PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx

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