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