Fix coverity defect
[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
e54bf126 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)
5707a463 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"
b19dcb40 41#include "Riostream.h"
5707a463 42
43// STEER includes
44#include "AliRun.h"
5707a463 45#include "AliRunLoader.h"
46#include "AliHeader.h"
47#include "AliLoader.h"
48#include "AliStack.h"
8ee6f794 49#include "AliCDBManager.h"
5707a463 50
51// MUON includes
04b9b6ba 52#include "AliMUONDataInterface.h"
53#include "AliMUONMCDataInterface.h"
54#include "AliMUONVHitStore.h"
55#include "AliMUONVTriggerStore.h"
5707a463 56#include "AliMUONHit.h"
5707a463 57#include "AliMUONDigit.h"
5707a463 58#include "AliMUONGlobalTrigger.h"
5707a463 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
04b9b6ba 74void MUONTriggerEfficiencyPt(const char* filenameSim="galice_sim.root",
75 const char* filenameRec="galice.root",
b19dcb40 76 Bool_t readFromRP = 0)
5707a463 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
07d83718 103// gROOT->ForceStyle();
5707a463 104 //TGaxis::SetMaxDigits(3);
105
106// beginning of macro
5707a463 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;
5707a463 112
b19dcb40 113 Int_t nevents;
5707a463 114 Double_t coincmuon=0;
5707a463 115 Double_t lptmuon=0;
116 Double_t hptmuon=0;
522b4724 117 Double_t ptmu=0;
5707a463 118
8ee6f794 119 AliCDBManager* cdbManager = AliCDBManager::Instance();
120 cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
121 cdbManager->SetRun(0);
122
04b9b6ba 123 AliMUONMCDataInterface diSim(filenameSim);
124 AliMUONDataInterface diRec(filenameRec);
b19dcb40 125
04b9b6ba 126 if (!diSim.IsValid())
127 {
128 cerr << "Cannot access " << filenameSim << endl;
b19dcb40 129 return;
130 }
5707a463 131
04b9b6ba 132 if (!diRec.IsValid())
133 {
134 cerr << "Cannot access " << filenameRec << endl;
5707a463 135 return;
136 }
137
04b9b6ba 138 FILE* fdat = fopen("MUONTriggerEfficiencyPt.out","w");
139 if (!fdat)
140 {
141 cerr << "Cannot create output file" << endl;
142 return;
143 }
b19dcb40 144
04b9b6ba 145 nevents = diSim.NumberOfEvents();
b19dcb40 146
04b9b6ba 147 AliRunLoader* runLoader = AliRunLoader::Open(filenameSim);
5707a463 148
5707a463 149 for (Int_t ievent=0; ievent<nevents; ievent++) { // Event loop
150
5707a463 151 if (ievent%500==0) printf("ievent = %d \n",ievent);
5707a463 152// kine
04b9b6ba 153
8ee6f794 154 runLoader->GetEvent(ievent);
155 runLoader->LoadKinematics();
156 AliStack* stack = runLoader->Stack();
04b9b6ba 157
158 Int_t nparticles = (Int_t) stack->GetNtrack();
159
160 for (Int_t iparticle=0; iparticle<nparticles; iparticle++) {
5707a463 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
04b9b6ba 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);
07d83718 181 }
5707a463 182
5707a463 183// Hits
5707a463 184 Int_t itrack, ntracks, NbHits[4];
185 Int_t SumNbHits;
04b9b6ba 186
5707a463 187 for (Int_t j=0; j<4; j++) NbHits[j]=0;
188
04b9b6ba 189 ntracks = (Int_t) diSim.NumberOfTracks(ievent);
b19dcb40 190
5707a463 191 for (itrack=0; itrack<ntracks; itrack++) { // track loop
04b9b6ba 192 AliMUONVHitStore* hitStore = diSim.HitStore(ievent,itrack);
193 AliMUONHit* mHit;
194 TIter next(hitStore->CreateIterator());
5707a463 195
04b9b6ba 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 }
522b4724 209 } // end track loop
07d83718 210
04b9b6ba 211
5707a463 212// 3/4 coincidence
213 SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
522b4724 214
5707a463 215 if (SumNbHits==3 || SumNbHits==4) {
216 coincmuon++;
217 ptcoinc->Fill(ptmu);
218 }
219
04b9b6ba 220
b19dcb40 221 } // end loop on event
04b9b6ba 222
b19dcb40 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 }
5707a463 231
522b4724 232
04b9b6ba 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
5707a463 335}