fd09fca508f227de4cf7b0cb65a5c062cf323275
[u/mrichter/AliRoot.git] / PWGPP / MUON / dep / AliAnalysisTaskMuonPerformance.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 //-----------------------------------------------------------------------------
17 /// \class AliAnalysisTaskMuonPerformance
18 /// Analysis task to chek the tracker and trigger reconstruction
19 /// The output is a list of histograms.
20 ///
21 /// \author Diego Stocco and Philippe Pillot
22 //-----------------------------------------------------------------------------
23
24 // ROOT includes
25 #include "TROOT.h"
26 #include "TH1.h"
27 #include "TH2.h"
28 #include "TH3.h"
29 #include "TF1.h"
30 #include "TLegend.h"
31 #include "TGraphErrors.h"
32 #include "TGraphAsymmErrors.h"
33 #include "TAxis.h"
34 #include "TObjArray.h"
35 #include "TCanvas.h"
36 #include "TMath.h"
37 #include "TMCProcess.h"
38 #include "TGeoGlobalMagField.h"
39 #include "TGeoManager.h"
40
41 // STEER includes
42 #include "AliMCEvent.h"
43 #include "AliMCParticle.h"
44 #include "AliMCEventHandler.h"
45 #include "AliESDEvent.h"
46 #include "AliESDMuonTrack.h"
47 #include "AliCDBManager.h"
48 #include "AliGeomManager.h"
49
50 // ANALYSIS includes
51 #include "AliAnalysisManager.h"
52 #include "AliAnalysisDataSlot.h"
53 #include "AliAnalysisDataContainer.h"
54 #include "AliCentrality.h"
55
56 // CORRFW includes
57 #include "AliCFContainer.h"
58 #include "AliCFGridSparse.h"
59 #include "AliCFEffGrid.h"
60
61 // MUON includes
62 #include "AliMUONConstants.h"
63 #include "AliMUONTrack.h"
64 #include "AliMUONTriggerTrack.h"
65 #include "AliMUONLocalTrigger.h"
66 #include "AliMUONVTrackStore.h"
67 #include "AliMUONVTriggerTrackStore.h"
68 #include "AliMUONESDInterface.h"
69 #include "AliMUONRecoParam.h"
70 #include "AliMUONCDB.h"
71 #include "AliMUONTrackExtrap.h"
72 #include "AliMUONTrackParam.h"
73 #include "AliMUONRecoCheck.h"
74 #include "AliMUONVCluster.h"
75 #include "AliMUONVTrackReconstructor.h"
76
77 // MUON mapping includes
78 #include "AliMpSegmentation.h"
79 #include "AliMpDEIterator.h"
80
81 #include "AliAnalysisTaskMuonPerformance.h"
82
83
84 /// \cond CLASSIMP
85 ClassImp(AliAnalysisTaskMuonPerformance) // Class implementation in ROOT context
86 /// \endcond
87
88
89 //________________________________________________________________________
90 AliAnalysisTaskMuonPerformance::AliAnalysisTaskMuonPerformance() :
91 AliAnalysisTaskSE(),
92 fDefaultStorage(""),
93 fNPBins(30),
94 fCorrectForSystematics(kFALSE),
95 fFitResiduals(kFALSE),
96 fEnforceTrkCriteria(kFALSE),
97 fUseMCKinematics(kFALSE),
98 fMCTrigLevelFromMatchTrk(kFALSE),
99 fRequestedStationMask(0),
100 fRequest2ChInSameSt45(0),
101 fSigmaCutTrig(-1.),
102 fNDE(0),
103 fCFContainer(0x0),
104 fEfficiencyList(0x0),
105 fTriggerList(0x0),
106 fTrackerList(0x0),
107 fPAtVtxList(0x0),
108 fSlopeAtVtxList(0x0),
109 fEtaAtVtxList(0x0),
110 fPhiAtVtxList(0x0),
111 fPAt1stClList(0x0),
112 fSlopeAt1stClList(0x0),
113 fDCAList(0x0),
114 fClusterList(0x0)
115 {
116   //
117   /// Default Constructor.
118   //
119   fPRange[0] = 0.;
120   fPRange[1] = 0.;
121   fClusterMaxRes[0] = 0.;
122   fClusterMaxRes[1] = 0.;
123   for (Int_t i = 0; i < 1100; i++) fDEIndices[i] = 0;
124   for (Int_t i = 0; i < 200; i++) fDEIds[i] = 0;
125 }
126
127
128 //________________________________________________________________________
129 AliAnalysisTaskMuonPerformance::AliAnalysisTaskMuonPerformance(const char *name) :
130 AliAnalysisTaskSE(name),
131 fDefaultStorage("raw://"),
132 fNPBins(30),
133 fCorrectForSystematics(kFALSE),
134 fFitResiduals(kFALSE),
135 fEnforceTrkCriteria(kFALSE),
136 fUseMCKinematics(kFALSE),
137 fMCTrigLevelFromMatchTrk(kFALSE),
138 fRequestedStationMask(0),
139 fRequest2ChInSameSt45(0),
140 fSigmaCutTrig(-1.),
141 fNDE(0),
142 fCFContainer(0x0),
143 fEfficiencyList(0x0),
144 fTriggerList(0x0),
145 fTrackerList(0x0),
146 fPAtVtxList(0x0),
147 fSlopeAtVtxList(0x0),
148 fEtaAtVtxList(0x0),
149 fPhiAtVtxList(0x0),
150 fPAt1stClList(0x0),
151 fSlopeAt1stClList(0x0),
152 fDCAList(0x0),
153 fClusterList(0x0)
154 {
155   //
156   /// Constructor.
157   //
158   fPRange[0] = 0.;
159   fPRange[1] = 300.;
160   fClusterMaxRes[0] = 0.;
161   fClusterMaxRes[1] = 0.;
162   for (Int_t i = 0; i < 1100; i++) fDEIndices[i] = 0;
163   for (Int_t i = 0; i < 200; i++) fDEIds[i] = 0;
164   
165   DefineOutput(1, AliCFContainer::Class());
166   DefineOutput(2, TObjArray::Class());
167   DefineOutput(3, TObjArray::Class());
168   DefineOutput(4, TObjArray::Class());
169   DefineOutput(5, TObjArray::Class());
170   DefineOutput(6, TObjArray::Class());
171   DefineOutput(7, TObjArray::Class());
172   DefineOutput(8, TObjArray::Class());
173   DefineOutput(9, TObjArray::Class());
174   DefineOutput(10, TObjArray::Class());
175   DefineOutput(11, TObjArray::Class());
176   DefineOutput(12, TObjArray::Class());
177 }
178
179
180 //________________________________________________________________________
181 AliAnalysisTaskMuonPerformance::~AliAnalysisTaskMuonPerformance()
182 {
183   //
184   /// Destructor
185   //
186   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
187     delete fCFContainer;
188     delete fEfficiencyList;
189     delete fTriggerList;
190     delete fTrackerList;
191     delete fPAtVtxList;
192     delete fSlopeAtVtxList;
193     delete fEtaAtVtxList;
194     delete fPhiAtVtxList;
195     delete fPAt1stClList;
196     delete fSlopeAt1stClList;
197     delete fDCAList;
198     delete fClusterList;
199   }
200 }
201
202
203 //___________________________________________________________________________
204 void AliAnalysisTaskMuonPerformance::NotifyRun()
205 {
206   //
207   /// Use the event handler information to correctly fill the analysis flags:
208   /// - check if Monte Carlo information is present
209   //
210   
211   // load OCDB objects only once
212   if (fSigmaCutTrig > 0) return;
213   
214   // set OCDB location
215   AliCDBManager* cdbm = AliCDBManager::Instance();
216   if (cdbm->IsDefaultStorageSet()) printf("PerformanceTask: CDB default storage already set!\n");
217   else cdbm->SetDefaultStorage(fDefaultStorage.Data());
218   if (cdbm->GetRun() > -1) printf("PerformanceTask: run number already set!\n");
219   else cdbm->SetRun(fCurrentRunNumber);
220   
221   // load magnetic field for track extrapolation
222   if (!TGeoGlobalMagField::Instance()->GetField()) {
223     if (!AliMUONCDB::LoadField()) return;
224   }
225   
226   // load mapping
227   if (!AliMpSegmentation::Instance(kFALSE)) {
228     if (!AliMUONCDB::LoadMapping(kTRUE)) return;
229   }
230
231   // load geometry for track extrapolation to vertex and for checking hits are under pads in reconstructible tracks
232   if (!AliGeomManager::GetGeometry()) {
233     AliGeomManager::LoadGeometry();
234     if (!AliGeomManager::GetGeometry()) return;  
235     if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
236   }
237   
238   // load recoParam
239   const AliMUONRecoParam* recoParam = (AliMUONESDInterface::GetTracker())
240   ? AliMUONESDInterface::GetTracker()->GetRecoParam()
241   : AliMUONCDB::LoadRecoParam();
242   
243   if (!recoParam) {
244     fRequestedStationMask = 0;
245     fRequest2ChInSameSt45 = kFALSE;
246     fSigmaCutTrig = -1.;
247     AliError("--> skip this run");
248     return;
249   }
250   
251   // compute the mask of requested stations from recoParam
252   fRequestedStationMask = 0;
253   for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) fRequestedStationMask |= ( 1 << i );
254   
255   // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
256   fRequest2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
257   
258   // get sigma cut from recoParam to associate trigger track to triggerable track
259   fSigmaCutTrig = recoParam->GetSigmaCutForTrigger();
260   
261   if (!AliMUONESDInterface::GetTracker()) AliMUONESDInterface::ResetTracker(recoParam);
262   
263   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
264     
265     // find the highest chamber resolution and set histogram bins
266     if (recoParam->GetDefaultNonBendingReso(i) > fClusterMaxRes[0]) fClusterMaxRes[0] = recoParam->GetDefaultNonBendingReso(i);
267     if (recoParam->GetDefaultBendingReso(i) > fClusterMaxRes[1]) fClusterMaxRes[1] = recoParam->GetDefaultBendingReso(i);
268     
269     // fill correspondence between DEId and indices for histo (starting from 1)
270     AliMpDEIterator it;
271     it.First(i);
272     while (!it.IsDone()) {
273       fNDE++;
274       fDEIndices[it.CurrentDEId()] = fNDE;
275       fDEIds[fNDE] = it.CurrentDEId();
276       it.Next();
277     }
278     
279   }
280   
281   UserCreateOutputObjects();
282 }
283
284
285 //___________________________________________________________________________
286 void AliAnalysisTaskMuonPerformance::UserCreateOutputObjects() 
287 {
288   //
289   /// Create output objects
290   //
291
292   // do it once the OCDB has been loaded (i.e. from NotifyRun())
293   if (fSigmaCutTrig < 0) return;
294   
295   // ------ CFContainer ------
296   
297   // define axes of particle container
298   Int_t nPtBins = 30;
299   Double_t ptMin = 0., ptMax = 15.;
300   TString ptTitle("p_{t}"), ptUnits("GeV/c");
301
302   Int_t nEtaBins = 15;
303   Double_t etaMin = -4., etaMax = -2.5;
304   TString etaTitle("#eta"), etaUnits("a.u.");
305
306   Int_t nPhiBins = 15;
307   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
308   TString phiTitle("#phi"), phiUnits("rad");
309
310   Int_t nThetaAbsEndBins = 4;
311   Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 3.5;
312   TString thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");  
313
314   Int_t nChargeBins = 2;
315   Double_t chargeMin = -2, chargeMax = 2.;
316   TString chargeTitle("charge"), chargeUnits("e");
317
318   Int_t nHasTrackerBins = 2;
319   Double_t hasTrackerMin = -0.5, hasTrackerMax = (Double_t)nHasTrackerBins - 0.5;
320   TString hasTrackerTitle("Has tracker"), hasTrackerUnits("");
321   
322   Int_t nTriggerBins = kNtrigCuts;
323   Double_t triggerMin = -0.5, triggerMax = (Double_t)nTriggerBins - 0.5;
324   TString triggerTitle("Trigger match"), triggerUnits("");
325
326   Int_t nMotherTypeBins = kNtrackSources;
327   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
328   TString motherTypeTitle("motherType"), motherTypeUnits("");
329
330   Int_t nMatchMCBins = kNMatchMC;
331   Double_t matchMCMin = -0.5, matchMCMax = (Double_t)kNMatchMC - 0.5;
332   TString matchMCTitle("MatchMC"), matchMCUnits("");
333   
334   Int_t nMCTriggerBins = kNtrigCuts;
335   Double_t mcTriggerMin = -0.5, mcTriggerMax = (Double_t)nMCTriggerBins - 0.5;
336   TString mcTriggerTitle("MC Trigger match"), mcTriggerUnits("");
337   
338   Int_t nCentBins = 22;
339   Double_t centMin = -5., centMax = 105.;
340   TString centTitle("centrality"), centUnits("%");
341   
342   Int_t nDupliTrgBins = 2;
343   Double_t dupliTrgMin = -0.5, dupliTrgMax = 1.5;
344   TString dupliTrgTitle("duplicate trigger"), dupliTrgUnits("");
345   
346   Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nThetaAbsEndBins, nChargeBins, nHasTrackerBins, nTriggerBins, nMotherTypeBins, nMatchMCBins, nMCTriggerBins, nCentBins, nDupliTrgBins};
347   Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, thetaAbsEndMin, chargeMin, hasTrackerMin, triggerMin, motherTypeMin, matchMCMin, mcTriggerMin, centMin, dupliTrgMin};
348   Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, thetaAbsEndMax, chargeMax, hasTrackerMax, triggerMax, motherTypeMax, matchMCMax, mcTriggerMax, centMax, dupliTrgMax};
349   TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, thetaAbsEndTitle, chargeTitle, hasTrackerTitle, triggerTitle, motherTypeTitle, matchMCTitle, mcTriggerTitle, centTitle, dupliTrgTitle};
350   TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, thetaAbsEndUnits, chargeUnits, hasTrackerUnits, triggerUnits, motherTypeUnits, matchMCUnits, mcTriggerUnits, centUnits, dupliTrgUnits};
351   
352   // create particle container
353   fCFContainer = new AliCFContainer(GetOutputSlot(1)->GetContainer()->GetName(),"container for tracks",kNsteps,kNvars,nbins);
354   
355   // set axes title and limits
356   for ( Int_t idim = 0; idim<kNvars; idim++) {
357     TString histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
358     histoTitle.ReplaceAll("()","");
359     fCFContainer->SetVarTitle(idim, histoTitle.Data());
360     fCFContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
361   }
362   
363   // define histo names
364   TString stepTitle[kNsteps] = {"reconstructed", "generated"};
365   
366   // define axes labels if any
367   TString trigName[kNtrigCuts];
368   trigName[kNoMatchTrig] = "NoMatch";
369   trigName[kOtherTrig]   = "Other";
370   trigName[kAllPtTrig]   = "AllPt";
371   trigName[kLowPtTrig]   = "LowPt";
372   trigName[kHighPtTrig]  = "HighPt";
373   
374   TString srcName[kNtrackSources];
375   srcName[kCharmMu]     = "Charm";
376   srcName[kBeautyMu]    = "Beauty";
377   srcName[kPrimaryMu]   = "Decay";
378   srcName[kSecondaryMu] = "Secondary";
379   srcName[kRecoHadron]  = "Hadrons";
380   srcName[kUnknownPart] = "Fakes";
381   
382   TString mMCName[kNMatchMC];
383   mMCName[kNoMatch]     = "NoMatch";
384   mMCName[kTrackerOnly] = "TrackerOnly";
385   mMCName[kMatchedSame] = "MatchedSame";
386   mMCName[kMatchedDiff] = "MatchedDiff";
387   mMCName[kTriggerOnly] = "TriggerOnly";
388   
389   // set histo name and axis labels if any
390   for (Int_t istep=0; istep<kNsteps; istep++) {
391     
392     // histo name
393     fCFContainer->SetStepTitle(istep, stepTitle[istep].Data());
394     AliCFGridSparse* gridSparse = fCFContainer->GetGrid(istep);
395     
396     // trigger info
397     TAxis* triggerAxis =  gridSparse->GetAxis(kVarTrigger);
398     for ( Int_t ibin=0; ibin<kNtrigCuts; ibin++ ) {
399       triggerAxis->SetBinLabel(ibin+1,trigName[ibin]);
400     }
401     
402     // mother type
403     TAxis* motherTypeAxis = gridSparse->GetAxis(kVarMotherType);
404     for ( Int_t ibin=0; ibin<kNtrackSources; ibin++ ) {
405       motherTypeAxis->SetBinLabel(ibin+1,srcName[ibin]);
406     }
407     
408     // MC matching flag
409     TAxis* matchMCAxis = gridSparse->GetAxis(kVarMatchMC);
410     for ( Int_t ibin=0; ibin<kNMatchMC; ibin++ ) {
411       matchMCAxis->SetBinLabel(ibin+1,mMCName[ibin]);
412     }
413     
414     // MC trigger info
415     TAxis* mcTriggerAxis =  gridSparse->GetAxis(kVarMCTrigger);
416     for ( Int_t ibin=0; ibin<kNtrigCuts; ibin++ ) {
417       mcTriggerAxis->SetBinLabel(ibin+1,trigName[ibin]);
418     }
419     
420   }
421   
422   // ------ trigger histograms ------
423   
424   fTriggerList = new TObjArray(100);
425   fTriggerList->SetOwner();
426   
427   TH1F* h1 = new TH1F("hResTrigX11", "Residual X11;X11_{reco} - X11_{MC} (cm)", 100, -10., 10.);
428   fTriggerList->AddAt(h1, kResTrigX11);
429   h1 = new TH1F("hResTrigY11", "Residual Y11;Y11_{reco} - Y11_{MC} (cm)", 100, -10., 10.);
430   fTriggerList->AddAt(h1, kResTrigY11);
431   h1 = new TH1F("hResTrigSlopeY", "Residual slope y;ySlope_{reco} - ySlope_{MC} (rad)", 100, -0.1, 0.1);
432   fTriggerList->AddAt(h1, kResTrigSlopeY);
433   
434   // ------ tracker histograms ------
435   
436   fTrackerList = new TObjArray(100);
437   fTrackerList->SetOwner();
438   
439   // momentum resolution at vertex
440   const Int_t deltaPAtVtxNBins = 250;
441   Double_t deltaPAtVtxEdges[2];
442   deltaPAtVtxEdges[0] = -20. - 0.05 * fPRange[1];
443   deltaPAtVtxEdges[1] =   5. + 0.05 * fPRange[1];
444   
445   h1 = new TH1F("hResPAtVtx"," delta P at vertex;#Delta_{p} (GeV/c)",deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
446   fTrackerList->AddAt(h1, kResPAtVtx);
447   TH2F *h2 = new TH2F("hResPAtVtxVsP","#Delta_{p} at vertex versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
448   fTrackerList->AddAt(h2, kResPAtVtxVsP);
449   h2 = new TH2F("hResPAtVtxVsPIn23deg","#Delta_{p} at vertex versus p for tracks between 2 and 3 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
450   fTrackerList->AddAt(h2, kResPAtVtxVsPIn23deg);
451   h2 = new TH2F("hResPAtVtxVsPIn310deg","#Delta_{p} at vertex versus p for tracks between 3 and 10 degrees at absorber end;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
452   fTrackerList->AddAt(h2, kResPAtVtxVsPIn310deg);
453   h2 = new TH2F("hResPAtVtxVsPIn02degMC","#Delta_{p} at vertex versus p for tracks with MC angle below 2 degrees;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtVtxNBins/10,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
454   fTrackerList->AddAt(h2, kResPAtVtxVsPIn02degMC);
455   h2 = new TH2F("hResPAtVtxVsPosAbsEndIn02degMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
456   fTrackerList->AddAt(h2, kResPAtVtxVsPosAbsEndIn02degMC);
457   h2 = new TH2F("hResPAtVtxVsPosAbsEndIn23degMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
458   fTrackerList->AddAt(h2, kResPAtVtxVsPosAbsEndIn23degMC);
459   h2 = new TH2F("hResPAtVtxVsPosAbsEndIn310degMC","#Delta_{p} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{p} (GeV/c)",1000,0.,100.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
460   fTrackerList->AddAt(h2, kResPAtVtxVsPosAbsEndIn310degMC);
461   h2 = new TH2F("hResPAtVtxVsAngleAtAbsEnd","#Delta_{p} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
462   fTrackerList->AddAt(h2, kResPAtVtxVsAngleAtAbsEnd);
463   h2 = new TH2F("hResPAtVtxVsMCAngle","#Delta_{p} at vertex versus MC angle;MC angle (Deg);#Delta_{p} (GeV/c)",10,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
464   fTrackerList->AddAt(h2, kResPAtVtxVsMCAngle);
465   TH3F *h3 = new TH3F("hResPAtVtxVsAngleAtAbsEndVsP","#Delta_{p} at vertex versus track position at absorber end converted to degrees versus momentum;p (GeV/c);angle (Deg);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],100,0.,10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0],deltaPAtVtxEdges[1]);
466   fTrackerList->AddAt(h3, kResPAtVtxVsAngleAtAbsEndVsP);
467   
468   // transverse momentum resolution at vertex
469   h2 = new TH2F("hResPtAtVtxVsPt","#Delta_{p_{t}} at vertex versus p_{t};p_{t} (GeV/c);#Delta_{p_{t}} (GeV/c)",2*fNPBins,fPRange[0]/10.,fPRange[1]/10.,deltaPAtVtxNBins,deltaPAtVtxEdges[0]/10.,deltaPAtVtxEdges[1]/10.);
470   fTrackerList->AddAt(h2, kResPtAtVtxVsPt);
471   
472   // momentum resolution at first cluster
473   const Int_t deltaPAtFirstClNBins = 500;
474   Double_t deltaPAtFirstClEdges[2];
475   deltaPAtFirstClEdges[0] = -5. - 0.05 * fPRange[1];
476   deltaPAtFirstClEdges[1] =  5. + 0.05 * fPRange[1];
477   
478   h1 = new TH1F("hResPAt1stCl"," delta P at first cluster;#Delta_{p} (GeV/c)",deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
479   fTrackerList->AddAt(h1, kResPAt1stCl);
480   h2 = new TH2F("hResPAt1stClVsP","#Delta_{p} at first cluster versus p;p (GeV/c);#Delta_{p} (GeV/c)",2*fNPBins,fPRange[0],fPRange[1],deltaPAtFirstClNBins,deltaPAtFirstClEdges[0],deltaPAtFirstClEdges[1]);
481   fTrackerList->AddAt(h2, kResPAt1stClVsP);
482   
483   // transverse momentum resolution at vertex
484   h2 = new TH2F("hResPtAt1stClVsPt","#Delta_{p_{t}} at first cluster versus p_{t};p_{t} (GeV/c);#Delta_{p_{t}} (GeV/c)",2*fNPBins,fPRange[0]/10.,fPRange[1]/10.,deltaPAtFirstClNBins,deltaPAtFirstClEdges[0]/10.,deltaPAtFirstClEdges[1]/10.);
485   fTrackerList->AddAt(h2, kResPtAt1stClVsPt);
486   
487   // angular resolution at vertex
488   const Int_t deltaSlopeAtVtxNBins = 500;
489   const Double_t deltaSlopeAtVtxEdges[2] = {-0.05, 0.05};
490   
491   h1 = new TH1F("hResSlopeXAtVtx","#Delta_{slope_{X}} at vertex;#Delta_{slope_{X}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
492   fTrackerList->AddAt(h1, kResSlopeXAtVtx);
493   h1 = new TH1F("hResSlopeYAtVtx","#Delta_{slope_{Y}} at vertex;#Delta_{slope_{Y}}", deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
494   fTrackerList->AddAt(h1, kResSlopeYAtVtx);
495   h2 = new TH2F("hResSlopeXAtVtxVsP","#Delta_{slope_{X}} at vertex versus p;p (GeV/c);#Delta_{slope_{X}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
496   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsP);
497   h2 = new TH2F("hResSlopeYAtVtxVsP","#Delta_{slope_{Y}} at vertex versus p;p (GeV/c);#Delta_{slope_{Y}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
498   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsP);
499   h2 = new TH2F("hResSlopeXAtVtxVsPosAbsEndIn02degMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
500   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsPosAbsEndIn02degMC);
501   h2 = new TH2F("hResSlopeYAtVtxVsPosAbsEndIn02degMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
502   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsPosAbsEndIn02degMC);
503   h2 = new TH2F("hResSlopeXAtVtxVsPosAbsEndIn23degMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
504   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsPosAbsEndIn23degMC);
505   h2 = new TH2F("hResSlopeYAtVtxVsPosAbsEndIn23degMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
506   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsPosAbsEndIn23degMC);
507   h2 = new TH2F("hResSlopeXAtVtxVsPosAbsEndIn310degMC","#Delta_{slope_{X}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{X}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
508   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsPosAbsEndIn310degMC);
509   h2 = new TH2F("hResSlopeYAtVtxVsPosAbsEndIn310degMC","#Delta_{slope_{Y}} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{slope_{Y}}",1000,0.,100.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
510   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsPosAbsEndIn310degMC);
511   h2 = new TH2F("hResSlopeXAtVtxVsAngleAtAbsEnd","#Delta_{slope_{X}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
512   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsAngleAtAbsEnd);
513   h2 = new TH2F("hResSlopeYAtVtxVsAngleAtAbsEnd","#Delta_{slope_{Y}} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
514   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsAngleAtAbsEnd);
515   h2 = new TH2F("hResSlopeXAtVtxVsMCAngle","#Delta_{slope_{X}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{X}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
516   fTrackerList->AddAt(h2, kResSlopeXAtVtxVsMCAngle);
517   h2 = new TH2F("hResSlopeYAtVtxVsMCAngle","#Delta_{slope_{Y}} at vertex versus MC angle;MC angle (Deg);#Delta_{slope_{Y}}",10,0.,10.,deltaSlopeAtVtxNBins, deltaSlopeAtVtxEdges[0], deltaSlopeAtVtxEdges[1]);
518   fTrackerList->AddAt(h2, kResSlopeYAtVtxVsMCAngle);
519   
520   // angular resolution at first cluster
521   const Int_t deltaSlopeAtFirstClNBins = 500;
522   const Double_t deltaSlopeAtFirstClEdges[2] = {-0.01, 0.01};
523   
524   h1 = new TH1F("hResSlopeXAt1stCl","#Delta_{slope_{X}} at first cluster;#Delta_{slope_{X}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
525   fTrackerList->AddAt(h1, kResSlopeXAt1stCl);
526   h1 = new TH1F("hResSlopeYAt1stCl","#Delta_{slope_{Y}} at first cluster;#Delta_{slope_{Y}}", deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
527   fTrackerList->AddAt(h1, kResSlopeYAt1stCl);
528   h2 = new TH2F("hResSlopeXAt1stClVsP","#Delta_{slope_{X}} at first cluster versus p;p (GeV/c);#Delta_{slope_{X}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
529   fTrackerList->AddAt(h2, kResSlopeXAt1stClVsP);
530   h2 = new TH2F("hResSlopeYAt1stClVsP","#Delta_{slope_{Y}} at first cluster versus p;p (GeV/c);#Delta_{slope_{Y}}",2*fNPBins,fPRange[0],fPRange[1], deltaSlopeAtFirstClNBins, deltaSlopeAtFirstClEdges[0], deltaSlopeAtFirstClEdges[1]);
531   fTrackerList->AddAt(h2, kResSlopeYAt1stClVsP);
532   
533   // eta resolution at vertex
534   const Int_t deltaEtaAtVtxNBins = 500;
535   const Double_t deltaEtaAtVtxEdges[2] = {-0.5, 0.5};
536   
537   h1 = new TH1F("hResEtaAtVtx","#Delta_{eta} at vertex;#Delta_{eta}", deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
538   fTrackerList->AddAt(h1, kResEtaAtVtx);
539   h2 = new TH2F("hResEtaAtVtxVsP","#Delta_{eta} at vertex versus p;p (GeV/c);#Delta_{eta}",2*fNPBins,fPRange[0],fPRange[1], deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
540   fTrackerList->AddAt(h2, kResEtaAtVtxVsP);
541   h2 = new TH2F("hResEtaAtVtxVsPosAbsEndIn02degMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
542   fTrackerList->AddAt(h2, kResEtaAtVtxVsPosAbsEndIn02degMC);
543   h2 = new TH2F("hResEtaAtVtxVsPosAbsEndIn23degMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
544   fTrackerList->AddAt(h2, kResEtaAtVtxVsPosAbsEndIn23degMC);
545   h2 = new TH2F("hResEtaAtVtxVsPosAbsEndIn310degMC","#Delta_{eta} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{eta}",1000,0.,100.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
546   fTrackerList->AddAt(h2, kResEtaAtVtxVsPosAbsEndIn310degMC);
547   h2 = new TH2F("hResEtaAtVtxVsAngleAtAbsEnd","#Delta_{eta} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
548   fTrackerList->AddAt(h2, kResEtaAtVtxVsAngleAtAbsEnd);
549   h2 = new TH2F("hResEtaAtVtxVsMCAngle","#Delta_{eta} at vertex versus MC angle;MC angle (Deg);#Delta_{eta}",10,0.,10.,deltaEtaAtVtxNBins, deltaEtaAtVtxEdges[0], deltaEtaAtVtxEdges[1]);
550   fTrackerList->AddAt(h2, kResEtaAtVtxVsMCAngle);
551   
552   // phi resolution at vertex
553   const Int_t deltaPhiAtVtxNBins = 500;
554   const Double_t deltaPhiAtVtxEdges[2] = {-0.5, 0.5};
555   
556   h1 = new TH1F("hResPhiAtVtx","#Delta_{phi} at vertex;#Delta_{phi}", deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
557   fTrackerList->AddAt(h1, kResPhiAtVtx);
558   h2 = new TH2F("hResPhiAtVtxVsP","#Delta_{phi} at vertex versus p;p (GeV/c);#Delta_{phi}",2*fNPBins,fPRange[0],fPRange[1], deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
559   fTrackerList->AddAt(h2, kResPhiAtVtxVsP);
560   h2 = new TH2F("hResPhiAtVtxVsPosAbsEndIn02degMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
561   fTrackerList->AddAt(h2, kResPhiAtVtxVsPosAbsEndIn02degMC);
562   h2 = new TH2F("hResPhiAtVtxVsPosAbsEndIn23degMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
563   fTrackerList->AddAt(h2, kResPhiAtVtxVsPosAbsEndIn23degMC);
564   h2 = new TH2F("hResPhiAtVtxVsPosAbsEndIn310degMC","#Delta_{phi} at vertex versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);#Delta_{phi}",1000,0.,100.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
565   fTrackerList->AddAt(h2, kResPhiAtVtxVsPosAbsEndIn310degMC);
566   h2 = new TH2F("hResPhiAtVtxVsAngleAtAbsEnd","#Delta_{phi} at vertex versus track position at absorber end converted to degrees;angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
567   fTrackerList->AddAt(h2, kResPhiAtVtxVsAngleAtAbsEnd);
568   h2 = new TH2F("hResPhiAtVtxVsMCAngle","#Delta_{phi} at vertex versus MC angle;MC angle (Deg);#Delta_{phi}",10,0.,10.,deltaPhiAtVtxNBins, deltaPhiAtVtxEdges[0], deltaPhiAtVtxEdges[1]);
569   fTrackerList->AddAt(h2, kResPhiAtVtxVsMCAngle);
570   
571   // DCA resolution
572   const Int_t deltaPDCANBins = 500;
573   const Double_t deltaPDCAEdges[2] = {0., 1000.};
574   const Double_t deltaPMCSAngEdges[2] = {-1.5, 1.5};
575   
576   h1 = new TH1F("hPDCA","p #times DCA at vertex;p #times DCA (GeV #times cm)", deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
577   fTrackerList->AddAt(h1, kPDCA);
578   h2 = new TH2F("hPDCAVsPIn23deg","p #times DCA versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
579   fTrackerList->AddAt(h2, kPDCAVsPIn23deg);
580   h2 = new TH2F("hPDCAVsPIn310deg","p #times DCA versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);p #times DCA (GeV #times cm)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
581   fTrackerList->AddAt(h2, kPDCAVsPIn310deg);
582   h2 = new TH2F("hPDCAVsPosAbsEndIn02degMC","p #times DCA versus track position at absorber end for tracks with MC angle < 2 degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
583   fTrackerList->AddAt(h2, kPDCAVsPosAbsEndIn02degMC);
584   h2 = new TH2F("hPDCAVsPosAbsEndIn23degMC","p #times DCA}versus track position at absorber end for tracks with MC angle in [2,3[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
585   fTrackerList->AddAt(h2, kPDCAVsPosAbsEndIn23degMC);
586   h2 = new TH2F("hPDCAVsPosAbsEndIn310degMC","p #times DCA versus track position at absorber end for tracks with MC angle in [3,10[ degrees;position (cm);p #times DCA (GeV #times cm)",1000,0.,100.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
587   fTrackerList->AddAt(h2, kPDCAVsPosAbsEndIn310degMC);
588   h2 = new TH2F("hPDCAVsAngleAtAbsEnd","p #times DCA versus track position at absorber end converted to degrees;angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
589   fTrackerList->AddAt(h2, kPDCAVsAngleAtAbsEnd);
590   h2 = new TH2F("hPDCAVsMCAngle","p #times DCA versus MC angle;MC angle (Deg);p #times DCA (GeV #times cm)",10,0.,10.,deltaPDCANBins, deltaPDCAEdges[0], deltaPDCAEdges[1]);
591   fTrackerList->AddAt(h2, kPDCAVsMCAngle);
592   
593   // MCS angular dispersion
594   h2 = new TH2F("hPMCSAngVsPIn23deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
595   fTrackerList->AddAt(h2, kPMCSAngVsPIn23deg);
596   h2 = new TH2F("hPMCSAngVsPIn310deg","p #times #Delta#theta_{MCS} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);p #times #Delta#theta_{MCS} (GeV)",2*fNPBins,fPRange[0],fPRange[1], deltaPDCANBins, deltaPMCSAngEdges[0], deltaPMCSAngEdges[1]);
597   fTrackerList->AddAt(h2, kPMCSAngVsPIn310deg);
598   
599   // cluster resolution
600   Int_t nCh = AliMUONConstants::NTrackingCh();
601   const Int_t clusterResNBins = 5000;
602   Double_t clusterResMaxX = 10.*fClusterMaxRes[0];
603   Double_t clusterResMaxY = 10.*fClusterMaxRes[1];
604   
605   h2 = new TH2F("hResClXVsCh", "cluster-track residual-X distribution per chamber;chamber ID;#Delta_{X} (cm)", nCh, 0.5, nCh+0.5, clusterResNBins, -clusterResMaxX, clusterResMaxX);
606   fTrackerList->AddAt(h2, kResClXVsCh);
607   h2 = new TH2F("hResClYVsCh", "cluster-track residual-Y distribution per chamber;chamber ID;#Delta_{Y} (cm)", nCh, 0.5, nCh+0.5, clusterResNBins, -clusterResMaxY, clusterResMaxY);
608   fTrackerList->AddAt(h2, kResClYVsCh);
609   h2 = new TH2F("hResClXVsDE", "cluster-track residual-X distribution per DE;DE ID;#Delta_{X} (cm)", fNDE, 0.5, fNDE+0.5, clusterResNBins, -clusterResMaxX, clusterResMaxX);
610   for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
611   fTrackerList->AddAt(h2, kResClXVsDE);
612   h2 = new TH2F("hResClYVsDE", "cluster-track residual-Y distribution per DE;DE ID;#Delta_{Y} (cm)", fNDE, 0.5, fNDE+0.5, clusterResNBins, -clusterResMaxY, clusterResMaxY);
613   for (Int_t i = 1; i <= fNDE; i++) h2->GetXaxis()->SetBinLabel(i, Form("%d",fDEIds[i]));
614   fTrackerList->AddAt(h2, kResClYVsDE);
615   
616   // Disable printout of AliMCEvent
617   AliLog::SetClassDebugLevel("AliMCEvent",-1);
618   
619   PostData(1, fCFContainer);
620   PostData(2, fTriggerList);
621   PostData(3, fTrackerList);
622 }
623
624 //________________________________________________________________________
625 void AliAnalysisTaskMuonPerformance::UserExec(Option_t * /*option*/) 
626 {
627   //
628   /// Main loop
629   /// Called for each event
630   //
631
632   // check that OCDB objects have been properly set
633   if (fSigmaCutTrig < 0) return;
634   
635   // Load ESD event
636   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
637   if ( ! esd ) {
638     AliError ("ESD event not found. Nothing done!");
639     return;
640   }
641
642   // Load MC event 
643   AliMCEventHandler* mcH = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
644   if ( ! mcH ) {
645     AliError ("MCH event handler not found. Nothing done!");
646     return;
647   }
648   
649   // get centrality
650   Double_t centrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
651   
652   // Get reference tracks
653   AliMUONRecoCheck rc(esd,mcH);
654   AliMUONVTriggerTrackStore* triggerTrackRefStore = rc.TriggerableTracks(-1);
655   AliMUONVTrackStore* trackRefStore = rc.TrackRefs(-1);
656   AliMUONVTrackStore* reconstructibleStore = rc.ReconstructibleTracks(-1, fRequestedStationMask, fRequest2ChInSameSt45);
657   
658   Double_t containerInput[kNvars];
659   containerInput[kVarCent] = centrality;
660   AliMUONTrackParam *trackParam;
661   Double_t x1,y1,z1,slopex1,slopey1,pX1,pY1,pZ1,p1,pT1,eta1,phi1;
662   Double_t x2,y2,z2,slopex2,slopey2,pX2,pY2,pZ2,p2,pT2,eta2,phi2;
663   Double_t dPhi;
664   Double_t xAbs,yAbs,dAbs,aAbs,aMCS,aMC;
665   Double_t xDCA,yDCA,dca,pU;
666   
667   // ------ Loop over reconstructed tracks ------
668   AliESDMuonTrack *esdTrack = 0x0;
669   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
670   Int_t *loCircuit = new Int_t[nMuTracks];
671   Int_t nTrgTracks = 0;
672   for (Int_t iMuTrack = 0; iMuTrack < nMuTracks; ++iMuTrack) {
673     
674     esdTrack = esd->GetMuonTrack(iMuTrack);
675
676     // Tracker
677     AliMUONTrack* matchedTrackRef = 0x0;
678     containerInput[kVarMatchMC] = static_cast<Double_t>(kNoMatch);
679     Bool_t isValid = kFALSE;
680     if (esdTrack->ContainTrackerData()) {
681     
682       // convert ESD track to MUON track (without recomputing track parameters at each clusters)
683       AliMUONTrack muonTrack;
684       AliMUONESDInterface::ESDToMUON(*esdTrack, muonTrack, kFALSE);
685       
686       // decide whether the track is valid for efficiency calculations
687       isValid = (!fEnforceTrkCriteria || muonTrack.IsValid(fRequestedStationMask, fRequest2ChInSameSt45));
688       
689       // get the associated simulated track (discard decays)
690       Int_t mcLabel = esdTrack->GetLabel();
691       if (mcLabel >= 0 && !esdTrack->TestBit(BIT(22)))
692         matchedTrackRef = static_cast<AliMUONTrack*>(trackRefStore->FindObject(mcLabel));
693       if (matchedTrackRef && isValid) containerInput[kVarMatchMC] = static_cast<Double_t>(kTrackerOnly);
694       
695       // compute track resolution (discard not-reconstructible trackRef)
696       if (matchedTrackRef && !esdTrack->TestBit(BIT(23))) {
697         
698         // simulated track parameters at vertex
699         trackParam = matchedTrackRef->GetTrackParamAtVertex();
700         x1 = trackParam->GetNonBendingCoor();
701         y1 = trackParam->GetBendingCoor();
702         z1 = trackParam->GetZ();
703         slopex1 = trackParam->GetNonBendingSlope();
704         slopey1 = trackParam->GetBendingSlope();
705         pX1 = trackParam->Px();
706         pY1 = trackParam->Py();
707         pZ1 = trackParam->Pz();
708         p1  = trackParam->P();
709         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
710         aMC = TMath::ATan(-pT1/pZ1) * TMath::RadToDeg();
711         eta1 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT1/pZ1)));
712         phi1 = TMath::Pi()+TMath::ATan2(-pY1, -pX1);
713         
714         // reconstructed track parameters at the end of the absorber
715         AliMUONTrackParam trackParamAtAbsEnd(*((AliMUONTrackParam*)muonTrack.GetTrackParamAtCluster()->First()));
716         AliMUONTrackExtrap::ExtrapToZ(&trackParamAtAbsEnd, AliMUONConstants::AbsZEnd());
717         xAbs = trackParamAtAbsEnd.GetNonBendingCoor();
718         yAbs = trackParamAtAbsEnd.GetBendingCoor();
719         dAbs = TMath::Sqrt(xAbs*xAbs + yAbs*yAbs);
720         aAbs = TMath::ATan(-dAbs/AliMUONConstants::AbsZEnd()) * TMath::RadToDeg();
721         pX2 = trackParamAtAbsEnd.Px();
722         pY2 = trackParamAtAbsEnd.Py();
723         pZ2 = trackParamAtAbsEnd.Pz();
724         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
725         aMCS = TMath::ATan(-pT2/pZ2) * TMath::RadToDeg();
726         
727         // reconstructed track parameters at vertex
728         trackParam = muonTrack.GetTrackParamAtVertex();
729         x2 = trackParam->GetNonBendingCoor();
730         y2 = trackParam->GetBendingCoor();
731         z2 = trackParam->GetZ();
732         slopex2 = trackParam->GetNonBendingSlope();
733         slopey2 = trackParam->GetBendingSlope();
734         pX2 = trackParam->Px();
735         pY2 = trackParam->Py();
736         pZ2 = trackParam->Pz();
737         p2  = trackParam->P();
738         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
739         eta2 = TMath::Log(TMath::Tan(0.5*TMath::ATan(-pT2/pZ2)));
740         phi2 = TMath::Pi()+TMath::ATan2(-pY2, -pX2);
741         
742         // reconstructed track parameters at DCA
743         AliMUONTrackParam trackParamAtDCA(*((AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->First()));
744         pU = trackParamAtDCA.P();
745         AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&trackParamAtDCA, z2);
746         xDCA = trackParamAtDCA.GetNonBendingCoor();
747         yDCA = trackParamAtDCA.GetBendingCoor();
748         dca = TMath::Sqrt(xDCA*xDCA + yDCA*yDCA);
749         
750         dPhi = phi2-phi1;
751         if (dPhi < -TMath::Pi()) dPhi += 2.*TMath::Pi();
752         else if (dPhi > TMath::Pi()) dPhi -= 2.*TMath::Pi();
753         
754         // fill histograms
755         static_cast<TH1*>(fTrackerList->At(kResPAtVtx))->Fill(p2-p1);
756         static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsP))->Fill(p1,p2-p1);
757         static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsAngleAtAbsEnd))->Fill(aAbs,p2-p1);
758         static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsMCAngle))->Fill(aMC,p2-p1);
759         static_cast<TH3*>(fTrackerList->At(kResPAtVtxVsAngleAtAbsEndVsP))->Fill(p1,aAbs,p2-p1);
760         static_cast<TH2*>(fTrackerList->At(kResPtAtVtxVsPt))->Fill(pT1,pT2-pT1);
761         
762         static_cast<TH1*>(fTrackerList->At(kResSlopeXAtVtx))->Fill(slopex2-slopex1);
763         static_cast<TH1*>(fTrackerList->At(kResSlopeYAtVtx))->Fill(slopey2-slopey1);
764         static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsP))->Fill(p1,slopex2-slopex1);
765         static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsP))->Fill(p1,slopey2-slopey1);
766         static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsAngleAtAbsEnd))->Fill(aAbs,slopex2-slopex1);
767         static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsAngleAtAbsEnd))->Fill(aAbs,slopey2-slopey1);
768         static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsMCAngle))->Fill(aMC,slopex2-slopex1);
769         static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsMCAngle))->Fill(aMC,slopey2-slopey1);
770         
771         static_cast<TH1*>(fTrackerList->At(kResEtaAtVtx))->Fill(eta2-eta1);
772         static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsP))->Fill(p1,eta2-eta1);
773         static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsAngleAtAbsEnd))->Fill(aAbs,eta2-eta1);
774         static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsMCAngle))->Fill(aMC,eta2-eta1);
775         
776         static_cast<TH1*>(fTrackerList->At(kResPhiAtVtx))->Fill(dPhi);
777         static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsP))->Fill(p1,dPhi);
778         static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsAngleAtAbsEnd))->Fill(aAbs,dPhi);
779         static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsMCAngle))->Fill(aMC,dPhi);
780         
781         static_cast<TH1*>(fTrackerList->At(kPDCA))->Fill(0.5*(p2+pU)*dca);
782         static_cast<TH2*>(fTrackerList->At(kPDCAVsAngleAtAbsEnd))->Fill(aAbs,0.5*(p2+pU)*dca);
783         static_cast<TH2*>(fTrackerList->At(kPDCAVsMCAngle))->Fill(aMC,0.5*(p2+pU)*dca);
784         
785         if (aAbs > 2. && aAbs < 3.) {
786           
787           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPIn23deg))->Fill(p1,p2-p1);
788           static_cast<TH2*>(fTrackerList->At(kPDCAVsPIn23deg))->Fill(p1,0.5*(p2+pU)*dca);
789           static_cast<TH2*>(fTrackerList->At(kPMCSAngVsPIn23deg))->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
790           
791         } else if (aAbs >= 3. && aAbs < 10.) {
792           
793           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPIn310deg))->Fill(p1,p2-p1);
794           static_cast<TH2*>(fTrackerList->At(kPDCAVsPIn310deg))->Fill(p1,0.5*(p2+pU)*dca);
795           static_cast<TH2*>(fTrackerList->At(kPMCSAngVsPIn310deg))->Fill(p1,0.5*(p2+pU)*(aMCS-aMC)*TMath::DegToRad());
796           
797         }
798         
799         if (aMC < 2.) {
800           
801           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPIn02degMC))->Fill(p1,p2-p1);
802           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,p2-p1);
803           static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,slopex2-slopex1);
804           static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,slopey2-slopey1);
805           static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,eta2-eta1);
806           static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsPosAbsEndIn02degMC))->Fill(dAbs,dPhi);
807           static_cast<TH2*>(fTrackerList->At(kPDCAVsPosAbsEndIn02degMC))->Fill(dAbs,0.5*(p2+pU)*dca);
808           
809         } else if (aMC >= 2. && aMC < 3) {
810           
811           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,p2-p1);
812           static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,slopex2-slopex1);
813           static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,slopey2-slopey1);
814           static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,eta2-eta1);
815           static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsPosAbsEndIn23degMC))->Fill(dAbs,dPhi);
816           static_cast<TH2*>(fTrackerList->At(kPDCAVsPosAbsEndIn23degMC))->Fill(dAbs,0.5*(p2+pU)*dca);
817           
818         } else if (aMC >= 3. && aMC < 10.) {
819           
820           static_cast<TH2*>(fTrackerList->At(kResPAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,p2-p1);
821           static_cast<TH2*>(fTrackerList->At(kResSlopeXAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,slopex2-slopex1);
822           static_cast<TH2*>(fTrackerList->At(kResSlopeYAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,slopey2-slopey1);
823           static_cast<TH2*>(fTrackerList->At(kResEtaAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,eta2-eta1);
824           static_cast<TH2*>(fTrackerList->At(kResPhiAtVtxVsPosAbsEndIn310degMC))->Fill(dAbs,dPhi);
825           static_cast<TH2*>(fTrackerList->At(kPDCAVsPosAbsEndIn310degMC))->Fill(dAbs,0.5*(p2+pU)*dca);
826           
827         }
828         
829         // simulated track parameters at vertex
830         trackParam = (AliMUONTrackParam*) matchedTrackRef->GetTrackParamAtCluster()->First();
831         x1 = trackParam->GetNonBendingCoor();
832         y1 = trackParam->GetBendingCoor();
833         z1 = trackParam->GetZ();
834         slopex1 = trackParam->GetNonBendingSlope();
835         slopey1 = trackParam->GetBendingSlope();
836         pX1 = trackParam->Px();
837         pY1 = trackParam->Py();
838         pZ1 = trackParam->Pz();
839         p1  = trackParam->P();
840         pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1);
841         
842         // reconstructed track parameters at vertex
843         trackParam = (AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->First();
844         x2 = trackParam->GetNonBendingCoor();
845         y2 = trackParam->GetBendingCoor();
846         z2 = trackParam->GetZ();
847         slopex2 = trackParam->GetNonBendingSlope();
848         slopey2 = trackParam->GetBendingSlope();
849         pX2 = trackParam->Px();
850         pY2 = trackParam->Py();
851         pZ2 = trackParam->Pz();
852         p2  = trackParam->P();
853         pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2);
854         
855         // fill histograms
856         static_cast<TH1*>(fTrackerList->At(kResPAt1stCl))->Fill(p2-p1);
857         static_cast<TH2*>(fTrackerList->At(kResPAt1stClVsP))->Fill(p1,p2-p1);
858         static_cast<TH2*>(fTrackerList->At(kResPtAt1stClVsPt))->Fill(pT1,pT2-pT1);
859         static_cast<TH1*>(fTrackerList->At(kResSlopeXAt1stCl))->Fill(slopex2-slopex1);
860         static_cast<TH1*>(fTrackerList->At(kResSlopeYAt1stCl))->Fill(slopey2-slopey1);
861         static_cast<TH2*>(fTrackerList->At(kResSlopeXAt1stClVsP))->Fill(p1,slopex2-slopex1);
862         static_cast<TH2*>(fTrackerList->At(kResSlopeYAt1stClVsP))->Fill(p1,slopey2-slopey1);
863         
864         // clusters residuals
865         for (Int_t iCl1 = 0; iCl1 < muonTrack.GetNClusters(); iCl1++) {
866           
867           AliMUONVCluster* cluster1 = static_cast<AliMUONTrackParam*>(muonTrack.GetTrackParamAtCluster()->UncheckedAt(iCl1))->GetClusterPtr();
868           Int_t chId = cluster1->GetChamberId();
869           Int_t deId = cluster1->GetDetElemId();
870           
871           for (Int_t iCl2 = 0; iCl2 < matchedTrackRef->GetNClusters(); iCl2++) {
872             
873             AliMUONVCluster* cluster2 = static_cast<AliMUONTrackParam*>(matchedTrackRef->GetTrackParamAtCluster()->UncheckedAt(iCl2))->GetClusterPtr();
874             
875             if (cluster2->GetDetElemId() == deId) {
876               
877               // fill histograms
878               static_cast<TH2*>(fTrackerList->At(kResClXVsCh))->Fill(chId+1, cluster1->GetX()-cluster2->GetX());
879               static_cast<TH2*>(fTrackerList->At(kResClYVsCh))->Fill(chId+1, cluster1->GetY()-cluster2->GetY());
880               static_cast<TH2*>(fTrackerList->At(kResClXVsDE))->Fill(fDEIndices[deId], cluster1->GetX()-cluster2->GetX());
881               static_cast<TH2*>(fTrackerList->At(kResClYVsDE))->Fill(fDEIndices[deId], cluster1->GetY()-cluster2->GetY());
882               
883               break;
884               
885             }
886             
887           }
888           
889         }
890         
891       } // end if (matchedTrackRef)
892       
893     } // end if (esdTrack->ContainTrackerData())
894     
895     // look for MC trigger associated to MC track
896     AliMUONTriggerTrack* triggerTrackRef = (isValid && matchedTrackRef)
897     ? static_cast<AliMUONTriggerTrack*>(triggerTrackRefStore->FindObject(matchedTrackRef->GetUniqueID()))
898     : 0x0;
899     
900     // Trigger
901     AliMUONTriggerTrack* matchedTrigTrackRef = 0x0;
902     containerInput[kVarDupliTrg] = 0.;
903     if (esdTrack->ContainTriggerData()) {
904
905       // check if this track is not already accounted for and record it if not
906       Bool_t trackExist = kFALSE;
907       for (Int_t i=0; i<nTrgTracks; i++) if (esdTrack->LoCircuit() == loCircuit[i]) trackExist = kTRUE;
908       if (trackExist) containerInput[kVarDupliTrg] = 1.;
909       else loCircuit[nTrgTracks++] = esdTrack->LoCircuit();
910       
911       // Convert ESD track to trigger track
912       AliMUONLocalTrigger locTrg;
913       AliMUONESDInterface::ESDToMUON(*esdTrack, locTrg);
914       AliMUONTriggerTrack trigTrack;
915       rc.TriggerToTrack(locTrg, trigTrack);
916       
917       // try to match the trigger track with the same MC track if any
918       if (triggerTrackRef && trigTrack.Match(*triggerTrackRef, fSigmaCutTrig)) matchedTrigTrackRef = triggerTrackRef;
919       if (matchedTrigTrackRef) containerInput[kVarMatchMC] = static_cast<Double_t>(kMatchedSame);
920       else { // or try to match with any triggerable track
921         matchedTrigTrackRef = rc.FindCompatibleTrack(trigTrack, *triggerTrackRefStore, fSigmaCutTrig);
922         if (matchedTrigTrackRef) {
923           if (isValid && matchedTrackRef) {
924             containerInput[kVarMatchMC] = static_cast<Double_t>(kMatchedDiff);
925             if (!fMCTrigLevelFromMatchTrk) triggerTrackRef = matchedTrigTrackRef;
926           } else {
927             containerInput[kVarMatchMC] = static_cast<Double_t>(kTriggerOnly);
928             triggerTrackRef = matchedTrigTrackRef;
929           }
930         }
931       }
932       
933       // fill histograms
934       if (matchedTrigTrackRef) {
935         static_cast<TH1*>(fTriggerList->At(kResTrigX11))->Fill(trigTrack.GetX11() - matchedTrigTrackRef->GetX11());
936         static_cast<TH1*>(fTriggerList->At(kResTrigY11))->Fill(trigTrack.GetY11() - matchedTrigTrackRef->GetY11());
937         static_cast<TH1*>(fTriggerList->At(kResTrigSlopeY))->Fill(trigTrack.GetSlopeY() - matchedTrigTrackRef->GetSlopeY());
938       }
939       
940     }
941     
942     // fill container
943     if (triggerTrackRef) {
944       if (triggerTrackRef->GetPtCutLevel() == 0) containerInput[kVarMCTrigger] = static_cast<Double_t>(kOtherTrig);
945       else if (triggerTrackRef->GetPtCutLevel() == 1) containerInput[kVarMCTrigger] = static_cast<Double_t>(kAllPtTrig);
946       else if (triggerTrackRef->GetPtCutLevel() == 2) containerInput[kVarMCTrigger] = static_cast<Double_t>(kLowPtTrig);
947       else if (triggerTrackRef->GetPtCutLevel() == 3) containerInput[kVarMCTrigger] = static_cast<Double_t>(kHighPtTrig);
948     } else containerInput[kVarMCTrigger] = kNoMatchTrig;
949     
950     // get MC particle ID
951     Int_t mcID = -1;
952     if (isValid && matchedTrackRef) mcID = static_cast<Int_t>(matchedTrackRef->GetUniqueID());
953     else if (matchedTrigTrackRef) mcID = static_cast<Int_t>(matchedTrigTrackRef->GetUniqueID());
954     
955     // fill particle container
956     FillContainerInfoReco(containerInput, esdTrack, isValid, mcID);
957     fCFContainer->Fill(containerInput, kStepReconstructed);
958
959   }
960
961   // clean memory
962   delete[] loCircuit;
963   containerInput[kVarDupliTrg] = 0.;
964   
965   // ------ Loop over reconstructible tracks ------
966   AliMUONTrack* trackRef = 0x0;
967   TIter next(reconstructibleStore->CreateIterator());
968   while ((trackRef = static_cast<AliMUONTrack*>(next()))) {
969     
970     // find the corresponding triggerable track if any
971     UInt_t mcID = trackRef->GetUniqueID();
972     AliMUONTriggerTrack* trigTrackRef = static_cast<AliMUONTriggerTrack*>(triggerTrackRefStore->FindObject(mcID));
973     
974     // fill particle container
975     FillContainerInfoMC(containerInput, static_cast<AliMCParticle*>(fMCEvent->GetTrack(static_cast<Int_t>(mcID))));
976     containerInput[kVarHasTracker] = 1.;
977     if (trigTrackRef) {
978       if (trigTrackRef->GetPtCutLevel() == 0) containerInput[kVarTrigger] = static_cast<Double_t>(kOtherTrig);
979       else if (trigTrackRef->GetPtCutLevel() == 1) containerInput[kVarTrigger] = static_cast<Double_t>(kAllPtTrig);
980       else if (trigTrackRef->GetPtCutLevel() == 2) containerInput[kVarTrigger] = static_cast<Double_t>(kLowPtTrig);
981       else if (trigTrackRef->GetPtCutLevel() == 3) containerInput[kVarTrigger] = static_cast<Double_t>(kHighPtTrig);
982     } else containerInput[kVarTrigger] = static_cast<Double_t>(kNoMatchTrig);
983     containerInput[kVarMatchMC] = static_cast<Double_t>(kNoMatch);
984     containerInput[kVarMCTrigger] = static_cast<Double_t>(kNoMatchTrig);
985     fCFContainer->Fill(containerInput, kStepGeneratedMC);
986     
987   }
988
989   // ------ Loop over triggerable ghosts ------
990   AliMUONTriggerTrack* trigTrackRef = 0x0;
991   TIter nextTrig(triggerTrackRefStore->CreateIterator());
992   while ((trigTrackRef = static_cast<AliMUONTriggerTrack*>(nextTrig()))) {
993     
994     // discard tracks also reconstructible
995     UInt_t mcID = trigTrackRef->GetUniqueID();
996     if (reconstructibleStore->FindObject(mcID)) continue;
997
998     // fill particle container
999     FillContainerInfoMC(containerInput, static_cast<AliMCParticle*>(fMCEvent->GetTrack(static_cast<Int_t>(mcID))));
1000     containerInput[kVarHasTracker] = 0.;
1001     if (trigTrackRef->GetPtCutLevel() == 0) containerInput[kVarTrigger] = static_cast<Double_t>(kOtherTrig);
1002     else if (trigTrackRef->GetPtCutLevel() == 1) containerInput[kVarTrigger] = static_cast<Double_t>(kAllPtTrig);
1003     else if (trigTrackRef->GetPtCutLevel() == 2) containerInput[kVarTrigger] = static_cast<Double_t>(kLowPtTrig);
1004     else if (trigTrackRef->GetPtCutLevel() == 3) containerInput[kVarTrigger] = static_cast<Double_t>(kHighPtTrig);
1005     containerInput[kVarMatchMC] = static_cast<Double_t>(kNoMatch);
1006     containerInput[kVarMCTrigger] = static_cast<Double_t>(kNoMatchTrig);
1007     fCFContainer->Fill(containerInput, kStepGeneratedMC);
1008     
1009   }
1010   
1011   // Post final data
1012   PostData(1, fCFContainer);
1013   PostData(2, fTriggerList);
1014   PostData(3, fTrackerList);
1015 }
1016
1017 //________________________________________________________________________
1018 void AliAnalysisTaskMuonPerformance::Terminate(Option_t *)
1019 {
1020   //
1021   /// Draw some histogram at the end.
1022   //
1023
1024   // do it only locally
1025   //if (gROOT->IsBatch()) return;
1026   
1027   // get output containers
1028   fCFContainer = dynamic_cast<AliCFContainer*>(GetOutputData(1));
1029   fTriggerList = dynamic_cast<TObjArray*>(GetOutputData(2));
1030   fTrackerList = dynamic_cast<TObjArray*>(GetOutputData(3));
1031   if (!fCFContainer || !fTriggerList || !fTrackerList) {
1032     AliWarning("Output containers not found: summary histograms are not created"); 
1033     return;
1034   }
1035   
1036   // --- compute efficiencies ---
1037   fEfficiencyList = new TObjArray(100);
1038   fEfficiencyList->SetOwner();
1039   
1040   TObjArray* effAnyPt = new TObjArray(100);
1041   effAnyPt->SetName("effAnyPt");
1042   effAnyPt->SetOwner();
1043   fEfficiencyList->AddLast(effAnyPt);
1044   
1045   TObjArray* effAllPt = new TObjArray(100);
1046   effAllPt->SetName("effAllPt");
1047   effAllPt->SetOwner();
1048   fEfficiencyList->AddLast(effAllPt);
1049   
1050   TObjArray* effLowPt = new TObjArray(100);
1051   effLowPt->SetName("effLowPt");
1052   effLowPt->SetOwner();
1053   fEfficiencyList->AddLast(effLowPt);
1054   
1055   TObjArray* effHighPt = new TObjArray(100);
1056   effHighPt->SetName("effHighPt");
1057   effHighPt->SetOwner();
1058   fEfficiencyList->AddLast(effHighPt);
1059   
1060   TObjArray* notTrgable = new TObjArray(100);
1061   notTrgable->SetName("notTrgable");
1062   notTrgable->SetOwner();
1063   fEfficiencyList->AddLast(notTrgable);
1064   
1065   TObjArray* trgableNoPtOnly = new TObjArray(100);
1066   trgableNoPtOnly->SetName("trgableNoPtOnly");
1067   trgableNoPtOnly->SetOwner();
1068   fEfficiencyList->AddLast(trgableNoPtOnly);
1069   
1070   TObjArray* trgableAPtOnly = new TObjArray(100);
1071   trgableAPtOnly->SetName("trgableAPtOnly");
1072   trgableAPtOnly->SetOwner();
1073   fEfficiencyList->AddLast(trgableAPtOnly);
1074   
1075   TObjArray* trgableLPtOnly = new TObjArray(100);
1076   trgableLPtOnly->SetName("trgableLPtOnly");
1077   trgableLPtOnly->SetOwner();
1078   fEfficiencyList->AddLast(trgableLPtOnly);
1079   
1080   TObjArray* trgableHPtOnly = new TObjArray(100);
1081   trgableHPtOnly->SetName("trgableHPtOnly");
1082   trgableHPtOnly->SetOwner();
1083   fEfficiencyList->AddLast(trgableHPtOnly);
1084   
1085   AliCFEffGrid* efficiency = new AliCFEffGrid("eff","",*fCFContainer);
1086   efficiency->CalculateEfficiency(kStepReconstructed, kStepGeneratedMC);
1087   Double_t totalEff = 0., totalEffErr = 0.;
1088   Int_t sumEffBin = 0;
1089   
1090   // add histogram summarizing global efficiencies
1091   TH1D* effSummary = new TH1D("effSummary", "Efficiency summary", 1, 0., 0.);
1092   effSummary->GetYaxis()->SetTitle("Efficiency");
1093   fEfficiencyList->AddLast(effSummary);
1094   
1095   // ------ Tracker only ------
1096   
1097   // Tracker efficiency using all reconstructed tracks
1098   efficiency->GetNum()->SetRangeUser(kVarHasTracker, 1., 1.);
1099   efficiency->GetDen()->SetRangeUser(kVarHasTracker, 1., 1.);
1100   efficiency->GetNum()->GetAxis(kVarTrigger)->SetRange();
1101   efficiency->GetDen()->GetAxis(kVarTrigger)->SetRange();
1102   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1103   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1104   FillEffHistos(efficiency, "trackerTracks", fEfficiencyList);
1105   GetEfficiency(efficiency, totalEff, totalEffErr);
1106   effSummary->Fill("Tracker_all", totalEff);
1107   effSummary->SetBinError(++sumEffBin, totalEffErr);
1108   printf("Tracker efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1109   
1110   // Tracker efficiency using tracks matched with reconstructible ones
1111   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kMatchedDiff);
1112   FillEffHistos(efficiency, "trackerTracksMatchMC", fEfficiencyList);
1113   GetEfficiency(efficiency, totalEff, totalEffErr);
1114   effSummary->Fill("Tracker_MCId", totalEff);
1115   effSummary->SetBinError(++sumEffBin, totalEffErr);
1116   printf("Tracker efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1117   
1118   // ------ Tracker matched with any trigger ------
1119   
1120   // Matched efficiency using all reconstructed tracks
1121   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kHighPtTrig);
1122   efficiency->GetDen()->SetRangeUser(kVarTrigger, kOtherTrig, kHighPtTrig);
1123   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1124   FillEffHistos(efficiency, "matchedTracks", effAnyPt);
1125   GetEfficiency(efficiency, totalEff, totalEffErr);
1126   effSummary->Fill("Matched_all", totalEff);
1127   effSummary->SetBinError(++sumEffBin, totalEffErr);
1128   printf("Matched efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1129   
1130   // Matched efficiency using tracks matched with reconstructible ones, triggerable or not
1131   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kMatchedDiff);
1132   FillEffHistos(efficiency, "matchedTracksMatchMC", effAnyPt);
1133   GetEfficiency(efficiency, totalEff, totalEffErr);
1134   effSummary->Fill("Matched_MCId", totalEff);
1135   effSummary->SetBinError(++sumEffBin, totalEffErr);
1136   printf("Matched efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1137   
1138   // Matched efficiency using tracks matched with reconstructible ones triggerable
1139   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kHighPtTrig);
1140   FillEffHistos(efficiency, "matchedTracksMatchMCAnypt", effAnyPt);
1141   GetEfficiency(efficiency, totalEff, totalEffErr);
1142   effSummary->Fill("Matched_MCIdAnypt", totalEff);
1143   effSummary->SetBinError(++sumEffBin, totalEffErr);
1144   printf("Matched efficiency using reconstructed tracks matching MC-anyPt = %f +- %f\n", totalEff, totalEffErr);
1145   
1146   // Matched efficiency using tracks matched with reconstructible ones not triggerable
1147   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kNoMatchTrig);
1148   FillEffHistos(efficiency, "matchedTracksMatchMCNoTrig", effAnyPt);
1149   GetEfficiency(efficiency, totalEff, totalEffErr);
1150   effSummary->Fill("Matched_MCIdNoTrig", totalEff);
1151   effSummary->SetBinError(++sumEffBin, totalEffErr);
1152   printf("Matched efficiency using reconstructed tracks matching MC-noTrig = %f +- %f\n", totalEff, totalEffErr);
1153   
1154   // Matched efficiency using tracks matched with same reconstructible & triggerable ones
1155   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kMatchedSame);
1156   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1157   FillEffHistos(efficiency, "matchedTracksMatchSameMC", effAnyPt);
1158   GetEfficiency(efficiency, totalEff, totalEffErr);
1159   effSummary->Fill("Matched_SameMCId", totalEff);
1160   effSummary->SetBinError(++sumEffBin, totalEffErr);
1161   printf("Matched efficiency using reconstructed tracks matching same MC = %f +- %f\n", totalEff, totalEffErr);
1162   
1163   // Matched efficiency using tracks matched with different reconstructible & triggerable ones
1164   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedDiff, kMatchedDiff);
1165   FillEffHistos(efficiency, "matchedTracksMatchDiffMC", effAnyPt);
1166   GetEfficiency(efficiency, totalEff, totalEffErr);
1167   effSummary->Fill("Matched_DiffMCId", totalEff);
1168   effSummary->SetBinError(++sumEffBin, totalEffErr);
1169   printf("Matched efficiency using reconstructed tracks matching different MC = %f +- %f\n", totalEff, totalEffErr);
1170   
1171   // Matched efficiency using tracks matched with reconstructible ones only
1172   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kTrackerOnly);
1173   FillEffHistos(efficiency, "matchedTracksMatchTrkMC", effAnyPt);
1174   GetEfficiency(efficiency, totalEff, totalEffErr);
1175   effSummary->Fill("Matched_TrkMCId", totalEff);
1176   effSummary->SetBinError(++sumEffBin, totalEffErr);
1177   printf("Matched efficiency using reconstructed tracks matching tracker MC = %f +- %f\n", totalEff, totalEffErr);
1178   
1179   // ------ Tracker matched with all pt trigger ------
1180   
1181   // Matched all pt efficiency using all reconstructed tracks
1182   efficiency->GetDen()->SetRangeUser(kVarTrigger, kAllPtTrig, kHighPtTrig);
1183   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1184   FillEffHistos(efficiency, "matchedTracks", effAllPt);
1185   GetEfficiency(efficiency, totalEff, totalEffErr);
1186   printf("Matched Apt efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1187   
1188   // Matched all pt efficiency using tracks matched with reconstructible ones, triggerable or not
1189   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kMatchedDiff);
1190   FillEffHistos(efficiency, "matchedTracksMatchMC", effAllPt);
1191   GetEfficiency(efficiency, totalEff, totalEffErr);
1192   printf("Matched Apt efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1193   
1194   // Matched all pt efficiency using tracks matched with reconstructible ones triggerable Apt
1195   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kAllPtTrig, kHighPtTrig);
1196   FillEffHistos(efficiency, "matchedTracksMatchMCApt", effAllPt);
1197   GetEfficiency(efficiency, totalEff, totalEffErr);
1198   printf("Matched Apt efficiency using reconstructed tracks matching MC-Apt = %f +- %f\n", totalEff, totalEffErr);
1199   
1200   // Matched all pt efficiency using tracks matched with reconstructible ones triggerable other pt
1201   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kOtherTrig);
1202   FillEffHistos(efficiency, "matchedTracksMatchMCOther", effAllPt);
1203   GetEfficiency(efficiency, totalEff, totalEffErr);
1204   printf("Matched Apt efficiency using reconstructed tracks matching MC-other = %f +- %f\n", totalEff, totalEffErr);
1205   
1206   // Matched all pt efficiency using tracks matched with reconstructible ones not triggerable
1207   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kNoMatchTrig);
1208   FillEffHistos(efficiency, "matchedTracksMatchMCNoTrig", effAllPt);
1209   GetEfficiency(efficiency, totalEff, totalEffErr);
1210   printf("Matched Apt efficiency using reconstructed tracks matching MC-noTrig = %f +- %f\n", totalEff, totalEffErr);
1211   
1212   // Matched all pt efficiency using tracks matched with same reconstructible & triggerable ones (all pt MC trig)
1213   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kMatchedSame);
1214   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kAllPtTrig, kHighPtTrig);
1215   FillEffHistos(efficiency, "matchedTracksMatchSameMCApt", effAllPt);
1216   
1217   // Matched all pt efficiency using tracks matched with same reconstructible & triggerable ones (other MC trig)
1218   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kOtherTrig);
1219   FillEffHistos(efficiency, "matchedTracksMatchSameMCOther", effAllPt);
1220   
1221   // Matched all pt efficiency using tracks matched with different reconstructible & triggerable ones (all pt MC trig)
1222   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedDiff, kMatchedDiff);
1223   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kAllPtTrig, kHighPtTrig);
1224   FillEffHistos(efficiency, "matchedTracksMatchDiffMCApt", effAllPt);
1225   
1226   // Matched all pt efficiency using tracks matched with different reconstructible & triggerable ones (other MC trig)
1227   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kOtherTrig);
1228   FillEffHistos(efficiency, "matchedTracksMatchDiffMCOther", effAllPt);
1229   
1230   // Matched efficiency using tracks matched with reconstructible ones only (all pt MC trig)
1231   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kTrackerOnly);
1232   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kAllPtTrig, kHighPtTrig);
1233   FillEffHistos(efficiency, "matchedTracksMatchTrkMCApt", effAllPt);
1234   
1235   // Matched efficiency using tracks matched with reconstructible ones only (other MC trig)
1236   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kOtherTrig);
1237   FillEffHistos(efficiency, "matchedTracksMatchTrkMCOther", effAllPt);
1238   
1239   // ------ Tracker matched with low pt trigger ------
1240   
1241   // Matched low pt efficiency using all reconstructed tracks
1242   efficiency->GetNum()->SetRangeUser(kVarTrigger, kLowPtTrig, kHighPtTrig);
1243   efficiency->GetDen()->SetRangeUser(kVarTrigger, kLowPtTrig, kHighPtTrig);
1244   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1245   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1246   FillEffHistos(efficiency, "matchedTracks", effLowPt);
1247   GetEfficiency(efficiency, totalEff, totalEffErr);
1248   printf("Matched Lpt efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1249   
1250   // Matched low pt efficiency using tracks matched with reconstructible ones, triggerable or not
1251   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kMatchedDiff);
1252   FillEffHistos(efficiency, "matchedTracksMatchMC", effLowPt);
1253   GetEfficiency(efficiency, totalEff, totalEffErr);
1254   printf("Matched Lpt efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1255   
1256   // Matched low pt efficiency using tracks matched with reconstructible ones triggerable Lpt
1257   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kLowPtTrig, kHighPtTrig);
1258   FillEffHistos(efficiency, "matchedTracksMatchMCLpt", effLowPt);
1259   GetEfficiency(efficiency, totalEff, totalEffErr);
1260   printf("Matched Lpt efficiency using reconstructed tracks matching MC-Lpt = %f +- %f\n", totalEff, totalEffErr);
1261   
1262   // Matched low pt efficiency using tracks matched with reconstructible ones triggerable other pt
1263   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kAllPtTrig);
1264   FillEffHistos(efficiency, "matchedTracksMatchMCOther", effLowPt);
1265   GetEfficiency(efficiency, totalEff, totalEffErr);
1266   printf("Matched Lpt efficiency using reconstructed tracks matching MC-other = %f +- %f\n", totalEff, totalEffErr);
1267   
1268   // Matched low pt efficiency using tracks matched with reconstructible ones not triggerable
1269   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kNoMatchTrig);
1270   FillEffHistos(efficiency, "matchedTracksMatchMCNoTrig", effLowPt);
1271   GetEfficiency(efficiency, totalEff, totalEffErr);
1272   printf("Matched Lpt efficiency using reconstructed tracks matching MC-noTrig = %f +- %f\n", totalEff, totalEffErr);
1273   
1274   // Matched low pt efficiency using tracks matched with same reconstructible & triggerable ones (low pt MC trig)
1275   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kMatchedSame);
1276   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kLowPtTrig, kHighPtTrig);
1277   FillEffHistos(efficiency, "matchedTracksMatchSameMCLpt", effLowPt);
1278   
1279   // Matched low pt efficiency using tracks matched with same reconstructible & triggerable ones (other MC trig)
1280   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kAllPtTrig);
1281   FillEffHistos(efficiency, "matchedTracksMatchSameMCOther", effLowPt);
1282   
1283   // Matched low pt efficiency using tracks matched with different reconstructible & triggerable ones (low pt MC trig)
1284   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedDiff, kMatchedDiff);
1285   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kLowPtTrig, kHighPtTrig);
1286   FillEffHistos(efficiency, "matchedTracksMatchDiffMCLpt", effLowPt);
1287   
1288   // Matched low pt efficiency using tracks matched with different reconstructible & triggerable ones (other MC trig)
1289   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kAllPtTrig);
1290   FillEffHistos(efficiency, "matchedTracksMatchDiffMCOther", effLowPt);
1291   
1292   // Matched efficiency using tracks matched with reconstructible ones only (low pt MC trig)
1293   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kTrackerOnly);
1294   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kLowPtTrig, kHighPtTrig);
1295   FillEffHistos(efficiency, "matchedTracksMatchTrkMCLpt", effLowPt);
1296   
1297   // Matched efficiency using tracks matched with reconstructible ones only (other MC trig)
1298   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kAllPtTrig);
1299   FillEffHistos(efficiency, "matchedTracksMatchTrkMCOther", effLowPt);
1300   
1301   // ------ Tracker matched with high pt trigger ------
1302   
1303   // Matched high pt efficiency using all reconstructed tracks
1304   efficiency->GetNum()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1305   efficiency->GetDen()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1306   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1307   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1308   FillEffHistos(efficiency, "matchedTracks", effHighPt);
1309   GetEfficiency(efficiency, totalEff, totalEffErr);
1310   printf("Matched Hpt efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1311   
1312   // Matched high pt efficiency using tracks matched with reconstructible ones, triggerable or not
1313   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kMatchedDiff);
1314   FillEffHistos(efficiency, "matchedTracksMatchMC", effHighPt);
1315   GetEfficiency(efficiency, totalEff, totalEffErr);
1316   printf("Matched Hpt efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1317   
1318   // Matched high pt efficiency using tracks matched with reconstructible ones triggerable Hpt
1319   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kHighPtTrig, kHighPtTrig);
1320   FillEffHistos(efficiency, "matchedTracksMatchMCHpt", effHighPt);
1321   GetEfficiency(efficiency, totalEff, totalEffErr);
1322   printf("Matched Hpt efficiency using reconstructed tracks matching MC-Hpt = %f +- %f\n", totalEff, totalEffErr);
1323   
1324   // Matched high pt efficiency using tracks matched with reconstructible ones triggerable other pt
1325   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kLowPtTrig);
1326   FillEffHistos(efficiency, "matchedTracksMatchMCOther", effHighPt);
1327   GetEfficiency(efficiency, totalEff, totalEffErr);
1328   printf("Matched Hpt efficiency using reconstructed tracks matching MC-other = %f +- %f\n", totalEff, totalEffErr);
1329   
1330   // Matched high pt efficiency using tracks matched with reconstructible ones not triggerable
1331   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kNoMatchTrig);
1332   FillEffHistos(efficiency, "matchedTracksMatchMCNoTrig", effHighPt);
1333   GetEfficiency(efficiency, totalEff, totalEffErr);
1334   printf("Matched Hpt efficiency using reconstructed tracks matching MC-noTrig = %f +- %f\n", totalEff, totalEffErr);
1335   
1336   // Matched high pt efficiency using tracks matched with same reconstructible & triggerable ones (high pt MC trig)
1337   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kMatchedSame);
1338   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kHighPtTrig, kHighPtTrig);
1339   FillEffHistos(efficiency, "matchedTracksMatchSameMCHpt", effHighPt);
1340   
1341   // Matched high pt efficiency using tracks matched with same reconstructible & triggerable ones (other MC trig)
1342   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kLowPtTrig);
1343   FillEffHistos(efficiency, "matchedTracksMatchSameMCOther", effHighPt);
1344   
1345   // Matched high pt efficiency using tracks matched with different reconstructible & triggerable ones (high pt MC trig)
1346   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedDiff, kMatchedDiff);
1347   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kHighPtTrig, kHighPtTrig);
1348   FillEffHistos(efficiency, "matchedTracksMatchDiffMCHpt", effHighPt);
1349   
1350   // Matched high pt efficiency using tracks matched with different reconstructible & triggerable ones (other MC trig)
1351   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kLowPtTrig);
1352   FillEffHistos(efficiency, "matchedTracksMatchDiffMCOther", effHighPt);
1353   
1354   // Matched efficiency using tracks matched with reconstructible ones only (high pt MC trig)
1355   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kTrackerOnly);
1356   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kHighPtTrig, kHighPtTrig);
1357   FillEffHistos(efficiency, "matchedTracksMatchTrkMCHpt", effHighPt);
1358   
1359   // Matched efficiency using tracks matched with reconstructible ones only (other MC trig)
1360   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kLowPtTrig);
1361   FillEffHistos(efficiency, "matchedTracksMatchTrkMCOther", effHighPt);
1362   
1363   // ------ Trigger only ------
1364   
1365   // Trigger efficiency using all reconstructed tracks
1366   efficiency->GetNum()->GetAxis(kVarHasTracker)->SetRange();
1367   efficiency->GetDen()->GetAxis(kVarHasTracker)->SetRange();
1368   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kHighPtTrig);
1369   efficiency->GetDen()->SetRangeUser(kVarTrigger, kOtherTrig, kHighPtTrig);
1370   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1371   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1372   efficiency->GetNum()->SetRangeUser(kVarDupliTrg, 0., 0.);
1373   FillEffHistos(efficiency, "triggerTracks", effAnyPt);
1374   GetEfficiency(efficiency, totalEff, totalEffErr);
1375   effSummary->Fill("Trigger_all", totalEff);
1376   effSummary->SetBinError(++sumEffBin, totalEffErr);
1377   printf("Trigger efficiency using all reconstructed tracks = %f +- %f\n", totalEff, totalEffErr);
1378   
1379   // Trigger efficiency using tracks matched with triggerable ones
1380   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kTriggerOnly);
1381   FillEffHistos(efficiency, "triggerTracksMatchMC", effAnyPt);
1382   GetEfficiency(efficiency, totalEff, totalEffErr);
1383   effSummary->Fill("Trigger_MCId", totalEff);
1384   effSummary->SetBinError(++sumEffBin, totalEffErr);
1385   printf("Trigger efficiency using reconstructed tracks matching MC = %f +- %f\n", totalEff, totalEffErr);
1386   
1387   // ------ All pt trigger only ------
1388   
1389   // All pt trigger efficiency using all reconstructed tracks
1390   efficiency->GetDen()->SetRangeUser(kVarTrigger, kAllPtTrig, kHighPtTrig);
1391   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1392   FillEffHistos(efficiency, "triggerTracks", effAllPt);
1393   
1394   // All pt trigger efficiency using tracks matched with all pt triggerable ones
1395   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kTriggerOnly);
1396   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kAllPtTrig, kHighPtTrig);
1397   FillEffHistos(efficiency, "triggerTracksMatchMCApt", effAllPt);
1398   
1399   // All pt trigger efficiency using tracks matched with other triggerable ones
1400   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kOtherTrig);
1401   FillEffHistos(efficiency, "triggerTracksMatchMCOther", effAllPt);
1402   
1403   // ------ Low pt trigger only ------
1404   
1405   // Low pt trigger efficiency using all reconstructed tracks
1406   efficiency->GetNum()->SetRangeUser(kVarTrigger, kLowPtTrig, kHighPtTrig);
1407   efficiency->GetDen()->SetRangeUser(kVarTrigger, kLowPtTrig, kHighPtTrig);
1408   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1409   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1410   FillEffHistos(efficiency, "triggerTracks", effLowPt);
1411   
1412   // Low pt trigger efficiency using tracks matched with Low pt triggerable ones
1413   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kTriggerOnly);
1414   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kLowPtTrig, kHighPtTrig);
1415   FillEffHistos(efficiency, "triggerTracksMatchMCLpt", effLowPt);
1416   
1417   // Low pt trigger efficiency using tracks matched with other triggerable ones
1418   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kAllPtTrig);
1419   FillEffHistos(efficiency, "triggerTracksMatchMCOther", effLowPt);
1420   
1421   // ------ High pt trigger only ------
1422   
1423   // High pt trigger efficiency using all reconstructed tracks
1424   efficiency->GetNum()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1425   efficiency->GetDen()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1426   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1427   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1428   FillEffHistos(efficiency, "triggerTracks", effHighPt);
1429   
1430   // High pt trigger efficiency using tracks matched with High pt triggerable ones
1431   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kMatchedSame, kTriggerOnly);
1432   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kHighPtTrig, kHighPtTrig);
1433   FillEffHistos(efficiency, "triggerTracksMatchMCHpt", effHighPt);
1434   
1435   // All pt trigger efficiency using tracks matched with other triggerable ones
1436   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kLowPtTrig);
1437   FillEffHistos(efficiency, "triggerTracksMatchMCOther", effHighPt);
1438   
1439   // ------ Tracker reconstructible not triggerable ------
1440   
1441   // all tracker tracks
1442   efficiency->GetNum()->SetRangeUser(kVarHasTracker, 1., 1.);
1443   efficiency->GetDen()->SetRangeUser(kVarHasTracker, 1., 1.);
1444   efficiency->GetNum()->GetAxis(kVarTrigger)->SetRange();
1445   efficiency->GetDen()->SetRangeUser(kVarTrigger, kNoMatchTrig, kNoMatchTrig);
1446   efficiency->GetNum()->SetRangeUser(kVarMatchMC, kTrackerOnly, kMatchedDiff);
1447   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kNoMatchTrig, kNoMatchTrig);
1448   efficiency->GetNum()->GetAxis(kVarDupliTrg)->SetRange();
1449   FillEffHistos(efficiency, "allTracks", notTrgable);
1450   
1451   // tracker not matched
1452   efficiency->GetNum()->SetRangeUser(kVarTrigger, kNoMatchTrig, kNoMatchTrig);
1453   FillEffHistos(efficiency, "notMatched", notTrgable);
1454   
1455   // tracker matched with all pt trigger
1456   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kAllPtTrig);
1457   FillEffHistos(efficiency, "MatchedApt", notTrgable);
1458   
1459   // tracker matched with low pt trigger
1460   efficiency->GetNum()->SetRangeUser(kVarTrigger, kLowPtTrig, kLowPtTrig);
1461   FillEffHistos(efficiency, "MatchedLpt", notTrgable);
1462   
1463   // tracker matched with high pt trigger
1464   efficiency->GetNum()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1465   FillEffHistos(efficiency, "MatchedHpt", notTrgable);
1466   
1467   // ------ Tracker reconstructible triggerable no pt ------
1468   
1469   // all tracker tracks
1470   efficiency->GetNum()->GetAxis(kVarTrigger)->SetRange();
1471   efficiency->GetDen()->SetRangeUser(kVarTrigger, kOtherTrig, kOtherTrig);
1472   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kOtherTrig, kOtherTrig);
1473   FillEffHistos(efficiency, "allTracks", trgableNoPtOnly);
1474   
1475   // tracker not matched
1476   efficiency->GetNum()->SetRangeUser(kVarTrigger, kNoMatchTrig, kNoMatchTrig);
1477   FillEffHistos(efficiency, "notMatched", trgableNoPtOnly);
1478   
1479   // tracker matched with all pt trigger
1480   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kAllPtTrig);
1481   FillEffHistos(efficiency, "MatchedApt", trgableNoPtOnly);
1482   
1483   // tracker matched with low pt trigger
1484   efficiency->GetNum()->SetRangeUser(kVarTrigger, kLowPtTrig, kLowPtTrig);
1485   FillEffHistos(efficiency, "MatchedLpt", trgableNoPtOnly);
1486   
1487   // tracker matched with high pt trigger
1488   efficiency->GetNum()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1489   FillEffHistos(efficiency, "MatchedHpt", trgableNoPtOnly);
1490   
1491   // ------ Tracker reconstructible triggerable Apt ------
1492   
1493   // all tracker tracks
1494   efficiency->GetNum()->GetAxis(kVarTrigger)->SetRange();
1495   efficiency->GetDen()->SetRangeUser(kVarTrigger, kAllPtTrig, kAllPtTrig);
1496   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kAllPtTrig, kAllPtTrig);
1497   FillEffHistos(efficiency, "allTracks", trgableAPtOnly);
1498   
1499   // tracker not matched
1500   efficiency->GetNum()->SetRangeUser(kVarTrigger, kNoMatchTrig, kNoMatchTrig);
1501   FillEffHistos(efficiency, "notMatched", trgableAPtOnly);
1502   
1503   // tracker matched with all pt trigger
1504   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kAllPtTrig);
1505   FillEffHistos(efficiency, "MatchedApt", trgableAPtOnly);
1506   
1507   // tracker matched with low pt trigger
1508   efficiency->GetNum()->SetRangeUser(kVarTrigger, kLowPtTrig, kLowPtTrig);
1509   FillEffHistos(efficiency, "MatchedLpt", trgableAPtOnly);
1510   
1511   // tracker matched with high pt trigger
1512   efficiency->GetNum()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1513   FillEffHistos(efficiency, "MatchedHpt", trgableAPtOnly);
1514   
1515   // ------ Tracker reconstructible triggerable Lpt ------
1516   
1517   // all tracker tracks
1518   efficiency->GetNum()->GetAxis(kVarTrigger)->SetRange();
1519   efficiency->GetDen()->SetRangeUser(kVarTrigger, kLowPtTrig, kLowPtTrig);
1520   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kLowPtTrig, kLowPtTrig);
1521   FillEffHistos(efficiency, "allTracks", trgableLPtOnly);
1522   
1523   // tracker not matched
1524   efficiency->GetNum()->SetRangeUser(kVarTrigger, kNoMatchTrig, kNoMatchTrig);
1525   FillEffHistos(efficiency, "notMatched", trgableLPtOnly);
1526   
1527   // tracker matched with all pt trigger
1528   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kAllPtTrig);
1529   FillEffHistos(efficiency, "MatchedApt", trgableLPtOnly);
1530   
1531   // tracker matched with low pt trigger
1532   efficiency->GetNum()->SetRangeUser(kVarTrigger, kLowPtTrig, kLowPtTrig);
1533   FillEffHistos(efficiency, "MatchedLpt", trgableLPtOnly);
1534   
1535   // tracker matched with high pt trigger
1536   efficiency->GetNum()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1537   FillEffHistos(efficiency, "MatchedHpt", trgableLPtOnly);
1538   
1539   // ------ Tracker reconstructible triggerable Hpt ------
1540   
1541   // all tracker tracks
1542   efficiency->GetNum()->GetAxis(kVarTrigger)->SetRange();
1543   efficiency->GetDen()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1544   efficiency->GetNum()->SetRangeUser(kVarMCTrigger, kHighPtTrig, kHighPtTrig);
1545   FillEffHistos(efficiency, "allTracks", trgableHPtOnly);
1546   
1547   // tracker not matched
1548   efficiency->GetNum()->SetRangeUser(kVarTrigger, kNoMatchTrig, kNoMatchTrig);
1549   FillEffHistos(efficiency, "notMatched", trgableHPtOnly);
1550   
1551   // tracker matched with all pt trigger
1552   efficiency->GetNum()->SetRangeUser(kVarTrigger, kAllPtTrig, kAllPtTrig);
1553   FillEffHistos(efficiency, "MatchedApt", trgableHPtOnly);
1554   
1555   // tracker matched with low pt trigger
1556   efficiency->GetNum()->SetRangeUser(kVarTrigger, kLowPtTrig, kLowPtTrig);
1557   FillEffHistos(efficiency, "MatchedLpt", trgableHPtOnly);
1558   
1559   // tracker matched with high pt trigger
1560   efficiency->GetNum()->SetRangeUser(kVarTrigger, kHighPtTrig, kHighPtTrig);
1561   FillEffHistos(efficiency, "MatchedHpt", trgableHPtOnly);
1562   
1563   // ------ reset ranges before saving CF containers ------
1564   efficiency->GetNum()->GetAxis(kVarHasTracker)->SetRange();
1565   efficiency->GetDen()->GetAxis(kVarHasTracker)->SetRange();
1566   efficiency->GetNum()->GetAxis(kVarTrigger)->SetRange();
1567   efficiency->GetDen()->GetAxis(kVarTrigger)->SetRange();
1568   efficiency->GetNum()->GetAxis(kVarMatchMC)->SetRange();
1569   efficiency->GetNum()->GetAxis(kVarMCTrigger)->SetRange();
1570   efficiency->GetNum()->GetAxis(kVarDupliTrg)->SetRange();
1571   
1572   // ------ Plot summary ------
1573   
1574   // plot histogram summarizing global efficiencies
1575   TCanvas* cEffSummary = new TCanvas("cEffSummary","Efficiency summary",20,20,310,310);
1576   cEffSummary->SetFillColor(10); cEffSummary->SetHighLightColor(10);
1577   cEffSummary->SetLeftMargin(0.15); cEffSummary->SetBottomMargin(0.15);
1578   effSummary->DrawCopy("etext");
1579   
1580   // --- plot trigger resolution ---
1581   TCanvas* cTriggerResolution = new TCanvas("cTriggerResolution","Trigger resolution",10,10,310,310);
1582   cTriggerResolution->SetFillColor(10); cTriggerResolution->SetHighLightColor(10);
1583   cTriggerResolution->SetLeftMargin(0.15); cTriggerResolution->SetBottomMargin(0.15);
1584   cTriggerResolution->Divide(2,2);
1585   cTriggerResolution->cd(1);
1586   static_cast<TH1*>(fTriggerList->UncheckedAt(kResTrigX11))->DrawCopy();
1587   cTriggerResolution->cd(2);
1588   static_cast<TH1*>(fTriggerList->UncheckedAt(kResTrigY11))->DrawCopy();
1589   cTriggerResolution->cd(3);
1590   static_cast<TH1*>(fTriggerList->UncheckedAt(kResTrigSlopeY))->DrawCopy();
1591   
1592   // --- compute momentum resolution at vertex ---
1593   fPAtVtxList = new TObjArray(100);
1594   fPAtVtxList->SetOwner();
1595   
1596   // define graphs
1597   TGraphAsymmErrors* gMeanResPAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1598   gMeanResPAtVtxVsP->SetName("gMeanResPAtVtxVsP");
1599   gMeanResPAtVtxVsP->SetTitle("<#Delta_{p}> at vertex versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
1600   fPAtVtxList->AddAt(gMeanResPAtVtxVsP, kMeanResPAtVtxVsP);
1601   TGraphAsymmErrors* gMostProbResPAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1602   gMostProbResPAtVtxVsP->SetName("gMostProbResPAtVtxVsP");
1603   gMostProbResPAtVtxVsP->SetTitle("Most probable #Delta_{p} at vertex versus p;p (GeV/c);Most prob. #Delta_{p} (GeV/c)");
1604   fPAtVtxList->AddAt(gMostProbResPAtVtxVsP, kMostProbResPAtVtxVsP);
1605   TGraphAsymmErrors* gSigmaResPAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1606   gSigmaResPAtVtxVsP->SetName("gSigmaResPAtVtxVsP");
1607   gSigmaResPAtVtxVsP->SetTitle("#sigma_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)");
1608   fPAtVtxList->AddAt(gSigmaResPAtVtxVsP, kSigmaResPAtVtxVsP);
1609   
1610   // fit histo and fill graphs
1611   TH2* h = static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsP));
1612   FitLandauGausResVsP(h, "momentum residuals at vertex", gMeanResPAtVtxVsP, gMostProbResPAtVtxVsP, gSigmaResPAtVtxVsP);
1613   
1614   // convert resolution into relative resolution
1615   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
1616   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1617     Double_t x,y;
1618     gSigmaResPAtVtxVsP->GetPoint(i/rebinFactorX-1, x, y);
1619     gSigmaResPAtVtxVsP->SetPoint(i/rebinFactorX-1, x, 100.*y/x);
1620     gSigmaResPAtVtxVsP->SetPointEYlow(i/rebinFactorX-1, 100.*gSigmaResPAtVtxVsP->GetErrorYlow(i/rebinFactorX-1)/x);
1621     gSigmaResPAtVtxVsP->SetPointEYhigh(i/rebinFactorX-1, 100.*gSigmaResPAtVtxVsP->GetErrorYhigh(i/rebinFactorX-1)/x);
1622   }
1623   
1624   // --- compute momentum resolution at first cluster ---
1625   fPAt1stClList = new TObjArray(100);
1626   fPAt1stClList->SetOwner();
1627   
1628   // define graphs
1629   TGraphAsymmErrors* gMeanResPAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1630   gMeanResPAt1stClVsP->SetName("gMeanResPAt1stClVsP");
1631   gMeanResPAt1stClVsP->SetTitle("<#Delta_{p}> at first cluster versus p;p (GeV/c);<#Delta_{p}> (GeV/c)");
1632   fPAt1stClList->AddAt(gMeanResPAt1stClVsP, kMeanResPAt1stClVsP);
1633   TGraphAsymmErrors* gSigmaResPAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1634   gSigmaResPAt1stClVsP->SetName("gSigmaResPAt1stClVsP");
1635   gSigmaResPAt1stClVsP->SetTitle("#sigma_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)");
1636   fPAt1stClList->AddAt(gSigmaResPAt1stClVsP, kSigmaResPAt1stClVsP);
1637   
1638   // fit histo and fill graphs
1639   h = static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAt1stClVsP));
1640   FitGausResVsMom(h, 0., 1., "momentum residuals at first cluster", gMeanResPAt1stClVsP, gSigmaResPAt1stClVsP);
1641   
1642   // convert resolution into relative resolution
1643   rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
1644   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
1645     Double_t x,y;
1646     gSigmaResPAt1stClVsP->GetPoint(i/rebinFactorX-1, x, y);
1647     gSigmaResPAt1stClVsP->SetPoint(i/rebinFactorX-1, x, 100.*y/x);
1648     gSigmaResPAt1stClVsP->SetPointEYlow(i/rebinFactorX-1, 100.*gSigmaResPAt1stClVsP->GetErrorYlow(i/rebinFactorX-1)/x);
1649     gSigmaResPAt1stClVsP->SetPointEYhigh(i/rebinFactorX-1, 100.*gSigmaResPAt1stClVsP->GetErrorYhigh(i/rebinFactorX-1)/x);
1650   }
1651   
1652   // --- compute slope resolution at vertex ---
1653   fSlopeAtVtxList = new TObjArray(100);
1654   fSlopeAtVtxList->SetOwner();
1655   
1656   // define graphs
1657   TGraphAsymmErrors* gMeanResSlopeXAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1658   gMeanResSlopeXAtVtxVsP->SetName("gMeanResSlopeXAtVtxVsP");
1659   gMeanResSlopeXAtVtxVsP->SetTitle("<#Delta_{slope_{X}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{X}}>");
1660   fSlopeAtVtxList->AddAt(gMeanResSlopeXAtVtxVsP, kMeanResSlopeXAtVtxVsP);
1661   TGraphAsymmErrors* gMeanResSlopeYAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1662   gMeanResSlopeYAtVtxVsP->SetName("gMeanResSlopeYAtVtxVsP");
1663   gMeanResSlopeYAtVtxVsP->SetTitle("<#Delta_{slope_{Y}}> at vertex versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
1664   fSlopeAtVtxList->AddAt(gMeanResSlopeYAtVtxVsP, kMeanResSlopeYAtVtxVsP);
1665   TGraphAsymmErrors* gSigmaResSlopeXAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1666   gSigmaResSlopeXAtVtxVsP->SetName("gSigmaResSlopeXAtVtxVsP");
1667   gSigmaResSlopeXAtVtxVsP->SetTitle("#sigma_{slope_{X}} at vertex versus p;p (GeV/c);#sigma_{slope_{X}}");
1668   fSlopeAtVtxList->AddAt(gSigmaResSlopeXAtVtxVsP, kSigmaResSlopeXAtVtxVsP);
1669   TGraphAsymmErrors* gSigmaResSlopeYAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1670   gSigmaResSlopeYAtVtxVsP->SetName("gSigmaResSlopeYAtVtxVsP");
1671   gSigmaResSlopeYAtVtxVsP->SetTitle("#sigma_{slope_{Y}} at vertex versus p;p (GeV/c);#sigma_{slope_{Y}}");
1672   fSlopeAtVtxList->AddAt(gSigmaResSlopeYAtVtxVsP, kSigmaResSlopeYAtVtxVsP);
1673   
1674   // fit histo and fill graphs
1675   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsP)), 0., 2.e-3,
1676                   "slopeX residuals at vertex", gMeanResSlopeXAtVtxVsP, gSigmaResSlopeXAtVtxVsP);
1677   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsP)), 0., 2.e-3,
1678                   "slopeY residuals at vertex", gMeanResSlopeYAtVtxVsP, gSigmaResSlopeYAtVtxVsP);
1679   
1680   // --- compute slope resolution at first cluster ---
1681   fSlopeAt1stClList = new TObjArray(100);
1682   fSlopeAt1stClList->SetOwner();
1683   
1684   // define graphs
1685   TGraphAsymmErrors* gMeanResSlopeXAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1686   gMeanResSlopeXAt1stClVsP->SetName("gMeanResSlopeXAt1stClVsP");
1687   gMeanResSlopeXAt1stClVsP->SetTitle("<#Delta_{slope_{X}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{X}}>");
1688   fSlopeAt1stClList->AddAt(gMeanResSlopeXAt1stClVsP, kMeanResSlopeXAt1stClVsP);
1689   TGraphAsymmErrors* gMeanResSlopeYAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1690   gMeanResSlopeYAt1stClVsP->SetName("gMeanResSlopeYAt1stClVsP");
1691   gMeanResSlopeYAt1stClVsP->SetTitle("<#Delta_{slope_{Y}}> at first cluster versus p;p (GeV/c);<#Delta_{slope_{Y}}>");
1692   fSlopeAt1stClList->AddAt(gMeanResSlopeYAt1stClVsP, kMeanResSlopeYAt1stClVsP);
1693   TGraphAsymmErrors* gSigmaResSlopeXAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1694   gSigmaResSlopeXAt1stClVsP->SetName("gSigmaResSlopeXAt1stClVsP");
1695   gSigmaResSlopeXAt1stClVsP->SetTitle("#sigma_{slope_{X}} at first cluster versus p;p (GeV/c);#sigma_{slope_{X}}");
1696   fSlopeAt1stClList->AddAt(gSigmaResSlopeXAt1stClVsP, kSigmaResSlopeXAt1stClVsP);
1697   TGraphAsymmErrors* gSigmaResSlopeYAt1stClVsP = new TGraphAsymmErrors(fNPBins);
1698   gSigmaResSlopeYAt1stClVsP->SetName("gSigmaResSlopeYAt1stClVsP");
1699   gSigmaResSlopeYAt1stClVsP->SetTitle("#sigma_{slope_{Y}} at first cluster versus p;p (GeV/c);#sigma_{slope_{Y}}");
1700   fSlopeAt1stClList->AddAt(gSigmaResSlopeYAt1stClVsP, kSigmaResSlopeYAt1stClVsP);
1701   
1702   // fit histo and fill graphs
1703   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAt1stClVsP)), 0., 3.e-4,
1704                   "slopeX residuals at first cluster", gMeanResSlopeXAt1stClVsP, gSigmaResSlopeXAt1stClVsP);
1705   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAt1stClVsP)), 0., 3.e-4,
1706                   "slopeY residuals at first cluster", gMeanResSlopeYAt1stClVsP, gSigmaResSlopeYAt1stClVsP);
1707   
1708   // --- compute eta resolution at vertex ---
1709   fEtaAtVtxList = new TObjArray(100);
1710   fEtaAtVtxList->SetOwner();
1711   
1712   // define graphs
1713   TGraphAsymmErrors* gMeanResEtaAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1714   gMeanResEtaAtVtxVsP->SetName("gMeanResEtaAtVtxVsP");
1715   gMeanResEtaAtVtxVsP->SetTitle("<#Delta_{eta}> at vertex versus p;p (GeV/c);<#Delta_{eta}>");
1716   fEtaAtVtxList->AddAt(gMeanResEtaAtVtxVsP, kMeanResEtaAtVtxVsP);
1717   TGraphAsymmErrors* gSigmaResEtaAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1718   gSigmaResEtaAtVtxVsP->SetName("gSigmaResEtaAtVtxVsP");
1719   gSigmaResEtaAtVtxVsP->SetTitle("#sigma_{eta} at vertex versus p;p (GeV/c);#sigma_{eta}");
1720   fEtaAtVtxList->AddAt(gSigmaResEtaAtVtxVsP, kSigmaResEtaAtVtxVsP);
1721   
1722   // fit histo and fill graphs
1723   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsP)), 0., 0.1,
1724                   "eta residuals at vertex", gMeanResEtaAtVtxVsP, gSigmaResEtaAtVtxVsP);
1725   
1726   // --- compute phi resolution at vertex ---
1727   fPhiAtVtxList = new TObjArray(100);
1728   fPhiAtVtxList->SetOwner();
1729   
1730   // define graphs
1731   TGraphAsymmErrors* gMeanResPhiAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1732   gMeanResPhiAtVtxVsP->SetName("gMeanResPhiAtVtxVsP");
1733   gMeanResPhiAtVtxVsP->SetTitle("<#Delta_{phi}> at vertex versus p;p (GeV/c);<#Delta_{phi}>");
1734   fPhiAtVtxList->AddAt(gMeanResPhiAtVtxVsP, kMeanResPhiAtVtxVsP);
1735   TGraphAsymmErrors* gSigmaResPhiAtVtxVsP = new TGraphAsymmErrors(fNPBins);
1736   gSigmaResPhiAtVtxVsP->SetName("gSigmaResPhiAtVtxVsP");
1737   gSigmaResPhiAtVtxVsP->SetTitle("#sigma_{phi} at vertex versus p;p (GeV/c);#sigma_{phi}");
1738   fPhiAtVtxList->AddAt(gSigmaResPhiAtVtxVsP, kSigmaResPhiAtVtxVsP);
1739   
1740   // fit histo and fill graphs
1741   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsP)), 0., 0.01,
1742                   "phi residuals at vertex", gMeanResPhiAtVtxVsP, gSigmaResPhiAtVtxVsP);
1743   
1744   // --- compute DCA resolution and MCS dispersion ---
1745   fDCAList = new TObjArray(100);
1746   fDCAList->SetOwner();
1747   
1748   // define graphs
1749   TGraphAsymmErrors* gMeanPDCAVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1750   gMeanPDCAVsPIn23deg->SetName("gMeanPDCAVsPIn23deg");
1751   gMeanPDCAVsPIn23deg->SetTitle("<p #times DCA> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
1752   fDCAList->AddAt(gMeanPDCAVsPIn23deg, kMeanPDCAVsPIn23deg);
1753   TGraphAsymmErrors* gSigmaPDCAVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1754   gSigmaPDCAVsPIn23deg->SetName("gSigmaPDCAVsPIn23deg");
1755   gSigmaPDCAVsPIn23deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
1756   fDCAList->AddAt(gSigmaPDCAVsPIn23deg, kSigmaPDCAVsPIn23deg);
1757   TGraphAsymmErrors* gMeanPDCAVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1758   gMeanPDCAVsPIn310deg->SetName("gMeanPDCAVsPIn310deg");
1759   gMeanPDCAVsPIn310deg->SetTitle("<p #times DCA> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times DCA> (GeV #times cm)");
1760   fDCAList->AddAt(gMeanPDCAVsPIn310deg, kMeanPDCAVsPIn310deg);
1761   TGraphAsymmErrors* gSigmaPDCAVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1762   gSigmaPDCAVsPIn310deg->SetName("gSigmaPDCAVsPIn310deg");
1763   gSigmaPDCAVsPIn310deg->SetTitle("#sigma_{p #times DCA} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times DCA} (GeV #times cm)");
1764   fDCAList->AddAt(gSigmaPDCAVsPIn310deg, kSigmaPDCAVsPIn310deg);
1765   
1766   TGraphAsymmErrors* gMeanPMCSAngVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1767   gMeanPMCSAngVsPIn23deg->SetName("gMeanPMCSAngVsPIn23deg");
1768   gMeanPMCSAngVsPIn23deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
1769   fDCAList->AddAt(gMeanPMCSAngVsPIn23deg, kMeanPMCSAngVsPIn23deg);
1770   TGraphAsymmErrors* gSigmaPMCSAngVsPIn23deg = new TGraphAsymmErrors(fNPBins);
1771   gSigmaPMCSAngVsPIn23deg->SetName("gSigmaPMCSAngVsPIn23deg");
1772   gSigmaPMCSAngVsPIn23deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [2,3[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
1773   fDCAList->AddAt(gSigmaPMCSAngVsPIn23deg, kSigmaPMCSAngVsPIn23deg);
1774   TGraphAsymmErrors* gMeanPMCSAngVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1775   gMeanPMCSAngVsPIn310deg->SetName("gMeanPMCSAngVsPIn310deg");
1776   gMeanPMCSAngVsPIn310deg->SetTitle("<p #times #Delta#theta_{MCS}> versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);<p #times #Delta#theta_{MCS}> (GeV)");
1777   fDCAList->AddAt(gMeanPMCSAngVsPIn310deg, kMeanPMCSAngVsPIn310deg);
1778   TGraphAsymmErrors* gSigmaPMCSAngVsPIn310deg = new TGraphAsymmErrors(fNPBins);
1779   gSigmaPMCSAngVsPIn310deg->SetName("gSigmaPMCSAngVsPIn310deg");
1780   gSigmaPMCSAngVsPIn310deg->SetTitle("#sigma_{p #times #Delta#theta_{MCS}} versus p for tracks within [3,10[ degrees at absorber end;p (GeV/c);#sigma_{p #times #Delta#theta_{MCS}} (GeV)");
1781   fDCAList->AddAt(gSigmaPMCSAngVsPIn310deg, kSigmaPMCSAngVsPIn310deg);
1782   
1783   // fit histo and fill graphs
1784   FitPDCAVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPIn23deg)),
1785                "p*DCA (tracks in [2,3] deg.)", gMeanPDCAVsPIn23deg, gSigmaPDCAVsPIn23deg);
1786   FitPDCAVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPIn310deg)),
1787                "p*DCA (tracks in [3,10] deg.)", gMeanPDCAVsPIn310deg, gSigmaPDCAVsPIn310deg);
1788   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPMCSAngVsPIn23deg)), 0., 2.e-3,
1789                   "p*MCSAngle (tracks in [2,3] deg.)", gMeanPMCSAngVsPIn23deg, gSigmaPMCSAngVsPIn23deg);
1790   FitGausResVsMom(static_cast<TH2*>(fTrackerList->UncheckedAt(kPMCSAngVsPIn310deg)), 0., 2.e-3,
1791                   "p*MCSAngle (tracks in [3,10] deg.)", gMeanPMCSAngVsPIn310deg, gSigmaPMCSAngVsPIn310deg);
1792   
1793   // --- compute cluster resolution ---
1794   fClusterList = new TObjArray(100);
1795   fClusterList->SetOwner();
1796   
1797   // define graphs per chamber
1798   TGraphErrors* gMeanResClXVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1799   gMeanResClXVsCh->SetName("gMeanResClXVsCh");
1800   gMeanResClXVsCh->SetTitle("cluster-trackRef residual-X per Ch: mean;chamber ID;<#Delta_{X}> (cm)");
1801   gMeanResClXVsCh->SetMarkerStyle(kFullDotLarge);
1802   fClusterList->AddAt(gMeanResClXVsCh, kMeanResClXVsCh);
1803   TGraphErrors* gMeanResClYVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1804   gMeanResClYVsCh->SetName("gMeanResClYVsCh");
1805   gMeanResClYVsCh->SetTitle("cluster-trackRef residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)");
1806   gMeanResClYVsCh->SetMarkerStyle(kFullDotLarge);
1807   fClusterList->AddAt(gMeanResClYVsCh, kMeanResClYVsCh);
1808   TGraphErrors* gSigmaResClXVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1809   gSigmaResClXVsCh->SetName("gSigmaResClXVsCh");
1810   gSigmaResClXVsCh->SetTitle("cluster-trackRef residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)");
1811   gSigmaResClXVsCh->SetMarkerStyle(kFullDotLarge);
1812   fClusterList->AddAt(gSigmaResClXVsCh, kSigmaResClXVsCh);
1813   TGraphErrors* gSigmaResClYVsCh = new TGraphErrors(AliMUONConstants::NTrackingCh());
1814   gSigmaResClYVsCh->SetName("gSigmaResClYVsCh");
1815   gSigmaResClYVsCh->SetTitle("cluster-trackRef residual-Y per Ch: sigma;chamber ID;#sigma_{Y} (cm)");
1816   gSigmaResClYVsCh->SetMarkerStyle(kFullDotLarge);
1817   fClusterList->AddAt(gSigmaResClYVsCh, kSigmaResClYVsCh);
1818   
1819   // define graphs per DE
1820   TGraphErrors* gMeanResClXVsDE = new TGraphErrors(fNDE);
1821   gMeanResClXVsDE->SetName("gMeanResClXVsDE");
1822   gMeanResClXVsDE->SetTitle("cluster-trackRef residual-X per DE: mean;DE ID;<#Delta_{X}> (cm)");
1823   gMeanResClXVsDE->SetMarkerStyle(kFullDotLarge);
1824   fClusterList->AddAt(gMeanResClXVsDE, kMeanResClXVsDE);
1825   TGraphErrors* gMeanResClYVsDE = new TGraphErrors(fNDE);
1826   gMeanResClYVsDE->SetName("gMeanResClYVsDE");
1827   gMeanResClYVsDE->SetTitle("cluster-trackRef residual-Y per dE: mean;DE ID;<#Delta_{Y}> (cm)");
1828   gMeanResClYVsDE->SetMarkerStyle(kFullDotLarge);
1829   fClusterList->AddAt(gMeanResClYVsDE, kMeanResClYVsDE);
1830   TGraphErrors* gSigmaResClXVsDE = new TGraphErrors(fNDE);
1831   gSigmaResClXVsDE->SetName("gSigmaResClXVsDE");
1832   gSigmaResClXVsDE->SetTitle("cluster-trackRef residual-X per DE: sigma;DE ID;#sigma_{X} (cm)");
1833   gSigmaResClXVsDE->SetMarkerStyle(kFullDotLarge);
1834   fClusterList->AddAt(gSigmaResClXVsDE, kSigmaResClXVsDE);
1835   TGraphErrors* gSigmaResClYVsDE = new TGraphErrors(fNDE);
1836   gSigmaResClYVsDE->SetName("gSigmaResClYVsDE");
1837   gSigmaResClYVsDE->SetTitle("cluster-trackRef residual-Y per DE: sigma;DE ID;#sigma_{Y} (cm)");
1838   gSigmaResClYVsDE->SetMarkerStyle(kFullDotLarge);
1839   fClusterList->AddAt(gSigmaResClYVsDE, kSigmaResClYVsDE);
1840   
1841   // fit histo and fill graphs per chamber
1842   Double_t clusterResPerCh[10][2];
1843   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
1844     TH1D *tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClXVsCh))->ProjectionY("tmp",i+1,i+1,"e");
1845     FitClusterResidual(tmp, i, clusterResPerCh[i][0], gMeanResClXVsCh, gSigmaResClXVsCh);
1846     delete tmp;
1847     tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClYVsCh))->ProjectionY("tmp",i+1,i+1,"e");
1848     FitClusterResidual(tmp, i, clusterResPerCh[i][1], gMeanResClYVsCh, gSigmaResClYVsCh);
1849     delete tmp;
1850   }
1851   
1852   // fit histo and fill graphs per DE
1853   Double_t clusterResPerDE[200][2];
1854   for (Int_t i = 0; i < fNDE; i++) {
1855     TH1D *tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClXVsDE))->ProjectionY("tmp",i+1,i+1,"e");
1856     FitClusterResidual(tmp, i, clusterResPerDE[i][0], gMeanResClXVsDE, gSigmaResClXVsDE);
1857     delete tmp;
1858     tmp = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClYVsDE))->ProjectionY("tmp",i+1,i+1,"e");
1859     FitClusterResidual(tmp, i, clusterResPerDE[i][1], gMeanResClYVsDE, gSigmaResClYVsDE);
1860     delete tmp;
1861   }
1862   
1863   // set DE graph labels
1864   TAxis* xAxis = static_cast<TH2*>(fTrackerList->UncheckedAt(kResClXVsDE))->GetXaxis();
1865   gMeanResClXVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1866   gMeanResClYVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1867   gSigmaResClXVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1868   gSigmaResClYVsDE->GetXaxis()->Set(fNDE, 0.5, fNDE+0.5);
1869   for (Int_t i = 1; i <= fNDE; i++) {
1870     const char* label = xAxis->GetBinLabel(i);
1871     gMeanResClXVsDE->GetXaxis()->SetBinLabel(i, label);
1872     gMeanResClYVsDE->GetXaxis()->SetBinLabel(i, label);
1873     gSigmaResClXVsDE->GetXaxis()->SetBinLabel(i, label);
1874     gSigmaResClYVsDE->GetXaxis()->SetBinLabel(i, label);
1875   }
1876   
1877   // --- diplay momentum residuals at vertex ---
1878   TCanvas* cResPAtVtx = DrawVsAng("cResPAtVtx", "momentum residual at vertex in 3 angular regions",
1879                                   static_cast<TH1*>(fTrackerList->UncheckedAt(kResPAtVtx)),
1880                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsAngleAtAbsEnd)));
1881   fPAtVtxList->AddAt(cResPAtVtx, kcResPAtVtx);
1882   TCanvas* cResPAtVtxMC = DrawVsAng("cResPAtVtxMC", "momentum residual at vertex in 3 MC angular regions",
1883                                     static_cast<TH1*>(fTrackerList->UncheckedAt(kResPAtVtx)),
1884                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsMCAngle)));
1885   fPAtVtxList->AddAt(cResPAtVtxMC, kcResPAtVtxMC);
1886   TCanvas* cResPAtVtxVsPosAbsEndMC = DrawVsPos("cResPAtVtxVsPosAbsEndMC", "momentum residual at vertex versus position at absorber end in 3 MC angular regions",
1887                                                static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPosAbsEndIn02degMC)),
1888                                                static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPosAbsEndIn23degMC)),
1889                                                static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPosAbsEndIn310degMC)));
1890   fPAtVtxList->AddAt(cResPAtVtxVsPosAbsEndMC, kcResPAtVtxVsPosAbsEndMC);
1891   TCanvas* cResPAtVtxVsPIn23deg = DrawFitLandauGausResPVsP("cResPAtVtxVsPIn23deg", "momentum residual for tracks between 2 and 3 degrees",
1892                                                       static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPIn23deg)),
1893                                                       10, "momentum residuals at vertex (tracks in [2,3] deg.)");
1894   fPAtVtxList->AddAt(cResPAtVtxVsPIn23deg, kcResPAtVtxVsPIn23deg);
1895   TCanvas* cResPAtVtxVsPIn310deg = DrawFitLandauGausResPVsP("cResPAtVtxVsPIn310deg", "momentum residual for tracks between 3 and 10 degrees",
1896                                                        static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPIn310deg)),
1897                                                        10, "momentum residuals at vertex (tracks in [3,10] deg.)");
1898   fPAtVtxList->AddAt(cResPAtVtxVsPIn310deg, kcResPAtVtxVsPIn310deg);
1899   TCanvas* cResPAtVtxVsPIn02degMC = DrawResPVsP("cResPAtVtxVsPIn02degMC", "momentum residuals for tracks with MC angle < 2 degrees",
1900                                            static_cast<TH2*>(fTrackerList->UncheckedAt(kResPAtVtxVsPIn02degMC)), 5);
1901   fPAtVtxList->AddAt(cResPAtVtxVsPIn02degMC, kcResPAtVtxVsPIn02degMC);
1902   
1903   // --- diplay slopeX residuals at vertex ---
1904   TCanvas* cResSlopeXAtVtx = DrawVsAng("cResSlopeXAtVtx", "slope_{X} residual at vertex in 3 angular regions",
1905                                        static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeXAtVtx)),
1906                                        static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsAngleAtAbsEnd)));
1907   fSlopeAtVtxList->AddAt(cResSlopeXAtVtx, kcResSlopeXAtVtx);
1908   TCanvas* cResSlopeYAtVtx = DrawVsAng("cResSlopeYAtVtx", "slope_{Y} residual at vertex in 3 angular regions",
1909                                        static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeYAtVtx)),
1910                                        static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsAngleAtAbsEnd)));
1911   fSlopeAtVtxList->AddAt(cResSlopeYAtVtx, kcResSlopeYAtVtx);
1912   TCanvas* cResSlopeXAtVtxMC = DrawVsAng("cResSlopeXAtVtxMC", "slope_{X} residual at vertex in 3 MC angular regions",
1913                                          static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeXAtVtx)),
1914                                          static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsMCAngle)));
1915   fSlopeAtVtxList->AddAt(cResSlopeXAtVtxMC, kcResSlopeXAtVtxMC);
1916   TCanvas* cResSlopeYAtVtxMC = DrawVsAng("cResSlopeYAtVtxMC", "slope_{Y} residual at vertex in 3 MC angular regions",
1917                                          static_cast<TH1*>(fTrackerList->UncheckedAt(kResSlopeYAtVtx)),
1918                                          static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsMCAngle)));
1919   fSlopeAtVtxList->AddAt(cResSlopeYAtVtxMC, kcResSlopeYAtVtxMC);
1920   TCanvas* cResSlopeXAtVtxVsPosAbsEndMC = DrawVsPos("cResSlopeXAtVtxVsPosAbsEndMC", "slope_{X} residual at vertex versus position at absorber end in 3 MC angular regions",
1921                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsPosAbsEndIn02degMC)),
1922                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsPosAbsEndIn23degMC)),
1923                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeXAtVtxVsPosAbsEndIn310degMC)));
1924   fSlopeAtVtxList->AddAt(cResSlopeXAtVtxVsPosAbsEndMC, kcResSlopeXAtVtxVsPosAbsEndMC);
1925   TCanvas* cResSlopeYAtVtxVsPosAbsEndMC = DrawVsPos("cResSlopeYAtVtxVsPosAbsEndMC", "slope_{Y} residual at vertex versus position at absorber end in 3 MC angular regions",
1926                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsPosAbsEndIn02degMC)),
1927                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsPosAbsEndIn23degMC)),
1928                                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResSlopeYAtVtxVsPosAbsEndIn310degMC)));
1929   fSlopeAtVtxList->AddAt(cResSlopeYAtVtxVsPosAbsEndMC, kcResSlopeYAtVtxVsPosAbsEndMC);
1930   
1931   // --- diplay eta residuals at vertex ---
1932   TCanvas* cResEtaAtVtx = DrawVsAng("cResEtaAtVtx", "eta residual at vertex in 3 angular regions",
1933                                     static_cast<TH1*>(fTrackerList->UncheckedAt(kResEtaAtVtx)),
1934                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsAngleAtAbsEnd)));
1935   fEtaAtVtxList->AddAt(cResEtaAtVtx, kcResEtaAtVtx);
1936   TCanvas* cResEtaAtVtxMC = DrawVsAng("cResEtaAtVtxMC", "eta residual at vertex in 3 MC angular regions",
1937                                       static_cast<TH1*>(fTrackerList->UncheckedAt(kResEtaAtVtx)),
1938                                       static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsMCAngle)));
1939   fEtaAtVtxList->AddAt(cResEtaAtVtxMC, kcResEtaAtVtxMC);
1940   TCanvas* cResEtaAtVtxVsPosAbsEndMC = DrawVsPos("cResEtaAtVtxVsPosAbsEndMC", "eta residual at vertex versus position at absorber end in 3 MC angular regions",
1941                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsPosAbsEndIn02degMC)),
1942                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsPosAbsEndIn23degMC)),
1943                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResEtaAtVtxVsPosAbsEndIn310degMC)));
1944   fEtaAtVtxList->AddAt(cResEtaAtVtxVsPosAbsEndMC, kcResEtaAtVtxVsPosAbsEndMC);
1945   
1946   // --- diplay phi residuals at vertex ---
1947   TCanvas* cResPhiAtVtx = DrawVsAng("cResPhiAtVtx", "phi residual at vertex in 3 angular regions",
1948                                     static_cast<TH1*>(fTrackerList->UncheckedAt(kResPhiAtVtx)),
1949                                     static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsAngleAtAbsEnd)));
1950   fPhiAtVtxList->AddAt(cResPhiAtVtx, kcResPhiAtVtx);
1951   TCanvas* cResPhiAtVtxMC = DrawVsAng("cResPhiAtVtxMC", "phi residual at vertex in 3 MC angular regions",
1952                                       static_cast<TH1*>(fTrackerList->UncheckedAt(kResPhiAtVtx)),
1953                                       static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsMCAngle)));
1954   fPhiAtVtxList->AddAt(cResPhiAtVtxMC, kcResPhiAtVtxMC);
1955   TCanvas* cResPhiAtVtxVsPosAbsEndMC = DrawVsPos("cResPhiAtVtxVsPosAbsEndMC", "phi residual at vertex versus position at absorber end in 3 MC angular regions",
1956                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsPosAbsEndIn02degMC)),
1957                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsPosAbsEndIn23degMC)),
1958                                                  static_cast<TH2*>(fTrackerList->UncheckedAt(kResPhiAtVtxVsPosAbsEndIn310degMC)));
1959   fPhiAtVtxList->AddAt(cResPhiAtVtxVsPosAbsEndMC, kcResPhiAtVtxVsPosAbsEndMC);
1960   
1961   // --- diplay P*DCA ---
1962   TCanvas* cPDCA = DrawVsAng("cPDCA", "p #times DCA in 3 angular regions",
1963                              static_cast<TH1*>(fTrackerList->UncheckedAt(kPDCA)),
1964                              static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsAngleAtAbsEnd)));
1965   fDCAList->AddAt(cPDCA, kcPDCA);
1966   TCanvas* cPDCAMC = DrawVsAng("cPDCAMC", "p #times DCA in 3 MC angular regions",
1967                                static_cast<TH1*>(fTrackerList->UncheckedAt(kPDCA)),
1968                                static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsMCAngle)));
1969   fDCAList->AddAt(cPDCAMC, kcPDCAMC);
1970   TCanvas* cPDCAVsPosAbsEndMC = DrawVsPos("cPDCAVsPosAbsEndMC", "p #times DCA versus position at absorber end in 3 MC angular regions",
1971                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPosAbsEndIn02degMC)),
1972                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPosAbsEndIn23degMC)),
1973                                   static_cast<TH2*>(fTrackerList->UncheckedAt(kPDCAVsPosAbsEndIn310degMC)));
1974   fDCAList->AddAt(cPDCAVsPosAbsEndMC, kcPDCAVsPosAbsEndMC);
1975   
1976   // post param containers
1977   PostData(4, fEfficiencyList);
1978   PostData(5, fPAtVtxList);
1979   PostData(6, fSlopeAtVtxList);
1980   PostData(7, fEtaAtVtxList);
1981   PostData(8, fPhiAtVtxList);
1982   PostData(9, fPAt1stClList);
1983   PostData(10, fSlopeAt1stClList);
1984   PostData(11, fDCAList);
1985   PostData(12, fClusterList);
1986 }
1987
1988
1989 //________________________________________________________________________
1990 Bool_t AliAnalysisTaskMuonPerformance::GetEfficiency(AliCFEffGrid* efficiency, Double_t& calcEff, Double_t& calcEffErr)
1991 {
1992   //
1993   /// Calculate the efficiency when cuts on the THnSparse are applied
1994   //
1995
1996   Bool_t isGood = kTRUE;
1997
1998   TH1* histo = 0x0;
1999   Double_t sum[2] = {0., 0.};
2000   for ( Int_t ihisto=0; ihisto<2; ihisto++ ) {
2001     histo = ( ihisto==0 ) ? efficiency->GetNum()->Project(kVarCharge) : efficiency->GetDen()->Project(kVarCharge);
2002     sum[ihisto] = histo->Integral();
2003     delete histo;
2004   }
2005
2006   if ( sum[1] == 0. ) isGood = kFALSE;
2007   
2008   calcEff = ( isGood ) ? sum[0]/sum[1] : 0.;
2009   if ( calcEff > 1. ) isGood = kFALSE;
2010
2011   calcEffErr = ( isGood ) ? TMath::Sqrt(calcEff*(1-calcEff)/sum[1]) : 0.;
2012
2013   return isGood;
2014 }
2015
2016 //________________________________________________________________________
2017 Int_t AliAnalysisTaskMuonPerformance::RecoTrackMother(AliMCParticle* mcParticle)
2018 {
2019   //
2020   /// Find track mother from kinematics
2021   //
2022
2023   if ( ! mcParticle )
2024     return kUnknownPart;
2025
2026   Int_t recoPdg = mcParticle->PdgCode();
2027
2028   // Track is not a muon
2029   if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
2030
2031   Int_t imother = mcParticle->GetMother();
2032
2033   Bool_t isFirstMotherHF = kFALSE;
2034   Int_t step = 0;
2035
2036   //printf("\n"); // REMEMBER TO CUT
2037   while ( imother >= 0 ) {
2038     TParticle* part = static_cast<AliMCParticle*>(fMCEvent->GetTrack(imother))->Particle();
2039     //printf("Particle %s  -> %s\n", part->GetName(), TMCProcessName[part->GetUniqueID()]); // REMEMBER TO CUT
2040
2041     if ( part->GetUniqueID() == kPHadronic ) 
2042       return kSecondaryMu;
2043
2044     Int_t absPdg = TMath::Abs(part->GetPdgCode());
2045
2046     step++;
2047     if ( step == 1 ) // only for first mother
2048       isFirstMotherHF = ( ( absPdg >= 400  && absPdg < 600 ) || 
2049                           ( absPdg >= 4000 && absPdg < 6000 ) );
2050       
2051     if ( absPdg < 4 )
2052       return kPrimaryMu;
2053
2054     if ( isFirstMotherHF) {
2055       if ( absPdg == 4 )
2056         return kCharmMu;
2057       else if ( absPdg == 5 )
2058         return kBeautyMu;
2059     }
2060
2061     imother = part->GetFirstMother();
2062   }
2063
2064   return kPrimaryMu;
2065
2066 }
2067
2068 //________________________________________________________________________
2069 Float_t AliAnalysisTaskMuonPerformance::GetBinThetaAbsEnd(Float_t RAtAbsEnd, Bool_t isTheta)
2070 {
2071   //
2072   /// Get bin of theta at absorber end region
2073   //
2074   Float_t thetaDeg = ( isTheta ) ? RAtAbsEnd : TMath::ATan( RAtAbsEnd / 505. );
2075   thetaDeg *= TMath::RadToDeg();
2076   if ( thetaDeg < 2. )
2077     return 0.;
2078   else if ( thetaDeg < 3. )
2079     return 1.;
2080   else if ( thetaDeg < 10. )
2081     return 2.;
2082
2083   return 3.;
2084 }
2085
2086 //________________________________________________________________________
2087 void AliAnalysisTaskMuonPerformance::FillContainerInfoReco(Double_t* containerInput, AliESDMuonTrack* esdTrack,
2088                                                            Bool_t isValid, Int_t mcID)
2089 {
2090   //
2091   /// Fill container info (except kVarMatchMC, kVarMCTrigger, kVarCent, kVarDupliTrg) for reconstructed tracks
2092   //
2093   
2094   AliMCParticle* mcPart = (mcID >= 0) ? static_cast<AliMCParticle*>(fMCEvent->GetTrack(mcID)) : 0x0;
2095   
2096   if (fUseMCKinematics && mcPart) {
2097     containerInput[kVarPt] = mcPart->Pt();
2098     containerInput[kVarEta] = mcPart->Eta();
2099     containerInput[kVarPhi] = mcPart->Phi();
2100   } else {
2101     containerInput[kVarPt] = esdTrack->Pt();
2102     containerInput[kVarEta] = esdTrack->Eta();
2103     containerInput[kVarPhi] = esdTrack->Phi();
2104   }
2105   containerInput[kVarThetaZones] = GetBinThetaAbsEnd(esdTrack->GetRAtAbsorberEnd());
2106   containerInput[kVarCharge] = static_cast<Double_t>(esdTrack->Charge());
2107   containerInput[kVarHasTracker] = static_cast<Double_t>(esdTrack->ContainTrackerData() && isValid);
2108   if (esdTrack->GetMatchTrigger() == 0) containerInput[kVarTrigger] = static_cast<Double_t>(kNoMatchTrig);
2109   else if (esdTrack->GetMatchTrigger() == 1) containerInput[kVarTrigger] = static_cast<Double_t>(kAllPtTrig);
2110   else if (esdTrack->GetMatchTrigger() == 2) containerInput[kVarTrigger] = static_cast<Double_t>(kLowPtTrig);
2111   else if (esdTrack->GetMatchTrigger() == 3) containerInput[kVarTrigger] = static_cast<Double_t>(kHighPtTrig);
2112   containerInput[kVarMotherType] = static_cast<Double_t>(RecoTrackMother(mcPart));
2113
2114 }
2115
2116 //________________________________________________________________________
2117 void AliAnalysisTaskMuonPerformance::FillContainerInfoMC(Double_t* containerInput, AliMCParticle* mcPart)
2118 {
2119   //
2120   /// Fill container info (except kVarMatchMC, kVarMCTrigger, kVarHasTracker and kVarTrigger, kVarCent, kVarDupliTrg) for MC tracks
2121   //
2122   
2123   containerInput[kVarPt] = mcPart->Pt();
2124   containerInput[kVarEta] = mcPart->Eta();
2125   containerInput[kVarPhi] = mcPart->Phi();
2126   containerInput[kVarThetaZones] = GetBinThetaAbsEnd(TMath::Pi()-mcPart->Theta(),kTRUE);
2127   containerInput[kVarCharge] = static_cast<Double_t>(mcPart->Charge())/3.;
2128   containerInput[kVarMotherType] = static_cast<Double_t>(RecoTrackMother(mcPart));
2129   
2130 }
2131
2132 //________________________________________________________________________
2133 Double_t langaufun(Double_t *x, Double_t *par) {
2134   
2135   //Fit parameters:
2136   //par[0]=Width (scale) parameter of Landau density
2137   //par[1]=Most Probable (MP, location) parameter of Landau density
2138   //par[2]=Total area (integral -inf to inf, normalization constant)
2139   //par[3]=Width (sigma) of convoluted Gaussian function
2140   //
2141   //In the Landau distribution (represented by the CERNLIB approximation), 
2142   //the maximum is located at x=-0.22278298 with the location parameter=0.
2143   //This shift is corrected within this function, so that the actual
2144   //maximum is identical to the MP parameter.
2145   
2146   // Numeric constants
2147   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
2148   Double_t mpshift  = -0.22278298;       // Landau maximum location
2149   
2150   // Control constants
2151   Double_t np = 100.0; // number of convolution steps
2152   Double_t sc = 5.0;   // convolution extends to +-sc Gaussian sigmas
2153   
2154   // Variables
2155   Double_t xx;
2156   Double_t mpc;
2157   Double_t fland;
2158   Double_t sum = 0.0;
2159   Double_t xlow,xupp;
2160   Double_t step;
2161   Double_t i;
2162   
2163   
2164   // MP shift correction
2165   mpc = par[1] - mpshift * par[0]; 
2166   
2167   // Range of convolution integral
2168   xlow = x[0] - sc * par[3];
2169   xupp = x[0] + sc * par[3];
2170   
2171   step = (xupp-xlow) / np;
2172   
2173   // Convolution integral of Landau and Gaussian by sum
2174   for(i=1.0; i<=np/2; i++) {
2175     xx = xlow + (i-.5) * step;
2176     //change x -> -x because the tail of the Landau is at the left here...
2177     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
2178     sum += fland * TMath::Gaus(x[0],xx,par[3]);
2179     
2180     //change x -> -x because the tail of the Landau is at the left here...
2181     xx = xupp - (i-.5) * step;
2182     fland = TMath::Landau(-xx,mpc,par[0]) / par[0];
2183     sum += fland * TMath::Gaus(x[0],xx,par[3]);
2184   }
2185   
2186   return (par[2] * step * sum * invsq2pi / par[3]);
2187 }
2188
2189 //________________________________________________________________________
2190 void AliAnalysisTaskMuonPerformance::FitLandauGausResVsP(TH2* h, const char* fitting, TGraphAsymmErrors* gMean,
2191                                                          TGraphAsymmErrors* gMostProb, TGraphAsymmErrors* gSigma)
2192 {
2193   /// generic function to fit residuals versus momentum with a landau convoluted with a gaussian
2194   
2195   static TF1* fLandauGaus = 0x0;
2196   if (!fLandauGaus) fLandauGaus = new TF1("fLandauGaus",langaufun,h->GetYaxis()->GetBinLowEdge(1),h->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1),4);
2197   
2198   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
2199   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
2200     
2201     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,fNPBins)<<flush;
2202     
2203     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
2204     
2205     // first fit
2206     fLandauGaus->SetParameters(0.2,0.,(Double_t)tmp->GetEntries(),1.);
2207     tmp->Fit("fLandauGaus","WWNQ");
2208     
2209     // rebin histo
2210     Double_t fwhm = fLandauGaus->GetParameter(0);
2211     Double_t sigma = fLandauGaus->GetParameter(3);
2212     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
2213     Int_t rebin = TMath::Max(static_cast<Int_t>(0.5*sigmaP/tmp->GetBinWidth(1)),1);
2214     while (tmp->GetNbinsX()%rebin!=0) rebin--;
2215     tmp->Rebin(rebin);
2216     
2217     // second fit
2218     tmp->Fit("fLandauGaus","NQ");
2219     
2220     // get fit results and fill histograms
2221     fwhm = fLandauGaus->GetParameter(0);
2222     sigma = fLandauGaus->GetParameter(3);
2223     sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
2224     Double_t fwhmErr = fLandauGaus->GetParError(0);
2225     Double_t sigmaErr = fLandauGaus->GetParError(3);
2226     Double_t sigmaPErr = TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + fwhm*fwhm*fwhmErr*fwhmErr/(64.*log(2.)*log(2.))) / sigmaP;
2227     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
2228     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
2229     h->GetXaxis()->SetRange();
2230     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
2231     gMean->SetPoint(i/rebinFactorX-1, p, tmp->GetMean());
2232     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], tmp->GetMeanError(), tmp->GetMeanError());
2233     gMostProb->SetPoint(i/rebinFactorX-1, p, -fLandauGaus->GetParameter(1));
2234     gMostProb->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fLandauGaus->GetParError(1), fLandauGaus->GetParError(1));
2235     gSigma->SetPoint(i/rebinFactorX-1, p, sigmaP);
2236     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], sigmaPErr, sigmaPErr);
2237     
2238     // clean memory
2239     delete tmp;
2240   }
2241   
2242   cout<<Form("\rFitting %s... %d/%d",fitting,fNPBins,fNPBins)<<endl;
2243   
2244 }
2245
2246 //________________________________________________________________________
2247 void AliAnalysisTaskMuonPerformance::FitGausResVsMom(TH2* h, const Double_t mean0,
2248                                                      const Double_t sigma0, const char* fitting,
2249                                                      TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
2250 {
2251   /// generic function to fit residuals versus momentum with a gaussian
2252   
2253   static TF1* fGaus = 0x0;
2254   if (!fGaus) fGaus = new TF1("fGaus","gaus");
2255   
2256   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
2257   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
2258     
2259     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,fNPBins)<<flush;
2260     
2261     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
2262     
2263     // first fit
2264     fGaus->SetParameters(tmp->GetEntries(), mean0, sigma0);
2265     tmp->Fit("fGaus","WWNQ");
2266     
2267     // rebin histo
2268     Int_t rebin = TMath::Max(static_cast<Int_t>(0.5*fGaus->GetParameter(2)/tmp->GetBinWidth(1)),1);
2269     while (tmp->GetNbinsX()%rebin!=0) rebin--;
2270     tmp->Rebin(rebin);
2271     
2272     // second fit
2273     tmp->Fit("fGaus","NQ");
2274     
2275     // get fit results and fill histograms
2276     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
2277     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
2278     h->GetXaxis()->SetRange();
2279     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
2280     gMean->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(1));
2281     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(1), fGaus->GetParError(1));
2282     gSigma->SetPoint(i/rebinFactorX-1, p, fGaus->GetParameter(2));
2283     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fGaus->GetParError(2), fGaus->GetParError(2));
2284     
2285     // clean memory
2286     delete tmp;
2287   }
2288   
2289   cout<<Form("\rFitting %s... %d/%d",fitting,fNPBins,fNPBins)<<endl;
2290   
2291 }
2292
2293 //________________________________________________________________________
2294 void AliAnalysisTaskMuonPerformance::FitPDCAVsMom(TH2* h, const char* fitting,
2295                                                   TGraphAsymmErrors* gMean, TGraphAsymmErrors* gSigma)
2296 {
2297   /// generic function to fit p*DCA distributions
2298   
2299   static TF1* fPGaus = 0x0;
2300   if (!fPGaus) fPGaus = new TF1("fPGaus","x*gaus");
2301   
2302   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/fNPBins, 1);
2303   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
2304     
2305     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,fNPBins)<<flush;
2306     
2307     TH1D *tmp = h->ProjectionY("tmp",i-rebinFactorX+1,i,"e");
2308     
2309     // rebin histo
2310     Int_t rebin = static_cast<Int_t>(25*(tmp->GetNbinsX()/(tmp->GetBinLowEdge(tmp->GetNbinsX()+1)-tmp->GetBinLowEdge(1))));
2311     while (tmp->GetNbinsX()%rebin!=0) rebin--;
2312     tmp->Rebin(rebin);
2313     
2314     // fit
2315     fPGaus->SetParameters(1.,0.,80.);
2316     tmp->Fit("fPGaus","NQ");
2317     
2318     // get fit results and fill histograms
2319     h->GetXaxis()->SetRange(i-rebinFactorX+1,i);
2320     Double_t p = (tmp->GetEntries() > 0) ? h->GetMean() : 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
2321     h->GetXaxis()->SetRange();
2322     Double_t pErr[2] = {p-h->GetBinLowEdge(i-rebinFactorX+1), h->GetBinLowEdge(i+1)-p};
2323     gMean->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(1));
2324     gMean->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(1), fPGaus->GetParError(1));
2325     gSigma->SetPoint(i/rebinFactorX-1, p, fPGaus->GetParameter(2));
2326     gSigma->SetPointError(i/rebinFactorX-1, pErr[0], pErr[1], fPGaus->GetParError(2), fPGaus->GetParError(2));
2327     
2328     // clean memory
2329     delete tmp;
2330   }
2331   
2332   cout<<Form("\rFitting %s... %d/%d",fitting,fNPBins,fNPBins)<<endl;
2333   
2334 }
2335
2336 //________________________________________________________________________
2337 void AliAnalysisTaskMuonPerformance::FitClusterResidual(TH1* h, Int_t i, Double_t& sigma,
2338                                                         TGraphErrors* gMean, TGraphErrors* gSigma)
2339 {
2340   /// fill graphs with residual mean and sigma
2341   
2342   static TF1* fRGaus = 0x0;
2343   Double_t mean, meanErr, sigmaErr;
2344   
2345   if (fFitResiduals) {
2346     
2347     if (!fRGaus) fRGaus = new TF1("fRGaus","gaus");
2348     
2349     // first fit
2350     Double_t xMin = h->GetXaxis()->GetXmin();
2351     Double_t xMax = h->GetXaxis()->GetXmax();
2352     fRGaus->SetRange(xMin, xMax);
2353     fRGaus->SetParameters(h->GetEntries(), 0., 0.1);
2354     fRGaus->SetParLimits(1, xMin, xMax);
2355     h->Fit("fRGaus", "WWNQ");
2356     
2357     // rebin histo
2358     Int_t rebin = TMath::Max(static_cast<Int_t>(0.3*fRGaus->GetParameter(2)/h->GetBinWidth(1)),1);
2359     while (h->GetNbinsX()%rebin!=0) rebin--;
2360     h->Rebin(rebin);
2361     
2362     // second fit
2363     xMin = TMath::Max(fRGaus->GetParameter(1)-10.*fRGaus->GetParameter(2), h->GetXaxis()->GetXmin());
2364     xMax = TMath::Min(fRGaus->GetParameter(1)+10.*fRGaus->GetParameter(2), h->GetXaxis()->GetXmax());
2365     fRGaus->SetRange(xMin, xMax);
2366     fRGaus->SetParLimits(1, xMin, xMax);
2367     h->Fit("fRGaus","NQR");
2368     
2369     mean = fRGaus->GetParameter(1);
2370     meanErr = fRGaus->GetParError(1);
2371     sigma = fRGaus->GetParameter(2);
2372     sigmaErr = fRGaus->GetParError(2);
2373     
2374   } else {
2375     
2376     Zoom(h);
2377     mean = h->GetMean();
2378     meanErr = h->GetMeanError();
2379     sigma = h->GetRMS();
2380     sigmaErr = h->GetRMSError();
2381     h->GetXaxis()->SetRange(0,0);
2382     
2383   }
2384   
2385   gMean->SetPoint(i, i+1, mean);
2386   gMean->SetPointError(i, 0., meanErr);
2387   
2388   if (fCorrectForSystematics) {
2389     Double_t s = TMath::Sqrt(mean*mean + sigma*sigma);
2390     sigmaErr = (s>0.) ? TMath::Sqrt(sigma*sigma*sigmaErr*sigmaErr + mean*mean*meanErr*meanErr) / s : 0.;
2391     sigma = s;
2392   }
2393   
2394   gSigma->SetPoint(i, i+1, sigma);
2395   gSigma->SetPointError(i, 0., sigmaErr);
2396   
2397 }
2398
2399 //________________________________________________________________________
2400 TCanvas* AliAnalysisTaskMuonPerformance::DrawVsAng(const char* name, const char* title, TH1* h1, TH2* h2)
2401 {
2402   /// generic function to draw histograms versus absorber angular region
2403   TCanvas* c = new TCanvas(name, title);
2404   c->cd();
2405   h1->DrawCopy();
2406   TH1D *proj1 = h2->ProjectionY(Form("%s_proj_0_2",h2->GetName()),1,2);
2407   proj1->SetLineColor(2);
2408   proj1->Draw("sames");
2409   TH1D *proj2 = h2->ProjectionY(Form("%s_proj_2_3",h2->GetName()),3,3);
2410   proj2->SetLineColor(4);
2411   proj2->Draw("sames");
2412   TH1D *proj3 = h2->ProjectionY(Form("%s__proj_3_10",h2->GetName()),4,10);
2413   proj3->SetLineColor(3);
2414   proj3->Draw("sames");
2415   return c;
2416 }
2417
2418 //________________________________________________________________________
2419 TCanvas* AliAnalysisTaskMuonPerformance::DrawVsPos(const char* name, const char* title, TH2* h1, TH2* h2, TH2* h3)
2420 {
2421   /// generic function to draw histograms versus position at absorber end
2422   TCanvas* c = new TCanvas(name, title);
2423   c->cd();
2424   h1->SetMarkerColor(2);
2425   h1->DrawCopy();
2426   h2->SetMarkerColor(4);
2427   h2->DrawCopy("sames");
2428   h3->SetMarkerColor(3);
2429   h3->DrawCopy("sames");
2430   return c;
2431 }
2432
2433 //________________________________________________________________________
2434 TCanvas* AliAnalysisTaskMuonPerformance::DrawFitLandauGausResPVsP(const char* name, const char* title,
2435                                                                   TH2* h, const Int_t nBins, const char* fitting)
2436 {
2437   /// generic function to draw and fit momentum residuals versus momentum
2438   
2439   static TF1* fLandauGaus = 0x0;
2440   if (!fLandauGaus) fLandauGaus = new TF1("fLandauGaus",langaufun,h->GetYaxis()->GetBinLowEdge(1),h->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1),4);
2441   
2442   TLegend* l = new TLegend(0.15,0.25,0.3,0.85);
2443   TCanvas* c = new TCanvas(name, title);
2444   c->cd();
2445   
2446   h->Sumw2();
2447   
2448   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
2449   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
2450     
2451     cout<<Form("\rFitting %s... %d/%d",fitting,i/rebinFactorX,nBins)<<flush;
2452     
2453     // draw projection
2454     TH1D* proj = h->ProjectionY(Form("%s_%d",h->GetName(),i/rebinFactorX),i-rebinFactorX+1,i);
2455     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
2456     proj->Draw((i==rebinFactorX)?"hist":"histsames");
2457     proj->SetLineColor(i/rebinFactorX);
2458     
2459     // first fit
2460     fLandauGaus->SetParameters(0.2,0.,1.,1.);
2461     fLandauGaus->SetLineColor(i/rebinFactorX);
2462     proj->Fit("fLandauGaus","WWNQ","sames");
2463     
2464     // rebin histo
2465     Double_t fwhm = fLandauGaus->GetParameter(0);
2466     Double_t sigma = fLandauGaus->GetParameter(3);
2467     Double_t sigmaP = TMath::Sqrt(sigma*sigma + fwhm*fwhm/(8.*log(2.)));
2468     Int_t rebin = TMath::Max(static_cast<Int_t>(0.5*sigmaP/proj->GetBinWidth(1)),1);
2469     while (proj->GetNbinsX()%rebin!=0) rebin--;
2470     proj->Rebin(rebin);
2471     proj->Scale(1./rebin);
2472     
2473     // second fit
2474     proj->Fit("fLandauGaus","Q","sames");
2475     
2476     // set label
2477     Double_t p = 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
2478     l->AddEntry(proj,Form("%5.1f GeV",p));
2479     
2480   }
2481   
2482   cout<<Form("\rFitting %s... %d/%d",fitting,nBins,nBins)<<endl;
2483   
2484   l->Draw("same");
2485   
2486   return c;
2487 }
2488
2489 //________________________________________________________________________
2490 TCanvas* AliAnalysisTaskMuonPerformance::DrawResPVsP(const char* name, const char* title, TH2* h, const Int_t nBins)
2491 {
2492   /// generic function to draw momentum residuals versus momentum
2493   
2494   TLegend* l = new TLegend(0.15,0.25,0.3,0.85);
2495   TCanvas* c = new TCanvas(name, title);
2496   c->cd();
2497   
2498   h->Sumw2();
2499   
2500   Int_t rebinFactorX = TMath::Max(h->GetNbinsX()/nBins, 1);
2501   for (Int_t i = rebinFactorX; i <= h->GetNbinsX(); i+=rebinFactorX) {
2502     
2503     // draw projection
2504     TH1D* proj = h->ProjectionY(Form("%s_%d",h->GetName(),i/rebinFactorX),i-rebinFactorX+1,i);
2505     if (proj->GetEntries() > 0) proj->Scale(1./proj->GetEntries());
2506     proj->Draw((i==rebinFactorX)?"hist":"histsames");
2507     proj->SetLineColor(i/rebinFactorX);
2508     proj->SetLineWidth(2);
2509     
2510     // set label
2511     Double_t p = 0.5 * (h->GetBinLowEdge(i-rebinFactorX+1) + h->GetBinLowEdge(i+1));
2512     l->AddEntry(proj,Form("%5.1f GeV",p));
2513     
2514   }
2515   
2516   l->Draw("same");
2517   
2518   return c;
2519 }
2520
2521 //________________________________________________________________________
2522 void AliAnalysisTaskMuonPerformance::Zoom(TH1* h, Double_t fractionCut)
2523 {
2524   /// Reduce the range of the histogram by removing a given fration of the statistic at each edge
2525   
2526   Double_t maxEventsCut = fractionCut * h->GetEntries();
2527   Int_t nBins = h->GetNbinsX();
2528   
2529   // set low edge  
2530   Int_t minBin;
2531   Double_t eventsCut = 0.;
2532   for (minBin = 1; minBin <= nBins; minBin++) {
2533     eventsCut += h->GetBinContent(minBin);
2534     if (eventsCut > maxEventsCut) break;
2535   }
2536   
2537   // set high edge
2538   Int_t maxBin;
2539   eventsCut = 0.;
2540   for (maxBin = nBins; maxBin >= 1; maxBin--) {
2541     eventsCut += h->GetBinContent(maxBin);
2542     if (eventsCut > maxEventsCut) break;
2543   }
2544   
2545   // set new axis range
2546   h->GetXaxis()->SetRange(--minBin, ++maxBin);
2547 }
2548
2549 //________________________________________________________________________
2550 void AliAnalysisTaskMuonPerformance::FillEffHistos(AliCFEffGrid* efficiency, const char* suffix, TObjArray* list)
2551 {
2552   /// Compute efficiency histograms and save them to the given list
2553   
2554   TH1* auxHisto = efficiency->Project(kVarPt);
2555   auxHisto->SetName(Form("effVsPt_%s",suffix));
2556   list->AddLast(auxHisto);
2557   
2558   auxHisto = efficiency->Project(kVarEta);
2559   auxHisto->SetName(Form("effVsEta_%s",suffix));
2560   list->AddLast(auxHisto);
2561   
2562   auxHisto = efficiency->Project(kVarPhi);
2563   auxHisto->SetName(Form("effVsPhi_%s",suffix));
2564   list->AddLast(auxHisto);
2565   
2566   auxHisto = efficiency->Project(kVarCent);
2567   auxHisto->SetName(Form("effVsCent_%s",suffix));
2568   list->AddLast(auxHisto);
2569   
2570 }
2571