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