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