Micromegas LAPP reconstruction

int main	(	int	argc,
		char **	argv
	)
Definition at line 48 of file profile4layers_simu.cpp.
00048                               {
00049 
00050   Arguments arg(argc,argv," inputrootfilename outputrootfilename hitcutstart");
00051 
00052   if ( arg.getNbOfArgs() < 1) {
00053     arg.usage();
00054     exit(1);
00055   }
00056 
00057 
00058 
00059   int Nsample = 3;
00060   int z_sample[3] = {14,24,45};
00061   float calib[3]  = {1,1,1};
00062 
00063   TString infilename = argv[1];
00064   TString outfilename = argv[2];
00065 
00066   int hitcutstart = atoi(argv[3]);
00067 
00068 
00069 
00070   float calib_layer[50];
00071   float layer_response = 1.0;
00072   float layer_response_rel_rms = 0.25;
00073 
00074 
00075 
00076   int Nlayer = 50;
00077 
00078   int nhit_calo_cut = 60;
00079 
00080   int nshower = 0;
00081   int nmuon = 0;
00082   int npenetrating_muon = 0;
00083 
00084   int nstart = 0;
00085 
00086   int nhit_cumul = 0;
00087 
00088   int zstart = -1;
00089 
00090   int nevt_tot = 0;
00091 
00092 
00093 
00094 
00095 
00096 
00097   TH1F * hdensity = new TH1F("hdensity","",100,0,50);
00098   TH2F * hnhit_density = new TH2F("hnhit_density","",1000,0,2000,200,0,100);
00099 
00100   TH1F * hrhadron = new TH1F("hrhadron","",1000,0,10);
00101   TH2F * hnhit_rhadron = new TH2F("hnhit_rhadron","",1000,0,2000,200,0,10);
00102 
00103   TProfile * hlong1 = new TProfile("hlong1","",50,0,50,0,150,"");
00104   TProfile * hlong2 = new TProfile("hlong2","",50,0,50,0,150,"");
00105   TProfile * hlong3 = new TProfile("hlong3","",50,0,50,0,150,"");
00106 
00107   TGraphErrors * tgmip = new TGraphErrors();
00108   TGraphErrors * tgshower = new TGraphErrors();
00109 
00110   tgmip->SetName("tgmip");
00111   tgshower->SetName("tgshower");
00112 
00113   TH1F * hlayer_mult = new TH1F("hlayer_mult","",50,0,50);
00114 
00115   TH1I * hnhit_calo = new TH1I("hnhit_calo","",500,0,5000);
00116   TH1I * hnhit_evt  = new TH1I("hnhit_evt","all showers",500,0,5000);
00117   TH1I * hnhit_evt2 = new TH1I("hnhit_evt2","all showers starting wthin 5 first planes",500,0,5000);
00118   TH1I * hnhit_evt3 = new TH1I("hnhit_evt3","all showers without MIP part",500,0,5000);
00119   TH1I * hnhit_evt4 = new TH1I("hnhit_evt4","all showers without MIP part starting wthin 5 first planes",500,0,5000);
00120   TH1I * hnhit_evt5  = new TH1I("hnhit_evt5","all showers with multiplicity per layer",500,0,5000);
00121   TH1I * hnhit_evt6 = new TH1I("hnhit_evt6","all showers without MIP part starting wthin 5 first planes with multiplicity per layer",500,0,5000);
00122 
00123   TH1I * hstart = new TH1I("hstart","",Nlayer,0,Nlayer);
00124 
00125   TH1I * hnhit_layer = new TH1I("hnhit_layer","",200,0,200);
00126 
00127   TH1I * hz_nhit = new TH1I("hz_nhit","",50,0,50);
00128   TH1I * hz_nhit_cumul = new TH1I("hz_nhit_cumul","",50,0,50);
00129 
00130   TH2I * hnhit_zstart = new TH2I("hnhit_zstart","Nhit for pions vs shower start",50,0,50,400,0,2000);
00131   TH2I * hnhit2_zstart = new TH2I("hnhit2_zstart","Nhit after removal of MIP hits for pions vs shower start",50,0,50,400,0,2000);
00132 
00133 
00134   TString name,title;
00135   TH2F * hxy[Nlayer];
00136   TH1F * hx[Nlayer];
00137   TH1F * hy[Nlayer];
00138   TH1I * hnhit_mip[Nlayer];
00139   TH1I * hnhit_shower[Nlayer];
00140   TH1I * hnhit_shower2[Nlayer];
00141   TH1I * hnhit_noise[Nlayer];
00142   TH1I * hdt[Nlayer];
00143   int nhit[Nlayer];
00144   int nnoise[Nlayer];
00145 
00146   for (int i=0;i<Nlayer;i++)
00147     {
00148       nhit[i] = 0;
00149       nnoise[i] = 0;
00150 
00151       name = Form("hnhit_mip_%i",i);
00152       title = Form("Number of CaloHit - chb %i",i);
00153       hnhit_mip[i]= new TH1I(name,title,100,0,100);
00154 
00155       name = Form("hnhit_shower_%i",i);
00156       title = Form("Number of CaloHit - chb %i",i);
00157       hnhit_shower[i]= new TH1I(name,title,1000,0,1000);
00158 
00159       name = Form("hnhit_shower2_%i",i);
00160       title = Form("Number of CaloHit - chb %i",i);
00161       hnhit_shower2[i]= new TH1I(name,title,1000,0,1000);
00162 
00163       name = Form("hnhit_noise_%i",i);
00164       title = Form("Number of NoiseHit - chb %i",i);
00165       hnhit_noise[i]= new TH1I(name,title,100,0,100);
00166 
00167       name = Form("hdt_%i",i);
00168       title = Form("Delta_t - chb %i",i);
00169       hdt[i]= new TH1I(name,title,21,-10,10);
00170 
00171       name = Form("hxy_%i",i);
00172       title = Form("CaloHit position - chb %i",i);
00173 
00174       if (i<48)
00175         {
00176           hxy[i] = new TH2F(name,title,96,0,96*1.04,96,0,96*1.04);
00177 
00178           name = Form("hx_%i",i);
00179           title = Form("X - chb %i",i);
00180           hx[i] = new TH1F(name,title,96,0,96*1.04);
00181 
00182           name = Form("hy_%i",i);
00183           title = Form("Y - chb %i",i);
00184           hy[i] = new TH1F(name,title,96,0,96*1.04);
00185         }
00186       else
00187         {
00188           hxy[i] = new TH2F(name,title,194,0,97,193,0,96.5);
00189 
00190           name = Form("hx_%i",i);
00191           title = Form("X - chb %i",i);
00192           hx[i] = new TH1F(name,title,194,0,97);
00193 
00194           name = Form("hy_%i",i);
00195           title = Form("Y - chb %i",i);
00196           hy[i] = new TH1F(name,title,193,0,96.5);
00197         }
00198     }
00199 
00200 
00201 
00202 
00203   for (int i=0 ; i<50 ; i++)
00204     {
00205       calib_layer[i] = layer_response * fabs(gRandom->Gaus(1,layer_response * layer_response_rel_rms));
00206       hlayer_mult->Fill(i,calib_layer[i]);
00207     }
00208   
00209 
00210 
00211 
00212 
00213 
00214 
00215   TFile f(infilename);
00216 
00217 
00218   TIter nextkey(f.GetListOfKeys());
00219   TKey *key;
00220 
00221   int NbTrees = 0;
00222 
00223   while (key = (TKey*)nextkey()) 
00224     {
00225       TString keyname = key->GetName();
00226       TTree *tree = (TTree*)key->ReadObj();                
00227       cout << "---------------------------- NEW TTREE -----------------"<< endl;
00228 
00229       int nEvts = tree->GetEntries();
00230 
00231       int nHits;
00232       int HitX[maxHits];
00233       int HitY[maxHits];
00234       int HitZ[maxHits];
00235       float HitEng[maxHits];
00236       tree->SetBranchAddress("nHits", &nHits);
00237       tree->SetBranchAddress("HitX", &HitX);
00238       tree->SetBranchAddress("HitY", &HitY);
00239       tree->SetBranchAddress("HitZ", &HitZ);
00240       tree->SetBranchAddress("HitEng", &HitEng);
00241 
00242 
00243       // Loop over events                                                                                                                                                                                                
00244 
00245       for (int evt = 0; evt < nEvts; evt++)
00246         {
00247           if (evt >= maxEvents) continue;
00248 
00249           tree->GetEntry(evt);
00250 
00251           int nhit_seperate1 = 0;
00252           int nhit_seperate2 = 0;
00253           float rhadron = 0;
00254 
00255           int nhit_evt = 0;
00256           int nhit_evt2 = 0;
00257           int nhit_evt3 = 0;
00258           int nhit_evt4 = 0;
00259           float nhit_evt5 = 0;
00260           float nhit_evt6 = 0;
00261 
00262           int nhotlayer = 0;
00263           float density = 0;
00264 
00265           int nhot_rear  = 0;
00266           int nhot_front = 0;
00267 
00268           for(int i=0 ; i < Nlayer ; i++)
00269             {
00270               nhit[i] = 0;
00271             }
00272                 
00273           zstart = -1;
00274 
00275 
00276           cout<<"evt "<<evt<<" -> nhit = "<<nHits<<" \r"<<flush;
00277 
00278           hnhit_calo->Fill(nHits);
00279 
00280           
00281 
00282 
00283           /********** read of hits ************/
00284                                                                                                                                                                             
00285           for (int hit = 0; hit < nHits; hit++)
00286             {
00287               int xpos = HitX[hit] + 48;
00288               int ypos = HitY[hit] + 48;
00289               int zSlot = HitZ[hit];
00290               float Edep = HitEng[hit];
00291 
00292               if (Edep !=0 )
00293                 {
00294                   hxy[zSlot]->Fill(xpos,ypos);
00295                         
00296                   nhit[zSlot]++;
00297                 }
00298             }   
00299 
00300 
00301 
00302 
00303           /********** sum of hits ************/
00304 
00305           for(int i=0 ; i < Nlayer ; i++)
00306             {
00307               nhit_evt += nhit[i];
00308               nhit_evt5 += (nhit[i] * calib_layer[i]);
00309               
00310               if (nhit[i] != 0)
00311                 {
00312                   density += nhit[i];
00313                   nhotlayer++;
00314 
00315                   if (i<5){nhit_seperate1++;}
00316                   if (i>=15){nhit_seperate2++;}
00317 
00318                   if (i<10){nhot_front++;}
00319                   if (i>=19 && i<30){nhot_rear++;}
00320                 }
00321             }
00322 
00323 
00324 
00325 
00326 
00327 
00328           /************ calculate density ****************/
00329 
00330           if (nhotlayer != 0)
00331             {
00332               density /= (1.*nhotlayer);
00333               hdensity->Fill(density);
00334             }
00335 
00336           if (nhit_seperate2 != 0)
00337             {
00338               rhadron = 1. * nhit_seperate1 / nhit_seperate2;
00339             }
00340 
00341 
00342 
00343 
00344           /******** select muons *******/
00345 
00346           if (density !=0 && density < 2)
00347             {
00348               /********** select penetrating ************/
00349 
00350               nmuon++;
00351 
00352               if (nhot_front >= 6 && nhot_rear >= 6)
00353                 {
00354                   
00355                   npenetrating_muon++;
00356                   
00357                   for(int i=0 ; i < Nlayer ; i++)
00358                     {
00359                       hnhit_mip[i]->Fill(nhit[i]);
00360                     }
00361                 }
00362             }
00363 
00364 
00365 
00366 
00367           /******** select showers *******/
00368 
00369           if (density >= 5 && nhit_evt >= nhit_calo_cut)
00370             {          
00371               
00372               hrhadron->Fill(rhadron);
00373               hnhit_rhadron->Fill(nhit_evt,rhadron);
00374 
00375               nshower++;
00376               
00377               /******** find shower start *********/
00378               
00379               for(int i=0 ; i < Nlayer-2 ; i++)
00380                 {
00381                   
00382                   /******** rising number of hits in 3 consecutive layer *********/
00383                   
00384                   if ( (nhit[i+2] > nhit[i+1]) && (nhit[i+1] > nhit[i]) )
00385                     {
00386                       for (int k=0 ; k<3 ; k++)
00387                         {                         
00388 
00389                           /******* shower begins when nhit >= 4 ***********/
00390 
00391                           if (nhit[i+k] >= hitcutstart)
00392                             {
00393                               nstart++;
00394                               zstart = i+k;
00395                               
00396                               hstart->Fill(zstart);
00397                               k = 3;
00398                               i = Nlayer;
00399                             }
00400                         }
00401                     }
00402                 }
00403               
00404               
00405               
00406               /********* fill profiles *********/
00407               
00408               if (zstart != -1)
00409                 {
00410 
00411                   hnhit_zstart->Fill(zstart,nHits);
00412                   
00413                   for(int i=0 ; i < Nlayer ; i++)
00414                     {
00415                       hnhit_shower[i]->Fill(nhit[i]);
00416                       
00417                       hz_nhit->Fill(i,nhit[i]);
00418                       hlong1->Fill(i,nhit[i]);
00419 
00420                       //all showers starting wthin 5 first planes
00421                       if (zstart < 5)
00422                         {
00423                           nhit_evt2 += nhit[i];
00424                         }
00425                       
00426                       //all showers without MIP part
00427                       if (i >= zstart && i <= 45)
00428                         {
00429                           nhit_evt3 += nhit[i];
00430 
00431                           //all showers without MIP part starting wthin 5 first planes
00432                           if (zstart < 5)
00433                             {
00434                               nhit_evt4 += nhit[i];
00435                               nhit_evt6 += (nhit[i] * calib_layer[i]);
00436                             }
00437                         }
00438                     }                                               
00439 
00440 
00441                   //all showers
00442                   hnhit_evt->Fill(nhit_evt);
00443                   hnhit_evt5->Fill(nhit_evt5);
00444 
00445                   //all showers starting wthin 5 first planes
00446                   if (nhit_evt2 != 0){hnhit_evt2->Fill(nhit_evt2);}
00447 
00448                   //all showers without MIP part
00449                   hnhit_evt3->Fill(nhit_evt3);
00450 
00451                   //all showers without MIP part starting wthin 5 first planes
00452                   if (nhit_evt4 != 0){hnhit_evt4->Fill(nhit_evt4);hnhit_evt6->Fill(nhit_evt6);}
00453 
00454                   
00455 
00456 
00457                   for (int k=0 ; k<Nsample ; k++)
00458                     {
00459                       if (zstart <= z_sample[k])
00460                         {
00461                           hnhit_shower2[z_sample[k] - zstart]->Fill(nhit[z_sample[k]]);
00462                           hlong2->Fill(z_sample[k] - zstart,nhit[z_sample[k]]);
00463                           hlong3->Fill(z_sample[k] - zstart,nhit[z_sample[k]] * calib[k]);
00464                         }
00465                     }
00466                 }
00467             }
00468         }
00469     }
00470 
00471 
00472 
00473   cout<<endl;
00474   cout<<endl;
00475   cout<<"Nevt_tot = "<<nevt_tot<<endl;
00476   cout<<endl;
00477   cout<<"Nshower = "<<nshower<<endl;
00478   cout<<"Nstart = "<<nstart<<endl;
00479   cout<<endl;
00480   cout<<"Nmuon = "<<nmuon<<endl;
00481   cout<<"Npenetrating_muon = "<<npenetrating_muon<<endl;
00482   cout<<endl;
00483 
00484 
00485 
00486 
00487 
00488   cout<<endl;
00489   cout<<"Write to TFile "<<outfilename<<endl;
00490   cout<<endl;
00491   
00492   TFile * tf = new TFile(outfilename,"RECREATE");
00493   
00494   for (int i=0;i<Nlayer;i++)
00495     {
00496       tgmip->SetPoint(i,i,hnhit_mip[i]->GetMean());
00497       tgmip->SetPointError(i,0,hnhit_mip[i]->GetRMS());
00498       
00499       tgshower->SetPoint(i,i,hnhit_shower[i]->GetMean());
00500       tgshower->SetPointError(i,0,hnhit_shower[i]->GetRMS());
00501       
00502       hxy[i]->Write();
00503       hdt[i]->Write();
00504       hnhit_mip[i]->Write();
00505       hnhit_shower[i]->Write();
00506       hnhit_shower2[i]->Write();
00507     }
00508   
00509   tgmip->SetMarkerStyle(20);
00510   tgmip->SetMarkerColor(2);
00511   
00512   tgshower->SetMarkerStyle(20);
00513   tgshower->SetMarkerColor(3);
00514   
00515   tgmip->Write();
00516   tgshower->Write();
00517 
00518   hlayer_mult->Write();
00519 
00520   hrhadron->Write();
00521   hnhit_rhadron->Write();
00522   hnhit_zstart->Write();
00523   hdensity->Write();
00524   hnhit_density->Write();
00525   hstart->Write();
00526   hlong1->Write();
00527   hlong2->Write();
00528   hlong3->Write();
00529   hnhit_evt->Write();
00530   hnhit_evt2->Write();
00531   hnhit_evt3->Write();
00532   hnhit_evt4->Write();
00533   hnhit_evt5->Write();
00534   hnhit_evt6->Write();
00535   hnhit_layer->Write();
00536   hz_nhit->Write();
00537   hz_nhit_cumul->Write();
00538   hnhit_calo->Write();
00539   
00540   tf->Close();
00541 
00542   
00543 }

Functions
int	main (int argc, char **argv)
Variables
const int	maxHits = 10000
const int	maxEvents = 10000
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Iro/profile4layers_simu.cpp File Reference

Functions

Variables

Function Documentation

Variable Documentation
