Micromegas LAPP reconstruction

int main	(	int	argc,
		char **	argv
	)
Definition at line 43 of file myfirstprogam.cpp.
00043                               {
00044 
00045   Arguments arg(argc,argv," run1 [run2] max [run20] -o outputrootfilename");
00046 
00047   if ( arg.getNbOfArgs() < 1) {
00048     arg.usage();
00049     exit(1);
00050   }
00051 
00052 
00053 
00054   int Nsample = 3;
00055 
00056   int z_sample[3] = {14,24,39};
00057 
00058   float calib[3] = {1.8/1.641,1.8/1.587,1.8/1.611};
00059 
00060   int nhitcutstart = 9;
00061 
00062   int nhit_calo_cut = 220;  
00063 
00064 
00065   int NbOfRuns = arg.getNbOfArgs();
00066 
00067   TString runnumber[20];
00068   TString outfilename = "temp.root";
00069 
00070   for (int nr=0 ; nr < NbOfRuns ; nr ++ )
00071     {
00072       runnumber[nr] = arg.getArgument(nr+1); 
00073     }
00074 
00075   if ( arg.getOption("-o").size() != 0 )
00076     {
00077       outfilename  = arg.getOption("-o");
00078     }
00079 
00080   if ( arg.getOption("-c").size() != 0 )
00081     {
00082       nhitcutstart  = atoi(arg.getOption("-c").c_str());
00083 
00084       TString temp = arg.getOption("-c");
00085       for (int i=0;i<5;i++){outfilename.Chop();}
00086       outfilename += "_cut" + temp + ".root";
00087     }
00088 
00089 
00090   if ( arg.getOption("-n").size() != 0 )
00091     {
00092       nhit_calo_cut = atoi(arg.getOption("-n").c_str());
00093       cout<<"Nhit_calo_cut = "<<nhit_calo_cut<<endl;
00094     }
00095 
00096    
00097 
00098 
00099   int Nlayer = 50;
00100 
00101   int nshower = 0;
00102   int nhshower = 0;
00103   int nmuon = 0;
00104   int npenetrating_muon = 0;
00105 
00106   int nstart = 0;
00107 
00108   int nhit_evt = 0;
00109   int nhit_cumul = 0;
00110 
00111   int zstart = -1;
00112 
00113   int nevt_tot = 0;
00114 
00115   int ndouble = 0;
00116 
00117   int nbHits_prev = -1;
00118 
00119 
00120 
00121 
00122 
00123   TH1F * hdensity = new TH1F("hdensity","",200,0,100);
00124   TH2F * hnhit_density = new TH2F("hnhit_density","",1000,0,2000,200,0,100);
00125 
00126   TH1F * hrhadron = new TH1F("hrhadron","",1000,0,10);
00127   TH2F * hnhit_rhadron = new TH2F("hnhit_rhadron","",1000,0,2000,200,0,10);
00128 
00129   TH2I * hbug = new TH2I("hbug","",100,0,1000,100,0,1000);
00130 
00131   TProfile * hlong1 = new TProfile("hlong1","",50,0,50,0,150,"");
00132   TProfile * hlong2 = new TProfile("hlong2","",50,0,50,0,150,"");
00133   TProfile * hlong3 = new TProfile("hlong3","",50,0,50,0,150,"");
00134 
00135   TGraphErrors * tgmip = new TGraphErrors();
00136   TGraphErrors * tgshower = new TGraphErrors();
00137   TGraphErrors * tgnoise = new TGraphErrors();
00138 
00139   tgmip->SetName("tgmip");
00140   tgshower->SetName("tgshower");
00141   tgnoise->SetName("tgnoise");
00142 
00143   TH1I * hnhit_calo = new TH1I("hnhit_calo","",1000,0,10000);
00144   TH1I * hnhit_evt  = new TH1I("hnhit_evt","all showers",500,0,5000);
00145   TH1I * hnhit_evt2 = new TH1I("hnhit_evt2","all showers starting wthin 5 first planes",500,0,5000);
00146   TH1I * hnhit_evt3 = new TH1I("hnhit_evt3","all showers without MIP part",500,0,5000);
00147   TH1I * hnhit_evt4 = new TH1I("hnhit_evt4","all showers without MIP part starting wthin 5 first planes",500,0,5000);
00148 
00149   TH1I * hstart = new TH1I("hstart","",Nlayer,0,Nlayer);
00150 
00151   TH1I * hnhit_layer = new TH1I("hnhit_layer","",200,0,200);
00152 
00153   TH1I * hz_nhit = new TH1I("hz_nhit","",50,0,50);
00154   TH1I * hz_nhit_cumul = new TH1I("hz_nhit_cumul","",50,0,50);
00155 
00156   TH1F * hz_rmsx = new TH1F("hz_rmsx","",50,0,50);
00157   TH1F * hz_rmsy = new TH1F("hz_rmsy","",50,0,50);
00158 
00159   TH1F * hx_evt = new TH1F("hx_evt","",96,0,96*1.04);
00160   TH1F * hy_evt = new TH1F("hy_evt","",96,0,96*1.04);
00161 
00162   TH2I * hnhit_zstart = new TH2I("hnhit_zstart","",50,0,50,400,0,2000);
00163 
00164 
00165   TString name,title;
00166   TH2F * hxy[Nlayer];
00167   TH1F * hx[Nlayer];
00168   TH1F * hy[Nlayer];
00169   TH1I * hnhit_mip[Nlayer];
00170   TH1I * hnhit_shower[Nlayer];
00171   TH1I * hnhit_shower2[Nlayer];
00172   TH1I * hnhit_noise[Nlayer];
00173   TH1I * hdt[Nlayer];
00174   int nhit[Nlayer];
00175   int nnoise[Nlayer];
00176 
00177   for (int i=0;i<Nlayer;i++)
00178     {
00179       nhit[i] = 0;
00180       nnoise[i] = 0;
00181 
00182       name = Form("hnhit_mip_%i",i);
00183       title = Form("Number of CaloHit - chb %i",i);
00184       hnhit_mip[i]= new TH1I(name,title,100,0,100);
00185 
00186       name = Form("hnhit_shower_%i",i);
00187       title = Form("Number of CaloHit - chb %i",i);
00188       hnhit_shower[i]= new TH1I(name,title,1000,0,1000);
00189 
00190       name = Form("hnhit_shower2_%i",i);
00191       title = Form("Number of CaloHit - chb %i",i);
00192       hnhit_shower2[i]= new TH1I(name,title,1000,0,1000);
00193 
00194       name = Form("hnhit_noise_%i",i);
00195       title = Form("Number of NoiseHit - chb %i",i);
00196       hnhit_noise[i]= new TH1I(name,title,100,0,100);
00197 
00198       name = Form("hdt_%i",i);
00199       title = Form("Delta_t - chb %i",i);
00200       hdt[i]= new TH1I(name,title,21,-10,10);
00201 
00202       name = Form("hxy_%i",i);
00203       title = Form("CaloHit position - chb %i",i);
00204 
00205       if (i<48)
00206         {
00207           hxy[i] = new TH2F(name,title,96,0,96*1.04,96,0,96*1.04);
00208 
00209           name = Form("hx_%i",i);
00210           title = Form("X - chb %i",i);
00211           hx[i] = new TH1F(name,title,96,0,96*1.04);
00212 
00213           name = Form("hy_%i",i);
00214           title = Form("Y - chb %i",i);
00215           hy[i] = new TH1F(name,title,96,0,96*1.04);
00216         }
00217       else
00218         {
00219           hxy[i] = new TH2F(name,title,194,0,97,193,0,96.5);
00220 
00221           name = Form("hx_%i",i);
00222           title = Form("X - chb %i",i);
00223           hx[i] = new TH1F(name,title,194,0,97);
00224 
00225           name = Form("hy_%i",i);
00226           title = Form("Y - chb %i",i);
00227           hy[i] = new TH1F(name,title,193,0,96.5);
00228         }
00229     }
00230 
00231 
00232 
00233 
00234 
00235 
00236 
00237   TChain chain("T");
00238 
00239   for ( int nr=0; nr< NbOfRuns; nr++){
00240 
00241     for (int i=0 ; i<4 ; i++)
00242       {
00243         TString SubFile = Form("_I%i",i);
00244         TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_" + runnumber[nr] + SubFile + "_0.slcio.root");
00245         chain.Add(fname);
00246         cout << fname << endl;
00247       }
00248   }
00249 
00250   TObjArray *fileElements=chain.GetListOfFiles();
00251   TIter next(fileElements);
00252   TChainElement *chEl=0;
00253   int NchainFiles=0;
00254 
00255   while (( chEl=(TChainElement*)next() ))
00256     {
00257       TFile f(chEl->GetTitle());
00258       TIter nextkey(f.GetListOfKeys());
00259       TKey *key;
00260 
00261       int NbTrees = 0;
00262 
00263       while (key = (TKey*)nextkey()) 
00264         {
00265 
00266           TString keyname = key->GetName();
00267           //if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
00268           //else{
00269 
00270             TTree *tree = (TTree*)key->ReadObj();                
00271             cout << "---------------------------- NEW TTREE -----------------"<< endl;
00272             int nbEvt =0;
00273             int nbHit = 0;
00274             NbTrees++;
00275             CaloRun *run=NULL;
00276 
00277             try
00278               {
00279                 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
00280                 if( run != NULL )
00281                   {
00282                     cout << "Run[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ;
00283                   }
00284               }
00285 
00286             catch ( ... ) {}
00287 
00288             int RunId=1;
00289             CaloEvent *evt =  new CaloEvent();
00290             TBranch *branch= tree->GetBranch("CaloEvent");
00291             branch->SetAddress(&evt);
00292 
00293             CaloHit* hit =NULL;
00294             int nbEntries = tree->GetEntries();
00295 
00296 
00297 
00298 
00299             //nbEntries = 1000;
00300             for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
00301               {
00302                 tree->GetEntry(evtNum);
00303                 nbEvt++;
00304 
00305                 int nhit_seperate1 = 0;
00306                 int nhit_seperate2 = 0;
00307                 float rhadron = 0;
00308 
00309                 int nhit_evt = 0;
00310                 int nhit_evt2 = 0;
00311                 int nhit_evt3 = 0;
00312                 int nhit_evt4 = 0;
00313                 
00314                 int nbHits = evt->GetNHits();
00315                 int maxDeltat = evt->GetDeltaTmax();
00316 
00317                 int nhotlayer = 0;
00318                 float density = 0;
00319 
00320                 int nhot_rear  = 0;
00321                 int nhot_front = 0;
00322 
00323                 for(int i=0 ; i < Nlayer ; i++)
00324                   {
00325                     nhit[i] = 0;
00326                     nnoise[i] = 0;
00327                   }
00328                 
00329                 zstart = -1;
00330 
00331 
00332                 if (evtNum!=0){hbug->Fill(nbHits,nbHits_prev);}
00333 
00334                 nevt_tot++;
00335 
00336                     
00337                     //cout<<"Evt "<<evtNum<<" Nhit = "<<nbHits<<" - Nhit_prev = "<<nbHits_prev<<endl;
00338                     cout<<"evt "<<evtNum<<" -> nhit = "<<nbHits<<" \r"<<flush;
00339 
00340                     hnhit_calo->Fill(nbHits);
00341 
00342 
00343 
00344                     /******** time cut ******/
00345 
00346                     for(int i=0 ; i < nbHits ; i++)
00347                       {
00348                         hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
00349                         
00350                         int deltat = hit->GetDeltaT();
00351                         int zSlot = hit->GetLayer();
00352                         
00353                         float xpos =  hit->GetX()*1e-4;
00354                         float ypos =  hit->GetY()*1e-4;
00355                         
00356                         hxy[zSlot]->Fill(xpos,ypos);
00357                         
00358                         hx[zSlot]->Fill(xpos);
00359                         hy[zSlot]->Fill(ypos);
00360                         
00361                         hx_evt->Fill(xpos);
00362                         hy_evt->Fill(ypos);
00363 
00364                         hdt[zSlot]->Fill(maxDeltat-deltat);
00365                         
00366                         if (zSlot >= 48)
00367                           {
00368                             if  (fabs(maxDeltat-deltat) <= 2)
00369                               {
00370                                 nhit[zSlot]++;
00371                               }
00372                             else
00373                               {
00374                                 nnoise[zSlot]++;
00375                               }
00376                           }
00377                         else
00378                           {
00379                             if  (fabs(maxDeltat-deltat) <= 1)
00380                               {
00381                                 nhit[zSlot]++;
00382                               }
00383                             else
00384                               {
00385                                 nnoise[zSlot]++;
00386                               }
00387                           }
00388                       }
00389 
00390                 
00391 
00392 
00393                     /********** sum of hits ************/
00394 
00395 
00396                     for(int i=0 ; i < Nlayer ; i++)
00397                       {
00398                         hnhit_noise[i]->Fill(nnoise[i]);
00399 
00400                         nhit_evt += nhit[i];
00401                         
00402                         if (nhit[i] != 0)
00403                           {
00404                             hnhit_layer->Fill(nhit[i]);
00405 
00406                             density += nhit[i];
00407                             nhotlayer++;
00408 
00409                             if (i<5){nhit_seperate1++;}
00410                             if (i>=15){nhit_seperate2++;}
00411 
00412                             if (i<10){nhot_front++;}
00413                             if (i>=19 && i<30){nhot_rear++;}
00414                           }
00415                       }
00416 
00417 
00418 
00419 
00420 
00421 
00422                     /************ calculate density ****************/
00423 
00424                     if (nhotlayer != 0)
00425                       {
00426                         density /= (1.*nhotlayer);
00427                         hdensity->Fill(density);
00428                         hnhit_density->Fill(nhit_evt,density);
00429                       }
00430 
00431                     if (nhit_seperate2 != 0)
00432                       {
00433                         rhadron = 1. * nhit_seperate1 / nhit_seperate2;
00434                       }
00435 
00436 
00437 
00438 
00439                     /******** select muons *******/
00440 
00441                     if (density !=0 && density < 5)
00442                       {
00443                         /********** select penetrating ************/
00444 
00445                         nmuon++;
00446 
00447                         if (nhot_front >= 6 && nhot_rear >= 6)
00448                           {
00449 
00450                             npenetrating_muon++;
00451 
00452                             for(int i=0 ; i < Nlayer ; i++)
00453                               {
00454                                 hnhit_mip[i]->Fill(nhit[i]);
00455                               }
00456                           }
00457                       }
00458 
00459 
00460 
00461 
00462                     /******** select showers *******/
00463 
00464                     if (density >= 5 && nhit_evt >= nhit_calo_cut)
00465                       {        
00466 
00467                         hrhadron->Fill(rhadron);
00468                         hnhit_rhadron->Fill(nhit_evt,rhadron);
00469 
00470                         nshower++;
00471 
00472                         /********** select pions - reject electrons ************/
00473 
00474                         if (rhadron <= 0.4)
00475                           {                     
00476                             nhshower++;
00477                         
00478                             /******** find shower start *********/
00479                             
00480                             for(int i=0 ; i < Nlayer-2 ; i++)
00481                               {
00482                                 
00483                                 /******** rising number of hits in 3 consecutive layer *********/
00484                                 
00485                                 if ( (nhit[i+2] > nhit[i+1]) && (nhit[i+1] > nhit[i]) )
00486                                   {
00487                                     for (int k=0 ; k<3 ; k++)
00488                                       {
00489                                         
00490                                         /******* shower begins when nhit >= 4 ***********/
00491                                         
00492                                         if (nhit[i+k] >= nhitcutstart)
00493                                           {
00494                                             nstart++;
00495                                             zstart = i+k;
00496                                             
00497                                             hstart->Fill(zstart);
00498                                             k = 3;
00499                                             i = Nlayer;
00500                                           }
00501                                       }
00502                                   }
00503                               }
00504                             
00505 
00506 
00507                             /********* fill profiles *********/
00508                             
00509                             if (zstart != -1)
00510                               {
00511                                 
00512                                 hnhit_zstart->Fill(zstart,nhit_evt);
00513                                 
00514                                 
00515                                 for(int i=0 ; i < Nlayer ; i++)
00516                                   {
00517                                     hnhit_shower[i]->Fill(nhit[i]);
00518                                     
00519                                     hz_nhit->Fill(i,nhit[i]);
00520                                     hlong1->Fill(i,nhit[i]);
00521                                     
00522                                     //all showers starting wthin 5 first planes
00523                                     if (zstart < 5)
00524                                       {
00525                                         nhit_evt2 += nhit[i];
00526                                       }
00527                                     
00528                                     //all showers without MIP part
00529                                     if (i >= zstart && i <= 45)
00530                                       {
00531                                         nhit_evt3 += nhit[i];
00532                                         
00533                                         //all showers without MIP part starting wthin 5 first planes
00534                                         if (zstart < 5)
00535                                           {
00536                                             nhit_evt4 += nhit[i];
00537                                           }
00538                                       }
00539                                   }
00540                                 
00541                                 
00542                                 //all showers
00543                                 hnhit_evt->Fill(nhit_evt);
00544                                 
00545                                 //all showers starting wthin 5 first planes
00546                                 if (nhit_evt2 != 0){hnhit_evt2->Fill(nhit_evt2);}
00547                                 
00548                                 //all showers without MIP part
00549                                 hnhit_evt3->Fill(nhit_evt3);
00550                                 
00551                                 //all showers without MIP part starting wthin 5 first planes
00552                                 if (nhit_evt4 != 0){hnhit_evt4->Fill(nhit_evt4);}
00553 
00554                   
00555 
00556                                 for (int k=0 ; k<Nsample ; k++)
00557                                   {
00558                                     if (zstart <= z_sample[k])
00559                                       {
00560                                         hnhit_shower2[z_sample[k] - zstart]->Fill(nhit[z_sample[k]]);
00561                                         hlong2->Fill(z_sample[k] - zstart,nhit[z_sample[k]]);
00562                                         hlong3->Fill(z_sample[k] - zstart,nhit[z_sample[k]] * calib[k]);
00563                                       }
00564                                   }
00565                               }
00566                           }
00567                       }
00568 
00569                     if (nbHits_prev == nbHits){ndouble++;}
00570 
00571                 nbHits_prev = nbHits;
00572               }          
00573         }
00574     }
00575 
00576 
00577 
00578   cout<<endl;
00579   cout<<endl;
00580   cout<<"Nevt_tot = "<<nevt_tot<<endl;
00581   cout<<"Ndouble = "<<ndouble<<endl;
00582   cout<<endl;
00583   cout<<"Nshower = "<<nshower<<endl;
00584   cout<<"Nhadronshower = "<<nhshower<<endl;
00585   cout<<"Nstart = "<<nstart<<endl;
00586   cout<<endl;
00587   cout<<"Nmuon = "<<nmuon<<endl;
00588   cout<<"Npenetrating_muon = "<<npenetrating_muon<<endl;
00589   cout<<endl;
00590 
00591 
00592 
00593 
00594 
00595   cout<<endl;
00596   cout<<"Write to TFile "<<outfilename<<endl;
00597   cout<<endl;
00598   
00599   TFile * tf = new TFile(outfilename,"RECREATE");
00600   
00601   for (int i=0;i<Nlayer;i++)
00602     {
00603       tgmip->SetPoint(i,i,hnhit_mip[i]->GetMean());
00604       tgmip->SetPointError(i,0,hnhit_mip[i]->GetRMS());
00605       
00606       tgshower->SetPoint(i,i,hnhit_shower[i]->GetMean());
00607       tgshower->SetPointError(i,0,hnhit_shower[i]->GetRMS());
00608       
00609       tgnoise->SetPoint(i,i,hnhit_noise[i]->GetMean());
00610       tgnoise->SetPointError(i,0,hnhit_noise[i]->GetRMS());
00611       
00612       hxy[i]->Write();
00613       hdt[i]->Write();
00614       hnhit_mip[i]->Write();
00615       hnhit_shower[i]->Write();
00616       hnhit_shower2[i]->Write();
00617       hnhit_noise[i]->Write();
00618     }
00619   
00620   tgmip->SetMarkerStyle(20);
00621   tgmip->SetMarkerColor(2);
00622   
00623   tgshower->SetMarkerStyle(20);
00624   tgshower->SetMarkerColor(3);
00625   
00626   tgnoise->SetMarkerStyle(20);
00627   tgnoise->SetMarkerColor(4);
00628   
00629   tgmip->Write();
00630   tgshower->Write();
00631   tgnoise->Write();
00632 
00633   hrhadron->Write();
00634   hnhit_rhadron->Write();
00635   hnhit_zstart->Write();
00636   hbug->Write();
00637   hdensity->Write();
00638   hnhit_density->Write();
00639   hstart->Write();
00640   hlong1->Write();
00641   hlong2->Write();
00642   hlong3->Write();
00643   hx_evt->Write();
00644   hy_evt->Write();
00645   hnhit_evt->Write();
00646   hnhit_evt2->Write();
00647   hnhit_evt3->Write();
00648   hnhit_evt4->Write();
00649   hnhit_layer->Write();
00650   hz_nhit->Write();
00651   hz_nhit_cumul->Write();
00652   hz_rmsx->Write();
00653   hz_rmsy->Write();
00654   hnhit_calo->Write();
00655   
00656   tf->Close();
00657 
00658   
00659 }
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Iro/myfirstprogam.cpp File Reference

Functions

Function Documentation
