/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/bookkeeping/Complete_detector/Beam/complete_detector_profiles.C

Go to the documentation of this file.

    {
   gROOT->Reset();
   //gROOT->SetStyle("Bold");
   gStyle->SetPalette(1);
   gStyle->SetOptStat(0);
   gSystem->Load("/LC/Detecteurs/MicroMegas/Offline/micromegasFrameWork/trunk/lib/libMicro.so");


   /*Set rotated telescope*/
   bool rotate90 = false;

   /*m2 prototype time cut*/
   int dt_min = 6;
   int dt_max = 8;

   /*telescope energy cut*/
   int adc_cut = 30;



   int nhit_max = 1e9;




   bool display = 0;
   int nshower = 250;

   if (display)
     {
       TCanvas * c2 = new TCanvas("c2","",10,10,500,500);
       //c2->Divide(3,1);
     }




   //READ INFO FROM XML
   std::ifstream xml;
   xml.open("complete_detector.xml",ifstream::in);
   std::string line = "";
   std::string lineb = "";
   while(line.find("path=")==string::npos || lineb.find("<output")==string::npos)
     {
        lineb=line;
        xml>>line;
     }
   int start = line.find("\"")+1;
   int end =  line.rfind("\"")-6;
   TString file = line.substr(start,end);

   for (int i=0;i<5;i++){file.Chop();}
   file += "_MGT.root";




   bool noisy[9][32][32];
   for (int i=0;i<32;i++){
     for (int j=0;j<32;j++){
       for (int k=0;k<9;k++){
         noisy[k][i][j] = 0;}}}



     if (rotate90)
  {
     /*HV pad*/
     noisy[6][5][5] = 1;
     noisy[7][5][5] = 1;
     noisy[8][5][5] = 1;

     /*floating pad*/
     noisy[6][5][17] = 1;
     noisy[8][1][15] = 1;
  }
  else
  {
     /*HV pad*/
/*
     noisy[6][5][5-2] = 1;
     noisy[7][5][5-2] = 1;
     noisy[8][5][5-2] = 1;
*/
     /*floating pad*/
/*
     noisy[6][17][5-2] = 1;
     noisy[8][15][9-2] = 1;
*/
   }





   int nchannel = 0;
   int hardid = 0;
   int chipid = 0;
   int nhit = 0;
   int nhit_cut = 0;
   int nhit_cut_dac1 = 0;
   int nhit_cut_dac2 = 0;
   int neff = 0;
   float xpos = 0;  
   float ypos = 0;  
   float zpos = 0;  
   int adc = 0;
   int chbid = 0;
   int ntree = 0;

   Long64_t t,t0,t1,t2,t3,dt;
   
   TH2I * hxy_m2 = new TH2I("hxy_m2","Hit position distribution;y (cm);x (cm)",96,0,96,96,0,96);
   TH2I * hxy_m2_cut = new TH2I("hxy_m2_cut","Hit position distribution - time cut;y (cm);x (cm)",96,0,96,96,0,96);

   TH2I * hxy_m2_dac3 = new TH2I("hxy_m2_dac3","Hit position distribution;y (pad);x (pad)",96,0,96,96,0,96);
   TH2I * hxy_m2_cut_dac3 = new TH2I("hxy_m2_cut_dac3","Hit position distribution - time cut;y (pad);x (pad)",96,0,96,96,0,96);

   TH2I * hxy_m2_dac0 = new TH2I("hxy_m2_dac0","Hit position distribution;y (pad);x (pad)",96,0,96,96,0,96);
   TH2I * hxy_m2_cut_dac0 = new TH2I("hxy_m2_cut_dac0","Hit position distribution - time cut;y (pad);x (pad)",96,0,96,96,0,96);

   TH2I * hxy_m2_dac1 = new TH2I("hxy_m2_dac1","Hit position distribution;y (pad);x (pad)",96,0,96,96,0,96);
   TH2I * hxy_m2_cut_dac1 = new TH2I("hxy_m2_cut_dac1","Hit position distribution - time cut;y (pad);x (pad)",96,0,96,96,0,96);

   TH2I * hxy_m2_dac2 = new TH2I("hxy_m2_dac2","Hit position distribution;y (pad);x (pad)",96,0,96,96,0,96);
   TH2I * hxy_m2_cut_dac2 = new TH2I("hxy_m2_cut_dac2","Hit position distribution - time cut;y (pad);x (pad)",96,0,96,96,0,96);

   TH2F * hxy_ratio_10 = new TH2F("hxy_ratio_10","Ratio of Nhit med/low;y (pad);x (pad)",96,0,96,96,0,96);
   TH2F * hxy_ratio_20 = new TH2F("hxy_ratio_20","Ratio of Nhit high/low;y (pad);x (pad)",96,0,96,96,0,96);

   TH1F * hy_ratio_10 = new TH1F("hy_ratio_10","Ratio of Nhit med/low;y (pad);ratio",96,0,96);
   TH1F * hy_ratio_20 = new TH1F("hy_ratio_20","Ratio of Nhit high/low;y (pad);ratio",96,0,96);

   TH1I * hnhit = new TH1I("hnhit","Number of hit distribution;Nhits",10000,0,10000);

   TH1I * hnhit_cut = new TH1I("hnhit_cut","Number of hits above low threshold;Nhits",10000,0,10000);
   TH1I * hnhit_cut_dac1 = new TH1I("hnhit_cut_dac1","Number of hits above medium threshold;Nhits",10000,0,10000);
   TH1I * hnhit_cut_dac2 = new TH1I("hnhit_cut_dac2","Number of hits above high threshold;Nhits",10000,0,10000);

   TH1I * hdt = new TH1I("hdt","Time to readout disitribution;time (200 ns clock cycle)",200,0,200);
   TH2I * hxy_tel[9];
   TString title,name;

   for (int i=0;i<9;i++)
     {
       name = Form("hxy_tel_chamber_%i",i+1);
       title = Form("Hit pattern (adc>%i) chamber %i;y (cm),x (cm)",adc_cut,i+1);
       hxy_tel[i] = new TH2I(name,title,32,0,32,32,0,32);
     }



   
   TFile *f = new TFile(file);
   cout<<endl;
   cout<<"Reading file : "<<file<<endl;
   cout<<endl;

   TIter nextkey(f.GetListOfKeys());
   TKey *key;




   while (key = (TKey*)nextkey())
     {

       TTree *tree = (TTree*)key->ReadObj();                
       MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");

       int nEvent = tree.GetEntries();
       ntree++;
       cout<<"  Tree "<<ntree<<" contains "<<nEvent<<" events"<<endl;

       MTEvent *evt = new MTEvent();

       TBranch *branch= tree->GetBranch("MTEvent");
       branch->SetAddress(&evt);

       //for (int i=0;i<nEvent-1;i++)
         for (int i=0;i<20;i++)
         {

cout<<endl;
cout<<"event "<<i<<endl;

           //if (i%100==0){cout<<"Progress: "<<i/1./nEvent*100<<" \r"<<flush;}

           tree.GetEntry(i);
           nchannel = evt->GetNchannel();
           hnhit->Fill(nchannel);

           if (nchannel < nhit_max)
             {
               if (i>0)
                 {
                   hnhit_cut->Fill(nhit_cut);
                   hnhit_cut_dac1->Fill(nhit_cut_dac1);
                   hnhit_cut_dac2->Fill(nhit_cut_dac2);
                   
                   nhit_cut = 0;
                   nhit_cut_dac1 = 0;
                   nhit_cut_dac2 = 0;
                 }

               if (display)
                 {
                   hxy_m2_cut_dac0->Reset();
                   hxy_m2_cut_dac1->Reset();
                   hxy_m2_cut_dac2->Reset();
                   hxy_m2_cut_dac3->Reset();
                 }

               for (int j=0;j<nchannel;j++)
                 {
                   
                   MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
                   
                   chbid = channel->GetChamberId();
                   if(rotate90 && chbid!=10)
                     {
                       //xpos = channel->GetY()*1e-4;
                       //ypos = 10 - channel->GetX()*1e-4;
                       xpos = channel->GetY()*1e-4;
                       ypos = 10 - channel->GetX()*1e-4;
                     }
                   else
                     {
                       xpos = channel->GetX()*1e-4;
                       ypos = channel->GetY()*1e-4;
                     }             
                   
                   /*if m2*/
                   if (chbid == 10)
                     {
                       
                       hxy_m2->Fill(ypos,xpos);
                       t2 = channel->GetBcId_Dif();
                       t3 = channel->GetBcId_Hit();
                       dt = t2 - t3;
                       hdt->Fill(dt);
                       
                       int digit = channel->GetDigitalValue();
                       if (digit>=1)
                         {
                           hxy_m2_dac0->Fill(ypos,xpos);
                           hxy_m2_dac3->SetBinContent(ypos+1,xpos+1,1);
                           if (digit>=2)
                             {
                               hxy_m2_dac1->Fill(ypos,xpos);
                               hxy_m2_dac3->SetBinContent(ypos+1,xpos+1,2);
                               if (digit==3)
                                 {
                                   hxy_m2_dac2->Fill(ypos,xpos);
                                   hxy_m2_dac3->SetBinContent(ypos+1,xpos+1,3);
                                 }
                             }
                         }
                       
                       if (dt>=dt_min && dt<dt_max)
                         {
                           nhit_cut++;
                           hxy_m2_cut->Fill(ypos,xpos);
                           
                           if (digit>=1)
                             {
                               hxy_m2_cut_dac0->Fill(ypos,xpos);
                               hxy_m2_cut_dac3->SetBinContent(ypos+1,xpos+1,1);
                               if (digit>=2)
                                 {
                                   nhit_cut_dac1++;
                                   hxy_m2_cut_dac1->Fill(ypos,xpos);
                                   hxy_m2_cut_dac3->SetBinContent(ypos+1,xpos+1,2);
                                   if (digit==3)
                                     {
                                       nhit_cut_dac2++;
                                       hxy_m2_cut_dac2->Fill(ypos,xpos);
                                       hxy_m2_cut_dac3->SetBinContent(ypos+1,xpos+1,3);
                                     }
                                 }
                             }
                           
                           /*estimate of efficiency*/
                           if (nhit_cut == 1)
                             {
                               neff++;
                             }
                           
                         }
                     }
               
                   /*telescope*/
                   else
                     {                 
                       if (!noisy[chbid-1][ypos][xpos])
                         {
                           
                           adc = channel->GetAnalogValue();
                           
                           if (adc>=adc_cut)
                             {

                                if (chbid<=6){xpos=xpos*1e4/250;ypos=ypos*1e4/250;}

                                cout<<chbid<<" "<<xpos<<" "<<ypos<<" "<<adc<<endl;
        
                               hxy_tel[chbid-1]->Fill(ypos,xpos);
                             }
                         }
                     }
                 }

               if (nhit_cut>=nshower){
                 
                 if (display){
                   cout<<endl;
                   cout<<i<<" "<<nhit_cut<<endl;
                   cout<<endl;}
                 
                 if (display)
                   {
                     title = Form("Event %i",i+1);
                     hxy_m2_cut_dac3->SetTitle(title.Data());
                     hxy_m2_cut_dac3->Draw("zcol");
                     /*
                       c2->cd(1);
                       hxy_m2_cut_dac0->Draw("zcol");
                       c2->cd(2);
                       hxy_m2_cut_dac1->Draw("zcol");
                       c2->cd(3);
                       hxy_m2_cut_dac2->Draw("zcol");
                     */
                     c2->Update();
                     
                     title = Form("~/hxy_shower_375v_event_%i.png",i+1);
                     c2->Print(title.Data());
                     
                     title = Form("~/hxy_shower_375v_event_%i.C",i+1);
                     c2->Print(title.Data());
                     
                     /*
                       int temp = 0;
                       cin>>temp;
                       if (temp==9){break;}
                     */
                   }
               }
               
             }
         }
     }
    


   cout<<endl;
   cout<<"Nhot_m2 / Nevent = "<<neff<<" / "<<nEvent<<" = "<<neff*100./nEvent<<" %"<<endl;
   cout<<endl;




   //Ratios

   for (int i=0;i<96;i++)
     {
       int content0 = hxy_m2_cut_dac0->ProjectionX()->GetBinContent(i+1);
       
       if(content0!=0)
         {
           int content1 = hxy_m2_cut_dac1->ProjectionX()->GetBinContent(i+1);
           
           if(content1>10)
             {
               float ratio10 = 1.*content1/content0;
               float ratio10_err = ratio10 * sqrt(1./content0 + 1./content1);
               
               hy_ratio_10->SetBinContent(i+1,ratio10);
               hy_ratio_10->SetBinError(i+1,ratio10_err);
               
               int content2 = hxy_m2_cut_dac2->ProjectionX()->GetBinContent(i+1);
               
               if(content2>10)
                 {
                   float ratio20 = 1.*content2/content0;
                   float ratio20_err = ratio20 * sqrt(1./content0 + 1./content2);                  

                   hy_ratio_20->SetBinContent(i+1,ratio20);
                   hy_ratio_20->SetBinError(i+1,ratio20_err);
                 }
             }
         }
     }
   

for (int i=0;i<96;i++)
     {
       for (int j=0;j<96;j++)
         {
           int content0 = hxy_m2_cut_dac0->GetBinContent(i+1,j+1);
           
           if(content0!=0)
             {
               int content1 = hxy_m2_cut_dac1->GetBinContent(i+1,j+1);

               if(content1>10)
                 {
                   float ratio10 = 1.*content1/content0;

                   hxy_ratio_10->SetBinContent(i+1,j+1,ratio10);

                   int content2 = hxy_m2_cut_dac2->GetBinContent(i+1,j+1);

                   if(content2>10)
                     {
                       float ratio20 = 1.*content2/content0;

                       hxy_ratio_20->SetBinContent(i+1,j+1,ratio20);
                     }
                 }
             }
         }
     }







   int xmax_tel[9];
   int ymax_tel[9];
   int xmax_m2 = 0;
   int ymax_m2 = 0;
   int zmax = 0;


   hxy_m2_cut->GetMaximumBin(ymax_m2,xmax_m2,zmax);


   cout<<endl;
   cout<<endl;
   cout<<"Maximum of hit profile in m2 prototype"<<endl;
   cout<<endl;
   cout<<"    y-horiz. (cm) = "<<ymax_m2-1<<endl;
   cout<<"    x-verti. (cm) = "<<xmax_m2-1<<endl;
   cout<<endl;
   cout<<endl;

   cout<<"Maximum of hit profile in telescope pad chambers"<<endl;
   cout<<endl;

   for (int i=6;i<9;i++)
     {
 
       hxy_tel[i]->GetMaximumBin(ymax_tel[i],xmax_tel[i],zmax);

       cout<<"  Chamber "<<i+1<<endl;
       cout<<"    y-horiz. (cm) = "<<ymax_tel[i]-1<<endl;
       cout<<"    x-verti. (cm) = "<<xmax_tel[i]-1<<endl;
       cout<<endl;

     }




    


if (0){


   TCanvas * c0 = new TCanvas("c0","Prototype m2",10,10,1200,800);
   c0->Divide(3,2);

   c0->cd(1);
   hxy_m2->Draw("zcol");

   c0->cd(2);
   hdt->Draw("zcol");

   c0->cd(3);
   hxy_m2_cut->Draw("zcol");

   c0->cd(4);
   hxy_tel[6]->Draw("zcol");

   c0->cd(5);
   hxy_tel[7]->Draw("zcol");

   c0->cd(6);
   hxy_tel[8]->Draw("zcol");

   c0->Print("Plots/complete_detector_beam_profile_2d.png");


   TCanvas * c1 = new TCanvas("c1","Prototype m2",10,10,1200,600);
   c1->Divide(4,2);

   c1->cd(1);
   hxy_m2_cut->ProjectionX()->SetFillColor(2);
   hxy_m2_cut->ProjectionX()->SetLineColor(1);
   hxy_m2_cut->ProjectionX()->Draw();

   c1->cd(1+4);
   hxy_m2_cut->ProjectionY()->SetFillColor(3);
   hxy_m2_cut->ProjectionY()->SetLineColor(1);
   hxy_m2_cut->ProjectionY()->Draw();

   c1->cd(2);
   hxy_tel[6]->ProjectionX()->Draw();
   hxy_tel[7]->ProjectionX()->SetLineColor(1);
   hxy_tel[6]->ProjectionX()->SetFillColor(2);

   c1->cd(2+4);
   hxy_tel[6]->ProjectionY()->Draw();
   hxy_tel[6]->ProjectionY()->SetLineColor(1);
   hxy_tel[6]->ProjectionY()->SetFillColor(3);

   c1->cd(3);
   hxy_tel[7]->ProjectionX()->Draw();
   hxy_tel[7]->ProjectionX()->SetLineColor(1);
   hxy_tel[7]->ProjectionX()->SetFillColor(2);

   c1->cd(3+4);
   hxy_tel[7]->ProjectionY()->Draw();
   hxy_tel[7]->ProjectionY()->SetLineColor(1);
   hxy_tel[7]->ProjectionY()->SetFillColor(3);

   c1->cd(4);
   hxy_tel[8]->ProjectionX()->Draw();
   hxy_tel[8]->ProjectionX()->SetLineColor(1);
   hxy_tel[8]->ProjectionX()->SetFillColor(2);

   c1->cd(4+4);
   hxy_tel[8]->ProjectionY()->Draw();
   hxy_tel[8]->ProjectionY()->SetLineColor(1);
   hxy_tel[8]->ProjectionY()->SetFillColor(3);


   c1->Print("Plots/complete_detector_beam_profile_1d.png");


   cout<<endl;}

}

}

---