updated dndeta analysis
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaCorrectionTask.cxx
CommitLineData
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
37ClassImp(AlidNdEtaCorrectionTask)
38
1c15d51a 39AlidNdEtaCorrectionTask::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 94AlidNdEtaCorrectionTask::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
153AlidNdEtaCorrectionTask::~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//________________________________________________________________________
169void 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
198void 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
341void 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
1124void 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
1289Bool_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