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