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