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