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