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

Go to the documentation of this file.

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


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

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

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



   //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] = 1;
     noisy[7][5][5] = 1;
     noisy[8][5][5] = 1;

     /*floating pad*/
     noisy[6][17][5] = 1;
     noisy[8][15][9] = 1;
   }





   int nchannel = 0;
   int hardid = 0;
   int chipid = 0;
   int nhit = 0;
   int nhit_cut = 0;
   int neff = 0;
   int xpos = 0;  
   int ypos = 0;  
   int 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);
   TH1I * hnhit = new TH1I("hnhit","Number of hit distribution;Nhits",10000,0,10000);
   TH1I * hnhit_cut = new TH1I("hnhit_cut","Number of hit distribution - time cut;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++)
         {

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

           tree.GetEntry(i);
           nchannel = evt->GetNchannel();
           hnhit->Fill(nchannel);
           if (i>0){hnhit_cut->Fill(nhit_cut);nhit_cut = 0;}


           for (int j=0;j<nchannel;j++)
             {

               MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);

               chbid = channel->GetChamberId();
               if(rotate90 && chbid!=10)
                 {
                   //xpos = channel->GetY();
                   //ypos = 10 - channel->GetX();
                   xpos = channel->GetY();
                   ypos = 10 - channel->GetX();
                 }
               else
                 {
                   xpos = channel->GetX();
                   ypos = channel->GetY();
                 }
               //cout<<"chamber id = "<<chbid<<"x="<<xpos<<" y="<<ypos<<endl;


               /*if m2*/
               if (chbid == 10)
                 {
                        if (xpos==35 && ypos==65) continue;
                        if (xpos==82 && ypos==63) continue;
                        if (xpos==2  && ypos==89) continue;
                   hxy_m2->Fill(ypos,xpos);
                   t2 = channel->GetBcId_Dif();
                   t3 = channel->GetBcId_Hit();
                   dt = t2 - t3;
                   hdt->Fill(dt);
         
                   if (dt>=dt_min && dt<dt_max)
                     {
                       nhit_cut++;
                       hxy_m2_cut->Fill(ypos,xpos);

                       /*estimate of efficiency*/
                       if (nhit_cut == 1)
                         {

                           neff++;

                         }

                     }
                 }
               
               /*telescope*/
               else
                 {

                   if (!noisy[chbid-1][ypos][xpos])
                     {

                       adc = channel->GetAnalogValue();
                       
                       if (adc>=adc_cut)
                         {
                           hxy_tel[chbid-1]->Fill(ypos,xpos);
                         }
                     }
                 }
             }
         }
     }
    


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




   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;

     }

   cout<<endl;



    }





   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;}

---