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