#include <map>#include "TKey.h"#include <TH2F.h>#include <TProfile.h>#include <TROOT.h>#include <TRint.h>#include <TList.h>#include <TTree.h>#include <TFile.h>#include <TChain.h>#include "TChainElement.h"#include <TMath.h>#include <stdlib.h>#include <TSystem.h>#include <iostream>#include <sstream>#include "Log.hh"#include "MicroException.hh"#include <string>#include "mysql/Mysql.hh"#include "root/CaloRun.hh"#include "root/CaloEvent.hh"#include "root/CaloHit.hh"#include "root/MTRun.hh"#include "root/MTEvent.hh"#include "root/MTChannel.hh"#include "TGraphErrors.h"#include "TGraphAsymmErrors.h"#include <TF1.h>Include dependency graph for getChamberEfficiency.cpp:
Go to the source code of this file.
Functions | |
| int | main (int argc, char **argv) |
Variables | |
| const int | maxEvents = 3000000 |
| const int | maxTrees = 100 |
| const int | maxHits = 200 |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 45 of file getChamberEfficiency.cpp.
00045 { 00046 00047 if ( argc !=2 ) { 00048 FILE_LOG(logERROR) << "usage: example rootFile " << endl; 00049 exit(1); 00050 } 00051 00052 int search = 1; 00053 00054 TH2F * hxy = new TH2F("hxy","",100,0,100,100,0,100); 00055 TH1I * hnhit = new TH1I("hnhit","Number of all hits",200,0,1000); 00056 TH1F * hDensity = new TH1F("hDensity","Hits Density ",200,0,50.); 00057 TH1I * hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",130,0,130); 00058 TH1I * hnHitForw = new TH1I("hnhitForw","MIP Number of hits Forw",50,0,50); 00059 TProfile * tp_nhit = new TProfile("tp_nhit","",60,0,60,0,200,""); 00060 TH1F * hxdir = new TH1F("hxdir"," x Angle between hit and all the next ", 200, -3., 3.); 00061 TH1F * hydir = new TH1F("hydir"," y Angle between hit and all the next ", 200, -3., 3.); 00062 TProfile * hLayerMult = new TProfile("hLayerMult","Average Multiplicity around track per layer Y",50,0,50, 0.,1000,""); 00063 TH1F * hChi2X = new TH1F("hChi2X"," Chi2 X ", 20, 0., 20.); 00064 TH1F * hChi2Y = new TH1F("hChi2Y"," Chi2 Y ", 20, 0., 20.); 00065 TH1F * hResidY = new TH1F("hResidY"," Residual Y ", 60, -30., 30.); 00066 TH1F * hResidX = new TH1F("hResidX"," Residual X ", 60, -30., 30.); 00067 TProfile * pResidX = new TProfile("pResidX","Average Residual per layer X",50,0,50,-1000.,1000,""); 00068 TProfile * pResidY = new TProfile("pResidY","Average Residual per layer Y",50,0,50,-1000.,1000,""); 00069 TH1F * hPenteY = new TH1F("hPentY"," Slope Y ", 200, -100.,100.); 00070 TH1F * hPenteX = new TH1F("hPentX"," Slope X ", 200, -100, 100.); 00071 TH1F * hOffsetX = new TH1F("hOffsetX"," Offset X ", 200, 0, 1000.); 00072 TH1F * hOffsetY = new TH1F("hOffsetY"," Offset Y ", 200, 0, 1000.); 00073 00074 TH2F * hxyLayer0 = new TH2F("hxyLayer0","X verus Y fitted track. Layer 0",100,0,1000,100,0,1000); 00075 TH1I * hz = new TH1I("hz","",100,0,100); 00076 TH1I * hGlNbHtsFound = new TH1I("hGlNbHtsFound"," Nb of Hits found in Chamber", 50, 0, 50); 00077 TH1I * hGlNbHtsExtrap = new TH1I("hGlNbHtsExtrap"," Nb of Hits extrapolated in Chamber", 50, 0, 50); 00078 TH1I * hNbHitsFit = new TH1I("hNbHitsFit","Nb of Hits Used in Fit", 100, 0, 200 ); 00079 00080 // Book Layer Hit maps and efficiencies 00081 TH2I * hxyMap[50]; 00082 TH2I * hTrackxyMap[50]; 00083 TH2F * hEffMap[50]; 00084 TH2I * hRegxyMap[50] ; 00085 TH2I * hxyEvtMap[50] ; 00086 00087 TString name,title; 00088 int NbinsEffx=12; 00089 int NbinsEffy=12; 00090 int hbin=96; 00091 for (int i=0;i<50;i++) 00092 { 00093 name = Form("hxyMap_%i",i); 00094 title = Form("Hit distribution in Layer %i",i); 00095 hxyMap[i]= new TH2I(name,title,96,0,hbin,96,0,hbin); 00096 name = Form("hTrackxyMap_%i",i); 00097 hTrackxyMap[i]= new TH2I(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin); 00098 name = Form("hEffMap_%i",i); 00099 hEffMap[i]= new TH2F(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin); 00100 name = Form("hRegxyMap_%i",i); 00101 hRegxyMap[i]= new TH2I(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin); 00102 } 00103 00104 string rootName; 00105 rootName.assign(argv[1]); 00106 int nbHit = 0; 00107 int nbEvt = 0; 00108 00109 int NbHitLayer[100]; 00110 int IsPenetratingMIP; 00111 // float Hits[2][maxHits]; 00112 Double_t Zlayer[maxHits]; // Before it was Float 00113 Double_t Xlayer[maxHits]; 00114 Double_t Ylayer[maxHits]; 00115 Double_t ErPos[maxHits]; 00116 TString OutFilename = "Muons"+rootName+".root"; 00117 TFile tf(OutFilename,"RECREATE"); 00118 int NbTrees = 0; 00119 int RunId=1; 00120 00121 TChain chain("T"); 00122 for (int i =0; i<5 ; i++){ 00123 // if(i==2)continue; 00124 TString SubFile = Form("_I%i",i); 00125 TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root"); 00126 chain.Add(fname); 00127 cout << fname << endl; 00128 } 00129 00130 TObjArray *fileElements=chain.GetListOfFiles(); 00131 TIter next(fileElements); 00132 TChainElement *chEl=0; 00133 int NchainFiles=0; 00134 while (( chEl=(TChainElement*)next() )) 00135 { 00136 TFile f(chEl->GetTitle()); 00137 // TFile f(rootName.c_str()); 00138 TIter nextkey(f.GetListOfKeys()); 00139 TKey *key; 00140 NchainFiles++; 00141 cout << " Reading a new file in the chain " << NchainFiles << endl; 00142 while (key = (TKey*)nextkey()) 00143 { 00144 00145 TString keyname = key->GetName(); 00146 if ((keyname == "hResidualX") || (keyname == "hResidualY")){} 00147 else 00148 { 00149 00150 TTree *tree = (TTree*)key->ReadObj(); 00151 cout << "---------------------------- NEW TTREE -----------------"<< endl; 00152 nbHit = 0; 00153 NbTrees++; 00154 if( NbTrees > maxTrees ) continue; 00155 CaloRun *run=NULL; 00156 00157 try { 00158 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First()); 00159 if( run != NULL ){cout << "Run[" << run->GetRunId() << "], temperature[" << run->GetTemperature() << "] pression[" << run->GetPressure() << "] " << endl ; } 00160 00161 } 00162 catch ( ... ) {} 00163 CaloEvent *evt = new CaloEvent(); 00164 TBranch *branch= tree->GetBranch("CaloEvent"); 00165 branch->SetAddress(&evt); 00166 00167 CaloHit* hit =NULL; 00168 int nbEntries = tree->GetEntries(); 00169 cout << " Nb of events to read " << nbEntries << endl; 00170 TString name = Form("deltaTRun_%i",RunId); 00171 TString title = Form("delta distribution Run %i",RunId); 00172 TH1I * deltaTRun = new TH1I(name,title,20,-10,10); 00173 00174 for ( int evtNum = 0; evtNum < nbEntries ; evtNum++) 00175 { 00176 tree->GetEntry(evtNum); 00177 nbEvt++; 00178 // if(nbEvt%2000==0){cout << nbEvt << endl;} 00179 if(nbEvt>maxEvents) continue; 00180 /* 00181 cout << "**** Event Information "<< dec << evtNum+1 << "/" << nbEntries<< " ****" << endl; 00182 cout << " event num in TTree:" << evtNum <<endl; 00183 cout << " event id:" << evt->fEventId <<endl; 00184 00185 cout << " number of channel hit:" << evt->GetNHits() <<endl; 00186 */ 00187 int nbHits = evt->GetNHits(); 00188 // cout << "nb hits:" << nbHits << endl; 00189 00190 hnhit->Fill(nbHits); 00191 00192 if (nbHits > 120) continue; 00193 00194 //number of hits in region of interest 00195 int nhit_roi = 0; 00196 int nhit_fit = 0; 00197 float xpos = 0; 00198 float ypos = 0; 00199 float zpos = 0; 00200 unsigned int ncol = 0; 00201 unsigned int nrow = 0; 00202 int zSlot = 0; 00203 int deltat=0; 00204 int maxDeltat = evt->GetDeltaTmax(); 00205 00206 for (int layer=0 ; layer <100 ; layer++){ 00207 NbHitLayer[layer]=0;} 00208 for (int i=0 ; i<maxHits ;i++){ 00209 Zlayer[i]=0; 00210 Xlayer[i]=0; 00211 Ylayer[i]=0; 00212 double theta = (0.0136/100.)*sqrt((i+1)*2/1.757)*(1+0.038*log((i+1)*2/1.757)); 00213 ErPos[i] = sqrt(pow(10./sqrt(12.),2)+pow(2*28.*theta,2)); 00214 } 00215 00216 for (int i=0; i<50; i++) // Fill Maps for this event 00217 { 00218 name = Form("hxyEvtMap_%i",i); 00219 title = Form("Hit distribution in Layer %i",i); 00220 hxyEvtMap[i]= new TH2I(name,title,96,0,96,96,0,96); 00221 } 00222 00223 for(int i=0 ; i < nbHits ; i++) 00224 { 00225 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00226 deltat=hit->GetDeltaT(); 00227 deltaTRun->Fill(maxDeltat-deltat); 00228 00229 xpos = hit->GetX()/1000 ; // unit = mm 00230 ypos = hit->GetY()/1000 ; 00231 zSlot = hit->GetLayer() ; 00232 zpos = hit->GetZ()/1000; 00233 ncol = hit->GetColInChamber(); 00234 nrow = hit->GetRowInChamber(); 00235 00236 if(zSlot <48 && abs(deltat-maxDeltat)>2) continue; // This is for RPCs 00237 if(zSlot >=48 && abs(deltat-maxDeltat)>5) continue; // For Micromegas 00238 00239 //Eliminate Bad zones 00240 // if(zSlot==33 && xpos>64 ) continue; 00241 // if(zSlot==49 && (xpos>=32 &&xpos<40) && ( ypos>=48 && ypos <56) ) continue; 00242 00243 hz->Fill(zSlot); 00244 00245 NbHitLayer[zSlot]++; 00246 00247 hxy->Fill(ypos,xpos); 00248 00249 nhit_roi++; // Nb of good hits 00250 hxyEvtMap[zSlot]->Fill(ncol,nrow,1); 00251 hxyMap[zSlot]->Fill(ncol,nrow,1); 00252 00253 // Prepare for fitting a muon 00254 // Choose layers to make the fit 00255 if(nhit_fit==0){ 00256 Zlayer[nhit_fit] = zpos; // Convert coordinates to mm 00257 Xlayer[nhit_fit] = xpos; 00258 Ylayer[nhit_fit] = ypos; 00259 nhit_fit++; 00260 } 00261 else{ 00262 if (nhit_fit < maxHits){ 00263 if(NbHitLayer[zSlot] >1 ) continue ; // keep only 1st hit in layers 00264 // Select hits that have almost the same slope as the previous one 00265 Double_t slopy = TMath::ATan2( (ypos-Ylayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) ); 00266 Double_t slopx = TMath::ATan2( (xpos-Xlayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) ); 00267 hxdir->Fill(slopx); hydir->Fill(slopy); 00268 // hSlopyY->Fill(slopy); hSlopx->Fill(slopx); 00269 if ( TMath::Abs(slopx)> 1. || TMath::Abs(slopy)>1. ) continue; 00270 Zlayer[nhit_fit] = zpos; 00271 Xlayer[nhit_fit] = xpos; 00272 Ylayer[nhit_fit] = ypos; 00273 nhit_fit++; 00274 } 00275 } 00276 } // end hit 00277 00278 hnhit_roi->Fill(nhit_roi); 00279 00280 // Compute Hit Density to find MIP events 00281 int NbLayers=0; 00282 for (int layer=0 ; layer < 50 ; layer++){ 00283 if(NbHitLayer[layer]>0) 00284 { 00285 NbLayers++; 00286 tp_nhit->Fill(layer,NbHitLayer[layer]); 00287 } 00288 } 00289 00290 float HitDensity = 1. * nhit_roi/NbLayers; 00291 00292 hDensity->Fill(HitDensity); 00293 00294 //Look for penetrating MIPS 00295 00296 int NbLayersForw=0; 00297 int NbLayersBack=0; 00298 00299 IsPenetratingMIP = 0; 00300 if(1 < HitDensity < 3){ 00301 00302 for (int layer=0 ; layer < 10 ; layer++){ 00303 if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersForw++;} 00304 } 00305 for (int layer=38 ; layer < 48 ; layer++){ 00306 if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersBack++;} 00307 } 00308 00309 hnHitForw->Fill(NbLayersForw); 00310 hnHitForw->Fill(NbLayersBack+20); 00311 00312 //Define penetrating MIP 00313 if(NbLayersForw>6 && NbLayersBack > 6) { IsPenetratingMIP = 1;} 00314 } // end Penetrating MIPS 00315 00316 00317 // Fit a straight line for penetrating MIPs 00318 if(IsPenetratingMIP==1 && (nhit_fit > 25) ){ 00319 hNbHitsFit->Fill(nhit_fit); // How many hits are used for the fit 00320 float zmax = 1372.; 00321 TGraphErrors *grx = new TGraphErrors(nhit_fit,Zlayer,Xlayer,0,ErPos); 00322 TF1 *f1 = new TF1("f1","pol1",0,zmax); 00323 grx->Fit("f1","Q"); 00324 double betaX = f1->GetParameter(0); 00325 double alfaX = f1->GetParameter(1); 00326 double chi2X = f1->GetChisquare(); 00327 delete f1; delete grx; 00328 TGraphErrors *gry = new TGraphErrors(nhit_fit,Zlayer,Ylayer,0,ErPos); 00329 TF1 *f2 = new TF1("f2","pol1",0,zmax); 00330 gry->Fit("f2","Q"); 00331 double betaY = f2->GetParameter(0); 00332 double alfaY = f2->GetParameter(1); 00333 double chi2Y = f2->GetChisquare(); 00334 // if(nbEvt<100) gry->Write(); 00335 delete f2; delete gry; 00336 hPenteX->Fill(alfaX*1000.); 00337 hPenteY->Fill(alfaY*1000.); 00338 hOffsetX->Fill(betaX); 00339 hOffsetY->Fill(betaY); 00340 hChi2X->Fill(chi2X/(nhit_fit-2)); 00341 hChi2Y->Fill(chi2Y/(nhit_fit-2)); 00342 hxyLayer0->Fill(betaX,betaY); // Make X Y Map at the first layer 00343 00344 // Compute residuals to fitted track 00345 for (int j=1; j<nhit_fit ; j++){ 00346 float ResX = Xlayer[j]-( alfaX*Zlayer[j]+betaX ); 00347 float ResY = Ylayer[j]-( alfaY*Zlayer[j]+betaY ); 00348 hResidY->Fill(ResY); 00349 hResidX->Fill(ResX); 00350 pResidX->Fill(Zlayer[j]/28,ResX); 00351 pResidY->Fill(Zlayer[j]/28,ResY); 00352 } 00353 00354 // Go to compute Efficiency with good tracks 00355 if( chi2X/(nhit_fit-2) < 6. && chi2Y/(nhit_fit-2) < 6 ) { 00356 for(int i=0; i<50 ; i++){ 00357 int NpadsLook = 3; 00358 float Xextrapol = alfaX*i*28. + betaX ; // X extrapolate to chamber i in mm 00359 float Yextrapol = alfaY*i*28. + betaY ; // Y extrapolate to chamber i in mm 00360 00361 // Find Central pad from computed x,y position 00362 Int_t xpadCentr = -1000; 00363 Int_t ypadCentr = -1000; 00364 if( i>47 && ( fabs(Yextrapol-322.5) < 20 || fabs(Yextrapol-650) < 20 ) ) continue; 00365 if( i>47 && ( fabs(Xextrapol-485.0) < 20 ) ) continue; 00366 if(i<48){ 00367 xpadCentr = (int)(Xextrapol/10.4); 00368 ypadCentr = (int)(Yextrapol/10.4); 00369 } 00370 else{ 00371 if(Xextrapol<=480) {xpadCentr = (int)(Xextrapol/10.);} 00372 else {xpadCentr = (int)( (Xextrapol-5)/10 );} 00373 if(Yextrapol<=320.) {ypadCentr = (int)( (Yextrapol/10.) );} 00374 if(Yextrapol>325 && Yextrapol<=645.) {ypadCentr = (int)( (Yextrapol-5.)/10.);} 00375 if(Yextrapol>=646) {ypadCentr = (int)((Yextrapol-10.)/10.);} 00376 } 00377 // check that we are inside the chamber 00378 if(xpadCentr < 0 || xpadCentr > 95 ) continue; 00379 if(ypadCentr < 0 || ypadCentr > 95 ) continue; 00380 // define fiducial volume 00381 if(xpadCentr <= NpadsLook || xpadCentr>95-NpadsLook) continue ; 00382 if(ypadCentr <= NpadsLook || ypadCentr>95-NpadsLook) continue ; 00383 00384 hTrackxyMap[i]->Fill(ypadCentr,xpadCentr,1); //Fill extrapolated point 00385 hGlNbHtsExtrap->Fill(i); 00386 //Check now if there is a hit around this point in a region of NpadsLook 00387 int xbinL = max(xpadCentr - NpadsLook,0); // For the PS data we must increase the search zone in x 00388 int xbinH = min(xpadCentr + NpadsLook,95); 00389 int ybinL = max(ypadCentr - NpadsLook,0); 00390 int ybinH = min(ypadCentr + NpadsLook,95); 00391 int NhitsRegSearch = hxyEvtMap[i]->Integral(ybinL,ybinH,xbinL,xbinH); 00392 if( NhitsRegSearch >0 ) { 00393 hRegxyMap[i]->Fill(ypadCentr,xpadCentr,1); //Fill extrapolated point 00394 hGlNbHtsFound->Fill(i) ; 00395 hLayerMult->Fill(i, NhitsRegSearch) ; } 00396 00397 } // end of efficiency coputing end layer loop 00398 } // end of Chi2 selection 00399 } // end of penetrating MIP 00400 00401 if (nbEvt%10000 == 0) { 00402 cout << "Processing event: " << nbEvt 00403 << endl; 00404 } 00405 for (int i=0; i<50 ; i++){ 00406 delete hxyEvtMap[i]; 00407 } 00408 tf.cd(0); 00409 } // end event 00410 cout << " Finished with events in this TREE " << endl; 00411 deltaTRun->Write(); 00412 delete deltaTRun; 00413 } // else 00414 } // end of tree or run 00415 cout << " Ending file in chain " << endl; 00416 } // end of chain 00417 cout << " Going to finsh the work after" << nbEvt << " events "<< endl; 00418 tf.cd(); 00419 hPenteY->Write(); 00420 hPenteX->Write(); 00421 hOffsetX->Write(); 00422 hOffsetY->Write(); 00423 hxdir->Write(); 00424 hydir->Write(); 00425 tp_nhit->Write(); 00426 hLayerMult->Write(); 00427 hxy->Write(); 00428 hxyMap[33]->Write(); 00429 hxyMap[31]->Write(); 00430 hxyMap[48]->Write(); 00431 hxyMap[49]->Write(); 00432 hz->Write(); 00433 hnhit->Write(); 00434 hnhit_roi->Write(); 00435 hDensity->Write(); 00436 hnHitForw->Write(); 00437 hNbHitsFit->Write(); 00438 hChi2X->Write(); 00439 hChi2Y->Write(); 00440 hResidY->Write(); 00441 hResidX->Write(); 00442 pResidX->Write(); 00443 pResidY->Write(); 00444 hxyLayer0->Write(); 00445 for (int i=0 ; i<50 ; i++){ 00446 hRegxyMap[i]->Write(); 00447 hTrackxyMap[i]->Write(); 00448 hEffMap[i]->Divide(hRegxyMap[i], hTrackxyMap[i], 100., 1.,"b"); // Binomial errors 00449 hEffMap[i]->Write(); 00450 } 00451 TH1F * hEffRPC = new TH1F("hEffRPC"," RPC Efficiencies for every chip ", 100, 0., 100.); 00452 TH1F * hEffmMegas = new TH1F("hEffmMegas"," mMegas Efficiencies for every chip ", 100, 0., 100.); 00453 for (int i = 0; i< 50; i++){ 00454 for (int j = 1; j<13 ; j++){ 00455 for (int k=1; k<13; k++) { 00456 float eff = hEffMap[i]->GetBinContent(j,k); 00457 if(i<48) { hEffRPC->Fill(eff); } 00458 if(i>=48){ hEffmMegas->Fill(eff); } 00459 }// end columns 00460 } //end rows 00461 }// end Chamber 00462 00463 hEffRPC->Write(); 00464 hEffmMegas->Write(); 00465 00466 // Compute Global efficiency 00467 TH1F * hGlEffChamber = new TH1F("hGlEffChamber"," Chamber Global Efficiency ", 50, 0, 50 ); 00468 hGlEffChamber->Divide(hGlNbHtsFound, hGlNbHtsExtrap, 100. , 1., "b"); // Binomial errors 00469 hGlNbHtsFound->Write(); 00470 hGlNbHtsExtrap->Write(); 00471 hGlEffChamber->Write(); 00472 for (int j=1; j<=50 ; j++){ 00473 double nbTrFound= hGlNbHtsFound->GetBinContent(j); 00474 double nbTrExtr=hGlNbHtsExtrap->GetBinContent(j); 00475 double ef=hGlEffChamber->GetBinContent(j); 00476 double erEff=hGlEffChamber->GetBinError(j); 00477 double w = nbTrFound/nbTrExtr; 00478 double binErr=w*(1-w)/nbTrExtr; 00479 cout << "chamber " << j << " Nb tracks found and extrapol " << nbTrFound << " " << nbTrExtr << " Efficiency " << ef << " +- " << binErr << endl; 00480 } 00481 TGraphAsymmErrors * gr = new TGraphAsymmErrors(hGlNbHtsFound, hGlNbHtsExtrap,""); 00482 gr->Write(); 00483 }
| const int maxEvents = 3000000 |
| const int maxTrees = 100 |
| const int maxHits = 200 |
1.4.7