]> 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
Fix bug in trigger class determination when running on MC productions (AliTimeStamp...
[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
909 // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality)  \r
910 \r
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
927 \r
928 if (vCharge > 0) {\r
929 correction = fHistMatrixCorrectionPlus[gCentralityInt-1]->GetBinContent(fHistMatrixCorrectionPlus[gCentralityInt-1]->GetBin(binEta, binPt, binPhi));\r
930 }\r
931 if (vCharge < 0) {\r
932 //correction = fHistMatrixCorrectionMinus->GetBinContent(fHistMatrixCorrectionMinus->GetBin(dimBin));\r
933 correction = fHistMatrixCorrectionMinus[gCentralityInt-1]->GetBinContent(fHistMatrixCorrectionMinus[gCentralityInt-1]->GetBin(binEta, binPt, binPhi));\r
934 }\r
935 if (correction == 0.) { \r
936 AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
937 correction = 1.; \r
938 } \r
939\r
940 return correction;\r
941}\r
942//========================correction=============================//\r
943\r
f06d59b3 944//________________________________________________________________________\r
35aff0f3 945TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
f06d59b3 946 // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
947 // Fills QA histograms\r
0879e280 948\r
f06d59b3 949 TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
0879e280 950\r
f06d59b3 951 //output TObjArray holding all good tracks\r
952 TObjArray* tracksAccepted = new TObjArray;\r
953 tracksAccepted->SetOwner(kTRUE);\r
0879e280 954\r
35aff0f3 955 Short_t vCharge;\r
9fd4b54e 956 Double_t vEta;\r
957 Double_t vY;\r
958 Double_t vPhi;\r
959 Double_t vPt;\r
0879e280 960\r
f06d59b3 961\r
962 if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
f06d59b3 963 // Loop over tracks in event\r
964 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
965 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
966 if (!aodTrack) {\r
967 AliError(Form("Could not receive track %d", iTracks));\r
968 continue;\r
969 }\r
970 \r
971 // AOD track cuts\r
972 \r
973 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
974 // take only TPC only tracks \r
c1847400 975 //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
976 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
977 fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
978 }\r
f06d59b3 979 if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
980 \r
9fd4b54e 981 vCharge = aodTrack->Charge();\r
982 vEta = aodTrack->Eta();\r
983 vY = aodTrack->Y();\r
c8ad9635 984 vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r
9fd4b54e 985 vPt = aodTrack->Pt();\r
f06d59b3 986 \r
9fd4b54e 987 Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r
988 Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r
f06d59b3 989 \r
990 \r
991 // Kinematics cuts from ESD track cuts\r
9fd4b54e 992 if( vPt < fPtMin || vPt > fPtMax) continue;\r
993 if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
f06d59b3 994 \r
995 // Extra DCA cuts (for systematic studies [!= -1])\r
996 if( fDCAxyCut != -1 && fDCAzCut != -1){\r
9fd4b54e 997 if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
f06d59b3 998 continue; // 2D cut\r
999 }\r
1000 }\r
1001 \r
1002 // Extra TPC cuts (for systematic studies [!= -1])\r
1003 if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1004 continue;\r
0879e280 1005 }\r
f06d59b3 1006 if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1007 continue;\r
1008 }\r
1009 \r
1010 // fill QA histograms\r
1011 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
9fd4b54e 1012 fHistDCA->Fill(dcaZ,dcaXY);\r
35aff0f3 1013 fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1014 fHistPt->Fill(vPt,gCentrality);\r
1015 fHistEta->Fill(vEta,gCentrality);\r
1016 fHistRapidity->Fill(vY,gCentrality);\r
1017 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1018 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1019 fHistPhi->Fill(vPhi,gCentrality);\r
1020 \r
1021 //=======================================correction\r
1022 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
f06d59b3 1023 \r
1024 // add the track to the TObjArray\r
35aff0f3 1025 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
f06d59b3 1026 }//track loop\r
1027 }// AOD analysis\r
0879e280 1028\r
f06d59b3 1029\r
1030 else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1031\r
9fd4b54e 1032 AliESDtrack *trackTPC = NULL;\r
f06d59b3 1033 AliMCParticle *mcTrack = NULL;\r
1034\r
1035 AliMCEvent* mcEvent = NULL;\r
0879e280 1036 \r
f06d59b3 1037 // for MC ESDs use also MC information\r
1038 if(gAnalysisLevel == "MCESD"){\r
1039 mcEvent = MCEvent(); \r
1040 if (!mcEvent) {\r
1041 AliError("mcEvent not available");\r
1042 return tracksAccepted;\r
1043 }\r
1044 }\r
0879e280 1045 \r
f06d59b3 1046 // Loop over tracks in event\r
1047 for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1048 AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1049 if (!track) {\r
1050 AliError(Form("Could not receive track %d", iTracks));\r
1051 continue;\r
1052 } \r
1053\r
1054 // for MC ESDs use also MC information --> MC track not used further???\r
1055 if(gAnalysisLevel == "MCESD"){\r
1056 Int_t label = TMath::Abs(track->GetLabel());\r
1057 if(label > mcEvent->GetNumberOfTracks()) continue;\r
1058 if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1059 \r
1060 mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1061 if(!mcTrack) continue;\r
1062 }\r
0879e280 1063\r
f06d59b3 1064 // take only TPC only tracks\r
9fd4b54e 1065 trackTPC = new AliESDtrack();\r
1066 if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
f06d59b3 1067 \r
1068 //ESD track cuts\r
1069 if(fESDtrackCuts) \r
9fd4b54e 1070 if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
f06d59b3 1071 \r
1072 // fill QA histograms\r
1073 Float_t b[2];\r
1074 Float_t bCov[3];\r
9fd4b54e 1075 trackTPC->GetImpactParameters(b,bCov);\r
f06d59b3 1076 if (bCov[0]<=0 || bCov[2]<=0) {\r
1077 AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1078 bCov[0]=0; bCov[2]=0;\r
1079 }\r
1080 \r
1081 Int_t nClustersTPC = -1;\r
9fd4b54e 1082 nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone\r
f06d59b3 1083 //nClustersTPC = track->GetTPCclusters(0); // global track\r
1084 Float_t chi2PerClusterTPC = -1;\r
1085 if (nClustersTPC!=0) {\r
9fd4b54e 1086 chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r
f06d59b3 1087 //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r
1088 }\r
1089 \r
1090 //===========================PID===============================// \r
1091 if(fUsePID) {\r
1092 Double_t prob[AliPID::kSPECIES]={0.};\r
1093 Double_t probTPC[AliPID::kSPECIES]={0.};\r
1094 Double_t probTOF[AliPID::kSPECIES]={0.};\r
1095 Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1096 \r
1097 Double_t nSigma = 0.;\r
1098 UInt_t detUsedTPC = 0;\r
1099 UInt_t detUsedTOF = 0;\r
1100 UInt_t detUsedTPCTOF = 0;\r
1101 \r
1102 //Decide what detector configuration we want to use\r
1103 switch(fPidDetectorConfig) {\r
1104 case kTPCpid:\r
1105 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1106 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1107 detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1108 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1109 prob[iSpecies] = probTPC[iSpecies];\r
1110 break;\r
1111 case kTOFpid:\r
1112 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1113 nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1114 detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1115 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1116 prob[iSpecies] = probTOF[iSpecies];\r
1117 break;\r
1118 case kTPCTOF:\r
1119 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1120 detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1121 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1122 prob[iSpecies] = probTPCTOF[iSpecies];\r
1123 break;\r
1124 default:\r
1125 break;\r
1126 }//end switch: define detector mask\r
1127 \r
1128 //Filling the PID QA\r
1129 Double_t tofTime = -999., length = 999., tof = -999.;\r
1130 Double_t c = TMath::C()*1.E-9;// m/ns\r
1131 Double_t beta = -999.;\r
1132 Double_t nSigmaTOFForParticleOfInterest = -999.;\r
1133 if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1134 (track->IsOn(AliESDtrack::kTIME)) ) { \r
1135 tofTime = track->GetTOFsignal();//in ps\r
1136 length = track->GetIntegratedLength();\r
1137 tof = tofTime*1E-3; // ns \r
1138 \r
1139 if (tof <= 0) {\r
1140 //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1141 continue;\r
1142 }\r
1143 if (length <= 0){\r
1144 //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1145 continue;\r
1146 }\r
1147 \r
1148 length = length*0.01; // in meters\r
1149 tof = tof*c;\r
1150 beta = length/tof;\r
1151 \r
1152 nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1153 fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1154 fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1155 fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1156 }//TOF signal \r
1157 \r
1158 \r
1159 Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1160 fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1161 fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1162 fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1163 fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1164 //end of QA-before pid\r
1165 \r
1166 if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1167 //Make the decision based on the n-sigma\r
1168 if(fUsePIDnSigma) {\r
1169 if(nSigma > fPIDNSigma) continue;}\r
1170 \r
1171 //Make the decision based on the bayesian\r
1172 else if(fUsePIDPropabilities) {\r
1173 if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1174 if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r
1175 }\r
1176 \r
1177 //Fill QA after the PID\r
1178 fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1179 fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1180 fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1181 \r
1182 fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1183 fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1184 fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1185 fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1186 }\r
1187 }\r
1188 //===========================PID===============================//\r
9fd4b54e 1189 vCharge = trackTPC->Charge();\r
1190 vY = trackTPC->Y();\r
1191 vEta = trackTPC->Eta();\r
c8ad9635 1192 vPhi = trackTPC->Phi();// * TMath::RadToDeg();\r
9fd4b54e 1193 vPt = trackTPC->Pt();\r
1194 fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
f06d59b3 1195 fHistDCA->Fill(b[1],b[0]);\r
35aff0f3 1196 fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1197 fHistPt->Fill(vPt,gCentrality);\r
1198 fHistEta->Fill(vEta,gCentrality);\r
1199 fHistPhi->Fill(vPhi,gCentrality);\r
1200 fHistRapidity->Fill(vY,gCentrality);\r
1201 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1202 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
f06d59b3 1203 \r
35aff0f3 1204 //=======================================correction\r
1205 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
1206\r
f06d59b3 1207 // add the track to the TObjArray\r
35aff0f3 1208 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
244c1a3d 1209\r
9fd4b54e 1210 delete trackTPC;\r
f06d59b3 1211 }//track loop\r
1212 }// ESD analysis\r
244c1a3d 1213\r
73609a48 1214 else if(gAnalysisLevel == "MC"){\r
f06d59b3 1215 \r
1216 // Loop over tracks in event\r
1217 for (Int_t iTracks = 0; iTracks < dynamic_cast<AliMCEvent*>(event)->GetNumberOfPrimaries(); iTracks++) {\r
1218 AliMCParticle* track = dynamic_cast<AliMCParticle *>(event->GetTrack(iTracks));\r
1219 if (!track) {\r
1220 AliError(Form("Could not receive particle %d", iTracks));\r
1221 continue;\r
1222 }\r
0879e280 1223 \r
f06d59b3 1224 //exclude non stable particles\r
1225 if(!(dynamic_cast<AliMCEvent*>(event)->IsPhysicalPrimary(iTracks))) continue;\r
0879e280 1226\r
9fd4b54e 1227 vEta = track->Eta();\r
1228 vPt = track->Pt();\r
1229 vY = track->Y();\r
f06d59b3 1230 \r
9fd4b54e 1231 if( vPt < fPtMin || vPt > fPtMax) \r
f06d59b3 1232 continue;\r
1233 if (!fUsePID) {\r
9fd4b54e 1234 if( vEta < fEtaMin || vEta > fEtaMax) continue;\r
f06d59b3 1235 }\r
1236 else if (fUsePID){\r
9fd4b54e 1237 if( vY < fEtaMin || vY > fEtaMax) continue;\r
f06d59b3 1238 }\r
1239 \r
1240 //analyze one set of particles\r
1241 if(fUseMCPdgCode) {\r
1242 TParticle *particle = track->Particle();\r
1243 if(!particle) continue;\r
1244 \r
1245 Int_t gPdgCode = particle->GetPdgCode();\r
1246 if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1247 continue;\r
1248 }\r
1249 \r
1250 //Use the acceptance parameterization\r
1251 if(fAcceptanceParameterization) {\r
1252 Double_t gRandomNumber = gRandom->Rndm();\r
1253 if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1254 continue;\r
1255 }\r
1256 \r
1257 //Exclude resonances\r
1258 if(fExcludeResonancesInMC) {\r
1259 TParticle *particle = track->Particle();\r
1260 if(!particle) continue;\r
1261 \r
1262 Bool_t kExcludeParticle = kFALSE;\r
1263 Int_t gMotherIndex = particle->GetFirstMother();\r
1264 if(gMotherIndex != -1) {\r
1265 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1266 if(motherTrack) {\r
1267 TParticle *motherParticle = motherTrack->Particle();\r
1268 if(motherParticle) {\r
1269 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1270 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1271 if(pdgCodeOfMother == 113) {\r
1272 kExcludeParticle = kTRUE;\r
0879e280 1273 }\r
0879e280 1274 }\r
f06d59b3 1275 }\r
1276 }\r
1277 \r
1278 //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1279 if(kExcludeParticle) continue;\r
1280 }\r
1281 \r
9fd4b54e 1282 vCharge = track->Charge();\r
1283 vPhi = track->Phi();\r
1284 //Printf("phi (before): %lf",vPhi);\r
f06d59b3 1285 \r
35aff0f3 1286 fHistPt->Fill(vPt,gCentrality);\r
1287 fHistEta->Fill(vEta,gCentrality);\r
1288 fHistPhi->Fill(vPhi,gCentrality);\r
1289 //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1290 fHistRapidity->Fill(vY,gCentrality);\r
1291 //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1292 //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1293 if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1294 else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
f06d59b3 1295 \r
1296 //Flow after burner\r
1297 if(fUseFlowAfterBurner) {\r
1298 Double_t precisionPhi = 0.001;\r
1299 Int_t maxNumberOfIterations = 100;\r
1300 \r
9fd4b54e 1301 Double_t phi0 = vPhi;\r
1302 Double_t gV2 = fDifferentialV2->Eval(vPt);\r
f06d59b3 1303 \r
1304 for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
9fd4b54e 1305 Double_t phiprev = vPhi;\r
648f1a5a 1306 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
1307 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
9fd4b54e 1308 vPhi -= fl/fp;\r
1309 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
f06d59b3 1310 }\r
9fd4b54e 1311 //Printf("phi (after): %lf\n",vPhi);\r
648f1a5a 1312 Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
9fd4b54e 1313 if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
35aff0f3 1314 fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
f06d59b3 1315 \r
648f1a5a 1316 Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
9fd4b54e 1317 if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
35aff0f3 1318 fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
1319 \r
f06d59b3 1320 }\r
1321 \r
c8ad9635 1322 //vPhi *= TMath::RadToDeg();\r
35aff0f3 1323\r
1324 //=======================================correction\r
1325 Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r
1326 \r
1327 tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
f06d59b3 1328 \r
1329 } //track loop\r
1330 }//MC\r
0879e280 1331 \r
f06d59b3 1332 return tracksAccepted; \r
1333}\r
1334//________________________________________________________________________\r
35aff0f3 1335TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
f06d59b3 1336 // Clones TObjArray and returns it with tracks after shuffling the charges\r
0879e280 1337\r
f06d59b3 1338 TObjArray* tracksShuffled = new TObjArray;\r
1339 tracksShuffled->SetOwner(kTRUE);\r
1340\r
1341 vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks \r
1342\r
1343 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
1344 {\r
1345 AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1346 chargeVector->push_back(track->Charge());\r
1347 } \r
1348 \r
1349 random_shuffle(chargeVector->begin(), chargeVector->end());\r
1350 \r
1351 for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
1352 AliVParticle* track = (AliVParticle*) tracks->At(i);\r
35aff0f3 1353 //==============================correction\r
1354 Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
1355 tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
0879e280 1356 }\r
f06d59b3 1357\r
1358 delete chargeVector;\r
1359 \r
1360 return tracksShuffled;\r
1361}\r
1362\r
0879e280 1363\r
1364//________________________________________________________________________\r
1365void AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1366 //Printf("END BF");\r
1367\r
1368 if (!fBalance) {\r
f06d59b3 1369 AliError("fBalance not available");\r
0879e280 1370 return;\r
1371 } \r
1372 if(fRunShuffling) {\r
1373 if (!fShuffledBalance) {\r
f06d59b3 1374 AliError("fShuffledBalance not available");\r
0879e280 1375 return;\r
1376 }\r
1377 }\r
1378\r
1379}\r
1380\r
1381//________________________________________________________________________\r
1382void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1383 // Draw result to the screen\r
1384 // Called once at the end of the query\r
1385\r
1386 // not implemented ...\r
1387\r
1388}\r