fixing psi in MC from the header
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskBFPsi.cxx
CommitLineData
0879e280 1#include "TChain.h"\r
2#include "TList.h"\r
3#include "TCanvas.h"\r
4#include "TLorentzVector.h"\r
5#include "TGraphErrors.h"\r
6#include "TH1F.h"\r
7#include "TH2F.h"\r
8#include "TH2D.h" \r
9#include "TH3D.h"\r
10#include "TArrayF.h"\r
11#include "TF1.h"\r
12#include "TRandom.h"\r
13\r
14#include "AliAnalysisTaskSE.h"\r
15#include "AliAnalysisManager.h"\r
16\r
17#include "AliESDVertex.h"\r
18#include "AliESDEvent.h"\r
19#include "AliESDInputHandler.h"\r
20#include "AliAODEvent.h"\r
21#include "AliAODTrack.h"\r
22#include "AliAODInputHandler.h"\r
23#include "AliGenEventHeader.h"\r
24#include "AliGenHijingEventHeader.h"\r
25#include "AliMCEventHandler.h"\r
26#include "AliMCEvent.h"\r
27#include "AliStack.h"\r
28#include "AliESDtrackCuts.h"\r
29#include "AliEventplane.h"\r
30\r
31#include "AliPID.h" \r
32#include "AliPIDResponse.h" \r
33#include "AliPIDCombined.h" \r
34\r
35#include "AliAnalysisTaskBFPsi.h"\r
36#include "AliBalancePsi.h"\r
37\r
38\r
39// Analysis task for the BF vs Psi code\r
40// Authors: Panos.Christakoglou@nikhef.nl\r
41\r
42ClassImp(AliAnalysisTaskBFPsi)\r
43\r
44//________________________________________________________________________\r
45AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
46: AliAnalysisTaskSE(name), \r
47 fBalance(0),\r
48 fRunShuffling(kFALSE),\r
49 fShuffledBalance(0),\r
50 fList(0),\r
51 fListBF(0),\r
52 fListBFS(0),\r
53 fHistListPIDQA(0),\r
54 fHistEventStats(0),\r
55 fHistCentStats(0),\r
56 fHistTriggerStats(0),\r
57 fHistTrackStats(0),\r
58 fHistVx(0),\r
59 fHistVy(0),\r
60 fHistVz(0),\r
61 fHistEventPlane(0),\r
62 fHistClus(0),\r
63 fHistDCA(0),\r
64 fHistChi2(0),\r
65 fHistPt(0),\r
66 fHistEta(0),\r
67 fHistPhi(0),\r
68 fHistPhiBefore(0),\r
69 fHistPhiAfter(0),\r
70 fHistV0M(0),\r
71 fHistRefTracks(0),\r
72 fHistdEdxVsPTPCbeforePID(NULL),\r
73 fHistBetavsPTOFbeforePID(NULL), \r
74 fHistProbTPCvsPtbeforePID(NULL), \r
75 fHistProbTOFvsPtbeforePID(NULL), \r
76 fHistProbTPCTOFvsPtbeforePID(NULL),\r
77 fHistNSigmaTPCvsPtbeforePID(NULL), \r
78 fHistNSigmaTOFvsPtbeforePID(NULL), \r
79 fHistdEdxVsPTPCafterPID(NULL),\r
80 fHistBetavsPTOFafterPID(NULL), \r
81 fHistProbTPCvsPtafterPID(NULL), \r
82 fHistProbTOFvsPtafterPID(NULL), \r
83 fHistProbTPCTOFvsPtafterPID(NULL),\r
84 fHistNSigmaTPCvsPtafterPID(NULL), \r
85 fHistNSigmaTOFvsPtafterPID(NULL), \r
86 fPIDResponse(0x0),\r
87 fPIDCombined(0x0),\r
88 fParticleOfInterest(kPion),\r
89 fPidDetectorConfig(kTPCTOF),\r
90 fUsePID(kFALSE),\r
91 fUsePIDnSigma(kTRUE),\r
92 fUsePIDPropabilities(kFALSE), \r
93 fPIDNSigma(3.),\r
94 fMinAcceptedPIDProbability(0.8),\r
95 fESDtrackCuts(0),\r
96 fCentralityEstimator("V0M"),\r
97 fUseCentrality(kFALSE),\r
98 fCentralityPercentileMin(0.), \r
99 fCentralityPercentileMax(5.),\r
100 fImpactParameterMin(0.),\r
101 fImpactParameterMax(20.),\r
102 fUseMultiplicity(kFALSE),\r
103 fNumberOfAcceptedTracksMin(0),\r
104 fNumberOfAcceptedTracksMax(10000),\r
105 fHistNumberOfAcceptedTracks(0),\r
106 fUseOfflineTrigger(kFALSE),\r
107 fVxMax(0.3),\r
108 fVyMax(0.3),\r
109 fVzMax(10.),\r
110 nAODtrackCutBit(128),\r
111 fPtMin(0.3),\r
112 fPtMax(1.5),\r
113 fEtaMin(-0.8),\r
114 fEtaMax(-0.8),\r
115 fDCAxyCut(-1),\r
116 fDCAzCut(-1),\r
117 fTPCchi2Cut(-1),\r
118 fNClustersTPCCut(-1),\r
119 fAcceptanceParameterization(0),\r
120 fDifferentialV2(0),\r
121 fUseFlowAfterBurner(kFALSE),\r
122 fExcludeResonancesInMC(kFALSE),\r
123 fUseMCPdgCode(kFALSE),\r
124 fPDGCodeToBeAnalyzed(-1) {\r
125 // Constructor\r
126 // Define input and output slots here\r
127 // Input slot #0 works with a TChain\r
128 DefineInput(0, TChain::Class());\r
129 // Output slot #0 writes into a TH1 container\r
130 DefineOutput(1, TList::Class());\r
131 DefineOutput(2, TList::Class());\r
132 DefineOutput(3, TList::Class());\r
133 DefineOutput(4, TList::Class());\r
134}\r
135\r
136//________________________________________________________________________\r
137AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r
138\r
139 // delete fBalance; \r
140 // delete fShuffledBalance; \r
141 // delete fList;\r
142 // delete fListBF; \r
143 // delete fListBFS;\r
144\r
145 // delete fHistEventStats; \r
146 // delete fHistTrackStats; \r
147 // delete fHistVx; \r
148 // delete fHistVy; \r
149 // delete fHistVz; \r
150\r
151 // delete fHistClus;\r
152 // delete fHistDCA;\r
153 // delete fHistChi2;\r
154 // delete fHistPt;\r
155 // delete fHistEta;\r
156 // delete fHistPhi;\r
157 // delete fHistV0M;\r
158}\r
159\r
160//________________________________________________________________________\r
161void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
162 // Create histograms\r
163 // Called once\r
164 if(!fBalance) {\r
165 fBalance = new AliBalancePsi();\r
166 fBalance->SetAnalysisLevel("ESD");\r
167 //fBalance->SetNumberOfBins(-1,16);\r
168 fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
169 }\r
170 if(fRunShuffling) {\r
171 if(!fShuffledBalance) {\r
172 fShuffledBalance = new AliBalancePsi();\r
173 fShuffledBalance->SetAnalysisLevel("ESD");\r
174 //fShuffledBalance->SetNumberOfBins(-1,16);\r
175 fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
176 }\r
177 }\r
178\r
179 //QA list\r
180 fList = new TList();\r
181 fList->SetName("listQA");\r
182 fList->SetOwner();\r
183\r
184 //Balance Function list\r
185 fListBF = new TList();\r
186 fListBF->SetName("listBF");\r
187 fListBF->SetOwner();\r
188\r
189 if(fRunShuffling) {\r
190 fListBFS = new TList();\r
191 fListBFS->SetName("listBFShuffled");\r
192 fListBFS->SetOwner();\r
193 }\r
194\r
195 //PID QA list\r
196 if(fUsePID) {\r
197 fHistListPIDQA = new TList();\r
198 fHistListPIDQA->SetName("listQAPID");\r
199 fHistListPIDQA->SetOwner();\r
200 }\r
201\r
202 //Event stats.\r
203 TString gCutName[4] = {"Total","Offline trigger",\r
204 "Vertex","Analyzed"};\r
205 fHistEventStats = new TH1F("fHistEventStats",\r
206 "Event statistics;;N_{events}",\r
207 4,0.5,4.5);\r
208 for(Int_t i = 1; i <= 4; i++)\r
209 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
210 fList->Add(fHistEventStats);\r
211\r
212 TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
213 fHistCentStats = new TH2F("fHistCentStats",\r
214 "Centrality statistics;;Cent percentile",\r
215 9,-0.5,8.5,220,-5,105);\r
216 for(Int_t i = 1; i <= 9; i++)\r
217 fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
218 fList->Add(fHistCentStats);\r
219\r
220 fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
221 fList->Add(fHistTriggerStats);\r
222\r
223 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
224 fList->Add(fHistTrackStats);\r
225\r
226 fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
227 fList->Add(fHistNumberOfAcceptedTracks);\r
228\r
229 // Vertex distributions\r
230 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
231 fList->Add(fHistVx);\r
232 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
233 fList->Add(fHistVy);\r
234 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
235 fList->Add(fHistVz);\r
236\r
237 //Event plane\r
f2583ca7 238 fHistEventPlane = new TH1F("fHistEventPlane",";#Psi_{2} [deg.];Counts",100,0,360.);\r
0879e280 239 fList->Add(fHistEventPlane);\r
240\r
241 // QA histograms\r
242 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
243 fList->Add(fHistClus);\r
244 fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
245 fList->Add(fHistChi2);\r
246 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
247 fList->Add(fHistDCA);\r
248 fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
249 fList->Add(fHistPt);\r
250 fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
251 fList->Add(fHistEta);\r
252 fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
253 fList->Add(fHistPhi);\r
254 fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
255 fList->Add(fHistPhiBefore);\r
256 fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
257 fList->Add(fHistPhiAfter);\r
258 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
259 fList->Add(fHistV0M);\r
260 TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
261 fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
262 for(Int_t i = 1; i <= 6; i++)\r
263 fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
264 fList->Add(fHistRefTracks);\r
265\r
266 // Balance function histograms\r
267 // Initialize histograms if not done yet\r
268 if(!fBalance->GetHistNp(0)){\r
269 AliWarning("Histograms not yet initialized! --> Will be done now");\r
270 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
271 fBalance->InitHistograms();\r
272 }\r
273\r
274 if(fRunShuffling) {\r
275 if(!fShuffledBalance->GetHistNp(0)) {\r
276 AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
277 AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
278 fShuffledBalance->InitHistograms();\r
279 }\r
280 }\r
281\r
282 for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
283 fListBF->Add(fBalance->GetHistNp(a));\r
284 fListBF->Add(fBalance->GetHistNn(a));\r
285 fListBF->Add(fBalance->GetHistNpn(a));\r
286 fListBF->Add(fBalance->GetHistNnn(a));\r
287 fListBF->Add(fBalance->GetHistNpp(a));\r
288 fListBF->Add(fBalance->GetHistNnp(a));\r
289\r
290 if(fRunShuffling) {\r
291 fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
292 fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
293 fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
294 fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
295 fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
296 fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
297 } \r
298 }\r
299\r
300 if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
301\r
302 //====================PID========================//\r
303 if(fUsePID) {\r
304 fPIDCombined = new AliPIDCombined();\r
305 fPIDCombined->SetDefaultTPCPriors();\r
306\r
307 fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
308 fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
309 \r
310 fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
311 fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
312 \r
313 fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
314 fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
315 \r
316 fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
317 fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
318\r
319 fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
320 fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
321 \r
322 fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
323 fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
324 \r
325 fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
326 fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
327 \r
328 fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
329 fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
330 \r
331 fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
332 fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
333 \r
334 fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
335 fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
336 \r
337 fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
338 fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition \r
339 \r
340 fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
341 fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
342\r
343 fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
344 fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition \r
345 \r
346 fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
347 fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
348 }\r
349 //====================PID========================//\r
350\r
351 // Post output data.\r
352 PostData(1, fList);\r
353 PostData(2, fListBF);\r
354 if(fRunShuffling) PostData(3, fListBFS);\r
355 if(fUsePID) PostData(4, fHistListPIDQA); //PID\r
356}\r
357\r
358//________________________________________________________________________\r
359void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
360 // Main loop\r
361 // Called for each event\r
362 TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
363\r
364 AliESDtrack *track_TPC = NULL;\r
365\r
366 Int_t gNumberOfAcceptedTracks = 0;\r
367 Float_t fCentrality = 0.;\r
368 Double_t gReactionPlane = 0.;\r
369\r
370 // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
371 vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled\r
372 vector<Double_t> *chargeVector[9]; // original charge\r
373 for(Int_t i = 0; i < 9; i++){\r
374 chargeVectorShuffle[i] = new vector<Double_t>;\r
375 chargeVector[i] = new vector<Double_t>;\r
376 }\r
377\r
378 Double_t v_charge;\r
379 Double_t v_y;\r
380 Double_t v_eta;\r
381 Double_t v_phi;\r
382 Double_t v_p[3];\r
383 Double_t v_pt;\r
384 Double_t v_E;\r
385\r
386 if(fUsePID) {\r
387 fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
388 if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
389 }\r
390 \r
391 //ESD analysis\r
392 if(gAnalysisLevel == "ESD") {\r
393 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
394 if (!gESD) {\r
395 Printf("ERROR: gESD not available");\r
396 return;\r
397 }\r
398\r
399 // store offline trigger bits\r
400 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
401\r
402 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
403 fHistEventStats->Fill(1); //all events\r
404 Bool_t isSelected = kTRUE;\r
405 if(fUseOfflineTrigger)\r
406 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
407 if(isSelected) {\r
408 fHistEventStats->Fill(2); //triggered events\r
409\r
410 if(fUseCentrality) {\r
411 //Centrality stuff\r
412 AliCentrality *centrality = gESD->GetCentrality();\r
413 \r
414 fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
415 \r
416 // take only events inside centrality class\r
417 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
418 fCentralityPercentileMax,\r
419 fCentralityEstimator.Data()))\r
420 return;\r
421 \r
422 // centrality QA (V0M)\r
423 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
424 }\r
425 \r
426 const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
427 if(vertex) {\r
428 if(vertex->GetNContributors() > 0) {\r
429 if(vertex->GetZRes() != 0) {\r
430 fHistEventStats->Fill(3); //events with a proper vertex\r
431 if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
432 if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
433 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
434 fHistEventStats->Fill(4); //analayzed events\r
435 fHistVx->Fill(vertex->GetXv());\r
436 fHistVy->Fill(vertex->GetYv());\r
437 fHistVz->Fill(vertex->GetZv());\r
438 \r
439 //========Get the VZERO event plane========//\r
440 Double_t gVZEROEventPlane = -10.0;\r
441 Double_t qxTot = 0.0, qyTot = 0.0;\r
442 AliEventplane *ep = gESD->GetEventplane();\r
443 if(ep) \r
444 gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
445 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
446 fHistEventPlane->Fill(gVZEROEventPlane);\r
f2583ca7 447 gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
0879e280 448 //========Get the VZERO event plane========//\r
449\r
450 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
451 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
452 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
453 if (!track) {\r
454 Printf("ERROR: Could not receive track %d", iTracks);\r
455 continue;\r
456 } \r
457 \r
458 // take only TPC only tracks\r
459 track_TPC = new AliESDtrack();\r
460 if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
461 \r
462 //ESD track cuts\r
463 if(fESDtrackCuts) \r
464 if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
465 \r
466 // fill QA histograms\r
467 Float_t b[2];\r
468 Float_t bCov[3];\r
469 track_TPC->GetImpactParameters(b,bCov);\r
470 if (bCov[0]<=0 || bCov[2]<=0) {\r
471 AliDebug(1, "Estimated b resolution lower or equal zero!");\r
472 bCov[0]=0; bCov[2]=0;\r
473 }\r
474 \r
475 Int_t nClustersTPC = -1;\r
476 nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone\r
477 //nClustersTPC = track->GetTPCclusters(0); // global track\r
478 Float_t chi2PerClusterTPC = -1;\r
479 if (nClustersTPC!=0) {\r
480 chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
481 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
482 }\r
483\r
484 //===========================PID===============================// \r
485 if(fUsePID) {\r
486 Double_t prob[AliPID::kSPECIES]={0.};\r
487 Double_t probTPC[AliPID::kSPECIES]={0.};\r
488 Double_t probTOF[AliPID::kSPECIES]={0.};\r
489 Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
490\r
491 Double_t nSigma = 0.;\r
492 UInt_t detUsedTPC = 0;\r
493 UInt_t detUsedTOF = 0;\r
494 UInt_t detUsedTPCTOF = 0;\r
495\r
496 //Decide what detector configuration we want to use\r
497 switch(fPidDetectorConfig) {\r
498 case kTPCpid:\r
499 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
500 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
501 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
502 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
503 prob[iSpecies] = probTPC[iSpecies];\r
504 break;\r
505 case kTOFpid:\r
506 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
507 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
508 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
509 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
510 prob[iSpecies] = probTOF[iSpecies];\r
511 break;\r
512 case kTPCTOF:\r
513 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
514 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
515 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
516 prob[iSpecies] = probTPCTOF[iSpecies];\r
517 break;\r
518 default:\r
519 break;\r
520 }//end switch: define detector mask\r
521 \r
522 //Filling the PID QA\r
523 Double_t tofTime = -999., length = 999., tof = -999.;\r
524 Double_t c = TMath::C()*1.E-9;// m/ns\r
525 Double_t beta = -999.;\r
526 Double_t nSigmaTOFForParticleOfInterest = -999.;\r
527 if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
528 (track->IsOn(AliESDtrack::kTIME)) ) { \r
529 tofTime = track->GetTOFsignal();//in ps\r
530 length = track->GetIntegratedLength();\r
531 tof = tofTime*1E-3; // ns \r
532 \r
533 if (tof <= 0) {\r
534 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
535 continue;\r
536 }\r
537 if (length <= 0){\r
538 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
539 continue;\r
540 }\r
541 \r
542 length = length*0.01; // in meters\r
543 tof = tof*c;\r
544 beta = length/tof;\r
545 \r
546 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
547 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
548 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
549 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
550 }//TOF signal \r
551 \r
552 \r
553 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
554 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
555 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
556 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
557 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
558 //end of QA-before pid\r
559 \r
560 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
561 //Make the decision based on the n-sigma\r
562 if(fUsePIDnSigma) {\r
563 if(nSigma > fPIDNSigma) continue;}\r
564 \r
565 //Make the decision based on the bayesian\r
566 else if(fUsePIDPropabilities) {\r
567 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
568 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
569 }\r
570 \r
571 //Fill QA after the PID\r
572 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
573 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
574 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
575 \r
576 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
577 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
578 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
579 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
580 }\r
581 \r
582 PostData(4, fHistListPIDQA);\r
583 }\r
584 //===========================PID===============================//\r
585 v_charge = track_TPC->Charge();\r
586 v_y = track_TPC->Y();\r
587 v_eta = track_TPC->Eta();\r
588 v_phi = track_TPC->Phi() * TMath::RadToDeg();\r
589 v_E = track_TPC->E();\r
590 v_pt = track_TPC->Pt();\r
591 track_TPC->PxPyPz(v_p);\r
592 fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
593 fHistDCA->Fill(b[1],b[0]);\r
594 fHistChi2->Fill(chi2PerClusterTPC);\r
595 fHistPt->Fill(v_pt);\r
596 fHistEta->Fill(v_eta);\r
597 fHistPhi->Fill(v_phi);\r
598\r
599 // fill charge vector\r
600 chargeVector[0]->push_back(v_charge);\r
601 chargeVector[1]->push_back(v_y);\r
602 chargeVector[2]->push_back(v_eta);\r
603 chargeVector[3]->push_back(v_phi);\r
604 chargeVector[4]->push_back(v_p[0]);\r
605 chargeVector[5]->push_back(v_p[1]);\r
606 chargeVector[6]->push_back(v_p[2]);\r
607 chargeVector[7]->push_back(v_pt);\r
608 chargeVector[8]->push_back(v_E);\r
609\r
610 if(fRunShuffling) {\r
611 chargeVectorShuffle[0]->push_back(v_charge);\r
612 chargeVectorShuffle[1]->push_back(v_y);\r
613 chargeVectorShuffle[2]->push_back(v_eta);\r
614 chargeVectorShuffle[3]->push_back(v_phi);\r
615 chargeVectorShuffle[4]->push_back(v_p[0]);\r
616 chargeVectorShuffle[5]->push_back(v_p[1]);\r
617 chargeVectorShuffle[6]->push_back(v_p[2]);\r
618 chargeVectorShuffle[7]->push_back(v_pt);\r
619 chargeVectorShuffle[8]->push_back(v_E);\r
620 }\r
621 \r
622 delete track_TPC;\r
623 \r
624 } //track loop\r
625 }//Vz cut\r
626 }//Vy cut\r
627 }//Vx cut\r
628 }//proper vertex resolution\r
629 }//proper number of contributors\r
630 }//vertex object valid\r
631 }//triggered event \r
632 }//ESD analysis\r
633 \r
634 //AOD analysis (vertex and track cuts also here!!!!)\r
635 else if(gAnalysisLevel == "AOD") {\r
636 AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
637 if(!gAOD) {\r
638 Printf("ERROR: gAOD not available");\r
639 return;\r
640 }\r
641\r
642 AliAODHeader *aodHeader = gAOD->GetHeader();\r
643\r
644 // store offline trigger bits\r
645 fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
646 \r
647 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
648 fHistEventStats->Fill(1); //all events\r
649 Bool_t isSelected = kTRUE;\r
650 if(fUseOfflineTrigger)\r
651 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
652 if(isSelected) {\r
653 fHistEventStats->Fill(2); //triggered events\r
654 \r
655 //Centrality stuff (centrality in AOD header)\r
656 if(fUseCentrality) {\r
657 fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
658 \r
659 // QA for centrality estimators\r
660 fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
661 fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
662 fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
663 fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
664 fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
665 fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
666 fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
667 fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
668 fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
669 \r
670 // take only events inside centrality class\r
671 if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
672 return;\r
673 \r
674 // centrality QA (V0M)\r
675 fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
676 \r
677 // centrality QA (reference tracks)\r
678 fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
679 fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
680 fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());\r
681 fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());\r
682 fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));\r
683 fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));\r
684 fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
685 fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
686 fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
687 }\r
688\r
689 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
690 \r
691 if(vertex) {\r
692 Double32_t fCov[6];\r
693 vertex->GetCovarianceMatrix(fCov);\r
694 \r
695 if(vertex->GetNContributors() > 0) {\r
696 if(fCov[5] != 0) {\r
697 fHistEventStats->Fill(3); //events with a proper vertex\r
698 if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
699 if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
700 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
701 fHistEventStats->Fill(4); //analyzed events\r
702 fHistVx->Fill(vertex->GetX());\r
703 fHistVy->Fill(vertex->GetY());\r
704 fHistVz->Fill(vertex->GetZ());\r
705 \r
706 //========Get the VZERO event plane========//\r
707 Double_t gVZEROEventPlane = -10.0;\r
708 Double_t qxTot = 0.0, qyTot = 0.0;\r
709 AliEventplane *ep = gAOD->GetEventplane();\r
710 if(ep) \r
711 gVZEROEventPlane = ep->CalculateVZEROEventPlane(gAOD,10,2,qxTot,qyTot);\r
712 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
713 fHistEventPlane->Fill(gVZEROEventPlane);\r
f2583ca7 714 gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
0879e280 715 //========Get the VZERO event plane========//\r
716\r
717 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
718 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
719 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
720 if (!aodTrack) {\r
721 Printf("ERROR: Could not receive track %d", iTracks);\r
722 continue;\r
723 }\r
724 \r
725 // AOD track cuts\r
726 \r
727 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
728 // take only TPC only tracks \r
729 fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
730 if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
731 \r
732 v_charge = aodTrack->Charge();\r
733 v_y = aodTrack->Y();\r
734 v_eta = aodTrack->Eta();\r
735 v_phi = aodTrack->Phi() * TMath::RadToDeg();\r
736 v_E = aodTrack->E();\r
737 v_pt = aodTrack->Pt();\r
738 aodTrack->PxPyPz(v_p);\r
739 \r
740 Float_t DCAxy = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
741 Float_t DCAz = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
742 \r
743 \r
744 // Kinematics cuts from ESD track cuts\r
745 if( v_pt < fPtMin || v_pt > fPtMax) continue;\r
746 if( v_eta < fEtaMin || v_eta > fEtaMax) continue;\r
747 \r
748 // Extra DCA cuts (for systematic studies [!= -1])\r
749 if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
750 if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
751 continue; // 2D cut\r
752 }\r
753 }\r
754 \r
755 // Extra TPC cuts (for systematic studies [!= -1])\r
756 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
757 continue;\r
758 }\r
759 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
760 continue;\r
761 }\r
762 \r
763 // fill QA histograms\r
764 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
765 fHistDCA->Fill(DCAz,DCAxy);\r
766 fHistChi2->Fill(aodTrack->Chi2perNDF());\r
767 fHistPt->Fill(v_pt);\r
768 fHistEta->Fill(v_eta);\r
769 fHistPhi->Fill(v_phi);\r
770\r
771 // fill charge vector\r
772 chargeVector[0]->push_back(v_charge);\r
773 chargeVector[1]->push_back(v_y);\r
774 chargeVector[2]->push_back(v_eta);\r
775 chargeVector[3]->push_back(v_phi);\r
776 chargeVector[4]->push_back(v_p[0]);\r
777 chargeVector[5]->push_back(v_p[1]);\r
778 chargeVector[6]->push_back(v_p[2]);\r
779 chargeVector[7]->push_back(v_pt);\r
780 chargeVector[8]->push_back(v_E);\r
781\r
782 if(fRunShuffling) {\r
783 chargeVectorShuffle[0]->push_back(v_charge);\r
784 chargeVectorShuffle[1]->push_back(v_y);\r
785 chargeVectorShuffle[2]->push_back(v_eta);\r
786 chargeVectorShuffle[3]->push_back(v_phi);\r
787 chargeVectorShuffle[4]->push_back(v_p[0]);\r
788 chargeVectorShuffle[5]->push_back(v_p[1]);\r
789 chargeVectorShuffle[6]->push_back(v_p[2]);\r
790 chargeVectorShuffle[7]->push_back(v_pt);\r
791 chargeVectorShuffle[8]->push_back(v_E);\r
792 }\r
793 \r
794 gNumberOfAcceptedTracks += 1;\r
795 \r
796 } //track loop\r
797 }//Vz cut\r
798 }//Vy cut\r
799 }//Vx cut\r
800 }//proper vertex resolution\r
801 }//proper number of contributors\r
802 }//vertex object valid\r
803 }//triggered event \r
804 }//AOD analysis\r
805\r
806 //MC-ESD analysis\r
807 if(gAnalysisLevel == "MCESD") {\r
808 AliMCEvent* mcEvent = MCEvent(); \r
809 if (!mcEvent) {\r
810 Printf("ERROR: mcEvent not available");\r
811 return;\r
812 }\r
813\r
814 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
815 if (!gESD) {\r
816 Printf("ERROR: gESD not available");\r
817 return;\r
818 }\r
819\r
820 // store offline trigger bits\r
821 fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
822\r
823 // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
824 fHistEventStats->Fill(1); //all events\r
825 Bool_t isSelected = kTRUE;\r
826 if(fUseOfflineTrigger)\r
827 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
828 if(isSelected) {\r
829 fHistEventStats->Fill(2); //triggered events\r
830\r
831 if(fUseCentrality) {\r
832 //Centrality stuff\r
833 AliCentrality *centrality = gESD->GetCentrality();\r
834\r
835 fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
836 \r
837 // take only events inside centrality class\r
838 if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
839 fCentralityPercentileMax,\r
840 fCentralityEstimator.Data()))\r
841 return;\r
842 \r
843 // centrality QA (V0M)\r
844 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
845 }\r
846 \r
847 const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
848 if(vertex) {\r
849 if(vertex->GetNContributors() > 0) {\r
850 if(vertex->GetZRes() != 0) {\r
851 fHistEventStats->Fill(3); //events with a proper vertex\r
852 if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
853 if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
854 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
855 fHistEventStats->Fill(4); //analayzed events\r
856 fHistVx->Fill(vertex->GetXv());\r
857 fHistVy->Fill(vertex->GetYv());\r
858 fHistVz->Fill(vertex->GetZv());\r
859 \r
860 //========Get the VZERO event plane========//\r
861 Double_t gVZEROEventPlane = -10.0;\r
862 Double_t qxTot = 0.0, qyTot = 0.0;\r
863 AliEventplane *ep = gESD->GetEventplane();\r
864 if(ep) \r
865 gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
866 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
867 fHistEventPlane->Fill(gVZEROEventPlane);\r
f2583ca7 868 gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
0879e280 869 //========Get the VZERO event plane========//\r
870\r
871 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
872 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
873 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
874 if (!track) {\r
875 Printf("ERROR: Could not receive track %d", iTracks);\r
876 continue;\r
877 } \r
878 \r
879 Int_t label = TMath::Abs(track->GetLabel());\r
880 if(label > mcEvent->GetNumberOfTracks()) continue;\r
881 if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
882 \r
883 AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
884 if(!mcTrack) continue;\r
885\r
886 // take only TPC only tracks\r
887 track_TPC = new AliESDtrack();\r
888 if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
889 \r
890 //ESD track cuts\r
891 if(fESDtrackCuts) \r
892 if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
893 \r
894 // fill QA histograms\r
895 Float_t b[2];\r
896 Float_t bCov[3];\r
897 track_TPC->GetImpactParameters(b,bCov);\r
898 if (bCov[0]<=0 || bCov[2]<=0) {\r
899 AliDebug(1, "Estimated b resolution lower or equal zero!");\r
900 bCov[0]=0; bCov[2]=0;\r
901 }\r
902 \r
903 Int_t nClustersTPC = -1;\r
904 nClustersTPC = track_TPC->GetTPCNclsIter1(); // TPC standalone\r
905 //nClustersTPC = track->GetTPCclusters(0); // global track\r
906 Float_t chi2PerClusterTPC = -1;\r
907 if (nClustersTPC!=0) {\r
908 chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
909 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
910 }\r
911 \r
912 v_charge = track_TPC->Charge();\r
913 v_y = track_TPC->Y();\r
914 v_eta = track_TPC->Eta();\r
915 v_phi = track_TPC->Phi() * TMath::RadToDeg();\r
916 v_E = track_TPC->E();\r
917 v_pt = track_TPC->Pt();\r
918 track_TPC->PxPyPz(v_p);\r
919\r
920 fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
921 fHistDCA->Fill(b[1],b[0]);\r
922 fHistChi2->Fill(chi2PerClusterTPC);\r
923 fHistPt->Fill(v_pt);\r
924 fHistEta->Fill(v_eta);\r
925 fHistPhi->Fill(v_phi);\r
926\r
927 // fill charge vector\r
928 chargeVector[0]->push_back(v_charge);\r
929 chargeVector[1]->push_back(v_y);\r
930 chargeVector[2]->push_back(v_eta);\r
931 chargeVector[3]->push_back(v_phi);\r
932 chargeVector[4]->push_back(v_p[0]);\r
933 chargeVector[5]->push_back(v_p[1]);\r
934 chargeVector[6]->push_back(v_p[2]);\r
935 chargeVector[7]->push_back(v_pt);\r
936 chargeVector[8]->push_back(v_E);\r
937\r
938 if(fRunShuffling) {\r
939 chargeVectorShuffle[0]->push_back(v_charge);\r
940 chargeVectorShuffle[1]->push_back(v_y);\r
941 chargeVectorShuffle[2]->push_back(v_eta);\r
942 chargeVectorShuffle[3]->push_back(v_phi);\r
943 chargeVectorShuffle[4]->push_back(v_p[0]);\r
944 chargeVectorShuffle[5]->push_back(v_p[1]);\r
945 chargeVectorShuffle[6]->push_back(v_p[2]);\r
946 chargeVectorShuffle[7]->push_back(v_pt);\r
947 chargeVectorShuffle[8]->push_back(v_E);\r
948 }\r
949 \r
950 delete track_TPC;\r
951 gNumberOfAcceptedTracks += 1;\r
952 \r
953 } //track loop\r
954 }//Vz cut\r
955 }//Vy cut\r
956 }//Vx cut\r
957 }//proper vertex resolution\r
958 }//proper number of contributors\r
959 }//vertex object valid\r
960 }//triggered event \r
961 }//MC-ESD analysis\r
962\r
963 //MC analysis\r
964 else if(gAnalysisLevel == "MC") {\r
965 AliMCEvent* mcEvent = MCEvent(); \r
966 if (!mcEvent) {\r
967 Printf("ERROR: mcEvent not available");\r
968 return;\r
969 }\r
970 fHistEventStats->Fill(1); //total events\r
971 fHistEventStats->Fill(2); //offline trigger\r
972\r
973 Double_t gImpactParameter = 0.;\r
974 if(fUseCentrality) {\r
975 //Get the MC header\r
976 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
977 if (headerH) {\r
978 //Printf("=====================================================");\r
979 //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());\r
980 //Printf("Impact parameter: %lf",headerH->ImpactParameter());\r
981 //Printf("=====================================================");\r
982 gReactionPlane = headerH->ReactionPlaneAngle();\r
983 gImpactParameter = headerH->ImpactParameter();\r
984 fCentrality = gImpactParameter;\r
985 }\r
986 fCentrality = gImpactParameter;\r
f2583ca7 987 fHistEventPlane->Fill(gReactionPlane*TMath::RadToDeg());\r
0879e280 988\r
989 // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)\r
990 if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))\r
991 return;\r
992 }\r
993 \r
994 AliGenEventHeader *header = mcEvent->GenEventHeader();\r
995 if(!header) return;\r
996 \r
997 TArrayF gVertexArray;\r
998 header->PrimaryVertex(gVertexArray);\r
999 //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
1000 //gVertexArray.At(0),\r
1001 //gVertexArray.At(1),\r
1002 //gVertexArray.At(2));\r
1003 fHistEventStats->Fill(3); //events with a proper vertex\r
1004 if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
1005 if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
1006 if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
1007 fHistEventStats->Fill(4); //analayzed events\r
1008 fHistVx->Fill(gVertexArray.At(0));\r
1009 fHistVy->Fill(gVertexArray.At(1));\r
1010 fHistVz->Fill(gVertexArray.At(2));\r
1011 \r
1012 Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
1013 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
1014 AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
1015 if (!track) {\r
1016 Printf("ERROR: Could not receive particle %d", iTracks);\r
1017 continue;\r
1018 }\r
1019 \r
1020 //exclude non stable particles\r
1021 if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;\r
1022\r
1023 v_eta = track->Eta();\r
1024 v_pt = track->Pt();\r
1025 \r
1026 if( v_pt < fPtMin || v_pt > fPtMax) \r
1027 continue;\r
1028 if( v_eta < fEtaMin || v_eta > fEtaMax) \r
1029 continue;\r
1030 \r
1031 //analyze one set of particles\r
1032 if(fUseMCPdgCode) {\r
1033 TParticle *particle = track->Particle();\r
1034 if(!particle) continue;\r
1035 \r
1036 Int_t gPdgCode = particle->GetPdgCode();\r
1037 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1038 continue;\r
1039 }\r
1040 \r
1041 //Use the acceptance parameterization\r
1042 if(fAcceptanceParameterization) {\r
1043 Double_t gRandomNumber = gRandom->Rndm();\r
1044 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1045 continue;\r
1046 }\r
1047 \r
1048 //Exclude resonances\r
1049 if(fExcludeResonancesInMC) {\r
1050 TParticle *particle = track->Particle();\r
1051 if(!particle) continue;\r
1052 \r
1053 Bool_t kExcludeParticle = kFALSE;\r
1054 Int_t gMotherIndex = particle->GetFirstMother();\r
1055 if(gMotherIndex != -1) {\r
1056 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1057 if(motherTrack) {\r
1058 TParticle *motherParticle = motherTrack->Particle();\r
1059 if(motherParticle) {\r
1060 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1061 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1062 if(pdgCodeOfMother == 113) {\r
1063 kExcludeParticle = kTRUE;\r
1064 }\r
1065 }\r
1066 }\r
1067 }\r
1068 \r
1069 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1070 if(kExcludeParticle) continue;\r
1071 }\r
1072\r
1073 v_charge = track->Charge();\r
1074 v_y = track->Y();\r
1075 v_phi = track->Phi();\r
1076 v_E = track->E();\r
1077 track->PxPyPz(v_p);\r
1078 //Printf("phi (before): %lf",v_phi);\r
1079\r
1080 fHistPt->Fill(v_pt);\r
1081 fHistEta->Fill(v_eta);\r
1082 fHistPhi->Fill(v_phi);\r
1083\r
1084 //Flow after burner\r
1085 if(fUseFlowAfterBurner) {\r
1086 Double_t precisionPhi = 0.001;\r
1087 Int_t maxNumberOfIterations = 100;\r
1088\r
1089 Double_t phi0 = v_phi;\r
1090 Double_t gV2 = fDifferentialV2->Eval(v_pt);\r
1091\r
1092 for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1093 Double_t phiprev = v_phi;\r
1094 Double_t fl = v_phi - phi0 + gV2*TMath::Sin(2.*(v_phi - gReactionPlane));\r
1095 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(v_phi - gReactionPlane)); \r
1096 v_phi -= fl/fp;\r
1097 if (TMath::AreEqualAbs(phiprev,v_phi,precisionPhi)) break;\r
1098 }\r
1099 //Printf("phi (after): %lf\n",v_phi);\r
1100 Double_t v_DeltaphiBefore = phi0 - gReactionPlane;\r
1101 if(v_DeltaphiBefore < 0) v_DeltaphiBefore += 2*TMath::Pi();\r
1102 fHistPhiBefore->Fill(v_DeltaphiBefore);\r
1103\r
1104 Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
1105 if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
1106 fHistPhiAfter->Fill(v_DeltaphiAfter);\r
1107 }\r
1108 \r
1109 v_phi *= TMath::RadToDeg();\r
1110\r
1111 // fill charge vector\r
1112 chargeVector[0]->push_back(v_charge);\r
1113 chargeVector[1]->push_back(v_y);\r
1114 chargeVector[2]->push_back(v_eta);\r
1115 chargeVector[3]->push_back(v_phi);\r
1116 chargeVector[4]->push_back(v_p[0]);\r
1117 chargeVector[5]->push_back(v_p[1]);\r
1118 chargeVector[6]->push_back(v_p[2]);\r
1119 chargeVector[7]->push_back(v_pt);\r
1120 chargeVector[8]->push_back(v_E);\r
1121 \r
1122 if(fRunShuffling) {\r
1123 chargeVectorShuffle[0]->push_back(v_charge);\r
1124 chargeVectorShuffle[1]->push_back(v_y);\r
1125 chargeVectorShuffle[2]->push_back(v_eta);\r
1126 chargeVectorShuffle[3]->push_back(v_phi);\r
1127 chargeVectorShuffle[4]->push_back(v_p[0]);\r
1128 chargeVectorShuffle[5]->push_back(v_p[1]);\r
1129 chargeVectorShuffle[6]->push_back(v_p[2]);\r
1130 chargeVectorShuffle[7]->push_back(v_pt);\r
1131 chargeVectorShuffle[8]->push_back(v_E);\r
1132 }\r
1133 gNumberOfAcceptedTracks += 1;\r
1134 \r
1135 } //track loop\r
1136 }//Vz cut\r
1137 }//Vy cut\r
1138 }//Vx cut\r
1139 }//MC analysis\r
1140 \r
1141 //multiplicity cut (used in pp)\r
1142 if(fUseMultiplicity) {\r
1143 if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
1144 return;\r
1145 }\r
1146 fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks);\r
1147 \r
1148 // calculate balance function\r
1149 if(fUseMultiplicity) \r
f2583ca7 1150 fBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane*TMath::RadToDeg(),chargeVector);\r
0879e280 1151 else \r
f2583ca7 1152 fBalance->CalculateBalance(fCentrality,gReactionPlane*TMath::RadToDeg(),chargeVector);\r
0879e280 1153\r
1154 if(fRunShuffling) {\r
1155 // shuffle charges\r
1156 random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
1157 if(fUseMultiplicity) \r
1158 fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane,chargeVectorShuffle);\r
1159 else \r
1160 fShuffledBalance->CalculateBalance(fCentrality,gReactionPlane,chargeVectorShuffle);\r
1161 }\r
1162} \r
1163\r
1164//________________________________________________________________________\r
1165void AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1166 //Printf("END BF");\r
1167\r
1168 if (!fBalance) {\r
1169 Printf("ERROR: fBalance not available");\r
1170 return;\r
1171 } \r
1172 if(fRunShuffling) {\r
1173 if (!fShuffledBalance) {\r
1174 Printf("ERROR: fShuffledBalance not available");\r
1175 return;\r
1176 }\r
1177 }\r
1178\r
1179}\r
1180\r
1181//________________________________________________________________________\r
1182void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1183 // Draw result to the screen\r
1184 // Called once at the end of the query\r
1185\r
1186 // not implemented ...\r
1187\r
1188}\r