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