Transition PWG0 -> PWGUD
[u/mrichter/AliRoot.git] / PWGUD / selectors / dNdEta / AlidNdEtaCorrectionTask.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaCorrectionTask.h"
4
5 #include <TCanvas.h>
6 #include <TChain.h>
7 #include <TFile.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <TProfile.h>
12 #include <TParticle.h>
13 #include <TParticlePDG.h>
14 #include <TDatabasePDG.h>
15 #include <TRandom.h>
16
17 #include <AliLog.h>
18 #include <AliESDVertex.h>
19 #include <AliESDEvent.h>
20 #include <AliESDRun.h>
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
30 #include "AliESDtrackCuts.h"
31 #include "AliPWG0Helper.h"
32 #include "dNdEtaAnalysis.h"
33 #include "AlidNdEtaCorrection.h"
34 #include "AliTriggerAnalysis.h"
35 #include "AliPhysicsSelection.h"
36
37 ClassImp(AlidNdEtaCorrectionTask)
38
39 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
40   AliAnalysisTask(),
41   fESD(0),
42   fOutput(0),
43   fOption(),
44   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
45   fTrigger(AliTriggerAnalysis::kMB1),
46   fFillPhi(kFALSE),
47   fDeltaPhiCut(-1),
48   fSymmetrize(kFALSE),
49   fMultAxisEta1(kFALSE),
50   fDiffTreatment(AliPWG0Helper::kMCFlags),
51   fSignMode(0),
52   fOnlyPrimaries(kFALSE),
53   fStatError(0),
54   fSystSkipParticles(kFALSE),
55   fEsdTrackCuts(0),
56   fdNdEtaCorrection(0),
57   fdNdEtaAnalysisMC(0),
58   fdNdEtaAnalysisESD(0),
59   fPIDParticles(0),
60   fPIDTracks(0),
61   fVertexCorrelation(0),
62   fVertexCorrelationShift(0),
63   fVertexProfile(0),
64   fVertexShift(0),
65   fVertexShiftNorm(0),
66   fEtaCorrelation(0),
67   fEtaCorrelationShift(0),
68   fEtaProfile(0),
69   fEtaResolution(0),
70   fDeltaPhiCorrelation(0),
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   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)
96 {
97   //
98   // Constructor. Initialization of pointers
99   //
100
101   for (Int_t i=0; i<4; i++)
102     fdNdEtaCorrectionSpecial[i] = 0;
103   
104   for (Int_t i=0; i<8; i++)
105     fDeltaPhi[i] = 0;
106
107   AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
108 }
109
110 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
111   AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
112   fESD(0),
113   fOutput(0),
114   fOption(opt),
115   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
116   fTrigger(AliTriggerAnalysis::kMB1),
117   fFillPhi(kFALSE),
118   fDeltaPhiCut(0),
119   fSymmetrize(kFALSE),
120   fMultAxisEta1(kFALSE),
121   fDiffTreatment(AliPWG0Helper::kMCFlags),
122   fSignMode(0),
123   fOnlyPrimaries(kFALSE),
124   fStatError(0),
125   fSystSkipParticles(kFALSE),
126   fEsdTrackCuts(0),
127   fdNdEtaCorrection(0),
128   fdNdEtaAnalysisMC(0),
129   fdNdEtaAnalysisESD(0),
130   fPIDParticles(0),
131   fPIDTracks(0),
132   fVertexCorrelation(0),
133   fVertexCorrelationShift(0),
134   fVertexProfile(0),
135   fVertexShift(0),
136   fVertexShiftNorm(0),
137   fEtaCorrelation(0),
138   fEtaCorrelationShift(0),
139   fEtaProfile(0),
140   fEtaResolution(0),
141   fDeltaPhiCorrelation(0),
142   fpTResolution(0),
143   fEsdTrackCutsPrim(0),
144   fEsdTrackCutsSec(0),
145   fTemp1(0),
146   fTemp2(0),
147   fMultAll(0),
148   fMultTr(0),
149   fMultVtx(0),
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)
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());
175
176   for (Int_t i=0; i<4; i++)
177     fdNdEtaCorrectionSpecial[i] = 0;
178   
179   for (Int_t i=0; i<8; i++)
180     fDeltaPhi[i] = 0;
181
182   AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning);
183 }
184
185 AlidNdEtaCorrectionTask::~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
193         
194   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
195     delete fOutput;
196     fOutput = 0;
197   }
198         
199 }
200
201 //________________________________________________________________________
202 void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
203 {
204   // Connect ESD
205   // Called once
206
207   Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
208
209   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
210
211   if (!esdH) {
212     Printf("ERROR: Could not get ESDInputHandler");
213   } else {
214     fESD = esdH->GetEvent();
215
216     // Enable only the needed branches
217     esdH->SetActiveBranches("AliESDHeader Vertex");
218
219     if (fAnalysisMode & AliPWG0Helper::kSPD)
220       esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
221
222     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
223       esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
224     }
225   }
226
227   // disable info messages of AliMCEvent (per event)
228   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
229 }
230
231 void AlidNdEtaCorrectionTask::CreateOutputObjects()
232 {
233   // create result objects and add to output list
234
235   AliDebug(2,Form("*********************************** fOption = %s", fOption.Data()));
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   }
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   }
257
258   fOutput = new TList;
259   fOutput->SetOwner();
260
261   fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode);
262   fOutput->Add(fdNdEtaCorrection);
263
264   fPIDParticles = new TH1F("fPIDParticles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
265   fOutput->Add(fPIDParticles);
266
267   fPIDTracks = new TH1F("fPIDTracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
268   fOutput->Add(fPIDTracks);
269
270   fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode);
271   fOutput->Add(fdNdEtaAnalysisMC);
272
273   fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode);
274   fOutput->Add(fdNdEtaAnalysisESD);
275
276   if (fEsdTrackCuts)
277   {
278     fEsdTrackCutsPrim = static_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim"));
279     fEsdTrackCutsSec  = static_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec"));
280     fEsdTrackCutsCheck  = static_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsCheck"));
281     fEsdTrackCutsCheck->SetPtRange(0.15);
282     fEsdTrackCutsCheck->SetEtaRange(-1.,1.);
283     fOutput->Add(fEsdTrackCutsPrim);
284     fOutput->Add(fEsdTrackCutsSec);
285   }
286
287   if (fOption.Contains("process-types")) {
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);
291
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]);
305   }
306
307   
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); 
310   fOutput->Add(fTemp1);
311   
312   fTemp2 = new TH2F("fTemp2", "fTemp2", 300, -15, 15, 80, -2.0, 2.0); 
313   fOutput->Add(fTemp2);
314
315   fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
316   fOutput->Add(fVertexCorrelation);
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);
318   fOutput->Add(fVertexCorrelationShift);
319   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
320   fOutput->Add(fVertexProfile);
321   fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
322   fOutput->Add(fVertexShift);
323   fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx);rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5);
324   fOutput->Add(fVertexShiftNorm);
325
326   fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
327   fOutput->Add(fEtaCorrelation);
328   fEtaCorrelationAllESD = new TH2F("fEtaCorrelationAllESD", "fEtaCorrelationAllESD;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
329   fOutput->Add(fEtaCorrelationAllESD);
330   fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
331   fOutput->Add(fEtaCorrelationShift);
332   fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
333   fOutput->Add(fEtaProfile);
334   fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2);
335   fOutput->Add(fEtaResolution);
336
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);
338   fOutput->Add(fpTResolution);
339
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);
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
349   fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
350   fOutput->Add(fMultAll);
351   fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
352   fOutput->Add(fMultVtx);
353   fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
354   fOutput->Add(fMultTr);
355
356   for (Int_t i=0; i<8; i++)
357   {
358     fDeltaPhi[i] = new TH2F(Form("fDeltaPhi_%d", i), ";#Delta phi;#phi;Entries", 2000, -0.1, 0.1, 18 * 5, 0, TMath::TwoPi());
359     fOutput->Add(fDeltaPhi[i]);
360   }
361
362   fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 109, -6.5, 102.5, 4, -0.5, 3.5);
363   fOutput->Add(fEventStats);
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
369
370   fEventStats->GetXaxis()->SetBinLabel(108, "INEL=0");   // x = -101
371   fEventStats->GetXaxis()->SetBinLabel(109, "INEL>0");   // x = -102
372   
373   for (Int_t i=-1; i<100; i++)
374     fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i)); 
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");
380
381   if (fEsdTrackCuts)
382   {
383     fEsdTrackCuts->SetName("fEsdTrackCuts");
384     // TODO like this we send an empty object back...
385     fOutput->Add(fEsdTrackCuts->Clone());
386   }
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);
404 }
405
406 void AlidNdEtaCorrectionTask::Exec(Option_t*)
407 {
408   // process the event
409
410   fEventNumber++;
411   // Check prerequisites
412   if (!fESD)
413   {
414     AliDebug(AliLog::kError, "ESD branch not available");
415     return;
416   }
417
418   if (fOnlyPrimaries)
419     Printf("WARNING: Processing only primaries. For systematical studies only!");
420     
421   if (fStatError > 0)
422     Printf("WARNING: Statistical error evaluation active!");
423     
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();
432
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     
444   if (!eventTriggered)
445     Printf("No trigger");
446
447   // post the data already here
448   PostData(0, fOutput);
449   
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
477   // get process type
478   Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
479   
480   AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
481
482   if (processType == AliPWG0Helper::kInvalidProcess)
483   {
484     AliDebug(AliLog::kWarning, "Unknown process type. Setting to ND");
485     processType = AliPWG0Helper::kND;
486   }
487
488   // get the MC vertex
489   AliGenEventHeader* genHeader = header->GenEventHeader();
490   TArrayF vtxMC(3);
491   genHeader->PrimaryVertex(vtxMC);
492   fVtxMC->Fill(vtxMC[2]);
493   AliDebug(2,Form("MC vtx: x = %.6f, y = %.6f, z = %.6f",vtxMC[0],vtxMC[1],vtxMC[2]));
494
495   // get the ESD vertex
496   Bool_t eventVertex = kFALSE;
497   Double_t vtx[3];
498   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
499   if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
500   {
501     vtxESD->GetXYZ(vtx);
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     }
510   }
511   else
512     vtxESD = 0;
513     
514   // create list of (label, eta, pt) tuples
515   Int_t inputCount = 0;
516   Int_t* labelArr = 0;
517   Int_t* labelArr2 = 0; // only for case of SPD
518   Float_t* etaArr = 0;
519   Float_t* thirdDimArr = 0;
520   Float_t* deltaPhiArr = 0;
521   if (fAnalysisMode & AliPWG0Helper::kSPD)
522   {
523     if (vtxESD)
524     {
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))
559           continue;
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       }
584       
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;
598     }
599   }
600   else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
601   {
602     if (!fEsdTrackCuts)
603     {
604       AliDebug(AliLog::kError, "fESDTrackCuts not available");
605       return;
606     }
607     
608     Bool_t foundInEta10 = kFALSE;
609     
610     if (vtxESD)
611     {
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
650       // get multiplicity from ESD tracks
651       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
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++)
664       {
665         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
666         if (!esdTrack)
667         {
668           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
669           continue;
670         }
671         
672         AliDebug(3,Form("particle %d: pt = %.6f, eta = %.6f",i,esdTrack->Pt(), esdTrack->Eta())); 
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
678         if (esdTrack->Pt() < fPtMin)
679           continue;
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;
689         }
690         
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)
711           foundInEta10 = kTRUE;
712
713         etaArr[inputCount] = esdTrack->Eta();
714         if (fSymmetrize)
715           etaArr[inputCount] = TMath::Abs(etaArr[inputCount]);
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       }
722
723       delete list;
724
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)
727       {
728         // collect values for primaries and secondaries
729         for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
730         {
731           AliESDtrack* track = 0;
732   
733           if (fAnalysisMode & AliPWG0Helper::kTPC)
734             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack);
735           else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
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   
753           if (TMath::Abs(track->Eta()) < 0.8 && track->Pt() > 0.15)
754           {
755             if (stack->IsPhysicalPrimary(label))
756             {
757               // primary
758               if (fEsdTrackCutsPrim->AcceptTrack(track)) 
759               {
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 //                 }
771               }
772             }
773             else
774             {
775               // secondary
776               fEsdTrackCutsSec->AcceptTrack(track);
777             }
778           }
779   
780           // TODO mem leak in the continue statements above
781           if (fAnalysisMode & AliPWG0Helper::kTPC)
782             delete track;
783         }
784       }
785     }
786     
787     if (!foundInEta10)
788       eventTriggered = kFALSE;
789     else{
790      //Printf("The event triggered. Its number in file is %d",fESD->GetEventNumberInFile());
791       fNumberEvent->Fill(fESD->GetEventNumberInFile());
792     }
793   }
794   else
795     return;
796
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   
812   // loop over mc particles
813   Int_t nPrim  = stack->GetNprimary();
814   Int_t nAccepted = 0;
815
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     
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)
833     {
834       oneParticleEvent = kTRUE;
835       fNumberEventMC->Fill(fESD->GetEventNumberInFile());
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   
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)
857     {
858       //if (TMath::Abs(particle->GetPdgCode()) > 3000 && TMath::Abs(particle->Eta()) < 1.0)
859       //  fPIDParticles->Fill(particle->GetPdgCode());
860       continue;
861     }
862
863     if (SignOK(particle->GetPDG()) == kFALSE)
864       continue;
865     
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)
871       fPIDParticles->Fill(particle->GetPdgCode());
872
873     Float_t eta = particle->Eta();
874     if (fSymmetrize)
875       eta = TMath::Abs(eta);
876     
877     Float_t thirdDim = -1;
878     if (fAnalysisMode & AliPWG0Helper::kSPD)
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);
894
895     Int_t processType2 = processType;
896     if (oneParticleEvent)
897       processType2 |= AliPWG0Helper::kOnePart;
898
899     fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
900
901     if (fOption.Contains("process-types"))
902     {
903       // non diffractive
904       if (processType==AliPWG0Helper::kND)
905         fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
906
907       // single diffractive
908       if (processType==AliPWG0Helper::kSD)
909         fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
910
911       // double diffractive
912       if (processType==AliPWG0Helper::kDD)
913         fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
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       }
926       fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2);
927     }
928
929     if (eventTriggered)
930       if (eventVertex)
931         fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, thirdDim);
932
933     // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
934     if (TMath::Abs(eta) < 1.5) // && pt > 0.2)
935       nAccepted++;
936   }
937
938   if (AliPWG0Helper::GetLastProcessType() >= -1)
939     fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
940
941   fMultAll->Fill(nAccepted);
942   if (eventTriggered) {
943     fMultTr->Fill(nAccepted);
944     if (eventVertex)
945       fMultVtx->Fill(nAccepted);
946   }
947
948   Int_t processed = 0;
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   }
957
958   Int_t nEta05 = 0;
959   Int_t nEta10 = 0;
960   for (Int_t i=0; i<inputCount; ++i)
961   {
962     if (TMath::Abs(etaArr[i]) < 0.5)
963       nEta05++;
964     if (TMath::Abs(etaArr[i]) < 1.0)
965       nEta10++;
966   }
967   
968   for (Int_t i=0; i<inputCount; ++i)
969   {
970     Int_t label = labelArr[i];
971     Int_t label2 = labelArr2[i];
972
973     if (label < 0)
974     {
975       Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label);
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     }
985
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)
991       fPIDTracks->Fill(particle->GetPdgCode());
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     }
1015
1016     if (eventTriggered && eventVertex)
1017     {
1018       if (SignOK(particle->GetPDG()) == kFALSE)
1019         continue;
1020
1021       processed++;
1022
1023       // resolutions
1024       if (fSymmetrize)
1025         fEtaResolution->Fill(TMath::Abs(particle->Eta()) - etaArr[i]);
1026       else
1027         fEtaResolution->Fill(particle->Eta() - etaArr[i]);
1028
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){
1035           fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
1036           fpTCorrelation->Fill(particle->Pt(),thirdDimArr[i]);
1037           fpTCorrelationShift->Fill(particle->Pt(),particle->Pt()-thirdDimArr[i]);
1038         }
1039       }
1040
1041       Float_t eta = -999;
1042       Float_t thirdDim = -1;
1043
1044       Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
1045       // in case of same label the MC values are filled, otherwise (background) the reconstructed values
1046       if (label == label2)
1047       {
1048         eta = particle->Eta();
1049         if (fSymmetrize)
1050           eta = TMath::Abs(eta);
1051         
1052         if (fAnalysisMode & AliPWG0Helper::kSPD)
1053         {
1054           if (fFillPhi)
1055           {
1056             thirdDim = particle->Phi();
1057           }
1058           else
1059             thirdDim = inputCount;
1060         }
1061         else
1062           thirdDim = particle->Pt();
1063       }
1064       else
1065       {
1066         if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi)
1067         {
1068           thirdDim = (fMultAxisEta1) ? nEta10 : inputCount;
1069         }
1070         else
1071           thirdDim = thirdDimArr[i];
1072
1073         eta = etaArr[i];
1074       }
1075
1076       Bool_t fillTrack = kTRUE;
1077
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)
1100       {
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);
1108         fTemp2->Fill(vtxMC[2], eta);
1109       }
1110       
1111       // eta comparison for tracklets with the same label (others are background)
1112       if (label == label2)
1113       {
1114         Float_t eta2 = particle->Eta();
1115         if (fSymmetrize)
1116           eta2 = TMath::Abs(eta2);
1117         
1118         fEtaProfile->Fill(eta2, eta2 - etaArr[i]);
1119         fEtaCorrelation->Fill(eta2, etaArr[i]);
1120         fEtaCorrelationShift->Fill(eta2, eta2 - etaArr[i]);
1121       }
1122
1123       if (fSymmetrize)
1124         fdNdEtaAnalysisESD->FillTrack(vtxMC[2], TMath::Abs(particle->Eta()), thirdDim);
1125       else
1126         fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
1127
1128       if (fOption.Contains("process-types"))
1129       {
1130         // non diffractive
1131         if (processType == AliPWG0Helper::kND)
1132           fdNdEtaCorrectionSpecial[0]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
1133
1134         // single diffractive
1135         if (processType == AliPWG0Helper::kSD)
1136           fdNdEtaCorrectionSpecial[1]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
1137
1138         // double diffractive
1139         if (processType == AliPWG0Helper::kDD)
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);
1157       }
1158
1159       // control histograms
1160       Int_t hist = -1;
1161       if (label == label2)
1162       {
1163         if (firstIsPrim)
1164         {
1165           hist = 0;
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;
1180
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
1208       }
1209       else
1210         hist = 7;
1211       
1212       fDeltaPhi[hist]->Fill(deltaPhiArr[i], thirdDimArr[i]);
1213     }
1214   }
1215   
1216   if (primCount)
1217   {
1218     delete[] primCount;
1219     primCount = 0;
1220   }
1221
1222   if (processed < inputCount)
1223     Printf("Only %d out of %d track(let)s were used", processed, inputCount); 
1224
1225   if (eventTriggered && eventVertex)
1226   {
1227     Double_t diff = vtxMC[2] - vtx[2];
1228     fVertexShift->Fill(diff);
1229     
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]);
1233   
1234     if (vtxESD->IsFromVertexerZ() && inputCount > 0)
1235       fVertexShiftNorm->Fill(diff, vtxESD->GetNContributors());
1236   }
1237
1238   Int_t multAxis = inputCount;
1239   if (fMultAxisEta1)
1240     multAxis = nEta10;
1241
1242   if (eventTriggered && eventVertex)
1243   {
1244     fdNdEtaAnalysisMC->FillEvent(vtxMC[2], multAxis);
1245     fdNdEtaAnalysisESD->FillEvent(vtxMC[2], multAxis);
1246   }
1247
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);
1254
1255   if (fOption.Contains("process-types"))
1256   {
1257     // non diffractive
1258     if (processType == AliPWG0Helper::kND)
1259       fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1260
1261     // single diffractive
1262     if (processType == AliPWG0Helper::kSD)
1263       fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1264
1265     // double diffractive
1266     if (processType == AliPWG0Helper::kDD)
1267       fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
1268   }
1269   
1270   if (fOption.Contains("particle-species"))
1271     for (Int_t id=0; id<4; id++)
1272       fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2);
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;
1284 }
1285
1286 void 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
1313   fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
1314   if (fEsdTrackCuts)
1315     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
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
1325   fdNdEtaCorrection->SaveHistograms();
1326   fdNdEtaAnalysisMC->SaveHistograms();
1327   fdNdEtaAnalysisESD->SaveHistograms();
1328
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();
1345
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();
1352
1353   fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
1354   if (fVertexCorrelation)
1355     fVertexCorrelation->Write();
1356   fVertexCorrelationShift = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexCorrelationShift"));
1357   if (fVertexCorrelationShift)
1358   {
1359     ((TH2*) fVertexCorrelationShift->Project3D("yx"))->FitSlicesY();
1360     fVertexCorrelationShift->Write();
1361   }
1362   fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
1363   if (fVertexProfile)
1364     fVertexProfile->Write();
1365   fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift"));
1366   if (fVertexShift)
1367     fVertexShift->Write();
1368   fVertexShiftNorm = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexShiftNorm"));
1369   if (fVertexShiftNorm)
1370   {
1371     fVertexShiftNorm->ProjectionX();
1372     fVertexShiftNorm->Write();
1373   }  
1374
1375   fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
1376   if (fEtaCorrelation)
1377     fEtaCorrelation->Write();
1378   fEtaCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationAllESD"));
1379   if (fEtaCorrelationAllESD)
1380     fEtaCorrelationAllESD->Write();
1381   fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
1382   if (fEtaCorrelationShift)
1383   {
1384     fEtaCorrelationShift->FitSlicesY();
1385     fEtaCorrelationShift->Write();
1386   }
1387   fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
1388   if (fEtaProfile)
1389     fEtaProfile->Write();
1390   fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution"));
1391   if (fEtaResolution)
1392     fEtaResolution->Write();
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   }
1402   fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution"));
1403   if (fpTResolution)
1404   {
1405     fpTResolution->FitSlicesY();
1406     fpTResolution->Write();
1407   }
1408
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   }
1418   fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
1419   if (fMultAll)
1420     fMultAll->Write();
1421
1422   fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr"));
1423   if (fMultTr)
1424     fMultTr->Write();
1425
1426   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
1427   if (fMultVtx)
1428     fMultVtx->Write();
1429
1430   for (Int_t i=0; i<8; ++i)
1431   {
1432     fDeltaPhi[i] = dynamic_cast<TH2*> (fOutput->FindObject(Form("fDeltaPhi_%d", i)));
1433     if (fDeltaPhi[i])
1434       fDeltaPhi[i]->Write();
1435   }
1436
1437   fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats"));
1438   if (fEventStats)
1439     fEventStats->Write();
1440
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();
1448
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
1472  //fdNdEtaCorrection->DrawHistograms();
1473
1474   Printf("Writing result to %s", fileName.Data());
1475
1476   if (fPIDParticles && fPIDTracks)
1477   {
1478     TDatabasePDG* pdgDB = new TDatabasePDG;
1479
1480     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
1481       if (fPIDParticles->GetBinContent(i) > 0)
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       }
1486
1487     delete pdgDB;
1488     pdgDB = 0;
1489   }
1490
1491   fout->Write();
1492   fout->Close();
1493 }
1494
1495 Bool_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 }
1523
1524 Double_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
1538 Double_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