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 <TH2F.h>
  12 #include <TH3F.h>
  13 #include <TParticle.h>
  14 #include <TRandom.h>
  15 #include <TNtuple.h>
  16 #include <TObjString.h>
  17 #include <TF1.h>
  18
  19 #include <AliLog.h>
  20 #include <AliESDVertex.h>
  21 #include <AliESDEvent.h>
  22 #include <AliStack.h>
  23 #include <AliHeader.h>
  24 #include <AliGenEventHeader.h>
  25 #include <AliMultiplicity.h>
  26 #include <AliAnalysisManager.h>
  27 #include <AliMCEventHandler.h>
  28 #include <AliMCEvent.h>
  29 #include <AliESDInputHandler.h>
  30
  31 #include "AliESDtrackCuts.h"
  32 #include "AliPWG0Helper.h"
  33 #include "dNdEta/AliMultiplicityCorrection.h"
  34 #include "AliCorrection.h"
  35 #include "AliCorrectionMatrix3D.h"
  36
  37 ClassImp(AliMultiplicityTask)
  38
  39 AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
  40   AliAnalysisTask("AliMultiplicityTask", ""),
  41   fESD(0),
  42   fOption(opt),
  43   fAnalysisMode(AliPWG0Helper::kSPD),
  44   fReadMC(kFALSE),
  45   fMultiplicity(0),
  46   fEsdTrackCuts(0),
  47   fSystSkipParticles(kFALSE),
  48   fSelectProcessType(0),
  49   fParticleSpecies(0),
  50   fPtSpectrum(0),
  51   fOutput(0)
  52 {
  53   //
  54   // Constructor. Initialization of pointers
  55   //
  56
  57   for (Int_t i = 0; i<4; i++)
  58     fParticleCorrection[i] = 0;
  59
  60   // Define input and output slots here
  61   DefineInput(0, TChain::Class());
  62   DefineOutput(0, TList::Class());
  63 }
  64
  65 AliMultiplicityTask::~AliMultiplicityTask()
  66 {
  67   //
  68   // Destructor
  69   //
  70
  71   // histograms are in the output list and deleted when the output
  72   // list is deleted by the TSelector dtor
  73
  74   if (fOutput) {
  75     delete fOutput;
  76     fOutput = 0;
  77   }
  78 }
  79
  80 //________________________________________________________________________
  81 void AliMultiplicityTask::ConnectInputData(Option_t *)
  82 {
  83   // Connect ESD
  84   // Called once
  85
  86   Printf("AliMultiplicityTask::ConnectInputData called");
  87
  88   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  89   if (!tree) {
  90     Printf("ERROR: Could not read tree from input slot 0");
  91   } else {
  92     // Disable all branches and enable only the needed ones
  93     //tree->SetBranchStatus("*", 0);
  94
  95     tree->SetBranchStatus("fTriggerMask", 1);
  96     tree->SetBranchStatus("fSPDVertex*", 1);
  97
  98     if (fAnalysisMode == AliPWG0Helper::kSPD)
  99       tree->SetBranchStatus("fSPDMult*", 1);
 100
 101     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
 102       AliESDtrackCuts::EnableNeededBranches(tree);
 103       tree->SetBranchStatus("fTracks.fLabel", 1);
 104     }
 105
 106     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 107
 108     if (!esdH) {
 109       Printf("ERROR: Could not get ESDInputHandler");
 110     } else
 111       fESD = esdH->GetEvent();
 112   }
 113 }
 114
 115 void AliMultiplicityTask::CreateOutputObjects()
 116 {
 117   // create result objects and add to output list
 118
 119   fOutput = new TList;
 120   fOutput->SetOwner();
 121
 122   fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 123   fOutput->Add(fMultiplicity);
 124
 125   if (fOption.Contains("skip-particles"))
 126   {
 127     fSystSkipParticles = kTRUE;
 128     AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
 129   }
 130
 131   if (fOption.Contains("particle-efficiency"))
 132     for (Int_t i = 0; i<4; i++)
 133     {
 134       fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
 135       fOutput->Add(fParticleCorrection[i]);
 136     }
 137
 138   if (fOption.Contains("only-process-type-nd"))
 139     fSelectProcessType = 1;
 140
 141   if (fOption.Contains("only-process-type-sd"))
 142     fSelectProcessType = 2;
 143
 144   if (fOption.Contains("only-process-type-dd"))
 145     fSelectProcessType = 3;
 146
 147   if (fSelectProcessType != 0)
 148     AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
 149
 150   if (fOption.Contains("pt-spectrum-hist"))
 151   {
 152     TFile* file = TFile::Open("ptspectrum_fit.root");
 153     if (file)
 154     {
 155       TString subStr(fOption(fOption.Index("pt-spectrum")+17, 3));
 156       TString histName(Form("ptspectrum_%s", subStr.Data()));
 157       AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
 158       fPtSpectrum = (TH1*) file->Get(histName);
 159       if (!fPtSpectrum)
 160         AliError("Histogram not found");
 161     }
 162     else
 163       AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
 164
 165     if (fPtSpectrum)
 166       AliInfo("WARNING: Systematic study enabled. Pt spectrum will be modified");
 167   }
 168
 169   if (fOption.Contains("pt-spectrum-func"))
 170   {
 171     if (fPtSpectrum)
 172     {
 173       Printf("Using function from input list for systematic p_t study");
 174     }
 175     else
 176     {
 177       fPtSpectrum = new TH1F("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
 178       fPtSpectrum->SetBinContent(1, 1);
 179     }
 180
 181     if (fPtSpectrum)
 182       AliInfo("WARNING: Systematic study enabled. Pt spectrum will be modified");
 183   }
 184
 185   if (fOption.Contains("particle-species")) {
 186     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");
 187     fOutput->Add(fParticleSpecies);
 188   }
 189
 190   // TODO set seed for random generator
 191 }
 192
 193 void AliMultiplicityTask::Exec(Option_t*)
 194 {
 195   // process the event
 196
 197   // Check prerequisites
 198   if (!fESD)
 199   {
 200     AliDebug(AliLog::kError, "ESD not available");
 201     return;
 202   }
 203
 204   //MB1 definition
 205   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB1);
 206   // only FASTOR
 207   //Bool_t eventTriggered = fESD->GetTriggerMask() & 32;
 208
 209   Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD->GetVertex());
 210
 211   // get the ESD vertex
 212   const AliESDVertex* vtxESD = fESD->GetVertex();
 213   Double_t vtx[3];
 214   vtxESD->GetXYZ(vtx);
 215
 216   // post the data already here
 217   PostData(0, fOutput);
 218
 219   //const Float_t kPtCut = 0.3;
 220
 221   // create list of (label, eta) tuples
 222   Int_t inputCount = 0;
 223   Int_t* labelArr = 0;
 224   Float_t* etaArr = 0;
 225   if (fAnalysisMode == AliPWG0Helper::kSPD)
 226   {
 227     // get tracklets
 228     const AliMultiplicity* mult = fESD->GetMultiplicity();
 229     if (!mult)
 230     {
 231       AliDebug(AliLog::kError, "AliMultiplicity not available");
 232       return;
 233     }
 234
 235     labelArr = new Int_t[mult->GetNumberOfTracklets()];
 236     etaArr = new Float_t[mult->GetNumberOfTracklets()];
 237
 238     // get multiplicity from ITS tracklets
 239     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 240     {
 241       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 242
 243       // This removes non-tracklets in PDC06 data. Very bad solution. New solution is implemented for newer data. Keeping this for compatibility.
 244       if (mult->GetDeltaPhi(i) < -1000)
 245         continue;
 246
 247       etaArr[inputCount] = mult->GetEta(i);
 248       // TODO add second label array
 249       labelArr[inputCount] = mult->GetLabel(i, 0);
 250       ++inputCount;
 251     }
 252   }
 253   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
 254   {
 255     if (!fEsdTrackCuts)
 256     {
 257       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 258       return;
 259     }
 260
 261     // get multiplicity from ESD tracks
 262     TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
 263     Int_t nGoodTracks = list->GetEntries();
 264
 265     labelArr = new Int_t[nGoodTracks];
 266     etaArr = new Float_t[nGoodTracks];
 267
 268     // loop over esd tracks
 269     for (Int_t i=0; i<nGoodTracks; i++)
 270     {
 271       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 272       if (!esdTrack)
 273       {
 274         AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 275         continue;
 276       }
 277
 278       etaArr[inputCount] = esdTrack->Eta();
 279       labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 280       ++inputCount;
 281     }
 282
 283     delete list;
 284   }
 285
 286   // eta range for nMCTracksSpecies, nESDTracksSpecies
 287   Float_t etaRange = -1;
 288   switch (fAnalysisMode) {
 289     case AliPWG0Helper::kInvalid: break;
 290     case AliPWG0Helper::kTPC:
 291     case AliPWG0Helper::kTPCITS:
 292         etaRange = 0.9; break;
 293     case AliPWG0Helper::kSPD: etaRange = 2.0; break;
 294   }
 295
 296   if (!fReadMC) // Processing of ESD information
 297   {
 298     if (eventTriggered && eventVertex)
 299     {
 300       Int_t nESDTracks05 = 0;
 301       Int_t nESDTracks09 = 0;
 302       Int_t nESDTracks15 = 0;
 303       Int_t nESDTracks20 = 0;
 304
 305       for (Int_t i=0; i<inputCount; ++i)
 306       {
 307         Float_t eta = etaArr[i];
 308
 309         if (TMath::Abs(eta) < 0.5)
 310           nESDTracks05++;
 311
 312         if (TMath::Abs(eta) < 0.9)
 313           nESDTracks09++;
 314
 315         if (TMath::Abs(eta) < 1.5)
 316           nESDTracks15++;
 317
 318         if (TMath::Abs(eta) < 2.0)
 319           nESDTracks20++;
 320       }
 321
 322       fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
 323     }
 324   }
 325   else if (fReadMC)   // Processing of MC information
 326   {
 327     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 328     if (!eventHandler) {
 329       Printf("ERROR: Could not retrieve MC event handler");
 330       return;
 331     }
 332
 333     AliMCEvent* mcEvent = eventHandler->MCEvent();
 334     if (!mcEvent) {
 335       Printf("ERROR: Could not retrieve MC event");
 336       return;
 337     }
 338
 339     AliStack* stack = mcEvent->Stack();
 340     if (!stack)
 341     {
 342       AliDebug(AliLog::kError, "Stack not available");
 343       return;
 344     }
 345
 346     AliHeader* header = mcEvent->Header();
 347     if (!header)
 348     {
 349       AliDebug(AliLog::kError, "Header not available");
 350       return;
 351     }
 352
 353     Bool_t processEvent = kTRUE;
 354     if (fSelectProcessType > 0)
 355     {
 356       // getting process information; NB: this only works for Pythia
 357       Int_t processtype = AliPWG0Helper::GetPythiaEventProcessType(header);
 358
 359       processEvent = kFALSE;
 360
 361       // non diffractive
 362       if (fSelectProcessType == 1 && processtype !=92 && processtype !=93 && processtype != 94)
 363         processEvent = kTRUE;
 364
 365       // single diffractive
 366       if (fSelectProcessType == 2 && (processtype == 92 || processtype == 93))
 367         processEvent = kTRUE;
 368
 369       // double diffractive
 370       if (fSelectProcessType == 3 && processtype == 94)
 371         processEvent = kTRUE;
 372
 373       if (!processEvent)
 374         AliDebug(AliLog::kDebug, Form("Skipping this event, because it is not of the requested process type (%d)", processtype));
 375     }
 376
 377     // systematic study: 10% lower efficiency
 378     if (fSystSkipParticles && (gRandom->Uniform() < 0.1))
 379       processEvent = kFALSE;
 380
 381     if (processEvent)
 382     {
 383       // get the MC vertex
 384       AliGenEventHeader* genHeader = header->GenEventHeader();
 385       TArrayF vtxMC(3);
 386       genHeader->PrimaryVertex(vtxMC);
 387
 388       // get number of tracks from MC
 389       // loop over mc particles
 390       Int_t nPrim  = stack->GetNprimary();
 391       Int_t nMCPart = stack->GetNtrack();
 392
 393       // tracks in different eta ranges
 394       Int_t nMCTracks05 = 0;
 395       Int_t nMCTracks09 = 0;
 396       Int_t nMCTracks15 = 0;
 397       Int_t nMCTracks20 = 0;
 398       Int_t nMCTracksAll = 0;
 399
 400       // tracks per particle species, in |eta| < 2 (systematic study)
 401       Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
 402       for (Int_t i = 0; i<4; ++i)
 403         nMCTracksSpecies[i] = 0;
 404
 405       for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 406       {
 407         AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 408
 409         TParticle* particle = stack->Particle(iMc);
 410
 411         if (!particle)
 412         {
 413           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 414           continue;
 415         }
 416
 417         Bool_t debug = kFALSE;
 418         if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
 419         {
 420           //printf("%d) DROPPED ", iMC);
 421           //particle->Print();
 422           continue;
 423         }
 424
 425         //printf("%d) OK      ", iMC);
 426         //particle->Print();
 427
 428         //if (particle->Pt() < kPtCut)
 429         //  continue;
 430
 431         Int_t particleWeight = 1;
 432
 433         Float_t pt = particle->Pt();
 434
 435         // in case of systematic study, weight according to the change of the pt spectrum
 436         // it cannot be just multiplied because we cannot count "half of a particle"
 437         // instead a random generator decides if the particle is counted twice (if value > 1)
 438         // or not (if value < 0)
 439         if (fPtSpectrum)
 440         {
 441           Int_t bin = fPtSpectrum->FindBin(pt);
 442           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 443           {
 444             Float_t factor = fPtSpectrum->GetBinContent(bin);
 445             if (factor > 0)
 446             {
 447               Float_t random = gRandom->Uniform();
 448               if (factor > 1 && random < factor - 1)
 449               {
 450                 particleWeight = 2;
 451               }
 452               else if (factor < 1 && random < 1 - factor)
 453                 particleWeight = 0;
 454             }
 455           }
 456         }
 457
 458         //Printf("MC weight is: %d", particleWeight);
 459
 460         if (TMath::Abs(particle->Eta()) < 0.5)
 461           nMCTracks05 += particleWeight;
 462
 463         if (TMath::Abs(particle->Eta()) < 0.9)
 464           nMCTracks09 += particleWeight;
 465
 466         if (TMath::Abs(particle->Eta()) < 1.5)
 467           nMCTracks15 += particleWeight;
 468
 469         if (TMath::Abs(particle->Eta()) < 2.0)
 470           nMCTracks20 += particleWeight;
 471
 472         nMCTracksAll += particleWeight;
 473
 474         if (fParticleCorrection[0] || fParticleSpecies)
 475         {
 476           Int_t id = -1;
 477           switch (TMath::Abs(particle->GetPdgCode()))
 478           {
 479             case 211: id = 0; break;
 480             case 321: id = 1; break;
 481             case 2212: id = 2; break;
 482             default: id = 3; break;
 483           }
 484
 485           if (TMath::Abs(particle->Eta()) < etaRange)
 486             nMCTracksSpecies[id]++;
 487
 488           if (fParticleCorrection[id])
 489             fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
 490         }
 491       } // end of mc particle
 492
 493       fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, (Int_t) nMCTracks05, (Int_t) nMCTracks09, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
 494
 495       if (eventTriggered && eventVertex)
 496       {
 497         Int_t nESDTracks05 = 0;
 498         Int_t nESDTracks09 = 0;
 499         Int_t nESDTracks15 = 0;
 500         Int_t nESDTracks20 = 0;
 501
 502         // tracks per particle species, in |eta| < 2 (systematic study)
 503         Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
 504         for (Int_t i = 0; i<7; ++i)
 505           nESDTracksSpecies[i] = 0;
 506
 507         Bool_t* foundPrimaries = new Bool_t[nPrim];   // to prevent double counting
 508         for (Int_t i=0; i<nPrim; i++)
 509           foundPrimaries[i] = kFALSE;
 510
 511         Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
 512         for (Int_t i=0; i<nMCPart; i++)
 513           foundTracks[i] = kFALSE;
 514
 515         for (Int_t i=0; i<inputCount; ++i)
 516         {
 517           Float_t eta = etaArr[i];
 518           Int_t label = labelArr[i];
 519
 520           Int_t particleWeight = 1;
 521
 522           // in case of systematic study, weight according to the change of the pt spectrum
 523           if (fPtSpectrum)
 524           {
 525             TParticle* mother = 0;
 526
 527             // preserve label for later
 528             Int_t labelCopy = label;
 529             if (labelCopy >= 0)
 530               labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
 531             if (labelCopy >= 0)
 532               mother = stack->Particle(labelCopy);
 533
 534             // in case of pt study we do not count particles w/o label, because they cannot be scaled
 535             if (!mother)
 536               continue;
 537
 538             // it cannot be just multiplied because we cannot count "half of a particle"
 539             // instead a random generator decides if the particle is counted twice (if value > 1)
 540             // or not (if value < 0)
 541             Int_t bin = fPtSpectrum->FindBin(mother->Pt());
 542             if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 543             {
 544               Float_t factor = fPtSpectrum->GetBinContent(bin);
 545               if (factor > 0)
 546               {
 547                 Float_t random = gRandom->Uniform();
 548                 if (factor > 1 && random < factor - 1)
 549                 {
 550                   particleWeight = 2;
 551                 }
 552                 else if (factor < 1 && random < 1 - factor)
 553                   particleWeight = 0;
 554               }
 555             }
 556           }
 557
 558           //Printf("ESD weight is: %d", particleWeight);
 559
 560           if (TMath::Abs(eta) < 0.5)
 561             nESDTracks05 += particleWeight;
 562
 563           if (TMath::Abs(eta) < 0.9)
 564             nESDTracks09 += particleWeight;
 565
 566           if (TMath::Abs(eta) < 1.5)
 567             nESDTracks15 += particleWeight;
 568
 569           if (TMath::Abs(eta) < 2.0)
 570             nESDTracks20 += particleWeight;
 571
 572
 573           if (fParticleCorrection[0] || fParticleSpecies)
 574           {
 575             Int_t motherLabel = -1;
 576             TParticle* mother = 0;
 577
 578             // find mother
 579             if (label >= 0)
 580               motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
 581             if (motherLabel >= 0)
 582               mother = stack->Particle(motherLabel);
 583
 584             if (!mother)
 585             {
 586               // count tracks that did not have a label
 587               if (TMath::Abs(eta) < etaRange)
 588                 nESDTracksSpecies[4]++;
 589               continue;
 590             }
 591
 592             // get particle type (pion, proton, kaon, other)
 593             Int_t id = -1;
 594             switch (TMath::Abs(mother->GetPdgCode()))
 595             {
 596               case 211: id = 0; break;
 597               case 321: id = 1; break;
 598               case 2212: id = 2; break;
 599               default: id = 3; break;
 600             }
 601
 602             // double counting is ok for particle ratio study
 603             if (TMath::Abs(eta) < etaRange)
 604               nESDTracksSpecies[id]++;
 605
 606             // double counting is not ok for efficiency study
 607
 608             // 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)
 609             if (foundTracks[label])
 610             {
 611               if (TMath::Abs(eta) < etaRange)
 612                 nESDTracksSpecies[6]++;
 613               continue;
 614             }
 615             foundTracks[label] = kTRUE;
 616
 617             // particle (primary) already counted?
 618             if (foundPrimaries[motherLabel])
 619             {
 620               if (TMath::Abs(eta) < etaRange)
 621                 nESDTracksSpecies[5]++;
 622               continue;
 623             }
 624             foundPrimaries[motherLabel] = kTRUE;
 625
 626             if (fParticleCorrection[id])
 627               fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
 628           }
 629         }
 630
 631         delete[] foundTracks;
 632         delete[] foundPrimaries;
 633
 634         if ((Int_t) nMCTracks20 == 0 && nESDTracks20 > 3)
 635             printf("WARNING: Event has %d generated and %d reconstructed...\n", nMCTracks20, nESDTracks20);
 636
 637         // fill response matrix using vtxMC (best guess)
 638         fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks09,  nMCTracks15,  nMCTracks20,  nMCTracksAll,  nESDTracks05,  nESDTracks09,  nESDTracks15,  nESDTracks20);
 639
 640         fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
 641
 642         if (fParticleSpecies)
 643           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]);
 644       }
 645     }
 646   }
 647
 648   delete etaArr;
 649   delete labelArr;
 650 }
 651
 652 void AliMultiplicityTask::Terminate(Option_t *)
 653 {
 654   // The Terminate() function is the last function to be called during
 655   // a query. It always runs on the client, it can be used to present
 656   // the results graphically or save the results to file.
 657
 658   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 659   if (!fOutput) {
 660     Printf("ERROR: fOutput not available");
 661     return;
 662   }
 663
 664   fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
 665   for (Int_t i = 0; i < 4; ++i)
 666     fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
 667   fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
 668
 669   if (!fMultiplicity)
 670   {
 671     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
 672     return;
 673   }
 674
 675   TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
 676
 677   fMultiplicity->SaveHistograms();
 678   for (Int_t i = 0; i < 4; ++i)
 679     if (fParticleCorrection[i])
 680       fParticleCorrection[i]->SaveHistograms();
 681   if (fParticleSpecies)
 682     fParticleSpecies->Write();
 683
 684   TObjString option(fOption);
 685   option.Write();
 686
 687   file->Close();
 688
 689   Printf("Writting result to multiplicityMC.root");
 690 }