]>
Commit | Line | Data |
---|---|---|
0f67a57c | 1 | /* $Id$ */ |
2 | ||
3 | #include "AlidNdEtaCorrectionTask.h" | |
4 | ||
5 | #include <TCanvas.h> | |
6 | #include <TChain.h> | |
7 | #include <TFile.h> | |
8 | #include <TH1F.h> | |
54b096ef | 9 | #include <TH2F.h> |
1c15d51a | 10 | #include <TH3F.h> |
54b096ef | 11 | #include <TProfile.h> |
0f67a57c | 12 | #include <TParticle.h> |
1c15d51a | 13 | #include <TParticlePDG.h> |
c180f65d | 14 | #include <TDatabasePDG.h> |
81be4ee8 | 15 | #include <TRandom.h> |
0f67a57c | 16 | |
17 | #include <AliLog.h> | |
18 | #include <AliESDVertex.h> | |
19 | #include <AliESDEvent.h> | |
81be4ee8 | 20 | #include <AliESDRun.h> |
0f67a57c | 21 | #include <AliStack.h> |
22 | #include <AliHeader.h> | |
23 | #include <AliGenEventHeader.h> | |
24 | #include <AliMultiplicity.h> | |
25 | #include <AliAnalysisManager.h> | |
26 | #include <AliMCEventHandler.h> | |
27 | #include <AliMCEvent.h> | |
28 | #include <AliESDInputHandler.h> | |
29 | ||
745d6088 | 30 | #include "AliESDtrackCuts.h" |
0f67a57c | 31 | #include "AliPWG0Helper.h" |
0f67a57c | 32 | #include "dNdEta/dNdEtaAnalysis.h" |
33 | #include "dNdEta/AlidNdEtaCorrection.h" | |
70fdd197 | 34 | #include "AliTriggerAnalysis.h" |
81be4ee8 | 35 | #include "AliPhysicsSelection.h" |
0f67a57c | 36 | |
37 | ClassImp(AlidNdEtaCorrectionTask) | |
38 | ||
1c15d51a | 39 | AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() : |
40 | AliAnalysisTask(), | |
41 | fESD(0), | |
42 | fOutput(0), | |
43 | fOption(), | |
a7f69e56 | 44 | fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)), |
70fdd197 | 45 | fTrigger(AliTriggerAnalysis::kMB1), |
69b09e3b | 46 | fFillPhi(kFALSE), |
47 | fDeltaPhiCut(-1), | |
1d107532 | 48 | fSymmetrize(kFALSE), |
81be4ee8 | 49 | fMultAxisEta1(kFALSE), |
50 | fDiffTreatment(AliPWG0Helper::kMCFlags), | |
1c15d51a | 51 | fSignMode(0), |
52 | fOnlyPrimaries(kFALSE), | |
69b09e3b | 53 | fStatError(0), |
81be4ee8 | 54 | fSystSkipParticles(kFALSE), |
1c15d51a | 55 | fEsdTrackCuts(0), |
56 | fdNdEtaCorrection(0), | |
57 | fdNdEtaAnalysisMC(0), | |
58 | fdNdEtaAnalysisESD(0), | |
59 | fPIDParticles(0), | |
60 | fPIDTracks(0), | |
61 | fVertexCorrelation(0), | |
51f6de65 | 62 | fVertexCorrelationShift(0), |
1c15d51a | 63 | fVertexProfile(0), |
64 | fVertexShift(0), | |
65 | fVertexShiftNorm(0), | |
66 | fEtaCorrelation(0), | |
51f6de65 | 67 | fEtaCorrelationShift(0), |
1c15d51a | 68 | fEtaProfile(0), |
69 | fEtaResolution(0), | |
69b09e3b | 70 | fDeltaPhiCorrelation(0), |
1c15d51a | 71 | fpTResolution(0), |
72 | fEsdTrackCutsPrim(0), | |
73 | fEsdTrackCutsSec(0), | |
74 | fTemp1(0), | |
75 | fTemp2(0), | |
76 | fMultAll(0), | |
77 | fMultTr(0), | |
78 | fMultVtx(0), | |
79 | fEventStats(0) | |
80 | { | |
81 | // | |
82 | // Constructor. Initialization of pointers | |
83 | // | |
84 | ||
69b09e3b | 85 | for (Int_t i=0; i<4; i++) |
86 | fdNdEtaCorrectionSpecial[i] = 0; | |
1c15d51a | 87 | |
88 | for (Int_t i=0; i<8; i++) | |
89 | fDeltaPhi[i] = 0; | |
81be4ee8 | 90 | |
91 | AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning); | |
1c15d51a | 92 | } |
93 | ||
0f67a57c | 94 | AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) : |
95 | AliAnalysisTask("AlidNdEtaCorrectionTask", ""), | |
96 | fESD(0), | |
97 | fOutput(0), | |
98 | fOption(opt), | |
a7f69e56 | 99 | fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)), |
70fdd197 | 100 | fTrigger(AliTriggerAnalysis::kMB1), |
69b09e3b | 101 | fFillPhi(kFALSE), |
102 | fDeltaPhiCut(0), | |
1d107532 | 103 | fSymmetrize(kFALSE), |
81be4ee8 | 104 | fMultAxisEta1(kFALSE), |
105 | fDiffTreatment(AliPWG0Helper::kMCFlags), | |
0f67a57c | 106 | fSignMode(0), |
3d7758c1 | 107 | fOnlyPrimaries(kFALSE), |
69b09e3b | 108 | fStatError(0), |
81be4ee8 | 109 | fSystSkipParticles(kFALSE), |
0f67a57c | 110 | fEsdTrackCuts(0), |
111 | fdNdEtaCorrection(0), | |
112 | fdNdEtaAnalysisMC(0), | |
113 | fdNdEtaAnalysisESD(0), | |
114 | fPIDParticles(0), | |
96fcc8a7 | 115 | fPIDTracks(0), |
54b096ef | 116 | fVertexCorrelation(0), |
51f6de65 | 117 | fVertexCorrelationShift(0), |
54b096ef | 118 | fVertexProfile(0), |
1c15d51a | 119 | fVertexShift(0), |
54b096ef | 120 | fVertexShiftNorm(0), |
3d7758c1 | 121 | fEtaCorrelation(0), |
51f6de65 | 122 | fEtaCorrelationShift(0), |
3d7758c1 | 123 | fEtaProfile(0), |
1c15d51a | 124 | fEtaResolution(0), |
69b09e3b | 125 | fDeltaPhiCorrelation(0), |
1c15d51a | 126 | fpTResolution(0), |
127 | fEsdTrackCutsPrim(0), | |
128 | fEsdTrackCutsSec(0), | |
129 | fTemp1(0), | |
130 | fTemp2(0), | |
4ef701a0 | 131 | fMultAll(0), |
132 | fMultTr(0), | |
50ec344d | 133 | fMultVtx(0), |
134 | fEventStats(0) | |
0f67a57c | 135 | { |
136 | // | |
137 | // Constructor. Initialization of pointers | |
138 | // | |
139 | ||
140 | // Define input and output slots here | |
141 | DefineInput(0, TChain::Class()); | |
142 | DefineOutput(0, TList::Class()); | |
96fcc8a7 | 143 | |
69b09e3b | 144 | for (Int_t i=0; i<4; i++) |
145 | fdNdEtaCorrectionSpecial[i] = 0; | |
4ef701a0 | 146 | |
3d7758c1 | 147 | for (Int_t i=0; i<8; i++) |
4ef701a0 | 148 | fDeltaPhi[i] = 0; |
81be4ee8 | 149 | |
150 | AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning); | |
0f67a57c | 151 | } |
152 | ||
153 | AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask() | |
154 | { | |
155 | // | |
156 | // Destructor | |
157 | // | |
158 | ||
159 | // histograms are in the output list and deleted when the output | |
160 | // list is deleted by the TSelector dtor | |
161 | ||
162 | if (fOutput) { | |
163 | delete fOutput; | |
164 | fOutput = 0; | |
165 | } | |
166 | } | |
167 | ||
168 | //________________________________________________________________________ | |
169 | void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *) | |
170 | { | |
171 | // Connect ESD | |
172 | // Called once | |
173 | ||
174 | Printf("AlidNdEtaCorrectionTask::ConnectInputData called"); | |
175 | ||
c17301f3 | 176 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
177 | ||
178 | if (!esdH) { | |
179 | Printf("ERROR: Could not get ESDInputHandler"); | |
0f67a57c | 180 | } else { |
c17301f3 | 181 | fESD = esdH->GetEvent(); |
0f67a57c | 182 | |
c17301f3 | 183 | // Enable only the needed branches |
184 | esdH->SetActiveBranches("AliESDHeader Vertex"); | |
0f67a57c | 185 | |
a7f69e56 | 186 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
c17301f3 | 187 | esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity"); |
0f67a57c | 188 | |
a7f69e56 | 189 | if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) { |
c17301f3 | 190 | esdH->SetActiveBranches("AliESDHeader Vertex Tracks"); |
0f67a57c | 191 | } |
0f67a57c | 192 | } |
96fcc8a7 | 193 | |
194 | // disable info messages of AliMCEvent (per event) | |
195 | AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1); | |
0f67a57c | 196 | } |
197 | ||
198 | void AlidNdEtaCorrectionTask::CreateOutputObjects() | |
199 | { | |
200 | // create result objects and add to output list | |
201 | ||
202 | if (fOption.Contains("only-positive")) | |
203 | { | |
204 | Printf("INFO: Processing only positive particles."); | |
205 | fSignMode = 1; | |
206 | } | |
207 | else if (fOption.Contains("only-negative")) | |
208 | { | |
209 | Printf("INFO: Processing only negative particles."); | |
210 | fSignMode = -1; | |
211 | } | |
69b09e3b | 212 | |
213 | if (fOption.Contains("stat-error-1")) | |
214 | { | |
215 | Printf("INFO: Evaluation statistical errors. Mode: 1."); | |
216 | fStatError = 1; | |
217 | } | |
218 | else if (fOption.Contains("stat-error-2")) | |
219 | { | |
220 | Printf("INFO: Evaluation statistical errors. Mode: 2."); | |
221 | fStatError = 2; | |
222 | } | |
0f67a57c | 223 | |
0f67a57c | 224 | fOutput = new TList; |
225 | fOutput->SetOwner(); | |
226 | ||
770a1f1d | 227 | fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode); |
0f67a57c | 228 | fOutput->Add(fdNdEtaCorrection); |
229 | ||
1c15d51a | 230 | fPIDParticles = new TH1F("fPIDParticles", "PID of generated primary particles", 10001, -5000.5, 5000.5); |
0f67a57c | 231 | fOutput->Add(fPIDParticles); |
232 | ||
1c15d51a | 233 | fPIDTracks = new TH1F("fPIDTracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5); |
0f67a57c | 234 | fOutput->Add(fPIDTracks); |
235 | ||
770a1f1d | 236 | fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode); |
0f67a57c | 237 | fOutput->Add(fdNdEtaAnalysisMC); |
238 | ||
770a1f1d | 239 | fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode); |
0f67a57c | 240 | fOutput->Add(fdNdEtaAnalysisESD); |
96fcc8a7 | 241 | |
1c15d51a | 242 | if (fEsdTrackCuts) |
243 | { | |
244 | fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim")); | |
245 | fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec")); | |
246 | fOutput->Add(fEsdTrackCutsPrim); | |
247 | fOutput->Add(fEsdTrackCutsSec); | |
248 | } | |
249 | ||
96fcc8a7 | 250 | if (fOption.Contains("process-types")) { |
69b09e3b | 251 | fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode); |
252 | fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode); | |
253 | fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode); | |
96fcc8a7 | 254 | |
69b09e3b | 255 | fOutput->Add(fdNdEtaCorrectionSpecial[0]); |
256 | fOutput->Add(fdNdEtaCorrectionSpecial[1]); | |
257 | fOutput->Add(fdNdEtaCorrectionSpecial[2]); | |
258 | } | |
259 | ||
260 | if (fOption.Contains("particle-species")) { | |
261 | fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_pi", "dndeta_correction_pi", fAnalysisMode); | |
262 | fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_K", "dndeta_correction_K", fAnalysisMode); | |
263 | fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_p", "dndeta_correction_p", fAnalysisMode); | |
264 | fdNdEtaCorrectionSpecial[3] = new AlidNdEtaCorrection("dndeta_correction_other", "dndeta_correction_other", fAnalysisMode); | |
265 | ||
266 | for (Int_t i=0; i<4; i++) | |
267 | fOutput->Add(fdNdEtaCorrectionSpecial[i]); | |
96fcc8a7 | 268 | } |
269 | ||
69b09e3b | 270 | |
81be4ee8 | 271 | //fTemp1 = new TH2F("fTemp1", "fTemp1", 4, 0.5, 4.5, 101, -1.5, 99.5); // nsd study |
272 | fTemp1 = new TH2F("fTemp1", "fTemp1", 300, -15, 15, 80, -2.0, 2.0); | |
1c15d51a | 273 | fOutput->Add(fTemp1); |
81be4ee8 | 274 | |
275 | fTemp2 = new TH2F("fTemp2", "fTemp2", 300, -15, 15, 80, -2.0, 2.0); | |
1c15d51a | 276 | fOutput->Add(fTemp2); |
54b096ef | 277 | |
0fc41645 | 278 | fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30); |
279 | fOutput->Add(fVertexCorrelation); | |
a7f69e56 | 280 | fVertexCorrelationShift = new TH3F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx;rec. tracks", 120, -30, 30, 100, -1, 1, 100, -0.5, 99.5); |
51f6de65 | 281 | fOutput->Add(fVertexCorrelationShift); |
54b096ef | 282 | fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20); |
0fc41645 | 283 | fOutput->Add(fVertexProfile); |
1c15d51a | 284 | fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2); |
285 | fOutput->Add(fVertexShift); | |
81be4ee8 | 286 | fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx);rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5); |
0fc41645 | 287 | fOutput->Add(fVertexShiftNorm); |
288 | ||
3d7758c1 | 289 | fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3); |
0fc41645 | 290 | fOutput->Add(fEtaCorrelation); |
69b09e3b | 291 | fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1); |
51f6de65 | 292 | fOutput->Add(fEtaCorrelationShift); |
3d7758c1 | 293 | fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3); |
0fc41645 | 294 | fOutput->Add(fEtaProfile); |
1c15d51a | 295 | fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2); |
296 | fOutput->Add(fEtaResolution); | |
297 | ||
69b09e3b | 298 | fpTResolution = new TH2F("fpTResolution", ";MC p_{T} (GeV/c);(MC p_{T} - ESD p_{T}) / MC p_{T}", 160, 0, 20, 201, -0.2, 0.2); |
1c15d51a | 299 | fOutput->Add(fpTResolution); |
3d7758c1 | 300 | |
4ef701a0 | 301 | fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5); |
0fc41645 | 302 | fOutput->Add(fMultAll); |
4ef701a0 | 303 | fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5); |
0fc41645 | 304 | fOutput->Add(fMultVtx); |
4ef701a0 | 305 | fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5); |
0fc41645 | 306 | fOutput->Add(fMultTr); |
4ef701a0 | 307 | |
3d7758c1 | 308 | for (Int_t i=0; i<8; i++) |
0fc41645 | 309 | { |
69b09e3b | 310 | fDeltaPhi[i] = new TH2F(Form("fDeltaPhi_%d", i), ";#Delta phi;#phi;Entries", 2000, -0.1, 0.1, 18 * 5, 0, TMath::TwoPi()); |
0fc41645 | 311 | fOutput->Add(fDeltaPhi[i]); |
312 | } | |
50ec344d | 313 | |
81be4ee8 | 314 | fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 109, -6.5, 102.5, 4, -0.5, 3.5); |
0fc41645 | 315 | fOutput->Add(fEventStats); |
69b09e3b | 316 | fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -5 |
317 | fEventStats->GetXaxis()->SetBinLabel(2, "NSD"); // x = -4 | |
318 | fEventStats->GetXaxis()->SetBinLabel(3, "ND"); // x = -3 | |
319 | fEventStats->GetXaxis()->SetBinLabel(4, "SD"); // x = -2 | |
320 | fEventStats->GetXaxis()->SetBinLabel(5, "DD"); // x = -1 | |
81be4ee8 | 321 | |
322 | fEventStats->GetXaxis()->SetBinLabel(108, "INEL=0"); // x = -101 | |
323 | fEventStats->GetXaxis()->SetBinLabel(109, "INEL>0"); // x = -102 | |
69b09e3b | 324 | |
81be4ee8 | 325 | for (Int_t i=-1; i<100; i++) |
69b09e3b | 326 | fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i)); |
50ec344d | 327 | |
328 | fEventStats->GetYaxis()->SetBinLabel(1, "nothing"); | |
329 | fEventStats->GetYaxis()->SetBinLabel(2, "trg"); | |
330 | fEventStats->GetYaxis()->SetBinLabel(3, "vtx"); | |
331 | fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx"); | |
1c15d51a | 332 | |
333 | if (fEsdTrackCuts) | |
334 | { | |
335 | fEsdTrackCuts->SetName("fEsdTrackCuts"); | |
81be4ee8 | 336 | // TODO like this we send an empty object back... |
337 | fOutput->Add(fEsdTrackCuts->Clone()); | |
1c15d51a | 338 | } |
0f67a57c | 339 | } |
340 | ||
341 | void AlidNdEtaCorrectionTask::Exec(Option_t*) | |
342 | { | |
343 | // process the event | |
344 | ||
345 | // Check prerequisites | |
346 | if (!fESD) | |
347 | { | |
348 | AliDebug(AliLog::kError, "ESD branch not available"); | |
349 | return; | |
350 | } | |
351 | ||
3d7758c1 | 352 | if (fOnlyPrimaries) |
353 | Printf("WARNING: Processing only primaries. For systematical studies only!"); | |
69b09e3b | 354 | |
355 | if (fStatError > 0) | |
356 | Printf("WARNING: Statistical error evaluation active!"); | |
a7f69e56 | 357 | |
81be4ee8 | 358 | AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
359 | if (!inputHandler) | |
360 | { | |
361 | Printf("ERROR: Could not receive input handler"); | |
362 | return; | |
363 | } | |
364 | ||
365 | Bool_t eventTriggered = inputHandler->IsEventSelected(); | |
54b096ef | 366 | |
81be4ee8 | 367 | static AliTriggerAnalysis* triggerAnalysis = 0; |
368 | if (!triggerAnalysis) | |
369 | { | |
370 | AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection()); | |
371 | if (physicsSelection) | |
372 | triggerAnalysis = physicsSelection->GetTriggerAnalysis(); | |
373 | } | |
374 | ||
375 | if (eventTriggered) | |
376 | eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger); | |
377 | ||
3d7758c1 | 378 | if (!eventTriggered) |
379 | Printf("No trigger"); | |
380 | ||
54b096ef | 381 | // post the data already here |
382 | PostData(0, fOutput); | |
69b09e3b | 383 | |
54b096ef | 384 | // MC info |
385 | AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
386 | if (!eventHandler) { | |
387 | Printf("ERROR: Could not retrieve MC event handler"); | |
388 | return; | |
389 | } | |
390 | ||
391 | AliMCEvent* mcEvent = eventHandler->MCEvent(); | |
392 | if (!mcEvent) { | |
393 | Printf("ERROR: Could not retrieve MC event"); | |
394 | return; | |
395 | } | |
396 | ||
397 | AliStack* stack = mcEvent->Stack(); | |
398 | if (!stack) | |
399 | { | |
400 | AliDebug(AliLog::kError, "Stack not available"); | |
401 | return; | |
402 | } | |
403 | ||
404 | AliHeader* header = mcEvent->Header(); | |
405 | if (!header) | |
406 | { | |
407 | AliDebug(AliLog::kError, "Header not available"); | |
408 | return; | |
409 | } | |
410 | ||
81be4ee8 | 411 | // get process type |
412 | Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment); | |
413 | ||
6b62a9c7 | 414 | AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType)); |
54b096ef | 415 | |
6b62a9c7 | 416 | if (processType == AliPWG0Helper::kInvalidProcess) |
81be4ee8 | 417 | { |
418 | AliDebug(AliLog::kWarning, "Unknown process type. Setting to ND"); | |
419 | processType = AliPWG0Helper::kND; | |
420 | } | |
54b096ef | 421 | |
422 | // get the MC vertex | |
423 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
424 | TArrayF vtxMC(3); | |
425 | genHeader->PrimaryVertex(vtxMC); | |
426 | ||
427 | // get the ESD vertex | |
770a1f1d | 428 | Bool_t eventVertex = kFALSE; |
a7f69e56 | 429 | Double_t vtx[3]; |
430 | const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode); | |
431 | if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode)) | |
54b096ef | 432 | { |
433 | vtxESD->GetXYZ(vtx); | |
81be4ee8 | 434 | eventVertex = kTRUE; |
435 | ||
436 | // remove vertices outside +- 15 cm | |
437 | if (TMath::Abs(vtx[2]) > 15) | |
438 | { | |
439 | eventVertex = kFALSE; | |
440 | vtxESD = 0; | |
441 | } | |
567160d6 | 442 | } |
a7f69e56 | 443 | else |
444 | vtxESD = 0; | |
445 | ||
0f67a57c | 446 | // create list of (label, eta, pt) tuples |
447 | Int_t inputCount = 0; | |
448 | Int_t* labelArr = 0; | |
3d7758c1 | 449 | Int_t* labelArr2 = 0; // only for case of SPD |
0f67a57c | 450 | Float_t* etaArr = 0; |
69b09e3b | 451 | Float_t* thirdDimArr = 0; |
4ef701a0 | 452 | Float_t* deltaPhiArr = 0; |
a7f69e56 | 453 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
0f67a57c | 454 | { |
81be4ee8 | 455 | if (vtxESD) |
0f67a57c | 456 | { |
81be4ee8 | 457 | // get tracklets |
458 | const AliMultiplicity* mult = fESD->GetMultiplicity(); | |
459 | if (!mult) | |
460 | { | |
461 | AliDebug(AliLog::kError, "AliMultiplicity not available"); | |
462 | return; | |
463 | } | |
464 | ||
465 | labelArr = new Int_t[mult->GetNumberOfTracklets()]; | |
466 | labelArr2 = new Int_t[mult->GetNumberOfTracklets()]; | |
467 | etaArr = new Float_t[mult->GetNumberOfTracklets()]; | |
468 | thirdDimArr = new Float_t[mult->GetNumberOfTracklets()]; | |
469 | deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()]; | |
470 | ||
471 | Bool_t foundInEta10 = kFALSE; | |
472 | ||
473 | // get multiplicity from SPD tracklets | |
474 | for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i) | |
475 | { | |
476 | //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i)); | |
477 | ||
478 | Float_t phi = mult->GetPhi(i); | |
479 | if (phi < 0) | |
480 | phi += TMath::Pi() * 2; | |
481 | Float_t deltaPhi = mult->GetDeltaPhi(i); | |
482 | ||
483 | if (TMath::Abs(deltaPhi) > 1) | |
484 | printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi); | |
485 | ||
486 | if (fOnlyPrimaries) | |
487 | if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0))) | |
488 | continue; | |
489 | ||
490 | if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025)) | |
3d7758c1 | 491 | continue; |
81be4ee8 | 492 | |
493 | if (fSystSkipParticles && gRandom->Uniform() < 0.0153) | |
494 | { | |
495 | Printf("Skipped tracklet!"); | |
496 | continue; | |
497 | } | |
498 | ||
499 | // TEST exclude potentially inefficient phi region | |
500 | //if (phi > 5.70 || phi < 0.06) | |
501 | // continue; | |
502 | ||
503 | // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles | |
504 | if (TMath::Abs(mult->GetEta(i)) < 1) | |
505 | foundInEta10 = kTRUE; | |
506 | ||
507 | etaArr[inputCount] = mult->GetEta(i); | |
508 | if (fSymmetrize) | |
509 | etaArr[inputCount] = TMath::Abs(etaArr[inputCount]); | |
510 | labelArr[inputCount] = mult->GetLabel(i, 0); | |
511 | labelArr2[inputCount] = mult->GetLabel(i, 1); | |
512 | thirdDimArr[inputCount] = phi; | |
513 | deltaPhiArr[inputCount] = deltaPhi; | |
514 | ++inputCount; | |
515 | } | |
69b09e3b | 516 | |
81be4ee8 | 517 | /* |
518 | for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i) | |
519 | { | |
520 | if (TMath::Abs(TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.))) < 1); | |
521 | { | |
522 | foundInEta10 = kTRUE; | |
523 | break; | |
524 | } | |
525 | } | |
526 | */ | |
527 | ||
528 | if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10) | |
529 | eventTriggered = kFALSE; | |
0f67a57c | 530 | } |
531 | } | |
a7f69e56 | 532 | else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) |
0f67a57c | 533 | { |
534 | if (!fEsdTrackCuts) | |
535 | { | |
536 | AliDebug(AliLog::kError, "fESDTrackCuts not available"); | |
537 | return; | |
538 | } | |
a7f69e56 | 539 | |
81be4ee8 | 540 | Bool_t foundInEta10 = kFALSE; |
541 | ||
69b09e3b | 542 | if (vtxESD) |
1c15d51a | 543 | { |
69b09e3b | 544 | // get multiplicity from ESD tracks |
a7f69e56 | 545 | TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC)); |
69b09e3b | 546 | Int_t nGoodTracks = list->GetEntries(); |
547 | ||
548 | Printf("Accepted %d tracks", nGoodTracks); | |
549 | ||
550 | labelArr = new Int_t[nGoodTracks]; | |
551 | labelArr2 = new Int_t[nGoodTracks]; | |
552 | etaArr = new Float_t[nGoodTracks]; | |
553 | thirdDimArr = new Float_t[nGoodTracks]; | |
554 | deltaPhiArr = new Float_t[nGoodTracks]; | |
555 | ||
556 | // loop over esd tracks | |
557 | for (Int_t i=0; i<nGoodTracks; i++) | |
1c15d51a | 558 | { |
69b09e3b | 559 | AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i)); |
560 | if (!esdTrack) | |
1c15d51a | 561 | { |
69b09e3b | 562 | AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i)); |
1c15d51a | 563 | continue; |
564 | } | |
69b09e3b | 565 | |
81be4ee8 | 566 | if (esdTrack->Pt() < 0.15) |
567 | continue; | |
a7f69e56 | 568 | |
569 | if (fOnlyPrimaries) | |
570 | { | |
571 | Int_t label = TMath::Abs(esdTrack->GetLabel()); | |
572 | if (label == 0) | |
573 | continue; | |
574 | ||
575 | if (stack->IsPhysicalPrimary(label) == kFALSE) | |
576 | continue; | |
81be4ee8 | 577 | } |
578 | ||
579 | // INEL>0 trigger | |
580 | if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15) | |
581 | foundInEta10 = kTRUE; | |
69b09e3b | 582 | |
583 | etaArr[inputCount] = esdTrack->Eta(); | |
1d107532 | 584 | if (fSymmetrize) |
585 | etaArr[inputCount] = TMath::Abs(etaArr[inputCount]); | |
69b09e3b | 586 | labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel()); |
587 | labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks | |
588 | thirdDimArr[inputCount] = esdTrack->Pt(); | |
589 | deltaPhiArr[inputCount] = 0; // no delta phi for tracks | |
590 | ++inputCount; | |
591 | } | |
1c15d51a | 592 | |
69b09e3b | 593 | delete list; |
1c15d51a | 594 | |
81be4ee8 | 595 | // TODO this code crashes for TPCITS because particles are requested from the stack for some labels that are out of bound |
596 | if (0 && eventTriggered) | |
69b09e3b | 597 | { |
598 | // collect values for primaries and secondaries | |
599 | for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++) | |
1c15d51a | 600 | { |
69b09e3b | 601 | AliESDtrack* track = 0; |
602 | ||
a7f69e56 | 603 | if (fAnalysisMode & AliPWG0Helper::kTPC) |
69b09e3b | 604 | track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack); |
a7f69e56 | 605 | else if (fAnalysisMode & AliPWG0Helper::kTPCITS) |
69b09e3b | 606 | track = fESD->GetTrack(iTrack); |
607 | ||
608 | if (!track) | |
609 | continue; | |
610 | ||
611 | Int_t label = TMath::Abs(track->GetLabel()); | |
612 | if (!stack->Particle(label) || !stack->Particle(label)->GetPDG()) | |
613 | { | |
614 | Printf("WARNING: No particle for %d", label); | |
615 | if (stack->Particle(label)) | |
616 | stack->Particle(label)->Print(); | |
617 | continue; | |
618 | } | |
619 | ||
620 | if (stack->Particle(label)->GetPDG()->Charge() == 0) | |
621 | continue; | |
622 | ||
a7f69e56 | 623 | if (TMath::Abs(track->Eta()) < 0.8 && track->Pt() > 0.15) |
69b09e3b | 624 | { |
625 | if (stack->IsPhysicalPrimary(label)) | |
1c15d51a | 626 | { |
69b09e3b | 627 | // primary |
628 | if (fEsdTrackCutsPrim->AcceptTrack(track)) | |
629 | { | |
a7f69e56 | 630 | // if (AliESDtrackCuts::GetSigmaToVertex(track) > 900) |
631 | // { | |
632 | // Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track)); | |
633 | // Float_t b[2]; | |
634 | // Float_t r[3]; | |
635 | // track->GetImpactParameters(b, r); | |
636 | // Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]); | |
637 | // track->Print(""); | |
638 | // if (vtxESD) | |
639 | // vtxESD->Print(); | |
640 | // } | |
69b09e3b | 641 | } |
1c15d51a | 642 | } |
69b09e3b | 643 | else |
644 | { | |
645 | // secondary | |
646 | fEsdTrackCutsSec->AcceptTrack(track); | |
1c15d51a | 647 | } |
648 | } | |
69b09e3b | 649 | |
650 | // TODO mem leak in the continue statements above | |
a7f69e56 | 651 | if (fAnalysisMode & AliPWG0Helper::kTPC) |
69b09e3b | 652 | delete track; |
1c15d51a | 653 | } |
1c15d51a | 654 | } |
655 | } | |
81be4ee8 | 656 | |
657 | if (!foundInEta10) | |
658 | eventTriggered = kFALSE; | |
0f67a57c | 659 | } |
54b096ef | 660 | else |
0f67a57c | 661 | return; |
0f67a57c | 662 | |
81be4ee8 | 663 | // fill process type |
664 | Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex; | |
665 | // INEL | |
666 | fEventStats->Fill(-6, biny); | |
667 | // NSD | |
668 | if (processType != AliPWG0Helper::kSD) | |
669 | fEventStats->Fill(-5, biny); | |
670 | // SD, ND, DD | |
671 | if (processType == AliPWG0Helper::kND) | |
672 | fEventStats->Fill(-4, biny); | |
673 | if (processType == AliPWG0Helper::kSD) | |
674 | fEventStats->Fill(-3, biny); | |
675 | if (processType == AliPWG0Helper::kDD) | |
676 | fEventStats->Fill(-2, biny); | |
677 | ||
0f67a57c | 678 | // loop over mc particles |
679 | Int_t nPrim = stack->GetNprimary(); | |
4ef701a0 | 680 | Int_t nAccepted = 0; |
0f67a57c | 681 | |
81be4ee8 | 682 | Bool_t oneParticleEvent = kFALSE; |
683 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) | |
684 | { | |
685 | //Printf("Getting particle %d", iMc); | |
686 | TParticle* particle = stack->Particle(iMc); | |
687 | ||
688 | if (!particle) | |
689 | continue; | |
690 | ||
691 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
692 | continue; | |
693 | ||
694 | if (TMath::Abs(particle->Eta()) < 1.0) | |
695 | { | |
696 | oneParticleEvent = kTRUE; | |
697 | break; | |
698 | } | |
699 | } | |
700 | ||
701 | if (TMath::Abs(vtxMC[2]) < 5.5) | |
702 | { | |
703 | if (oneParticleEvent) | |
704 | fEventStats->Fill(102, biny); | |
705 | else | |
706 | fEventStats->Fill(101, biny); | |
707 | } | |
708 | ||
0f67a57c | 709 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) |
710 | { | |
711 | //Printf("Getting particle %d", iMc); | |
712 | TParticle* particle = stack->Particle(iMc); | |
713 | ||
714 | if (!particle) | |
715 | continue; | |
716 | ||
717 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
1c15d51a | 718 | { |
719 | //if (TMath::Abs(particle->GetPdgCode()) > 3000 && TMath::Abs(particle->Eta()) < 1.0) | |
720 | // fPIDParticles->Fill(particle->GetPdgCode()); | |
0f67a57c | 721 | continue; |
1c15d51a | 722 | } |
0f67a57c | 723 | |
724 | if (SignOK(particle->GetPDG()) == kFALSE) | |
725 | continue; | |
81be4ee8 | 726 | |
1c15d51a | 727 | if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0) |
728 | fPIDParticles->Fill(particle->GetPdgCode()); | |
729 | ||
0f67a57c | 730 | Float_t eta = particle->Eta(); |
1d107532 | 731 | if (fSymmetrize) |
732 | eta = TMath::Abs(eta); | |
69b09e3b | 733 | |
734 | Float_t thirdDim = -1; | |
a7f69e56 | 735 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
69b09e3b | 736 | { |
737 | if (fFillPhi) | |
738 | { | |
739 | thirdDim = particle->Phi(); | |
740 | } | |
741 | else | |
742 | thirdDim = inputCount; | |
743 | } | |
744 | else | |
745 | thirdDim = particle->Pt(); | |
746 | ||
747 | // calculate y | |
748 | //Float_t y = 0.5 * TMath::Log((particle->Energy() + particle->Pz()) / (particle->Energy() - particle->Pz())); | |
749 | //fTemp1->Fill(eta); | |
750 | //fTemp2->Fill(y); | |
0f67a57c | 751 | |
81be4ee8 | 752 | Int_t processType2 = processType; |
753 | if (oneParticleEvent) | |
754 | processType2 |= AliPWG0Helper::kOnePart; | |
755 | fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); | |
0f67a57c | 756 | |
69b09e3b | 757 | if (fOption.Contains("process-types")) |
96fcc8a7 | 758 | { |
759 | // non diffractive | |
6b62a9c7 | 760 | if (processType==AliPWG0Helper::kND) |
81be4ee8 | 761 | fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 762 | |
763 | // single diffractive | |
6b62a9c7 | 764 | if (processType==AliPWG0Helper::kSD) |
81be4ee8 | 765 | fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 766 | |
767 | // double diffractive | |
6b62a9c7 | 768 | if (processType==AliPWG0Helper::kDD) |
81be4ee8 | 769 | fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
69b09e3b | 770 | } |
771 | ||
772 | if (fOption.Contains("particle-species")) | |
773 | { | |
774 | Int_t id = -1; | |
775 | switch (TMath::Abs(particle->GetPdgCode())) | |
776 | { | |
777 | case 211: id = 0; break; | |
778 | case 321: id = 1; break; | |
779 | case 2212: id = 2; break; | |
780 | default: id = 3; break; | |
781 | } | |
81be4ee8 | 782 | fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 783 | } |
784 | ||
0f67a57c | 785 | if (eventTriggered) |
786 | if (eventVertex) | |
51f6de65 | 787 | fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, thirdDim); |
4ef701a0 | 788 | |
0fc41645 | 789 | // TODO this value might be needed lower for the SPD study (only used for control histograms anyway) |
ea441adf | 790 | if (TMath::Abs(eta) < 1.5) // && pt > 0.2) |
4ef701a0 | 791 | nAccepted++; |
792 | } | |
793 | ||
81be4ee8 | 794 | if (AliPWG0Helper::GetLastProcessType() >= -1) |
795 | fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny); | |
ea441adf | 796 | |
4ef701a0 | 797 | fMultAll->Fill(nAccepted); |
798 | if (eventTriggered) { | |
799 | fMultTr->Fill(nAccepted); | |
800 | if (eventVertex) | |
801 | fMultVtx->Fill(nAccepted); | |
0f67a57c | 802 | } |
803 | ||
3d7758c1 | 804 | Int_t processed = 0; |
69b09e3b | 805 | |
806 | Bool_t* primCount = 0; | |
807 | if (fStatError > 0) | |
808 | { | |
809 | primCount = new Bool_t[nPrim]; | |
810 | for (Int_t i=0; i<nPrim; ++i) | |
811 | primCount[i] = kFALSE; | |
812 | } | |
3d7758c1 | 813 | |
70fdd197 | 814 | Int_t nEta05 = 0; |
81be4ee8 | 815 | Int_t nEta10 = 0; |
0f67a57c | 816 | for (Int_t i=0; i<inputCount; ++i) |
817 | { | |
70fdd197 | 818 | if (TMath::Abs(etaArr[i]) < 0.5) |
819 | nEta05++; | |
81be4ee8 | 820 | if (TMath::Abs(etaArr[i]) < 1.0) |
821 | nEta10++; | |
822 | } | |
70fdd197 | 823 | |
81be4ee8 | 824 | for (Int_t i=0; i<inputCount; ++i) |
825 | { | |
0f67a57c | 826 | Int_t label = labelArr[i]; |
3d7758c1 | 827 | Int_t label2 = labelArr2[i]; |
0f67a57c | 828 | |
829 | if (label < 0) | |
830 | { | |
0fc41645 | 831 | Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label); |
0f67a57c | 832 | continue; |
833 | } | |
834 | ||
835 | TParticle* particle = stack->Particle(label); | |
836 | if (!particle) | |
837 | { | |
838 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label)); | |
839 | continue; | |
840 | } | |
1c15d51a | 841 | |
81be4ee8 | 842 | if (TMath::Abs(particle->Eta()) < 1.0) |
843 | fPIDTracks->Fill(particle->GetPdgCode()); | |
3d7758c1 | 844 | |
845 | // find particle that is filled in the correction map | |
846 | // this should be the particle that has been reconstructed | |
847 | // for TPC this is clear and always identified by the track's label | |
848 | // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken: | |
849 | // L1 L2 (P = primary, S = secondary) | |
850 | // P P' : --> P | |
851 | // P S : --> P | |
852 | // S P : --> P | |
853 | // S S' : --> S | |
854 | ||
855 | if (label != label2 && label2 >= 0) | |
856 | { | |
857 | if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2)) | |
858 | { | |
859 | particle = stack->Particle(label2); | |
860 | if (!particle) | |
861 | { | |
862 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2)); | |
863 | continue; | |
864 | } | |
865 | } | |
866 | } | |
0f67a57c | 867 | |
868 | if (eventTriggered && eventVertex) | |
869 | { | |
3d7758c1 | 870 | if (SignOK(particle->GetPDG()) == kFALSE) |
871 | continue; | |
872 | ||
873 | processed++; | |
874 | ||
1c15d51a | 875 | // resolutions |
1d107532 | 876 | if (fSymmetrize) |
877 | fEtaResolution->Fill(TMath::Abs(particle->Eta()) - etaArr[i]); | |
878 | else | |
879 | fEtaResolution->Fill(particle->Eta() - etaArr[i]); | |
69b09e3b | 880 | |
a7f69e56 | 881 | if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) |
69b09e3b | 882 | if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0) |
883 | fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt()); | |
1c15d51a | 884 | |
51f6de65 | 885 | Float_t eta = -999; |
886 | Float_t thirdDim = -1; | |
887 | ||
3d7758c1 | 888 | Bool_t firstIsPrim = stack->IsPhysicalPrimary(label); |
a7f69e56 | 889 | // in case of same label the MC values are filled, otherwise (background) the reconstructed values |
890 | if (label == label2) | |
3d7758c1 | 891 | { |
51f6de65 | 892 | eta = particle->Eta(); |
1d107532 | 893 | if (fSymmetrize) |
894 | eta = TMath::Abs(eta); | |
69b09e3b | 895 | |
a7f69e56 | 896 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
69b09e3b | 897 | { |
898 | if (fFillPhi) | |
899 | { | |
900 | thirdDim = particle->Phi(); | |
901 | } | |
902 | else | |
903 | thirdDim = inputCount; | |
904 | } | |
905 | else | |
906 | thirdDim = particle->Pt(); | |
3d7758c1 | 907 | } |
908 | else | |
909 | { | |
a7f69e56 | 910 | if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi) |
69b09e3b | 911 | { |
81be4ee8 | 912 | thirdDim = (fMultAxisEta1) ? nEta10 : inputCount; |
69b09e3b | 913 | } |
914 | else | |
915 | thirdDim = thirdDimArr[i]; | |
916 | ||
51f6de65 | 917 | eta = etaArr[i]; |
3d7758c1 | 918 | } |
745d6088 | 919 | |
69b09e3b | 920 | Bool_t fillTrack = kTRUE; |
51f6de65 | 921 | |
69b09e3b | 922 | // statistical error evaluation active? |
923 | if (fStatError > 0) | |
924 | { | |
925 | Bool_t statErrorDecision = kFALSE; | |
926 | ||
927 | // primary and not yet counted | |
928 | if (label == label2 && firstIsPrim && !primCount[label]) | |
929 | { | |
930 | statErrorDecision = kTRUE; | |
931 | primCount[label] = kTRUE; | |
932 | } | |
933 | ||
934 | // in case of 1 we count only unique primaries, in case of 2 all the rest | |
935 | if (fStatError == 1) | |
936 | { | |
937 | fillTrack = statErrorDecision; | |
938 | } | |
939 | else if (fStatError == 2) | |
940 | fillTrack = !statErrorDecision; | |
941 | } | |
942 | ||
943 | if (fillTrack) | |
81be4ee8 | 944 | { |
69b09e3b | 945 | fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim); |
81be4ee8 | 946 | fTemp2->Fill(vtxMC[2], eta); |
947 | } | |
948 | ||
51f6de65 | 949 | // eta comparison for tracklets with the same label (others are background) |
950 | if (label == label2) | |
951 | { | |
1d107532 | 952 | Float_t eta2 = particle->Eta(); |
953 | if (fSymmetrize) | |
954 | eta2 = TMath::Abs(eta2); | |
955 | ||
956 | fEtaProfile->Fill(eta2, eta2 - etaArr[i]); | |
957 | fEtaCorrelation->Fill(etaArr[i], eta2); | |
958 | fEtaCorrelationShift->Fill(eta2, eta2 - etaArr[i]); | |
51f6de65 | 959 | } |
3d7758c1 | 960 | |
1d107532 | 961 | if (fSymmetrize) |
962 | fdNdEtaAnalysisESD->FillTrack(vtxMC[2], TMath::Abs(particle->Eta()), thirdDim); | |
963 | else | |
964 | fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim); | |
3d7758c1 | 965 | |
69b09e3b | 966 | if (fOption.Contains("process-types")) |
96fcc8a7 | 967 | { |
968 | // non diffractive | |
6b62a9c7 | 969 | if (processType == AliPWG0Helper::kND) |
69b09e3b | 970 | fdNdEtaCorrectionSpecial[0]->FillTrackedParticle(vtxMC[2], eta, thirdDim); |
96fcc8a7 | 971 | |
972 | // single diffractive | |
1c15d51a | 973 | if (processType == AliPWG0Helper::kSD) |
69b09e3b | 974 | fdNdEtaCorrectionSpecial[1]->FillTrackedParticle(vtxMC[2], eta, thirdDim); |
96fcc8a7 | 975 | |
976 | // double diffractive | |
6b62a9c7 | 977 | if (processType == AliPWG0Helper::kDD) |
69b09e3b | 978 | fdNdEtaCorrectionSpecial[2]->FillTrackedParticle(vtxMC[2], eta, thirdDim); |
979 | } | |
980 | ||
981 | if (fOption.Contains("particle-species")) | |
982 | { | |
983 | // find mother first | |
984 | TParticle* mother = AliPWG0Helper::FindPrimaryMother(stack, label); | |
985 | ||
986 | Int_t id = -1; | |
987 | switch (TMath::Abs(mother->GetPdgCode())) | |
988 | { | |
989 | case 211: id = 0; break; | |
990 | case 321: id = 1; break; | |
991 | case 2212: id = 2; break; | |
992 | default: id = 3; break; | |
993 | } | |
994 | fdNdEtaCorrectionSpecial[id]->FillTrackedParticle(vtxMC[2], eta, thirdDim); | |
96fcc8a7 | 995 | } |
4ef701a0 | 996 | |
745d6088 | 997 | // control histograms |
3d7758c1 | 998 | Int_t hist = -1; |
999 | if (label == label2) | |
4ef701a0 | 1000 | { |
69b09e3b | 1001 | if (firstIsPrim) |
3d7758c1 | 1002 | { |
69b09e3b | 1003 | hist = 0; |
3d7758c1 | 1004 | } |
1005 | else | |
1006 | hist = 1; | |
1007 | } | |
1008 | else if (label2 >= 0) | |
1009 | { | |
1010 | Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2); | |
1011 | if (firstIsPrim && secondIsPrim) | |
1012 | { | |
1013 | hist = 2; | |
1014 | } | |
1015 | else if (firstIsPrim && !secondIsPrim) | |
1016 | { | |
1017 | hist = 3; | |
69b09e3b | 1018 | |
3d7758c1 | 1019 | // check if secondary is caused by the primary or it is a fake combination |
1020 | //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2); | |
1021 | Int_t mother = label2; | |
1022 | while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0) | |
1023 | { | |
1024 | //Printf(" %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1)); | |
1025 | if (stack->Particle(mother)->GetMother(0) == label) | |
1026 | { | |
1027 | /*Printf("The untraceable decay was:"); | |
1028 | Printf(" from:"); | |
1029 | particle->Print(); | |
1030 | Printf(" to (status code %d):", stack->Particle(mother)->GetStatusCode()); | |
1031 | stack->Particle(mother)->Print();*/ | |
1032 | hist = 4; | |
1033 | } | |
1034 | mother = stack->Particle(mother)->GetMother(0); | |
1035 | } | |
1036 | } | |
1037 | else if (!firstIsPrim && secondIsPrim) | |
1038 | { | |
1039 | hist = 5; | |
1040 | } | |
1041 | else if (!firstIsPrim && !secondIsPrim) | |
1042 | { | |
1043 | hist = 6; | |
1044 | } | |
1045 | ||
4ef701a0 | 1046 | } |
1047 | else | |
3d7758c1 | 1048 | hist = 7; |
1049 | ||
69b09e3b | 1050 | fDeltaPhi[hist]->Fill(deltaPhiArr[i], thirdDimArr[i]); |
0f67a57c | 1051 | } |
1052 | } | |
69b09e3b | 1053 | |
1054 | if (primCount) | |
1055 | { | |
1056 | delete[] primCount; | |
1057 | primCount = 0; | |
1058 | } | |
0f67a57c | 1059 | |
3d7758c1 | 1060 | if (processed < inputCount) |
1061 | Printf("Only %d out of %d track(let)s were used", processed, inputCount); | |
1062 | ||
a7f69e56 | 1063 | if (eventTriggered && eventVertex) |
1064 | { | |
1065 | Double_t diff = vtxMC[2] - vtx[2]; | |
1066 | fVertexShift->Fill(diff); | |
81be4ee8 | 1067 | |
a7f69e56 | 1068 | fVertexCorrelation->Fill(vtxMC[2], vtx[2]); |
1069 | fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2], inputCount); | |
1070 | fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]); | |
81be4ee8 | 1071 | |
1072 | if (vtxESD->IsFromVertexerZ() && inputCount > 0) | |
1073 | fVertexShiftNorm->Fill(diff, vtxESD->GetNContributors()); | |
a7f69e56 | 1074 | } |
1075 | ||
81be4ee8 | 1076 | Int_t multAxis = inputCount; |
1077 | if (fMultAxisEta1) | |
1078 | multAxis = nEta10; | |
1079 | ||
0f67a57c | 1080 | if (eventTriggered && eventVertex) |
96fcc8a7 | 1081 | { |
81be4ee8 | 1082 | fdNdEtaAnalysisMC->FillEvent(vtxMC[2], multAxis); |
1083 | fdNdEtaAnalysisESD->FillEvent(vtxMC[2], multAxis); | |
96fcc8a7 | 1084 | } |
0f67a57c | 1085 | |
81be4ee8 | 1086 | Int_t processType2 = processType; |
1087 | if (oneParticleEvent) | |
1088 | processType2 |= AliPWG0Helper::kOnePart; | |
1089 | ||
1090 | // stuff regarding the vertex reco correction and trigger bias correction | |
1091 | fdNdEtaCorrection->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); | |
96fcc8a7 | 1092 | |
69b09e3b | 1093 | if (fOption.Contains("process-types")) |
96fcc8a7 | 1094 | { |
1095 | // non diffractive | |
a7f69e56 | 1096 | if (processType == AliPWG0Helper::kND) |
81be4ee8 | 1097 | fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 1098 | |
1099 | // single diffractive | |
6b62a9c7 | 1100 | if (processType == AliPWG0Helper::kSD) |
81be4ee8 | 1101 | fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 1102 | |
1103 | // double diffractive | |
6b62a9c7 | 1104 | if (processType == AliPWG0Helper::kDD) |
81be4ee8 | 1105 | fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
0f67a57c | 1106 | } |
69b09e3b | 1107 | |
1108 | if (fOption.Contains("particle-species")) | |
1109 | for (Int_t id=0; id<4; id++) | |
81be4ee8 | 1110 | fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
69b09e3b | 1111 | |
1112 | if (etaArr) | |
1113 | delete[] etaArr; | |
1114 | if (labelArr) | |
1115 | delete[] labelArr; | |
1116 | if (labelArr2) | |
1117 | delete[] labelArr2; | |
1118 | if (thirdDimArr) | |
1119 | delete[] thirdDimArr; | |
1120 | if (deltaPhiArr) | |
1121 | delete[] deltaPhiArr; | |
0f67a57c | 1122 | } |
1123 | ||
1124 | void AlidNdEtaCorrectionTask::Terminate(Option_t *) | |
1125 | { | |
1126 | // The Terminate() function is the last function to be called during | |
1127 | // a query. It always runs on the client, it can be used to present | |
1128 | // the results graphically or save the results to file. | |
1129 | ||
1130 | fOutput = dynamic_cast<TList*> (GetOutputData(0)); | |
1131 | if (!fOutput) { | |
1132 | Printf("ERROR: fOutput not available"); | |
1133 | return; | |
1134 | } | |
1135 | ||
1136 | fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction")); | |
1137 | fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC")); | |
1138 | fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD")); | |
1139 | if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD) | |
1140 | { | |
1141 | AliDebug(AliLog::kError, "Could not read object from output list"); | |
1142 | return; | |
1143 | } | |
1144 | ||
1145 | fdNdEtaCorrection->Finish(); | |
1146 | ||
1147 | TString fileName; | |
1148 | fileName.Form("correction_map%s.root", fOption.Data()); | |
1149 | TFile* fout = new TFile(fileName, "RECREATE"); | |
1150 | ||
1c15d51a | 1151 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts")); |
0f67a57c | 1152 | if (fEsdTrackCuts) |
1153 | fEsdTrackCuts->SaveHistograms("esd_track_cuts"); | |
1c15d51a | 1154 | |
1155 | fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPrim")); | |
1156 | if (fEsdTrackCutsPrim) | |
1157 | fEsdTrackCutsPrim->SaveHistograms("esd_track_cuts_primaries"); | |
1158 | ||
1159 | fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsSec")); | |
1160 | if (fEsdTrackCutsSec) | |
1161 | fEsdTrackCutsSec->SaveHistograms("esd_track_cuts_secondaries"); | |
1162 | ||
0f67a57c | 1163 | fdNdEtaCorrection->SaveHistograms(); |
1164 | fdNdEtaAnalysisMC->SaveHistograms(); | |
1165 | fdNdEtaAnalysisESD->SaveHistograms(); | |
1166 | ||
69b09e3b | 1167 | if (fOutput->FindObject("dndeta_correction_ND")) |
1168 | { | |
1169 | fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_ND")); | |
1170 | fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_SD")); | |
1171 | fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_DD")); | |
1172 | } | |
1173 | else | |
1174 | { | |
1175 | fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_pi")); | |
1176 | fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_K")); | |
1177 | fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_p")); | |
1178 | fdNdEtaCorrectionSpecial[3] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_other")); | |
1179 | } | |
1180 | for (Int_t i=0; i<4; ++i) | |
1181 | if (fdNdEtaCorrectionSpecial[i]) | |
1182 | fdNdEtaCorrectionSpecial[i]->SaveHistograms(); | |
96fcc8a7 | 1183 | |
1c15d51a | 1184 | fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1")); |
1185 | if (fTemp1) | |
1186 | fTemp1->Write(); | |
1187 | fTemp2 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp2")); | |
1188 | if (fTemp2) | |
1189 | fTemp2->Write(); | |
96fcc8a7 | 1190 | |
0fc41645 | 1191 | fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation")); |
54b096ef | 1192 | if (fVertexCorrelation) |
1193 | fVertexCorrelation->Write(); | |
a7f69e56 | 1194 | fVertexCorrelationShift = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexCorrelationShift")); |
51f6de65 | 1195 | if (fVertexCorrelationShift) |
a7f69e56 | 1196 | { |
1197 | ((TH2*) fVertexCorrelationShift->Project3D("yx"))->FitSlicesY(); | |
51f6de65 | 1198 | fVertexCorrelationShift->Write(); |
a7f69e56 | 1199 | } |
0fc41645 | 1200 | fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile")); |
54b096ef | 1201 | if (fVertexProfile) |
1202 | fVertexProfile->Write(); | |
1c15d51a | 1203 | fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift")); |
1204 | if (fVertexShift) | |
1205 | fVertexShift->Write(); | |
a7f69e56 | 1206 | fVertexShiftNorm = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexShiftNorm")); |
54b096ef | 1207 | if (fVertexShiftNorm) |
a7f69e56 | 1208 | { |
1209 | fVertexShiftNorm->ProjectionX(); | |
54b096ef | 1210 | fVertexShiftNorm->Write(); |
a7f69e56 | 1211 | } |
3d7758c1 | 1212 | |
0fc41645 | 1213 | fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation")); |
3d7758c1 | 1214 | if (fEtaCorrelation) |
1215 | fEtaCorrelation->Write(); | |
51f6de65 | 1216 | fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift")); |
1217 | if (fEtaCorrelationShift) | |
a7f69e56 | 1218 | { |
1219 | fEtaCorrelationShift->FitSlicesY(); | |
51f6de65 | 1220 | fEtaCorrelationShift->Write(); |
a7f69e56 | 1221 | } |
0fc41645 | 1222 | fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile")); |
3d7758c1 | 1223 | if (fEtaProfile) |
1224 | fEtaProfile->Write(); | |
1c15d51a | 1225 | fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution")); |
1226 | if (fEtaResolution) | |
1227 | fEtaResolution->Write(); | |
69b09e3b | 1228 | fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution")); |
1c15d51a | 1229 | if (fpTResolution) |
69b09e3b | 1230 | { |
1231 | fpTResolution->FitSlicesY(); | |
1c15d51a | 1232 | fpTResolution->Write(); |
69b09e3b | 1233 | } |
3d7758c1 | 1234 | |
0fc41645 | 1235 | fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll")); |
4ef701a0 | 1236 | if (fMultAll) |
1237 | fMultAll->Write(); | |
1238 | ||
0fc41645 | 1239 | fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr")); |
4ef701a0 | 1240 | if (fMultTr) |
1241 | fMultTr->Write(); | |
0fc41645 | 1242 | |
1243 | fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx")); | |
4ef701a0 | 1244 | if (fMultVtx) |
1245 | fMultVtx->Write(); | |
1246 | ||
3d7758c1 | 1247 | for (Int_t i=0; i<8; ++i) |
0fc41645 | 1248 | { |
69b09e3b | 1249 | fDeltaPhi[i] = dynamic_cast<TH2*> (fOutput->FindObject(Form("fDeltaPhi_%d", i))); |
4ef701a0 | 1250 | if (fDeltaPhi[i]) |
1251 | fDeltaPhi[i]->Write(); | |
0fc41645 | 1252 | } |
54b096ef | 1253 | |
0fc41645 | 1254 | fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats")); |
50ec344d | 1255 | if (fEventStats) |
1256 | fEventStats->Write(); | |
1257 | ||
1c15d51a | 1258 | fPIDParticles = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDParticles")); |
1259 | if (fPIDParticles) | |
1260 | fPIDParticles->Write(); | |
1261 | ||
1262 | fPIDTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDTracks")); | |
1263 | if (fPIDTracks) | |
1264 | fPIDTracks->Write(); | |
0f67a57c | 1265 | |
1c15d51a | 1266 | //fdNdEtaCorrection->DrawHistograms(); |
0f67a57c | 1267 | |
96fcc8a7 | 1268 | Printf("Writting result to %s", fileName.Data()); |
1269 | ||
0f67a57c | 1270 | if (fPIDParticles && fPIDTracks) |
1271 | { | |
0f67a57c | 1272 | TDatabasePDG* pdgDB = new TDatabasePDG; |
1273 | ||
1274 | for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i) | |
1275 | if (fPIDParticles->GetBinContent(i) > 0) | |
1c15d51a | 1276 | { |
1277 | TObject* pdgParticle = pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i)); | |
1278 | printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), (pdgParticle) ? pdgParticle->GetName() : "not found", (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1)); | |
1279 | } | |
0f67a57c | 1280 | |
1281 | delete pdgDB; | |
1282 | pdgDB = 0; | |
1283 | } | |
1c15d51a | 1284 | |
1285 | fout->Write(); | |
1286 | fout->Close(); | |
0f67a57c | 1287 | } |
1288 | ||
1289 | Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle) | |
1290 | { | |
1291 | // returns if a particle with this sign should be counted | |
1292 | // this is determined by the value of fSignMode, which should have the same sign | |
1293 | // as the charge | |
1294 | // if fSignMode is 0 all particles are counted | |
1295 | ||
1296 | if (fSignMode == 0) | |
1297 | return kTRUE; | |
1298 | ||
1299 | if (!particle) | |
1300 | { | |
1301 | printf("WARNING: not counting a particle that does not have a pdg particle\n"); | |
1302 | return kFALSE; | |
1303 | } | |
1304 | ||
1305 | Double_t charge = particle->Charge(); | |
1306 | ||
1307 | if (fSignMode > 0) | |
1308 | if (charge < 0) | |
1309 | return kFALSE; | |
1310 | ||
1311 | if (fSignMode < 0) | |
1312 | if (charge > 0) | |
1313 | return kFALSE; | |
1314 | ||
1315 | return kTRUE; | |
1316 | } | |
96fcc8a7 | 1317 |