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