macro to get single muon trigger effiency versus pt (Fabien)
[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 /* Macro to produce trigger single muon efficiency versus pt plots for the 
20    3 pt cuts. Results are compared to the reference (red curves).   
21    To be used with (at least) 10000 events as follows
22    AliGenBox * gener = new AliGenBox(1);
23    gener->SetPtRange(0.,10.);
24    gener->SetPhiRange(0., 360.);         
25    gener->SetThetaRange(171.000,178.001);
26    gener->SetPart(13);           // or -13
27    gener->SetOrigin(0.,0., 0.);  
28    gener->SetSigma(0.0, 0.0, 0.0);     
29
30    Author: Fabien Guerin, LPC Clermont-Ferrand, Jan. 2006
31 */
32
33
34 // ROOT includes
35 #include "TBranch.h"
36 #include "TClonesArray.h"
37 #include "TFile.h"
38 #include "TF1.h"
39 #include "TH1.h"
40 #include "TH2.h"
41 #include "TParticle.h"
42 #include "TTree.h"
43 #include "TMath.h"
44 #include "TStyle.h"
45 #include "TCanvas.h"
46 #include "TLegend.h"
47
48
49 // STEER includes
50 #include "AliRun.h"
51 #include "AliDetector.h"
52 #include "AliRunLoader.h"
53 #include "AliHeader.h"
54 #include "AliLoader.h"
55 #include "AliStack.h"
56 #include "AliSegmentation.h"
57 #include "AliMC.h"
58
59 // MUON includes
60 #include "AliMUON.h"
61 #include "AliMUONData.h"
62 #include "AliMUONHit.h"
63 #include "AliMUONChamber.h"
64 #include "AliMUONConstants.h"
65 #include "AliMUONDigit.h"
66 #include "AliMUONRawCluster.h"
67 #include "AliMUONGlobalTrigger.h"
68 #include "AliMUONLocalTrigger.h"
69 #include "AliMUONTrack.h"
70
71 Double_t fitArc(Double_t *x,Double_t *par) 
72 {
73   Double_t h1=par[1]*(x[0]+par[2]);
74   return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
75 }
76
77 Double_t fitArch(Double_t *x,Double_t *par) 
78 {
79   Double_t h0=2*par[0]/(TMath::Pi());
80   Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0];
81   return h0*TMath::TanH(h1)+par[3];
82 }
83
84 void MUONTriggerEfficiencyPt(char filename[10]="galice.root")
85 {
86
87 // define style
88     TStyle *st1 = new TStyle("st1","My STYLE");
89     
90     st1->SetOptStat(0);
91     st1->SetOptFit(111);  
92     st1->SetFrameFillColor(10);
93     st1->SetCanvasColor(10);
94     st1->SetFillColor(10);
95     st1->SetTitleBorderSize(0);
96     st1->SetTitleOffset(1.2,"XY");
97     st1->SetTitleSize(0.05,"XY"); 
98     st1->SetLabelSize(0.045,"XY");  
99     st1->SetLabelFont(22,"XY");  
100     st1->SetTitleFont(22,"XY");     
101     st1->SetOptTitle(0);   
102     st1->SetStatColor(10);
103     st1->SetStatFont(62);     
104     st1->SetFitFormat("4.4g");
105     st1->SetPadLeftMargin(0.15);
106     st1->SetPadRightMargin(0.06); 
107     st1->SetPadTopMargin(0.06);
108     st1->SetPadBottomMargin(0.15); 
109     st1->cd();
110     
111     gROOT->ForceStyle();
112     //TGaxis::SetMaxDigits(3);
113     
114 // beginning of macro    
115     TH1F *ptapt    = new TH1F("ptapt","",50,0.15,5.); 
116     TH1F *ptlpt    = new TH1F("ptlpt","",50,0.15,5.); 
117     TH1F *pthpt    = new TH1F("pthpt","",50,0.15,5.);
118     TH1F *ptcoinc  = new TH1F("ptcoinc","",50,0.15,5.);  
119         
120     TParticle *particle;
121     AliStack* stack; 
122     
123     Double_t coincmuon=0;
124     Double_t aptmuon=0;
125     Double_t lptmuon=0;
126     Double_t hptmuon=0;
127     Double_t ptmu;
128
129 // output file
130     char digitdat[100];
131     sprintf(digitdat,"MUONTriggerEfficiencyPt.out");  
132     FILE *fdat=fopen(digitdat,"w");
133
134 // Creating Run Loader and openning file containing Hits
135     AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
136     
137     if (RunLoader ==0x0) {
138         printf(">>> Error : Error Opening %s file \n",filename);
139         return;
140     }
141     
142     Int_t nevents = RunLoader->GetNumberOfEvents();
143     
144 // loading hits
145     AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");     
146     MUONLoader->LoadHits("READ");
147     AliMUONData data_hits(MUONLoader,"MUON","MUON");
148     
149     TClonesArray * globalTrigger;
150     AliMUONGlobalTrigger * gloTrg;
151
152 // loading trigger
153     MUONLoader->LoadDigits("READ");  
154     AliMUONData muondata(MUONLoader,"MUON","MUON");
155
156 // Loading kine
157     RunLoader->LoadKinematics("READ");
158     
159     for (Int_t ievent=0; ievent<nevents; ievent++) {  // Event loop
160         
161         RunLoader->GetEvent(ievent);
162         
163         if (ievent%500==0) printf("ievent = %d \n",ievent);
164
165
166 // kine
167         Int_t iparticle, nparticles;
168         stack = RunLoader->Stack();
169         nparticles = (Int_t) stack->GetNtrack();        
170         for (iparticle=0; iparticle<nparticles; iparticle++) {             
171             particle = stack->Particle(iparticle);   
172             Float_t pt = particle->Pt();            
173             Int_t pdgcode = particle->GetPdgCode();                        
174             if (pdgcode==-13 || pdgcode==13) ptmu = pt;             
175         }
176
177 // trigger 
178         muondata.SetTreeAddress("D,GLT");
179         muondata.GetTriggerD();
180     
181         globalTrigger = muondata.GlobalTrigger();
182
183         Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1 
184         for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
185             gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
186             
187             if (gloTrg->SinglePlusApt()>=1 || gloTrg->SingleMinusApt()>=1 || gloTrg->SingleUndefApt()) {
188                 aptmuon++;
189                 ptapt->Fill(ptmu);
190             }
191             if (gloTrg->SinglePlusLpt()>=1 || gloTrg->SingleMinusLpt()>=1 || gloTrg->SingleUndefLpt()) {
192                 lptmuon++;
193                 ptlpt->Fill(ptmu);
194             }       
195             if (gloTrg->SinglePlusHpt()>=1 || gloTrg->SingleMinusHpt()>=1 || gloTrg->SingleUndefHpt()) {
196                 hptmuon++;
197                 pthpt->Fill(ptmu);
198             }
199         } // end of loop on Global Trigger      
200         muondata.ResetTrigger();  
201
202 // Hits
203         RunLoader->GetEvent(ievent);  
204         data_hits.SetTreeAddress("H");    
205         
206         Int_t itrack, ntracks, NbHits[4];
207         Int_t SumNbHits;
208     
209         for (Int_t j=0; j<4; j++) NbHits[j]=0;
210  
211         ntracks = (Int_t) data_hits.GetNtracks();      
212      
213         for (itrack=0; itrack<ntracks; itrack++) { // track loop
214          data_hits.GetTrack(itrack); 
215
216          Int_t ihit, nhits;
217          nhits = (Int_t) data_hits.Hits()->GetEntriesFast();   
218          AliMUONHit* mHit;
219     
220           for (ihit=0; ihit<nhits; ihit++) {
221             mHit = static_cast<AliMUONHit*>(data_hits.Hits()->At(ihit));
222             Int_t Nch        = mHit->Chamber(); 
223             Int_t hittrack   = mHit->Track();
224             Int_t IdPart     = mHit->Particle();
225                                         
226             for (Int_t j=0;j<4;j++) {
227               Int_t kch=11+j;
228               if (hittrack==0) {
229                 if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++; 
230               }               
231             }    
232           }
233           data_hits.ResetHits();
234         } // end track loop     
235
236         
237 // 3/4 coincidence 
238         SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
239         
240         if (SumNbHits==3 || SumNbHits==4) {
241             coincmuon++;
242             ptcoinc->Fill(ptmu);
243         }            
244                 
245       } // end loop on event  
246       
247       MUONLoader->UnloadHits();
248       MUONLoader->UnloadDigits();
249       MUONLoader->UnloadRecPoints();
250       RunLoader->UnloadKinematics();   
251
252       delete RunLoader;
253      
254       fprintf(fdat,"\n");    
255       fprintf(fdat,"\n");
256       fprintf(fdat," Number of events = %d \n",nevents);  
257       fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
258       fprintf(fdat," Number of dimuons with 3/4 coinc Apt cut = %d \n",(Int_t)aptmuon);
259       fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon); 
260       fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);  
261       fprintf(fdat,"\n"); 
262       
263       Double_t efficiency,error;  
264       
265       efficiency=aptmuon/coincmuon;
266       error=efficiency*TMath::Sqrt((aptmuon+coincmuon)/(aptmuon*coincmuon));
267       fprintf(fdat," Efficiency Apt cut = %4.4f +/- %4.4f\n",efficiency,error);
268       
269       efficiency=lptmuon/coincmuon;  
270       error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));  
271       fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
272       
273       efficiency=hptmuon/coincmuon; 
274       error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon)); 
275       fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
276       
277       fclose(fdat);
278       
279       Double_t x1,x2,xval,xerr;
280       
281       for (Int_t i=0;i<50;i++) {
282           x1=ptcoinc->GetBinContent(i+1);
283           
284           x2=ptapt->GetBinContent(i+1);
285           if (x1!=0 && x2!=0) {
286               xval=x2/x1;
287               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
288               ptapt->SetBinContent(i+1,xval);
289               ptapt->SetBinError(i+1,0);
290           } else {
291               ptapt->SetBinContent(i+1,0.);
292               ptapt->SetBinError(i+1,0.);     
293           }
294           
295           x2=ptlpt->GetBinContent(i+1);
296           if (x1!=0 && x2!=0) {
297               xval=x2/x1;
298               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
299               ptlpt->SetBinContent(i+1,xval);
300               ptlpt->SetBinError(i+1,0);
301           } else {
302               ptlpt->SetBinContent(i+1,0.);
303               ptlpt->SetBinError(i+1,0.);     
304           } 
305           
306           x2=pthpt->GetBinContent(i+1);
307           if (x1!=0 && x2!=0) {
308               xval=x2/x1;
309               xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
310               pthpt->SetBinContent(i+1,xval);
311               pthpt->SetBinError(i+1,0);
312           } else {
313               pthpt->SetBinContent(i+1,0.);
314               pthpt->SetBinError(i+1,0.);     
315           }              
316       }
317       
318       
319       TF1 *fitapt = new TF1("fitapt",fitArch,0.,5.,4);     
320       TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4); 
321       TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);      
322       
323       TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400); 
324       c1->Divide(3,1);      
325       
326       c1->cd(1);
327       ptapt->SetMinimum(0.);
328       ptapt->SetMaximum(1.05);   
329       ptapt->SetTitle("");
330       ptapt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); 
331       ptapt->GetYaxis()->SetTitle("Efficiency");     
332       //ptapt->GetXaxis()->SetRange(3);
333       ptapt->Draw("");
334       fitapt->SetParameters(1.853,1.57,-0.0203,-0.178);   
335       fitapt->SetLineColor(2);
336       fitapt->SetLineWidth(2);
337       fitapt->Draw("SAME");    
338       TLegend * leg = new TLegend(0.5,0.38,0.85,0.53); 
339       leg->SetBorderSize(0);
340       leg->SetFillColor(0);
341       leg->SetTextSize(0.05);
342       leg->SetTextFont(22);
343       leg->SetHeader("Apt trigger pt cut");
344       leg->AddEntry(fitapt,"Ref.","l");
345       leg->AddEntry(ptapt,"New","l");                 
346       leg->Draw("SAME");     
347       
348       c1->cd(2);
349       ptlpt->SetMinimum(0.);
350       ptlpt->SetMaximum(1.05);   
351       ptlpt->SetTitle("");
352       ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
353       ptlpt->GetYaxis()->SetTitle("Efficiency");     
354       //ptlpt->GetXaxis()->SetRange(3);
355       ptlpt->Draw(""); 
356       fitlpt->SetParameters(0.602,0.774,0.273,0.048);  
357       fitlpt->SetLineColor(2);
358       fitlpt->SetLineWidth(2);
359       fitlpt->Draw("SAME");    
360       leg = new TLegend(0.5,0.38,0.85,0.53); 
361       leg->SetBorderSize(0);
362       leg->SetFillColor(0);
363       leg->SetTextSize(0.05);
364       leg->SetTextFont(22);
365       leg->SetHeader("Lpt trigger pt cut");
366       leg->AddEntry(fitlpt,"Ref.","l");
367       leg->AddEntry(ptlpt,"New","l");                 
368       leg->Draw("SAME");
369       
370       c1->cd(3);
371       pthpt->SetMinimum(0.);
372       pthpt->SetMaximum(1.05);   
373       pthpt->SetTitle("");
374       pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");  
375       pthpt->GetYaxis()->SetTitle("Efficiency");     
376       //pthpt->GetXaxis()->SetRange(3);
377       pthpt->Draw(""); 
378       fithpt->SetParameters(0.627,0.393,0.855,0.0081); 
379       fithpt->SetLineColor(2);
380       fithpt->SetLineWidth(2);
381       fithpt->Draw("SAME");    
382       leg = new TLegend(0.5,0.38,0.85,0.53); 
383       leg->SetBorderSize(0);
384       leg->SetFillColor(0);
385       leg->SetTextSize(0.05);
386       leg->SetTextFont(22);
387       leg->SetHeader("Hpt trigger pt cut");
388       leg->AddEntry(fithpt,"Ref.","l");
389       leg->AddEntry(pthpt,"New","l");                 
390       leg->Draw("SAME");
391       
392       c1->SaveAs("MUONTriggerEfficiencyPt.gif");
393       c1->SaveAs("MUONTriggerEfficiencyPt.eps");
394       
395 }