]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/NetParticle/AliAnalysisNetParticleHelper.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes for DQM QA histos
[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"
9be43c8e 8#include "TDatabasePDG.h"
cb68eb1d 9
10#include "AliStack.h"
11#include "AliMCEvent.h"
12#include "AliMCParticle.h"
13#include "AliESDtrackCuts.h"
14#include "AliESDInputHandler.h"
15#include "AliESDpid.h"
9be43c8e 16#include "AliAODInputHandler.h"
17#include "AliAODEvent.h"
18#include "AliAODMCParticle.h"
cb68eb1d 19#include "AliCentrality.h"
20#include "AliTracker.h"
21
22#include "AliAnalysisNetParticleHelper.h"
23
24using namespace std;
25
26// Task for NetParticle checks
27// Author: Jochen Thaeder <jochen@thaeder.de>
28
29ClassImp(AliAnalysisNetParticleHelper)
30
3aa68b92 31/*
32 * ---------------------------------------------------------------------------------
33 * particle names
34 * ---------------------------------------------------------------------------------
35 */
36
37 const Char_t* aPartNames[AliPID::kSPECIES][2] = {
38 {"ele", "posi"},
39 {"mubar", "mu"},
40 {"pibar", "pi"},
41 {"kbar", "k"},
42 {"pbar", "p"}
43 };
44
45 const Char_t* aPartTitles[AliPID::kSPECIES][2] = {
46 {"Electron", "Positron"},
47 {"Anti-Muon", "Muon"},
48 {"Anti-Pion", "Proton"},
49 {"Anti-Kaon", "Kaon"},
50 {"Anti-Proton", "Proton"}
51 };
52
53 const Char_t* aPartTitlesLatex[AliPID::kSPECIES][2] = {
54 {"e^{-}", "e^{+}" },
55 {"#mu^{-}", "#mu^{+}" },
56 {"#pi^{-}", "#pi^{+}" },
57 {"K^{-}", "K^{+}" },
58 {"#bar{p}", "p"}
59 };
60
61 const Char_t* aControlPartNames[2] = {"lambdabar", "lambda"};
62 const Char_t* aControlPartTitles[2] = {"Anti-Lambda", "Lambda"};
63
cb68eb1d 64/*
65 * ---------------------------------------------------------------------------------
66 * Constructor / Destructor
67 * ---------------------------------------------------------------------------------
68 */
69
70//________________________________________________________________________
71AliAnalysisNetParticleHelper::AliAnalysisNetParticleHelper() :
72 fESDHandler(NULL),
73 fPIDResponse(NULL),
74 fESD(NULL),
9be43c8e 75 fAODHandler(NULL),
76 fAOD(NULL),
cb68eb1d 77 fMCEvent(NULL),
78 fStack(NULL),
79
80 fCentralityBin(-1),
81 fCentralityPercentile(-1.),
82
83 fCentralityBinMax(7),
84 fVertexZMax(10.),
85 fRapidityMax(0.5),
86 fMinTrackLengthMC(70.),
87 fNSigmaMaxCdd(3.),
88 fNSigmaMaxCzz(3.),
89
90 fParticleSpecies(AliPID::kProton),
91 fControlParticleCode(3122),
92 fControlParticleIsNeutral(kTRUE),
93 fControlParticleName("Lambda"),
94
95 fNSigmaMaxTPC(2.5),
96 fNSigmaMaxTOF(2.5),
97 fMinPtForTOFRequired(0.8),
98
99 fHEventStat0(NULL),
100 fHEventStat1(NULL),
101 fHEventStatMax(6),
102
103 fHTriggerStat(NULL),
104 fNTriggers(5),
105
106 fHCentralityStat(NULL),
107 fNCentralityBins(10),
108
109 fEtaCorrFunc(NULL),
110 fCorr0(NULL),
111 fCorr1(NULL) {
112 // Constructor
113
114 AliLog::SetClassDebugLevel("AliAnalysisNetParticleHelper",10);
115}
116
117//________________________________________________________________________
118AliAnalysisNetParticleHelper::~AliAnalysisNetParticleHelper() {
119 // Destructor
120
121 for (Int_t jj = 0; jj < 2; ++jj) {
122 if (fCorr0[jj]) delete[] fCorr0[jj];
123 if (fCorr1[jj]) delete[] fCorr1[jj];
124 }
125 if (fCorr0) delete[] fCorr0;
126 if (fCorr1) delete[] fCorr1;
127
128 return;
129}
130
3aa68b92 131/*
132 * ---------------------------------------------------------------------------------
133 * Setter
134 * ---------------------------------------------------------------------------------
135 */
136
137//________________________________________________________________________
138void AliAnalysisNetParticleHelper::SetParticleSpecies(AliPID::EParticleType pid) {
139 // -- Set particle species (ID, Name, Title, Title LATEX)
140
141 if( (Int_t)pid < 0 || (Int_t)pid > AliPID::kSPECIES){
142 AliWarning("Particle ID not in AliPID::kSPECIES --> Set to protons");
143 pid = AliPID::kProton;
144 }
145
146 fParticleSpecies = pid;
147
148 for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
149 fPartName[idxPart] = aPartNames[fParticleSpecies][idxPart];
150 fPartTitle[idxPart] = aPartTitles[fParticleSpecies][idxPart];
151 fPartTitleLatex[idxPart] = aPartTitlesLatex[fParticleSpecies][idxPart];
152 }
153}
154
155
156/*
157 * ---------------------------------------------------------------------------------
158 * Getter
159 * ---------------------------------------------------------------------------------
160 */
161
162//________________________________________________________________________
163TString AliAnalysisNetParticleHelper::GetParticleName(Int_t idxPart) {
164 // -- Get particle Name
165
166 if( idxPart != 0 && idxPart != 1){
167 AliWarning("Particle type not known --> Set to antiparticles");
168 idxPart = 0;
169 }
170
171 return fPartName[idxPart];
172}
173
174//________________________________________________________________________
175TString AliAnalysisNetParticleHelper::GetParticleTitle(Int_t idxPart) {
176 // -- Get particle Title
177
178 if( idxPart != 0 && idxPart != 1){
179 AliWarning("Particle type not known --> Set to antiparticles");
180 idxPart = 0;
181 }
182
183 return fPartTitle[idxPart];
184}
185
186//________________________________________________________________________
187TString AliAnalysisNetParticleHelper::GetParticleTitleLatex(Int_t idxPart) {
188 // -- Get particle Title LATEX
189
190 if( idxPart != 0 && idxPart != 1){
191 AliWarning("Particle type not known --> Set to antiparticles");
192 idxPart = 0;
193 }
194
195 return fPartTitleLatex[idxPart];
196}
197
198//________________________________________________________________________
199TString AliAnalysisNetParticleHelper::GetControlParticleName(Int_t idxPart) {
200 // -- Get particle Title LATEX
201
202 if( idxPart != 0 && idxPart != 1){
203 AliWarning("Particle type not known --> Set to antiparticles");
204 idxPart = 0;
205 }
206
207 return aControlPartNames[idxPart];
208}
209
210//________________________________________________________________________
211TString AliAnalysisNetParticleHelper::GetControlParticleTitle(Int_t idxPart) {
212 // -- Get particle Title LATEX
213
214 if( idxPart != 0 && idxPart != 1){
215 AliWarning("Particle type not known --> Set to antiparticles");
216 idxPart = 0;
217 }
218
219 return aControlPartTitles[idxPart];
220}
221
cb68eb1d 222/*
223 * ---------------------------------------------------------------------------------
224 * Public Methods
225 * ---------------------------------------------------------------------------------
226 */
227
228//________________________________________________________________________
229Int_t AliAnalysisNetParticleHelper::Initialize(Bool_t isMC) {
230 // -- Initialize helper
231
232 Int_t iResult = 0;
233
234 // -- Setup event cut statistics
235 InitializeEventStats();
236
237 // -- Setup trigger statistics
238 InitializeTriggerStats();
239
240 // -- Setup centrality statistics
241 InitializeCentralityStats();
242
243 // -- Load Eta correction function
244 iResult = InitializeEtaCorrection(isMC);
245
3aa68b92 246 // -- Load track by track correction function
cb68eb1d 247 iResult = InitializeTrackbyTrackCorrection();
248
249 return iResult;
250}
251
252//________________________________________________________________________
9be43c8e 253Int_t AliAnalysisNetParticleHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
cb68eb1d 254 // -- Setup Event
255
256 // -- Get ESD objects
9be43c8e 257 if(esdHandler){
258 fESDHandler = esdHandler;
259 fPIDResponse = esdHandler->GetPIDResponse();
260 fESD = fESDHandler->GetEvent();
261 }
262
263 // -- Get AOD objects
264 else if(aodHandler){
265 fAODHandler = aodHandler;
266 fPIDResponse = aodHandler->GetPIDResponse();
267 fAOD = fAODHandler->GetEvent();
268 }
cb68eb1d 269
270 // -- Get MC objects
271 fMCEvent = mcEvent;
272 if (fMCEvent)
273 fStack = fMCEvent->Stack();
274
275 // -- Get event centrality
276 // > 0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
277 // > 0 1 2 3 4 5 6 7 8 9
278
9be43c8e 279 AliCentrality *centrality = NULL;
280
281 if(esdHandler)
282 centrality = fESD->GetCentrality();
283 else if(aodHandler)
284 centrality = fAOD->GetHeader()->GetCentralityP();
285
cb68eb1d 286 Int_t centBin = centrality->GetCentralityClass10("V0M");
287 if (centBin == 0)
288 fCentralityBin = centrality->GetCentralityClass5("V0M");
289 else if (centBin == 10 || centBin == -1.)
290 fCentralityBin = -1;
291 else if (centBin > 0 && centBin < fNCentralityBins)
292 fCentralityBin = centBin + 1;
293 else
294 fCentralityBin = -2;
9be43c8e 295
cb68eb1d 296 // -- Stay within the max centrality bin
297 if (fCentralityBin >= fCentralityBinMax)
298 fCentralityBin = -2;
9be43c8e 299
cb68eb1d 300 fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
301
9be43c8e 302 if(esdHandler){
303 // -- Update TPC pid with eta correction (only for ESDs?) XXX
304 UpdateEtaCorrectedTPCPid();
305 }
cb68eb1d 306
307 return 0;
308}
309
310/*
311 * ---------------------------------------------------------------------------------
312 * Event / Trigger Statistics
313 * ---------------------------------------------------------------------------------
314 */
315
316//________________________________________________________________________
317Bool_t AliAnalysisNetParticleHelper::IsEventTriggered() {
318 // -- Check if Event is triggered and fill Trigger Histogram
319
320 Bool_t *aTriggerFired = new Bool_t[fNTriggers];
321 for (Int_t ii = 0; ii < fNTriggers; ++ii)
322 aTriggerFired[ii] = kFALSE;
323
9be43c8e 324 if(fESDHandler){
325 if ((fESDHandler->IsEventSelected() & AliVEvent::kMB)) aTriggerFired[0] = kTRUE;
326 if ((fESDHandler->IsEventSelected() & AliVEvent::kCentral)) aTriggerFired[1] = kTRUE;
327 if ((fESDHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
328 if ((fESDHandler->IsEventSelected() & AliVEvent::kEMCEJE)) aTriggerFired[3] = kTRUE;
329 if ((fESDHandler->IsEventSelected() & AliVEvent::kEMCEGA)) aTriggerFired[4] = kTRUE;
330 }
331 else if(fAODHandler){
332 if ((fAODHandler->IsEventSelected() & AliVEvent::kMB)) aTriggerFired[0] = kTRUE;
333 if ((fAODHandler->IsEventSelected() & AliVEvent::kCentral)) aTriggerFired[1] = kTRUE;
334 if ((fAODHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
335 if ((fAODHandler->IsEventSelected() & AliVEvent::kEMCEJE)) aTriggerFired[3] = kTRUE;
336 if ((fAODHandler->IsEventSelected() & AliVEvent::kEMCEGA)) aTriggerFired[4] = kTRUE;
337 }
338
cb68eb1d 339 Bool_t isTriggered = kFALSE;
340
341 for (Int_t ii=0; ii<fNTriggers; ++ii) {
342 if(aTriggerFired[ii]) {
343 isTriggered = kTRUE;
344 fHTriggerStat->Fill(ii);
345 }
346 }
347
348 delete[] aTriggerFired;
349
350 return isTriggered;
351}
352
353//________________________________________________________________________
354Bool_t AliAnalysisNetParticleHelper::IsEventRejected() {
355 // -- Evaluate event statistics histograms
356
357 Int_t *aEventCuts = new Int_t[fHEventStatMax];
358 // set aEventCuts[ii] to 1 in case of reject
359
360 for (Int_t ii=0;ii<fHEventStatMax; ++ii)
361 aEventCuts[ii] = 0;
362
363 Int_t iCut = 0;
364
365 // -- 0 - Before Physics Selection
366 aEventCuts[iCut] = 0;
367
368 // -- 1 - No Trigger fired
369 ++iCut;
370 if (!IsEventTriggered())
371 aEventCuts[iCut] = 1;
372
373 // -- 2 - No Vertex
374 ++iCut;
9be43c8e 375 const AliESDVertex* vtxESD = NULL;
376 const AliAODVertex* vtxAOD = NULL;
377 if (fESD){
378 vtxESD = fESD->GetPrimaryVertexTracks();
379 if (!vtxESD)
380 aEventCuts[iCut] = 1;
381 }
382 else if (fAOD){
383 vtxAOD = fAOD->GetPrimaryVertex();
384 if (!vtxAOD)
385 aEventCuts[iCut] = 1;
386 }
cb68eb1d 387
388 // -- 3 - Vertex z outside cut window
389 ++iCut;
390 if (vtxESD){
391 if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax)
392 aEventCuts[iCut] = 1;
393 }
9be43c8e 394 else if(vtxAOD){
395 if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax)
396 aEventCuts[iCut] = 1;
397 }
cb68eb1d 398 else
399 aEventCuts[iCut] = 1;
9be43c8e 400
cb68eb1d 401 // -- 4 - Centrality = -1 (no centrality or not hadronic)
402 ++iCut;
403 if(fCentralityBin == -1.)
404 aEventCuts[iCut] = 1;
405
406 // -- 5 - Centrality < fCentralityMax
407 ++iCut;
408 if(fCentralityBin == -2.)
409 aEventCuts[iCut] = 1;
410
411 // -- Fill statistics / reject event
412 Bool_t isRejected = FillEventStats(aEventCuts);
413
414 // -- Cleanup
415 delete[] aEventCuts;
416
417 return isRejected;
418}
419
420/*
421 * ---------------------------------------------------------------------------------
422 * Accept Particle Methods - private
423 * ---------------------------------------------------------------------------------
424 */
425
426//________________________________________________________________________
427Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(TParticle *particle, Int_t idxMC) {
428 // -- Check if MC particle is accepted for basic parameters
429
430 if (!particle)
431 return kFALSE;
432
433 // -- check if PDF code exists
434 if (!particle->GetPDG())
435 return kFALSE;
436
437 // -- check if charged
438 if (particle->GetPDG()->Charge() == 0.0)
439 return kFALSE;
440
441 // -- check if physical primary
442 if(!fStack->IsPhysicalPrimary(idxMC))
443 return kFALSE;
444
445 return kTRUE;
446}
447//________________________________________________________________________
9be43c8e 448Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(AliAODMCParticle *particle) {
449 // -- Check if MC particle is accepted for basic parameters
450
451 if (!particle)
452 return kFALSE;
453
454 // -- check if PDF code exists
455 if (!(TDatabasePDG::Instance()->GetParticle(particle->PdgCode())))
456 return kFALSE;
457
458 // -- check if charged
459 if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() == 0.0)
460 return kFALSE;
461
462 // -- check if physical primary
463 if(!particle->IsPhysicalPrimary())
464 return kFALSE;
465
466 return kTRUE;
467}
468//________________________________________________________________________
cb68eb1d 469Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(TParticle *particle, Int_t idxMC) {
470 // -- Check if MC particle is accepted for basic parameters
471
472 if (!particle)
473 return kFALSE;
474
475 // -- check if PDF code exists
476 if (!particle->GetPDG())
477 return kFALSE;
478
479 // -- check if neutral
480 if (particle->GetPDG()->Charge() != 0.0)
481 return kFALSE;
482
483 // -- check if physical primary
484 if(!fStack->IsPhysicalPrimary(idxMC))
485 return kFALSE;
486
487 return kTRUE;
488}
9be43c8e 489//________________________________________________________________________
490Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(AliAODMCParticle *particle) {
491 // -- Check if MC particle is accepted for basic parameters
492
493 if (!particle)
494 return kFALSE;
cb68eb1d 495
9be43c8e 496 // -- check if PDF code exists
497 if (!(TDatabasePDG::Instance()->GetParticle(particle->PdgCode())))
498 return kFALSE;
499
500 // -- check if charged
501 if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() != 0.0)
502 return kFALSE;
503
504 // -- check if physical primary
505 if(!particle->IsPhysicalPrimary())
506 return kFALSE;
507
508 return kTRUE;
509}
cb68eb1d 510//________________________________________________________________________
511Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedRapidity(TParticle *particle, Double_t &yP) {
512 // -- Check if particle is accepted
513 // > in rapidity
514 // > return 0 if not accepted
515
516 Double_t mP = AliPID::ParticleMass(fParticleSpecies);
517
518 // -- Calculate rapidities and kinematics
519 Double_t p = particle->P();
520 Double_t pz = particle->Pz();
521
522 Double_t eP = TMath::Sqrt(p*p + mP*mP);
523 yP = 0.5 * TMath::Log((eP + pz) / (eP - pz));
524
525 // -- Check Rapidity window
526 if (TMath::Abs(yP) > fRapidityMax)
527 return kFALSE;
528
529 return kTRUE;
530}
531
532//_____________________________________________________________________________
533Bool_t AliAnalysisNetParticleHelper::IsParticleFindable(Int_t label) {
534 // -- Check if MC particle is findable tracks
535
536 AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
537 if(!mcParticle)
538 return kFALSE;
539
540 Int_t counter;
541 Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
542
543 return (tpcTrackLength > fMinTrackLengthMC);
544}
545
546/*
547 * ---------------------------------------------------------------------------------
548 * Accept Track Methods - public
549 * ---------------------------------------------------------------------------------
550 */
9be43c8e 551
cb68eb1d 552//________________________________________________________________________
553Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliESDtrack *track) {
554 // -- Check if track is accepted
555 // > for basic parameters
9be43c8e 556
cb68eb1d 557 if (!track)
558 return kFALSE;
559
560 if (track->Charge() == 0)
561 return kFALSE;
562
563 // -- Get momentum for dEdx
564 if (!track->GetInnerParam())
565 return kFALSE;
566
9be43c8e 567 return kTRUE;
568}
569
570//________________________________________________________________________
571Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliAODTrack *track) {
572 // -- Check if track is accepted
573 // > for basic parameters
574
575 if (!track)
576 return kFALSE;
577
578 if (track->Charge() == 0)
579 return kFALSE;
580
581 //// -- Get momentum for dEdx --> returns always ZERO XXX
582 //if (!track->GetInnerParam())
583 // return kFALSE;
584
cb68eb1d 585 return kTRUE;
586}
587
588//________________________________________________________________________
9be43c8e 589Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP) {
cb68eb1d 590 // -- Check if track is accepted
591 // > in rapidity
592 // > return 0 if not accepted
593
594 Double_t mP = AliPID::ParticleMass(fParticleSpecies);
595
596 // -- Calculate rapidities and kinematics
597 Double_t pvec[3];
598 track->GetPxPyPz(pvec);
599
9be43c8e 600 Double_t p = track->P();
cb68eb1d 601 Double_t eP = TMath::Sqrt(p*p + mP*mP);
602 yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
603
604 // -- Check Rapidity window
605 if (TMath::Abs(yP) > fRapidityMax)
606 return kFALSE;
607
608 return kTRUE;
609}
610
611//________________________________________________________________________
9be43c8e 612Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedDCA(AliESDtrack *track) { //ONLY FOR ESDs so far XXX
cb68eb1d 613 // -- Check if track is accepted
614 // > for DCA, if both SPD layers have hits
615
616 Bool_t isAccepted = kTRUE;
617
618 // -- Get nHits SPD
619 if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
620
621 // -- Get DCA nSigmas
622 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
623 track->GetImpactParameters(dca,cov);
624
625 Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
626 Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
627
628 if (fNSigmaMaxCdd != 0.) {
629 if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
630 isAccepted = kFALSE;
631 }
632
633 if (fNSigmaMaxCzz != 0.) {
634 if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
635 isAccepted = kFALSE;
636 }
637 }
638
639 return isAccepted;
640}
641
642//________________________________________________________________________
9be43c8e 643Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid) {
cb68eb1d 644 // -- Check if track is accepted
645 // > provides TPC and TOF nSigmas to the argument
646
647 Bool_t isAcceptedTPC = kFALSE;
648 Bool_t isAcceptedTOF = kTRUE;
649
650 // -- Get PID
651 pid[0] = fPIDResponse->NumberOfSigmasTPC(track, fParticleSpecies);
652 pid[1] = fPIDResponse->NumberOfSigmasTOF(track, fParticleSpecies);
653
654 // -- Check PID with TPC
655 if (TMath::Abs(pid[0]) < fNSigmaMaxTPC)
656 isAcceptedTPC = kTRUE;
657
658 // -- Check PID with TOF
659 if (isAcceptedTPC) {
660
661 // Has TOF PID availible
662 if (track->GetStatus() & AliVTrack::kTOFpid) {
663 if (TMath::Abs(pid[1]) < fNSigmaMaxTOF)
664 isAcceptedTOF = kTRUE;
665 else
666 isAcceptedTOF = kFALSE;
667 }
668
669 // Has no TOF PID availible
670 else {
671 if (track->Pt() > fMinPtForTOFRequired)
672 isAcceptedTOF = kFALSE;
673 else
674 isAcceptedTOF = kTRUE;
675 }
676 } // if (isAcceptedTPC) {
677
678 return (isAcceptedTPC && isAcceptedTOF);
679}
680
681/*
682 * ---------------------------------------------------------------------------------
683 * Helper Methods
684 * ---------------------------------------------------------------------------------
685 */
686
687//________________________________________________________________________
688void AliAnalysisNetParticleHelper::UpdateEtaCorrectedTPCPid() {
689 // -- Update eta corrected TPC pid
690
691 if (!fEtaCorrFunc)
692 return;
693
694 for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
695 AliESDtrack *track = fESD->GetTrack(idxTrack);
696
697 // -- Check if charged track is accepted for basic parameters
698 if (!IsTrackAcceptedBasicCharged(track))
699 continue;
700
701 Double_t newTPCSignal = track->GetTPCsignal() / fEtaCorrFunc->Eval(track->Eta());
702 track->SetTPCsignal(newTPCSignal, track->GetTPCsignalSigma(), track->GetTPCsignalN());
703 }
704}
705
706//________________________________________________________________________
707Double_t AliAnalysisNetParticleHelper::GetTrackbyTrackCorrectionFactor(Double_t *aTrack, Int_t flag) {
708 // -- Get efficiency correctionf of particle dependent on (eta, phi, pt, centrality)
709
710 Int_t idxPart = (aTrack[2] < 0) ? 0 : 1;
711 THnSparseF* corrMatrix = (flag == 0) ? fCorr0[idxPart][fCentralityBin] : fCorr1[idxPart][fCentralityBin];
712
713 Double_t dimBin[3] = {aTrack[3], aTrack[4], aTrack[1]}; // eta, phi, pt
714
715 Double_t corr = corrMatrix->GetBinContent(corrMatrix->GetBin(dimBin));
716 if (corr == 0.) {
717 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",
718 aTrack[3], aTrack[4], aTrack[1], fCentralityBin));
719 corr = 1.;
720 }
721
722 return corr;
723}
724
725//________________________________________________________________________
726void AliAnalysisNetParticleHelper::BinLogAxis(const THnSparseF *h, Int_t axisNumber) {
727 // -- Method for the correct logarithmic binning of histograms
728
729 TAxis *axis = h->GetAxis(axisNumber);
730 Int_t bins = axis->GetNbins();
731
732 Double_t from = axis->GetXmin();
733 Double_t to = axis->GetXmax();
734 Double_t *newBins = new Double_t[bins + 1];
735
736 newBins[0] = from;
737 Double_t factor = pow(to/from, 1./bins);
738
739 for (int i = 1; i <= bins; i++) {
740 newBins[i] = factor * newBins[i-1];
741 }
742 axis->Set(bins, newBins);
743 delete [] newBins;
744}
745
746///////////////////////////////////////////////////////////////////////////////////
747
748/*
749 * ---------------------------------------------------------------------------------
750 * Initialize - Private
751 * ---------------------------------------------------------------------------------
752 */
753
754//________________________________________________________________________
755void AliAnalysisNetParticleHelper::InitializeEventStats() {
756 // -- Initialize event statistics histograms
757
758 const Char_t* aEventNames[] = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,90]%"};
759 const Char_t* aCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
760
761 fHEventStat0 = new TH1F("fHEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
762 fHEventStat1 = new TH1F("fHEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
763
764 for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
765 fHEventStat0->GetXaxis()->SetBinLabel(ii+1, aEventNames[ii]);
766 fHEventStat1->GetXaxis()->SetBinLabel(ii+1, aEventNames[ii]);
767 }
768 fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", aCentralityMaxNames[fCentralityBinMax-1]));
769 fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", aCentralityMaxNames[fCentralityBinMax-1]));
770}
771
772//________________________________________________________________________
773void AliAnalysisNetParticleHelper::InitializeTriggerStats() {
774 // -- Initialize trigger statistics histograms
775
776 const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
777
778 fHTriggerStat = new TH1F("fHTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
779
780 for ( Int_t ii=0; ii < fNTriggers; ii++ )
781 fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
782}
783
784//________________________________________________________________________
785void AliAnalysisNetParticleHelper::InitializeCentralityStats() {
786 // -- Initialize trigger statistics histograms
787
788 const Char_t* aCentralityNames[] = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%",
789 "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
790
791 fHCentralityStat = new TH1F("fHCentralityStat","Centrality statistics;Centrality Bins;Events",
792 fNCentralityBins,-0.5,fNCentralityBins-0.5);
793
794 for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
795 fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, aCentralityNames[ii]);
796}
797
798//________________________________________________________________________
799Int_t AliAnalysisNetParticleHelper::InitializeEtaCorrection(Bool_t isMC) {
800 // -- Initialize eta correction maps for TPC pid
801
802 TFile fileEtaCorrMaps("${ALICE_ROOT}/PWGDQ/dielectron/files/EtaCorrMaps.root");
803 if (!fileEtaCorrMaps.IsOpen())
804 return -1;
805
806 TList *keyList = fileEtaCorrMaps.GetListOfKeys();
807
808 TString sList("");
809 sList = (isMC) ? "LHC11a10" : "LHC10h.pass2";
810
811 for (Int_t idx = 0; idx < keyList->GetEntries(); ++idx) {
812 TString keyName = keyList->At(idx)->GetName();
813 TPRegexp reg(keyName);
814 if (reg.MatchB(sList)){
815 AliInfo(Form("Using Eta Correction Function: %s",keyName.Data()));
816 fEtaCorrFunc = static_cast<TF1*>(fileEtaCorrMaps.Get(keyName.Data()));
817 return 0;
818 }
819 }
820
821 return -2;
822}
823
824//________________________________________________________________________
825Int_t AliAnalysisNetParticleHelper::InitializeTrackbyTrackCorrection() {
826 // -- Initialize track by track correction matrices
827
3aa68b92 828 if( fParticleSpecies != AliPID::kProton){
829 AliWarning("Only proton corrections (preliminary available) --> Will correct with these");
830 }
831
832 // In future --> !!!!!!!!!!!!!!!!!!!!!!!
833 //TFile* corrFile = TFile::Open("${ALICE_ROOT}/PWGCF/EBYE/NetParticle/eff/effectiveCorrection.root");
cb68eb1d 834 TFile* corrFile = TFile::Open("${ALICE_ROOT}/PWGCF/EBYE/NetParticle/eff/effectiveCorrectionProton.root");
835 if (!corrFile) {
836 AliError("Track-by-track correction file can not be opened!");
837 return -1;
838 }
839
840 // -- correction - not cross section corrected
841 fCorr0 = new THnSparseF**[2];
842 fCorr0[0] = new THnSparseF*[fCentralityBinMax];
843 fCorr0[1] = new THnSparseF*[fCentralityBinMax];
844
845 for (Int_t idxCent = 0; idxCent < fCentralityBinMax; ++idxCent) {
3aa68b92 846
847 // In future --> !!!!!!!!!!!!!!!!!!!!!!!
848 // THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr0_Cent_%d", fPartName[0].Data(), idxCent)));
849 // THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr0_Cent_%d", fPartName[1].Data(), idxCent)));
cb68eb1d 850 THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("pbar_Corr0_Cent_%d", idxCent)));
851 THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("p_Corr0_Cent_%d", idxCent)));
852
853 if (!sp0 || !sp1) {
854 AliError(Form("Effective correction objects 0 - idxCent %d can not be retrieved!",idxCent));
855 return -1;
856 }
857
858 fCorr0[0][idxCent] = static_cast<THnSparseF*>(sp0->Clone());
859 fCorr0[1][idxCent] = static_cast<THnSparseF*>(sp1->Clone());
860 }
861
862 // -- correction - ross section corrected
863 fCorr1 = new THnSparseF**[2];
864 fCorr1[0] = new THnSparseF*[fCentralityBinMax];
865 fCorr1[1] = new THnSparseF*[fCentralityBinMax];
866
867 for (Int_t idxCent = 0; idxCent < fCentralityBinMax; ++idxCent) {
3aa68b92 868
869 // In future --> !!!!!!!!!!!!!!!!!!!!!!!
870 //THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr1_Cent_%d", fPartName[0].Data(), idxCent)));
871 //THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr1_Cent_%d", fPartName[1].Data(), idxCent)));
cb68eb1d 872 THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("pbar_Corr1_Cent_%d", idxCent)));
873 THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("p_Corr1_Cent_%d", idxCent)));
874
875 if (!sp0 || !sp1) {
876 AliError(Form("Effective correction objects 1 - idxCent %d can not be retrieved!",idxCent));
877 return -1;
878 }
879
880 fCorr1[0][idxCent] = static_cast<THnSparseF*>(sp0->Clone());
881 fCorr1[1][idxCent] = static_cast<THnSparseF*>(sp1->Clone());
882 }
883
884 return 0;
885}
886
887/*
888 * ---------------------------------------------------------------------------------
889 * Event / Trigger Statistics - private
890 * ---------------------------------------------------------------------------------
891 */
892
893//________________________________________________________________________
894Bool_t AliAnalysisNetParticleHelper::FillEventStats(Int_t *aEventCuts) {
895 // -- Fill event / centrality statistics
896
897 Bool_t isRejected = kFALSE;
898
899 // -- Fill event statistics
900 for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
9be43c8e 901
cb68eb1d 902 if (aEventCuts[idx])
903 isRejected = kTRUE;
904 else
905 fHEventStat0->Fill(idx);
906 }
907
908 for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
909 if (aEventCuts[idx])
910 break;
911 fHEventStat1->Fill(idx);
912 }
913
914 // -- Fill centrality statistics of accepted events
915 if (!isRejected)
916 fHCentralityStat->Fill(fCentralityBin);
917
918 return isRejected;
919}