renaming function to avoid library conflict (discovered in test/generators/TUHKMgen)
[u/mrichter/AliRoot.git] / MUON / MUONTriggerEfficiencyPt.C
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18
19 /// \ingroup macros
20 /// \file MUONTriggerEfficiencyPt.C
21 /// \brief Macro to produce trigger single muon efficiency versus pt plots for the 
22 /// 2 pt cuts.
23 /// 
24 /// See full description on the \ref README_trigger page.
25 ///
26 /// \author Fabien Guerin (LPC)
27
28 // ROOT includes
29 #include "TBranch.h"
30 #include "TClonesArray.h"
31 #include "TFile.h"
32 #include "TF1.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "TParticle.h"
36 #include "TTree.h"
37 #include "TMath.h"
38 #include "TStyle.h"
39 #include "TCanvas.h"
40 #include "TLegend.h"
41 #include "Riostream.h"
42
43 // STEER includes
44 #include "AliRun.h"
45 #include "AliRunLoader.h"
46 #include "AliHeader.h"
47 #include "AliLoader.h"
48 #include "AliStack.h"
49 #include "AliCDBManager.h"
50
51 // MUON includes
52 #include "AliMUONDataInterface.h"
53 #include "AliMUONMCDataInterface.h"
54 #include "AliMUONVHitStore.h"
55 #include "AliMUONVTriggerStore.h"
56 #include "AliMUONHit.h"
57 #include "AliMUONDigit.h"
58 #include "AliMUONGlobalTrigger.h"
59 #include "AliMUONTrack.h"
60
61 Double_t fitArc(Double_t *x,Double_t *par) 
62 {
63   Double_t h1=par[1]*(x[0]+par[2]);
64   return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
65 }
66
67 Double_t fitArch(Double_t *x,Double_t *par) 
68 {
69   Double_t h0=2*par[0]/(TMath::Pi());
70   Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0];
71   return h0*TMath::TanH(h1)+par[3];
72 }
73
74 void MUONTriggerEfficiencyPt(const char* filenameSim="galice_sim.root", 
75                              const char* filenameRec="galice.root",  
76                              Bool_t readFromRP = 0)
77 {
78
79 // define style
80     TStyle *st1 = new TStyle("st1","My STYLE");
81     
82     st1->SetOptStat(0);
83     st1->SetOptFit(111);  
84     st1->SetFrameFillColor(10);
85     st1->SetCanvasColor(10);
86     st1->SetFillColor(10);
87     st1->SetTitleBorderSize(0);
88     st1->SetTitleOffset(1.2,"XY");
89     st1->SetTitleSize(0.05,"XY"); 
90     st1->SetLabelSize(0.045,"XY");  
91     st1->SetLabelFont(22,"XY");  
92     st1->SetTitleFont(22,"XY");     
93     st1->SetOptTitle(0);   
94     st1->SetStatColor(10);
95     st1->SetStatFont(62);     
96     st1->SetFitFormat("4.4g");
97     st1->SetPadLeftMargin(0.15);
98     st1->SetPadRightMargin(0.06); 
99     st1->SetPadTopMargin(0.06);
100     st1->SetPadBottomMargin(0.15); 
101     st1->cd();
102     
103 //    gROOT->ForceStyle();
104     //TGaxis::SetMaxDigits(3);
105     
106 // beginning of macro    
107     TH1F *ptlpt    = new TH1F("ptlpt","",50,0.15,5.); 
108     TH1F *pthpt    = new TH1F("pthpt","",50,0.15,5.);
109     TH1F *ptcoinc  = new TH1F("ptcoinc","",50,0.15,5.);  
110         
111     TParticle *particle;
112     
113     Int_t nevents;
114     Double_t coincmuon=0;
115     Double_t lptmuon=0;
116     Double_t hptmuon=0;
117     Double_t ptmu=0;
118
119     AliCDBManager* cdbManager = AliCDBManager::Instance();
120     cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
121     cdbManager->SetRun(0);
122
123     AliMUONMCDataInterface diSim(filenameSim);
124     AliMUONDataInterface diRec(filenameRec);
125     
126     if (!diSim.IsValid())
127     {
128         cerr << "Cannot access " << filenameSim << endl;
129         return;
130     }
131     
132     if (!diRec.IsValid())
133     {
134         cerr << "Cannot access " << filenameRec << endl;
135         return;
136     }
137     
138     FILE* fdat = fopen("MUONTriggerEfficiencyPt.out","w");          
139     if (!fdat)
140     {
141         cerr << "Cannot create output file" << endl;
142         return;
143     }
144
145     nevents = diSim.NumberOfEvents();
146
147     AliRunLoader* runLoader = AliRunLoader::Open(filenameSim);    
148
149     for (Int_t ievent=0; ievent<nevents; ievent++) {  // Event loop
150         
151         if (ievent%500==0) printf("ievent = %d \n",ievent);
152 // kine
153         
154         runLoader->GetEvent(ievent);
155         runLoader->LoadKinematics();
156         AliStack* stack = runLoader->Stack();
157         
158         Int_t nparticles = (Int_t) stack->GetNtrack();        
159
160         for (Int_t iparticle=0; iparticle<nparticles; iparticle++) {
161             particle = stack->Particle(iparticle);   
162             Float_t pt = particle->Pt();            
163             Int_t pdgcode = particle->GetPdgCode();                        
164             if (pdgcode==-13 || pdgcode==13) ptmu = pt;             
165         }
166
167 // global trigger
168         TString tree("D");
169         if ( readFromRP ) tree = "R";
170         
171         AliMUONVTriggerStore* triggerStore = diRec.TriggerStore(ievent,tree.Data());
172         AliMUONGlobalTrigger* gloTrg = triggerStore->Global();
173         
174         if (gloTrg->SingleLpt()>=1 ) {
175             lptmuon++;
176             ptlpt->Fill(ptmu);
177         }       
178         if (gloTrg->SingleHpt()>=1 ) {
179             hptmuon++;
180             pthpt->Fill(ptmu);
181         }
182
183 // Hits
184         Int_t itrack, ntracks, NbHits[4];
185         Int_t SumNbHits;
186         
187         for (Int_t j=0; j<4; j++) NbHits[j]=0;
188  
189         ntracks = (Int_t) diSim.NumberOfTracks(ievent);
190    
191         for (itrack=0; itrack<ntracks; itrack++) { // track loop
192             AliMUONVHitStore* hitStore = diSim.HitStore(ievent,itrack);      
193             AliMUONHit* mHit;
194             TIter next(hitStore->CreateIterator());
195
196             while ( ( mHit = static_cast<AliMUONHit*>(next()) ) )
197             {
198                 Int_t Nch        = mHit->Chamber(); 
199                 Int_t hittrack   = mHit->Track();
200                 Float_t IdPart   = mHit->Particle();        
201             
202                 for (Int_t j=0;j<4;j++) {
203                     Int_t kch=11+j;
204                     if (hittrack==0) {
205                         if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++; 
206                     }               
207                 }    
208             }
209         } // end track loop
210
211
212 // 3/4 coincidence 
213         SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
214
215         if (SumNbHits==3 || SumNbHits==4) {
216             coincmuon++;
217             ptcoinc->Fill(ptmu);
218         }            
219                 
220
221     } // end loop on event
222
223     if (coincmuon==0) {
224         cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
225         cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
226         cout << " >>> muon track per event as it should " << "\n";
227         cout << " >>> see README for further information " << "\n";
228         cout << " >>> exiting now ! " << "\n";
229         return;
230     }
231
232
233     fprintf(fdat,"\n");    
234     fprintf(fdat,"\n");
235     fprintf(fdat," Number of events = %d \n",nevents);  
236     fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
237     fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon); 
238     fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);  
239     fprintf(fdat,"\n"); 
240     
241     Double_t efficiency,error;
242     
243     efficiency=lptmuon/coincmuon;  
244     error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));  
245     fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
246     
247     efficiency=hptmuon/coincmuon; 
248     error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon)); 
249     fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
250     
251     fclose(fdat);
252     
253     Double_t x1,x2,xval,xerr;
254     
255     for (Int_t i=0;i<50;i++) {
256         x1=ptcoinc->GetBinContent(i+1);
257         
258         x2=ptlpt->GetBinContent(i+1);
259         if (x1!=0 && x2!=0) {
260             xval=x2/x1;
261             xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
262             ptlpt->SetBinContent(i+1,xval);
263             ptlpt->SetBinError(i+1,0);
264         } else {
265             ptlpt->SetBinContent(i+1,0.);
266             ptlpt->SetBinError(i+1,0.);     
267         } 
268         
269         x2=pthpt->GetBinContent(i+1);
270         if (x1!=0 && x2!=0) {
271             xval=x2/x1;
272             xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
273             pthpt->SetBinContent(i+1,xval);
274             pthpt->SetBinError(i+1,0);
275         } else {
276             pthpt->SetBinContent(i+1,0.);
277             pthpt->SetBinError(i+1,0.);     
278         }              
279     }      
280     
281     TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4); 
282     TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);      
283     
284     TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400); 
285     c1->Divide(2,1);      
286     
287     c1->cd(1);
288     ptlpt->SetMinimum(0.);
289     ptlpt->SetMaximum(1.05);   
290     ptlpt->SetTitle("");
291     ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
292     ptlpt->GetYaxis()->SetTitle("Efficiency");     
293     //ptlpt->GetXaxis()->SetRange(3);
294     ptlpt->Draw(""); 
295     fitlpt->SetParameters(0.602,0.774,0.273,0.048);  
296     fitlpt->SetLineColor(2);
297     fitlpt->SetLineWidth(2);
298     fitlpt->Draw("SAME");    
299     TLegend * leg = new TLegend(0.5,0.38,0.85,0.53); 
300     leg = new TLegend(0.5,0.38,0.85,0.53); 
301     leg->SetBorderSize(0);
302     leg->SetFillColor(0);
303     leg->SetTextSize(0.05);
304     leg->SetTextFont(22);
305     leg->SetHeader("Lpt trigger pt cut");
306     leg->AddEntry(fitlpt,"Ref.","l");
307     leg->AddEntry(ptlpt,"New","l");                 
308     leg->Draw("SAME");
309     
310     c1->cd(2);
311     pthpt->SetMinimum(0.);
312     pthpt->SetMaximum(1.05);   
313     pthpt->SetTitle("");
314     pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");  
315     pthpt->GetYaxis()->SetTitle("Efficiency");     
316     //pthpt->GetXaxis()->SetRange(3);
317     pthpt->Draw(""); 
318     fithpt->SetParameters(0.627,0.393,0.855,0.0081); 
319     fithpt->SetLineColor(2);
320     fithpt->SetLineWidth(2);
321     fithpt->Draw("SAME");    
322     leg = new TLegend(0.5,0.38,0.85,0.53); 
323     leg->SetBorderSize(0);
324     leg->SetFillColor(0);
325     leg->SetTextSize(0.05);
326     leg->SetTextFont(22);
327     leg->SetHeader("Hpt trigger pt cut");
328     leg->AddEntry(fithpt,"Ref.","l");
329     leg->AddEntry(pthpt,"New","l");                 
330     leg->Draw("SAME");
331     
332     c1->SaveAs("MUONTriggerEfficiencyPt.gif");
333     c1->SaveAs("MUONTriggerEfficiencyPt.eps");      
334
335 }