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