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