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