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