PWG0/multiplicity/AliMultiplicityTask.cxx

   1 /* $Id$ */
   2
   3 #include "AliMultiplicityTask.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TVector3.h>
   9 #include <TChain.h>
  10 #include <TFile.h>
  11 #include <TH1D.h>
  12 #include <TH2F.h>
  13 #include <TH3F.h>
  14 #include <TParticle.h>
  15 #include <TRandom.h>
  16 #include <TNtuple.h>
  17 #include <TObjString.h>
  18 #include <TF1.h>
  19
  20 #include <AliLog.h>
  21 #include <AliESDVertex.h>
  22 #include <AliESDEvent.h>
  23 #include <AliStack.h>
  24 #include <AliHeader.h>
  25 #include <AliGenEventHeader.h>
  26 #include <AliMultiplicity.h>
  27 #include <AliAnalysisManager.h>
  28 #include <AliMCEventHandler.h>
  29 #include <AliMCEvent.h>
  30 #include <AliESDInputHandler.h>
  31
  32 #include "AliESDtrackCuts.h"
  33 #include "AliPWG0Helper.h"
  34 #include "multiplicity/AliMultiplicityCorrection.h"
  35 #include "AliCorrection.h"
  36 #include "AliCorrectionMatrix3D.h"
  37
  38 ClassImp(AliMultiplicityTask)
  39
  40 AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
  41   AliAnalysisTask("AliMultiplicityTask", ""),
  42   fESD(0),
  43   fOption(opt),
  44   fAnalysisMode(AliPWG0Helper::kSPD),
  45   fTrigger(AliPWG0Helper::kMB1),
  46   fDeltaPhiCut(-1),
  47   fReadMC(kFALSE),
  48   fUseMCVertex(kFALSE),
  49   fMultiplicity(0),
  50   fEsdTrackCuts(0),
  51   fSystSkipParticles(kFALSE),
  52   fSelectProcessType(0),
  53   fParticleSpecies(0),
  54   fdNdpT(0),
  55   fPtSpectrum(0),
  56   fOutput(0)
  57 {
  58   //
  59   // Constructor. Initialization of pointers
  60   //
  61
  62   for (Int_t i = 0; i<8; i++)
  63     fParticleCorrection[i] = 0;
  64
  65   // Define input and output slots here
  66   DefineInput(0, TChain::Class());
  67   DefineOutput(0, TList::Class());
  68 }
  69
  70 AliMultiplicityTask::~AliMultiplicityTask()
  71 {
  72   //
  73   // Destructor
  74   //
  75
  76   // histograms are in the output list and deleted when the output
  77   // list is deleted by the TSelector dtor
  78
  79   if (fOutput) {
  80     delete fOutput;
  81     fOutput = 0;
  82   }
  83 }
  84
  85 //________________________________________________________________________
  86 void AliMultiplicityTask::ConnectInputData(Option_t *)
  87 {
  88   // Connect ESD
  89   // Called once
  90
  91   Printf("AliMultiplicityTask::ConnectInputData called");
  92
  93   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  94   if (!tree) {
  95     Printf("ERROR: Could not read tree from input slot 0");
  96   } else {
  97     // Disable all branches and enable only the needed ones
  98     /*
  99     tree->SetBranchStatus("*", 0);
 100
 101     tree->SetBranchStatus("AliESDHeader*", 1);
 102     tree->SetBranchStatus("*Vertex*", 1);
 103
 104     if (fAnalysisMode == AliPWG0Helper::kSPD) {
 105       tree->SetBranchStatus("AliMultiplicity*", 1);
 106     }
 107
 108     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
 109       //AliESDtrackCuts::EnableNeededBranches(tree);
 110       tree->SetBranchStatus("Tracks*", 1);
 111     }
 112     */
 113
 114     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 115
 116     if (!esdH) {
 117       Printf("ERROR: Could not get ESDInputHandler");
 118     } else
 119       fESD = esdH->GetEvent();
 120   }
 121
 122   // disable info messages of AliMCEvent (per event)
 123   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 124 }
 125
 126 void AliMultiplicityTask::CreateOutputObjects()
 127 {
 128   // create result objects and add to output list
 129
 130   fOutput = new TList;
 131   fOutput->SetOwner();
 132
 133   fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 134   fOutput->Add(fMultiplicity);
 135
 136   fdNdpT = new TH1F("fdNdpT", "fdNdpT", 1000, 0, 10);
 137   fdNdpT->Sumw2();
 138   fOutput->Add(fdNdpT);
 139
 140   if (fOption.Contains("skip-particles"))
 141   {
 142     fSystSkipParticles = kTRUE;
 143     AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
 144   }
 145
 146   if (fOption.Contains("particle-efficiency"))
 147     for (Int_t i = 0; i<8; i++)
 148     {
 149       fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
 150       fOutput->Add(fParticleCorrection[i]);
 151     }
 152
 153   if (fOption.Contains("only-process-type-nd"))
 154     fSelectProcessType = 1;
 155
 156   if (fOption.Contains("only-process-type-sd"))
 157     fSelectProcessType = 2;
 158
 159   if (fOption.Contains("only-process-type-dd"))
 160     fSelectProcessType = 3;
 161
 162   if (fSelectProcessType != 0)
 163     AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
 164
 165   if (fOption.Contains("pt-spectrum-hist"))
 166   {
 167     TFile* file = TFile::Open("ptspectrum_fit.root");
 168     if (file)
 169     {
 170       TString subStr(fOption(fOption.Index("pt-spectrum")+17, 3));
 171       TString histName(Form("ptspectrum_%s", subStr.Data()));
 172       AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
 173       fPtSpectrum = (TH1D*) file->Get(histName);
 174       if (!fPtSpectrum)
 175         AliError("Histogram not found");
 176     }
 177     else
 178       AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
 179   }
 180
 181   if (fOption.Contains("pt-spectrum-func"))
 182   {
 183     if (fPtSpectrum)
 184     {
 185       Printf("Using function for systematic p_t study");
 186     }
 187     else
 188     {
 189       Printf("ERROR: Could not find function for systematic p_t study");
 190       fPtSpectrum = new TH1D("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
 191       fPtSpectrum->SetBinContent(1, 1);
 192     }
 193   }
 194
 195   if (fPtSpectrum)
 196     Printf("WARNING: Systematic study enabled. Pt spectrum will be modified");
 197
 198   if (fOption.Contains("particle-species")) {
 199     fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "vtx:Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec:nolabel:doublePrim:doubleCount");
 200     fOutput->Add(fParticleSpecies);
 201   }
 202
 203   // TODO set seed for random generator
 204 }
 205
 206 void AliMultiplicityTask::Exec(Option_t*)
 207 {
 208   // process the event
 209
 210   // Check prerequisites
 211   if (!fESD)
 212   {
 213     AliDebug(AliLog::kError, "ESD not available");
 214     return;
 215   }
 216
 217   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
 218   //Printf("%lld", fESD->GetTriggerMask());
 219
 220   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 221   Bool_t eventVertex = (vtxESD != 0);
 222
 223   Double_t vtx[3];
 224   if (vtxESD)
 225     vtxESD->GetXYZ(vtx);
 226
 227   // post the data already here
 228   PostData(0, fOutput);
 229
 230   //const Float_t kPtCut = 0.3;
 231
 232   // create list of (label, eta) tuples
 233   Int_t inputCount = 0;
 234   Int_t* labelArr = 0;
 235   Float_t* etaArr = 0;
 236   if (fAnalysisMode == AliPWG0Helper::kSPD)
 237   {
 238     // get tracklets
 239     const AliMultiplicity* mult = fESD->GetMultiplicity();
 240     if (!mult)
 241     {
 242       AliDebug(AliLog::kError, "AliMultiplicity not available");
 243       return;
 244     }
 245
 246     labelArr = new Int_t[mult->GetNumberOfTracklets()];
 247     etaArr = new Float_t[mult->GetNumberOfTracklets()];
 248
 249     // get multiplicity from ITS tracklets
 250     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 251     {
 252       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 253
 254       Float_t deltaPhi = mult->GetDeltaPhi(i);
 255
 256       if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
 257         continue;
 258
 259       etaArr[inputCount] = mult->GetEta(i);
 260       if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
 261       {
 262         labelArr[inputCount] = mult->GetLabel(i, 0);
 263       }
 264       else
 265         labelArr[inputCount] = -1;
 266
 267       ++inputCount;
 268     }
 269   }
 270   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
 271   {
 272     if (!fEsdTrackCuts)
 273     {
 274       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 275       return;
 276     }
 277
 278     if (vtxESD)
 279     {
 280       // get multiplicity from ESD tracks
 281       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
 282       Int_t nGoodTracks = list->GetEntries();
 283
 284       labelArr = new Int_t[nGoodTracks];
 285       etaArr = new Float_t[nGoodTracks];
 286
 287       // loop over esd tracks
 288       for (Int_t i=0; i<nGoodTracks; i++)
 289       {
 290         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 291         if (!esdTrack)
 292         {
 293           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 294           continue;
 295         }
 296
 297         etaArr[inputCount] = esdTrack->Eta();
 298         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 299         ++inputCount;
 300       }
 301
 302       delete list;
 303     }
 304   }
 305
 306   // eta range for nMCTracksSpecies, nESDTracksSpecies
 307   Float_t etaRange = -1;
 308   switch (fAnalysisMode) {
 309     case AliPWG0Helper::kInvalid: break;
 310     case AliPWG0Helper::kTPC:
 311     case AliPWG0Helper::kTPCITS:
 312         etaRange = 1.0; break;
 313     case AliPWG0Helper::kSPD: etaRange = 1.0; break;
 314     default: break;
 315   }
 316
 317   if (!fReadMC) // Processing of ESD information
 318   {
 319     if (eventTriggered && eventVertex)
 320     {
 321       Int_t nESDTracks05 = 0;
 322       Int_t nESDTracks10 = 0;
 323       Int_t nESDTracks14 = 0;
 324
 325       for (Int_t i=0; i<inputCount; ++i)
 326       {
 327         Float_t eta = etaArr[i];
 328
 329         if (TMath::Abs(eta) < 0.5)
 330           nESDTracks05++;
 331
 332         if (TMath::Abs(eta) < 1.0)
 333           nESDTracks10++;
 334
 335         if (TMath::Abs(eta) < 1.4)
 336           nESDTracks14++;
 337       }
 338
 339       fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 340     }
 341   }
 342   else if (fReadMC)   // Processing of MC information
 343   {
 344     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 345     if (!eventHandler) {
 346       Printf("ERROR: Could not retrieve MC event handler");
 347       return;
 348     }
 349
 350     AliMCEvent* mcEvent = eventHandler->MCEvent();
 351     if (!mcEvent) {
 352       Printf("ERROR: Could not retrieve MC event");
 353       return;
 354     }
 355
 356     AliStack* stack = mcEvent->Stack();
 357     if (!stack)
 358     {
 359       AliDebug(AliLog::kError, "Stack not available");
 360       return;
 361     }
 362
 363     AliHeader* header = mcEvent->Header();
 364     if (!header)
 365     {
 366       AliDebug(AliLog::kError, "Header not available");
 367       return;
 368     }
 369
 370     if (fUseMCVertex)
 371     {
 372       Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 373       // get the MC vertex
 374       AliGenEventHeader* genHeader = header->GenEventHeader();
 375       TArrayF vtxMC(3);
 376       genHeader->PrimaryVertex(vtxMC);
 377
 378       vtx[2] = vtxMC[2];
 379     }
 380
 381     // get process information
 382     AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(header);
 383
 384     Bool_t processEvent = kTRUE;
 385     if (fSelectProcessType > 0)
 386     {
 387       processEvent = kFALSE;
 388
 389       // non diffractive
 390       if (fSelectProcessType == 1 && processType == AliPWG0Helper::kND)
 391         processEvent = kTRUE;
 392
 393       // single diffractive
 394       if (fSelectProcessType == 2 && processType == AliPWG0Helper::kSD)
 395         processEvent = kTRUE;
 396
 397       // double diffractive
 398       if (fSelectProcessType == 3 && processType == AliPWG0Helper::kDD)
 399         processEvent = kTRUE;
 400
 401       if (!processEvent)
 402         Printf("Skipping this event, because it is not of the requested process type (%d)", (Int_t) processType);
 403     }
 404
 405     if (processEvent)
 406     {
 407       // get the MC vertex
 408       AliGenEventHeader* genHeader = header->GenEventHeader();
 409       TArrayF vtxMC(3);
 410       genHeader->PrimaryVertex(vtxMC);
 411
 412       // get number of tracks from MC
 413       // loop over mc particles
 414       Int_t nPrim  = stack->GetNprimary();
 415       Int_t nMCPart = stack->GetNtrack();
 416
 417       // tracks in different eta ranges
 418       Int_t nMCTracks05 = 0;
 419       Int_t nMCTracks10 = 0;
 420       Int_t nMCTracks14 = 0;
 421       Int_t nMCTracksAll = 0;
 422
 423       // tracks per particle species, in |eta| < 2 (systematic study)
 424       Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
 425       for (Int_t i = 0; i<4; ++i)
 426         nMCTracksSpecies[i] = 0;
 427
 428       for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 429       {
 430         AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 431
 432         TParticle* particle = stack->Particle(iMc);
 433
 434         if (!particle)
 435         {
 436           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 437           continue;
 438         }
 439
 440         Bool_t debug = kFALSE;
 441         if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
 442         {
 443           //printf("%d) DROPPED ", iMC);
 444           //particle->Print();
 445           continue;
 446         }
 447
 448         //printf("%d) OK      ", iMC);
 449         //particle->Print();
 450
 451         //if (particle->Pt() < kPtCut)
 452         //  continue;
 453
 454         Int_t particleWeight = 1;
 455
 456         Float_t pt = particle->Pt();
 457
 458         // in case of systematic study, weight according to the change of the pt spectrum
 459         // it cannot be just multiplied because we cannot count "half of a particle"
 460         // instead a random generator decides if the particle is counted twice (if value > 1)
 461         // or not (if value < 0)
 462         if (fPtSpectrum)
 463         {
 464           Int_t bin = fPtSpectrum->FindBin(pt);
 465           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 466           {
 467             Float_t factor = fPtSpectrum->GetBinContent(bin);
 468             if (factor > 0)
 469             {
 470               Float_t random = gRandom->Uniform();
 471               if (factor > 1 && random < factor - 1)
 472               {
 473                 particleWeight = 2;
 474               }
 475               else if (factor < 1 && random < 1 - factor)
 476                 particleWeight = 0;
 477             }
 478           }
 479         }
 480
 481         //Printf("MC weight is: %d", particleWeight);
 482
 483         if (TMath::Abs(particle->Eta()) < 0.5)
 484           nMCTracks05 += particleWeight;
 485
 486         if (TMath::Abs(particle->Eta()) < 1.0)
 487           nMCTracks10 += particleWeight;
 488
 489         if (TMath::Abs(particle->Eta()) < 1.4)
 490           nMCTracks14 += particleWeight;
 491
 492         nMCTracksAll += particleWeight;
 493
 494         if (particle->Pt() > 0 && TMath::Abs(particle->Eta()) < 1.0)
 495           fdNdpT->Fill(particle->Pt(), 1.0 / particle->Pt());
 496
 497         if (fParticleCorrection[0] || fParticleSpecies)
 498         {
 499           Int_t id = -1;
 500           switch (TMath::Abs(particle->GetPdgCode()))
 501           {
 502             case 211: id = 0; break;
 503             case 321: id = 1; break;
 504             case 2212: id = 2; break;
 505             default: id = 3; break;
 506           }
 507
 508           if (TMath::Abs(particle->Eta()) < etaRange)
 509             nMCTracksSpecies[id]++;
 510
 511           if (fParticleCorrection[id])
 512             fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
 513         }
 514       } // end of mc particle
 515
 516       fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
 517
 518       if (eventTriggered && eventVertex)
 519       {
 520         Int_t nESDTracks05 = 0;
 521         Int_t nESDTracks10 = 0;
 522         Int_t nESDTracks14 = 0;
 523
 524         // tracks per particle species, in |eta| < 2 (systematic study)
 525         Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
 526         for (Int_t i = 0; i<7; ++i)
 527           nESDTracksSpecies[i] = 0;
 528
 529         Bool_t* foundPrimaries = new Bool_t[nPrim];   // to prevent double counting
 530         for (Int_t i=0; i<nPrim; i++)
 531           foundPrimaries[i] = kFALSE;
 532
 533         Bool_t* foundPrimaries2 = new Bool_t[nPrim];   // to prevent double counting
 534         for (Int_t i=0; i<nPrim; i++)
 535           foundPrimaries2[i] = kFALSE;
 536
 537         Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
 538         for (Int_t i=0; i<nMCPart; i++)
 539           foundTracks[i] = kFALSE;
 540
 541         for (Int_t i=0; i<inputCount; ++i)
 542         {
 543           Float_t eta = etaArr[i];
 544           Int_t label = labelArr[i];
 545
 546           Int_t particleWeight = 1;
 547
 548           // systematic study: 5% lower efficiency
 549           if (fSystSkipParticles && (gRandom->Uniform() < 0.05))
 550             continue;
 551
 552           // in case of systematic study, weight according to the change of the pt spectrum
 553           if (fPtSpectrum)
 554           {
 555             TParticle* mother = 0;
 556
 557             // preserve label for later
 558             Int_t labelCopy = label;
 559             if (labelCopy >= 0)
 560               labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
 561             if (labelCopy >= 0)
 562               mother = stack->Particle(labelCopy);
 563
 564             // in case of pt study we do not count particles w/o label, because they cannot be scaled
 565             if (!mother)
 566               continue;
 567
 568             // it cannot be just multiplied because we cannot count "half of a particle"
 569             // instead a random generator decides if the particle is counted twice (if value > 1)
 570             // or not (if value < 0)
 571             Int_t bin = fPtSpectrum->FindBin(mother->Pt());
 572             if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 573             {
 574               Float_t factor = fPtSpectrum->GetBinContent(bin);
 575               if (factor > 0)
 576               {
 577                 Float_t random = gRandom->Uniform();
 578                 if (factor > 1 && random < factor - 1)
 579                 {
 580                   particleWeight = 2;
 581                 }
 582                 else if (factor < 1 && random < 1 - factor)
 583                   particleWeight = 0;
 584               }
 585             }
 586           }
 587
 588           //Printf("ESD weight is: %d", particleWeight);
 589
 590           if (TMath::Abs(eta) < 0.5)
 591             nESDTracks05 += particleWeight;
 592
 593           if (TMath::Abs(eta) < 1.0)
 594             nESDTracks10 += particleWeight;
 595
 596           if (TMath::Abs(eta) < 1.4)
 597             nESDTracks14 += particleWeight;
 598
 599           if (fParticleSpecies)
 600           {
 601             Int_t motherLabel = -1;
 602             TParticle* mother = 0;
 603
 604             // find mother
 605             if (label >= 0)
 606               motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
 607             if (motherLabel >= 0)
 608               mother = stack->Particle(motherLabel);
 609
 610             if (!mother)
 611             {
 612               // count tracks that did not have a label
 613               if (TMath::Abs(eta) < etaRange)
 614                 nESDTracksSpecies[4]++;
 615             }
 616             else
 617             {
 618               // get particle type (pion, proton, kaon, other) of mother
 619               Int_t idMother = -1;
 620               switch (TMath::Abs(mother->GetPdgCode()))
 621               {
 622                 case 211: idMother = 0; break;
 623                 case 321: idMother = 1; break;
 624                 case 2212: idMother = 2; break;
 625                 default: idMother = 3; break;
 626               }
 627
 628               // double counting is ok for particle ratio study
 629               if (TMath::Abs(eta) < etaRange)
 630                 nESDTracksSpecies[idMother]++;
 631
 632               // double counting is not ok for efficiency study
 633
 634               // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
 635               if (foundTracks[label])
 636               {
 637                 if (TMath::Abs(eta) < etaRange)
 638                   nESDTracksSpecies[6]++;
 639               }
 640               else
 641               {
 642                 foundTracks[label] = kTRUE;
 643
 644                 // particle (primary) already counted?
 645                 if (foundPrimaries[motherLabel])
 646                 {
 647                   if (TMath::Abs(eta) < etaRange)
 648                     nESDTracksSpecies[5]++;
 649                 }
 650                 else
 651                   foundPrimaries[motherLabel] = kTRUE;
 652               }
 653             }
 654           }
 655
 656           if (fParticleCorrection[0])
 657           {
 658             if (label >= 0 && stack->IsPhysicalPrimary(label))
 659             {
 660               TParticle* particle = stack->Particle(label);
 661
 662               // get particle type (pion, proton, kaon, other)
 663               Int_t id = -1;
 664               switch (TMath::Abs(particle->GetPdgCode()))
 665               {
 666                 case 211: id = 0; break;
 667                 case 321: id = 1; break;
 668                 case 2212: id = 2; break;
 669                 default: id = 3; break;
 670               }
 671
 672               // todo check if values are not completely off??
 673
 674               // particle (primary) already counted?
 675               if (!foundPrimaries2[label])
 676               {
 677                 foundPrimaries2[label] = kTRUE;
 678                 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
 679               }
 680             }
 681           }
 682         }
 683
 684         if (fParticleCorrection[0])
 685         {
 686           // if the particle decays/stops before this radius we do not see it
 687           // 8cm larger than SPD layer 2
 688           // 123cm TPC radius where a track has about 50 clusters (cut limit)
 689           const Float_t endRadius = (fAnalysisMode == AliPWG0Helper::kSPD) ? 8. : 123;
 690
 691           // loop over all primaries that have not been found
 692           for (Int_t i=0; i<nPrim; i++)
 693           {
 694             // already found
 695             if (foundPrimaries2[i])
 696               continue;
 697
 698             TParticle* particle = 0;
 699             TClonesArray* trackrefs = 0;
 700             mcEvent->GetParticleAndTR(i, particle, trackrefs);
 701
 702             // true primary and charged
 703             if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
 704               continue;
 705
 706             //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
 707             if (TMath::Abs(particle->Eta()) > 3)
 708               continue;
 709
 710             // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
 711
 712             // get particle type (pion, proton, kaon, other)
 713             Int_t id = -1;
 714             switch (TMath::Abs(particle->GetPdgCode()))
 715             {
 716               case 211: id = 4; break;
 717               case 321: id = 5; break;
 718               case 2212: id = 6; break;
 719               default: id = 7; break;
 720             }
 721
 722             if (!fParticleCorrection[id])
 723               continue;
 724
 725             // get last track reference
 726             AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
 727
 728             if (!trackref)
 729             {
 730               Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
 731               particle->Print();
 732               continue;
 733             }
 734
 735             // particle in tracking volume long enough...
 736             if (trackref->R() > endRadius)
 737               continue;
 738
 739             if (particle->GetLastDaughter() >= 0)
 740             {
 741               Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
 742               //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
 743               if (uID == kPDecay)
 744               {
 745                 // decayed
 746
 747                 Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
 748                 particle->Print();
 749                 Printf("Daughers:");
 750                 for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
 751                   stack->Particle(d)->Print();
 752                 Printf("");
 753
 754                 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
 755                 continue;
 756               }
 757             }
 758
 759             if (trackref->DetectorId() == -1)
 760             {
 761               // stopped
 762               Printf("Particle %d stopped at %f:", i, trackref->R());
 763               particle->Print();
 764               Printf("");
 765
 766               fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
 767               continue;
 768             }
 769
 770             Printf("Particle %d simply not tracked", i);
 771             particle->Print();
 772             Printf("");
 773           }
 774         }
 775
 776         delete[] foundTracks;
 777         delete[] foundPrimaries;
 778         delete[] foundPrimaries2;
 779
 780         if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
 781         {
 782             TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 783             printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
 784         }
 785
 786         // fill response matrix using vtxMC (best guess)
 787         fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks14,  nMCTracksAll,  nESDTracks05,  nESDTracks10, nESDTracks14);
 788
 789         fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 790
 791         if (fParticleSpecies)
 792           fParticleSpecies->Fill(vtxMC[2], nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
 793       }
 794     }
 795   }
 796
 797   if (etaArr)
 798     delete[] etaArr;
 799   if (labelArr)
 800     delete[] labelArr;
 801 }
 802
 803 void AliMultiplicityTask::Terminate(Option_t *)
 804 {
 805   // The Terminate() function is the last function to be called during
 806   // a query. It always runs on the client, it can be used to present
 807   // the results graphically or save the results to file.
 808
 809   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 810   if (!fOutput) {
 811     Printf("ERROR: fOutput not available");
 812     return;
 813   }
 814
 815   fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
 816   for (Int_t i = 0; i < 8; ++i)
 817     fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
 818   fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
 819
 820   fdNdpT = dynamic_cast<TH1*> (fOutput->FindObject("fdNdpT"));
 821
 822   if (!fMultiplicity)
 823   {
 824     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
 825     return;
 826   }
 827
 828   TFile* file = TFile::Open("multiplicity.root", "RECREATE");
 829
 830   fMultiplicity->SaveHistograms();
 831   for (Int_t i = 0; i < 8; ++i)
 832     if (fParticleCorrection[i])
 833       fParticleCorrection[i]->SaveHistograms();
 834   if (fParticleSpecies)
 835     fParticleSpecies->Write();
 836   if (fdNdpT)
 837     fdNdpT->Write();
 838
 839   TObjString option(fOption);
 840   option.Write();
 841
 842   file->Close();
 843
 844   Printf("Written result to multiplicity.root");
 845 }