]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
changed default binning
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskBFPsi.cxx
CommitLineData
c683985a 1#include "TChain.h"
2#include "TList.h"
3#include "TCanvas.h"
4#include "TLorentzVector.h"
5#include "TGraphErrors.h"
6#include "TH1F.h"
7#include "TH2F.h"
8#include "TH3F.h"
9#include "TH2D.h"
10#include "TH3D.h"
11#include "TArrayF.h"
12#include "TF1.h"
13#include "TRandom.h"
3b408588 14#include "TROOT.h"
c683985a 15
16#include "AliAnalysisTaskSE.h"
17#include "AliAnalysisManager.h"
18
19#include "AliESDVertex.h"
20#include "AliESDEvent.h"
21#include "AliESDInputHandler.h"
22#include "AliAODEvent.h"
23#include "AliAODTrack.h"
24#include "AliAODInputHandler.h"
25#include "AliAODMCParticle.h"
26#include "AliCollisionGeometry.h"
27#include "AliGenEventHeader.h"
28#include "AliMCEventHandler.h"
29#include "AliMCEvent.h"
30#include "AliStack.h"
31#include "AliESDtrackCuts.h"
32#include "AliEventplane.h"
33#include "AliTHn.h"
34#include "AliLog.h"
35#include "AliAnalysisUtils.h"
36
37#include "AliEventPoolManager.h"
38
39#include "AliPID.h"
40#include "AliPIDResponse.h"
41#include "AliPIDCombined.h"
42
43#include "AliAnalysisTaskBFPsi.h"
44#include "AliBalancePsi.h"
45#include "AliAnalysisTaskTriggeredBF.h"
46
47
48// Analysis task for the BF vs Psi code
49// Authors: Panos.Christakoglou@nikhef.nl
50
51using std::cout;
52using std::endl;
53
54ClassImp(AliAnalysisTaskBFPsi)
55
56//________________________________________________________________________
57AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
58: AliAnalysisTaskSE(name),
59 fDebugLevel(kFALSE),
60 fArrayMC(0), //+++++++++++++
61 fBalance(0),
62 fRunShuffling(kFALSE),
63 fShuffledBalance(0),
64 fRunMixing(kFALSE),
65 fRunMixingEventPlane(kFALSE),
66 fMixingTracks(50000),
67 fMixedBalance(0),
68 fPoolMgr(0),
69 fList(0),
70 fListBF(0),
71 fListBFS(0),
72 fListBFM(0),
73 fHistListPIDQA(0),
74 fHistEventStats(0),
75 fHistCentStats(0),
76 fHistCentStatsUsed(0),
77 fHistTriggerStats(0),
78 fHistTrackStats(0),
79 fHistVx(0),
80 fHistVy(0),
81 fHistVz(0),
82 fHistTPCvsVZEROMultiplicity(0),
83 fHistVZEROSignal(0),
84 fHistEventPlane(0),
85 fHistClus(0),
86 fHistDCA(0),
87 fHistChi2(0),
88 fHistPt(0),
89 fHistEta(0),
90 fHistRapidity(0),
91 fHistPhi(0),
92 fHistEtaPhiPos(0),
93 fHistEtaPhiNeg(0),
94 fHistPhiBefore(0),
95 fHistPhiAfter(0),
96 fHistPhiPos(0),
97 fHistPhiNeg(0),
98 fHistV0M(0),
99 fHistRefTracks(0),
100 fHistdEdxVsPTPCbeforePID(NULL),
101 fHistBetavsPTOFbeforePID(NULL),
102 fHistProbTPCvsPtbeforePID(NULL),
103 fHistProbTOFvsPtbeforePID(NULL),
104 fHistProbTPCTOFvsPtbeforePID(NULL),
105 fHistNSigmaTPCvsPtbeforePID(NULL),
106 fHistNSigmaTOFvsPtbeforePID(NULL),
107 fHistBetaVsdEdXbeforePID(NULL), //+++++++
108 fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++
109 fHistdEdxVsPTPCafterPID(NULL),
110 fHistBetavsPTOFafterPID(NULL),
111 fHistProbTPCvsPtafterPID(NULL),
112 fHistProbTOFvsPtafterPID(NULL),
113 fHistProbTPCTOFvsPtafterPID(NULL),
114 fHistNSigmaTPCvsPtafterPID(NULL),
115 fHistNSigmaTOFvsPtafterPID(NULL),
116 fHistBetaVsdEdXafterPID(NULL), //+++++++
117 fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++
118 fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++
119 fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++
120 fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++
121 fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++
122 fCentralityArrayBinsForCorrections(kCENTRALITY),
123 fCentralityWeights(0x0),
124 fPIDResponse(0x0),
125 fPIDCombined(0x0),
126 fParticleOfInterest(kPion),
127 fPidDetectorConfig(kTPCTOF),
128 fUsePID(kFALSE),
129 fUsePIDnSigma(kTRUE),
130 fUsePIDPropabilities(kFALSE),
131 fPIDNSigma(3.),
132 fMinAcceptedPIDProbability(0.8),
133 fElectronRejection(kFALSE),
134 fElectronOnlyRejection(kFALSE),
135 fElectronRejectionNSigma(-1.),
136 fElectronRejectionMinPt(0.),
137 fElectronRejectionMaxPt(1000.),
138 fESDtrackCuts(0),
139 fCentralityEstimator("V0M"),
140 fUseCentrality(kFALSE),
141 fCentralityPercentileMin(0.),
142 fCentralityPercentileMax(5.),
143 fImpactParameterMin(0.),
144 fImpactParameterMax(20.),
145 fMultiplicityEstimator("V0A"),
146 fUseMultiplicity(kFALSE),
147 fNumberOfAcceptedTracksMin(0),
148 fNumberOfAcceptedTracksMax(10000),
149 fHistNumberOfAcceptedTracks(0),
150 fHistMultiplicity(0),
151 fUseOfflineTrigger(kFALSE),
152 fCheckFirstEventInChunk(kFALSE),
153 fCheckPileUp(kFALSE),
154 fCheckPrimaryFlagAOD(kFALSE),
155 fUseMCforKinematics(kFALSE),
156 fVxMax(0.3),
157 fVyMax(0.3),
158 fVzMax(10.),
159 fnAODtrackCutBit(128),
160 fPtMin(0.3),
161 fPtMax(1.5),
c683985a 162 fEtaMin(-0.8),
b66d1de3 163 fEtaMax(0.8),
c683985a 164 fPhiMin(0.),
165 fPhiMax(360.),
c683985a 166 fDCAxyCut(-1),
167 fDCAzCut(-1),
168 fTPCchi2Cut(-1),
169 fNClustersTPCCut(-1),
096141dd 170 fTPCsharedCut(-1),
c683985a 171 fAcceptanceParameterization(0),
172 fDifferentialV2(0),
173 fUseFlowAfterBurner(kFALSE),
174 fExcludeResonancesInMC(kFALSE),
175 fExcludeElectronsInMC(kFALSE),
176 fUseMCPdgCode(kFALSE),
177 fPDGCodeToBeAnalyzed(-1),
178 fEventClass("EventPlane"),
179 fCustomBinning(""),
180 fHistVZEROAGainEqualizationMap(0),
181 fHistVZEROCGainEqualizationMap(0),
182 fHistVZEROChannelGainEqualizationMap(0) {
183 // Constructor
184 // Define input and output slots here
185 // Input slot #0 works with a TChain
186
187 //======================================================correction
188 for (Int_t i=0; i<kCENTRALITY; i++){
189 fHistCorrectionPlus[i] = NULL;
190 fHistCorrectionMinus[i] = NULL;
191 fCentralityArrayForCorrections[i] = -1.;
192 }
193 //=====================================================correction
194
195 DefineInput(0, TChain::Class());
196 // Output slot #0 writes into a TH1 container
197 DefineOutput(1, TList::Class());
198 DefineOutput(2, TList::Class());
199 DefineOutput(3, TList::Class());
200 DefineOutput(4, TList::Class());
201 DefineOutput(5, TList::Class());
202}
203
204//________________________________________________________________________
205AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {
206
207 // delete fBalance;
208 // delete fShuffledBalance;
209 // delete fList;
210 // delete fListBF;
211 // delete fListBFS;
212
213 // delete fHistEventStats;
214 // delete fHistTrackStats;
215 // delete fHistVx;
216 // delete fHistVy;
217 // delete fHistVz;
218
219 // delete fHistClus;
220 // delete fHistDCA;
221 // delete fHistChi2;
222 // delete fHistPt;
223 // delete fHistEta;
224 // delete fHistPhi;
225 // delete fHistEtaPhiPos;
226 // delete fHistEtaPhiNeg;
227 // delete fHistV0M;
228}
229
230//________________________________________________________________________
231void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
232 // Create histograms
233 // Called once
234
235 // global switch disabling the reference
236 // (to avoid "Replacing existing TH1" if several wagons are created in train)
237 Bool_t oldStatus = TH1::AddDirectoryStatus();
238 TH1::AddDirectory(kFALSE);
239
240 if(!fBalance) {
241 fBalance = new AliBalancePsi();
242 fBalance->SetAnalysisLevel("ESD");
243 fBalance->SetEventClass(fEventClass);
244 //fBalance->SetNumberOfBins(-1,16);
245 //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
246 }
247 if(fRunShuffling) {
248 if(!fShuffledBalance) {
249 fShuffledBalance = new AliBalancePsi();
250 fShuffledBalance->SetAnalysisLevel("ESD");
251 fShuffledBalance->SetEventClass(fEventClass);
252 //fShuffledBalance->SetNumberOfBins(-1,16);
253 //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);
254 }
255 }
256 if(fRunMixing) {
257 if(!fMixedBalance) {
258 fMixedBalance = new AliBalancePsi();
259 fMixedBalance->SetAnalysisLevel("ESD");
260 fMixedBalance->SetEventClass(fEventClass);
261 }
262 }
263
264 //QA list
265 fList = new TList();
266 fList->SetName("listQA");
267 fList->SetOwner();
268
269 //Balance Function list
270 fListBF = new TList();
271 fListBF->SetName("listBF");
272 fListBF->SetOwner();
273
274 if(fRunShuffling) {
275 fListBFS = new TList();
276 fListBFS->SetName("listBFShuffled");
277 fListBFS->SetOwner();
278 }
279
280 if(fRunMixing) {
281 fListBFM = new TList();
282 fListBFM->SetName("listTriggeredBFMixed");
283 fListBFM->SetOwner();
284 }
285
286 //PID QA list
287 if(fUsePID || fElectronRejection) {
288 fHistListPIDQA = new TList();
289 fHistListPIDQA->SetName("listQAPID");
290 fHistListPIDQA->SetOwner();
291 }
292
293 //Event stats.
294 TString gCutName[7] = {"Total","Offline trigger",
295 "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};
296 fHistEventStats = new TH2F("fHistEventStats",
297 "Event statistics;;Centrality percentile;N_{events}",
298 7,0.5,7.5,220,-5,105);
299 for(Int_t i = 1; i <= 7; i++)
300 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
301 fList->Add(fHistEventStats);
302
303 TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
304 fHistCentStats = new TH2F("fHistCentStats",
305 "Centrality statistics;;Cent percentile",
306 13,-0.5,12.5,220,-5,105);
307 for(Int_t i = 1; i <= 13; i++){
308 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
309 //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
310 }
311 fList->Add(fHistCentStats);
312
313 fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);
314 fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());
315 fList->Add(fHistCentStatsUsed);
316
317 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);
318 fList->Add(fHistTriggerStats);
319
320 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);
321 fList->Add(fHistTrackStats);
322
323 fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);
324 fList->Add(fHistNumberOfAcceptedTracks);
325
326 fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);
327 fList->Add(fHistMultiplicity);
328
329 // Vertex distributions
330 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
331 fList->Add(fHistVx);
332 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
333 fList->Add(fHistVy);
334 fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);
335 fList->Add(fHistVz);
336
337 //TPC vs VZERO multiplicity
338 fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
339 if(fMultiplicityEstimator == "V0A")
340 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");
341 else if(fMultiplicityEstimator == "V0C")
342 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");
343 else
344 fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");
345 fList->Add(fHistTPCvsVZEROMultiplicity);
346
347 fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
348 fList->Add(fHistVZEROSignal);
349
350 //Event plane
351 fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);
352 fList->Add(fHistEventPlane);
353
354 // QA histograms
355 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
356 fList->Add(fHistClus);
357 fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);
358 fList->Add(fHistChi2);
359 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
360 fList->Add(fHistDCA);
361 fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);
362 fList->Add(fHistPt);
363 fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);
364 fList->Add(fHistEta);
365 fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);
366 fList->Add(fHistRapidity);
367 fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);
368 fList->Add(fHistPhi);
655d9d07 369 fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",40,-1.6,1.6,72,-TMath::Pi()/2.,1.5*TMath::Pi(),220,-5,105);
c683985a 370 fList->Add(fHistEtaPhiPos);
655d9d07 371 fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",40,-1.6,1.6,72,-TMath::Pi()/2.,1.5*TMath::Pi(),220,-5,105);
c683985a 372 fList->Add(fHistEtaPhiNeg);
373 fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
374 fList->Add(fHistPhiBefore);
375 fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
376 fList->Add(fHistPhiAfter);
377 fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);
378 fList->Add(fHistPhiPos);
379 fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);
380 fList->Add(fHistPhiNeg);
381 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
382 fList->Add(fHistV0M);
383 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
384 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
385 for(Int_t i = 1; i <= 6; i++)
386 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
387 fList->Add(fHistRefTracks);
388
389 // Balance function histograms
390 // Initialize histograms if not done yet (including the custom binning)
391 if(!fBalance->GetHistNp()){
392 AliInfo("Histograms not yet initialized! --> Will be done now");
393 fBalance->SetCustomBinning(fCustomBinning);
394 fBalance->InitHistograms();
395 }
396
397 if(fRunShuffling) {
398 if(!fShuffledBalance->GetHistNp()) {
399 AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");
400 fShuffledBalance->SetCustomBinning(fCustomBinning);
401 fShuffledBalance->InitHistograms();
402 }
403 }
404
405 if(fRunMixing) {
406 if(!fMixedBalance->GetHistNp()) {
407 AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");
408 fMixedBalance->SetCustomBinning(fCustomBinning);
409 fMixedBalance->InitHistograms();
410 }
411 }
412
413 // QA histograms for different cuts
414 fList->Add(fBalance->GetQAHistHBTbefore());
415 fList->Add(fBalance->GetQAHistHBTafter());
416 fList->Add(fBalance->GetQAHistConversionbefore());
417 fList->Add(fBalance->GetQAHistConversionafter());
418 fList->Add(fBalance->GetQAHistPsiMinusPhi());
419 fList->Add(fBalance->GetQAHistResonancesBefore());
420 fList->Add(fBalance->GetQAHistResonancesRho());
421 fList->Add(fBalance->GetQAHistResonancesK0());
422 fList->Add(fBalance->GetQAHistResonancesLambda());
423 fList->Add(fBalance->GetQAHistQbefore());
424 fList->Add(fBalance->GetQAHistQafter());
425
426 //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
427 fListBF->Add(fBalance->GetHistNp());
428 fListBF->Add(fBalance->GetHistNn());
429 fListBF->Add(fBalance->GetHistNpn());
430 fListBF->Add(fBalance->GetHistNnn());
431 fListBF->Add(fBalance->GetHistNpp());
432 fListBF->Add(fBalance->GetHistNnp());
433
434 if(fRunShuffling) {
435 fListBFS->Add(fShuffledBalance->GetHistNp());
436 fListBFS->Add(fShuffledBalance->GetHistNn());
437 fListBFS->Add(fShuffledBalance->GetHistNpn());
438 fListBFS->Add(fShuffledBalance->GetHistNnn());
439 fListBFS->Add(fShuffledBalance->GetHistNpp());
440 fListBFS->Add(fShuffledBalance->GetHistNnp());
441 }
442
443 if(fRunMixing) {
444 fListBFM->Add(fMixedBalance->GetHistNp());
445 fListBFM->Add(fMixedBalance->GetHistNn());
446 fListBFM->Add(fMixedBalance->GetHistNpn());
447 fListBFM->Add(fMixedBalance->GetHistNnn());
448 fListBFM->Add(fMixedBalance->GetHistNpp());
449 fListBFM->Add(fMixedBalance->GetHistNnp());
450 }
451 //}
452
453
454 // Event Mixing
455 if(fRunMixing){
456 Int_t trackDepth = fMixingTracks;
457 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
458
459 // centrality bins
460 Double_t* centbins = NULL;
461 Int_t nCentralityBins;
462 if(fBalance->IsUseVertexBinning()){
463 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);
464 }
465 else{
466 centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);
467 }
468
469 // multiplicity bins
470 Double_t* multbins = NULL;
471 Int_t nMultiplicityBins;
472 multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);
473
474 // Zvtx bins
475 Double_t* vtxbins = NULL;
476 Int_t nVertexBins;
477 if(fBalance->IsUseVertexBinning()){
478 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);
479 }
480 else{
481 vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);
482 }
483
484 // Event plane angle (Psi) bins
485 Double_t* psibins = NULL;
486 Int_t nPsiBins;
487 psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);
488
b78ae8b4 489
c683985a 490 // run the event mixing also in bins of event plane (statistics!)
491 if(fRunMixingEventPlane){
492 if(fEventClass=="Multiplicity"){
493 if(multbins && vtxbins && psibins){
494 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
495 }
496 }
497 else{
498 if(centbins && vtxbins && psibins){
499 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
500 }
501 }
502 }
503 else{
504 if(fEventClass=="Multiplicity"){
505 if(multbins && vtxbins){
506 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
507 }
508 }
509 else{
510 if(centbins && vtxbins){
511 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
512 }
513 }
514 }
b78ae8b4 515
516 if(centbins) delete [] centbins;
517 if(multbins) delete [] multbins;
518 if(vtxbins) delete [] vtxbins;
70de3ddf 519 if(psibins) delete [] psibins;
b78ae8b4 520
c683985a 521 // check pool manager
522 if(!fPoolMgr){
523 AliError("Event Mixing required, but Pool Manager not initialized...");
524 return;
525 }
526 }
527
528 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
529
530 //====================PID========================//
531 if(fUsePID) {
532 fPIDCombined = new AliPIDCombined();
533 fPIDCombined->SetDefaultTPCPriors();
534
535 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
536 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);
537
538 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
539 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID);
540
541 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
542 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID);
543
544 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
545 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);
546
547 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
548 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);
549
550 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
551 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);
552
553 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
554 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID);
555
556 fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2);
557 fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++
558
559 fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500);
560 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++
561
562 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
563 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);
564
565 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
566 fHistListPIDQA->Add(fHistBetavsPTOFafterPID);
567
568 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
569 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);
570
571 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2);
572 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID);
573
574 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
575 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);
576
577 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
578 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);
579
580 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
581 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);
582
583 fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2);
584 fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++
585
586 fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500);
587 fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++
588 }
589
590 // for electron rejection only TPC nsigma histograms
591 if(fElectronRejection) {
592
593 fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000);
594 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron);
595
596 fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500);
597 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron);
598
599 fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000);
600 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron);
601
602 fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500);
603 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron);
604 }
605 //====================PID========================//
606
607 // Post output data.
608 PostData(1, fList);
609 PostData(2, fListBF);
610 if(fRunShuffling) PostData(3, fListBFS);
611 if(fRunMixing) PostData(4, fListBFM);
612 if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID
613
614 AliInfo("Finished setting up the Output");
615
616 TH1::AddDirectory(oldStatus);
b78ae8b4 617
c683985a 618}
619
620
621//________________________________________________________________________
622void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename,
623 Int_t nCentralityBins,
624 Double_t *centralityArrayForCorrections) {
625 //Open files that will be used for correction
626 fCentralityArrayBinsForCorrections = nCentralityBins;
627 for (Int_t i=0; i<nCentralityBins; i++)
628 fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];
629
630 // No file specified -> run without corrections
631 if(!filename.Contains(".root")) {
632 AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));
633 return;
634 }
635
636 //Open the input file
637 TFile *f = TFile::Open(filename);
638 if(!f->IsOpen()) {
639 AliInfo(Form("File %s not found --> run without corrections",filename.Data()));
640 return;
641 }
642
643 //TString listEffName = "";
644 for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {
645 //Printf("iCent %d:",iCent);
646 TString histoName = "fHistCorrectionPlus";
647 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
648 fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));
649 if(!fHistCorrectionPlus[iCent]) {
650 AliError(Form("fHist %s not found!!!",histoName.Data()));
651 return;
652 }
653
654 histoName = "fHistCorrectionMinus";
655 histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));
656 fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data()));
657 if(!fHistCorrectionMinus[iCent]) {
658 AliError(Form("fHist %s not found!!!",histoName.Data()));
659 return;
660 }
661 }//loop over centralities: ONLY the PbPb case is covered
c683985a 662}
663
664//________________________________________________________________________
665void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
666 // Main loop
667 // Called for each event
668
669 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
670 Int_t gNumberOfAcceptedTracks = 0;
671 Double_t lMultiplicityVar = -999.; //-1
672 Double_t gReactionPlane = -1.;
673 Float_t bSign = 0.;
674
675 // get the event (for generator level: MCEvent())
676 AliVEvent* eventMain = NULL;
677 if(gAnalysisLevel == "MC") {
678 eventMain = dynamic_cast<AliVEvent*>(MCEvent());
679 }
680 else{
681 eventMain = dynamic_cast<AliVEvent*>(InputEvent());
682 // for HBT like cuts need magnetic field sign
683 bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;
684 }
685 if(!eventMain) {
686 AliError("eventMain not available");
687 return;
688 }
689
690 // PID Response task active?
691 if(fUsePID || fElectronRejection) {
692 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
693 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
694 }
695
696 // check event cuts and fill event histograms
697 if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){
698 return;
699 }
700 // get the reaction plane
3b408588 701 if(fEventClass != "Multiplicity" && gAnalysisLevel!="AODnano") {
c683985a 702 gReactionPlane = GetEventPlane(eventMain);
703 fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);
704 if(gReactionPlane < 0){
705 return;
706 }
707 }
708
709 // get the accepted tracks in main event
710 TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);
711 gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();
712
713 //multiplicity cut (used in pp)
714 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);
715
716 // store charges of all accepted tracks,shuffle and reassign(two extra loops)
717 TObjArray* tracksShuffled = NULL;
718 if(fRunShuffling){
719 tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);
720 }
721
722 // Event mixing
723 if (fRunMixing)
724 {
725 // 1. First get an event pool corresponding in mult (cent) and
726 // zvertex to the current event. Once initialized, the pool
727 // should contain nMix (reduced) events. This routine does not
728 // pre-scan the chain. The first several events of every chain
729 // will be skipped until the needed pools are filled to the
730 // specified depth. If the pool categories are not too rare, this
731 // should not be a problem. If they are rare, you could lose`
732 // statistics.
733
734 // 2. Collect the whole pool's content of tracks into one TObjArray
735 // (bgTracks), which is effectively a single background super-event.
736
737 // 3. The reduced and bgTracks arrays must both be passed into
738 // FillCorrelations(). Also nMix should be passed in, so a weight
739 // of 1./nMix can be applied.
740
741 AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);
742
743 if (!pool){
744 AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));
745 }
746 else{
747
748 //pool->SetDebug(1);
749
750 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
751
752
753 Int_t nMix = pool->GetCurrentNEvents();
754 //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
755
756 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
757 //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
758 //if (pool->IsReady())
759 //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
760
761 // Fill mixed-event histos here
762 for (Int_t jMix=0; jMix<nMix; jMix++)
763 {
764 TObjArray* tracksMixed = pool->GetEvent(jMix);
765 fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
766 }
767 }
768
769 // Update the Event pool
770 pool->UpdatePool(tracksMain);
771 //pool->PrintInfo();
772
773 }//pool NULL check
774 }//run mixing
775
776 // calculate balance function
777 fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
778
779 // calculate shuffled balance function
780 if(fRunShuffling && tracksShuffled != NULL) {
781 fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());
782 }
783}
784
785//________________________________________________________________________
786Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
787 // Checks the Event cuts
788 // Fills Event statistics histograms
789
790 Bool_t isSelectedMain = kTRUE;
791 Float_t gRefMultiplicity = -1.;
792 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
793
794 AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);
795
796 fHistEventStats->Fill(1,gRefMultiplicity); //all events
797
798 // check first event in chunk (is not needed for new reconstructions)
799 if(fCheckFirstEventInChunk){
800 AliAnalysisUtils ut;
801 if(ut.IsFirstEventInChunk(event))
802 return -1.;
803 fHistEventStats->Fill(6,gRefMultiplicity);
804 }
805 // check for pile-up event
806 if(fCheckPileUp){
807 AliAnalysisUtils ut;
808 ut.SetUseMVPlpSelection(kTRUE);
809 ut.SetUseOutOfBunchPileUp(kTRUE);
810 if(ut.IsPileUpEvent(event))
811 return -1.;
812 fHistEventStats->Fill(7,gRefMultiplicity);
813 }
814
815 // Event trigger bits
816 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
817 if(fUseOfflineTrigger)
818 isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
819
820 if(isSelectedMain) {
821 fHistEventStats->Fill(2,gRefMultiplicity); //triggered events
822
823 // Event Vertex MC
824 if(gAnalysisLevel == "MC") {
825 if(!event) {
826 AliError("mcEvent not available");
827 return 0x0;
828 }
829
830 if(mcevent){
831 AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());
832 if(header){
833 TArrayF gVertexArray;
834 header->PrimaryVertex(gVertexArray);
835 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
836 //gVertexArray.At(0),
837 //gVertexArray.At(1),
838 //gVertexArray.At(2));
839 fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex
840 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
841 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
842 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
843 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
844
845 // get the reference multiplicty or centrality
846 gRefMultiplicity = GetRefMultiOrCentrality(event);
847
848 fHistVx->Fill(gVertexArray.At(0));
849 fHistVy->Fill(gVertexArray.At(1));
850 fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);
851
852 // take only events inside centrality class
853 if(fUseCentrality) {
854 if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){
855 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
856 return gRefMultiplicity;
857 }//centrality class
858 }
859 // take events only within the same multiplicity class
860 else if(fUseMultiplicity) {
861 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
862 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
863 return gRefMultiplicity;
864 }
865 }//multiplicity range
866 }//Vz cut
867 }//Vy cut
868 }//Vx cut
869 }//header
870 }//MC event object
871 }//MC
872
873 // Event Vertex AOD, ESD, ESDMC
874 else{
875 const AliVVertex *vertex = event->GetPrimaryVertex();
876
877 if(vertex) {
878 Double32_t fCov[6];
879 vertex->GetCovarianceMatrix(fCov);
880 if(vertex->GetNContributors() > 0) {
881 if(fCov[5] != 0) {
882 fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex
883 if(TMath::Abs(vertex->GetX()) < fVxMax) {
884 if(TMath::Abs(vertex->GetY()) < fVyMax) {
885 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
886 fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events
887
888 // get the reference multiplicty or centrality
889 gRefMultiplicity = GetRefMultiOrCentrality(event);
890
891 fHistVx->Fill(vertex->GetX());
892 fHistVy->Fill(vertex->GetY());
893 fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);
894
895 // take only events inside centrality class
896 if(fUseCentrality) {
897 if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){
898
899 // centrality weighting (optional for 2011 if central and semicentral triggers are used)
900 if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){
901 AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
902 return -1;
903 }
904
905 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality
906 return gRefMultiplicity;
907 }//centrality class
908 }
909 // take events only within the same multiplicity class
910 else if(fUseMultiplicity) {
911 //if(fDebugLevel)
912 //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,
913 //fNumberOfAcceptedTracksMax,gRefMultiplicity);
914
915 if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {
916 fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity
917 return gRefMultiplicity;
918 }
919 }//multiplicity range
920 }//Vz cut
921 }//Vy cut
922 }//Vx cut
923 }//proper vertex resolution
924 }//proper number of contributors
925 }//vertex object valid
926 }//triggered event
927 }//AOD,ESD,ESDMC
928
929 // in all other cases return -1 (event not accepted)
930 return -1;
931}
932
933
934//________________________________________________________________________
935Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
936 // Checks the Event cuts
937 // Fills Event statistics histograms
3b408588 938
c683985a 939 Float_t gCentrality = -1.;
940 Double_t gMultiplicity = -1.;
941 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
942
943
944 // calculate centrality always (not only in centrality mode)
945 if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++
946 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
947 if(header){
948 gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
949
950 // QA for centrality estimators
951 fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
952 fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));
953 fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));
954 fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
955 fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
956 fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
957 fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
958 fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
959 fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));
960 fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));
961 fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
962 fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
963 fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
964
965 // Centrality estimator USED ++++++++++++++++++++++++++++++
966 fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));
967
968 // centrality QA (V0M)
969 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
970
971 // centrality QA (reference tracks)
972 fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
973 fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
974 fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
975 fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
976 fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
977 fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
978 fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
979 fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
980 fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
981
982 }//AOD header
983 }//AOD
3b408588 984
985 // calculate centrality always (not only in centrality mode)
986 else if(gAnalysisLevel == "AODnano" ) { //centrality via JF workaround
987
988 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
989 if(header){
990 gCentrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", header,fCentralityEstimator.Data())) / 100 - 1.0;
991
992 // QA histogram
993 fHistCentStatsUsed->Fill(0.,gCentrality);
994
995 }//AOD header
996 }//AODnano
c683985a 997
998 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs
999 AliCentrality *centrality = event->GetCentrality();
1000 gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
1001
1002 // QA for centrality estimators
1003 fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));
1004 fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));
1005 fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));
1006 fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));
1007 fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));
1008 fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));
1009 fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));
1010 fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));
1011 fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));
1012 fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));
1013 fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));
1014 fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));
1015 fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));
1016
1017 // Centrality estimator USED ++++++++++++++++++++++++++++++
1018 fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
1019
1020 // centrality QA (V0M)
1021 fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
1022 }//ESD
1023
1024 else if(gAnalysisLevel == "MC"){
1025 Double_t gImpactParameter = 0.;
1026 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1027 if(gMCEvent){
1028 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1029 if(headerH){
1030 gImpactParameter = headerH->ImpactParameter();
1031 gCentrality = gImpactParameter;
1032 }//MC header
1033 }//MC event cast
1034 }//MC
1035
1036 else{
1037 gCentrality = -1.;
1038 }
1039
1040 // calculate reference multiplicity always (not only in multiplicity mode)
1041 if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){
1042 AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);
1043 if(gESDEvent){
1044 gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1045 fHistMultiplicity->Fill(gMultiplicity);
1046 }//AliESDevent cast
1047 }//ESD mode
1048
1049 else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){
1050 AliAODHeader *header = (AliAODHeader*) event->GetHeader();
1051 if ((fMultiplicityEstimator == "V0M")||
1052 (fMultiplicityEstimator == "V0A")||
1053 (fMultiplicityEstimator == "V0C") ||
1054 (fMultiplicityEstimator == "TPC")) {
1055 gMultiplicity = GetReferenceMultiplicityFromAOD(event);
1056 if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);
1057 }
1058 else {
1059 if(header)
1060 gMultiplicity = header->GetRefMultiplicity();
1061 if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);
1062 }
1063 fHistMultiplicity->Fill(gMultiplicity);
1064 }//AOD mode
1065 else if(gAnalysisLevel == "MC") {
1066 AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);
1067 //Calculating the multiplicity as the number of charged primaries
1068 //within \pm 0.8 in eta and pT > 0.1 GeV/c
1069 for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {
1070 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));
1071 if (!track) {
1072 AliError(Form("Could not receive particle %d", iParticle));
1073 continue;
1074 }
1075
1076 //exclude non stable particles
1077 if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;
1078
1079 //++++++++++++++++
1080 if (fMultiplicityEstimator == "V0M") {
1081 if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7))
1082 continue;}
1083 else if (fMultiplicityEstimator == "V0A") {
1084 if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}
1085 else if (fMultiplicityEstimator == "V0C") {
1086 if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}
1087 else if (fMultiplicityEstimator == "TPC") {
1088 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1089 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1090 }
1091 else{
1092 if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;
1093 if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;
1094 }
1095 //++++++++++++++++
1096
1097 if(track->Charge() == 0) continue;
1098
1099 gMultiplicity += 1;
1100 }//loop over primaries
1101 fHistMultiplicity->Fill(gMultiplicity);
1102 }//MC mode
1103 else{
1104 gMultiplicity = -1;
1105 }
1106
1107
1108 // decide what should be returned only here
1109 Double_t lReturnVal = -100;
1110 if(fEventClass=="Multiplicity"){
1111 lReturnVal = gMultiplicity;
1112 }else if(fEventClass=="Centrality"){
1113 lReturnVal = gCentrality;
1114 }
1115 return lReturnVal;
1116}
1117
1118//________________________________________________________________________
1119Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){
1120 //Function that returns the reference multiplicity from AODs (data or reco MC)
1121 //Different ref. mult. implemented: V0M, V0A, V0C, TPC
1122 Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
1123 Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
1124
1125 AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
1126 if(!header) {
1127 Printf("ERROR: AOD header not available");
1128 return -999;
1129 }
1130 Int_t gRunNumber = header->GetRunNumber();
1131
1132 // Loop over tracks in event
1133 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1134 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1135 if (!aodTrack) {
1136 AliError(Form("Could not receive track %d", iTracks));
1137 continue;
1138 }
1139
1140 // AOD track cuts
1141 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1142
1143 if(aodTrack->Charge() == 0) continue;
1144 // Kinematics cuts from ESD track cuts
1145 if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;
1146 if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;
1147
1148 //=================PID (so far only for electron rejection)==========================//
1149 if(fElectronRejection) {
1150 // get the electron nsigma
1151 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1152
1153 // check only for given momentum range
1154 if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){
1155 //look only at electron nsigma
1156 if(!fElectronOnlyRejection) {
1157 //Make the decision based on the n-sigma of electrons
1158 if(nSigma < fElectronRejectionNSigma) continue;
1159 }
1160 else {
1161 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1162 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1163 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1164
1165 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1166 if(nSigma < fElectronRejectionNSigma
1167 && nSigmaPions > fElectronRejectionNSigma
1168 && nSigmaKaons > fElectronRejectionNSigma
1169 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1170 }
1171 }
1172 }//electron rejection
1173
1174 gRefMultiplicityTPC += 1.0;
1175 }// track loop
1176
1177 //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
1178 for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
1179 fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
1180
1181 if(iChannel < 32)
1182 gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
1183 else if(iChannel >= 32)
1184 gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
1185 }//loop over PMTs
1186
1187 //Equalization of gain
1188 Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
1189 if(gFactorA != 0)
1190 gRefMultiplicityVZEROA /= gFactorA;
1191 Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
1192 if(gFactorC != 0)
1193 gRefMultiplicityVZEROC /= gFactorC;
1194 if((gFactorA != 0)&&(gFactorC != 0))
1195 gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
1196
1197 if(fDebugLevel)
1198 Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);
1199
1200 fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
1201
1202 if(fMultiplicityEstimator == "TPC")
1203 gRefMultiplicity = gRefMultiplicityTPC;
1204 else if(fMultiplicityEstimator == "V0M")
1205 gRefMultiplicity = gRefMultiplicityVZERO;
1206 else if(fMultiplicityEstimator == "V0A")
1207 gRefMultiplicity = gRefMultiplicityVZEROA;
1208 else if(fMultiplicityEstimator == "V0C")
1209 gRefMultiplicity = gRefMultiplicityVZEROC;
1210
1211 return gRefMultiplicity;
1212}
1213
1214//________________________________________________________________________
1215Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){
1216 // Get the event plane
1217
1218 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1219
1220 Float_t gVZEROEventPlane = -10.;
1221 Float_t gReactionPlane = -10.;
1222 Double_t qxTot = 0.0, qyTot = 0.0;
1223
1224 //MC: from reaction plane
1225 if(gAnalysisLevel == "MC"){
1226 if(!event) {
1227 AliError("mcEvent not available");
1228 return 0x0;
1229 }
1230
1231 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
1232 if(gMCEvent){
1233 AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());
1234 if (headerH) {
1235 gReactionPlane = headerH->ReactionPlaneAngle();
1236 //gReactionPlane *= TMath::RadToDeg();
1237 }//MC header
1238 }//MC event cast
1239 }//MC
1240
1241 // AOD,ESD,ESDMC: from VZERO Event Plane
1242 else{
1243
1244 AliEventplane *ep = event->GetEventplane();
1245 if(ep){
1246 gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
1247 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
1248 //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
1249 gReactionPlane = gVZEROEventPlane;
1250 }
1251 }//AOD,ESD,ESDMC
1252
1253 return gReactionPlane;
1254}
1255
1256//________________________________________________________________________
1257Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta,
1258 Double_t vPhi,
1259 Double_t vPt,
1260 Short_t vCharge,
1261 Double_t gCentrality) {
1262 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)
1263
1264 Double_t correction = 1.;
c683985a 1265 Int_t gCentralityInt = -1;
967e4db5 1266
c683985a 1267 for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){
1268 if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){
1269 gCentralityInt = i;
1270 break;
1271 }
1272 }
1273
1274 // centrality not in array --> no correction
1275 if(gCentralityInt < 0){
1276 correction = 1.;
1277 }
1278 else{
1279
1280 //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);
1281
1282 if(fHistCorrectionPlus[gCentralityInt]){
1283 if (vCharge > 0) {
967e4db5 1284 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
c683985a 1285 //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);
1286 }
1287 if (vCharge < 0) {
967e4db5 1288 correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi));
1289 //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);
c683985a 1290 }
1291 }
1292 else {
1293 correction = 1.;
1294 }
1295 }//centrality in array
1296
1297 if (correction == 0.) {
1298 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt));
1299 correction = 1.;
1300 }
1301
1302 return correction;
1303}
1304
1305//________________________________________________________________________
1306TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){
1307 // Returns TObjArray with tracks after all track cuts (only for AOD!)
1308 // Fills QA histograms
1309
1310 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1311
1312 //output TObjArray holding all good tracks
1313 TObjArray* tracksAccepted = new TObjArray;
1314 tracksAccepted->SetOwner(kTRUE);
1315
1316 Short_t vCharge;
1317 Double_t vEta;
1318 Double_t vY;
1319 Double_t vPhi;
1320 Double_t vPt;
1321
1322
1323 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD
1324 // Loop over tracks in event
1325
1326 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1327 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1328 if (!aodTrack) {
1329 AliError(Form("Could not receive track %d", iTracks));
1330 continue;
1331 }
1332
1333 // AOD track cuts
1334
1335 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
1336 // take only TPC only tracks
1337 //fHistTrackStats->Fill(aodTrack->GetFilterMap());
1338 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1339 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1340 }
1341
1342 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1343
1344
1345 // additional check on kPrimary flag
1346 if(fCheckPrimaryFlagAOD){
1347 if(aodTrack->GetType() != AliAODTrack::kPrimary)
1348 continue;
1349 }
1350
1351
1352 vCharge = aodTrack->Charge();
1353 vEta = aodTrack->Eta();
1354 vY = aodTrack->Y();
1355 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1356 vPt = aodTrack->Pt();
1357
1358 //===========================PID (so far only for electron rejection)===============================//
1359 if(fElectronRejection) {
1360
1361 // get the electron nsigma
1362 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1363
1364 //Fill QA before the PID
1365 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1366 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1367 //end of QA-before pid
1368
1369 // check only for given momentum range
1370 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1371
1372 //look only at electron nsigma
1373 if(!fElectronOnlyRejection){
1374
1375 //Make the decision based on the n-sigma of electrons
1376 if(nSigma < fElectronRejectionNSigma) continue;
1377 }
1378 else{
1379
1380 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1381 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1382 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1383
1384 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1385 if(nSigma < fElectronRejectionNSigma
1386 && nSigmaPions > fElectronRejectionNSigma
1387 && nSigmaKaons > fElectronRejectionNSigma
1388 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1389 }
1390 }
1391
1392 //Fill QA after the PID
1393 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1394 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1395
1396 }
1397 //===========================end of PID (so far only for electron rejection)===============================//
1398
1399 //+++++++++++++++++++++++++++++//
1400 //===========================PID===============================//
1401 if(fUsePID) {
1402 Double_t prob[AliPID::kSPECIES]={0.};
1403 Double_t probTPC[AliPID::kSPECIES]={0.};
1404 Double_t probTOF[AliPID::kSPECIES]={0.};
1405 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1406
1407 Double_t nSigma = 0.;
1408 Double_t nSigmaTPC = 0.; //++++
1409 Double_t nSigmaTOF = 0.; //++++
1410 Double_t nSigmaTPCTOF = 0.; //++++
1411 UInt_t detUsedTPC = 0;
1412 UInt_t detUsedTOF = 0;
1413 UInt_t detUsedTPCTOF = 0;
1414
1415 nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1416 nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));
1417 nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);
1418
1419 //Decide what detector configuration we want to use
1420 switch(fPidDetectorConfig) {
1421 case kTPCpid:
1422 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1423 nSigma = nSigmaTPC;
1424 detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
1425 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1426 prob[iSpecies] = probTPC[iSpecies];
1427 break;
1428 case kTOFpid:
1429 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1430 nSigma = nSigmaTOF;
1431 detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
1432 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1433 prob[iSpecies] = probTOF[iSpecies];
1434 break;
1435 case kTPCTOF:
1436 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1437 nSigma = nSigmaTPCTOF;
1438 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);
1439 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1440 prob[iSpecies] = probTPCTOF[iSpecies];
1441 break;
1442 default:
1443 break;
1444 }//end switch: define detector mask
1445
1446 //Filling the PID QA
1447 Double_t tofTime = -999., length = 999., tof = -999.;
1448 Double_t c = TMath::C()*1.E-9;// m/ns
1449 Double_t beta = -999.;
1450 if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&
1451 (aodTrack->IsOn(AliAODTrack::kTIME)) ) {
1452 tofTime = aodTrack->GetTOFsignal();//in ps
1453 length = aodTrack->GetIntegratedLength();
1454 tof = tofTime*1E-3; // ns
1455
1456 if (tof <= 0) {
1457 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
1458 continue;
1459 }
1460 if (length <= 0){
1461 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
1462 continue;
1463 }
1464
1465 length = length*0.01; // in meters
1466 tof = tof*c;
1467 beta = length/tof;
1468
1469 fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1470 fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1471 fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF);
1472 }//TOF signal
1473
1474 fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1475 fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1476 fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC);
1477
1478 fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1479
1480 //combined TPC&TOF
1481 fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
1482 fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF);
1483 Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF);
1484
1485 //end of QA-before pid
1486
1487 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
1488 //Make the decision based on the n-sigma
1489 if(fUsePIDnSigma) {
1490 if(nSigma > fPIDNSigma) continue;
1491
1492 fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF);
1493 fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC);
1494 fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF);
1495 }
1496 //Make the decision based on the bayesian
1497 else if(fUsePIDPropabilities) {
1498 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
1499 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
1500
1501 fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);
1502 fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]);
1503 fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
1504
1505 }
1506
1507 //Fill QA after the PID
1508 fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
1509 fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1510 fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++
1511 }
1512 }
1513 //===========================PID===============================//
1514 //+++++++++++++++++++++++++++++//
1515
1516
1517 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1518 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1519
1520
1521 // Kinematics cuts from ESD track cuts
1522 if( vPt < fPtMin || vPt > fPtMax) continue;
1523 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1524
1525 // Extra DCA cuts (for systematic studies [!= -1])
1526 if( fDCAxyCut != -1 && fDCAzCut != -1){
1527 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1528 continue; // 2D cut
1529 }
1530 }
1531
1532 // Extra TPC cuts (for systematic studies [!= -1])
1533 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1534 continue;
1535 }
1536 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1537 continue;
1538 }
096141dd 1539
1540 // Extra cut on shared clusters
1541 if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
1542 continue;
1543 }
c683985a 1544
1545 // fill QA histograms
1546 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1547 fHistDCA->Fill(dcaZ,dcaXY);
1548 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1549 fHistPt->Fill(vPt,gCentrality);
1550 fHistEta->Fill(vEta,gCentrality);
1551 fHistRapidity->Fill(vY,gCentrality);
1552 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1553 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1554 fHistPhi->Fill(vPhi,gCentrality);
1555 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1556 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1557
1558 //=======================================correction
1559 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1560 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1561
1562 // add the track to the TObjArray
1563 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1564 }//track loop
1565 }// AOD analysis
1566
b66d1de3 1567
3b408588 1568 // nano AODs
1569 else if(gAnalysisLevel == "AODnano") { // not fully supported yet (PID missing)
1570 // Loop over tracks in event
1571
1572 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
b66d1de3 1573 AliVTrack* aodTrack = dynamic_cast<AliVTrack *>(event->GetTrack(iTracks));
3b408588 1574 if (!aodTrack) {
1575 AliError(Form("Could not receive track %d", iTracks));
1576 continue;
1577 }
1578
b66d1de3 1579 // AOD track cuts (not needed)
1580 //if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
3b408588 1581
1582 vCharge = aodTrack->Charge();
1583 vEta = aodTrack->Eta();
b66d1de3 1584 vY = -999.;
3b408588 1585 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1586 vPt = aodTrack->Pt();
b66d1de3 1587
3b408588 1588
1589 // Kinematics cuts from ESD track cuts
1590 if( vPt < fPtMin || vPt > fPtMax) continue;
1591 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1592
b66d1de3 1593
3b408588 1594 // fill QA histograms
3b408588 1595 fHistPt->Fill(vPt,gCentrality);
1596 fHistEta->Fill(vEta,gCentrality);
1597 fHistRapidity->Fill(vY,gCentrality);
1598 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1599 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1600 fHistPhi->Fill(vPhi,gCentrality);
1601 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1602 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1603
1604 //=======================================correction
1605 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1606 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1607
1608 // add the track to the TObjArray
1609 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1610 }//track loop
1611 }// AOD nano analysis
1612
1613
c683985a 1614 //==============================================================================================================
1615 else if(gAnalysisLevel == "MCAOD") {
1616
1617 AliMCEvent* mcEvent = MCEvent();
1618 if (!mcEvent) {
1619 AliError("ERROR: Could not retrieve MC event");
1620 }
1621 else{
1622
1623 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
1624 AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks);
1625 if (!aodTrack) {
1626 AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));
1627 continue;
1628 }
1629
1630 if(!aodTrack->IsPhysicalPrimary()) continue;
1631
1632 vCharge = aodTrack->Charge();
1633 vEta = aodTrack->Eta();
1634 vY = aodTrack->Y();
1635 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1636 vPt = aodTrack->Pt();
1637
1638 // Kinematics cuts from ESD track cuts
1639 if( vPt < fPtMin || vPt > fPtMax) continue;
1640 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1641
1642 // Remove neutral tracks
1643 if( vCharge == 0 ) continue;
1644
1645 //Exclude resonances
1646 if(fExcludeResonancesInMC) {
1647
1648 Bool_t kExcludeParticle = kFALSE;
1649 Int_t gMotherIndex = aodTrack->GetMother();
1650 if(gMotherIndex != -1) {
1651 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1652 if(motherTrack) {
1653 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1654 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1655 //if(pdgCodeOfMother == 113) {
1656 if(pdgCodeOfMother == 113 // rho0
1657 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1658 // || pdgCodeOfMother == 221 // eta
1659 // || pdgCodeOfMother == 331 // eta'
1660 // || pdgCodeOfMother == 223 // omega
1661 // || pdgCodeOfMother == 333 // phi
1662 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1663 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1664 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1665 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1666 || pdgCodeOfMother == 111 // pi0 Dalitz
1667 || pdgCodeOfMother == 22 // photon
1668 ) {
1669 kExcludeParticle = kTRUE;
1670 }
1671 }
1672 }
1673
1674 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1675 if(kExcludeParticle) continue;
1676 }
1677
1678 //Exclude electrons with PDG
1679 if(fExcludeElectronsInMC) {
1680
1681 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
1682
1683 }
1684
1685 // fill QA histograms
1686 fHistPt->Fill(vPt,gCentrality);
1687 fHistEta->Fill(vEta,gCentrality);
1688 fHistRapidity->Fill(vY,gCentrality);
1689 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1690 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1691 fHistPhi->Fill(vPhi,gCentrality);
1692 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1693 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1694
1695 //=======================================correction
1696 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1697 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1698
1699 // add the track to the TObjArray
1700 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1701 }//aodTracks
1702 }//MC event
1703 }//MCAOD
1704 //==============================================================================================================
1705
1706 //==============================================================================================================
1707 else if(gAnalysisLevel == "MCAODrec") {
1708
1709 /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
1710 if (!fAOD) {
1711 printf("ERROR: fAOD not available\n");
1712 return;
1713 }*/
1714
1715 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
1716 if (!fArrayMC) {
1717 AliError("No array of MC particles found !!!");
1718 }
1719
1720 AliMCEvent* mcEvent = MCEvent();
1721 if (!mcEvent) {
1722 AliError("ERROR: Could not retrieve MC event");
1723 return tracksAccepted;
1724 }
1725
1726 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1727 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1728 if (!aodTrack) {
1729 AliError(Form("Could not receive track %d", iTracks));
1730 continue;
1731 }
1732
1733 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1734 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1735 }
1736 if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;
1737
1738 vCharge = aodTrack->Charge();
1739 vEta = aodTrack->Eta();
1740 vY = aodTrack->Y();
1741 vPhi = aodTrack->Phi();// * TMath::RadToDeg();
1742 vPt = aodTrack->Pt();
1743
1744 //===========================use MC information for Kinematics===============================//
1745 if(fUseMCforKinematics){
1746
1747 Int_t label = TMath::Abs(aodTrack->GetLabel());
1748 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1749
1750 if(AODmcTrack){
1751 vCharge = AODmcTrack->Charge();
1752 vEta = AODmcTrack->Eta();
1753 vY = AODmcTrack->Y();
1754 vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();
1755 vPt = AODmcTrack->Pt();
1756 }
1757 else{
1758 AliDebug(1, "no MC particle for this track");
1759 continue;
1760 }
1761 }
1762 //===========================end of use MC information for Kinematics========================//
1763
1764
1765 //===========================PID (so far only for electron rejection)===============================//
1766 if(fElectronRejection) {
1767
1768 // get the electron nsigma
1769 Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));
1770
1771 //Fill QA before the PID
1772 fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1773 fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1774 //end of QA-before pid
1775
1776 // check only for given momentum range
1777 if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){
1778
1779 //look only at electron nsigma
1780 if(!fElectronOnlyRejection){
1781
1782 //Make the decision based on the n-sigma of electrons
1783 if(nSigma < fElectronRejectionNSigma) continue;
1784 }
1785 else{
1786
1787 Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));
1788 Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));
1789 Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));
1790
1791 //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)
1792 if(nSigma < fElectronRejectionNSigma
1793 && nSigmaPions > fElectronRejectionNSigma
1794 && nSigmaKaons > fElectronRejectionNSigma
1795 && nSigmaProtons > fElectronRejectionNSigma ) continue;
1796 }
1797 }
1798
1799 //Fill QA after the PID
1800 fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());
1801 fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma);
1802
1803 }
1804 //===========================end of PID (so far only for electron rejection)===============================//
1805
1806 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)
1807 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)
1808
1809 // Kinematics cuts from ESD track cuts
1810 if( vPt < fPtMin || vPt > fPtMax) continue;
1811 if( vEta < fEtaMin || vEta > fEtaMax) continue;
1812
1813 // Extra DCA cuts (for systematic studies [!= -1])
1814 if( fDCAxyCut != -1 && fDCAzCut != -1){
1815 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
1816 continue; // 2D cut
1817 }
1818 }
1819
1820 // Extra TPC cuts (for systematic studies [!= -1])
1821 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
1822 continue;
1823 }
1824 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
1825 continue;
1826 }
1827
096141dd 1828 // Extra cut on shared clusters
1829 if( fTPCsharedCut != -1 && aodTrack->GetTPCnclsS() > fTPCsharedCut){
1830 continue;
1831 }
1832
c683985a 1833 //Exclude resonances
1834 if(fExcludeResonancesInMC) {
1835
1836 Bool_t kExcludeParticle = kFALSE;
1837
1838 Int_t label = TMath::Abs(aodTrack->GetLabel());
1839 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1840
1841 if (AODmcTrack){
1842 //if (AODmcTrack->IsPhysicalPrimary()){
1843
1844 Int_t gMotherIndex = AODmcTrack->GetMother();
1845 if(gMotherIndex != -1) {
1846 AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1847 if(motherTrack) {
1848 Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
1849 if(pdgCodeOfMother == 113 // rho0
1850 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
1851 // || pdgCodeOfMother == 221 // eta
1852 // || pdgCodeOfMother == 331 // eta'
1853 // || pdgCodeOfMother == 223 // omega
1854 // || pdgCodeOfMother == 333 // phi
1855 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
1856 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
1857 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
1858 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
1859 || pdgCodeOfMother == 111 // pi0 Dalitz
1860 || pdgCodeOfMother == 22 // photon
1861 ) {
1862 kExcludeParticle = kTRUE;
1863 }
1864 }
1865 }
1866 }
1867 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1868 if(kExcludeParticle) continue;
1869 }
1870
1871 //Exclude electrons with PDG
1872 if(fExcludeElectronsInMC) {
1873
1874 Int_t label = TMath::Abs(aodTrack->GetLabel());
1875 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
1876
1877 if (AODmcTrack){
1878 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
1879 }
1880 }
1881
1882 // fill QA histograms
1883 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
1884 fHistDCA->Fill(dcaZ,dcaXY);
1885 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);
1886 fHistPt->Fill(vPt,gCentrality);
1887 fHistEta->Fill(vEta,gCentrality);
1888 fHistRapidity->Fill(vY,gCentrality);
1889 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
1890 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
1891 fHistPhi->Fill(vPhi,gCentrality);
1892 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
1893 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
1894
1895 //=======================================correction
1896 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
1897 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
1898
1899 // add the track to the TObjArray
1900 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
1901 }//track loop
1902 }//MCAODrec
1903 //==============================================================================================================
1904
1905 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD
1906
1907 AliESDtrack *trackTPC = NULL;
1908 AliMCParticle *mcTrack = NULL;
1909
1910 AliMCEvent* mcEvent = NULL;
1911
1912 // for MC ESDs use also MC information
1913 if(gAnalysisLevel == "MCESD"){
1914 mcEvent = MCEvent();
1915 if (!mcEvent) {
1916 AliError("mcEvent not available");
1917 return tracksAccepted;
1918 }
1919 }
1920
1921 // Loop over tracks in event
1922 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
1923 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));
1924 if (!track) {
1925 AliError(Form("Could not receive track %d", iTracks));
1926 continue;
1927 }
1928
1929 // for MC ESDs use also MC information --> MC track not used further???
1930 if(gAnalysisLevel == "MCESD"){
1931 Int_t label = TMath::Abs(track->GetLabel());
1932 if(label > mcEvent->GetNumberOfTracks()) continue;
1933 if(label > mcEvent->GetNumberOfPrimaries()) continue;
1934
1935 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1936 if(!mcTrack) continue;
1937 }
1938
1939 // take only TPC only tracks
1940 trackTPC = new AliESDtrack();
1941 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1942
1943 //ESD track cuts
1944 if(fESDtrackCuts)
1945 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1946
1947 // fill QA histograms
1948 Float_t b[2];
1949 Float_t bCov[3];
1950 trackTPC->GetImpactParameters(b,bCov);
1951 if (bCov[0]<=0 || bCov[2]<=0) {
1952 AliDebug(1, "Estimated b resolution lower or equal zero!");
1953 bCov[0]=0; bCov[2]=0;
1954 }
1955
1956 Int_t nClustersTPC = -1;
1957 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone
1958 //nClustersTPC = track->GetTPCclusters(0); // global track
1959 Float_t chi2PerClusterTPC = -1;
1960 if (nClustersTPC!=0) {
1961 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone
1962 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track
1963 }
1964
1965 //===========================PID===============================//
1966 if(fUsePID) {
1967 Double_t prob[AliPID::kSPECIES]={0.};
1968 Double_t probTPC[AliPID::kSPECIES]={0.};
1969 Double_t probTOF[AliPID::kSPECIES]={0.};
1970 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1971
1972 Double_t nSigma = 0.;
1973 UInt_t detUsedTPC = 0;
1974 UInt_t detUsedTOF = 0;
1975 UInt_t detUsedTPCTOF = 0;
1976
1977 //Decide what detector configuration we want to use
1978 switch(fPidDetectorConfig) {
1979 case kTPCpid:
1980 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1981 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));
1982 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);
1983 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1984 prob[iSpecies] = probTPC[iSpecies];
1985 break;
1986 case kTOFpid:
1987 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1988 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));
1989 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);
1990 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1991 prob[iSpecies] = probTOF[iSpecies];
1992 break;
1993 case kTPCTOF:
1994 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1995 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);
1996 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1997 prob[iSpecies] = probTPCTOF[iSpecies];
1998 break;
1999 default:
2000 break;
2001 }//end switch: define detector mask
2002
2003 //Filling the PID QA
2004 Double_t tofTime = -999., length = 999., tof = -999.;
2005 Double_t c = TMath::C()*1.E-9;// m/ns
2006 Double_t beta = -999.;
2007 Double_t nSigmaTOFForParticleOfInterest = -999.;
2008 if ( (track->IsOn(AliESDtrack::kTOFin)) &&
2009 (track->IsOn(AliESDtrack::kTIME)) ) {
2010 tofTime = track->GetTOFsignal();//in ps
2011 length = track->GetIntegratedLength();
2012 tof = tofTime*1E-3; // ns
2013
2014 if (tof <= 0) {
2015 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
2016 continue;
2017 }
2018 if (length <= 0){
2019 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
2020 continue;
2021 }
2022
2023 length = length*0.01; // in meters
2024 tof = tof*c;
2025 beta = length/tof;
2026
2027 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);
2028 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
2029 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
2030 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
2031 }//TOF signal
2032
2033
2034 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);
2035 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
2036 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
2037 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
2038 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
2039 //end of QA-before pid
2040
2041 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
2042 //Make the decision based on the n-sigma
2043 if(fUsePIDnSigma) {
2044 if(nSigma > fPIDNSigma) continue;}
2045
2046 //Make the decision based on the bayesian
2047 else if(fUsePIDPropabilities) {
2048 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
2049 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
2050 }
2051
2052 //Fill QA after the PID
2053 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);
2054 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);
2055 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);
2056
2057 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());
2058 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]);
2059 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);
2060 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest);
2061 }
2062 }
2063 //===========================PID===============================//
2064 vCharge = trackTPC->Charge();
2065 vY = trackTPC->Y();
2066 vEta = trackTPC->Eta();
2067 vPhi = trackTPC->Phi();// * TMath::RadToDeg();
2068 vPt = trackTPC->Pt();
2069 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
2070 fHistDCA->Fill(b[1],b[0]);
2071 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);
2072 fHistPt->Fill(vPt,gCentrality);
2073 fHistEta->Fill(vEta,gCentrality);
2074 fHistPhi->Fill(vPhi,gCentrality);
2075 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2076 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2077 fHistRapidity->Fill(vY,gCentrality);
2078 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2079 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2080
2081 //=======================================correction
2082 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2083 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2084
2085 // add the track to the TObjArray
2086 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2087
2088 delete trackTPC;
2089 }//track loop
2090 }// ESD analysis
2091
2092 else if(gAnalysisLevel == "MC"){
2093 if(!event) {
2094 AliError("mcEvent not available");
2095 return 0x0;
2096 }
2097
2098 AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);
2099 if(gMCEvent) {
2100 // Loop over tracks in event
2101 for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {
2102 AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));
2103 if (!track) {
2104 AliError(Form("Could not receive particle %d", iTracks));
2105 continue;
2106 }
2107
2108 //exclude non stable particles
2109 if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;
2110
2111 vCharge = track->Charge();
2112 vEta = track->Eta();
2113 vPt = track->Pt();
2114 vY = track->Y();
2115
2116 if( vPt < fPtMin || vPt > fPtMax)
2117 continue;
2118 if (!fUsePID) {
2119 if( vEta < fEtaMin || vEta > fEtaMax) continue;
2120 }
2121 else if (fUsePID){
2122 if( vY < fEtaMin || vY > fEtaMax) continue;
2123 }
2124
2125 // Remove neutral tracks
2126 if( vCharge == 0 ) continue;
2127
2128 //analyze one set of particles
2129 if(fUseMCPdgCode) {
2130 TParticle *particle = track->Particle();
2131 if(!particle) continue;
2132
2133 Int_t gPdgCode = particle->GetPdgCode();
2134 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
2135 continue;
2136 }
2137
2138 //Use the acceptance parameterization
2139 if(fAcceptanceParameterization) {
2140 Double_t gRandomNumber = gRandom->Rndm();
2141 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
2142 continue;
2143 }
2144
2145 //Exclude resonances
2146 if(fExcludeResonancesInMC) {
2147 TParticle *particle = track->Particle();
2148 if(!particle) continue;
2149
2150 Bool_t kExcludeParticle = kFALSE;
2151 Int_t gMotherIndex = particle->GetFirstMother();
2152 if(gMotherIndex != -1) {
2153 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));
2154 if(motherTrack) {
2155 TParticle *motherParticle = motherTrack->Particle();
2156 if(motherParticle) {
2157 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
2158 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
2159 if(pdgCodeOfMother == 113 // rho0
2160 || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+
2161 // || pdgCodeOfMother == 221 // eta
2162 // || pdgCodeOfMother == 331 // eta'
2163 // || pdgCodeOfMother == 223 // omega
2164 // || pdgCodeOfMother == 333 // phi
2165 || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0
2166 // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*
2167 // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*
2168 || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda
2169 || pdgCodeOfMother == 111 // pi0 Dalitz
2170 ) {
2171 kExcludeParticle = kTRUE;
2172 }
2173 }
2174 }
2175 }
2176
2177 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
2178 if(kExcludeParticle) continue;
2179 }
2180
2181 //Exclude electrons with PDG
2182 if(fExcludeElectronsInMC) {
2183
2184 TParticle *particle = track->Particle();
2185
2186 if (particle){
2187 if(TMath::Abs(particle->GetPdgCode()) == 11) continue;
2188 }
2189 }
2190
2191 vPhi = track->Phi();
2192 //Printf("phi (before): %lf",vPhi);
2193
2194 fHistPt->Fill(vPt,gCentrality);
2195 fHistEta->Fill(vEta,gCentrality);
2196 fHistPhi->Fill(vPhi,gCentrality);
2197 if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);
2198 else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);
2199 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2200 fHistRapidity->Fill(vY,gCentrality);
2201 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2202 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);
2203 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);
2204 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);
2205
2206 //Flow after burner
2207 if(fUseFlowAfterBurner) {
2208 Double_t precisionPhi = 0.001;
2209 Int_t maxNumberOfIterations = 100;
2210
2211 Double_t phi0 = vPhi;
2212 Double_t gV2 = fDifferentialV2->Eval(vPt);
2213
2214 for (Int_t j = 0; j < maxNumberOfIterations; j++) {
2215 Double_t phiprev = vPhi;
2216 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2217 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad()));
2218 vPhi -= fl/fp;
2219 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
2220 }
2221 //Printf("phi (after): %lf\n",vPhi);
2222 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();
2223 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
2224 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);
2225
2226 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();
2227 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
2228 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);
2229
2230 }
2231
2232 //vPhi *= TMath::RadToDeg();
2233
2234 //=======================================correction
2235 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);
2236 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2237
2238 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));
2239 } //track loop
2240 }//MC event object
2241 }//MC
2242
2243 return tracksAccepted;
2244}
2245
2246//________________________________________________________________________
2247TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){
2248 // Clones TObjArray and returns it with tracks after shuffling the charges
2249
2250 TObjArray* tracksShuffled = new TObjArray;
2251 tracksShuffled->SetOwner(kTRUE);
2252
2253 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks
2254
2255 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2256 {
2257 AliVParticle* track = (AliVParticle*) tracks->At(i);
2258 chargeVector->push_back(track->Charge());
2259 }
2260
2261 random_shuffle(chargeVector->begin(), chargeVector->end());
2262
2263 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
2264 AliVParticle* track = (AliVParticle*) tracks->At(i);
2265 //==============================correction
2266 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);
2267 //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);
2268 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));
2269 }
2270
2271 delete chargeVector;
2272
2273 return tracksShuffled;
2274}
2275
2276//________________________________________________________________________
2277void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,
2278 const char* lhcPeriod) {
2279 //Function to setup the VZERO gain equalization
2280 //============Get the equilization map============//
2281 TFile *calibrationFile = TFile::Open(filename);
2282 if((!calibrationFile)||(!calibrationFile->IsOpen())) {
2283 Printf("No calibration file found!!!");
2284 return;
2285 }
2286
2287 TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
2288 if(!list) {
2289 Printf("Calibration TList not found!!!");
2290 return;
2291 }
2292
2293 fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
2294 if(!fHistVZEROAGainEqualizationMap) {
2295 Printf("VZERO-A calibration object not found!!!");
2296 return;
2297 }
2298 fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
2299 if(!fHistVZEROCGainEqualizationMap) {
2300 Printf("VZERO-C calibration object not found!!!");
2301 return;
2302 }
2303
2304 fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
2305 if(!fHistVZEROChannelGainEqualizationMap) {
2306 Printf("VZERO channel calibration object not found!!!");
2307 return;
2308 }
2309}
2310
2311//________________________________________________________________________
2312Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run,
2313 Int_t channel) {
2314 //
2315 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2316
2317 for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
2318 Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2319 if(gRunNumber == run)
2320 return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
2321 }
2322
2323 return 1.0;
2324}
2325
2326//________________________________________________________________________
2327Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run,
2328 const char* side) {
2329 //
2330 if(!fHistVZEROAGainEqualizationMap) return 1.0;
2331
2332 TString gVZEROSide = side;
2333 for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
2334 Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
2335 //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
2336 if(gRunNumber == run) {
2337 if(gVZEROSide == "A")
2338 return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
2339 else if(gVZEROSide == "C")
2340 return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
2341 }
2342 }
2343
2344 return 1.0;
2345}
2346
2347//____________________________________________________________________
2348Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality)
2349{
2350 // copied from AliAnalysisTaskPhiCorrelations
2351 //
2352 // rejects "randomly" events such that the centrality gets flat
2353 // uses fCentralityWeights histogram
2354
2355 // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
2356
2357 Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
2358 Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
2359
2360 Bool_t result = kFALSE;
2361 if (centralityDigits < weight)
2362 result = kTRUE;
2363
2364 AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
2365
2366 return result;
2367}
2368
2369//________________________________________________________________________
2370void AliAnalysisTaskBFPsi::FinishTaskOutput(){
2371 //Printf("END BF");
2372
2373 if (!fBalance) {
2374 AliError("fBalance not available");
2375 return;
2376 }
2377 if(fRunShuffling) {
2378 if (!fShuffledBalance) {
2379 AliError("fShuffledBalance not available");
2380 return;
2381 }
2382 }
2383
2384}
2385
2386//________________________________________________________________________
2387void AliAnalysisTaskBFPsi::Terminate(Option_t *) {
2388 // Draw result to the screen
2389 // Called once at the end of the query
2390
2391 // not implemented ...
2392
2393}
2394