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