changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisNetParticleHelper.cxx
CommitLineData
cb68eb1d 1//-*- Mode: C++ -*-
2
3#include "TMath.h"
4#include "TAxis.h"
5#include "TSystem.h"
6#include "TFile.h"
7#include "TPRegexp.h"
8
9#include "AliStack.h"
10#include "AliMCEvent.h"
11#include "AliMCParticle.h"
12#include "AliESDtrackCuts.h"
13#include "AliESDInputHandler.h"
14#include "AliESDpid.h"
9be43c8e 15#include "AliAODInputHandler.h"
16#include "AliAODEvent.h"
17#include "AliAODMCParticle.h"
cb68eb1d 18#include "AliCentrality.h"
19#include "AliTracker.h"
20
21#include "AliAnalysisNetParticleHelper.h"
22
23using namespace std;
24
4918c45f 25/**
a2ddc3d0 26 * Class for NetParticle Distributions
4918c45f 27 * -- Helper class for net particle istributions
28 * Authors: Jochen Thaeder <jochen@thaeder.de>
29 * Michael Weber <m.weber@cern.ch>
30 */
cb68eb1d 31
32ClassImp(AliAnalysisNetParticleHelper)
33
34/*
35 * ---------------------------------------------------------------------------------
3aa68b92 36 * particle names
37 * ---------------------------------------------------------------------------------
38 */
39
478c95cf 40 // MW make fgk ... static const
3aa68b92 41 const Char_t* aPartNames[AliPID::kSPECIES][2] = {
42 {"ele", "posi"},
43 {"mubar", "mu"},
44 {"pibar", "pi"},
45 {"kbar", "k"},
46 {"pbar", "p"}
47 };
48
49 const Char_t* aPartTitles[AliPID::kSPECIES][2] = {
50 {"Electron", "Positron"},
51 {"Anti-Muon", "Muon"},
52 {"Anti-Pion", "Proton"},
53 {"Anti-Kaon", "Kaon"},
54 {"Anti-Proton", "Proton"}
55 };
56
57 const Char_t* aPartTitlesLatex[AliPID::kSPECIES][2] = {
58 {"e^{-}", "e^{+}" },
59 {"#mu^{-}", "#mu^{+}" },
60 {"#pi^{-}", "#pi^{+}" },
61 {"K^{-}", "K^{+}" },
62 {"#bar{p}", "p"}
63 };
64
3aa68b92 65/*
66 * ---------------------------------------------------------------------------------
cb68eb1d 67 * Constructor / Destructor
68 * ---------------------------------------------------------------------------------
69 */
70
71//________________________________________________________________________
72AliAnalysisNetParticleHelper::AliAnalysisNetParticleHelper() :
478c95cf 73 fModeDistCreation(0),
74
75 fInputEventHandler(NULL),
cb68eb1d 76 fPIDResponse(NULL),
77 fESD(NULL),
a2ddc3d0 78 fESDTrackCuts(NULL),
9be43c8e 79 fAOD(NULL),
a2ddc3d0 80 fAODtrackCutBit(1024),
81 fIsMC(kFALSE),
cb68eb1d 82 fMCEvent(NULL),
83 fStack(NULL),
84
85 fCentralityBin(-1),
86 fCentralityPercentile(-1.),
87
88 fCentralityBinMax(7),
89 fVertexZMax(10.),
90 fRapidityMax(0.5),
4918c45f 91 fPhiMin(0.),
92 fPhiMax(TMath::TwoPi()),
cb68eb1d 93 fMinTrackLengthMC(70.),
94 fNSigmaMaxCdd(3.),
95 fNSigmaMaxCzz(3.),
96
97 fParticleSpecies(AliPID::kProton),
cb68eb1d 98
478c95cf 99 fUsePID(kTRUE),
a2ddc3d0 100 fPIDStrategy(0),
101 fNSigmaMaxITS(4.),
102 fNSigmaMaxTPC(4.),
103 fNSigmaMaxTPClow(3.),
104 fNSigmaMaxTOF(4.),
105 fMinPtForTOFRequired(0.69),
106 fMaxPtForTPClow(0.69),
cb68eb1d 107
108 fHEventStat0(NULL),
109 fHEventStat1(NULL),
110 fHEventStatMax(6),
111
112 fHTriggerStat(NULL),
113 fNTriggers(5),
114
115 fHCentralityStat(NULL),
116 fNCentralityBins(10),
117
a2ddc3d0 118 fNSubSamples(20),
119 fSubSampleIdx(0),
120 fRandom(NULL) {
cb68eb1d 121 // Constructor
122
123 AliLog::SetClassDebugLevel("AliAnalysisNetParticleHelper",10);
124}
125
4918c45f 126const Float_t AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap = 0.075;
a2ddc3d0 127const Float_t AliAnalysisNetParticleHelper::fgkfHistBinWitdthPt = 0.3; // 0.08 // 300 MeV // was 80 MeV
4918c45f 128
129const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeCent[] = {-0.5, 8.5};
130const Int_t AliAnalysisNetParticleHelper::fgkfHistNBinsCent = 9 ;
131
132const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeEta[] = {-0.9, 0.9};
133const Int_t AliAnalysisNetParticleHelper::fgkfHistNBinsEta = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangeEta[1] -
134 AliAnalysisNetParticleHelper::fgkfHistRangeEta[0]) /
135 AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap) +1;
136
137const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeRap[] = {-0.5, 0.5};
138const Int_t AliAnalysisNetParticleHelper::fgkfHistNBinsRap = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangeRap[1] -
139 AliAnalysisNetParticleHelper::fgkfHistRangeRap[0]) /
140 AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap) +1;
141
71aa6041 142const Float_t AliAnalysisNetParticleHelper::fgkfHistRangePhi[] = {0.0, static_cast<Float_t>(TMath::TwoPi())};
4918c45f 143const Int_t AliAnalysisNetParticleHelper::fgkfHistNBinsPhi = 42 ;
144
a2ddc3d0 145const Float_t AliAnalysisNetParticleHelper::fgkfHistRangePt[] = {0.2, 2.9}; // {0.2, 5.}; // was {0.3, 2.22}
4918c45f 146const Int_t AliAnalysisNetParticleHelper::fgkfHistNBinsPt = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangePt[1] -
147 AliAnalysisNetParticleHelper::fgkfHistRangePt[0]) /
148 AliAnalysisNetParticleHelper::fgkfHistBinWitdthPt);
149
150const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeSign[] = {-1.5, 1.5};
151const Int_t AliAnalysisNetParticleHelper::fgkfHistNBinsSign = 3;
152
153const Char_t* AliAnalysisNetParticleHelper::fgkEventNames[] = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,90]%"};
154const Char_t* AliAnalysisNetParticleHelper::fgkCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
155const Char_t* AliAnalysisNetParticleHelper::fgkTriggerNames[] = {"kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
156const Char_t* AliAnalysisNetParticleHelper::fgkCentralityNames[] = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%",
157 "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
158
cb68eb1d 159//________________________________________________________________________
160AliAnalysisNetParticleHelper::~AliAnalysisNetParticleHelper() {
161 // Destructor
162
a2ddc3d0 163 if (fRandom)
164 delete fRandom;
165
cb68eb1d 166 return;
167}
168
169/*
170 * ---------------------------------------------------------------------------------
3aa68b92 171 * Setter
172 * ---------------------------------------------------------------------------------
173 */
174
175//________________________________________________________________________
5de3d4d1 176void AliAnalysisNetParticleHelper::SetPhiRange(Float_t f1, Float_t f2) {
177 // -- Set phi range and adopt to phi-histogram
178
179 fPhiMin = f1;
180 fPhiMax = (f1 < f2) ? f2 : f2+TMath::TwoPi();
a2ddc3d0 181
182 Float_t phiMin = fPhiMin;
183 Float_t phiMax = fPhiMax;
5de3d4d1 184
185 // -- Update Ranges
186 Float_t binWidth = (AliAnalysisNetParticleHelper::fgkfHistRangePhi[1] - AliAnalysisNetParticleHelper::fgkfHistRangePhi[0]) /
187 Float_t(AliAnalysisNetParticleHelper::fgkfHistNBinsPhi);
188
189 Float_t lowEdge = AliAnalysisNetParticleHelper::fgkfHistRangePhi[0] - binWidth;
190 Float_t highEdge = AliAnalysisNetParticleHelper::fgkfHistRangePhi[0];
191
192 for (Int_t ii = 1; ii <= AliAnalysisNetParticleHelper::fgkfHistNBinsPhi; ++ii) {
193 lowEdge += binWidth;
194 highEdge += binWidth;
195
a2ddc3d0 196 if (phiMin >= lowEdge && phiMin < highEdge )
197 phiMin = lowEdge;
198 if (phiMax > lowEdge && phiMax <= highEdge )
199 phiMax = highEdge;
5de3d4d1 200 }
201
a2ddc3d0 202 printf(">>>> Update Phi Range : [%f,%f] -> [%f,%f]\n", fPhiMin, fPhiMax, phiMin, phiMax);
203 fPhiMin = phiMin;
204 fPhiMax = phiMax;
5de3d4d1 205}
206
207
208//________________________________________________________________________
3aa68b92 209void AliAnalysisNetParticleHelper::SetParticleSpecies(AliPID::EParticleType pid) {
210 // -- Set particle species (ID, Name, Title, Title LATEX)
8c6128b7 211
212 if ( Int_t(pid) < 0 || Int_t(pid) >= AliPID::kSPECIES) {
3aa68b92 213 AliWarning("Particle ID not in AliPID::kSPECIES --> Set to protons");
214 pid = AliPID::kProton;
215 }
478c95cf 216
217 fParticleSpecies = pid;
218
3aa68b92 219 for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
220 fPartName[idxPart] = aPartNames[fParticleSpecies][idxPart];
221 fPartTitle[idxPart] = aPartTitles[fParticleSpecies][idxPart];
222 fPartTitleLatex[idxPart] = aPartTitlesLatex[fParticleSpecies][idxPart];
223 }
224}
5de3d4d1 225
478c95cf 226//________________________________________________________________________
227void AliAnalysisNetParticleHelper::SetUsePID(Bool_t usePID) {
228 // -- Set usage of PID
229 // > if turn off, set charge types (ID, Name, Title, Title LATEX)
230
231 fUsePID = usePID;
232
233 if (!usePID) {
234 fParticleSpecies = AliPID::kUnknown;
3aa68b92 235
478c95cf 236 fPartName[0] = "neg";
237 fPartName[1] = "pos";
238 fPartTitle[0] = "Negative";
239 fPartTitle[1] = "Positive";
240 fPartTitleLatex[0] = "Negative";
241 fPartTitleLatex[1] = "Positive";
242 }
243}
3aa68b92 244
245/*
246 * ---------------------------------------------------------------------------------
247 * Getter
248 * ---------------------------------------------------------------------------------
249 */
250
251//________________________________________________________________________
252TString AliAnalysisNetParticleHelper::GetParticleName(Int_t idxPart) {
253 // -- Get particle Name
254
255 if( idxPart != 0 && idxPart != 1){
256 AliWarning("Particle type not known --> Set to antiparticles");
257 idxPart = 0;
258 }
259
260 return fPartName[idxPart];
261}
262
263//________________________________________________________________________
264TString AliAnalysisNetParticleHelper::GetParticleTitle(Int_t idxPart) {
265 // -- Get particle Title
266
267 if( idxPart != 0 && idxPart != 1){
268 AliWarning("Particle type not known --> Set to antiparticles");
269 idxPart = 0;
270 }
271
272 return fPartTitle[idxPart];
273}
274
275//________________________________________________________________________
276TString AliAnalysisNetParticleHelper::GetParticleTitleLatex(Int_t idxPart) {
277 // -- Get particle Title LATEX
278
279 if( idxPart != 0 && idxPart != 1){
280 AliWarning("Particle type not known --> Set to antiparticles");
281 idxPart = 0;
282 }
283
284 return fPartTitleLatex[idxPart];
285}
286
3aa68b92 287/*
288 * ---------------------------------------------------------------------------------
cb68eb1d 289 * Public Methods
290 * ---------------------------------------------------------------------------------
291 */
292
293//________________________________________________________________________
a2ddc3d0 294Int_t AliAnalysisNetParticleHelper::Initialize(AliESDtrackCuts *cuts, Bool_t isMC, Int_t trackCutBit, Int_t modeDistCreation) {
cb68eb1d 295 // -- Initialize helper
296
297 Int_t iResult = 0;
a2ddc3d0 298
299 // -- ESD track cuts
300 fESDTrackCuts = cuts;
301
302 // -- Is MC
303 fIsMC = isMC;
304
305 // -- AOD track filter bit
306 fAODtrackCutBit = trackCutBit;
307
308 // -- mode Distribution creation
478c95cf 309 fModeDistCreation = modeDistCreation;
cb68eb1d 310
311 // -- Setup event cut statistics
312 InitializeEventStats();
313
314 // -- Setup trigger statistics
315 InitializeTriggerStats();
316
317 // -- Setup centrality statistics
318 InitializeCentralityStats();
319
a2ddc3d0 320 // -- PRINT PID Strategy
321 // 0 : TPC(TPClow+TPCHigh)
322 // 1 : ITS
323 // 2 : TOF
324 // 3 : ITS+TPC(TPClow+TPCHigh)
325 // 4 : TPC(TPClow+TPCHigh)+TOF
326 // 5 : TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
327 // 6 : TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
328 // 7 : TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
329 // 8 : TPC(TPClow+TPCHigh)+ITS+TOF
330 printf(">>>> USE PID %d || PID STRATEGY: %d || sigmaMax: ITS %.2f TPC %.2f TOF %.2f \n", fUsePID, fPIDStrategy, fNSigmaMaxITS, fNSigmaMaxTPC, fNSigmaMaxTOF);
331
332 // -- Initialize random number generator
333 fRandom = new TRandom3();
334 fRandom->SetSeed();
335
cb68eb1d 336 return iResult;
337}
338
339//________________________________________________________________________
9be43c8e 340Int_t AliAnalysisNetParticleHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
cb68eb1d 341 // -- Setup Event
478c95cf 342
cb68eb1d 343 // -- Get ESD objects
9be43c8e 344 if(esdHandler){
478c95cf 345 fInputEventHandler = static_cast<AliInputEventHandler*>(esdHandler);
346 fESD = dynamic_cast<AliESDEvent*>(fInputEventHandler->GetEvent());
36ef73d2 347 if (!fESD) {
348 AliError("ESD event handler not available");
349 return -1;
350 }
9be43c8e 351 }
352
353 // -- Get AOD objects
354 else if(aodHandler){
478c95cf 355 fInputEventHandler = static_cast<AliInputEventHandler*>(aodHandler);
356 fAOD = dynamic_cast<AliAODEvent*>(fInputEventHandler->GetEvent());
36ef73d2 357 if (!fAOD) {
358 AliError("AOD event handler not available");
359 return -1;
360 }
9be43c8e 361 }
cb68eb1d 362
478c95cf 363 // -- Get Common objects
364 fPIDResponse = fInputEventHandler->GetPIDResponse();
365
cb68eb1d 366 // -- Get MC objects
367 fMCEvent = mcEvent;
368 if (fMCEvent)
369 fStack = fMCEvent->Stack();
370
371 // -- Get event centrality
372 // > 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
373 // > 0 1 2 3 4 5 6 7 8 9
374
9be43c8e 375 AliCentrality *centrality = NULL;
376
377 if(esdHandler)
378 centrality = fESD->GetCentrality();
379 else if(aodHandler)
380 centrality = fAOD->GetHeader()->GetCentralityP();
381
36ef73d2 382 if (!centrality) {
383 AliError("Centrality not available");
384 return -1;
385 }
8c6128b7 386
cb68eb1d 387 Int_t centBin = centrality->GetCentralityClass10("V0M");
388 if (centBin == 0)
389 fCentralityBin = centrality->GetCentralityClass5("V0M");
390 else if (centBin == 10 || centBin == -1.)
391 fCentralityBin = -1;
392 else if (centBin > 0 && centBin < fNCentralityBins)
393 fCentralityBin = centBin + 1;
394 else
395 fCentralityBin = -2;
478c95cf 396
cb68eb1d 397 // -- Stay within the max centrality bin
398 if (fCentralityBin >= fCentralityBinMax)
399 fCentralityBin = -2;
478c95cf 400
cb68eb1d 401 fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
402
a2ddc3d0 403 // -- Get current subsample idx
404 fSubSampleIdx = fRandom->Integer(fNSubSamples);
cb68eb1d 405
406 return 0;
407}
408
409/*
410 * ---------------------------------------------------------------------------------
411 * Event / Trigger Statistics
412 * ---------------------------------------------------------------------------------
413 */
414
415//________________________________________________________________________
416Bool_t AliAnalysisNetParticleHelper::IsEventTriggered() {
417 // -- Check if Event is triggered and fill Trigger Histogram
418
419 Bool_t *aTriggerFired = new Bool_t[fNTriggers];
420 for (Int_t ii = 0; ii < fNTriggers; ++ii)
421 aTriggerFired[ii] = kFALSE;
422
478c95cf 423 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kMB)) aTriggerFired[0] = kTRUE;
424 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kCentral)) aTriggerFired[1] = kTRUE;
425 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
426 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEJE)) aTriggerFired[3] = kTRUE;
427 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEGA)) aTriggerFired[4] = kTRUE;
9be43c8e 428
cb68eb1d 429 Bool_t isTriggered = kFALSE;
430
431 for (Int_t ii=0; ii<fNTriggers; ++ii) {
432 if(aTriggerFired[ii]) {
433 isTriggered = kTRUE;
434 fHTriggerStat->Fill(ii);
435 }
436 }
437
438 delete[] aTriggerFired;
439
440 return isTriggered;
441}
442
443//________________________________________________________________________
444Bool_t AliAnalysisNetParticleHelper::IsEventRejected() {
445 // -- Evaluate event statistics histograms
446
447 Int_t *aEventCuts = new Int_t[fHEventStatMax];
448 // set aEventCuts[ii] to 1 in case of reject
449
450 for (Int_t ii=0;ii<fHEventStatMax; ++ii)
451 aEventCuts[ii] = 0;
452
453 Int_t iCut = 0;
454
455 // -- 0 - Before Physics Selection
456 aEventCuts[iCut] = 0;
457
458 // -- 1 - No Trigger fired
459 ++iCut;
460 if (!IsEventTriggered())
461 aEventCuts[iCut] = 1;
462
463 // -- 2 - No Vertex
464 ++iCut;
9be43c8e 465 const AliESDVertex* vtxESD = NULL;
466 const AliAODVertex* vtxAOD = NULL;
467 if (fESD){
468 vtxESD = fESD->GetPrimaryVertexTracks();
469 if (!vtxESD)
470 aEventCuts[iCut] = 1;
471 }
472 else if (fAOD){
473 vtxAOD = fAOD->GetPrimaryVertex();
474 if (!vtxAOD)
475 aEventCuts[iCut] = 1;
476 }
cb68eb1d 477
478 // -- 3 - Vertex z outside cut window
479 ++iCut;
480 if (vtxESD){
e690d4d0 481 if(TMath::Abs(vtxESD->GetZ()) > fVertexZMax)
cb68eb1d 482 aEventCuts[iCut] = 1;
483 }
9be43c8e 484 else if(vtxAOD){
485 if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax)
486 aEventCuts[iCut] = 1;
487 }
cb68eb1d 488 else
489 aEventCuts[iCut] = 1;
478c95cf 490
cb68eb1d 491 // -- 4 - Centrality = -1 (no centrality or not hadronic)
492 ++iCut;
493 if(fCentralityBin == -1.)
494 aEventCuts[iCut] = 1;
495
496 // -- 5 - Centrality < fCentralityMax
497 ++iCut;
498 if(fCentralityBin == -2.)
499 aEventCuts[iCut] = 1;
500
501 // -- Fill statistics / reject event
502 Bool_t isRejected = FillEventStats(aEventCuts);
503
504 // -- Cleanup
505 delete[] aEventCuts;
506
507 return isRejected;
508}
509
510/*
511 * ---------------------------------------------------------------------------------
512 * Accept Particle Methods - private
513 * ---------------------------------------------------------------------------------
514 */
515
516//________________________________________________________________________
fbb73c19 517Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC) {
9be43c8e 518 // -- Check if MC particle is accepted for basic parameters
519
520 if (!particle)
521 return kFALSE;
522
9be43c8e 523 // -- check if charged
0c47a116 524 if (particle->Charge() == 0.0)
9be43c8e 525 return kFALSE;
0c47a116 526
fbb73c19 527 // -- check if physical primary - ESD
528 if (fESD) {
529 if(!fStack->IsPhysicalPrimary(idxMC))
530 return kFALSE;
531 }
532 // -- check if physical primary - AOD
533 else {
534 if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
535 return kFALSE;
536 }
cb68eb1d 537
cb68eb1d 538 return kTRUE;
539}
fbb73c19 540
9be43c8e 541//________________________________________________________________________
fbb73c19 542Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC) {
9be43c8e 543 // -- Check if MC particle is accepted for basic parameters
544
545 if (!particle)
546 return kFALSE;
cb68eb1d 547
9be43c8e 548 // -- check if charged
0c47a116 549 if (particle->Charge() != 0.0)
9be43c8e 550 return kFALSE;
0c47a116 551
fbb73c19 552 // -- check if physical primary - ESD
553 if (fESD) {
554 if(!fStack->IsPhysicalPrimary(idxMC))
555 return kFALSE;
556 }
557 // -- check if physical primary - AOD
558 else {
559 if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
560 return kFALSE;
8c6128b7 561 }
8c6128b7 562
563 return kTRUE;
564}
565
566//________________________________________________________________________
fbb73c19 567Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP) {
568 // -- Check if particle is accepted
8c6128b7 569 // > in rapidity
4918c45f 570 // > if no pid : return kTRUE, yP = eta
cb68eb1d 571 // > return 0 if not accepted
572
478c95cf 573 if (!fUsePID) {
478c95cf 574 yP = particle->Eta();
4918c45f 575 return kTRUE;
478c95cf 576 }
577
cb68eb1d 578 Double_t mP = AliPID::ParticleMass(fParticleSpecies);
579
580 // -- Calculate rapidities and kinematics
581 Double_t p = particle->P();
582 Double_t pz = particle->Pz();
583
584 Double_t eP = TMath::Sqrt(p*p + mP*mP);
585 yP = 0.5 * TMath::Log((eP + pz) / (eP - pz));
586
587 // -- Check Rapidity window
588 if (TMath::Abs(yP) > fRapidityMax)
589 return kFALSE;
590
591 return kTRUE;
592}
593
4918c45f 594//________________________________________________________________________
fbb73c19 595Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedPhi(AliVParticle *particle) {
4918c45f 596 // -- Check if particle is accepted
597 // > in phi
598 // > return 0 if not accepted
599
600 if (particle->Phi() > fPhiMin && particle->Phi() <= fPhiMax)
601 return kTRUE;
602 else if (particle->Phi() < fPhiMin && (particle->Phi() + TMath::TwoPi()) <= fPhiMax)
603 return kTRUE;
604 else
605 return kFALSE;
606}
607
cb68eb1d 608//_____________________________________________________________________________
609Bool_t AliAnalysisNetParticleHelper::IsParticleFindable(Int_t label) {
610 // -- Check if MC particle is findable tracks
611
612 AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
613 if(!mcParticle)
614 return kFALSE;
615
616 Int_t counter;
617 Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
618
619 return (tpcTrackLength > fMinTrackLengthMC);
620}
621
622/*
623 * ---------------------------------------------------------------------------------
624 * Accept Track Methods - public
625 * ---------------------------------------------------------------------------------
626 */
9be43c8e 627
cb68eb1d 628//________________________________________________________________________
478c95cf 629Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliVTrack* track) {
cb68eb1d 630 // -- Check if track is accepted
631 // > for basic parameters
9be43c8e 632
cb68eb1d 633 if (!track)
634 return kFALSE;
635
636 if (track->Charge() == 0)
637 return kFALSE;
638
cb68eb1d 639 return kTRUE;
9be43c8e 640}
641
642//________________________________________________________________________
9be43c8e 643Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP) {
cb68eb1d 644 // -- Check if track is accepted
645 // > in rapidity
4918c45f 646 // > if no pid : return kTRUE
cb68eb1d 647 // > return 0 if not accepted
648
478c95cf 649 if (!fUsePID) {
650 yP = track->Eta();
651 return kTRUE;
652 }
653
cb68eb1d 654 Double_t mP = AliPID::ParticleMass(fParticleSpecies);
655
656 // -- Calculate rapidities and kinematics
657 Double_t pvec[3];
658 track->GetPxPyPz(pvec);
659
9be43c8e 660 Double_t p = track->P();
cb68eb1d 661 Double_t eP = TMath::Sqrt(p*p + mP*mP);
662 yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
663
664 // -- Check Rapidity window
665 if (TMath::Abs(yP) > fRapidityMax)
666 return kFALSE;
667
668 return kTRUE;
669}
670
671//________________________________________________________________________
478c95cf 672Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedDCA(AliVTrack *vTrack) {
673 // -- Check if track is accepted - ONLY FOR ESDs so far
cb68eb1d 674 // > for DCA, if both SPD layers have hits
478c95cf 675 // > For now only Implemented for ESDs
cb68eb1d 676
677 Bool_t isAccepted = kTRUE;
678
478c95cf 679 if (!fESD)
680 return isAccepted;
681
682 AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
683 if (!track)
684 return kFALSE;
685
cb68eb1d 686 // -- Get nHits SPD
687 if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
688
689 // -- Get DCA nSigmas
690 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
691 track->GetImpactParameters(dca,cov);
692
693 Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
694 Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
695
696 if (fNSigmaMaxCdd != 0.) {
697 if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
698 isAccepted = kFALSE;
699 }
700
701 if (fNSigmaMaxCzz != 0.) {
702 if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
703 isAccepted = kFALSE;
704 }
705 }
706
707 return isAccepted;
708}
709
710//________________________________________________________________________
9be43c8e 711Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid) {
cb68eb1d 712 // -- Check if track is accepted
a2ddc3d0 713 // > provides ITS, TPC and TOF nSigmas to the argument
cb68eb1d 714
a2ddc3d0 715 Bool_t isAcceptedITS = kFALSE;
716 Bool_t isAcceptedTPC = kFALSE;
717 Bool_t isAcceptedTPClow = kFALSE;
718 Bool_t isAcceptedTOF = kFALSE;
719 Bool_t isAccepted = kFALSE;
478c95cf 720
721 // -- In case not PID is used
722 if (!fUsePID) {
723 pid[0] = 10.;
724 pid[1] = 10.;
a2ddc3d0 725 pid[2] = 10.;
478c95cf 726 return kTRUE;
727 }
cb68eb1d 728
a2ddc3d0 729 // -- Get PID with ITS and check
730 if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kITS, track, fParticleSpecies, pid[0]) == AliPIDResponse::kDetPidOk) {
731 if (TMath::Abs(pid[0]) < fNSigmaMaxITS)
732 isAcceptedITS = kTRUE;
733 }
734
735 // -- Get PID with TPC and check
736 if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, fParticleSpecies, pid[1]) == AliPIDResponse::kDetPidOk) {
737 if (TMath::Abs(pid[1]) < fNSigmaMaxTPC)
738 isAcceptedTPC = kTRUE;
739 if (TMath::Abs(pid[1]) < fNSigmaMaxTPClow)
740 isAcceptedTPClow = kTRUE;
741 if (track->Pt() < fMaxPtForTPClow)
742 isAcceptedTPC = isAcceptedTPClow;
743 }
cb68eb1d 744
a2ddc3d0 745 // -- Get PID with TOF and check
746 Bool_t hasPIDTOF = kFALSE;
747 if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, fParticleSpecies, pid[2]) == AliPIDResponse::kDetPidOk) {
748 hasPIDTOF = kTRUE;
749 if (TMath::Abs(pid[2]) < fNSigmaMaxTOF)
750 isAcceptedTOF = kTRUE;
751 }
752
753 // -- Check TOF missmatch for MC
cb68eb1d 754
a2ddc3d0 755 //if (ESD)
756 if (fIsMC && isAcceptedTOF) {
757 Int_t tofLabel[3];
758 // AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
759 // TODO add code for AOD
cb68eb1d 760
a2ddc3d0 761 (dynamic_cast<AliESDtrack*>(track))->GetTOFLabel(tofLabel);
762
763 Bool_t hasMatchTOF = kTRUE;
764 if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0)
765 hasMatchTOF = kFALSE;
766
767 TParticle *matchedTrack = fStack->Particle(TMath::Abs(tofLabel[0]));
768 if (TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
769 hasMatchTOF = kTRUE;
770
771 isAcceptedTOF = hasMatchTOF;
772 }
773
774 // 0 : TPC(TPClow+TPCHigh)
775 // 1 : ITS
776 // 2 : TOF
777 // 3 : ITS+TPC(TPClow+TPCHigh)
778 // 4 : TPC(TPClow+TPCHigh)+TOF
779 // 5 : TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
780 // 6 : TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
781 // 7 : TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
782 // 8 : TPC(TPClow+TPCHigh)+ITS+TOF
783 if (fPIDStrategy == 0) { // TPC PID
784 isAccepted = isAcceptedTPC;
785 }
786 else if (fPIDStrategy == 1) { // ITS PID
787 isAccepted = isAcceptedITS;
788 }
789 else if (fPIDStrategy == 2) { // TOF PID
790 isAccepted = isAcceptedTOF;
791 }
792 else if (fPIDStrategy == 3) { // TPC+ITS PID
793 isAccepted = isAcceptedTPC && isAcceptedITS;
794 }
795 else if (fPIDStrategy == 4) { // TPC+TOF PID
796 isAccepted = isAcceptedTPC && isAcceptedTOF;
797 }
798 else if (fPIDStrategy == 5) { // TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
799 if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
800 isAcceptedTOF = kTRUE;
801 isAccepted = isAcceptedTPC && isAcceptedTOF;
802 }
803 else if (fPIDStrategy == 6) { // ITS+TPC+TOF PID -- TOF only for those tracks which have TOF information
804 isAccepted = isAcceptedTPC && isAcceptedITS;
805 if (hasPIDTOF)
806 isAccepted = isAccepted && isAcceptedTOF;
807 }
808 else if (fPIDStrategy == 7) { // ITS+TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
809 if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
810 isAcceptedTOF = kTRUE;
811 isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
812 }
813 else if (fPIDStrategy == 8) { // ITS+TPC+TOF PID
814 isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
815 }
cb68eb1d 816
a2ddc3d0 817 return isAccepted;
cb68eb1d 818}
819
4918c45f 820//________________________________________________________________________
821Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPhi(AliVTrack *track) {
822 // -- Check if track is accepted
823 // > in phi
824 // > return 0 if not accepted
825
826 if (track->Phi() > fPhiMin && track->Phi() <= fPhiMax)
827 return kTRUE;
828 else if (track->Phi() < fPhiMin && (track->Phi() + TMath::TwoPi()) <= fPhiMax)
829 return kTRUE;
830 else
831 return kFALSE;
832}
833
cb68eb1d 834/*
835 * ---------------------------------------------------------------------------------
836 * Helper Methods
837 * ---------------------------------------------------------------------------------
838 */
839
840//________________________________________________________________________
a2ddc3d0 841void AliAnalysisNetParticleHelper::BinLogAxis(const THnBase *hn, Int_t axisNumber, AliESDtrackCuts* cuts) {
cb68eb1d 842 // -- Method for the correct logarithmic binning of histograms
4918c45f 843 // -- and update fMinPtForTOFRequired using the logarithmic scale
a2ddc3d0 844
845 AliESDtrackCuts* esdTrackCuts = (cuts) ? cuts : fESDTrackCuts;
4918c45f 846
847 // -- Make logarithmic binning
848 TAxis *axis = hn->GetAxis(axisNumber);
849 Int_t nBins = axis->GetNbins();
cb68eb1d 850
851 Double_t from = axis->GetXmin();
852 Double_t to = axis->GetXmax();
4918c45f 853 Double_t *newBins = new Double_t[nBins + 1];
cb68eb1d 854
855 newBins[0] = from;
4918c45f 856 Double_t factor = TMath::Power(to/from, 1./nBins);
cb68eb1d 857
4918c45f 858 for (int ii = 1; ii <= nBins; ii++)
859 newBins[ii] = factor * newBins[ii-1];
478c95cf 860
4918c45f 861 axis->Set(nBins, newBins);
478c95cf 862
cb68eb1d 863 delete [] newBins;
4918c45f 864
a2ddc3d0 865 // -- Update Ranges
866 // ------------------
867 Float_t ptRange[2];
868 Float_t oldPtRange[2];
869 esdTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
870 esdTrackCuts->GetPtRange(oldPtRange[0],oldPtRange[1]);
871
872 Float_t minPtForTOFRequired = fMinPtForTOFRequired;
873 Float_t maxPtForTPClow = fMaxPtForTPClow;
874
875 // -- Update minPt Cut
876 Int_t bin = axis->FindBin(ptRange[0]+10e-6);
877 ptRange[0] = axis->GetBinLowEdge(bin);
878
879 // -- Update maxPt Cut
880 bin = axis->FindBin(ptRange[1]-10e-6);
881 ptRange[1] = axis->GetBinUpEdge(bin);
882
883 if (ptRange[0] != oldPtRange[0] || ptRange[1] != oldPtRange[1]) {
884 printf(">>>> Update Pt Range : [%f,%f] -> [%f,%f]\n", oldPtRange[0], oldPtRange[1], ptRange[0], ptRange[1]);
885 esdTrackCuts->SetPtRange(ptRange[0],ptRange[1]);
886 }
887
4918c45f 888 // -- Update MinPtForTOFRequired
a2ddc3d0 889 bin = axis->FindBin(minPtForTOFRequired-10e-6);
890 minPtForTOFRequired = axis->GetBinUpEdge(bin);
891
892 if (minPtForTOFRequired != fMinPtForTOFRequired) {
893 printf(">>>> Update Min Pt for TOF : %f -> %f\n", fMinPtForTOFRequired, minPtForTOFRequired);
894 fMinPtForTOFRequired = minPtForTOFRequired;
895 }
896
897 // -- Update MaxPtForTPClow
898 bin = axis->FindBin(maxPtForTPClow-10e-6);
899 maxPtForTPClow = axis->GetBinUpEdge(bin);
900
901 if (maxPtForTPClow != fMaxPtForTPClow) {
902 printf(">>>> Update Max Pt for TPC Low : %f -> %f\n", fMaxPtForTPClow, maxPtForTPClow);
903 fMaxPtForTPClow = maxPtForTPClow;
904 }
cb68eb1d 905}
906
907///////////////////////////////////////////////////////////////////////////////////
908
909/*
910 * ---------------------------------------------------------------------------------
911 * Initialize - Private
912 * ---------------------------------------------------------------------------------
913 */
914
915//________________________________________________________________________
916void AliAnalysisNetParticleHelper::InitializeEventStats() {
917 // -- Initialize event statistics histograms
918
a2ddc3d0 919 fHEventStat0 = new TH1F("hEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
920 fHEventStat1 = new TH1F("hEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
cb68eb1d 921
922 for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
4918c45f 923 fHEventStat0->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkEventNames[ii]);
924 fHEventStat1->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkEventNames[ii]);
cb68eb1d 925 }
4918c45f 926
927 fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliAnalysisNetParticleHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
928 fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliAnalysisNetParticleHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
cb68eb1d 929}
930
931//________________________________________________________________________
932void AliAnalysisNetParticleHelper::InitializeTriggerStats() {
933 // -- Initialize trigger statistics histograms
934
a2ddc3d0 935 fHTriggerStat = new TH1F("hTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
cb68eb1d 936
937 for ( Int_t ii=0; ii < fNTriggers; ii++ )
4918c45f 938 fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkTriggerNames[ii]);
cb68eb1d 939}
940
941//________________________________________________________________________
942void AliAnalysisNetParticleHelper::InitializeCentralityStats() {
943 // -- Initialize trigger statistics histograms
944
a2ddc3d0 945 fHCentralityStat = new TH1F("hCentralityStat","Centrality statistics;Centrality Bins;Events",
cb68eb1d 946 fNCentralityBins,-0.5,fNCentralityBins-0.5);
947
948 for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
4918c45f 949 fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkCentralityNames[ii]);
cb68eb1d 950}
951
cb68eb1d 952/*
953 * ---------------------------------------------------------------------------------
954 * Event / Trigger Statistics - private
955 * ---------------------------------------------------------------------------------
956 */
957
958//________________________________________________________________________
959Bool_t AliAnalysisNetParticleHelper::FillEventStats(Int_t *aEventCuts) {
960 // -- Fill event / centrality statistics
961
962 Bool_t isRejected = kFALSE;
963
964 // -- Fill event statistics
965 for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
9be43c8e 966
cb68eb1d 967 if (aEventCuts[idx])
968 isRejected = kTRUE;
969 else
970 fHEventStat0->Fill(idx);
971 }
972
973 for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
974 if (aEventCuts[idx])
975 break;
976 fHEventStat1->Fill(idx);
977 }
978
979 // -- Fill centrality statistics of accepted events
980 if (!isRejected)
981 fHCentralityStat->Fill(fCentralityBin);
982
983 return isRejected;
984}