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