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