/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/getChamberEfficiency.cpp File Reference

#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 "tools/Arguments.hh"

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 = 200000
const int maxTrees = 100
const int maxHits = 200


Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 46 of file getChamberEfficiency.cpp.

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


Variable Documentation

const int maxEvents = 200000

Definition at line 41 of file getChamberEfficiency.cpp.

const int maxTrees = 100

Definition at line 42 of file getChamberEfficiency.cpp.

Referenced by main().

const int maxHits = 200

Definition at line 43 of file getChamberEfficiency.cpp.


Generated on Mon Jan 7 13:17:15 2013 for MicromegasFramework by  doxygen 1.4.7