/Users/jacquem/NetBeansProjects/MicromegasFramework/trunk/src/analyse/Gassiplex/Shower_profile/Profile_energy_loop.C

Go to the documentation of this file.

{gROOT->Reset();
  
  gSystem->Load("/home1/chefdevi/micromegasFrameWork/trunk/libMicro.so");
  

  bool display = 0;
  bool write_root_file = 1;

  

  TString particle[2] = {"electron","hadron"};
  TString nabsorber[6] = {"2","4","6","8","10","12"};



  for (int p=0;p<2;p++){

    for (int a=0;a<6;a++){


      TString filename= particle[p] + "_2gev_" + nabsorber[a] + "absorber.root";
      
      TString path = "/lapp_data/LC/Detecteurs/MicroMegas/data/TB2009/May/Rootfiles/dev/";
      //TString path = "/home1/chefdevi/micromegasFrameWork/trunk/src/analyse/Gassiplex/Shower_profile/";
      
      TFile *f = new TFile(path+filename);
      cout<<"reading File : "<<path<<filename<<endl;
      TTree *tree;
      tree  = (TTree*) f.Get("t1");cout<<"tree built\n";
      int nEvent = tree.GetEntries();
      MTEvent *evt = new MTEvent();
      TBranch *branch= tree->GetBranch("MTEvent");
      branch->SetAddress(&evt);
      cout <<"nEvent="<<nEvent<<endl;
      int NC = 672;



      /*Rootfile to save the histos*/
      
      if (write_root_file){
        TString filename2= "Root_histo_files/adc_" + filename;
        TFile * tf = new TFile(filename2,"recreate");}

      
      
      /*Distributions of the energy (or ADC)*/
      
      TH1F * hadc1 = new TH1F("hadc1","",400,0,1024*5);
      TH1F * hadc2 = new TH1F("hadc2","",400,0,1024*5);
      TH1F * hadc3 = new TH1F("hadc3","",400,0,1024*5);
      TH1F * hadc4 = new TH1F("hadc4","",400,0,1024*10);
      
      TH1F * hadc1_clean = new TH1F("hadc1_clean","",400,0,1024*5);
      TH1F * hadc2_clean = new TH1F("hadc2_clean","",400,0,1024*5);
      TH1F * hadc3_clean = new TH1F("hadc3_clean","",400,0,1024*5);
      TH1F * hadc4_clean = new TH1F("hadc4_clean","",400,0,1024*10);
      
      TH1F * hadc1_no_sat = new TH1F("hadc1_no_sat","",400,0,1024*5);
      TH1F * hadc2_no_sat = new TH1F("hadc2_no_sat","",400,0,1024*5);
      TH1F * hadc3_no_sat = new TH1F("hadc3_no_sat","",400,0,1024*5);
      TH1F * hadc4_no_sat = new TH1F("hadc4_no_sat","",400,0,1024*10);
      
      TH1F * hpedestal = new TH1F("hpedestal","",2048,-1024,1024);
      TH1F * hmean_pedestal = new TH1F("hmean_pedestal","",2048,-1024,1024);
      
      
      
      
      /*Set Line Width*/
      
      hadc1->SetLineWidth(2);
      hadc2->SetLineWidth(2);
      hadc3->SetLineWidth(2);
      hadc4->SetLineWidth(2);
      
      hadc1_clean->SetLineWidth(2);
      hadc2_clean->SetLineWidth(2);
      hadc3_clean->SetLineWidth(2);
      hadc4_clean->SetLineWidth(2);
      
      hadc1_no_sat->SetLineWidth(2);
      hadc2_no_sat->SetLineWidth(2);
      hadc3_no_sat->SetLineWidth(2);
      hadc4_no_sat->SetLineWidth(2);
      
      /*Set Line Color*/
      
      hadc2->SetLineColor(2);
      hadc3->SetLineColor(4);
      hadc4->SetLineColor(8);
      
      hadc2_clean->SetLineColor(2);
      hadc3_clean->SetLineColor(4);
      hadc4_clean->SetLineColor(8);
      
      hadc2_no_sat->SetLineColor(2);
      hadc3_no_sat->SetLineColor(4);
      hadc4_no_sat->SetLineColor(8);
      
      
      /*Set Line Style*/
      
      hadc1_clean->SetLineStyle(2);
      hadc2_clean->SetLineStyle(2);
      hadc3_clean->SetLineStyle(2);
      hadc4_clean->SetLineStyle(2);
      
      hadc1_no_sat->SetLineStyle(2);
      hadc2_no_sat->SetLineStyle(2);
      hadc3_no_sat->SetLineStyle(2);
      hadc4_no_sat->SetLineStyle(2);
      

      
      /*Define variables*/
      
      int xpos1,ypos1,adc1,sum_adc1;
      int xpos2,ypos2,adc2,sum_adc2;
      int xpos3,ypos3,adc3,sum_adc3;
      int xpos4,ypos4,adc4,sum_adc4;
      
      int threshold = 35;
      float nhit1=0,nhit2=0,nhit3=0,nhit4=0;
      int nclean = 0;
      bool clean = 0;
      int adc_sat = 900;
      int pedestal = 0;
      int correction = 0;
      int nunsaturated = 0;
      bool no_sat1 = 1,no_sat2 = 1;
      bool no_sat3 = 1,no_sat4 = 1;
      
      int event = 0;
      int crap_event = 0;
      
      
      
      
      
      /*Loop over all events*/
      
      for(int i=0;i<nEvent;i++){
        //for( int i=0;i<1;i++){
        //for( int i=event;i<event+1;i++){
        
        if (display){cout<<"Event # "<<i<<endl;}
        
        tree.GetEntry(i);
        NC = evt->GetNchannel();
        
        if (NC == 672){
          
          
          clean=0;
          no_sat1 = 1;no_sat2 = 1;no_sat3 = 1;no_sat4 = 1;
          nhit1=0,nhit2=0,nhit3=0,nhit4=0;
          sum_adc1=0;sum_adc2=0;sum_adc3=0;sum_adc4=0;
          hpedestal->Reset();
          
          
          /*Chamber beta 2.1 D*/
          
          for(int j=0;j<96;j++){
            
            MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
            xpos1 = channel.GetX();
            ypos1 = channel.GetY();
            adc1  = channel.GetValue();
            
            if (adc1>=threshold){
              
              sum_adc1+=adc1;
              /*no saturation if all adcs in chamber are below 1000*/
              if (adc1>adc_sat){no_sat1 = 0;}
              if (display){cout<<" x1="<<xpos1<<"  y1="<<ypos1<<"  adc1="<<adc1<<endl;}
              nhit1++;}}
          
          
          
          /*Chamber beta 2.1 E*/
          
          for(int j=0+96;j<96+96;j++){
            
            MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
            xpos2 = channel.GetX();
            ypos2 = channel.GetY();
            adc2  = channel.GetValue();
            
            if (adc2>=threshold){
              
              sum_adc2+=adc2;
              /*no saturation if all adcs in chamber are below 1000*/
              if (adc2>adc_sat){no_sat2 = 0;}
              if (display){cout<<" x2="<<xpos2<<"  y2="<<ypos2<<"  adc2="<<adc2<<endl;}
              nhit2++;}}
          
          
          
          /*Chamber beta 2.1 A*/
          
          for(int j=0+96+96;j<96+96+96;j++){
            
            MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
            xpos3 = channel.GetX();
            ypos3 = channel.GetY();
            adc3  = channel.GetValue();
            
            if (adc3>=threshold){
              
              sum_adc3+=adc3;
              /*no saturation if all adcs in chamber are below 1000*/
              if (adc3>adc_sat){no_sat3 = 0;}
              if (display){cout<<" x3="<<xpos3<<"  y3="<<ypos3<<"  adc3="<<adc3<<endl;}
              nhit3++;}}
          
          
          
          
          /*clean event if one and only one hit in each beta 2.1 chamber*/
          clean = (nhit1==1) && (nhit2==1) && (nhit3==1);
          
          
          
          
          /*chamber beta 2.4*/
          
          for(int j=287;j<672;j++){
            
            MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
            xpos4 = channel.GetX();
            ypos4 = channel.GetY();
            adc4  = channel.GetValue();

            if (adc4>=threshold){
              
              sum_adc4+=adc4;
              /*no saturation if all adcs in chamber are below 1000*/
              if (adc4>adc_sat){no_sat4 = 0;}
              if (display){cout<<" x4="<<xpos4<<"  y4="<<ypos4<<"  adc4="<<adc4<<endl;}
              nhit4++;}
            
            else{hpedestal->Fill(adc4);}}
          
          
          pedestal = hpedestal->GetMean();
          hmean_pedestal->Fill(pedestal);
          correction = nhit4 * pedestal;
          
          
          
          if (sum_adc1>0){hadc1->Fill(sum_adc1-correction);}
          if (sum_adc2>0){hadc2->Fill(sum_adc2-correction);}
          if (sum_adc3>0){hadc3->Fill(sum_adc3-correction);}
          if (sum_adc4>0){hadc4->Fill(sum_adc4-correction);nunsaturated++;}
          
          
          if (clean){
            
            if (sum_adc1>0){hadc1_clean->Fill(sum_adc1-correction);}
            if (sum_adc2>0){hadc2_clean->Fill(sum_adc2-correction);}
            if (sum_adc3>0){hadc3_clean->Fill(sum_adc3-correction);}
            if (sum_adc4>0){hadc4_clean->Fill(sum_adc4-correction);}
            
            if (no_sat1 && (sum_adc1>0)){hadc1_no_sat->Fill(sum_adc1-correction);}
            if (no_sat2 && (sum_adc2>0)){hadc2_no_sat->Fill(sum_adc2-correction);}
            if (no_sat3 && (sum_adc3>0)){hadc3_no_sat->Fill(sum_adc3-correction);}
            if (no_sat4 && (sum_adc4>0)){hadc4_no_sat->Fill(sum_adc4-correction);}
            
            if (display){cout<<"clean event!"<<endl;}
            
            nclean++;}
          
          
          if (display){cout<<nhit1<<" "<<nhit2<<" "<<nhit3<<" "<<nhit4<<endl;}}
        
        
        else{crap_event++;}}
      
      
      
      
      
      cout<<"Run "<<filename<<endl;
      cout<<"ntotal = "<<nEvent<<"  nclean = "<<nclean<<"  ncrap = "<<crap_event<<"  nunsaturated-4 = "<<nunsaturated<<endl;
      cout<<endl;
      
      
      
      /*save histo in a rootfile*/
      if (write_root_file){
        tf->Write();
        tf->Close();}}}









}

---