1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 /// \file MUONTriggerEfficiencyPt.C
21 /// \brief Macro to produce trigger single muon efficiency versus pt plots for the
24 /// See full description on the \ref README_trigger page.
26 /// \author Fabien Guerin (LPC)
30 #include "TClonesArray.h"
35 #include "TParticle.h"
41 #include "Riostream.h"
45 #include "AliRunLoader.h"
46 #include "AliHeader.h"
47 #include "AliLoader.h"
51 #include "AliMUONDataInterface.h"
52 #include "AliMUONMCDataInterface.h"
53 #include "AliMUONVHitStore.h"
54 #include "AliMUONVTriggerStore.h"
55 #include "AliMUONHit.h"
56 #include "AliMUONDigit.h"
57 #include "AliMUONGlobalTrigger.h"
58 #include "AliMUONTrack.h"
60 Double_t fitArc(Double_t *x,Double_t *par)
62 Double_t h1=par[1]*(x[0]+par[2]);
63 return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
66 Double_t fitArch(Double_t *x,Double_t *par)
68 Double_t h0=2*par[0]/(TMath::Pi());
69 Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0];
70 return h0*TMath::TanH(h1)+par[3];
73 void MUONTriggerEfficiencyPt(const char* filenameSim="galice_sim.root",
74 const char* filenameRec="galice.root",
75 Bool_t readFromRP = 0)
79 TStyle *st1 = new TStyle("st1","My STYLE");
83 st1->SetFrameFillColor(10);
84 st1->SetCanvasColor(10);
85 st1->SetFillColor(10);
86 st1->SetTitleBorderSize(0);
87 st1->SetTitleOffset(1.2,"XY");
88 st1->SetTitleSize(0.05,"XY");
89 st1->SetLabelSize(0.045,"XY");
90 st1->SetLabelFont(22,"XY");
91 st1->SetTitleFont(22,"XY");
93 st1->SetStatColor(10);
95 st1->SetFitFormat("4.4g");
96 st1->SetPadLeftMargin(0.15);
97 st1->SetPadRightMargin(0.06);
98 st1->SetPadTopMargin(0.06);
99 st1->SetPadBottomMargin(0.15);
102 // gROOT->ForceStyle();
103 //TGaxis::SetMaxDigits(3);
105 // beginning of macro
106 TH1F *ptlpt = new TH1F("ptlpt","",50,0.15,5.);
107 TH1F *pthpt = new TH1F("pthpt","",50,0.15,5.);
108 TH1F *ptcoinc = new TH1F("ptcoinc","",50,0.15,5.);
113 Double_t coincmuon=0;
118 AliMUONMCDataInterface diSim(filenameSim);
119 AliMUONDataInterface diRec(filenameRec);
121 if (!diSim.IsValid())
123 cerr << "Cannot access " << filenameSim << endl;
127 if (!diRec.IsValid())
129 cerr << "Cannot access " << filenameRec << endl;
133 FILE* fdat = fopen("MUONTriggerEfficiencyPt.out","w");
136 cerr << "Cannot create output file" << endl;
140 nevents = diSim.NumberOfEvents();
142 AliRunLoader* runLoader = AliRunLoader::Open(filenameSim);
144 for (Int_t ievent=0; ievent<nevents; ievent++) { // Event loop
146 if (ievent%500==0) printf("ievent = %d \n",ievent);
149 runLoader->GetRunLoader()->GetEvent(ievent);
150 runLoader->GetRunLoader()->LoadKinematics();
151 AliStack* stack = runLoader->GetRunLoader()->Stack();
153 Int_t nparticles = (Int_t) stack->GetNtrack();
155 for (Int_t iparticle=0; iparticle<nparticles; iparticle++) {
156 particle = stack->Particle(iparticle);
157 Float_t pt = particle->Pt();
158 Int_t pdgcode = particle->GetPdgCode();
159 if (pdgcode==-13 || pdgcode==13) ptmu = pt;
164 if ( readFromRP ) tree = "R";
166 AliMUONVTriggerStore* triggerStore = diRec.TriggerStore(ievent,tree.Data());
167 AliMUONGlobalTrigger* gloTrg = triggerStore->Global();
169 if (gloTrg->SingleLpt()>=1 ) {
173 if (gloTrg->SingleHpt()>=1 ) {
179 Int_t itrack, ntracks, NbHits[4];
182 for (Int_t j=0; j<4; j++) NbHits[j]=0;
184 ntracks = (Int_t) diSim.NumberOfTracks(ievent);
186 for (itrack=0; itrack<ntracks; itrack++) { // track loop
187 AliMUONVHitStore* hitStore = diSim.HitStore(ievent,itrack);
189 TIter next(hitStore->CreateIterator());
191 while ( ( mHit = static_cast<AliMUONHit*>(next()) ) )
193 Int_t Nch = mHit->Chamber();
194 Int_t hittrack = mHit->Track();
195 Float_t IdPart = mHit->Particle();
197 for (Int_t j=0;j<4;j++) {
200 if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++;
208 SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
210 if (SumNbHits==3 || SumNbHits==4) {
216 } // end loop on event
219 cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
220 cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
221 cout << " >>> muon track per event as it should " << "\n";
222 cout << " >>> see README for further information " << "\n";
223 cout << " >>> exiting now ! " << "\n";
230 fprintf(fdat," Number of events = %d \n",nevents);
231 fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
232 fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon);
233 fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);
236 Double_t efficiency,error;
238 efficiency=lptmuon/coincmuon;
239 error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));
240 fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
242 efficiency=hptmuon/coincmuon;
243 error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon));
244 fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
248 Double_t x1,x2,xval,xerr;
250 for (Int_t i=0;i<50;i++) {
251 x1=ptcoinc->GetBinContent(i+1);
253 x2=ptlpt->GetBinContent(i+1);
254 if (x1!=0 && x2!=0) {
256 xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
257 ptlpt->SetBinContent(i+1,xval);
258 ptlpt->SetBinError(i+1,0);
260 ptlpt->SetBinContent(i+1,0.);
261 ptlpt->SetBinError(i+1,0.);
264 x2=pthpt->GetBinContent(i+1);
265 if (x1!=0 && x2!=0) {
267 xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);
268 pthpt->SetBinContent(i+1,xval);
269 pthpt->SetBinError(i+1,0);
271 pthpt->SetBinContent(i+1,0.);
272 pthpt->SetBinError(i+1,0.);
276 TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4);
277 TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);
279 TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400);
283 ptlpt->SetMinimum(0.);
284 ptlpt->SetMaximum(1.05);
286 ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
287 ptlpt->GetYaxis()->SetTitle("Efficiency");
288 //ptlpt->GetXaxis()->SetRange(3);
290 fitlpt->SetParameters(0.602,0.774,0.273,0.048);
291 fitlpt->SetLineColor(2);
292 fitlpt->SetLineWidth(2);
293 fitlpt->Draw("SAME");
294 TLegend * leg = new TLegend(0.5,0.38,0.85,0.53);
295 leg = new TLegend(0.5,0.38,0.85,0.53);
296 leg->SetBorderSize(0);
297 leg->SetFillColor(0);
298 leg->SetTextSize(0.05);
299 leg->SetTextFont(22);
300 leg->SetHeader("Lpt trigger pt cut");
301 leg->AddEntry(fitlpt,"Ref.","l");
302 leg->AddEntry(ptlpt,"New","l");
306 pthpt->SetMinimum(0.);
307 pthpt->SetMaximum(1.05);
309 pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
310 pthpt->GetYaxis()->SetTitle("Efficiency");
311 //pthpt->GetXaxis()->SetRange(3);
313 fithpt->SetParameters(0.627,0.393,0.855,0.0081);
314 fithpt->SetLineColor(2);
315 fithpt->SetLineWidth(2);
316 fithpt->Draw("SAME");
317 leg = new TLegend(0.5,0.38,0.85,0.53);
318 leg->SetBorderSize(0);
319 leg->SetFillColor(0);
320 leg->SetTextSize(0.05);
321 leg->SetTextFont(22);
322 leg->SetHeader("Hpt trigger pt cut");
323 leg->AddEntry(fithpt,"Ref.","l");
324 leg->AddEntry(pthpt,"New","l");
327 c1->SaveAs("MUONTriggerEfficiencyPt.gif");
328 c1->SaveAs("MUONTriggerEfficiencyPt.eps");