ifstream member should not be persistent. Thanks to Christian for spotting it
[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"
b19dcb40 37#include "Riostream.h"
5707a463 38
39// STEER includes
40#include "AliRun.h"
5707a463 41#include "AliRunLoader.h"
42#include "AliHeader.h"
43#include "AliLoader.h"
44#include "AliStack.h"
5707a463 45
46// MUON includes
b19dcb40 47#include "AliMUONSimData.h"
48#include "AliMUONRecData.h"
5707a463 49#include "AliMUONHit.h"
5707a463 50#include "AliMUONDigit.h"
5707a463 51#include "AliMUONGlobalTrigger.h"
52#include "AliMUONLocalTrigger.h"
53#include "AliMUONTrack.h"
54
55Double_t fitArc(Double_t *x,Double_t *par)
56{
57 Double_t h1=par[1]*(x[0]+par[2]);
58 return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
59}
60
61Double_t fitArch(Double_t *x,Double_t *par)
62{
63 Double_t h0=2*par[0]/(TMath::Pi());
64 Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0];
65 return h0*TMath::TanH(h1)+par[3];
66}
67
b19dcb40 68void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",
69 char filenameRec[10]="galice.root",
70 Bool_t readFromRP = 0)
5707a463 71{
72
73// define style
74 TStyle *st1 = new TStyle("st1","My STYLE");
75
76 st1->SetOptStat(0);
77 st1->SetOptFit(111);
78 st1->SetFrameFillColor(10);
79 st1->SetCanvasColor(10);
80 st1->SetFillColor(10);
81 st1->SetTitleBorderSize(0);
82 st1->SetTitleOffset(1.2,"XY");
83 st1->SetTitleSize(0.05,"XY");
84 st1->SetLabelSize(0.045,"XY");
85 st1->SetLabelFont(22,"XY");
86 st1->SetTitleFont(22,"XY");
87 st1->SetOptTitle(0);
88 st1->SetStatColor(10);
89 st1->SetStatFont(62);
90 st1->SetFitFormat("4.4g");
91 st1->SetPadLeftMargin(0.15);
92 st1->SetPadRightMargin(0.06);
93 st1->SetPadTopMargin(0.06);
94 st1->SetPadBottomMargin(0.15);
95 st1->cd();
96
07d83718 97// gROOT->ForceStyle();
5707a463 98 //TGaxis::SetMaxDigits(3);
99
100// beginning of macro
5707a463 101 TH1F *ptlpt = new TH1F("ptlpt","",50,0.15,5.);
102 TH1F *pthpt = new TH1F("pthpt","",50,0.15,5.);
103 TH1F *ptcoinc = new TH1F("ptcoinc","",50,0.15,5.);
104
105 TParticle *particle;
106 AliStack* stack;
107
b19dcb40 108 Int_t nevents;
5707a463 109 Double_t coincmuon=0;
5707a463 110 Double_t lptmuon=0;
111 Double_t hptmuon=0;
522b4724 112 Double_t ptmu=0;
5707a463 113
114// output file
115 char digitdat[100];
b19dcb40 116 char currentfile[100];
5707a463 117 sprintf(digitdat,"MUONTriggerEfficiencyPt.out");
118 FILE *fdat=fopen(digitdat,"w");
b19dcb40 119
120// Initialise AliRoot
5707a463 121// Creating Run Loader and openning file containing Hits
b19dcb40 122 AliRunLoader * RunLoaderSim = AliRunLoader::Open(filenameSim,"MUONFolderSim","READ");
123
124 if (RunLoaderSim ==0x0) {
125 printf(">>> Error : Error Opening %s file \n",currentfile);
126 return;
127 }
5707a463 128
b19dcb40 129 AliRunLoader * RunLoaderRec = AliRunLoader::Open(filenameRec,"MUONFolder","READ");
130
131 if (RunLoaderRec ==0x0) {
132 printf(">>> Error : Error Opening %s file \n",currentfile);
5707a463 133 return;
134 }
135
b19dcb40 136 nevents = RunLoaderSim->GetNumberOfEvents();
5707a463 137
b19dcb40 138 AliLoader * MUONLoaderSim = RunLoaderSim->GetLoader("MUONLoader");
139 AliLoader * MUONLoaderRec = RunLoaderRec->GetLoader("MUONLoader");
5707a463 140
b19dcb40 141 if (!readFromRP) {
142 cout << " reading from digits \n";
143 MUONLoaderSim->LoadDigits("READ");
144 } else {
145 cout << " reading from RecPoints \n";
146 MUONLoaderRec->LoadRecPoints("READ");
147 }
148 MUONLoaderSim->LoadHits("READ");
149 RunLoaderSim->LoadKinematics("READ");
150
151
152// Creating MUON data containers
153 AliMUONSimData muondataSim(MUONLoaderSim,"MUON","MUON");
154 AliMUONRecData muondataRec(MUONLoaderRec,"MUON","MUON");
155
5707a463 156 TClonesArray * globalTrigger;
157 AliMUONGlobalTrigger * gloTrg;
158
5707a463 159 for (Int_t ievent=0; ievent<nevents; ievent++) { // Event loop
160
b19dcb40 161 RunLoaderSim->GetEvent(ievent);
162 RunLoaderRec->GetEvent(ievent);
5707a463 163
164 if (ievent%500==0) printf("ievent = %d \n",ievent);
165
5707a463 166// kine
167 Int_t iparticle, nparticles;
b19dcb40 168 stack = RunLoaderSim->Stack();
5707a463 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
07d83718 178 if (!readFromRP) {
b19dcb40 179 muondataSim.SetTreeAddress("D,GLT");
180 muondataSim.GetTriggerD();
181 globalTrigger = muondataSim.GlobalTrigger();
07d83718 182 } else {
b19dcb40 183 muondataRec.SetTreeAddress("RC,TC");
184 muondataRec.GetTrigger();
185 globalTrigger = muondataRec.GlobalTrigger();
186
07d83718 187 }
5707a463 188
189 Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
190 for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
191 gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
192
522b4724 193 if (gloTrg->SingleLpt()>=1 ) {
5707a463 194 lptmuon++;
195 ptlpt->Fill(ptmu);
196 }
522b4724 197 if (gloTrg->SingleHpt()>=1 ) {
5707a463 198 hptmuon++;
199 pthpt->Fill(ptmu);
200 }
201 } // end of loop on Global Trigger
b19dcb40 202 muondataSim.ResetTrigger();
203 muondataRec.ResetTrigger();
5707a463 204
205// Hits
b19dcb40 206 muondataSim.SetTreeAddress("H");
207
5707a463 208 Int_t itrack, ntracks, NbHits[4];
209 Int_t SumNbHits;
522b4724 210
5707a463 211 for (Int_t j=0; j<4; j++) NbHits[j]=0;
212
b19dcb40 213 ntracks = (Int_t) muondataSim.GetNtracks();
214
5707a463 215 for (itrack=0; itrack<ntracks; itrack++) { // track loop
b19dcb40 216 muondataSim.GetTrack(itrack);
5707a463 217
218 Int_t ihit, nhits;
b19dcb40 219 nhits = (Int_t) muondataSim.Hits()->GetEntriesFast();
5707a463 220 AliMUONHit* mHit;
221
222 for (ihit=0; ihit<nhits; ihit++) {
b19dcb40 223 mHit = static_cast<AliMUONHit*>(muondataSim.Hits()->At(ihit));
5707a463 224 Int_t Nch = mHit->Chamber();
225 Int_t hittrack = mHit->Track();
09f9a8f1 226 Float_t IdPart = mHit->Particle();
522b4724 227
5707a463 228 for (Int_t j=0;j<4;j++) {
229 Int_t kch=11+j;
522b4724 230 if (hittrack==0) {
5707a463 231 if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++;
522b4724 232 }
5707a463 233 }
234 }
b19dcb40 235 muondataSim.ResetHits();
522b4724 236 } // end track loop
07d83718 237
5707a463 238// 3/4 coincidence
239 SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
522b4724 240
5707a463 241 if (SumNbHits==3 || SumNbHits==4) {
242 coincmuon++;
243 ptcoinc->Fill(ptmu);
244 }
245
b19dcb40 246 } // end loop on event
247
248 if (coincmuon==0) {
249 cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
250 cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
251 cout << " >>> muon track per event as it should " << "\n";
252 cout << " >>> see README for further information " << "\n";
253 cout << " >>> exiting now ! " << "\n";
254 return;
255 }
5707a463 256
b19dcb40 257
258 MUONLoaderSim->UnloadHits();
259 if (!readFromRP) {
260 MUONLoaderSim->UnloadDigits();
261 } else {
262 MUONLoaderRec->UnloadRecPoints();
263 }
264 RunLoaderSim->UnloadKinematics();
5707a463 265
266 fprintf(fdat,"\n");
267 fprintf(fdat,"\n");
268 fprintf(fdat," Number of events = %d \n",nevents);
269 fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
5707a463 270 fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon);
271 fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);
272 fprintf(fdat,"\n");
273
522b4724 274 Double_t efficiency,error;
275
5707a463 276 efficiency=lptmuon/coincmuon;
277 error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));
278 fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
279
280 efficiency=hptmuon/coincmuon;
281 error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon));
282 fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
283
284 fclose(fdat);
285
286 Double_t x1,x2,xval,xerr;
287
288 for (Int_t i=0;i<50;i++) {
289 x1=ptcoinc->GetBinContent(i+1);
290
5707a463 291 x2=ptlpt->GetBinContent(i+1);
292 if (x1!=0 && x2!=0) {
293 xval=x2/x1;
294 xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
295 ptlpt->SetBinContent(i+1,xval);
296 ptlpt->SetBinError(i+1,0);
297 } else {
298 ptlpt->SetBinContent(i+1,0.);
299 ptlpt->SetBinError(i+1,0.);
300 }
301
302 x2=pthpt->GetBinContent(i+1);
303 if (x1!=0 && x2!=0) {
304 xval=x2/x1;
305 xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
306 pthpt->SetBinContent(i+1,xval);
307 pthpt->SetBinError(i+1,0);
308 } else {
309 pthpt->SetBinContent(i+1,0.);
310 pthpt->SetBinError(i+1,0.);
311 }
b19dcb40 312 }
5707a463 313
5707a463 314 TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4);
315 TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);
316
317 TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400);
522b4724 318 c1->Divide(2,1);
5707a463 319
320 c1->cd(1);
5707a463 321 ptlpt->SetMinimum(0.);
322 ptlpt->SetMaximum(1.05);
323 ptlpt->SetTitle("");
324 ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
325 ptlpt->GetYaxis()->SetTitle("Efficiency");
326 //ptlpt->GetXaxis()->SetRange(3);
327 ptlpt->Draw("");
328 fitlpt->SetParameters(0.602,0.774,0.273,0.048);
329 fitlpt->SetLineColor(2);
330 fitlpt->SetLineWidth(2);
331 fitlpt->Draw("SAME");
522b4724 332 TLegend * leg = new TLegend(0.5,0.38,0.85,0.53);
5707a463 333 leg = new TLegend(0.5,0.38,0.85,0.53);
334 leg->SetBorderSize(0);
335 leg->SetFillColor(0);
336 leg->SetTextSize(0.05);
337 leg->SetTextFont(22);
338 leg->SetHeader("Lpt trigger pt cut");
339 leg->AddEntry(fitlpt,"Ref.","l");
340 leg->AddEntry(ptlpt,"New","l");
341 leg->Draw("SAME");
342
522b4724 343 c1->cd(2);
5707a463 344 pthpt->SetMinimum(0.);
345 pthpt->SetMaximum(1.05);
346 pthpt->SetTitle("");
347 pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
348 pthpt->GetYaxis()->SetTitle("Efficiency");
349 //pthpt->GetXaxis()->SetRange(3);
350 pthpt->Draw("");
351 fithpt->SetParameters(0.627,0.393,0.855,0.0081);
352 fithpt->SetLineColor(2);
353 fithpt->SetLineWidth(2);
354 fithpt->Draw("SAME");
355 leg = new TLegend(0.5,0.38,0.85,0.53);
356 leg->SetBorderSize(0);
357 leg->SetFillColor(0);
358 leg->SetTextSize(0.05);
359 leg->SetTextFont(22);
360 leg->SetHeader("Hpt trigger pt cut");
361 leg->AddEntry(fithpt,"Ref.","l");
362 leg->AddEntry(pthpt,"New","l");
363 leg->Draw("SAME");
364
365 c1->SaveAs("MUONTriggerEfficiencyPt.gif");
07d83718 366 c1->SaveAs("MUONTriggerEfficiencyPt.eps");
5707a463 367}