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