]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioHelper.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update in PR task: drathee
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidRatioHelper.cxx
CommitLineData
0a28d543 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: ALICE Offline. *
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// AliEbyE Analysis for Particle Ratio Fluctuation //
18// Deepika Rathee | Satyajit Jena //
19// drathee@cern.ch | sjena@cern.ch //
20// Date: Wed Jul 9 18:38:30 CEST 2014 //
21// New approch to find particle ratio to reduce memory //
22// (Test Only) //
23//=========================================================================//
24
25#include "TMath.h"
26#include "TAxis.h"
27#include "TSystem.h"
28#include "TFile.h"
29#include "TPRegexp.h"
30
31#include "AliStack.h"
32#include "AliMCEvent.h"
33#include "AliMCParticle.h"
34#include "AliESDtrackCuts.h"
35#include "AliESDInputHandler.h"
36#include "AliESDpid.h"
37#include "AliAODInputHandler.h"
38#include "AliAODEvent.h"
39#include "AliAODMCParticle.h"
40#include "AliCentrality.h"
41#include "AliTracker.h"
42
43#include "AliEbyEPidRatioHelper.h"
44
45using namespace std;
46
47
48ClassImp(AliEbyEPidRatioHelper)
49//________________________________________________________________________
50AliEbyEPidRatioHelper::AliEbyEPidRatioHelper() :
51 fModeDistCreation(0),
52
53 fInputEventHandler(NULL),
54 fPIDResponse(NULL),
55 fESD(NULL),
56 fESDTrackCuts(NULL),
57 fAOD(NULL),
58 fAODtrackCutBit(1024),
59 fIsMC(kFALSE),
60 fMCEvent(NULL),
61 fStack(NULL),
62
63 fCentralityBin(-1),
64 fCentralityPercentile(-1.),
65
bdb6c35e 66 fCentralityBinMax(11),
0a28d543 67 fVertexZMax(10.),
9b5bc5a1 68 fRapidityMax(0.9),
0a28d543 69 fPhiMin(0.),
70 fPhiMax(TMath::TwoPi()),
71 fMinTrackLengthMC(70.),
72 fNSigmaMaxCdd(3.),
73 fNSigmaMaxCzz(3.),
74
75 fPIDStrategy(0),
76 fNSigmaMaxITS(4.),
77 fNSigmaMaxTPC(4.),
78 fNSigmaMaxTPClow(3.),
79 fNSigmaMaxTOF(4.),
80 fMinPtForTOFRequired(0.69),
81 fMaxPtForTPClow(0.69),
82
83 fHEventStat0(NULL),
84 fHEventStat1(NULL),
85 fHEventStatMax(6),
86
87 fHTriggerStat(NULL),
88 fNTriggers(5),
89
90 fHCentralityStat(NULL),
15bb9247 91 fHCentralityPer(NULL),
92 fHCentralityPerAll(NULL),
bdb6c35e 93 fNCentralityBins(11),
0a28d543 94
c3bf3f97 95 fSubSamples(25),
f7ea34d2 96 fRandom(NULL),
c3bf3f97 97 fSubSampleIdx(1),
98
1648d22e 99 fIsRatio(kFALSE),
56a3fc5b 100 fIsPtBin(kFALSE), fIsDetectorWise(kFALSE) {
0a28d543 101 // Constructor
102
103 AliLog::SetClassDebugLevel("AliEbyEPidRatioHelper",10);
104}
105
f7ea34d2 106const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthRap = 0.1;
0a28d543 107const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthPt = 0.3; // 0.08 // 300 MeV // was 80 MeV
108
109const Float_t AliEbyEPidRatioHelper::fgkfHistRangeCent[] = {-0.5, 10.5};
110const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsCent = 11 ;
111
112const Float_t AliEbyEPidRatioHelper::fgkfHistRangeEta[] = {-0.9, 0.9};
113const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsEta = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeEta[1] -
114 AliEbyEPidRatioHelper::fgkfHistRangeEta[0]) /
115 AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
116
f7ea34d2 117const Float_t AliEbyEPidRatioHelper::fgkfHistRangeRap[] = {-0.8, 0.8};
0a28d543 118const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsRap = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeRap[1] - AliEbyEPidRatioHelper::fgkfHistRangeRap[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
119
120const Float_t AliEbyEPidRatioHelper::fgkfHistRangePhi[] = {0.0, TMath::TwoPi()};
f7ea34d2 121const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsPhi = 21 ;
0a28d543 122
123const Float_t AliEbyEPidRatioHelper::fgkfHistRangePt[] = {0.2, 2.9}; // {0.2, 5.}; // was {0.3, 2.22}
124const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsPt = Int_t((AliEbyEPidRatioHelper::fgkfHistRangePt[1] - AliEbyEPidRatioHelper::fgkfHistRangePt[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthPt);
125
126const Float_t AliEbyEPidRatioHelper::fgkfHistRangeSign[] = {-1.5, 1.5};
127const Int_t AliEbyEPidRatioHelper::fgkfHistNBinsSign = 3;
128
bdb6c35e 129const Char_t* AliEbyEPidRatioHelper::fgkEventNames[] = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,100]%"};
0a28d543 130const Char_t* AliEbyEPidRatioHelper::fgkCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
131const Char_t* AliEbyEPidRatioHelper::fgkTriggerNames[] = {"kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
132const Char_t* AliEbyEPidRatioHelper::fgkCentralityNames[] = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%","50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
133
134const Char_t* AliEbyEPidRatioHelper::fgkPidName[4] = {"Nch","Npi","Nka","Npr"};
135const Char_t* AliEbyEPidRatioHelper::fgkPidLatex[4][2] = {{"N_{-}","N_{+}"}, {"N_{#pi^{-}}","N_{#pi^{+}}"},{"N_{K^{-}}","N_{K^{+}}"}, {"N_{#bar{p}}","N_{p}"}};
136const Char_t* AliEbyEPidRatioHelper::fgkPidTitles[4][2] = {{"Negative","Positive"},{"Anti-Pions","Pions"},{"Anti-Kaons","Kaons"}, {"Anti-Protons","Protons"}};
137
138
139
140//________________________________________________________________________
141AliEbyEPidRatioHelper::~AliEbyEPidRatioHelper() {
142 // Destructor
143
144 if (fRandom)
145 delete fRandom;
146
147 return;
148}
149
150//________________________________________________________________________
151void AliEbyEPidRatioHelper::SetPhiRange(Float_t f1, Float_t f2) {
152 // -- Set phi range and adopt to phi-histogram
153
154 fPhiMin = f1;
155 fPhiMax = (f1 < f2) ? f2 : f2+TMath::TwoPi();
156
157 Float_t phiMin = fPhiMin;
158 Float_t phiMax = fPhiMax;
159
160 // -- Update Ranges
161 Float_t binWidth = (AliEbyEPidRatioHelper::fgkfHistRangePhi[1] - AliEbyEPidRatioHelper::fgkfHistRangePhi[0]) /
162 Float_t(AliEbyEPidRatioHelper::fgkfHistNBinsPhi);
163
164 Float_t lowEdge = AliEbyEPidRatioHelper::fgkfHistRangePhi[0] - binWidth;
165 Float_t highEdge = AliEbyEPidRatioHelper::fgkfHistRangePhi[0];
166
167 for (Int_t ii = 1; ii <= AliEbyEPidRatioHelper::fgkfHistNBinsPhi; ++ii) {
168 lowEdge += binWidth;
169 highEdge += binWidth;
170
171 if (phiMin >= lowEdge && phiMin < highEdge )
172 phiMin = lowEdge;
173 if (phiMax > lowEdge && phiMax <= highEdge )
174 phiMax = highEdge;
175 }
176
177 printf(">>>> Update Phi Range : [%f,%f] -> [%f,%f]\n", fPhiMin, fPhiMax, phiMin, phiMax);
178 fPhiMin = phiMin;
179 fPhiMax = phiMax;
180}
181
182
183
184//________________________________________________________________________
56a3fc5b 185Int_t AliEbyEPidRatioHelper::Initialize(AliESDtrackCuts *cuts, Bool_t isMC, Bool_t isRatio, Bool_t isPtBin, Bool_t isDetWise, Int_t trackCutBit, Int_t modeDistCreation) {
0a28d543 186 // -- Initialize helper
187
188 Int_t iResult = 0;
189
190 // -- ESD track cuts
191 fESDTrackCuts = cuts;
192
1648d22e 193
0a28d543 194 // -- Is MC
195 fIsMC = isMC;
f7ea34d2 196 fIsRatio = isRatio;
1648d22e 197 fIsPtBin = isPtBin;
56a3fc5b 198 fIsDetectorWise = isDetWise;
1648d22e 199
200
0a28d543 201
202 // -- AOD track filter bit
203 fAODtrackCutBit = trackCutBit;
204
205 // -- mode Distribution creation
206 fModeDistCreation = modeDistCreation;
207
208 // -- Setup event cut statistics
209 InitializeEventStats();
210
211 // -- Setup trigger statistics
212 InitializeTriggerStats();
213
214 // -- Setup centrality statistics
215 InitializeCentralityStats();
216
217 // -- PRINT PID Strategy
218 // 0 : TPC(TPClow+TPCHigh)
219 // 1 : ITS
220 // 2 : TOF
221 // 3 : ITS+TPC(TPClow+TPCHigh)
222 // 4 : TPC(TPClow+TPCHigh)+TOF
223 // 5 : TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
224 // 6 : TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
225 // 7 : TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
226 // 8 : TPC(TPClow+TPCHigh)+ITS+TOF
227 printf(">>>> PID STRATEGY: %d || sigmaMax: ITS %.2f TPC %.2f TOF %.2f \n", fPIDStrategy, fNSigmaMaxITS, fNSigmaMaxTPC, fNSigmaMaxTOF);
228
229 // -- Initialize random number generator
230 fRandom = new TRandom3();
231 fRandom->SetSeed();
232
233 return iResult;
234}
235
236//________________________________________________________________________
237Int_t AliEbyEPidRatioHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
238 // -- Setup Event
239
240 // -- Get ESD objects
241 if(esdHandler){
242 fInputEventHandler = static_cast<AliInputEventHandler*>(esdHandler);
243 fESD = dynamic_cast<AliESDEvent*>(fInputEventHandler->GetEvent());
244 if (!fESD) {
245 AliError("ESD event handler not available");
246 return -1;
247 }
248 }
249
250 // -- Get AOD objects
251 else if(aodHandler){
252 fInputEventHandler = static_cast<AliInputEventHandler*>(aodHandler);
253 fAOD = dynamic_cast<AliAODEvent*>(fInputEventHandler->GetEvent());
254 if (!fAOD) {
255 AliError("AOD event handler not available");
256 return -1;
257 }
258 }
259
260 // -- Get Common objects
261 fPIDResponse = fInputEventHandler->GetPIDResponse();
262
263 // -- Get MC objects
264 fMCEvent = mcEvent;
265 if (fMCEvent)
266 fStack = fMCEvent->Stack();
267
268 // -- Get event centrality
bdb6c35e 269 // > 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90|90-100 --> 11 bins
270 // > 0 1 2 3 4 5 6 7 8 9 10
0a28d543 271
272 AliCentrality *centrality = NULL;
273
274 if(esdHandler)
275 centrality = fESD->GetCentrality();
276 else if(aodHandler)
277 centrality = fAOD->GetHeader()->GetCentralityP();
278
279 if (!centrality) {
280 AliError("Centrality not available");
281 return -1;
282 }
283
e146cc53 284
285 // Int_t a = centrality->GetCentralityClass5("V0M");
286 // if (a < 0 || a >= 20 ) fCentralityBin = -2;
287 // else if (a <= 1) fCentralityBin = a;
288 // else fCentralityBin = 1 + centrality->GetCentralityClass10("V0M");
289
290
291
0a28d543 292 Int_t centBin = centrality->GetCentralityClass10("V0M");
bdb6c35e 293 if (centBin == 0) { fCentralityBin = centrality->GetCentralityClass5("V0M"); }
294 else if (centBin == 11 || centBin == -1.) { fCentralityBin = -1; }
295 else if (centBin > 0 && centBin < fNCentralityBins) { fCentralityBin = centBin + 1; }
296 else { fCentralityBin = -2; }
297
e146cc53 298
0a28d543 299 if (fCentralityBin >= fCentralityBinMax)
300 fCentralityBin = -2;
e146cc53 301
302
0a28d543 303
304 fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
0a28d543 305
c3bf3f97 306 fSubSampleIdx = fRandom->Integer(fSubSamples);
e146cc53 307
0a28d543 308 return 0;
309}
310//________________________________________________________________________
311Bool_t AliEbyEPidRatioHelper::IsEventTriggered() {
312 // -- Check if Event is triggered and fill Trigger Histogram
313
314 Bool_t *aTriggerFired = new Bool_t[fNTriggers];
315 for (Int_t ii = 0; ii < fNTriggers; ++ii)
316 aTriggerFired[ii] = kFALSE;
317
318 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kMB)) aTriggerFired[0] = kTRUE;
319 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kCentral)) aTriggerFired[1] = kTRUE;
320 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
321 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEJE)) aTriggerFired[3] = kTRUE;
322 if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEGA)) aTriggerFired[4] = kTRUE;
323
324 Bool_t isTriggered = kFALSE;
325
326 for (Int_t ii=0; ii<fNTriggers; ++ii) {
327 if(aTriggerFired[ii]) {
328 isTriggered = kTRUE;
329 fHTriggerStat->Fill(ii);
330 }
331 }
332
333 delete[] aTriggerFired;
334
335 return isTriggered;
336}
337
338//________________________________________________________________________
339Bool_t AliEbyEPidRatioHelper::IsEventRejected() {
340 // -- Evaluate event statistics histograms
341
342 Int_t *aEventCuts = new Int_t[fHEventStatMax];
343 // set aEventCuts[ii] to 1 in case of reject
344
345 for (Int_t ii=0;ii<fHEventStatMax; ++ii)
346 aEventCuts[ii] = 0;
347
348 Int_t iCut = 0;
349
350 // -- 0 - Before Physics Selection
351 aEventCuts[iCut] = 0;
352
353 // -- 1 - No Trigger fired
354 ++iCut;
355 if (!IsEventTriggered())
356 aEventCuts[iCut] = 1;
357
358 // -- 2 - No Vertex
359 ++iCut;
360 const AliESDVertex* vtxESD = NULL;
361 const AliAODVertex* vtxAOD = NULL;
362 if (fESD){
363 vtxESD = fESD->GetPrimaryVertexTracks();
364 if (!vtxESD)
365 aEventCuts[iCut] = 1;
366 }
367 else if (fAOD){
368 vtxAOD = fAOD->GetPrimaryVertex();
369 if (!vtxAOD)
370 aEventCuts[iCut] = 1;
371 }
372
373 // -- 3 - Vertex z outside cut window
374 ++iCut;
375 if (vtxESD){
376 if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax)
377 aEventCuts[iCut] = 1;
378 }
379 else if(vtxAOD){
380 if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax)
381 aEventCuts[iCut] = 1;
382 }
383 else
384 aEventCuts[iCut] = 1;
385
386 // -- 4 - Centrality = -1 (no centrality or not hadronic)
387 ++iCut;
388 if(fCentralityBin == -1.)
389 aEventCuts[iCut] = 1;
390
391 // -- 5 - Centrality < fCentralityMax
392 ++iCut;
393 if(fCentralityBin == -2.)
394 aEventCuts[iCut] = 1;
395
396 // -- Fill statistics / reject event
397 Bool_t isRejected = FillEventStats(aEventCuts);
398
399 // -- Cleanup
400 delete[] aEventCuts;
401
402 //cout << isRejected << endl;
403
404 return isRejected;
405}
406
407//________________________________________________________________________
408Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC) {
409 // -- Check if MC particle is accepted for basic parameters
410
411 if (!particle)
412 return kFALSE;
413
414 // -- check if charged
415 if (particle->Charge() == 0.0)
416 return kFALSE;
417
418 // -- check if physical primary - ESD
419 if (fESD) {
420 if(!fStack->IsPhysicalPrimary(idxMC))
421 return kFALSE;
422 }
423 // -- check if physical primary - AOD
424 else {
425 if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
426 return kFALSE;
427 }
428
429 return kTRUE;
430}
431
432//________________________________________________________________________
433Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC) {
434 // -- Check if MC particle is accepted for basic parameters
435
436 if (!particle)
437 return kFALSE;
438
439 // -- check if charged
440 if (particle->Charge() != 0.0)
441 return kFALSE;
442
443 // -- check if physical primary - ESD
444 if (fESD) {
445 if(!fStack->IsPhysicalPrimary(idxMC))
446 return kFALSE;
447 }
448 // -- check if physical primary - AOD
449 else {
450 if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
451 return kFALSE;
452 }
453
454 return kTRUE;
455}
456
457//________________________________________________________________________
458Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP, Int_t gCurPid) {
9b5bc5a1 459
0a28d543 460 if (gCurPid == 0) {
461 yP = particle->Eta();
462 return kTRUE;
463 }
464
9b5bc5a1 465 Double_t mP = AliPID::ParticleMass(AliPID::kPion);
0a28d543 466 if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
9b5bc5a1 467 else if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
468 else if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
0a28d543 469
470 // -- Calculate rapidities and kinematics
471 Double_t p = particle->P();
472 Double_t pz = particle->Pz();
473
474 Double_t eP = TMath::Sqrt(p*p + mP*mP);
475 yP = 0.5 * TMath::Log((eP + pz) / (eP - pz));
476
477 // -- Check Rapidity window
478 if (TMath::Abs(yP) > fRapidityMax)
479 return kFALSE;
480
481 return kTRUE;
482}
483
484//________________________________________________________________________
485Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedPhi(AliVParticle *particle) {
486 // -- Check if particle is accepted
487 // > in phi
488 // > return 0 if not accepted
489
490 if (particle->Phi() > fPhiMin && particle->Phi() <= fPhiMax)
491 return kTRUE;
492 else if (particle->Phi() < fPhiMin && (particle->Phi() + TMath::TwoPi()) <= fPhiMax)
493 return kTRUE;
494 else
495 return kFALSE;
496}
497
498//_____________________________________________________________________________
499Bool_t AliEbyEPidRatioHelper::IsParticleFindable(Int_t label) {
500 // -- Check if MC particle is findable tracks
501
502 AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
503 if(!mcParticle)
504 return kFALSE;
505
506 Int_t counter;
507 Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
508
509 return (tpcTrackLength > fMinTrackLengthMC);
510}
511//________________________________________________________________________
512Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedBasicCharged(AliVTrack* track) {
513 // -- Check if track is accepted
514 // > for basic parameters
515
516 if (!track)
517 return kFALSE;
518
519 if (track->Charge() == 0)
520 return kFALSE;
521
522 return kTRUE;
523}
524
525//________________________________________________________________________
526Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP, Int_t gCurPid) {
527 if (gCurPid == 0) {
528 yP = track->Eta();
529 return kTRUE;
530 }
531
9b5bc5a1 532 Double_t mP = AliPID::ParticleMass(AliPID::kPion);
0a28d543 533 if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
9b5bc5a1 534 else if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
535 else if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
0a28d543 536
537 // -- Calculate rapidities and kinematics
538 Double_t pvec[3];
539 track->GetPxPyPz(pvec);
540
541 Double_t p = track->P();
542 Double_t eP = TMath::Sqrt(p*p + mP*mP);
543 yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
544
545 // -- Check Rapidity window
546 if (TMath::Abs(yP) > fRapidityMax)
547 return kFALSE;
548
549 return kTRUE;
550}
551
552//________________________________________________________________________
553Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedDCA(AliVTrack *vTrack) {
554 Bool_t isAccepted = kTRUE;
555
556 if (!fESD)
557 return isAccepted;
558
559 AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
560 if (!track)
561 return kFALSE;
562
563 // -- Get nHits SPD
564 if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
565
566 // -- Get DCA nSigmas
567 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
568 track->GetImpactParameters(dca,cov);
569
570 Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
571 Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
572
573 if (fNSigmaMaxCdd != 0.) {
574 if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
575 isAccepted = kFALSE;
576 }
577
578 if (fNSigmaMaxCzz != 0.) {
579 if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
580 isAccepted = kFALSE;
581 }
582 }
583
584 return isAccepted;
585}
586
587//________________________________________________________________________
588Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid, AliPID::EParticleType gCurPid) {
589
590 Bool_t isAcceptedITS = kFALSE;
591 Bool_t isAcceptedTPC = kFALSE;
592 Bool_t isAcceptedTPClow = kFALSE;
593 Bool_t isAcceptedTOF = kFALSE;
594 Bool_t isAccepted = kFALSE;
595
596 // -- In case not PID is used
597 if (gCurPid == 0) {
598 pid[0] = 10.;
599 pid[1] = 10.;
600 pid[2] = 10.;
601 return kTRUE;
602 }
603
604 if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kITS, track, gCurPid, pid[0]) == AliPIDResponse::kDetPidOk) {
605 if (TMath::Abs(pid[0]) < fNSigmaMaxITS)
606 isAcceptedITS = kTRUE;
607 }
608
609 if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, gCurPid, pid[1]) == AliPIDResponse::kDetPidOk) {
610 if (TMath::Abs(pid[1]) < fNSigmaMaxTPC)
611 isAcceptedTPC = kTRUE;
612 if (TMath::Abs(pid[1]) < fNSigmaMaxTPClow)
613 isAcceptedTPClow = kTRUE;
614 if (track->Pt() < fMaxPtForTPClow)
615 isAcceptedTPC = isAcceptedTPClow;
616 }
617
618 Bool_t hasPIDTOF = kFALSE;
619 if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, gCurPid, pid[2]) == AliPIDResponse::kDetPidOk) {
620 hasPIDTOF = kTRUE;
621 if (TMath::Abs(pid[2]) < fNSigmaMaxTOF)
622 isAcceptedTOF = kTRUE;
623 }
0a28d543 624 // -- Check TOF missmatch for MC
625
626 //if (ESD)
627 if (fIsMC && isAcceptedTOF) {
628 Int_t tofLabel[3];
629 // AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
630 // TODO add code for AOD
631
632 (dynamic_cast<AliESDtrack*>(track))->GetTOFLabel(tofLabel);
633
634 Bool_t hasMatchTOF = kTRUE;
635 if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0)
636 hasMatchTOF = kFALSE;
637
638 TParticle *matchedTrack = fStack->Particle(TMath::Abs(tofLabel[0]));
639 if (TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
640 hasMatchTOF = kTRUE;
641
642 isAcceptedTOF = hasMatchTOF;
643 }
644
645 // 0 : TPC(TPClow+TPCHigh)
646 // 1 : ITS
647 // 2 : TOF
648 // 3 : ITS+TPC(TPClow+TPCHigh)
649 // 4 : TPC(TPClow+TPCHigh)+TOF
650 // 5 : TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
651 // 6 : TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
652 // 7 : TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
653 // 8 : TPC(TPClow+TPCHigh)+ITS+TOF
654 if (fPIDStrategy == 0) { // TPC PID
655 isAccepted = isAcceptedTPC;
656 }
657 else if (fPIDStrategy == 1) { // ITS PID
658 isAccepted = isAcceptedITS;
659 }
660 else if (fPIDStrategy == 2) { // TOF PID
661 isAccepted = isAcceptedTOF;
662 }
663 else if (fPIDStrategy == 3) { // TPC+ITS PID
664 isAccepted = isAcceptedTPC && isAcceptedITS;
665 }
666 else if (fPIDStrategy == 4) { // TPC+TOF PID
667 isAccepted = isAcceptedTPC && isAcceptedTOF;
668 }
669 else if (fPIDStrategy == 5) { // TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
670 if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
671 isAcceptedTOF = kTRUE;
672 isAccepted = isAcceptedTPC && isAcceptedTOF;
673 }
674 else if (fPIDStrategy == 6) { // ITS+TPC+TOF PID -- TOF only for those tracks which have TOF information
675 isAccepted = isAcceptedTPC && isAcceptedITS;
676 if (hasPIDTOF)
677 isAccepted = isAccepted && isAcceptedTOF;
678 }
679 else if (fPIDStrategy == 7) { // ITS+TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
680 if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
681 isAcceptedTOF = kTRUE;
682 isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
683 }
684 else if (fPIDStrategy == 8) { // ITS+TPC+TOF PID
685 isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
686 }
687
688 return isAccepted;
689}
690
691//________________________________________________________________________
692Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedPhi(AliVTrack *track) {
693 // -- Check if track is accepted
694 // > in phi
695 // > return 0 if not accepted
696
697 if (track->Phi() > fPhiMin && track->Phi() <= fPhiMax)
698 return kTRUE;
699 else if (track->Phi() < fPhiMin && (track->Phi() + TMath::TwoPi()) <= fPhiMax)
700 return kTRUE;
701 else
702 return kFALSE;
703}
704
705//________________________________________________________________________
706void AliEbyEPidRatioHelper::BinLogAxis(const THnBase *hn, Int_t axisNumber, AliESDtrackCuts* cuts) {
707 AliESDtrackCuts* esdTrackCuts = (cuts) ? cuts : fESDTrackCuts;
708
709 // -- Make logarithmic binning
710 TAxis *axis = hn->GetAxis(axisNumber);
711 Int_t nBins = axis->GetNbins();
712
713 Double_t from = axis->GetXmin();
714 Double_t to = axis->GetXmax();
715 Double_t *newBins = new Double_t[nBins + 1];
716
717 newBins[0] = from;
718 Double_t factor = TMath::Power(to/from, 1./nBins);
719
720 for (int ii = 1; ii <= nBins; ii++)
721 newBins[ii] = factor * newBins[ii-1];
722
723 axis->Set(nBins, newBins);
724
725 delete [] newBins;
726
727 // -- Update Ranges
728 // ------------------
729 Float_t ptRange[2];
730 Float_t oldPtRange[2];
731 esdTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
732 esdTrackCuts->GetPtRange(oldPtRange[0],oldPtRange[1]);
733
734 Float_t minPtForTOFRequired = fMinPtForTOFRequired;
735 Float_t maxPtForTPClow = fMaxPtForTPClow;
736
737 // -- Update minPt Cut
6ce4ad92 738 Int_t bin = axis->FindBin(ptRange[0]+10e-7);
0a28d543 739 ptRange[0] = axis->GetBinLowEdge(bin);
740
741 // -- Update maxPt Cut
6ce4ad92 742 bin = axis->FindBin(ptRange[1]-10e-7);
0a28d543 743 ptRange[1] = axis->GetBinUpEdge(bin);
744
745 if (ptRange[0] != oldPtRange[0] || ptRange[1] != oldPtRange[1]) {
746 printf(">>>> Update Pt Range : [%f,%f] -> [%f,%f]\n", oldPtRange[0], oldPtRange[1], ptRange[0], ptRange[1]);
747 esdTrackCuts->SetPtRange(ptRange[0],ptRange[1]);
748 }
749
750 // -- Update MinPtForTOFRequired
6ce4ad92 751 bin = axis->FindBin(minPtForTOFRequired-10e-7);
0a28d543 752 minPtForTOFRequired = axis->GetBinUpEdge(bin);
753
754 if (minPtForTOFRequired != fMinPtForTOFRequired) {
755 printf(">>>> Update Min Pt for TOF : %f -> %f\n", fMinPtForTOFRequired, minPtForTOFRequired);
756 fMinPtForTOFRequired = minPtForTOFRequired;
757 }
758
759 // -- Update MaxPtForTPClow
6ce4ad92 760 bin = axis->FindBin(maxPtForTPClow-10e-7);
0a28d543 761 maxPtForTPClow = axis->GetBinUpEdge(bin);
762
763 if (maxPtForTPClow != fMaxPtForTPClow) {
764 printf(">>>> Update Max Pt for TPC Low : %f -> %f\n", fMaxPtForTPClow, maxPtForTPClow);
765 fMaxPtForTPClow = maxPtForTPClow;
766 }
767}
768
769//________________________________________________________________________
770void AliEbyEPidRatioHelper::InitializeEventStats() {
771 // -- Initialize event statistics histograms
772
773 fHEventStat0 = new TH1F("hEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
774 fHEventStat1 = new TH1F("hEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
775
776 for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
777 fHEventStat0->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkEventNames[ii]);
778 fHEventStat1->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkEventNames[ii]);
779 }
780
781 fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
782 fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
783}
784
785//________________________________________________________________________
786void AliEbyEPidRatioHelper::InitializeTriggerStats() {
787 // -- Initialize trigger statistics histograms
788
789 fHTriggerStat = new TH1F("hTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
790
791 for ( Int_t ii=0; ii < fNTriggers; ii++ )
792 fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkTriggerNames[ii]);
793}
794
795//________________________________________________________________________
796void AliEbyEPidRatioHelper::InitializeCentralityStats() {
797 // -- Initialize trigger statistics histograms
798
799 fHCentralityStat = new TH1F("hCentralityStat","Centrality statistics;Centrality Bins;Events",
800 fNCentralityBins,-0.5,fNCentralityBins-0.5);
15bb9247 801
0a28d543 802 for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
803 fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkCentralityNames[ii]);
15bb9247 804
805 fHCentralityPer = new TH1F("hCentralityPercentileAccepted","Centrality Percentile statistics;Centrality Bins;Events",
806 100,-0.5,99.5);
807
808 fHCentralityPerAll = new TH1F("hCentralityPercentileAll","Centrality Percentile statistics;Centrality Bins;Events",
809 100,-0.5,99.5);
810
811
0a28d543 812}
813//________________________________________________________________________
814Bool_t AliEbyEPidRatioHelper::FillEventStats(Int_t *aEventCuts) {
815 // -- Fill event / centrality statistics
816
817 Bool_t isRejected = kFALSE;
818
e146cc53 819
0a28d543 820 // -- Fill event statistics
821 for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
822
823 if (aEventCuts[idx])
824 isRejected = kTRUE;
825 else
826 fHEventStat0->Fill(idx);
827 }
828
829 for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
830 if (aEventCuts[idx])
831 break;
832 fHEventStat1->Fill(idx);
833 }
834
835 // -- Fill centrality statistics of accepted events
15bb9247 836 if (!isRejected) {
0a28d543 837 fHCentralityStat->Fill(fCentralityBin);
15bb9247 838 fHCentralityPer->Fill(fCentralityPercentile);
839 }
840
841 fHCentralityPerAll->Fill(fCentralityPercentile);
842
843
0a28d543 844
845 return isRejected;
846}