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