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