AliRoot Core  edcc906 (edcc906)
MUONTriggerEfficiencyPt.C
Go to the documentation of this file.
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 
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 }
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
Float_t Particle() const
Return particle id.
Definition: AliMUONHit.h:48
Virtual store to hold digit.
Bool_t IsValid() const
Returns true if the data interface was able to open the root file correctly.
An easy to use interface to MUON data.
Double_t fitArc(Double_t *x, Double_t *par)
TTree * tree
Base class of a trigger information store.
Int_t Track() const
Definition: AliHit.h:24
Int_t SingleLpt() const
Return number of Single Low pt.
static AliRunLoader * Open(const char *filename="galice.root", const char *eventfoldername=AliConfig::GetDefaultEventFolderName(), Option_t *option="READ")
virtual Int_t GetNtrack() const
Definition: AliStack.h:134
AliMUONVTriggerStore * TriggerStore(Int_t event, const char *treeLetter="R")
Int_t Chamber() const
Definition: AliMUONHit.cxx:160
virtual AliMUONGlobalTrigger * Global() const =0
Return global trigger.
void SetRun(Int_t run)
void MUONTriggerEfficiencyPt(const char *filenameSim="galice_sim.root", const char *filenameRec="galice.root", Bool_t readFromRP=0)
void SetDefaultStorage(const char *dbString)
Int_t GetEvent(Int_t evno)
Bool_t IsValid() const
Returns true if the data interface was able to open the root file correctly.
TParticle * Particle(Int_t id, Bool_t useInEmbedding=kFALSE)
Definition: AliStack.cxx:688
Easy to use data access to MC information.
MonteCarlo hit.
Definition: AliMUONHit.h:24
virtual TIterator * CreateIterator() const =0
Return an iterator to loop over hits.
Global trigger object.
AliMUONVHitStore * HitStore(Int_t event, Int_t track)
static AliCDBManager * Instance(TMap *entryCache=NULL, Int_t run=-1)
Int_t SingleHpt() const
Return number of Single High pt.
Double_t fitArch(Double_t *x, Double_t *par)
Int_t LoadKinematics(Option_t *option="READ")
AliStack * Stack() const
Definition: AliRunLoader.h:95