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