Fix histo name
[u/mrichter/AliRoot.git] / MUON / MUONTriggerEfficiency.C
index 6af15c3..2a8c1d6 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*  */
+// $Id$
 
-//
-// Macro for checking trigger integrated efficiency for dimuons.
-// The efficiency is calculated with respect to the 3/4 coincidence.
-// Author: Fabien Guerin, LPC Clermont-Ferrand, Jan. 2006
-//
+/// \ingroup macros
+/// \file MUONTriggerEfficiency.C
+/// \brief Macro for checking trigger integrated efficiency for dimuons.
+///
+/// The efficiency is calculated with respect to the 3/4 coincidence.
+///
+/// \author Fabien Guerin, LPC Clermont-Ferrand, Jan. 2006
 
 #if !defined(__CINT__) || defined(__MAKECINT__)
 // ROOT includes
 #include "AliHeader.h"
 #include "AliLoader.h"
 #include "AliStack.h"
+#include "AliCDBManager.h"
 
 // MUON includes
 #include "AliMUON.h"
-#include "AliMUONData.h"
 #include "AliMUONHit.h"
 #include "AliMUONConstants.h"
 #include "AliMUONDigit.h"
-#include "AliMUONRawCluster.h"
 #include "AliMUONGlobalTrigger.h"
-#include "AliMUONLocalTrigger.h"
 #include "AliMUONTrack.h"
+
+#include "AliMUONDataInterface.h"
+#include "AliMUONMCDataInterface.h"
+#include "AliMUONVHitStore.h"
+#include "AliMUONVTriggerStore.h"
+
 #endif
 
 // Upsilon(1S)
 
 
-void MUONTriggerEfficiency (char filename[10]="galice.root"){
-// output file
-  char digitdat[100];
-  char currentfile[100]; 
-  sprintf(digitdat,"MUONTriggerEfficiency.out"); 
+void MUONTriggerEfficiency(const char* filenameSim="galice_sim.root", 
+                           const char* filenameRec="galice.root",  
+                           Bool_t readFromRP = 0)
+{
+  
+  // Set default CDB storage
+  AliCDBManager* man = AliCDBManager::Instance();
+  man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+
+  // output file
+  
+  AliMUONMCDataInterface diSim(filenameSim);
+  AliMUONDataInterface diRec(filenameRec);
+  
+  if (!diSim.IsValid())
+  {
+    cerr << "Cannot access " << filenameSim << endl;
+    return;
+  }
+  
+  if (!diRec.IsValid())
+  {
+    cerr << "Cannot access " << filenameRec << endl;
+    return;
+  }
   
-  FILE *fdat=fopen(digitdat,"w");          
+  FILE* fdat = fopen("MUONTriggerEfficiency.out","w");          
+  if (!fdat)
+  {
+    cerr << "Cannot create output file" << endl;
+    return;
+  }
   
-  Int_t nevents,coincmuon,muonapt,muonlpt,muonhpt;
+  Int_t coincmuon,muonlpt,muonhpt;
   Int_t CoincMuPlus,CoincMuMinus;
   coincmuon=0;
-  muonapt=0;
   muonlpt=0;  
-  muonhpt=0;         
+  muonhpt=0;
   
-// Initialise AliRoot
-   // Creating Run Loader and openning file containing Hits
-   AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
-    
-   if (RunLoader ==0x0) {
-       printf(">>> Error : Error Opening %s file \n",currentfile);
-       return;
-   }
-             
-   nevents = RunLoader->GetNumberOfEvents();          
-     
-   //Loading hits  
-   AliLoader * MUONHits = RunLoader->GetLoader("MUONLoader");     
-   MUONHits->LoadHits("READ");
-   AliMUONData data_hits(MUONHits,"MUON","MUON");
-   
-    //Loading Digits  
-   AliLoader * MUONDigits = RunLoader->GetLoader("MUONLoader");    
-   MUONDigits->LoadDigits("READ");
-   AliMUONData data_digits(MUONDigits,"MUON","MUON");
+  Int_t nevents = diSim.NumberOfEvents();          
   
-   TClonesArray * globalTrigger;
-   AliMUONGlobalTrigger * gloTrg; 
-   for (Int_t ievent=0; ievent<nevents; ievent++) {    // event loop
-       CoincMuPlus=0;
-       CoincMuMinus=0;
-       
-       
-     RunLoader->GetEvent(ievent);
-         
-     if (ievent%1000==0) printf("\t Event = %d\n",ievent);    
-
-// Hits
-    data_hits.SetTreeAddress("H");    
-        
-    Int_t itrack, ntracks, NbHits[2][4];
-    Int_t SumNbHits;
-    for (Int_t j=0; j<2; j++) {
-     for (Int_t jj=0; jj<4; jj++) {
-      NbHits[j][jj]=0;
-     }
+  for (Int_t ievent=0; ievent<nevents; ++ievent) 
+  {   
+    CoincMuPlus=0;
+    CoincMuMinus=0;
+    
+    if (ievent%1000==0) printf("\t Event = %d\n",ievent);    
+    
+    // Hits
+    
+    Int_t NbHits[2][4];
+    for (Int_t j=0; j<2; j++) 
+    {
+      for (Int_t jj=0; jj<4; jj++) 
+      {
+        NbHits[j][jj]=0;
+      }
     } 
-    ntracks = (Int_t) data_hits.GetNtracks();      
-    for (itrack=0; itrack<ntracks; itrack++) { // Track loop
-      data_hits.GetTrack(itrack); 
-
-      Int_t ihit, nhits;
-      nhits = (Int_t) data_hits.Hits()->GetEntriesFast();   
+    
+    Int_t ntracks = (Int_t) diSim.NumberOfTracks(ievent);
+    
+    for ( Int_t itrack=0; itrack<ntracks; ++itrack ) 
+    {
+      AliMUONVHitStore* hitStore = diSim.HitStore(ievent,itrack);      
       AliMUONHit* mHit;
-      for(ihit=0; ihit<nhits; ihit++) {
-        mHit = static_cast<AliMUONHit*>(data_hits.Hits()->At(ihit));
-        Int_t Nch        = mHit->Chamber(); 
-        Int_t hittrack   = mHit->Track();
-        Float_t IdPart     = mHit->Particle();
-
-        for (Int_t j=0;j<4;j++) {
-         Int_t kch=11+j;
-         if (hittrack==1 || hittrack==2) {
-          if (Nch==kch && IdPart==-13 && NbHits[0][j]==0) NbHits[0][j]++; 
-          if (Nch==kch && IdPart==13 && NbHits[1][j]==0) NbHits[1][j]++; 
+      TIter next(hitStore->CreateIterator());
+      
+      while ( ( mHit = static_cast<AliMUONHit*>(next()) ) )
+      {
+        Int_t Nch = mHit->Chamber(); 
+        Int_t hittrack = mHit->Track();
+        Float_t IdPart = mHit->Particle();
+        
+        for (Int_t j=0;j<4;j++) 
+        {
+          Int_t kch=11+j;
+          if (hittrack==1 || hittrack==2) 
+          {
+            if (Nch==kch && IdPart==-13 && NbHits[0][j]==0) NbHits[0][j]++; 
+            if (Nch==kch && IdPart==13 && NbHits[1][j]==0) NbHits[1][j]++; 
           }       
         }
-       }
-       
-      data_hits.ResetHits();
-      
+      }
     } // end track loop     
     
-// 3/4 coincidence
-    SumNbHits=NbHits[0][0]+NbHits[0][1]+NbHits[0][2]+NbHits[0][3];
+    // 3/4 coincidence
+    Int_t SumNbHits=NbHits[0][0]+NbHits[0][1]+NbHits[0][2]+NbHits[0][3];
     if (SumNbHits==3 || SumNbHits==4) CoincMuPlus=1;
-       
+    
     SumNbHits=NbHits[1][0]+NbHits[1][1]+NbHits[1][2]+NbHits[1][3];
     if (SumNbHits==3 || SumNbHits==4) CoincMuMinus=1;
     
     if (CoincMuPlus==1 && CoincMuMinus==1) coincmuon++;         
-           
-// Trigger
-    data_digits.SetTreeAddress("D,GLT"); 
-    data_digits.GetTriggerD();
-   
-    globalTrigger = data_digits.GlobalTrigger();
-
-    Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
+          
+    TString tree("D");
+    if ( readFromRP ) tree = "R";
+    
+    AliMUONVTriggerStore* triggerStore = diRec.TriggerStore(ievent,tree.Data());
 
-    for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
-      gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
-      
-       if (gloTrg->PairUnlikeApt()>=1) muonapt++;                     
-       if (gloTrg->PairUnlikeLpt()>=1) muonlpt++;
-       if (gloTrg->PairUnlikeHpt()>=1) muonhpt++;
-                               
-    } // end of loop on Global Trigger         
+    AliMUONGlobalTrigger* gloTrg = triggerStore->Global();
     
-    data_digits.ResetTrigger();    
+    if (gloTrg->PairUnlikeLpt()>=1) muonlpt++;
+    if (gloTrg->PairUnlikeHpt()>=1) muonhpt++;
 
   } // end loop on event  
-
-  MUONHits->UnloadHits();
-  MUONDigits->UnloadDigits();  
-  
  
   // calculate efficiency with as a ref. at least 3/4 planes fired
-  Float_t efficiencyapt,efficiencylpt,efficiencyhpt;
-  Double_t coincmu,aptmu,lptmu,hptmu;  
+  Float_t efficiencylpt,efficiencyhpt;
+  Double_t coincmu,lptmu,hptmu;  
   Float_t error;  
-  coincmu=Double_t(coincmuon);  
-  aptmu=Double_t(muonapt);
+  coincmu=Double_t(coincmuon);
   lptmu=Double_t(muonlpt); 
   hptmu=Double_t(muonhpt); 
 
-
   // output
   fprintf(fdat,"\n"); 
   fprintf(fdat," Number of events = %d \n",nevents);  
   fprintf(fdat," Number of events with 3/4 coinc = %d \n",coincmuon);  
-  fprintf(fdat," Number of dimuons with 3/4 coinc Apt cut = %d \n",muonapt);  
   fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",muonlpt);  
   fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",muonhpt);  
   fprintf(fdat,"\n");
   
-  efficiencyapt=aptmu/coincmu;
-  error=efficiencyapt*TMath::Sqrt((aptmu+coincmu)/(aptmu*coincmu));
-  fprintf(fdat," Efficiency Apt cut = %4.4f +/- %4.4f\n",efficiencyapt,error);
-               
   efficiencylpt=lptmu/coincmu;  
   error=efficiencylpt*TMath::Sqrt((lptmu+coincmu)/(lptmu*coincmu));  
   fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiencylpt,error);