/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/bookkeeping/Complete_detector/Beam3/complete_detector_efficiency.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");


   /*Offset m2 prototype - telescope*/
   int offset_y = 29;
   int offset_x = 33;

   /*1D half length of search area*/
   int search = 5;

   /*m2 prototype time cut*/
   int dt_min = 5;
   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;}}}

   /*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 nhot_m2 = 0;
   int nhot_tel = 0;
   int nsuper_hot_tel = 0;
   int nhit[10];
   bool hot[10];
   bool super_hot[10];
   for (int i=0;i<10;i++){nhit[i] = 0;hot[i] = 0;}

   int nchannel = 0;
   int hardid = 0;
   int chipid = 0;
   int xpos = 0;  
   int ypos = 0;  
   int zpos = 0;  
   int adc = 0;
   int chbid = 0;
   int ntree = 0;
   int mult = 0;
   int ntrack = 0;
   int neff = 0;
   int track_x = 0;
   int track_y = 0;
   int target_x = 0;
   int target_y = 0;
   int dt_tp_max = 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_target = new TH2I("hxy_m2_target","Extrapolated 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_cut_track = new TH2I("hxy_m2_cut_track","Hit position distribution when track in telescope- time cut;y (cm);x (cm)",96,0,96,96,0,96);
   TH2I * hxy_m2_cut_temp = new TH2I("hxy_m2_cut_temp","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 * hmult_m2 = new TH1I("hmult_m2","Number of hit in target area in m2;Nhit",10,0,10);
   TH1I * hnhot_tel = new TH1I("hnhot_tel","Number of hot telescope chambers (Nhit>0);Nhot",9,0,9);
   TH1I * hnsuper_hot_tel = new TH1I("hnsuper_hot_tel","Number of super hot telescope chambers (Nhit=1);Nhot",9,0,9);
   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);
   TH1I * hdt_first = new TH1I("hdt_first","Time to readout of first hit;time (200 ns clock cycle)",20,0,20);
   TH1I * hdt_width = new TH1I("hdt_width","Time to first hit disitribution;time (200 ns clock cycle)",20,0,20);
   TH1I * hdt_tp = new TH1I("hdt_tp","Time to readout disitribution;time (200 ns clock cycle)",20,0,20);
   TH1I * hx_res = new TH1I("hx_res","Vertical residuals Prototype - telescope;Xm2 - Xtarget (cm)",51,-25,25);
   TH1I * hy_res = new TH1I("hy_res","Horizontal residuals Prototype - telescope;Ym2 - Ytarget (cm)",51,-25,25);
   TH1F * hr_res = new TH1F("hr_res","Residuals Prototype - telescope;d(m2 - target) (cm)",100,0,25);
   TH2I * hx_cor = new TH2I("hx_cor","Correlation Prototype - telescope;Xtarget (cm);Xm2 (cm)",100,0,100,100,0,100);
   TH2I * hy_cor = new TH2I("hy_cor","Correlation Prototype - telescope;Ytarget (cm);Ym2 (cm)",100,0,100,100,0,100);
   TH2I * hxy_tel_pad_temp = new TH2I("hxy_tel_pad_temp","Hit position distribution in telescope pad chamber for super hot events (Nhit=1);y (cm);x (cm)",32,0,32,32,0,32);

   TH1I * hnhit_tel[9];
   TH2I * hxy_tel[9];
   TH2I * hxy_tel_temp[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);

       name = Form("hxy_tel_temp_chamber_%i",i+1);
       title = Form("Hit pattern (adc>%i) chamber %i;y (cm);x (cm)",adc_cut,i+1);
       hxy_tel_temp[i] = new TH2I(name,title,32,0,32,32,0,32);

       name = Form("hnhit_tel_chamber_%i",i+1);
       title = Form("Number of hits (adc>%i) in chamber %i;Nhit",adc_cut,i+1);
       hnhit_tel[i] = new TH1I(name,title,96,0,96);
     }



   
   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();

           nhot_m2 = 0;
           nhot_tel = 0;
           nsuper_hot_tel = 0;
           for (int k=0;k<10;k++){nhit[k] = 0;super_hot[k] = 0;}

           /*Look in telescope*/
           for (int j=0;j<nchannel;j++)
             {

               MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
               chbid = channel->GetChamberId();
               xpos = channel->GetX()*1e-4;
               ypos = channel->GetY()*1e-4;

               if (chbid < 10)
                 {
                   if (!noisy[chbid-1][ypos][xpos])
                     {
                       adc = channel->GetAnalogValue();
                       
                       if (adc>=adc_cut)
                         {                 
                           nhit[chbid-1]++;

                           hxy_tel[chbid-1]->Fill(ypos,xpos);
                           hxy_tel_temp[chbid-1]->Fill(ypos,xpos);

                           if (chbid == 7 || chbid == 8 || chbid == 9)
                             {

                               hxy_tel_pad_temp->Fill(ypos,xpos);

                             }

                           if (chbid == 8)
                             {

                               track_x = xpos;
                               track_y = ypos;

                             }

                         }
                     }
                 }
             }

           /*Find tracks in pad telescope*/

           for (int k=6;k<9;k++)
             {

               if (nhit[k] != 0)
                 {
                   hot[k] = 1;

                   if (nhit[k] == 1)
                     {
                       super_hot[k] = 1;
                     }
                 }

               hnhit_tel[k]->Fill(nhit[k]);

               if (hot[k])
                 {
                   nhot_tel++;
                 }

               if (super_hot[k])
                 {
                   nsuper_hot_tel++;
                 }
             }

           hnhot_tel->Fill(nhot_tel);
           hnsuper_hot_tel->Fill(nsuper_hot_tel);


           /*If track, look in m2 prototype*/
           if (nsuper_hot_tel == 3)
             {

               if (hxy_tel_pad_temp->GetRMS(1) == 0 && hxy_tel_pad_temp->GetRMS(2) == 0){

               ntrack++;

               /*Extrapolate track to m2 prototype*/
               target_x = track_x + offset_x;
               target_y = track_y + offset_y;
               hxy_m2_target->Fill(target_y,target_x);

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

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

                   if (chbid == 10)
                     {

                       t2 = channel->GetBcId_Dif();
                       t3 = channel->GetBcId_Hit();
                       dt = t2 - t3;
                       hdt->Fill(dt);
                       hdt_tp->Fill(dt);                       

                       xpos = channel->GetX()*1e-4;
                       ypos = channel->GetY()*1e-4;
                       hxy_m2->Fill(ypos,xpos);

                       if (dt>=dt_min && dt<dt_max)
                         {
                           nhit[chbid-1]++;
                           hxy_m2_cut->Fill(ypos,xpos);
                           hxy_m2_cut_temp->Fill(ypos,xpos);

                           hx_cor->Fill(target_x,xpos);
                           hy_cor->Fill(target_y,ypos);

                           hr_res->Fill(sqrt(pow(xpos-target_x,2) + pow(ypos-target_y,2)));
                           hx_res->Fill(xpos - target_x);
                           hy_res->Fill(ypos - target_y);

                           /*look in target area*/
                           if (fabs(xpos - target_x)<=search && fabs(ypos - target_y)<=search)
                             {

                               mult++;
                               hxy_m2_cut_track->Fill(ypos,xpos);

                             }
                         }
                     }
                 }

               for (k=10;k>0;k--){if (hdt_tp->GetBinContent(k+1)!=0){dt_tp_max = k;k = 0;}}
               hdt_first->Fill(dt_tp_max);
               for (k=10;k>0;k--){if (hdt_tp->GetBinContent(k+1)!=0){hdt_width->Fill(dt_tp_max-k);}}
               hdt_tp->Reset();

               if (mult != 0)
                 {
                   neff++;
                   hmult_m2->Fill(mult);
                   mult = 0;
                 }

               hxy_m2_cut_temp->Reset();

               }
             }

           hxy_tel_pad_temp->Reset();

         }
     }
    



   cout<<endl;
   cout<<endl;
   cout<<"Number of tracks in telescope = "<<ntrack<<endl;
   cout<<"Number of time the m2 prototype was efficienct (search area of "<<search*2+1<<"*"<<search*2+1<<") = "<<neff<<endl;
   cout<<endl;
   cout<<"Global efficiency = "<<neff<<" / "<<ntrack<<" = "<<neff*1./ntrack*100.<<" %"<<endl;
   cout<<endl;
   cout<<"Hit multiplicity  = "<<hmult_m2->GetMean()<<endl;
   cout<<endl;
   cout<<endl;





}

---