Modifications to work with Stores (Philippe C.)
authorpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jul 2007 14:59:54 +0000 (14:59 +0000)
committerpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jul 2007 14:59:54 +0000 (14:59 +0000)
MUON/MUONTriggerEfficiencyPt.C

index 74114e2..d7acf21 100644 (file)
 #include "AliStack.h"
 
 // MUON includes
-#include "AliMUONSimData.h"
-#include "AliMUONRecData.h"
+#include "AliMUONDataInterface.h"
+#include "AliMUONMCDataInterface.h"
+#include "AliMUONVHitStore.h"
+#include "AliMUONVTriggerStore.h"
 #include "AliMUONHit.h"
 #include "AliMUONDigit.h"
 #include "AliMUONGlobalTrigger.h"
-#include "AliMUONLocalTrigger.h"
 #include "AliMUONTrack.h"
 
 Double_t fitArc(Double_t *x,Double_t *par) 
@@ -65,8 +66,8 @@ Double_t fitArch(Double_t *x,Double_t *par)
   return h0*TMath::TanH(h1)+par[3];
 }
 
-void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",  
-                             char filenameRec[10]="galice.root",
+void MUONTriggerEfficiencyPt(const char* filenameSim="galice_sim.root", 
+                            const char* filenameRec="galice.root",  
                              Bool_t readFromRP = 0)
 {
 
@@ -111,130 +112,95 @@ void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",
     Double_t hptmuon=0;
     Double_t ptmu=0;
 
-// output file
-    char digitdat[100];
-    char currentfile[100];
-    sprintf(digitdat,"MUONTriggerEfficiencyPt.out");  
-    FILE *fdat=fopen(digitdat,"w");
+    AliMUONMCDataInterface diSim(filenameSim);
+    AliMUONDataInterface diRec(filenameRec);
     
-// Initialise AliRoot
-// Creating Run Loader and openning file containing Hits
-    AliRunLoader * RunLoaderSim = AliRunLoader::Open(filenameSim,"MUONFolderSim","READ");
-    
-    if (RunLoaderSim ==0x0) {
-       printf(">>> Error : Error Opening %s file \n",currentfile);
+    if (!diSim.IsValid())
+    {
+       cerr << "Cannot access " << filenameSim << endl;
        return;
     }
     
-    AliRunLoader * RunLoaderRec = AliRunLoader::Open(filenameRec,"MUONFolder","READ");
-    
-    if (RunLoaderRec ==0x0) {
-       printf(">>> Error : Error Opening %s file \n",currentfile);
+    if (!diRec.IsValid())
+    {
+       cerr << "Cannot access " << filenameRec << endl;
        return;
     }
     
-    nevents = RunLoaderSim->GetNumberOfEvents();
-    
-    AliLoader * MUONLoaderSim = RunLoaderSim->GetLoader("MUONLoader");
-    AliLoader * MUONLoaderRec = RunLoaderRec->GetLoader("MUONLoader");
-    
-    if (!readFromRP) {
-       cout << " reading from digits \n";
-       MUONLoaderSim->LoadDigits("READ");
-    } else {
-       cout << " reading from RecPoints \n";
-       MUONLoaderRec->LoadRecPoints("READ");
-   }
-    MUONLoaderSim->LoadHits("READ");
-    RunLoaderSim->LoadKinematics("READ");
-
+    FILE* fdat = fopen("MUONTriggerEfficiencyPt.out","w");          
+    if (!fdat)
+    {
+       cerr << "Cannot create output file" << endl;
+       return;
+    }
 
-// Creating MUON data containers
-    AliMUONSimData muondataSim(MUONLoaderSim,"MUON","MUON");
-    AliMUONRecData muondataRec(MUONLoaderRec,"MUON","MUON");
+    nevents = diSim.NumberOfEvents();
 
-    TClonesArray * globalTrigger;
-    AliMUONGlobalTrigger * gloTrg;
+    AliRunLoader* runLoader = AliRunLoader::Open(filenameSim);    
 
     for (Int_t ievent=0; ievent<nevents; ievent++) {  // Event loop
        
-       RunLoaderSim->GetEvent(ievent);
-       RunLoaderRec->GetEvent(ievent);
-       
        if (ievent%500==0) printf("ievent = %d \n",ievent);
-
 // kine
-       Int_t iparticle, nparticles;
-       stack = RunLoaderSim->Stack();
-       nparticles = (Int_t) stack->GetNtrack();        
-       for (iparticle=0; iparticle<nparticles; iparticle++) {             
+       
+       runLoader->GetRunLoader()->GetEvent(ievent);
+       runLoader->GetRunLoader()->LoadKinematics();
+       AliStack* stack = runLoader->GetRunLoader()->Stack();
+       
+       Int_t nparticles = (Int_t) stack->GetNtrack();        
+
+        for (Int_t iparticle=0; iparticle<nparticles; iparticle++) {
             particle = stack->Particle(iparticle);   
             Float_t pt = particle->Pt();           
             Int_t pdgcode = particle->GetPdgCode();                        
            if (pdgcode==-13 || pdgcode==13) ptmu = pt;             
         }
 
-// trigger 
-       if (!readFromRP) {
-           muondataSim.SetTreeAddress("D,GLT"); 
-           muondataSim.GetTriggerD();
-           globalTrigger = muondataSim.GlobalTrigger();
-       } else {    
-           muondataRec.SetTreeAddress("RC,TC"); 
-           muondataRec.GetTrigger();
-           globalTrigger = muondataRec.GlobalTrigger();
-
+// global trigger
+       TString tree("D");
+       if ( readFromRP ) tree = "R";
+       
+       AliMUONVTriggerStore* triggerStore = diRec.TriggerStore(ievent,tree.Data());
+       AliMUONGlobalTrigger* gloTrg = triggerStore->Global();
+       
+       if (gloTrg->SingleLpt()>=1 ) {
+           lptmuon++;
+           ptlpt->Fill(ptmu);
+       }       
+       if (gloTrg->SingleHpt()>=1 ) {
+           hptmuon++;
+           pthpt->Fill(ptmu);
        }
 
-        Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1 
-        for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
-           gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
-           
-           if (gloTrg->SingleLpt()>=1 ) {
-               lptmuon++;
-               ptlpt->Fill(ptmu);
-           }       
-           if (gloTrg->SingleHpt()>=1 ) {
-               hptmuon++;
-               pthpt->Fill(ptmu);
-           }
-        } // end of loop on Global Trigger      
-        muondataSim.ResetTrigger();  
-        muondataRec.ResetTrigger();  
-
 // Hits
-       muondataSim.SetTreeAddress("H");    
-
         Int_t itrack, ntracks, NbHits[4];
         Int_t SumNbHits;
-
+       
         for (Int_t j=0; j<4; j++) NbHits[j]=0;
  
-       ntracks = (Int_t) muondataSim.GetNtracks();
+       ntracks = (Int_t) diSim.NumberOfTracks(ievent);
    
         for (itrack=0; itrack<ntracks; itrack++) { // track loop
-           muondataSim.GetTrack(itrack);
+           AliMUONVHitStore* hitStore = diSim.HitStore(ievent,itrack);      
+           AliMUONHit* mHit;
+           TIter next(hitStore->CreateIterator());
 
-         Int_t ihit, nhits;
-        nhits = (Int_t) muondataSim.Hits()->GetEntriesFast();
-         AliMUONHit* mHit;
-    
-          for (ihit=0; ihit<nhits; ihit++) {
-             mHit = static_cast<AliMUONHit*>(muondataSim.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==0) {
-                if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++; 
-               }               
-            }    
-          }
-         muondataSim.ResetHits();
+           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==0) {
+                       if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++; 
+                   }               
+               }    
+           }
         } // end track loop
 
+
 // 3/4 coincidence 
         SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
 
@@ -243,8 +209,9 @@ void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",
             ptcoinc->Fill(ptmu);
         }            
                 
+
     } // end loop on event
-    
+
     if (coincmuon==0) {
        cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
        cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
@@ -254,114 +221,107 @@ void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",
        return;
     }
 
-      
-       MUONLoaderSim->UnloadHits();
-       if (!readFromRP) {
-          MUONLoaderSim->UnloadDigits();  
-       } else {    
-          MUONLoaderRec->UnloadRecPoints();
-       }
-       RunLoaderSim->UnloadKinematics();
-     
-      fprintf(fdat,"\n");    
-      fprintf(fdat,"\n");
-      fprintf(fdat," Number of events = %d \n",nevents);  
-      fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
-      fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon); 
-      fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);  
-      fprintf(fdat,"\n"); 
-      
-      Double_t efficiency,error;
 
-      efficiency=lptmuon/coincmuon;  
-      error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));  
-      fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
-      
-      efficiency=hptmuon/coincmuon; 
-      error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon)); 
-      fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
-      
-      fclose(fdat);
-      
-      Double_t x1,x2,xval,xerr;
-      
-      for (Int_t i=0;i<50;i++) {
-         x1=ptcoinc->GetBinContent(i+1);
-         
-         x2=ptlpt->GetBinContent(i+1);
-         if (x1!=0 && x2!=0) {
-             xval=x2/x1;
-             xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
-             ptlpt->SetBinContent(i+1,xval);
-             ptlpt->SetBinError(i+1,0);
-         } else {
-             ptlpt->SetBinContent(i+1,0.);
-             ptlpt->SetBinError(i+1,0.);     
-         } 
-         
-         x2=pthpt->GetBinContent(i+1);
-         if (x1!=0 && x2!=0) {
-             xval=x2/x1;
-             xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
-             pthpt->SetBinContent(i+1,xval);
-             pthpt->SetBinError(i+1,0);
-         } else {
-             pthpt->SetBinContent(i+1,0.);
-             pthpt->SetBinError(i+1,0.);     
-         }              
-      }      
-      
-      TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4); 
-      TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);      
-      
-      TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400); 
-      c1->Divide(2,1);      
-      
-      c1->cd(1);
-      ptlpt->SetMinimum(0.);
-      ptlpt->SetMaximum(1.05);   
-      ptlpt->SetTitle("");
-      ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
-      ptlpt->GetYaxis()->SetTitle("Efficiency");     
-      //ptlpt->GetXaxis()->SetRange(3);
-      ptlpt->Draw(""); 
-      fitlpt->SetParameters(0.602,0.774,0.273,0.048);  
-      fitlpt->SetLineColor(2);
-      fitlpt->SetLineWidth(2);
-      fitlpt->Draw("SAME");    
-      TLegend * leg = new TLegend(0.5,0.38,0.85,0.53); 
-      leg = new TLegend(0.5,0.38,0.85,0.53); 
-      leg->SetBorderSize(0);
-      leg->SetFillColor(0);
-      leg->SetTextSize(0.05);
-      leg->SetTextFont(22);
-      leg->SetHeader("Lpt trigger pt cut");
-      leg->AddEntry(fitlpt,"Ref.","l");
-      leg->AddEntry(ptlpt,"New","l");                 
-      leg->Draw("SAME");
-      
-      c1->cd(2);
-      pthpt->SetMinimum(0.);
-      pthpt->SetMaximum(1.05);   
-      pthpt->SetTitle("");
-      pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");  
-      pthpt->GetYaxis()->SetTitle("Efficiency");     
-      //pthpt->GetXaxis()->SetRange(3);
-      pthpt->Draw(""); 
-      fithpt->SetParameters(0.627,0.393,0.855,0.0081); 
-      fithpt->SetLineColor(2);
-      fithpt->SetLineWidth(2);
-      fithpt->Draw("SAME");    
-      leg = new TLegend(0.5,0.38,0.85,0.53); 
-      leg->SetBorderSize(0);
-      leg->SetFillColor(0);
-      leg->SetTextSize(0.05);
-      leg->SetTextFont(22);
-      leg->SetHeader("Hpt trigger pt cut");
-      leg->AddEntry(fithpt,"Ref.","l");
-      leg->AddEntry(pthpt,"New","l");                 
-      leg->Draw("SAME");
-      
-      c1->SaveAs("MUONTriggerEfficiencyPt.gif");
-      c1->SaveAs("MUONTriggerEfficiencyPt.eps");      
+    fprintf(fdat,"\n");    
+    fprintf(fdat,"\n");
+    fprintf(fdat," Number of events = %d \n",nevents);  
+    fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
+    fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon); 
+    fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);  
+    fprintf(fdat,"\n"); 
+    
+    Double_t efficiency,error;
+    
+    efficiency=lptmuon/coincmuon;  
+    error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));  
+    fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
+    
+    efficiency=hptmuon/coincmuon; 
+    error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon)); 
+    fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
+    
+    fclose(fdat);
+    
+    Double_t x1,x2,xval,xerr;
+    
+    for (Int_t i=0;i<50;i++) {
+       x1=ptcoinc->GetBinContent(i+1);
+       
+       x2=ptlpt->GetBinContent(i+1);
+       if (x1!=0 && x2!=0) {
+           xval=x2/x1;
+           xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
+           ptlpt->SetBinContent(i+1,xval);
+           ptlpt->SetBinError(i+1,0);
+       } else {
+           ptlpt->SetBinContent(i+1,0.);
+           ptlpt->SetBinError(i+1,0.);     
+       } 
+       
+       x2=pthpt->GetBinContent(i+1);
+       if (x1!=0 && x2!=0) {
+           xval=x2/x1;
+           xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
+           pthpt->SetBinContent(i+1,xval);
+           pthpt->SetBinError(i+1,0);
+       } else {
+           pthpt->SetBinContent(i+1,0.);
+           pthpt->SetBinError(i+1,0.);     
+       }              
+    }      
+    
+    TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4); 
+    TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);      
+    
+    TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400); 
+    c1->Divide(2,1);      
+    
+    c1->cd(1);
+    ptlpt->SetMinimum(0.);
+    ptlpt->SetMaximum(1.05);   
+    ptlpt->SetTitle("");
+    ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
+    ptlpt->GetYaxis()->SetTitle("Efficiency");     
+    //ptlpt->GetXaxis()->SetRange(3);
+    ptlpt->Draw(""); 
+    fitlpt->SetParameters(0.602,0.774,0.273,0.048);  
+    fitlpt->SetLineColor(2);
+    fitlpt->SetLineWidth(2);
+    fitlpt->Draw("SAME");    
+    TLegend * leg = new TLegend(0.5,0.38,0.85,0.53); 
+    leg = new TLegend(0.5,0.38,0.85,0.53); 
+    leg->SetBorderSize(0);
+    leg->SetFillColor(0);
+    leg->SetTextSize(0.05);
+    leg->SetTextFont(22);
+    leg->SetHeader("Lpt trigger pt cut");
+    leg->AddEntry(fitlpt,"Ref.","l");
+    leg->AddEntry(ptlpt,"New","l");                 
+    leg->Draw("SAME");
+    
+    c1->cd(2);
+    pthpt->SetMinimum(0.);
+    pthpt->SetMaximum(1.05);   
+    pthpt->SetTitle("");
+    pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");  
+    pthpt->GetYaxis()->SetTitle("Efficiency");     
+    //pthpt->GetXaxis()->SetRange(3);
+    pthpt->Draw(""); 
+    fithpt->SetParameters(0.627,0.393,0.855,0.0081); 
+    fithpt->SetLineColor(2);
+    fithpt->SetLineWidth(2);
+    fithpt->Draw("SAME");    
+    leg = new TLegend(0.5,0.38,0.85,0.53); 
+    leg->SetBorderSize(0);
+    leg->SetFillColor(0);
+    leg->SetTextSize(0.05);
+    leg->SetTextFont(22);
+    leg->SetHeader("Hpt trigger pt cut");
+    leg->AddEntry(fithpt,"Ref.","l");
+    leg->AddEntry(pthpt,"New","l");                 
+    leg->Draw("SAME");
+    
+    c1->SaveAs("MUONTriggerEfficiencyPt.gif");
+    c1->SaveAs("MUONTriggerEfficiencyPt.eps");      
+
 }