4a3b4678fc4cca0b17dd04ef602d3ecf99f8c688
[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 <TProfile.h>
12 #include <TParticle.h>
13
14 #include <AliLog.h>
15 #include <AliESDVertex.h>
16 #include <AliESDEvent.h>
17 #include <AliStack.h>
18 #include <AliHeader.h>
19 #include <AliGenEventHeader.h>
20 #include <AliMultiplicity.h>
21 #include <AliAnalysisManager.h>
22 #include <AliMCEventHandler.h>
23 #include <AliMCEvent.h>
24 #include <AliESDInputHandler.h>
25
26 #include "esdTrackCuts/AliESDtrackCuts.h"
27 #include "AliPWG0Helper.h"
28 //#include "AliCorrection.h"
29 //#include "AliCorrectionMatrix3D.h"
30 #include "dNdEta/dNdEtaAnalysis.h"
31 #include "dNdEta/AlidNdEtaCorrection.h"
32 #include "AliVertexerTracks.h"
33
34 ClassImp(AlidNdEtaCorrectionTask)
35
36 AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
37   AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
38   fESD(0),
39   fOutput(0),
40   fOption(opt),
41   fAnalysisMode(AliPWG0Helper::kTPC),
42   fSignMode(0),
43   fOnlyPrimaries(kFALSE),
44   fEsdTrackCuts(0),
45   fdNdEtaCorrection(0),
46   fdNdEtaAnalysisMC(0),
47   fdNdEtaAnalysisESD(0),
48   fPIDParticles(0),
49   fPIDTracks(0),
50   fVertexCorrelation(0),
51   fVertexProfile(0),
52   fVertexShiftNorm(0),
53   fEtaCorrelation(0),
54   fEtaProfile(0),
55   fSigmaVertexTracks(0),
56   fSigmaVertexPrim(0),
57   fMultAll(0),
58   fMultTr(0),
59   fMultVtx(0),
60   fEventStats(0)
61 {
62   //
63   // Constructor. Initialization of pointers
64   //
65
66   // Define input and output slots here
67   DefineInput(0, TChain::Class());
68   DefineOutput(0, TList::Class());
69
70   for (Int_t i=0; i<2; i++)
71     fdNdEtaCorrectionProcessType[i] = 0;
72   
73   for (Int_t i=0; i<8; i++)
74     fDeltaPhi[i] = 0;
75 }
76
77 AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
78 {
79   //
80   // Destructor
81   //
82
83   // histograms are in the output list and deleted when the output
84   // list is deleted by the TSelector dtor
85
86   if (fOutput) {
87     delete fOutput;
88     fOutput = 0;
89   }
90 }
91
92 //________________________________________________________________________
93 void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
94 {
95   // Connect ESD
96   // Called once
97
98   Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
99
100   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
101   if (!tree) {
102     Printf("ERROR: Could not read tree from input slot 0");
103   } else {
104     // Disable all branches and enable only the needed ones
105     //tree->SetBranchStatus("*", 0);
106
107     tree->SetBranchStatus("fTriggerMask", 1);
108     tree->SetBranchStatus("fSPDVertex*", 1);
109     // PrimaryVertex
110
111     if (fAnalysisMode == AliPWG0Helper::kSPD)
112       tree->SetBranchStatus("fSPDMult*", 1);
113
114     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
115       AliESDtrackCuts::EnableNeededBranches(tree);
116       tree->SetBranchStatus("fTracks.fLabel", 1);
117     }
118
119     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
120
121     if (!esdH) {
122       Printf("ERROR: Could not get ESDInputHandler");
123     } else
124       fESD = esdH->GetEvent();
125   }
126
127   // disable info messages of AliMCEvent (per event)
128   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
129 }
130
131 void AlidNdEtaCorrectionTask::CreateOutputObjects()
132 {
133   // create result objects and add to output list
134
135   if (fOption.Contains("only-positive"))
136   {
137     Printf("INFO: Processing only positive particles.");
138     fSignMode = 1;
139   }
140   else if (fOption.Contains("only-negative"))
141   {
142     Printf("INFO: Processing only negative particles.");
143     fSignMode = -1;
144   }
145
146   fOutput = new TList;
147   fOutput->SetOwner();
148
149   fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode);
150   fOutput->Add(fdNdEtaCorrection);
151
152   fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
153   fOutput->Add(fPIDParticles);
154
155   fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
156   fOutput->Add(fPIDTracks);
157
158   fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode);
159   fOutput->Add(fdNdEtaAnalysisMC);
160
161   fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode);
162   fOutput->Add(fdNdEtaAnalysisESD);
163
164   if (fOption.Contains("process-types")) {
165     fdNdEtaCorrectionProcessType[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
166     fdNdEtaCorrectionProcessType[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
167     fdNdEtaCorrectionProcessType[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode);
168
169     fOutput->Add(fdNdEtaCorrectionProcessType[0]);
170     fOutput->Add(fdNdEtaCorrectionProcessType[1]);
171     fOutput->Add(fdNdEtaCorrectionProcessType[2]);
172   }
173
174   if (fOption.Contains("sigma-vertex"))
175   {
176     fSigmaVertexTracks = new TH1F("fSigmaVertexTracks", "fSigmaVertexTracks;Nsigma2vertex;NacceptedTracks", 60, 0.05, 6.05);
177     fSigmaVertexPrim = new TH1F("fSigmaVertexPrim", "fSigmaVertexPrim;Nsigma2vertex;NacceptedPrimaries", 60, 0.05, 6.05);
178     fOutput->Add(fSigmaVertexTracks);
179     fOutput->Add(fSigmaVertexPrim);
180     Printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms");
181   }
182
183   fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 80, -20, 20, 80, -20, 20);
184   fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
185   fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
186   
187   fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
188   fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
189
190   fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
191   fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
192   fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
193
194   for (Int_t i=0; i<8; i++)
195     fDeltaPhi[i] = new TH1F(Form("fDeltaPhi_%d", i), ";#Delta phi;Entries", 2000, -0.1, 0.1);
196
197   fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 2, -0.5, 1.5, 4, -0.5, 3.5);
198   fEventStats->GetXaxis()->SetBinLabel(1, "INEL");
199   fEventStats->GetXaxis()->SetBinLabel(2, "NSD");
200
201   fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
202   fEventStats->GetYaxis()->SetBinLabel(2, "trg");
203   fEventStats->GetYaxis()->SetBinLabel(3, "vtx");
204   fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx");
205 }
206
207 void AlidNdEtaCorrectionTask::Exec(Option_t*)
208 {
209   // process the event
210
211   // Check prerequisites
212   if (!fESD)
213   {
214     AliDebug(AliLog::kError, "ESD branch not available");
215     return;
216   }
217
218   if (fOnlyPrimaries)
219     Printf("WARNING: Processing only primaries. For systematical studies only!");
220
221   // trigger definition
222   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
223
224   if (!eventTriggered)
225     Printf("No trigger");
226
227   // post the data already here
228   PostData(0, fOutput);
229
230   // MC info
231   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
232   if (!eventHandler) {
233     Printf("ERROR: Could not retrieve MC event handler");
234     return;
235   }
236
237   AliMCEvent* mcEvent = eventHandler->MCEvent();
238   if (!mcEvent) {
239     Printf("ERROR: Could not retrieve MC event");
240     return;
241   }
242
243   AliStack* stack = mcEvent->Stack();
244   if (!stack)
245   {
246     AliDebug(AliLog::kError, "Stack not available");
247     return;
248   }
249
250   AliHeader* header = mcEvent->Header();
251   if (!header)
252   {
253     AliDebug(AliLog::kError, "Header not available");
254     return;
255   }
256
257   // get process type; NB: this only works for Pythia
258   Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
259   AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));
260
261   if (processType<0)
262     AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
263
264   // get the MC vertex
265   AliGenEventHeader* genHeader = header->GenEventHeader();
266   TArrayF vtxMC(3);
267   genHeader->PrimaryVertex(vtxMC);
268
269   // get the ESD vertex
270   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
271
272   Bool_t eventVertex = kFALSE;
273   Double_t vtx[3];
274   if (vtxESD) 
275   {
276     vtxESD->GetXYZ(vtx);
277     eventVertex = kTRUE;
278     
279     Double_t diff = vtxMC[2] - vtx[2];
280     if (vtxESD->GetZRes() > 0) 
281       fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
282   } 
283   else
284     Printf("No vertex found");
285
286   // fill process type
287   Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
288   // INEL
289   fEventStats->Fill(0.0, biny);
290   // NDS
291   if (processType != 92 && processType != 93)
292     fEventStats->Fill(1.0, biny);
293   
294   // create list of (label, eta, pt) tuples
295   Int_t inputCount = 0;
296   Int_t* labelArr = 0;
297   Int_t* labelArr2 = 0; // only for case of SPD
298   Float_t* etaArr = 0;
299   Float_t* ptArr = 0;
300   Float_t* deltaPhiArr = 0;
301   if (fAnalysisMode == AliPWG0Helper::kSPD)
302   {
303     // get tracklets
304     const AliMultiplicity* mult = fESD->GetMultiplicity();
305     if (!mult)
306     {
307       AliDebug(AliLog::kError, "AliMultiplicity not available");
308       return;
309     }
310
311     labelArr = new Int_t[mult->GetNumberOfTracklets()];
312     labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
313     etaArr = new Float_t[mult->GetNumberOfTracklets()];
314     ptArr = new Float_t[mult->GetNumberOfTracklets()];
315     deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
316
317     // get multiplicity from ITS tracklets
318     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
319     {
320       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
321
322       Float_t deltaPhi = mult->GetDeltaPhi(i);
323       // prevent values to be shifted by 2 Pi()
324       if (deltaPhi < -TMath::Pi())
325         deltaPhi += TMath::Pi() * 2;
326       if (deltaPhi > TMath::Pi())
327         deltaPhi -= TMath::Pi() * 2;
328
329       if (TMath::Abs(deltaPhi) > 1)
330         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
331
332       if (fOnlyPrimaries)
333         if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
334           continue;
335
336       etaArr[inputCount] = mult->GetEta(i);
337       labelArr[inputCount] = mult->GetLabel(i, 0);
338       labelArr2[inputCount] = mult->GetLabel(i, 1);
339       ptArr[inputCount] = 0; // no pt for tracklets
340       deltaPhiArr[inputCount] = deltaPhi;
341       ++inputCount;
342     }
343   }
344   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
345   {
346     if (!fEsdTrackCuts)
347     {
348       AliDebug(AliLog::kError, "fESDTrackCuts not available");
349       return;
350     }
351
352     // get multiplicity from ESD tracks
353     TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
354     Int_t nGoodTracks = list->GetEntries();
355
356     Printf("Accepted %d tracks", nGoodTracks);
357
358     labelArr = new Int_t[nGoodTracks];
359     labelArr2 = new Int_t[nGoodTracks];
360     etaArr = new Float_t[nGoodTracks];
361     ptArr = new Float_t[nGoodTracks];
362     deltaPhiArr = new Float_t[nGoodTracks];
363
364     // loop over esd tracks
365     for (Int_t i=0; i<nGoodTracks; i++)
366     {
367       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
368       if (!esdTrack)
369       {
370         AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
371         continue;
372       }
373
374       etaArr[inputCount] = esdTrack->Eta();
375       labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
376       labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
377       ptArr[inputCount] = esdTrack->Pt();
378       deltaPhiArr[inputCount] = 0; // no delta phi for tracks
379       ++inputCount;
380     }
381
382     delete list;
383   }
384   else
385     return;
386
387   if (eventTriggered && eventVertex)
388   {
389     fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
390     fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
391   }
392
393
394   // loop over mc particles
395   Int_t nPrim  = stack->GetNprimary();
396   Int_t nAccepted = 0;
397
398   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
399   {
400     //Printf("Getting particle %d", iMc);
401     TParticle* particle = stack->Particle(iMc);
402
403     if (!particle)
404       continue;
405
406     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
407       continue;
408
409     if (SignOK(particle->GetPDG()) == kFALSE)
410       continue;
411
412     Float_t eta = particle->Eta();
413     Float_t pt = particle->Pt();
414
415     fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
416
417     if (fdNdEtaCorrectionProcessType[0])
418     {
419       // non diffractive
420       if (processType!=92 && processType!=93 && processType!=94)
421         fdNdEtaCorrectionProcessType[0]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
422
423       // single diffractive
424       if (processType==92 || processType==93)
425         fdNdEtaCorrectionProcessType[1]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
426
427       // double diffractive
428       if (processType==94)
429         fdNdEtaCorrectionProcessType[2]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
430     }
431
432     if (eventTriggered)
433       if (eventVertex)
434         fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
435
436     if (TMath::Abs(eta) < 1 && pt > 0.2)
437       nAccepted++;
438   }
439
440   fMultAll->Fill(nAccepted);
441   if (eventTriggered) {
442     fMultTr->Fill(nAccepted);
443     if (eventVertex)
444       fMultVtx->Fill(nAccepted);
445   }
446
447   Int_t processed = 0;
448
449   for (Int_t i=0; i<inputCount; ++i)
450   {
451     Int_t label = labelArr[i];
452     Int_t label2 = labelArr2[i];
453
454     if (label < 0)
455     {
456       Printf("WARNING: cannot find corresponding mc part for track(let) %d with label %d.", i, label);
457       continue;
458     }
459
460     TParticle* particle = stack->Particle(label);
461     if (!particle)
462     {
463       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
464       continue;
465     }
466     
467     // find particle that is filled in the correction map
468     // this should be the particle that has been reconstructed
469     // for TPC this is clear and always identified by the track's label
470     // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken: 
471     //  L1 L2 (P = primary, S = secondary)
472     //   P P' : --> P
473     //   P S  : --> P
474     //   S P  : --> P
475     //   S S' : --> S
476
477     if (label != label2 && label2 >= 0)
478     {
479       if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2))
480       {
481         particle = stack->Particle(label2);
482         if (!particle)
483         {
484           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2));
485           continue;
486         }
487       }
488     }
489
490     if (eventTriggered && eventVertex)
491     {
492       if (SignOK(particle->GetPDG()) == kFALSE)
493         continue;
494
495       processed++;
496
497       Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
498       // in case of primary the MC values are filled, otherwise (background) the reconstructed values
499       if (label == label2 && firstIsPrim)
500       {
501         fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
502       }
503       else
504       {
505         fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], etaArr[i], ptArr[i]);
506         fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
507       }
508               
509       fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
510
511       fEtaCorrelation->Fill(etaArr[i], particle->Eta());
512
513       if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
514       {
515         fPIDTracks->Fill(particle->GetPdgCode());
516       }
517
518       if (fdNdEtaCorrectionProcessType[0])
519       {
520         // non diffractive
521         if (processType!=92 && processType!=93 && processType!=94)
522           fdNdEtaCorrectionProcessType[0]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
523
524         // single diffractive
525         if (processType==92 || processType==93)
526           fdNdEtaCorrectionProcessType[1]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
527
528         // double diffractive
529         if (processType==94)
530           fdNdEtaCorrectionProcessType[2]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
531       }
532
533       Int_t hist = -1;
534       if (label == label2)
535       {
536         if (firstIsPrim)
537         {
538           hist = 0;
539         }
540         else
541           hist = 1; 
542       }
543       else if (label2 >= 0)
544       {
545         Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2);
546         if (firstIsPrim && secondIsPrim)
547         {
548           hist = 2;
549         }
550         else if (firstIsPrim && !secondIsPrim)
551         {
552           hist = 3;
553           
554           // check if secondary is caused by the primary or it is a fake combination
555           //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2);
556           Int_t mother = label2;
557           while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0)
558           {
559             //Printf("  %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1));
560             if (stack->Particle(mother)->GetMother(0) == label)
561             {
562               /*Printf("The untraceable decay was:");
563               Printf("   from:");
564               particle->Print();
565               Printf("   to (status code %d):", stack->Particle(mother)->GetStatusCode());
566               stack->Particle(mother)->Print();*/
567               hist = 4;
568             }
569             mother = stack->Particle(mother)->GetMother(0);
570           }
571         }
572         else if (!firstIsPrim && secondIsPrim)
573         {
574           hist = 5;
575         }
576         else if (!firstIsPrim && !secondIsPrim)
577         {
578           hist = 6;
579         }
580
581       }
582       else
583         hist = 7;
584       
585       fDeltaPhi[hist]->Fill(deltaPhiArr[i]);
586     }
587   }
588
589   if (processed < inputCount)
590     Printf("Only %d out of %d track(let)s were used", processed, inputCount); 
591
592   if (eventTriggered && eventVertex)
593   {
594     fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
595     fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
596   }
597
598    // stuff regarding the vertex reco correction and trigger bias correction
599   fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
600
601   if (fdNdEtaCorrectionProcessType[0])
602   {
603     // non diffractive
604     if (processType!=92 && processType!=93 && processType!=94)
605       fdNdEtaCorrectionProcessType[0]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
606
607     // single diffractive
608     if (processType==92 || processType==93)
609       fdNdEtaCorrectionProcessType[1]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
610
611     // double diffractive
612     if (processType==94)
613       fdNdEtaCorrectionProcessType[2]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
614   }
615
616   delete[] etaArr;
617   delete[] labelArr;
618   delete[] labelArr2;
619   delete[] ptArr;
620   delete[] deltaPhiArr;
621
622   // fills the fSigmaVertex histogram (systematic study)
623   if (fSigmaVertexTracks)
624   {
625     // save the old value
626     Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
627
628     // set to maximum
629     fEsdTrackCuts->SetMinNsigmaToVertex(6);
630
631     TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
632     Int_t nGoodTracks = list->GetEntries();
633
634     // loop over esd tracks
635     for (Int_t i=0; i<nGoodTracks; i++)
636     {
637       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
638       if (!esdTrack)
639       {
640         AliError(Form("ERROR: Could not retrieve track %d.", i));
641         continue;
642       }
643
644       Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
645
646       for (Double_t nSigma = 0.1; nSigma < 6.05; nSigma += 0.1)
647       {
648         if (sigma2Vertex < nSigma)
649         {
650           fSigmaVertexTracks->Fill(nSigma);
651
652           Int_t label = TMath::Abs(esdTrack->GetLabel());
653           TParticle* particle = stack->Particle(label);
654           if (!particle || label >= nPrim)
655             continue;
656
657           if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
658             fSigmaVertexPrim->Fill(nSigma);
659         }
660       }
661     }
662
663     delete list;
664     list = 0;
665
666     // set back the old value
667     fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
668   }
669 }
670
671 void AlidNdEtaCorrectionTask::Terminate(Option_t *)
672 {
673   // The Terminate() function is the last function to be called during
674   // a query. It always runs on the client, it can be used to present
675   // the results graphically or save the results to file.
676
677   fOutput = dynamic_cast<TList*> (GetOutputData(0));
678   if (!fOutput) {
679     Printf("ERROR: fOutput not available");
680     return;
681   }
682
683   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
684   fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
685   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
686   if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
687   {
688     AliDebug(AliLog::kError, "Could not read object from output list");
689     return;
690   }
691
692   fdNdEtaCorrection->Finish();
693
694   TString fileName;
695   fileName.Form("correction_map%s.root", fOption.Data());
696   TFile* fout = new TFile(fileName, "RECREATE");
697
698   if (fEsdTrackCuts)
699     fEsdTrackCuts->SaveHistograms("esd_track_cuts");
700   fdNdEtaCorrection->SaveHistograms();
701   fdNdEtaAnalysisMC->SaveHistograms();
702   fdNdEtaAnalysisESD->SaveHistograms();
703
704   fdNdEtaCorrectionProcessType[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_ND"));
705   fdNdEtaCorrectionProcessType[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_SD"));
706   fdNdEtaCorrectionProcessType[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_DD"));
707   for (Int_t i=0; i<3; ++i)
708     if (fdNdEtaCorrectionProcessType[i])
709       fdNdEtaCorrectionProcessType[i]->SaveHistograms();
710
711   fSigmaVertexTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertexTracks"));
712   if (fSigmaVertexTracks)
713     fSigmaVertexTracks->Write();
714   fSigmaVertexPrim = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertexPrim"));
715   if (fSigmaVertexPrim)
716     fSigmaVertexPrim->Write();
717
718   if (fVertexCorrelation)
719     fVertexCorrelation->Write();
720   if (fVertexProfile)
721     fVertexProfile->Write();
722   if (fVertexShiftNorm)
723     fVertexShiftNorm->Write();
724
725   if (fEtaCorrelation)
726     fEtaCorrelation->Write();
727   if (fEtaProfile)
728     fEtaProfile->Write();
729
730   if (fMultAll)
731     fMultAll->Write();
732
733   if (fMultTr)
734     fMultTr->Write();
735   
736   if (fMultVtx)
737     fMultVtx->Write();
738
739   for (Int_t i=0; i<8; ++i)
740     if (fDeltaPhi[i])
741       fDeltaPhi[i]->Write();
742
743   if (fEventStats)
744     fEventStats->Write();
745
746   fout->Write();
747   fout->Close();
748
749   //fdNdEtaCorrection->DrawHistograms();
750
751   Printf("Writting result to %s", fileName.Data());
752
753   if (fPIDParticles && fPIDTracks)
754   {
755     new TCanvas("pidcanvas", "pidcanvas", 500, 500);
756
757     fPIDParticles->Draw();
758     fPIDTracks->SetLineColor(2);
759     fPIDTracks->Draw("SAME");
760
761     TDatabasePDG* pdgDB = new TDatabasePDG;
762
763     for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
764       if (fPIDParticles->GetBinContent(i) > 0)
765         printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
766
767     delete pdgDB;
768     pdgDB = 0;
769   }
770 }
771
772 Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
773 {
774   // returns if a particle with this sign should be counted
775   // this is determined by the value of fSignMode, which should have the same sign
776   // as the charge
777   // if fSignMode is 0 all particles are counted
778
779   if (fSignMode == 0)
780     return kTRUE;
781
782   if (!particle)
783   {
784     printf("WARNING: not counting a particle that does not have a pdg particle\n");
785     return kFALSE;
786   }
787
788   Double_t charge = particle->Charge();
789
790   if (fSignMode > 0)
791     if (charge < 0)
792       return kFALSE;
793
794   if (fSignMode < 0)
795     if (charge > 0)
796       return kFALSE;
797
798   return kTRUE;
799 }
800