PWG0/dNdEta/AliMultiplicityMCSelector.cxx

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