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