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