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