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