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