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