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