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