]>
Commit | Line | Data |
---|---|---|
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 | |
1cb2a06e | 24 | #include "AliAODMCParticle.h" \r |
951cf00e | 25 | #include "AliCollisionGeometry.h"\r |
0879e280 | 26 | #include "AliGenEventHeader.h"\r |
0879e280 | 27 | #include "AliMCEventHandler.h"\r |
28 | #include "AliMCEvent.h"\r | |
29 | #include "AliStack.h"\r | |
30 | #include "AliESDtrackCuts.h"\r | |
31 | #include "AliEventplane.h"\r | |
f06d59b3 | 32 | #include "AliTHn.h" \r |
6d6eb2df | 33 | #include "AliLog.h"\r |
7556a0f0 | 34 | #include "AliAnalysisUtils.h"\r |
f06d59b3 | 35 | \r |
36 | #include "AliEventPoolManager.h" \r | |
0879e280 | 37 | \r |
38 | #include "AliPID.h" \r | |
39 | #include "AliPIDResponse.h" \r | |
40 | #include "AliPIDCombined.h" \r | |
41 | \r | |
42 | #include "AliAnalysisTaskBFPsi.h"\r | |
43 | #include "AliBalancePsi.h"\r | |
f06d59b3 | 44 | #include "AliAnalysisTaskTriggeredBF.h"\r |
0879e280 | 45 | \r |
46 | \r | |
47 | // Analysis task for the BF vs Psi code\r | |
48 | // Authors: Panos.Christakoglou@nikhef.nl\r | |
49 | \r | |
c64cb1f6 | 50 | using std::cout;\r |
51 | using std::endl;\r | |
52 | \r | |
0879e280 | 53 | ClassImp(AliAnalysisTaskBFPsi)\r |
54 | \r | |
55 | //________________________________________________________________________\r | |
56 | AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r | |
6f719a5b | 57 | : AliAnalysisTaskSE(name),\r |
a3ed6fc1 | 58 | fDebugLevel(kFALSE),\r |
6f719a5b | 59 | fArrayMC(0), //+++++++++++++\r |
0879e280 | 60 | fBalance(0),\r |
61 | fRunShuffling(kFALSE),\r | |
62 | fShuffledBalance(0),\r | |
f06d59b3 | 63 | fRunMixing(kFALSE),\r |
8cdfbf52 | 64 | fRunMixingEventPlane(kFALSE),\r |
f06d59b3 | 65 | fMixingTracks(50000),\r |
66 | fMixedBalance(0),\r | |
67 | fPoolMgr(0),\r | |
0879e280 | 68 | fList(0),\r |
69 | fListBF(0),\r | |
70 | fListBFS(0),\r | |
f06d59b3 | 71 | fListBFM(0),\r |
0879e280 | 72 | fHistListPIDQA(0),\r |
73 | fHistEventStats(0),\r | |
74 | fHistCentStats(0),\r | |
6f719a5b | 75 | fHistCentStatsUsed(0),\r |
0879e280 | 76 | fHistTriggerStats(0),\r |
77 | fHistTrackStats(0),\r | |
78 | fHistVx(0),\r | |
79 | fHistVy(0),\r | |
80 | fHistVz(0),\r | |
4d528efb | 81 | fHistTPCvsVZEROMultiplicity(0),\r |
82 | fHistVZEROSignal(0),\r | |
0879e280 | 83 | fHistEventPlane(0),\r |
84 | fHistClus(0),\r | |
85 | fHistDCA(0),\r | |
86 | fHistChi2(0),\r | |
87 | fHistPt(0),\r | |
88 | fHistEta(0),\r | |
244c1a3d | 89 | fHistRapidity(0),\r |
0879e280 | 90 | fHistPhi(0),\r |
f2e8af26 | 91 | fHistEtaPhiPos(0), \r |
92 | fHistEtaPhiNeg(0), \r | |
0879e280 | 93 | fHistPhiBefore(0),\r |
94 | fHistPhiAfter(0),\r | |
244c1a3d | 95 | fHistPhiPos(0),\r |
96 | fHistPhiNeg(0),\r | |
0879e280 | 97 | fHistV0M(0),\r |
98 | fHistRefTracks(0),\r | |
99 | fHistdEdxVsPTPCbeforePID(NULL),\r | |
100 | fHistBetavsPTOFbeforePID(NULL), \r | |
101 | fHistProbTPCvsPtbeforePID(NULL), \r | |
102 | fHistProbTOFvsPtbeforePID(NULL), \r | |
103 | fHistProbTPCTOFvsPtbeforePID(NULL),\r | |
104 | fHistNSigmaTPCvsPtbeforePID(NULL), \r | |
105 | fHistNSigmaTOFvsPtbeforePID(NULL), \r | |
90bec8c5 | 106 | fHistBetaVsdEdXbeforePID(NULL), //+++++++ \r |
107 | fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++\r | |
0879e280 | 108 | fHistdEdxVsPTPCafterPID(NULL),\r |
109 | fHistBetavsPTOFafterPID(NULL), \r | |
110 | fHistProbTPCvsPtafterPID(NULL), \r | |
111 | fHistProbTOFvsPtafterPID(NULL), \r | |
112 | fHistProbTPCTOFvsPtafterPID(NULL),\r | |
113 | fHistNSigmaTPCvsPtafterPID(NULL), \r | |
114 | fHistNSigmaTOFvsPtafterPID(NULL), \r | |
90bec8c5 | 115 | fHistBetaVsdEdXafterPID(NULL), //+++++++ \r |
116 | fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++\r | |
6fa567bd | 117 | fCentralityArrayBinsForCorrections(kCENTRALITY),\r |
0879e280 | 118 | fPIDResponse(0x0),\r |
119 | fPIDCombined(0x0),\r | |
120 | fParticleOfInterest(kPion),\r | |
121 | fPidDetectorConfig(kTPCTOF),\r | |
122 | fUsePID(kFALSE),\r | |
123 | fUsePIDnSigma(kTRUE),\r | |
124 | fUsePIDPropabilities(kFALSE), \r | |
125 | fPIDNSigma(3.),\r | |
126 | fMinAcceptedPIDProbability(0.8),\r | |
c443da34 | 127 | fElectronRejection(kFALSE),\r |
8b3da35d | 128 | fElectronOnlyRejection(kFALSE),\r |
c443da34 | 129 | fElectronRejectionNSigma(-1.),\r |
8b3da35d | 130 | fElectronRejectionMinPt(0.),\r |
131 | fElectronRejectionMaxPt(1000.),\r | |
0879e280 | 132 | fESDtrackCuts(0),\r |
133 | fCentralityEstimator("V0M"),\r | |
134 | fUseCentrality(kFALSE),\r | |
135 | fCentralityPercentileMin(0.), \r | |
136 | fCentralityPercentileMax(5.),\r | |
137 | fImpactParameterMin(0.),\r | |
138 | fImpactParameterMax(20.),\r | |
c16b246d | 139 | fMultiplicityEstimator("V0A"),\r |
0879e280 | 140 | fUseMultiplicity(kFALSE),\r |
141 | fNumberOfAcceptedTracksMin(0),\r | |
142 | fNumberOfAcceptedTracksMax(10000),\r | |
143 | fHistNumberOfAcceptedTracks(0),\r | |
6ec16c49 | 144 | fHistMultiplicity(0),\r |
0879e280 | 145 | fUseOfflineTrigger(kFALSE),\r |
7556a0f0 | 146 | fCheckFirstEventInChunk(kFALSE),\r |
147 | fCheckPileUp(kFALSE),\r | |
d9eb282c | 148 | fUseMCforKinematics(kFALSE),\r |
0879e280 | 149 | fVxMax(0.3),\r |
150 | fVyMax(0.3),\r | |
151 | fVzMax(10.),\r | |
4d528efb | 152 | fnAODtrackCutBit(128),\r |
0879e280 | 153 | fPtMin(0.3),\r |
154 | fPtMax(1.5),\r | |
ce8126c8 | 155 | fPtMinForCorrections(0.3),//=================================correction\r |
156 | fPtMaxForCorrections(1.5),//=================================correction\r | |
35aff0f3 | 157 | fPtBinForCorrections(36), //=================================correction\r |
0879e280 | 158 | fEtaMin(-0.8),\r |
159 | fEtaMax(-0.8),\r | |
ce8126c8 | 160 | fEtaMinForCorrections(-0.8),//=================================correction\r |
161 | fEtaMaxForCorrections(0.8),//=================================correction\r | |
162 | fEtaBinForCorrections(16), //=================================correction\r | |
163 | fPhiMin(0.),\r | |
164 | fPhiMax(360.),\r | |
35aff0f3 | 165 | fPhiMinForCorrections(0.),//=================================correction\r |
166 | fPhiMaxForCorrections(360.),//=================================correction\r | |
167 | fPhiBinForCorrections(100), //=================================correction\r | |
0879e280 | 168 | fDCAxyCut(-1),\r |
169 | fDCAzCut(-1),\r | |
170 | fTPCchi2Cut(-1),\r | |
171 | fNClustersTPCCut(-1),\r | |
172 | fAcceptanceParameterization(0),\r | |
173 | fDifferentialV2(0),\r | |
174 | fUseFlowAfterBurner(kFALSE),\r | |
175 | fExcludeResonancesInMC(kFALSE),\r | |
8b3da35d | 176 | fExcludeElectronsInMC(kFALSE),\r |
0879e280 | 177 | fUseMCPdgCode(kFALSE),\r |
f0fb4ac1 | 178 | fPDGCodeToBeAnalyzed(-1),\r |
ea9867da | 179 | fEventClass("EventPlane"), \r |
91e5e210 | 180 | fCustomBinning(""),\r |
181 | fHistVZEROAGainEqualizationMap(0),\r | |
182 | fHistVZEROCGainEqualizationMap(0),\r | |
183 | fHistVZEROChannelGainEqualizationMap(0) {\r | |
0879e280 | 184 | // Constructor\r |
185 | // Define input and output slots here\r | |
186 | // Input slot #0 works with a TChain\r | |
ce8126c8 | 187 | \r |
188 | //======================================================correction\r | |
9006a70a | 189 | for (Int_t i=0; i<kCENTRALITY; i++){\r |
1cb2a06e | 190 | fHistCorrectionPlus[i] = NULL; \r |
191 | fHistCorrectionMinus[i] = NULL; \r | |
6fa567bd | 192 | fCentralityArrayForCorrections[i] = -1.;\r |
35aff0f3 | 193 | }\r |
ce8126c8 | 194 | //=====================================================correction\r |
195 | \r | |
0879e280 | 196 | DefineInput(0, TChain::Class());\r |
197 | // Output slot #0 writes into a TH1 container\r | |
198 | DefineOutput(1, TList::Class());\r | |
199 | DefineOutput(2, TList::Class());\r | |
200 | DefineOutput(3, TList::Class());\r | |
201 | DefineOutput(4, TList::Class());\r | |
f06d59b3 | 202 | DefineOutput(5, TList::Class());\r |
0879e280 | 203 | }\r |
204 | \r | |
205 | //________________________________________________________________________\r | |
206 | AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r | |
207 | \r | |
208 | // delete fBalance; \r | |
209 | // delete fShuffledBalance; \r | |
210 | // delete fList;\r | |
211 | // delete fListBF; \r | |
212 | // delete fListBFS;\r | |
213 | \r | |
214 | // delete fHistEventStats; \r | |
215 | // delete fHistTrackStats; \r | |
216 | // delete fHistVx; \r | |
217 | // delete fHistVy; \r | |
218 | // delete fHistVz; \r | |
219 | \r | |
220 | // delete fHistClus;\r | |
221 | // delete fHistDCA;\r | |
222 | // delete fHistChi2;\r | |
223 | // delete fHistPt;\r | |
224 | // delete fHistEta;\r | |
225 | // delete fHistPhi;\r | |
f2e8af26 | 226 | // delete fHistEtaPhiPos; \r |
227 | // delete fHistEtaPhiNeg;\r | |
0879e280 | 228 | // delete fHistV0M;\r |
229 | }\r | |
230 | \r | |
231 | //________________________________________________________________________\r | |
232 | void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r | |
233 | // Create histograms\r | |
234 | // Called once\r | |
211b716d | 235 | \r |
236 | // global switch disabling the reference \r | |
237 | // (to avoid "Replacing existing TH1" if several wagons are created in train)\r | |
238 | Bool_t oldStatus = TH1::AddDirectoryStatus();\r | |
239 | TH1::AddDirectory(kFALSE);\r | |
240 | \r | |
0879e280 | 241 | if(!fBalance) {\r |
242 | fBalance = new AliBalancePsi();\r | |
243 | fBalance->SetAnalysisLevel("ESD");\r | |
f0fb4ac1 | 244 | fBalance->SetEventClass(fEventClass);\r |
0879e280 | 245 | //fBalance->SetNumberOfBins(-1,16);\r |
9afe3098 | 246 | //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r |
0879e280 | 247 | }\r |
248 | if(fRunShuffling) {\r | |
249 | if(!fShuffledBalance) {\r | |
250 | fShuffledBalance = new AliBalancePsi();\r | |
251 | fShuffledBalance->SetAnalysisLevel("ESD");\r | |
f0fb4ac1 | 252 | fShuffledBalance->SetEventClass(fEventClass);\r |
0879e280 | 253 | //fShuffledBalance->SetNumberOfBins(-1,16);\r |
9afe3098 | 254 | //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r |
0879e280 | 255 | }\r |
256 | }\r | |
f06d59b3 | 257 | if(fRunMixing) {\r |
258 | if(!fMixedBalance) {\r | |
259 | fMixedBalance = new AliBalancePsi();\r | |
260 | fMixedBalance->SetAnalysisLevel("ESD");\r | |
f0fb4ac1 | 261 | fMixedBalance->SetEventClass(fEventClass);\r |
f06d59b3 | 262 | }\r |
263 | }\r | |
0879e280 | 264 | \r |
265 | //QA list\r | |
266 | fList = new TList();\r | |
267 | fList->SetName("listQA");\r | |
268 | fList->SetOwner();\r | |
269 | \r | |
270 | //Balance Function list\r | |
271 | fListBF = new TList();\r | |
272 | fListBF->SetName("listBF");\r | |
273 | fListBF->SetOwner();\r | |
274 | \r | |
275 | if(fRunShuffling) {\r | |
276 | fListBFS = new TList();\r | |
277 | fListBFS->SetName("listBFShuffled");\r | |
278 | fListBFS->SetOwner();\r | |
279 | }\r | |
280 | \r | |
f06d59b3 | 281 | if(fRunMixing) {\r |
282 | fListBFM = new TList();\r | |
283 | fListBFM->SetName("listTriggeredBFMixed");\r | |
284 | fListBFM->SetOwner();\r | |
285 | }\r | |
286 | \r | |
0879e280 | 287 | //PID QA list\r |
c443da34 | 288 | if(fUsePID || fElectronRejection) {\r |
0879e280 | 289 | fHistListPIDQA = new TList();\r |
290 | fHistListPIDQA->SetName("listQAPID");\r | |
291 | fHistListPIDQA->SetOwner();\r | |
292 | }\r | |
293 | \r | |
294 | //Event stats.\r | |
7556a0f0 | 295 | TString gCutName[7] = {"Total","Offline trigger",\r |
296 | "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};\r | |
6b046844 | 297 | fHistEventStats = new TH2F("fHistEventStats",\r |
298 | "Event statistics;;Centrality percentile;N_{events}",\r | |
7556a0f0 | 299 | 7,0.5,7.5,220,-5,105);\r |
300 | for(Int_t i = 1; i <= 7; i++)\r | |
0879e280 | 301 | fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r |
302 | fList->Add(fHistEventStats);\r | |
303 | \r | |
2213e863 | 304 | TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r |
0879e280 | 305 | fHistCentStats = new TH2F("fHistCentStats",\r |
306 | "Centrality statistics;;Cent percentile",\r | |
2213e863 | 307 | 13,-0.5,12.5,220,-5,105);\r |
308 | for(Int_t i = 1; i <= 13; i++){\r | |
0879e280 | 309 | fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r |
90bec8c5 | 310 | //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r |
2213e863 | 311 | }\r |
0879e280 | 312 | fList->Add(fHistCentStats);\r |
313 | \r | |
90bec8c5 | 314 | fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105);\r |
315 | fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());\r | |
316 | fList->Add(fHistCentStatsUsed);\r | |
2213e863 | 317 | \r |
c1847400 | 318 | fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r |
0879e280 | 319 | fList->Add(fHistTriggerStats);\r |
320 | \r | |
c1847400 | 321 | fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);\r |
0879e280 | 322 | fList->Add(fHistTrackStats);\r |
323 | \r | |
6b046844 | 324 | fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r |
0879e280 | 325 | fList->Add(fHistNumberOfAcceptedTracks);\r |
326 | \r | |
526d35fd | 327 | fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5);\r |
6ec16c49 | 328 | fList->Add(fHistMultiplicity);\r |
329 | \r | |
0879e280 | 330 | // Vertex distributions\r |
331 | fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r | |
332 | fList->Add(fHistVx);\r | |
333 | fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r | |
334 | fList->Add(fHistVy);\r | |
6b046844 | 335 | fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);\r |
0879e280 | 336 | fList->Add(fHistVz);\r |
337 | \r | |
4d528efb | 338 | //TPC vs VZERO multiplicity\r |
339 | fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",3001,-0.5,30000.5,4001,-0.5,4000.5);\r | |
340 | if(fMultiplicityEstimator == "V0A") \r | |
c16b246d | 341 | fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)");\r |
4d528efb | 342 | else if(fMultiplicityEstimator == "V0C") \r |
c16b246d | 343 | fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)");\r |
4d528efb | 344 | else \r |
345 | fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)");\r | |
346 | fList->Add(fHistTPCvsVZEROMultiplicity);\r | |
347 | \r | |
348 | fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);\r | |
349 | fList->Add(fHistVZEROSignal);\r | |
350 | \r | |
0879e280 | 351 | //Event plane\r |
6b046844 | 352 | fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);\r |
0879e280 | 353 | fList->Add(fHistEventPlane);\r |
354 | \r | |
355 | // QA histograms\r | |
356 | fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r | |
357 | fList->Add(fHistClus);\r | |
6b046844 | 358 | fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);\r |
0879e280 | 359 | fList->Add(fHistChi2);\r |
360 | fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r | |
361 | fList->Add(fHistDCA);\r | |
6b046844 | 362 | fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);\r |
0879e280 | 363 | fList->Add(fHistPt);\r |
6b046844 | 364 | fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);\r |
0879e280 | 365 | fList->Add(fHistEta);\r |
6b046844 | 366 | fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);\r |
244c1a3d | 367 | fList->Add(fHistRapidity);\r |
c8ad9635 | 368 | fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r |
0879e280 | 369 | fList->Add(fHistPhi);\r |
f2e8af26 | 370 | 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 |
371 | fList->Add(fHistEtaPhiPos); \r | |
372 | 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 | |
373 | fList->Add(fHistEtaPhiNeg);\r | |
6b046844 | 374 | fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r |
0879e280 | 375 | fList->Add(fHistPhiBefore);\r |
6b046844 | 376 | fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r |
0879e280 | 377 | fList->Add(fHistPhiAfter);\r |
c8ad9635 | 378 | fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r |
244c1a3d | 379 | fList->Add(fHistPhiPos);\r |
c8ad9635 | 380 | fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);\r |
244c1a3d | 381 | fList->Add(fHistPhiNeg);\r |
0879e280 | 382 | fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r |
383 | fList->Add(fHistV0M);\r | |
384 | TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r | |
385 | fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r | |
386 | for(Int_t i = 1; i <= 6; i++)\r | |
387 | fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r | |
388 | fList->Add(fHistRefTracks);\r | |
389 | \r | |
ea9867da | 390 | // Balance function histograms\r |
391 | // Initialize histograms if not done yet (including the custom binning)\r | |
392 | if(!fBalance->GetHistNp()){\r | |
393 | AliInfo("Histograms not yet initialized! --> Will be done now");\r | |
394 | fBalance->SetCustomBinning(fCustomBinning);\r | |
395 | fBalance->InitHistograms();\r | |
396 | }\r | |
397 | \r | |
398 | if(fRunShuffling) {\r | |
399 | if(!fShuffledBalance->GetHistNp()) {\r | |
400 | AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now");\r | |
401 | fShuffledBalance->SetCustomBinning(fCustomBinning);\r | |
402 | fShuffledBalance->InitHistograms();\r | |
403 | }\r | |
404 | }\r | |
405 | \r | |
406 | if(fRunMixing) {\r | |
407 | if(!fMixedBalance->GetHistNp()) {\r | |
408 | AliInfo("Histograms (mixing) not yet initialized! --> Will be done now");\r | |
409 | fMixedBalance->SetCustomBinning(fCustomBinning);\r | |
410 | fMixedBalance->InitHistograms();\r | |
411 | }\r | |
412 | }\r | |
413 | \r | |
5a861dc6 | 414 | // QA histograms for different cuts\r |
73609a48 | 415 | fList->Add(fBalance->GetQAHistHBTbefore());\r |
416 | fList->Add(fBalance->GetQAHistHBTafter());\r | |
417 | fList->Add(fBalance->GetQAHistConversionbefore());\r | |
418 | fList->Add(fBalance->GetQAHistConversionafter());\r | |
419 | fList->Add(fBalance->GetQAHistPsiMinusPhi());\r | |
f0754617 | 420 | fList->Add(fBalance->GetQAHistResonancesBefore());\r |
421 | fList->Add(fBalance->GetQAHistResonancesRho());\r | |
422 | fList->Add(fBalance->GetQAHistResonancesK0());\r | |
423 | fList->Add(fBalance->GetQAHistResonancesLambda());\r | |
5a861dc6 | 424 | fList->Add(fBalance->GetQAHistQbefore());\r |
425 | fList->Add(fBalance->GetQAHistQafter());\r | |
73609a48 | 426 | \r |
9afe3098 | 427 | //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r |
428 | fListBF->Add(fBalance->GetHistNp());\r | |
429 | fListBF->Add(fBalance->GetHistNn());\r | |
430 | fListBF->Add(fBalance->GetHistNpn());\r | |
431 | fListBF->Add(fBalance->GetHistNnn());\r | |
432 | fListBF->Add(fBalance->GetHistNpp());\r | |
433 | fListBF->Add(fBalance->GetHistNnp());\r | |
434 | \r | |
435 | if(fRunShuffling) {\r | |
436 | fListBFS->Add(fShuffledBalance->GetHistNp());\r | |
437 | fListBFS->Add(fShuffledBalance->GetHistNn());\r | |
438 | fListBFS->Add(fShuffledBalance->GetHistNpn());\r | |
439 | fListBFS->Add(fShuffledBalance->GetHistNnn());\r | |
440 | fListBFS->Add(fShuffledBalance->GetHistNpp());\r | |
441 | fListBFS->Add(fShuffledBalance->GetHistNnp());\r | |
442 | } \r | |
f06d59b3 | 443 | \r |
444 | if(fRunMixing) {\r | |
445 | fListBFM->Add(fMixedBalance->GetHistNp());\r | |
446 | fListBFM->Add(fMixedBalance->GetHistNn());\r | |
447 | fListBFM->Add(fMixedBalance->GetHistNpn());\r | |
448 | fListBFM->Add(fMixedBalance->GetHistNnn());\r | |
449 | fListBFM->Add(fMixedBalance->GetHistNpp());\r | |
450 | fListBFM->Add(fMixedBalance->GetHistNnp());\r | |
f0fb4ac1 | 451 | }\r |
9afe3098 | 452 | //}\r |
0879e280 | 453 | \r |
f06d59b3 | 454 | \r |
455 | // Event Mixing\r | |
8cdfbf52 | 456 | if(fRunMixing){\r |
457 | Int_t trackDepth = fMixingTracks; \r | |
458 | Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r | |
459 | \r | |
460 | // centrality bins\r | |
ea9867da | 461 | Double_t* centbins;\r |
462 | Int_t nCentralityBins;\r | |
463 | if(fBalance->IsUseVertexBinning()){\r | |
464 | centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins);\r | |
465 | }\r | |
466 | else{\r | |
467 | centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins);\r | |
468 | }\r | |
469 | \r | |
a9f20288 | 470 | // multiplicity bins\r |
ea9867da | 471 | Double_t* multbins;\r |
472 | Int_t nMultiplicityBins;\r | |
473 | multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins);\r | |
8cdfbf52 | 474 | \r |
475 | // Zvtx bins\r | |
ea9867da | 476 | Double_t* vtxbins; \r |
477 | Int_t nVertexBins;\r | |
478 | if(fBalance->IsUseVertexBinning()){\r | |
479 | vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins);\r | |
480 | }\r | |
481 | else{\r | |
482 | vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins);\r | |
483 | }\r | |
484 | \r | |
8cdfbf52 | 485 | // Event plane angle (Psi) bins\r |
ea9867da | 486 | Double_t* psibins;\r |
487 | Int_t nPsiBins; \r | |
488 | psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins);\r | |
489 | \r | |
8cdfbf52 | 490 | \r |
491 | // run the event mixing also in bins of event plane (statistics!)\r | |
492 | if(fRunMixingEventPlane){\r | |
a9f20288 | 493 | if(fEventClass=="Multiplicity"){\r |
494 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r | |
495 | }\r | |
496 | else{\r | |
497 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r | |
498 | }\r | |
8cdfbf52 | 499 | }\r |
500 | else{\r | |
a9f20288 | 501 | if(fEventClass=="Multiplicity"){\r |
502 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r | |
503 | }\r | |
504 | else{\r | |
505 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r | |
506 | }\r | |
8cdfbf52 | 507 | }\r |
508 | }\r | |
f06d59b3 | 509 | \r |
0879e280 | 510 | if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r |
511 | \r | |
512 | //====================PID========================//\r | |
513 | if(fUsePID) {\r | |
514 | fPIDCombined = new AliPIDCombined();\r | |
515 | fPIDCombined->SetDefaultTPCPriors();\r | |
516 | \r | |
517 | fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r | |
90bec8c5 | 518 | fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);\r |
0879e280 | 519 | \r |
520 | fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r | |
90bec8c5 | 521 | fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); \r |
0879e280 | 522 | \r |
523 | fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r | |
90bec8c5 | 524 | fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); \r |
0879e280 | 525 | \r |
526 | fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r | |
90bec8c5 | 527 | fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID);\r |
0879e280 | 528 | \r |
529 | fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r | |
90bec8c5 | 530 | fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID);\r |
0879e280 | 531 | \r |
532 | fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r | |
90bec8c5 | 533 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID);\r |
0879e280 | 534 | \r |
535 | fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r | |
90bec8c5 | 536 | fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); \r |
537 | \r | |
538 | fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2); \r | |
539 | fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++\r | |
540 | \r | |
541 | fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500); \r | |
542 | fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++\r | |
0879e280 | 543 | \r |
544 | fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r | |
90bec8c5 | 545 | fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);\r |
0879e280 | 546 | \r |
547 | fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r | |
90bec8c5 | 548 | fHistListPIDQA->Add(fHistBetavsPTOFafterPID); \r |
0879e280 | 549 | \r |
550 | fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r | |
90bec8c5 | 551 | fHistListPIDQA->Add(fHistProbTPCvsPtafterPID);\r |
0879e280 | 552 | \r |
553 | fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); \r | |
90bec8c5 | 554 | fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); \r |
0879e280 | 555 | \r |
556 | fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r | |
90bec8c5 | 557 | fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID);\r |
0879e280 | 558 | \r |
559 | fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r | |
90bec8c5 | 560 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID);\r |
0879e280 | 561 | \r |
562 | fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r | |
90bec8c5 | 563 | fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID);\r |
564 | \r | |
565 | fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2); \r | |
566 | fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++\r | |
567 | \r | |
568 | fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500); \r | |
569 | fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++\r | |
0879e280 | 570 | }\r |
c443da34 | 571 | \r |
572 | // for electron rejection only TPC nsigma histograms\r | |
573 | if(!fUsePID && fElectronRejection) {\r | |
574 | \r | |
575 | fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r | |
90bec8c5 | 576 | fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID);\r |
c443da34 | 577 | \r |
578 | fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r | |
579 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r | |
580 | \r | |
581 | fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r | |
90bec8c5 | 582 | fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID);\r |
c443da34 | 583 | \r |
584 | fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r | |
90bec8c5 | 585 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); \r |
c443da34 | 586 | }\r |
0879e280 | 587 | //====================PID========================//\r |
588 | \r | |
589 | // Post output data.\r | |
590 | PostData(1, fList);\r | |
591 | PostData(2, fListBF);\r | |
592 | if(fRunShuffling) PostData(3, fListBFS);\r | |
f06d59b3 | 593 | if(fRunMixing) PostData(4, fListBFM);\r |
c443da34 | 594 | if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID\r |
211b716d | 595 | \r |
ea9867da | 596 | AliInfo("Finished setting up the Output");\r |
597 | \r | |
211b716d | 598 | TH1::AddDirectory(oldStatus);\r |
0879e280 | 599 | }\r |
600 | \r | |
ce8126c8 | 601 | \r |
35aff0f3 | 602 | //________________________________________________________________________\r |
87dd0926 | 603 | void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename, \r |
c3b1f1a5 | 604 | Int_t nCentralityBins, \r |
605 | Double_t *centralityArrayForCorrections) {\r | |
35aff0f3 | 606 | //Open files that will be used for correction\r |
c3b1f1a5 | 607 | fCentralityArrayBinsForCorrections = nCentralityBins;\r |
608 | for (Int_t i=0; i<nCentralityBins; i++)\r | |
609 | fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];\r | |
87dd0926 | 610 | \r |
611 | // No file specified -> run without corrections\r | |
612 | if(!filename.Contains(".root")) {\r | |
613 | AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));\r | |
614 | return;\r | |
615 | }\r | |
616 | \r | |
35aff0f3 | 617 | //Open the input file\r |
618 | TFile *f = TFile::Open(filename);\r | |
619 | if(!f->IsOpen()) {\r | |
87dd0926 | 620 | AliInfo(Form("File %s not found --> run without corrections",filename.Data()));\r |
ce8126c8 | 621 | return;\r |
35aff0f3 | 622 | }\r |
1cb2a06e | 623 | \r |
624 | //TString listEffName = "";\r | |
c3b1f1a5 | 625 | for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) { \r |
626 | //Printf("iCent %d:",iCent); \r | |
1cb2a06e | 627 | TString histoName = "fHistCorrectionPlus";\r |
6fa567bd | 628 | histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r |
d79d2934 | 629 | fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));\r |
1cb2a06e | 630 | if(!fHistCorrectionPlus[iCent]) {\r |
c443da34 | 631 | AliError(Form("fHist %s not found!!!",histoName.Data()));\r |
35aff0f3 | 632 | return;\r |
35aff0f3 | 633 | }\r |
634 | \r | |
1cb2a06e | 635 | histoName = "fHistCorrectionMinus";\r |
6fa567bd | 636 | histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r |
d79d2934 | 637 | fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data())); \r |
1cb2a06e | 638 | if(!fHistCorrectionMinus[iCent]) {\r |
c443da34 | 639 | AliError(Form("fHist %s not found!!!",histoName.Data()));\r |
35aff0f3 | 640 | return; \r |
35aff0f3 | 641 | }\r |
642 | }//loop over centralities: ONLY the PbPb case is covered\r | |
643 | \r | |
1cb2a06e | 644 | if(fHistCorrectionPlus[0]){\r |
645 | fEtaMinForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmin();\r | |
646 | fEtaMaxForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmax();\r | |
647 | fEtaBinForCorrections = fHistCorrectionPlus[0]->GetNbinsX();\r | |
ce8126c8 | 648 | \r |
1cb2a06e | 649 | fPtMinForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmin();\r |
650 | fPtMaxForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmax();\r | |
651 | fPtBinForCorrections = fHistCorrectionPlus[0]->GetNbinsY();\r | |
ce8126c8 | 652 | \r |
1cb2a06e | 653 | fPhiMinForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmin();\r |
654 | fPhiMaxForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmax();\r | |
655 | fPhiBinForCorrections = fHistCorrectionPlus[0]->GetNbinsZ();\r | |
ce8126c8 | 656 | }\r |
35aff0f3 | 657 | }\r |
35aff0f3 | 658 | \r |
0879e280 | 659 | //________________________________________________________________________\r |
660 | void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r | |
661 | // Main loop\r | |
662 | // Called for each event\r | |
0879e280 | 663 | \r |
f06d59b3 | 664 | TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r |
0879e280 | 665 | Int_t gNumberOfAcceptedTracks = 0;\r |
90bec8c5 | 666 | Double_t lMultiplicityVar = -999.; //-1\r |
f06d59b3 | 667 | Double_t gReactionPlane = -1.; \r |
73609a48 | 668 | Float_t bSign = 0.;\r |
35aff0f3 | 669 | \r |
f06d59b3 | 670 | // get the event (for generator level: MCEvent())\r |
671 | AliVEvent* eventMain = NULL;\r | |
672 | if(gAnalysisLevel == "MC") {\r | |
673 | eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r | |
674 | }\r | |
675 | else{\r | |
90bec8c5 | 676 | eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r |
73609a48 | 677 | // for HBT like cuts need magnetic field sign\r |
678 | bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r | |
f06d59b3 | 679 | }\r |
680 | if(!eventMain) {\r | |
681 | AliError("eventMain not available");\r | |
682 | return;\r | |
683 | }\r | |
73609a48 | 684 | \r |
f06d59b3 | 685 | // PID Response task active?\r |
c443da34 | 686 | if(fUsePID || fElectronRejection) {\r |
0879e280 | 687 | fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r |
688 | if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r | |
689 | }\r | |
90bec8c5 | 690 | \r |
f06d59b3 | 691 | // check event cuts and fill event histograms\r |
90bec8c5 | 692 | if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){ \r |
f06d59b3 | 693 | return;\r |
694 | }\r | |
f06d59b3 | 695 | // get the reaction plane\r |
9efce966 | 696 | if(fEventClass != "Multiplicity") {\r |
697 | gReactionPlane = GetEventPlane(eventMain);\r | |
698 | fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r | |
699 | if(gReactionPlane < 0){\r | |
700 | return;\r | |
701 | }\r | |
f06d59b3 | 702 | }\r |
703 | \r | |
704 | // get the accepted tracks in main event\r | |
a9f20288 | 705 | TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r |
f06d59b3 | 706 | gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r |
0879e280 | 707 | \r |
f06d59b3 | 708 | //multiplicity cut (used in pp)\r |
a9f20288 | 709 | fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r |
0879e280 | 710 | \r |
9efce966 | 711 | // store charges of all accepted tracks,shuffle and reassign(two extra loops)\r |
f06d59b3 | 712 | TObjArray* tracksShuffled = NULL;\r |
713 | if(fRunShuffling){\r | |
a9f20288 | 714 | tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r |
f06d59b3 | 715 | }\r |
716 | \r | |
717 | // Event mixing \r | |
718 | if (fRunMixing)\r | |
719 | {\r | |
720 | // 1. First get an event pool corresponding in mult (cent) and\r | |
721 | // zvertex to the current event. Once initialized, the pool\r | |
722 | // should contain nMix (reduced) events. This routine does not\r | |
723 | // pre-scan the chain. The first several events of every chain\r | |
724 | // will be skipped until the needed pools are filled to the\r | |
725 | // specified depth. If the pool categories are not too rare, this\r | |
726 | // should not be a problem. If they are rare, you could lose`\r | |
727 | // statistics.\r | |
728 | \r | |
729 | // 2. Collect the whole pool's content of tracks into one TObjArray\r | |
730 | // (bgTracks), which is effectively a single background super-event.\r | |
731 | \r | |
732 | // 3. The reduced and bgTracks arrays must both be passed into\r | |
733 | // FillCorrelations(). Also nMix should be passed in, so a weight\r | |
734 | // of 1./nMix can be applied.\r | |
735 | \r | |
a9f20288 | 736 | AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r |
f06d59b3 | 737 | \r |
738 | if (!pool){\r | |
a9f20288 | 739 | AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r |
f06d59b3 | 740 | }\r |
741 | else{\r | |
0879e280 | 742 | \r |
f06d59b3 | 743 | //pool->SetDebug(1);\r |
f0fb4ac1 | 744 | \r |
f06d59b3 | 745 | if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r |
746 | \r | |
747 | \r | |
748 | Int_t nMix = pool->GetCurrentNEvents();\r | |
749 | //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r | |
750 | \r | |
751 | //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r | |
752 | //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r | |
753 | //if (pool->IsReady())\r | |
754 | //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r | |
755 | \r | |
756 | // Fill mixed-event histos here \r | |
757 | for (Int_t jMix=0; jMix<nMix; jMix++) \r | |
758 | {\r | |
759 | TObjArray* tracksMixed = pool->GetEvent(jMix);\r | |
683cbfb5 | 760 | fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r |
f06d59b3 | 761 | }\r |
762 | }\r | |
0879e280 | 763 | \r |
f06d59b3 | 764 | // Update the Event pool\r |
765 | pool->UpdatePool(tracksMain);\r | |
766 | //pool->PrintInfo();\r | |
0879e280 | 767 | \r |
f06d59b3 | 768 | }//pool NULL check \r |
769 | }//run mixing\r | |
770 | \r | |
771 | // calculate balance function\r | |
683cbfb5 | 772 | fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r |
f06d59b3 | 773 | \r |
774 | // calculate shuffled balance function\r | |
775 | if(fRunShuffling && tracksShuffled != NULL) {\r | |
683cbfb5 | 776 | fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r |
f06d59b3 | 777 | }\r |
778 | } \r | |
779 | \r | |
780 | //________________________________________________________________________\r | |
781 | Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r | |
782 | // Checks the Event cuts\r | |
783 | // Fills Event statistics histograms\r | |
784 | \r | |
785 | Bool_t isSelectedMain = kTRUE;\r | |
9efce966 | 786 | Float_t gRefMultiplicity = -1.;\r |
f06d59b3 | 787 | TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r |
788 | \r | |
35909ee1 | 789 | AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);\r |
790 | \r | |
c16b246d | 791 | fHistEventStats->Fill(1,gRefMultiplicity); //all events\r |
f06d59b3 | 792 | \r |
7556a0f0 | 793 | // check first event in chunk (is not needed for new reconstructions)\r |
794 | if(fCheckFirstEventInChunk){\r | |
795 | AliAnalysisUtils ut;\r | |
796 | if(ut.IsFirstEventInChunk(event)) \r | |
797 | return -1.;\r | |
c16b246d | 798 | fHistEventStats->Fill(6,gRefMultiplicity); \r |
7556a0f0 | 799 | }\r |
7556a0f0 | 800 | // check for pile-up event\r |
801 | if(fCheckPileUp){\r | |
802 | AliAnalysisUtils ut;\r | |
803 | ut.SetUseMVPlpSelection(kTRUE);\r | |
804 | ut.SetUseOutOfBunchPileUp(kTRUE);\r | |
805 | if(ut.IsPileUpEvent(event))\r | |
806 | return -1.;\r | |
c16b246d | 807 | fHistEventStats->Fill(7,gRefMultiplicity); \r |
7556a0f0 | 808 | }\r |
809 | \r | |
f06d59b3 | 810 | // Event trigger bits\r |
811 | fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r | |
812 | if(fUseOfflineTrigger)\r | |
813 | isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r | |
814 | \r | |
815 | if(isSelectedMain) {\r | |
c16b246d | 816 | fHistEventStats->Fill(2,gRefMultiplicity); //triggered events\r |
817 | \r | |
f06d59b3 | 818 | // Event Vertex MC\r |
7ecd1b5e | 819 | if(gAnalysisLevel == "MC") {\r |
9006a70a | 820 | if(!event) {\r |
821 | AliError("mcEvent not available");\r | |
822 | return 0x0;\r | |
823 | }\r | |
824 | \r | |
35909ee1 | 825 | if(mcevent){\r |
826 | AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());\r | |
88f76d4b | 827 | if(header){ \r |
828 | TArrayF gVertexArray;\r | |
829 | header->PrimaryVertex(gVertexArray);\r | |
830 | //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r | |
831 | //gVertexArray.At(0),\r | |
832 | //gVertexArray.At(1),\r | |
833 | //gVertexArray.At(2));\r | |
c16b246d | 834 | fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r |
88f76d4b | 835 | if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r |
836 | if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r | |
837 | if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r | |
c16b246d | 838 | fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r |
839 | \r | |
840 | // get the reference multiplicty or centrality\r | |
841 | gRefMultiplicity = GetRefMultiOrCentrality(event);\r | |
842 | \r | |
88f76d4b | 843 | fHistVx->Fill(gVertexArray.At(0));\r |
844 | fHistVy->Fill(gVertexArray.At(1));\r | |
c16b246d | 845 | fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity);\r |
88f76d4b | 846 | \r |
847 | // take only events inside centrality class\r | |
9efce966 | 848 | if(fUseCentrality) {\r |
c16b246d | 849 | if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){\r |
850 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r | |
851 | return gRefMultiplicity; \r | |
9efce966 | 852 | }//centrality class\r |
853 | }\r | |
854 | // take events only within the same multiplicity class\r | |
855 | else if(fUseMultiplicity) {\r | |
58b76b17 | 856 | if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r |
9efce966 | 857 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r |
858 | return gRefMultiplicity;\r | |
859 | }\r | |
860 | }//multiplicity range\r | |
88f76d4b | 861 | }//Vz cut\r |
862 | }//Vy cut\r | |
863 | }//Vx cut\r | |
864 | }//header \r | |
865 | }//MC event object\r | |
f06d59b3 | 866 | }//MC\r |
867 | \r | |
868 | // Event Vertex AOD, ESD, ESDMC\r | |
869 | else{\r | |
870 | const AliVVertex *vertex = event->GetPrimaryVertex();\r | |
871 | \r | |
0879e280 | 872 | if(vertex) {\r |
f06d59b3 | 873 | Double32_t fCov[6];\r |
874 | vertex->GetCovarianceMatrix(fCov);\r | |
0879e280 | 875 | if(vertex->GetNContributors() > 0) {\r |
f06d59b3 | 876 | if(fCov[5] != 0) {\r |
c16b246d | 877 | fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r |
f06d59b3 | 878 | if(TMath::Abs(vertex->GetX()) < fVxMax) {\r |
879 | if(TMath::Abs(vertex->GetY()) < fVyMax) {\r | |
880 | if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r | |
c16b246d | 881 | fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r |
882 | \r | |
883 | // get the reference multiplicty or centrality\r | |
884 | gRefMultiplicity = GetRefMultiOrCentrality(event);\r | |
885 | \r | |
f06d59b3 | 886 | fHistVx->Fill(vertex->GetX());\r |
887 | fHistVy->Fill(vertex->GetY());\r | |
c16b246d | 888 | fHistVz->Fill(vertex->GetZ(),gRefMultiplicity);\r |
0879e280 | 889 | \r |
f06d59b3 | 890 | // take only events inside centrality class\r |
9efce966 | 891 | if(fUseCentrality) {\r |
c16b246d | 892 | if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){\r |
893 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality\r | |
894 | return gRefMultiplicity; \r | |
9efce966 | 895 | }//centrality class\r |
896 | }\r | |
897 | // take events only within the same multiplicity class\r | |
898 | else if(fUseMultiplicity) {\r | |
a3ed6fc1 | 899 | //if(fDebugLevel) \r |
900 | //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin,\r | |
901 | //fNumberOfAcceptedTracksMax,gRefMultiplicity);\r | |
4e33f021 | 902 | \r |
58b76b17 | 903 | if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r |
9efce966 | 904 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r |
905 | return gRefMultiplicity;\r | |
906 | }\r | |
907 | }//multiplicity range\r | |
0879e280 | 908 | }//Vz cut\r |
909 | }//Vy cut\r | |
910 | }//Vx cut\r | |
911 | }//proper vertex resolution\r | |
912 | }//proper number of contributors\r | |
913 | }//vertex object valid\r | |
914 | }//triggered event \r | |
f06d59b3 | 915 | }//AOD,ESD,ESDMC\r |
0879e280 | 916 | \r |
f06d59b3 | 917 | // in all other cases return -1 (event not accepted)\r |
918 | return -1;\r | |
919 | }\r | |
0879e280 | 920 | \r |
f0fb4ac1 | 921 | \r |
922 | //________________________________________________________________________\r | |
923 | Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r | |
924 | // Checks the Event cuts\r | |
925 | // Fills Event statistics histograms\r | |
926 | \r | |
35aff0f3 | 927 | Float_t gCentrality = -1.;\r |
c16b246d | 928 | Double_t gMultiplicity = -1.;\r |
f0fb4ac1 | 929 | TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r |
f0fb4ac1 | 930 | \r |
c16b246d | 931 | \r |
932 | // calculate centrality always (not only in centrality mode)\r | |
933 | if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++\r | |
934 | AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r | |
935 | if(header){\r | |
936 | gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r | |
937 | \r | |
938 | // QA for centrality estimators\r | |
939 | fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r | |
940 | fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));\r | |
941 | fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));\r | |
942 | fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r | |
943 | fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r | |
944 | fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); \r | |
945 | fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r | |
946 | fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r | |
947 | fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));\r | |
948 | fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));\r | |
949 | fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r | |
950 | fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r | |
951 | fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r | |
952 | \r | |
953 | // Centrality estimator USED ++++++++++++++++++++++++++++++\r | |
954 | fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));\r | |
955 | \r | |
956 | // centrality QA (V0M)\r | |
957 | fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r | |
958 | \r | |
959 | // centrality QA (reference tracks)\r | |
960 | fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r | |
961 | fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r | |
962 | fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r | |
963 | fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r | |
964 | fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r | |
965 | fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r | |
966 | fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r | |
967 | fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r | |
968 | fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r | |
969 | \r | |
970 | }//AOD header\r | |
971 | }//AOD\r | |
972 | \r | |
973 | else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r | |
974 | AliCentrality *centrality = event->GetCentrality();\r | |
975 | gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r | |
976 | \r | |
977 | // QA for centrality estimators\r | |
978 | fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r | |
979 | fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));\r | |
980 | fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));\r | |
981 | fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));\r | |
982 | fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));\r | |
983 | fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));\r | |
984 | fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));\r | |
985 | fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));\r | |
986 | fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));\r | |
987 | fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));\r | |
988 | fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));\r | |
989 | fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));\r | |
990 | fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r | |
f0fb4ac1 | 991 | \r |
c16b246d | 992 | // Centrality estimator USED ++++++++++++++++++++++++++++++\r |
993 | fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));\r | |
994 | \r | |
995 | // centrality QA (V0M)\r | |
996 | fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r | |
997 | }//ESD\r | |
998 | \r | |
999 | else if(gAnalysisLevel == "MC"){\r | |
1000 | Double_t gImpactParameter = 0.;\r | |
1001 | AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r | |
1002 | if(gMCEvent){\r | |
1003 | AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); \r | |
1004 | if(headerH){\r | |
1005 | gImpactParameter = headerH->ImpactParameter();\r | |
1006 | gCentrality = gImpactParameter;\r | |
1007 | }//MC header\r | |
1008 | }//MC event cast\r | |
1009 | }//MC\r | |
1010 | \r | |
1011 | else{\r | |
1012 | gCentrality = -1.;\r | |
1013 | }\r | |
1014 | \r | |
1015 | // calculate reference multiplicity always (not only in multiplicity mode)\r | |
1016 | if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){\r | |
88f76d4b | 1017 | AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r |
1018 | if(gESDEvent){\r | |
7ecd1b5e | 1019 | gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r |
526d35fd | 1020 | fHistMultiplicity->Fill(gMultiplicity);\r |
c286b45e | 1021 | }//AliESDevent cast\r |
c16b246d | 1022 | }//ESD mode\r |
1023 | \r | |
1024 | else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){\r | |
f0fb4ac1 | 1025 | AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r |
4d528efb | 1026 | if ((fMultiplicityEstimator == "V0M")||\r |
1027 | (fMultiplicityEstimator == "V0A")||\r | |
1028 | (fMultiplicityEstimator == "V0C") ||\r | |
a3ed6fc1 | 1029 | (fMultiplicityEstimator == "TPC")) {\r |
4d528efb | 1030 | gMultiplicity = GetReferenceMultiplicityFromAOD(event);\r |
a3ed6fc1 | 1031 | if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity);\r |
1032 | }\r | |
4d528efb | 1033 | else {\r |
1034 | if(header)\r | |
1035 | gMultiplicity = header->GetRefMultiplicity();\r | |
a3ed6fc1 | 1036 | if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity);\r |
4d528efb | 1037 | }\r |
526d35fd | 1038 | fHistMultiplicity->Fill(gMultiplicity);\r |
c16b246d | 1039 | }//AOD mode\r |
1040 | else if(gAnalysisLevel == "MC") {\r | |
9efce966 | 1041 | AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r |
7ecd1b5e | 1042 | //Calculating the multiplicity as the number of charged primaries\r |
1043 | //within \pm 0.8 in eta and pT > 0.1 GeV/c\r | |
1044 | for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {\r | |
1045 | AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));\r | |
1046 | if (!track) {\r | |
1047 | AliError(Form("Could not receive particle %d", iParticle));\r | |
1048 | continue;\r | |
1049 | }\r | |
1050 | \r | |
1051 | //exclude non stable particles\r | |
1052 | if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;\r | |
c16b246d | 1053 | \r |
6ec16c49 | 1054 | //++++++++++++++++\r |
4d528efb | 1055 | if (fMultiplicityEstimator == "V0M") {\r |
1056 | if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) \r | |
c16b246d | 1057 | continue;}\r |
4d528efb | 1058 | else if (fMultiplicityEstimator == "V0A") {\r |
1059 | if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;}\r | |
1060 | else if (fMultiplicityEstimator == "V0C") {\r | |
1061 | if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;}\r | |
1062 | else if (fMultiplicityEstimator == "TPC") {\r | |
82d26474 | 1063 | if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r |
1064 | if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;\r | |
1065 | }\r | |
6ec16c49 | 1066 | else{\r |
82d26474 | 1067 | if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue;\r |
6ec16c49 | 1068 | if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue;\r |
1069 | }\r | |
1070 | //++++++++++++++++\r | |
c16b246d | 1071 | \r |
7ecd1b5e | 1072 | if(track->Charge() == 0) continue;\r |
c16b246d | 1073 | \r |
7ecd1b5e | 1074 | gMultiplicity += 1;\r |
1075 | }//loop over primaries\r | |
6ec16c49 | 1076 | fHistMultiplicity->Fill(gMultiplicity);\r |
c16b246d | 1077 | }//MC mode\r |
1078 | else{\r | |
1079 | gMultiplicity = -1;\r | |
1080 | }\r | |
6ec16c49 | 1081 | \r |
c16b246d | 1082 | \r |
1083 | // decide what should be returned only here\r | |
f0fb4ac1 | 1084 | Double_t lReturnVal = -100;\r |
1085 | if(fEventClass=="Multiplicity"){\r | |
7ecd1b5e | 1086 | lReturnVal = gMultiplicity;\r |
f0fb4ac1 | 1087 | }else if(fEventClass=="Centrality"){\r |
35aff0f3 | 1088 | lReturnVal = gCentrality;\r |
f0fb4ac1 | 1089 | }\r |
1090 | return lReturnVal;\r | |
1091 | }\r | |
1092 | \r | |
4d528efb | 1093 | //________________________________________________________________________\r |
1094 | Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){\r | |
1095 | //Function that returns the reference multiplicity from AODs (data or reco MC)\r | |
1096 | //Different ref. mult. implemented: V0M, V0A, V0C, TPC\r | |
91e5e210 | 1097 | Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;\r |
1098 | Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;\r | |
1099 | \r | |
1100 | AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());\r | |
1101 | if(!header) {\r | |
1102 | Printf("ERROR: AOD header not available");\r | |
1103 | return -999;\r | |
1104 | }\r | |
1105 | Int_t gRunNumber = header->GetRunNumber();\r | |
b84c4097 | 1106 | \r |
1107 | // Loop over tracks in event\r | |
1108 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r | |
1109 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r | |
1110 | if (!aodTrack) {\r | |
1111 | AliError(Form("Could not receive track %d", iTracks));\r | |
1112 | continue;\r | |
1113 | }\r | |
1114 | \r | |
1115 | // AOD track cuts\r | |
1116 | if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r | |
4d528efb | 1117 | \r |
b84c4097 | 1118 | if(aodTrack->Charge() == 0) continue;\r |
1119 | // Kinematics cuts from ESD track cuts\r | |
1120 | if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue;\r | |
1121 | if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue;\r | |
1122 | \r | |
1123 | //=================PID (so far only for electron rejection)==========================//\r | |
1124 | if(fElectronRejection) {\r | |
1125 | // get the electron nsigma\r | |
1126 | Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r | |
4d528efb | 1127 | \r |
b84c4097 | 1128 | // check only for given momentum range\r |
1129 | if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){\r | |
1130 | //look only at electron nsigma\r | |
1131 | if(!fElectronOnlyRejection) {\r | |
1132 | //Make the decision based on the n-sigma of electrons\r | |
1133 | if(nSigma < fElectronRejectionNSigma) continue;\r | |
4d528efb | 1134 | }\r |
b84c4097 | 1135 | else {\r |
1136 | Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r | |
1137 | Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r | |
1138 | Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r | |
1139 | \r | |
1140 | //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r | |
1141 | if(nSigma < fElectronRejectionNSigma\r | |
1142 | && nSigmaPions > fElectronRejectionNSigma\r | |
1143 | && nSigmaKaons > fElectronRejectionNSigma\r | |
1144 | && nSigmaProtons > fElectronRejectionNSigma ) continue;\r | |
1145 | }\r | |
1146 | }\r | |
1147 | }//electron rejection\r | |
1148 | \r | |
1149 | gRefMultiplicityTPC += 1.0;\r | |
1150 | }// track loop\r | |
1151 | \r | |
1152 | //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)\r | |
91e5e210 | 1153 | for(Int_t iChannel = 0; iChannel < 64; iChannel++) {\r |
1154 | fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));\r | |
b84c4097 | 1155 | \r |
91e5e210 | 1156 | if(iChannel < 32) \r |
1157 | gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);\r | |
1158 | else if(iChannel >= 32) \r | |
1159 | gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);\r | |
b84c4097 | 1160 | }//loop over PMTs\r |
1161 | \r | |
91e5e210 | 1162 | //Equalization of gain\r |
1163 | Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");\r | |
1164 | if(gFactorA != 0)\r | |
1165 | gRefMultiplicityVZEROA /= gFactorA;\r | |
1166 | Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");\r | |
1167 | if(gFactorC != 0)\r | |
1168 | gRefMultiplicityVZEROC /= gFactorC;\r | |
1169 | if((gFactorA != 0)&&(gFactorC != 0)) \r | |
1170 | gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);\r | |
1171 | \r | |
a3ed6fc1 | 1172 | if(fDebugLevel) \r |
1173 | Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC);\r | |
1174 | \r | |
b84c4097 | 1175 | fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);\r |
1176 | \r | |
1177 | if(fMultiplicityEstimator == "TPC") \r | |
1178 | gRefMultiplicity = gRefMultiplicityTPC;\r | |
91e5e210 | 1179 | else if(fMultiplicityEstimator == "V0M")\r |
b84c4097 | 1180 | gRefMultiplicity = gRefMultiplicityVZERO;\r |
91e5e210 | 1181 | else if(fMultiplicityEstimator == "V0A")\r |
1182 | gRefMultiplicity = gRefMultiplicityVZEROA;\r | |
1183 | else if(fMultiplicityEstimator == "V0C")\r | |
1184 | gRefMultiplicity = gRefMultiplicityVZEROC;\r | |
4d528efb | 1185 | \r |
1186 | return gRefMultiplicity;\r | |
1187 | }\r | |
1188 | \r | |
f06d59b3 | 1189 | //________________________________________________________________________\r |
1190 | Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r | |
1191 | // Get the event plane\r | |
0879e280 | 1192 | \r |
f06d59b3 | 1193 | TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r |
0879e280 | 1194 | \r |
f06d59b3 | 1195 | Float_t gVZEROEventPlane = -10.;\r |
1196 | Float_t gReactionPlane = -10.;\r | |
1197 | Double_t qxTot = 0.0, qyTot = 0.0;\r | |
1198 | \r | |
1199 | //MC: from reaction plane\r | |
1200 | if(gAnalysisLevel == "MC"){\r | |
9006a70a | 1201 | if(!event) {\r |
1202 | AliError("mcEvent not available");\r | |
1203 | return 0x0;\r | |
1204 | }\r | |
1205 | \r | |
88f76d4b | 1206 | AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r |
1207 | if(gMCEvent){\r | |
1208 | AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); \r | |
951cf00e | 1209 | if (headerH) {\r |
1210 | gReactionPlane = headerH->ReactionPlaneAngle();\r | |
1211 | //gReactionPlane *= TMath::RadToDeg();\r | |
1212 | }//MC header\r | |
1213 | }//MC event cast\r | |
f06d59b3 | 1214 | }//MC\r |
1215 | \r | |
1216 | // AOD,ESD,ESDMC: from VZERO Event Plane\r | |
f0fb4ac1 | 1217 | else{\r |
f06d59b3 | 1218 | \r |
1219 | AliEventplane *ep = event->GetEventplane();\r | |
1220 | if(ep){ \r | |
1221 | gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r | |
1222 | if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r | |
c8ad9635 | 1223 | //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r |
1224 | gReactionPlane = gVZEROEventPlane;\r | |
0879e280 | 1225 | }\r |
c8ad9635 | 1226 | }//AOD,ESD,ESDMC\r |
0879e280 | 1227 | \r |
f06d59b3 | 1228 | return gReactionPlane;\r |
1229 | }\r | |
0879e280 | 1230 | \r |
ce8126c8 | 1231 | //________________________________________________________________________\r |
1232 | Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r | |
1233 | Double_t vPhi, \r | |
1234 | Double_t vPt, \r | |
1235 | Short_t vCharge, \r | |
1236 | Double_t gCentrality) {\r | |
7c04a4d2 | 1237 | // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r |
1238 | \r | |
35aff0f3 | 1239 | Double_t correction = 1.;\r |
ce8126c8 | 1240 | Int_t binEta = 0, binPt = 0, binPhi = 0;\r |
1241 | \r | |
1242 | //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r | |
1243 | if(fEtaBinForCorrections != 0) {\r | |
1244 | Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r | |
1245 | if(fEtaMaxForCorrections != fEtaMinForCorrections) \r | |
1cb2a06e | 1246 | binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;\r |
ce8126c8 | 1247 | }\r |
1248 | \r | |
1249 | if(fPtBinForCorrections != 0) {\r | |
1250 | Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r | |
1251 | if(fPtMaxForCorrections != fPtMinForCorrections) \r | |
1cb2a06e | 1252 | binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;\r |
ce8126c8 | 1253 | }\r |
35aff0f3 | 1254 | \r |
ce8126c8 | 1255 | if(fPhiBinForCorrections != 0) {\r |
1256 | Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r | |
1257 | if(fPhiMaxForCorrections != fPhiMinForCorrections) \r | |
1cb2a06e | 1258 | binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;\r |
ce8126c8 | 1259 | }\r |
35aff0f3 | 1260 | \r |
d79d2934 | 1261 | Int_t gCentralityInt = -1;\r |
6fa567bd | 1262 | for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){\r |
1263 | if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){\r | |
35aff0f3 | 1264 | gCentralityInt = i;\r |
1cb2a06e | 1265 | break;\r |
1266 | }\r | |
d79d2934 | 1267 | } \r |
1cb2a06e | 1268 | \r |
d79d2934 | 1269 | // centrality not in array --> no correction\r |
1270 | if(gCentralityInt < 0){\r | |
ce8126c8 | 1271 | correction = 1.;\r |
35aff0f3 | 1272 | }\r |
d79d2934 | 1273 | else{\r |
1274 | \r | |
1275 | //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);\r | |
1276 | \r | |
1277 | if(fHistCorrectionPlus[gCentralityInt]){\r | |
1278 | if (vCharge > 0) {\r | |
1279 | correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r | |
1280 | //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt); \r | |
1281 | }\r | |
1282 | if (vCharge < 0) {\r | |
1283 | correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r | |
1284 | //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);\r | |
1285 | }\r | |
1286 | }\r | |
1287 | else {\r | |
1288 | correction = 1.;\r | |
1289 | }\r | |
1290 | }//centrality in array\r | |
1291 | \r | |
35aff0f3 | 1292 | if (correction == 0.) { \r |
1293 | AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r | |
1294 | correction = 1.; \r | |
1295 | } \r | |
d79d2934 | 1296 | \r |
35aff0f3 | 1297 | return correction;\r |
1298 | }\r | |
35aff0f3 | 1299 | \r |
f06d59b3 | 1300 | //________________________________________________________________________\r |
35aff0f3 | 1301 | TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r |
f06d59b3 | 1302 | // Returns TObjArray with tracks after all track cuts (only for AOD!)\r |
1303 | // Fills QA histograms\r | |
0879e280 | 1304 | \r |
f06d59b3 | 1305 | TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r |
0879e280 | 1306 | \r |
f06d59b3 | 1307 | //output TObjArray holding all good tracks\r |
1308 | TObjArray* tracksAccepted = new TObjArray;\r | |
1309 | tracksAccepted->SetOwner(kTRUE);\r | |
0879e280 | 1310 | \r |
35aff0f3 | 1311 | Short_t vCharge;\r |
9fd4b54e | 1312 | Double_t vEta;\r |
1313 | Double_t vY;\r | |
1314 | Double_t vPhi;\r | |
1315 | Double_t vPt;\r | |
0879e280 | 1316 | \r |
f06d59b3 | 1317 | \r |
1318 | if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r | |
f06d59b3 | 1319 | // Loop over tracks in event\r |
90bec8c5 | 1320 | \r |
f06d59b3 | 1321 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r |
1322 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r | |
1323 | if (!aodTrack) {\r | |
1324 | AliError(Form("Could not receive track %d", iTracks));\r | |
1325 | continue;\r | |
1326 | }\r | |
1327 | \r | |
1328 | // AOD track cuts\r | |
1329 | \r | |
1330 | // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r | |
1331 | // take only TPC only tracks \r | |
c1847400 | 1332 | //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r |
1333 | for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r | |
1334 | fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r | |
1335 | }\r | |
90bec8c5 | 1336 | \r |
4d528efb | 1337 | if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r |
8b3da35d | 1338 | \r |
1339 | vCharge = aodTrack->Charge();\r | |
1340 | vEta = aodTrack->Eta();\r | |
1341 | vY = aodTrack->Y();\r | |
1342 | vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r | |
1343 | vPt = aodTrack->Pt();\r | |
1344 | \r | |
c443da34 | 1345 | //===========================PID (so far only for electron rejection)===============================// \r |
1346 | if(fElectronRejection) {\r | |
8b3da35d | 1347 | \r |
c443da34 | 1348 | // get the electron nsigma\r |
1349 | Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r | |
1350 | \r | |
1351 | //Fill QA before the PID\r | |
1352 | fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r | |
1353 | fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r | |
1354 | //end of QA-before pid\r | |
1355 | \r | |
8b3da35d | 1356 | // check only for given momentum range\r |
1357 | if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r | |
1358 | \r | |
1359 | //look only at electron nsigma\r | |
1360 | if(!fElectronOnlyRejection){\r | |
1361 | \r | |
1362 | //Make the decision based on the n-sigma of electrons\r | |
1363 | if(nSigma < fElectronRejectionNSigma) continue;\r | |
1364 | }\r | |
1365 | else{\r | |
1366 | \r | |
1367 | Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r | |
1368 | Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r | |
1369 | Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r | |
1370 | \r | |
1371 | //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r | |
1372 | if(nSigma < fElectronRejectionNSigma\r | |
1373 | && nSigmaPions > fElectronRejectionNSigma\r | |
1374 | && nSigmaKaons > fElectronRejectionNSigma\r | |
1375 | && nSigmaProtons > fElectronRejectionNSigma ) continue;\r | |
1376 | }\r | |
1377 | }\r | |
1378 | \r | |
c443da34 | 1379 | //Fill QA after the PID\r |
1380 | fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r | |
1381 | fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r | |
8b3da35d | 1382 | \r |
c443da34 | 1383 | }\r |
1384 | //===========================end of PID (so far only for electron rejection)===============================//\r | |
8b3da35d | 1385 | \r |
90bec8c5 | 1386 | //+++++++++++++++++++++++++++++//\r |
1387 | //===========================PID===============================// \r | |
1388 | if(fUsePID) {\r | |
1389 | Double_t prob[AliPID::kSPECIES]={0.};\r | |
1390 | Double_t probTPC[AliPID::kSPECIES]={0.};\r | |
1391 | Double_t probTOF[AliPID::kSPECIES]={0.};\r | |
1392 | Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r | |
1393 | \r | |
1394 | Double_t nSigma = 0.;\r | |
1395 | Double_t nSigmaTPC = 0.; //++++\r | |
1396 | Double_t nSigmaTOF = 0.; //++++\r | |
1397 | UInt_t detUsedTPC = 0;\r | |
1398 | UInt_t detUsedTOF = 0;\r | |
1399 | UInt_t detUsedTPCTOF = 0;\r | |
1400 | \r | |
1401 | //Decide what detector configuration we want to use\r | |
1402 | switch(fPidDetectorConfig) {\r | |
1403 | case kTPCpid:\r | |
1404 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r | |
1405 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest));\r | |
1406 | detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);\r | |
1407 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r | |
1408 | prob[iSpecies] = probTPC[iSpecies];\r | |
1409 | break;\r | |
1410 | case kTOFpid:\r | |
1411 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r | |
1412 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest));\r | |
1413 | detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);\r | |
1414 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r | |
1415 | prob[iSpecies] = probTOF[iSpecies];\r | |
1416 | break;\r | |
1417 | case kTPCTOF:\r | |
1418 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r | |
1419 | nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest)); //++++++\r | |
1420 | nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest)); //++++++\r | |
1421 | nSigma = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF);//++++++\r | |
1422 | detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF);\r | |
1423 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r | |
1424 | prob[iSpecies] = probTPCTOF[iSpecies];\r | |
1425 | break;\r | |
1426 | default:\r | |
1427 | break;\r | |
1428 | }//end switch: define detector mask\r | |
1429 | \r | |
1430 | //Filling the PID QA\r | |
1431 | Double_t tofTime = -999., length = 999., tof = -999.;\r | |
1432 | Double_t c = TMath::C()*1.E-9;// m/ns\r | |
1433 | Double_t beta = -999.;\r | |
1434 | Double_t nSigmaTOFForParticleOfInterest = -999.;\r | |
1435 | if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) &&\r | |
1436 | (aodTrack->IsOn(AliAODTrack::kTIME)) ) { \r | |
1437 | tofTime = aodTrack->GetTOFsignal();//in ps\r | |
1438 | length = aodTrack->GetIntegratedLength();\r | |
1439 | tof = tofTime*1E-3; // ns \r | |
1440 | \r | |
1441 | if (tof <= 0) {\r | |
1442 | //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r | |
1443 | continue;\r | |
1444 | }\r | |
1445 | if (length <= 0){\r | |
1446 | //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r | |
1447 | continue;\r | |
1448 | }\r | |
1449 | \r | |
1450 | length = length*0.01; // in meters\r | |
1451 | tof = tof*c;\r | |
1452 | beta = length/tof;\r | |
1453 | \r | |
1454 | nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest);\r | |
1455 | fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);\r | |
1456 | fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);\r | |
1457 | fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOFForParticleOfInterest);\r | |
1458 | }//TOF signal \r | |
1459 | \r | |
1460 | Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest);\r | |
1461 | fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r | |
1462 | fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); \r | |
1463 | fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCForParticleOfInterest); \r | |
1464 | fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);\r | |
1465 | \r | |
1466 | //combined TPC&TOF\r | |
1467 | fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ \r | |
1468 | Double_t nSigmaTPCTOFForParticleOfInterest = -999.;//++++++++\r | |
1469 | nSigmaTPCTOFForParticleOfInterest = TMath::Sqrt(nSigmaTPCForParticleOfInterest*nSigmaTPCForParticleOfInterest + nSigmaTOFForParticleOfInterest*nSigmaTOFForParticleOfInterest);//++++++++\r | |
1470 | fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOFForParticleOfInterest); //++++++++\r | |
1471 | \r | |
1472 | //Printf("nSigma %lf",nSigma); //++++++++++++\r | |
1473 | //Printf("nSigmaTPCTOF %lf",nSigmaTPCTOFForParticleOfInterest); //++++++++++++\r | |
1474 | //end of QA-before pid\r | |
1475 | \r | |
1476 | if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r | |
1477 | //Make the decision based on the n-sigma\r | |
1478 | if(fUsePIDnSigma) {\r | |
1479 | if(nSigma > fPIDNSigma) continue; \r | |
1480 | \r | |
1481 | fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigma);\r | |
1482 | fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigma); \r | |
1483 | fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigma); //++++++++\r | |
1484 | }\r | |
1485 | //Make the decision based on the bayesian\r | |
1486 | else if(fUsePIDPropabilities) {\r | |
1487 | if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r | |
1488 | if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r | |
1489 | \r | |
1490 | fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]);\r | |
1491 | fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); \r | |
1492 | fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]);\r | |
1493 | \r | |
1494 | }\r | |
1495 | \r | |
1496 | //Printf("nSigmaAFTER %lf", nSigma); //++++++++++++\r | |
1497 | //Printf("nSigmaTPCTOFAFTER %lf", nSigmaTPCTOFForParticleOfInterest); //++++++++++++ \r | |
1498 | //Fill QA after the PID\r | |
1499 | fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);\r | |
1500 | fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r | |
1501 | fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ \r | |
1502 | }\r | |
1503 | }\r | |
1504 | //===========================PID===============================//\r | |
1505 | //+++++++++++++++++++++++++++++//\r | |
1506 | \r | |
1507 | \r | |
9fd4b54e | 1508 | Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r |
1509 | Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r | |
f06d59b3 | 1510 | \r |
1511 | \r | |
1512 | // Kinematics cuts from ESD track cuts\r | |
9fd4b54e | 1513 | if( vPt < fPtMin || vPt > fPtMax) continue;\r |
1514 | if( vEta < fEtaMin || vEta > fEtaMax) continue;\r | |
f06d59b3 | 1515 | \r |
1516 | // Extra DCA cuts (for systematic studies [!= -1])\r | |
1517 | if( fDCAxyCut != -1 && fDCAzCut != -1){\r | |
9fd4b54e | 1518 | if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r |
f06d59b3 | 1519 | continue; // 2D cut\r |
1520 | }\r | |
1521 | }\r | |
1522 | \r | |
1523 | // Extra TPC cuts (for systematic studies [!= -1])\r | |
1524 | if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r | |
1525 | continue;\r | |
0879e280 | 1526 | }\r |
f06d59b3 | 1527 | if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r |
1528 | continue;\r | |
1529 | }\r | |
1530 | \r | |
1531 | // fill QA histograms\r | |
1532 | fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r | |
9fd4b54e | 1533 | fHistDCA->Fill(dcaZ,dcaXY);\r |
35aff0f3 | 1534 | fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r |
1535 | fHistPt->Fill(vPt,gCentrality);\r | |
1536 | fHistEta->Fill(vEta,gCentrality);\r | |
1537 | fHistRapidity->Fill(vY,gCentrality);\r | |
1538 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r | |
1539 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r | |
1540 | fHistPhi->Fill(vPhi,gCentrality);\r | |
f2e8af26 | 1541 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r |
1542 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r | |
35aff0f3 | 1543 | \r |
1544 | //=======================================correction\r | |
1545 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r | |
1cb2a06e | 1546 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r |
f06d59b3 | 1547 | \r |
1548 | // add the track to the TObjArray\r | |
35aff0f3 | 1549 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r |
f06d59b3 | 1550 | }//track loop\r |
1551 | }// AOD analysis\r | |
0879e280 | 1552 | \r |
1cb2a06e | 1553 | //==============================================================================================================\r |
1554 | else if(gAnalysisLevel == "MCAOD") {\r | |
1555 | \r | |
1556 | AliMCEvent* mcEvent = MCEvent();\r | |
1557 | if (!mcEvent) {\r | |
6d6eb2df | 1558 | AliError("ERROR: Could not retrieve MC event");\r |
1cb2a06e | 1559 | }\r |
6d6eb2df | 1560 | else{\r |
1cb2a06e | 1561 | \r |
6d6eb2df | 1562 | for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r |
1563 | AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); \r | |
1564 | if (!aodTrack) {\r | |
1565 | AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));\r | |
1566 | continue;\r | |
1567 | }\r | |
1568 | \r | |
1569 | if(!aodTrack->IsPhysicalPrimary()) continue; \r | |
1570 | \r | |
1571 | vCharge = aodTrack->Charge();\r | |
1572 | vEta = aodTrack->Eta();\r | |
1573 | vY = aodTrack->Y();\r | |
1574 | vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r | |
1575 | vPt = aodTrack->Pt();\r | |
1576 | \r | |
1577 | // Kinematics cuts from ESD track cuts\r | |
1578 | if( vPt < fPtMin || vPt > fPtMax) continue;\r | |
1579 | if( vEta < fEtaMin || vEta > fEtaMax) continue;\r | |
938d8888 | 1580 | \r |
6d6eb2df | 1581 | // Remove neutral tracks\r |
1582 | if( vCharge == 0 ) continue;\r | |
1583 | \r | |
1584 | //Exclude resonances\r | |
1585 | if(fExcludeResonancesInMC) {\r | |
1586 | \r | |
1587 | Bool_t kExcludeParticle = kFALSE;\r | |
1588 | Int_t gMotherIndex = aodTrack->GetMother();\r | |
1589 | if(gMotherIndex != -1) {\r | |
1590 | AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r | |
1591 | if(motherTrack) {\r | |
1592 | Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r | |
1593 | //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r | |
1594 | //if(pdgCodeOfMother == 113) {\r | |
1595 | if(pdgCodeOfMother == 113 // rho0\r | |
1596 | || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r | |
1597 | // || pdgCodeOfMother == 221 // eta\r | |
1598 | // || pdgCodeOfMother == 331 // eta'\r | |
1599 | // || pdgCodeOfMother == 223 // omega\r | |
1600 | // || pdgCodeOfMother == 333 // phi\r | |
1601 | || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0\r | |
1602 | // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*\r | |
1603 | // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*\r | |
1604 | || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r | |
1605 | || pdgCodeOfMother == 111 // pi0 Dalitz\r | |
32fb422c | 1606 | || pdgCodeOfMother == 22 // photon\r |
6d6eb2df | 1607 | ) {\r |
1608 | kExcludeParticle = kTRUE;\r | |
1609 | }\r | |
87dd0926 | 1610 | }\r |
1611 | }\r | |
6d6eb2df | 1612 | \r |
1613 | //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r | |
1614 | if(kExcludeParticle) continue;\r | |
87dd0926 | 1615 | }\r |
8b3da35d | 1616 | \r |
1617 | //Exclude electrons with PDG\r | |
1618 | if(fExcludeElectronsInMC) {\r | |
1619 | \r | |
1620 | if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;\r | |
1621 | \r | |
1622 | }\r | |
938d8888 | 1623 | \r |
6d6eb2df | 1624 | // fill QA histograms\r |
1625 | fHistPt->Fill(vPt,gCentrality);\r | |
1626 | fHistEta->Fill(vEta,gCentrality);\r | |
1627 | fHistRapidity->Fill(vY,gCentrality);\r | |
1628 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r | |
1629 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r | |
1630 | fHistPhi->Fill(vPhi,gCentrality);\r | |
1631 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r | |
1632 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r | |
1633 | \r | |
1634 | //=======================================correction\r | |
1635 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r | |
1636 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); \r | |
1637 | \r | |
1638 | // add the track to the TObjArray\r | |
1639 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r | |
1640 | }//aodTracks\r | |
1641 | }//MC event\r | |
1cb2a06e | 1642 | }//MCAOD\r |
1643 | //==============================================================================================================\r | |
1644 | \r | |
6f719a5b | 1645 | //==============================================================================================================\r |
1646 | else if(gAnalysisLevel == "MCAODrec") {\r | |
1647 | \r | |
1648 | /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());\r | |
1649 | if (!fAOD) {\r | |
1650 | printf("ERROR: fAOD not available\n");\r | |
1651 | return;\r | |
1652 | }*/\r | |
1653 | \r | |
1654 | fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); \r | |
1655 | if (!fArrayMC) { \r | |
1656 | AliError("No array of MC particles found !!!");\r | |
1657 | }\r | |
1658 | \r | |
1659 | AliMCEvent* mcEvent = MCEvent();\r | |
1660 | if (!mcEvent) {\r | |
1661 | AliError("ERROR: Could not retrieve MC event");\r | |
d9eb282c | 1662 | return tracksAccepted;\r |
6f719a5b | 1663 | }\r |
1664 | \r | |
1665 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r | |
1666 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r | |
1667 | if (!aodTrack) {\r | |
1668 | AliError(Form("Could not receive track %d", iTracks));\r | |
1669 | continue;\r | |
1670 | }\r | |
1671 | \r | |
1672 | for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r | |
1673 | fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r | |
1674 | }\r | |
4d528efb | 1675 | if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue;\r |
6f719a5b | 1676 | \r |
1677 | vCharge = aodTrack->Charge();\r | |
1678 | vEta = aodTrack->Eta();\r | |
1679 | vY = aodTrack->Y();\r | |
1680 | vPhi = aodTrack->Phi();// * TMath::RadToDeg();\r | |
1681 | vPt = aodTrack->Pt();\r | |
d9eb282c | 1682 | \r |
1683 | //===========================use MC information for Kinematics===============================// \r | |
1684 | if(fUseMCforKinematics){\r | |
1685 | \r | |
1686 | Int_t label = TMath::Abs(aodTrack->GetLabel());\r | |
1687 | AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r | |
1688 | \r | |
1689 | if(AODmcTrack){\r | |
1690 | vCharge = AODmcTrack->Charge();\r | |
1691 | vEta = AODmcTrack->Eta();\r | |
1692 | vY = AODmcTrack->Y();\r | |
1693 | vPhi = AODmcTrack->Phi();// * TMath::RadToDeg();\r | |
1694 | vPt = AODmcTrack->Pt();\r | |
1695 | }\r | |
1696 | else{\r | |
1697 | AliDebug(1, "no MC particle for this track"); \r | |
1698 | continue;\r | |
1699 | }\r | |
1700 | }\r | |
1701 | //===========================end of use MC information for Kinematics========================// \r | |
1702 | \r | |
8b3da35d | 1703 | \r |
1704 | //===========================PID (so far only for electron rejection)===============================// \r | |
1705 | if(fElectronRejection) {\r | |
1706 | \r | |
1707 | // get the electron nsigma\r | |
1708 | Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r | |
1709 | \r | |
1710 | //Fill QA before the PID\r | |
1711 | fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r | |
1712 | fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r | |
1713 | //end of QA-before pid\r | |
1714 | \r | |
1715 | // check only for given momentum range\r | |
1716 | if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){\r | |
1717 | \r | |
1718 | //look only at electron nsigma\r | |
1719 | if(!fElectronOnlyRejection){\r | |
1720 | \r | |
1721 | //Make the decision based on the n-sigma of electrons\r | |
1722 | if(nSigma < fElectronRejectionNSigma) continue;\r | |
1723 | }\r | |
1724 | else{\r | |
1725 | \r | |
1726 | Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion));\r | |
1727 | Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon));\r | |
1728 | Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton));\r | |
1729 | \r | |
1730 | //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species)\r | |
1731 | if(nSigma < fElectronRejectionNSigma\r | |
1732 | && nSigmaPions > fElectronRejectionNSigma\r | |
1733 | && nSigmaKaons > fElectronRejectionNSigma\r | |
1734 | && nSigmaProtons > fElectronRejectionNSigma ) continue;\r | |
1735 | }\r | |
1736 | }\r | |
1737 | \r | |
1738 | //Fill QA after the PID\r | |
1739 | fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r | |
1740 | fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r | |
1741 | \r | |
1742 | }\r | |
1743 | //===========================end of PID (so far only for electron rejection)===============================//\r | |
6f719a5b | 1744 | \r |
1745 | Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on)\r | |
1746 | Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on)\r | |
1747 | \r | |
1748 | // Kinematics cuts from ESD track cuts\r | |
1749 | if( vPt < fPtMin || vPt > fPtMax) continue;\r | |
1750 | if( vEta < fEtaMin || vEta > fEtaMax) continue;\r | |
1751 | \r | |
1752 | // Extra DCA cuts (for systematic studies [!= -1])\r | |
1753 | if( fDCAxyCut != -1 && fDCAzCut != -1){\r | |
1754 | if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r | |
1755 | continue; // 2D cut\r | |
1756 | }\r | |
1757 | }\r | |
1758 | \r | |
1759 | // Extra TPC cuts (for systematic studies [!= -1])\r | |
1760 | if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r | |
1761 | continue;\r | |
1762 | }\r | |
1763 | if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r | |
1764 | continue;\r | |
1765 | }\r | |
1766 | \r | |
1767 | //Exclude resonances\r | |
1768 | if(fExcludeResonancesInMC) {\r | |
1769 | \r | |
1770 | Bool_t kExcludeParticle = kFALSE;\r | |
1771 | \r | |
1772 | Int_t label = TMath::Abs(aodTrack->GetLabel());\r | |
1773 | AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r | |
1774 | \r | |
1775 | if (AODmcTrack){ \r | |
1776 | //if (AODmcTrack->IsPhysicalPrimary()){\r | |
1777 | \r | |
1778 | Int_t gMotherIndex = AODmcTrack->GetMother();\r | |
1779 | if(gMotherIndex != -1) {\r | |
1780 | AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r | |
1781 | if(motherTrack) {\r | |
1782 | Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r | |
1783 | if(pdgCodeOfMother == 113 // rho0\r | |
1784 | || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r | |
1785 | // || pdgCodeOfMother == 221 // eta\r | |
1786 | // || pdgCodeOfMother == 331 // eta'\r | |
1787 | // || pdgCodeOfMother == 223 // omega\r | |
1788 | // || pdgCodeOfMother == 333 // phi\r | |
1789 | || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0\r | |
1790 | // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*\r | |
1791 | // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*\r | |
1792 | || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r | |
1793 | || pdgCodeOfMother == 111 // pi0 Dalitz\r | |
6cf380f1 | 1794 | || pdgCodeOfMother == 22 // photon\r |
6f719a5b | 1795 | ) {\r |
1796 | kExcludeParticle = kTRUE;\r | |
1797 | }\r | |
1798 | }\r | |
1799 | }\r | |
1800 | } \r | |
1801 | //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r | |
1802 | if(kExcludeParticle) continue;\r | |
1803 | }\r | |
8b3da35d | 1804 | \r |
1805 | //Exclude electrons with PDG\r | |
1806 | if(fExcludeElectronsInMC) {\r | |
1807 | \r | |
1808 | Int_t label = TMath::Abs(aodTrack->GetLabel());\r | |
1809 | AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r | |
1810 | \r | |
1811 | if (AODmcTrack){ \r | |
1812 | if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;\r | |
1813 | }\r | |
1814 | }\r | |
6f719a5b | 1815 | \r |
1816 | // fill QA histograms\r | |
1817 | fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r | |
1818 | fHistDCA->Fill(dcaZ,dcaXY);\r | |
1819 | fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r | |
1820 | fHistPt->Fill(vPt,gCentrality);\r | |
1821 | fHistEta->Fill(vEta,gCentrality);\r | |
1822 | fHistRapidity->Fill(vY,gCentrality);\r | |
1823 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r | |
1824 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r | |
1825 | fHistPhi->Fill(vPhi,gCentrality);\r | |
1826 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); \r | |
1827 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r | |
1828 | \r | |
1829 | //=======================================correction\r | |
1830 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r | |
1831 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r | |
1832 | \r | |
1833 | // add the track to the TObjArray\r | |
1834 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r | |
1835 | }//track loop\r | |
1836 | }//MCAODrec\r | |
1837 | //==============================================================================================================\r | |
1838 | \r | |
f06d59b3 | 1839 | else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r |
1840 | \r | |
9fd4b54e | 1841 | AliESDtrack *trackTPC = NULL;\r |
f06d59b3 | 1842 | AliMCParticle *mcTrack = NULL;\r |
1843 | \r | |
1844 | AliMCEvent* mcEvent = NULL;\r | |
0879e280 | 1845 | \r |
f06d59b3 | 1846 | // for MC ESDs use also MC information\r |
1847 | if(gAnalysisLevel == "MCESD"){\r | |
1848 | mcEvent = MCEvent(); \r | |
1849 | if (!mcEvent) {\r | |
1850 | AliError("mcEvent not available");\r | |
1851 | return tracksAccepted;\r | |
1852 | }\r | |
1853 | }\r | |
0879e280 | 1854 | \r |
f06d59b3 | 1855 | // Loop over tracks in event\r |
1856 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r | |
1857 | AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r | |
1858 | if (!track) {\r | |
1859 | AliError(Form("Could not receive track %d", iTracks));\r | |
1860 | continue;\r | |
1861 | } \r | |
1862 | \r | |
1863 | // for MC ESDs use also MC information --> MC track not used further???\r | |
1864 | if(gAnalysisLevel == "MCESD"){\r | |
1865 | Int_t label = TMath::Abs(track->GetLabel());\r | |
1866 | if(label > mcEvent->GetNumberOfTracks()) continue;\r | |
1867 | if(label > mcEvent->GetNumberOfPrimaries()) continue;\r | |
1868 | \r | |
1869 | mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r | |
1870 | if(!mcTrack) continue;\r | |
1871 | }\r | |
0879e280 | 1872 | \r |
f06d59b3 | 1873 | // take only TPC only tracks\r |
9fd4b54e | 1874 | trackTPC = new AliESDtrack();\r |
1875 | if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r | |
f06d59b3 | 1876 | \r |
1877 | //ESD track cuts\r | |
1878 | if(fESDtrackCuts) \r | |
9fd4b54e | 1879 | if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r |
f06d59b3 | 1880 | \r |
1881 | // fill QA histograms\r | |
1882 | Float_t b[2];\r | |
1883 | Float_t bCov[3];\r | |
9fd4b54e | 1884 | trackTPC->GetImpactParameters(b,bCov);\r |
f06d59b3 | 1885 | if (bCov[0]<=0 || bCov[2]<=0) {\r |
1886 | AliDebug(1, "Estimated b resolution lower or equal zero!");\r | |
1887 | bCov[0]=0; bCov[2]=0;\r | |
1888 | }\r | |
1889 | \r | |
1890 | Int_t nClustersTPC = -1;\r | |
9fd4b54e | 1891 | nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone\r |
f06d59b3 | 1892 | //nClustersTPC = track->GetTPCclusters(0); // global track\r |
1893 | Float_t chi2PerClusterTPC = -1;\r | |
1894 | if (nClustersTPC!=0) {\r | |
9fd4b54e | 1895 | chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone\r |
f06d59b3 | 1896 | //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track\r |
1897 | }\r | |
1898 | \r | |
1899 | //===========================PID===============================// \r | |
1900 | if(fUsePID) {\r | |
1901 | Double_t prob[AliPID::kSPECIES]={0.};\r | |
1902 | Double_t probTPC[AliPID::kSPECIES]={0.};\r | |
1903 | Double_t probTOF[AliPID::kSPECIES]={0.};\r | |
1904 | Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r | |
1905 | \r | |
1906 | Double_t nSigma = 0.;\r | |
1907 | UInt_t detUsedTPC = 0;\r | |
1908 | UInt_t detUsedTOF = 0;\r | |
1909 | UInt_t detUsedTPCTOF = 0;\r | |
1910 | \r | |
1911 | //Decide what detector configuration we want to use\r | |
1912 | switch(fPidDetectorConfig) {\r | |
1913 | case kTPCpid:\r | |
1914 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r | |
1915 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r | |
1916 | detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r | |
1917 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r | |
1918 | prob[iSpecies] = probTPC[iSpecies];\r | |
1919 | break;\r | |
1920 | case kTOFpid:\r | |
1921 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r | |
1922 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r | |
1923 | detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r | |
1924 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r | |
1925 | prob[iSpecies] = probTOF[iSpecies];\r | |
1926 | break;\r | |
1927 | case kTPCTOF:\r | |
1928 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r | |
1929 | detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r | |
1930 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r | |
1931 | prob[iSpecies] = probTPCTOF[iSpecies];\r | |
1932 | break;\r | |
1933 | default:\r | |
1934 | break;\r | |
1935 | }//end switch: define detector mask\r | |
1936 | \r | |
1937 | //Filling the PID QA\r | |
1938 | Double_t tofTime = -999., length = 999., tof = -999.;\r | |
1939 | Double_t c = TMath::C()*1.E-9;// m/ns\r | |
1940 | Double_t beta = -999.;\r | |
1941 | Double_t nSigmaTOFForParticleOfInterest = -999.;\r | |
1942 | if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r | |
1943 | (track->IsOn(AliESDtrack::kTIME)) ) { \r | |
1944 | tofTime = track->GetTOFsignal();//in ps\r | |
1945 | length = track->GetIntegratedLength();\r | |
1946 | tof = tofTime*1E-3; // ns \r | |
1947 | \r | |
1948 | if (tof <= 0) {\r | |
1949 | //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r | |
1950 | continue;\r | |
1951 | }\r | |
1952 | if (length <= 0){\r | |
1953 | //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r | |
1954 | continue;\r | |
1955 | }\r | |
1956 | \r | |
1957 | length = length*0.01; // in meters\r | |
1958 | tof = tof*c;\r | |
1959 | beta = length/tof;\r | |
1960 | \r | |
1961 | nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r | |
1962 | fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r | |
1963 | fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r | |
1964 | fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r | |
1965 | }//TOF signal \r | |
1966 | \r | |
1967 | \r | |
1968 | Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r | |
1969 | fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r | |
1970 | fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r | |
1971 | fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r | |
1972 | fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r | |
1973 | //end of QA-before pid\r | |
1974 | \r | |
1975 | if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r | |
1976 | //Make the decision based on the n-sigma\r | |
1977 | if(fUsePIDnSigma) {\r | |
1978 | if(nSigma > fPIDNSigma) continue;}\r | |
1979 | \r | |
1980 | //Make the decision based on the bayesian\r | |
1981 | else if(fUsePIDPropabilities) {\r | |
1982 | if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r | |
1983 | if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; \r | |
1984 | }\r | |
1985 | \r | |
1986 | //Fill QA after the PID\r | |
1987 | fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r | |
1988 | fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r | |
1989 | fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r | |
1990 | \r | |
1991 | fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r | |
1992 | fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r | |
1993 | fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r | |
1994 | fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r | |
1995 | }\r | |
1996 | }\r | |
1997 | //===========================PID===============================//\r | |
9fd4b54e | 1998 | vCharge = trackTPC->Charge();\r |
1999 | vY = trackTPC->Y();\r | |
2000 | vEta = trackTPC->Eta();\r | |
c8ad9635 | 2001 | vPhi = trackTPC->Phi();// * TMath::RadToDeg();\r |
9fd4b54e | 2002 | vPt = trackTPC->Pt();\r |
2003 | fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r | |
f06d59b3 | 2004 | fHistDCA->Fill(b[1],b[0]);\r |
35aff0f3 | 2005 | fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r |
2006 | fHistPt->Fill(vPt,gCentrality);\r | |
2007 | fHistEta->Fill(vEta,gCentrality);\r | |
2008 | fHistPhi->Fill(vPhi,gCentrality);\r | |
f2e8af26 | 2009 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r |
2010 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r | |
35aff0f3 | 2011 | fHistRapidity->Fill(vY,gCentrality);\r |
2012 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r | |
2013 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r | |
f06d59b3 | 2014 | \r |
35aff0f3 | 2015 | //=======================================correction\r |
2016 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r | |
1cb2a06e | 2017 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r |
35aff0f3 | 2018 | \r |
f06d59b3 | 2019 | // add the track to the TObjArray\r |
35aff0f3 | 2020 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r |
244c1a3d | 2021 | \r |
9fd4b54e | 2022 | delete trackTPC;\r |
f06d59b3 | 2023 | }//track loop\r |
2024 | }// ESD analysis\r | |
244c1a3d | 2025 | \r |
73609a48 | 2026 | else if(gAnalysisLevel == "MC"){\r |
9006a70a | 2027 | if(!event) {\r |
2028 | AliError("mcEvent not available");\r | |
2029 | return 0x0;\r | |
2030 | }\r | |
2031 | \r | |
88f76d4b | 2032 | AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r |
2033 | if(gMCEvent) {\r | |
2034 | // Loop over tracks in event\r | |
2035 | for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r | |
2036 | AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r | |
2037 | if (!track) {\r | |
2038 | AliError(Form("Could not receive particle %d", iTracks));\r | |
f06d59b3 | 2039 | continue;\r |
88f76d4b | 2040 | }\r |
2041 | \r | |
2042 | //exclude non stable particles\r | |
2043 | if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r | |
2044 | \r | |
87dd0926 | 2045 | vCharge = track->Charge();\r |
88f76d4b | 2046 | vEta = track->Eta();\r |
2047 | vPt = track->Pt();\r | |
2048 | vY = track->Y();\r | |
2049 | \r | |
2050 | if( vPt < fPtMin || vPt > fPtMax) \r | |
f06d59b3 | 2051 | continue;\r |
88f76d4b | 2052 | if (!fUsePID) {\r |
2053 | if( vEta < fEtaMin || vEta > fEtaMax) continue;\r | |
2054 | }\r | |
2055 | else if (fUsePID){\r | |
2056 | if( vY < fEtaMin || vY > fEtaMax) continue;\r | |
2057 | }\r | |
87dd0926 | 2058 | \r |
2059 | // Remove neutral tracks\r | |
2060 | if( vCharge == 0 ) continue;\r | |
88f76d4b | 2061 | \r |
2062 | //analyze one set of particles\r | |
2063 | if(fUseMCPdgCode) {\r | |
2064 | TParticle *particle = track->Particle();\r | |
2065 | if(!particle) continue;\r | |
2066 | \r | |
2067 | Int_t gPdgCode = particle->GetPdgCode();\r | |
2068 | if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r | |
2069 | continue;\r | |
2070 | }\r | |
2071 | \r | |
2072 | //Use the acceptance parameterization\r | |
2073 | if(fAcceptanceParameterization) {\r | |
2074 | Double_t gRandomNumber = gRandom->Rndm();\r | |
2075 | if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r | |
2076 | continue;\r | |
2077 | }\r | |
f06d59b3 | 2078 | \r |
88f76d4b | 2079 | //Exclude resonances\r |
2080 | if(fExcludeResonancesInMC) {\r | |
2081 | TParticle *particle = track->Particle();\r | |
2082 | if(!particle) continue;\r | |
2083 | \r | |
2084 | Bool_t kExcludeParticle = kFALSE;\r | |
2085 | Int_t gMotherIndex = particle->GetFirstMother();\r | |
2086 | if(gMotherIndex != -1) {\r | |
2087 | AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r | |
2088 | if(motherTrack) {\r | |
2089 | TParticle *motherParticle = motherTrack->Particle();\r | |
2090 | if(motherParticle) {\r | |
2091 | Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r | |
2092 | //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r | |
97f46539 | 2093 | if(pdgCodeOfMother == 113 // rho0\r |
2094 | || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r | |
b7d8f442 | 2095 | // || pdgCodeOfMother == 221 // eta\r |
2096 | // || pdgCodeOfMother == 331 // eta'\r | |
2097 | // || pdgCodeOfMother == 223 // omega\r | |
2098 | // || pdgCodeOfMother == 333 // phi\r | |
97f46539 | 2099 | || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0\r |
b7d8f442 | 2100 | // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0*\r |
2101 | // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+*\r | |
2102 | || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r | |
97f46539 | 2103 | || pdgCodeOfMother == 111 // pi0 Dalitz\r |
2104 | ) {\r | |
88f76d4b | 2105 | kExcludeParticle = kTRUE;\r |
2106 | }\r | |
0879e280 | 2107 | }\r |
0879e280 | 2108 | }\r |
f06d59b3 | 2109 | }\r |
88f76d4b | 2110 | \r |
2111 | //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r | |
2112 | if(kExcludeParticle) continue;\r | |
f06d59b3 | 2113 | }\r |
8b3da35d | 2114 | \r |
2115 | //Exclude electrons with PDG\r | |
2116 | if(fExcludeElectronsInMC) {\r | |
2117 | \r | |
2118 | TParticle *particle = track->Particle();\r | |
2119 | \r | |
2120 | if (particle){ \r | |
2121 | if(TMath::Abs(particle->GetPdgCode()) == 11) continue;\r | |
2122 | }\r | |
2123 | }\r | |
2124 | \r | |
88f76d4b | 2125 | vPhi = track->Phi();\r |
2126 | //Printf("phi (before): %lf",vPhi);\r | |
f06d59b3 | 2127 | \r |
88f76d4b | 2128 | fHistPt->Fill(vPt,gCentrality);\r |
2129 | fHistEta->Fill(vEta,gCentrality);\r | |
2130 | fHistPhi->Fill(vPhi,gCentrality);\r | |
2131 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r | |
2132 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r | |
2133 | //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r | |
2134 | fHistRapidity->Fill(vY,gCentrality);\r | |
2135 | //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r | |
2136 | //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r | |
2137 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r | |
2138 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r | |
f06d59b3 | 2139 | \r |
88f76d4b | 2140 | //Flow after burner\r |
2141 | if(fUseFlowAfterBurner) {\r | |
2142 | Double_t precisionPhi = 0.001;\r | |
2143 | Int_t maxNumberOfIterations = 100;\r | |
2144 | \r | |
2145 | Double_t phi0 = vPhi;\r | |
2146 | Double_t gV2 = fDifferentialV2->Eval(vPt);\r | |
2147 | \r | |
2148 | for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r | |
2149 | Double_t phiprev = vPhi;\r | |
2150 | Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r | |
2151 | Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r | |
2152 | vPhi -= fl/fp;\r | |
2153 | if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r | |
2154 | }\r | |
2155 | //Printf("phi (after): %lf\n",vPhi);\r | |
2156 | Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r | |
2157 | if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r | |
2158 | fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r | |
2159 | \r | |
2160 | Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r | |
2161 | if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r | |
2162 | fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r | |
2163 | \r | |
f06d59b3 | 2164 | }\r |
f06d59b3 | 2165 | \r |
88f76d4b | 2166 | //vPhi *= TMath::RadToDeg();\r |
2167 | \r | |
87dd0926 | 2168 | //=======================================correction\r |
2169 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); \r | |
2170 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r | |
35aff0f3 | 2171 | \r |
88f76d4b | 2172 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r |
2173 | } //track loop\r | |
2174 | }//MC event object\r | |
f06d59b3 | 2175 | }//MC\r |
0879e280 | 2176 | \r |
f06d59b3 | 2177 | return tracksAccepted; \r |
2178 | }\r | |
88f76d4b | 2179 | \r |
f06d59b3 | 2180 | //________________________________________________________________________\r |
35aff0f3 | 2181 | TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r |
f06d59b3 | 2182 | // Clones TObjArray and returns it with tracks after shuffling the charges\r |
0879e280 | 2183 | \r |
f06d59b3 | 2184 | TObjArray* tracksShuffled = new TObjArray;\r |
2185 | tracksShuffled->SetOwner(kTRUE);\r | |
2186 | \r | |
2187 | vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks \r | |
2188 | \r | |
2189 | for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r | |
2190 | {\r | |
2191 | AliVParticle* track = (AliVParticle*) tracks->At(i);\r | |
2192 | chargeVector->push_back(track->Charge());\r | |
2193 | } \r | |
2194 | \r | |
2195 | random_shuffle(chargeVector->begin(), chargeVector->end());\r | |
2196 | \r | |
2197 | for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r | |
2198 | AliVParticle* track = (AliVParticle*) tracks->At(i);\r | |
35aff0f3 | 2199 | //==============================correction\r |
2200 | Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r | |
1cb2a06e | 2201 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r |
35aff0f3 | 2202 | tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r |
0879e280 | 2203 | }\r |
f06d59b3 | 2204 | \r |
2205 | delete chargeVector;\r | |
2206 | \r | |
2207 | return tracksShuffled;\r | |
2208 | }\r | |
2209 | \r | |
91e5e210 | 2210 | //________________________________________________________________________\r |
2211 | void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename,\r | |
2212 | const char* lhcPeriod) {\r | |
2213 | //Function to setup the VZERO gain equalization\r | |
2214 | //============Get the equilization map============//\r | |
2215 | TFile *calibrationFile = TFile::Open(filename);\r | |
2216 | if((!calibrationFile)||(!calibrationFile->IsOpen())) {\r | |
2217 | Printf("No calibration file found!!!");\r | |
2218 | return;\r | |
2219 | }\r | |
2220 | \r | |
2221 | TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));\r | |
2222 | if(!list) {\r | |
2223 | Printf("Calibration TList not found!!!");\r | |
2224 | return;\r | |
2225 | }\r | |
2226 | \r | |
2227 | fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));\r | |
2228 | if(!fHistVZEROAGainEqualizationMap) {\r | |
2229 | Printf("VZERO-A calibration object not found!!!");\r | |
2230 | return;\r | |
2231 | }\r | |
2232 | fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));\r | |
2233 | if(!fHistVZEROCGainEqualizationMap) {\r | |
2234 | Printf("VZERO-C calibration object not found!!!");\r | |
2235 | return;\r | |
2236 | }\r | |
2237 | \r | |
2238 | fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));\r | |
2239 | if(!fHistVZEROChannelGainEqualizationMap) {\r | |
2240 | Printf("VZERO channel calibration object not found!!!");\r | |
2241 | return;\r | |
2242 | }\r | |
2243 | }\r | |
2244 | \r | |
2245 | //________________________________________________________________________\r | |
2246 | Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run, \r | |
2247 | Int_t channel) {\r | |
2248 | //\r | |
46b3b28b | 2249 | if(!fHistVZEROAGainEqualizationMap) return 1.0;\r |
2250 | \r | |
91e5e210 | 2251 | for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {\r |
2252 | Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));\r | |
2253 | if(gRunNumber == run)\r | |
2254 | return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);\r | |
2255 | }\r | |
2256 | \r | |
2257 | return 1.0;\r | |
2258 | }\r | |
2259 | \r | |
2260 | //________________________________________________________________________\r | |
2261 | Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run, \r | |
2262 | const char* side) {\r | |
2263 | //\r | |
46b3b28b | 2264 | if(!fHistVZEROAGainEqualizationMap) return 1.0;\r |
2265 | \r | |
91e5e210 | 2266 | TString gVZEROSide = side;\r |
90bec8c5 | 2267 | for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {\r |
91e5e210 | 2268 | Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));\r |
2269 | //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;\r | |
2270 | if(gRunNumber == run) {\r | |
2271 | if(gVZEROSide == "A") \r | |
2272 | return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);\r | |
2273 | else if(gVZEROSide == "C") \r | |
2274 | return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);\r | |
2275 | }\r | |
2276 | }\r | |
2277 | \r | |
2278 | return 1.0;\r | |
2279 | }\r | |
0879e280 | 2280 | \r |
2281 | //________________________________________________________________________\r | |
2282 | void AliAnalysisTaskBFPsi::FinishTaskOutput(){\r | |
2283 | //Printf("END BF");\r | |
2284 | \r | |
2285 | if (!fBalance) {\r | |
f06d59b3 | 2286 | AliError("fBalance not available");\r |
0879e280 | 2287 | return;\r |
2288 | } \r | |
2289 | if(fRunShuffling) {\r | |
2290 | if (!fShuffledBalance) {\r | |
f06d59b3 | 2291 | AliError("fShuffledBalance not available");\r |
0879e280 | 2292 | return;\r |
2293 | }\r | |
2294 | }\r | |
2295 | \r | |
2296 | }\r | |
2297 | \r | |
2298 | //________________________________________________________________________\r | |
2299 | void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r | |
2300 | // Draw result to the screen\r | |
2301 | // Called once at the end of the query\r | |
2302 | \r | |
2303 | // not implemented ...\r | |
2304 | \r | |
2305 | }\r | |
1cb2a06e | 2306 | \r |