fixing warnings
[u/mrichter/AliRoot.git] / 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::kTPC:
290     case AliPWG0Helper::kTPCITS:
291         etaRange = 0.9; break;
292     case AliPWG0Helper::kSPD: etaRange = 2.0; break;
293   }
294
295   if (!fReadMC) // Processing of ESD information
296   {
297     if (eventTriggered && eventVertex)
298     {
299       Int_t nESDTracks05 = 0;
300       Int_t nESDTracks09 = 0;
301       Int_t nESDTracks15 = 0;
302       Int_t nESDTracks20 = 0;
303
304       for (Int_t i=0; i<inputCount; ++i)
305       {
306         Float_t eta = etaArr[i];
307
308         if (TMath::Abs(eta) < 0.5)
309           nESDTracks05++;
310
311         if (TMath::Abs(eta) < 0.9)
312           nESDTracks09++;
313
314         if (TMath::Abs(eta) < 1.5)
315           nESDTracks15++;
316
317         if (TMath::Abs(eta) < 2.0)
318           nESDTracks20++;
319       }
320
321       fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
322     }
323   }
324   else if (fReadMC)   // Processing of MC information
325   {
326     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
327     if (!eventHandler) {
328       Printf("ERROR: Could not retrieve MC event handler");
329       return;
330     }
331
332     AliMCEvent* mcEvent = eventHandler->MCEvent();
333     if (!mcEvent) {
334       Printf("ERROR: Could not retrieve MC event");
335       return;
336     }
337
338     AliStack* stack = mcEvent->Stack();
339     if (!stack)
340     {
341       AliDebug(AliLog::kError, "Stack not available");
342       return;
343     }
344
345     AliHeader* header = mcEvent->Header();
346     if (!header)
347     {
348       AliDebug(AliLog::kError, "Header not available");
349       return;
350     }
351
352     Bool_t processEvent = kTRUE;
353     if (fSelectProcessType > 0)
354     {
355       // getting process information; NB: this only works for Pythia
356       Int_t processtype = AliPWG0Helper::GetPythiaEventProcessType(header);
357
358       processEvent = kFALSE;
359
360       // non diffractive
361       if (fSelectProcessType == 1 && processtype !=92 && processtype !=93 && processtype != 94)
362         processEvent = kTRUE;
363
364       // single diffractive
365       if (fSelectProcessType == 2 && (processtype == 92 || processtype == 93))
366         processEvent = kTRUE;
367
368       // double diffractive
369       if (fSelectProcessType == 3 && processtype == 94)
370         processEvent = kTRUE;
371
372       if (!processEvent)
373         AliDebug(AliLog::kDebug, Form("Skipping this event, because it is not of the requested process type (%d)", processtype));
374     }
375
376     // systematic study: 10% lower efficiency
377     if (fSystSkipParticles && (gRandom->Uniform() < 0.1))
378       processEvent = kFALSE;
379
380     if (processEvent)
381     {
382       // get the MC vertex
383       AliGenEventHeader* genHeader = header->GenEventHeader();
384       TArrayF vtxMC(3);
385       genHeader->PrimaryVertex(vtxMC);
386
387       // get number of tracks from MC
388       // loop over mc particles
389       Int_t nPrim  = stack->GetNprimary();
390       Int_t nMCPart = stack->GetNtrack();
391
392       // tracks in different eta ranges
393       Int_t nMCTracks05 = 0;
394       Int_t nMCTracks09 = 0;
395       Int_t nMCTracks15 = 0;
396       Int_t nMCTracks20 = 0;
397       Int_t nMCTracksAll = 0;
398
399       // tracks per particle species, in |eta| < 2 (systematic study)
400       Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
401       for (Int_t i = 0; i<4; ++i)
402         nMCTracksSpecies[i] = 0;
403
404       for (Int_t iMc = 0; iMc < nPrim; ++iMc)
405       {
406         AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
407
408         TParticle* particle = stack->Particle(iMc);
409
410         if (!particle)
411         {
412           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
413           continue;
414         }
415
416         Bool_t debug = kFALSE;
417         if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
418         {
419           //printf("%d) DROPPED ", iMC);
420           //particle->Print();
421           continue;
422         }
423
424         //printf("%d) OK      ", iMC);
425         //particle->Print();
426
427         //if (particle->Pt() < kPtCut)
428         //  continue;
429
430         Int_t particleWeight = 1;
431
432         Float_t pt = particle->Pt();
433
434         // in case of systematic study, weight according to the change of the pt spectrum
435         // it cannot be just multiplied because we cannot count "half of a particle"
436         // instead a random generator decides if the particle is counted twice (if value > 1)
437         // or not (if value < 0)
438         if (fPtSpectrum)
439         {
440           Int_t bin = fPtSpectrum->FindBin(pt);
441           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
442           {
443             Float_t factor = fPtSpectrum->GetBinContent(bin);
444             if (factor > 0)
445             {
446               Float_t random = gRandom->Uniform();
447               if (factor > 1 && random < factor - 1)
448               {
449                 particleWeight = 2;
450               }
451               else if (factor < 1 && random < 1 - factor)
452                 particleWeight = 0;
453             }
454           }
455         }
456
457         //Printf("MC weight is: %d", particleWeight);
458
459         if (TMath::Abs(particle->Eta()) < 0.5)
460           nMCTracks05 += particleWeight;
461
462         if (TMath::Abs(particle->Eta()) < 0.9)
463           nMCTracks09 += particleWeight;
464
465         if (TMath::Abs(particle->Eta()) < 1.5)
466           nMCTracks15 += particleWeight;
467
468         if (TMath::Abs(particle->Eta()) < 2.0)
469           nMCTracks20 += particleWeight;
470
471         nMCTracksAll += particleWeight;
472
473         if (fParticleCorrection[0] || fParticleSpecies)
474         {
475           Int_t id = -1;
476           switch (TMath::Abs(particle->GetPdgCode()))
477           {
478             case 211: id = 0; break;
479             case 321: id = 1; break;
480             case 2212: id = 2; break;
481             default: id = 3; break;
482           }
483
484           if (TMath::Abs(particle->Eta()) < etaRange)
485             nMCTracksSpecies[id]++;
486             
487           if (fParticleCorrection[id])
488             fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
489         }
490       } // end of mc particle
491
492       fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, (Int_t) nMCTracks05, (Int_t) nMCTracks09, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
493
494       if (eventTriggered && eventVertex)
495       {
496         Int_t nESDTracks05 = 0;
497         Int_t nESDTracks09 = 0;
498         Int_t nESDTracks15 = 0;
499         Int_t nESDTracks20 = 0;
500
501         // tracks per particle species, in |eta| < 2 (systematic study)
502         Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
503         for (Int_t i = 0; i<7; ++i)
504           nESDTracksSpecies[i] = 0;
505
506         Bool_t* foundPrimaries = new Bool_t[nPrim];   // to prevent double counting
507         for (Int_t i=0; i<nPrim; i++)
508           foundPrimaries[i] = kFALSE;
509
510         Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
511         for (Int_t i=0; i<nMCPart; i++)
512           foundTracks[i] = kFALSE;
513
514         for (Int_t i=0; i<inputCount; ++i)
515         {
516           Float_t eta = etaArr[i];
517           Int_t label = labelArr[i];
518
519           Int_t particleWeight = 1;
520
521           // in case of systematic study, weight according to the change of the pt spectrum
522           if (fPtSpectrum)
523           {
524             TParticle* mother = 0;
525
526             // preserve label for later
527             Int_t labelCopy = label;
528             if (labelCopy >= 0)
529               labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
530             if (labelCopy >= 0)
531               mother = stack->Particle(labelCopy);
532
533             // in case of pt study we do not count particles w/o label, because they cannot be scaled
534             if (!mother)
535               continue;
536
537             // it cannot be just multiplied because we cannot count "half of a particle"
538             // instead a random generator decides if the particle is counted twice (if value > 1) 
539             // or not (if value < 0)
540             Int_t bin = fPtSpectrum->FindBin(mother->Pt());
541             if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
542             {
543               Float_t factor = fPtSpectrum->GetBinContent(bin);
544               if (factor > 0)
545               {
546                 Float_t random = gRandom->Uniform();
547                 if (factor > 1 && random < factor - 1)
548                 {
549                   particleWeight = 2;
550                 }
551                 else if (factor < 1 && random < 1 - factor)
552                   particleWeight = 0;
553               }
554             }
555           }
556
557           //Printf("ESD weight is: %d", particleWeight);
558
559           if (TMath::Abs(eta) < 0.5)
560             nESDTracks05 += particleWeight;
561
562           if (TMath::Abs(eta) < 0.9)
563             nESDTracks09 += particleWeight;
564
565           if (TMath::Abs(eta) < 1.5)
566             nESDTracks15 += particleWeight;
567
568           if (TMath::Abs(eta) < 2.0)
569             nESDTracks20 += particleWeight;
570
571
572           if (fParticleCorrection[0] || fParticleSpecies)
573           {
574             Int_t motherLabel = -1;
575             TParticle* mother = 0;
576
577             // find mother
578             if (label >= 0)
579               motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
580             if (motherLabel >= 0)
581               mother = stack->Particle(motherLabel);
582
583             if (!mother)
584             {
585               // count tracks that did not have a label
586               if (TMath::Abs(eta) < etaRange)
587                 nESDTracksSpecies[4]++;
588               continue;
589             }
590
591             // get particle type (pion, proton, kaon, other)
592             Int_t id = -1;
593             switch (TMath::Abs(mother->GetPdgCode()))
594             {
595               case 211: id = 0; break;
596               case 321: id = 1; break;
597               case 2212: id = 2; break;
598               default: id = 3; break;
599             }
600
601             // double counting is ok for particle ratio study
602             if (TMath::Abs(eta) < etaRange)
603               nESDTracksSpecies[id]++;
604
605             // double counting is not ok for efficiency study
606
607             // 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)
608             if (foundTracks[label])
609             {
610               if (TMath::Abs(eta) < etaRange)
611                 nESDTracksSpecies[6]++;
612               continue;
613             }
614             foundTracks[label] = kTRUE;
615
616             // particle (primary) already counted?
617             if (foundPrimaries[motherLabel])
618             {
619               if (TMath::Abs(eta) < etaRange)
620                 nESDTracksSpecies[5]++;
621               continue;
622             }
623             foundPrimaries[motherLabel] = kTRUE;
624
625             if (fParticleCorrection[id])
626               fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
627           }
628         }
629
630         delete[] foundTracks;
631         delete[] foundPrimaries;
632
633         if ((Int_t) nMCTracks20 == 0 && nESDTracks20 > 3)
634             printf("WARNING: Event has %d generated and %d reconstructed...\n", nMCTracks20, nESDTracks20);
635
636         // fill response matrix using vtxMC (best guess)
637         fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks09,  nMCTracks15,  nMCTracks20,  nMCTracksAll,  nESDTracks05,  nESDTracks09,  nESDTracks15,  nESDTracks20);
638
639         fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
640
641         if (fParticleSpecies)
642           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]);
643       }
644     }
645   }
646
647   delete etaArr;
648   delete labelArr;
649 }
650
651 void AliMultiplicityTask::Terminate(Option_t *)
652 {
653   // The Terminate() function is the last function to be called during
654   // a query. It always runs on the client, it can be used to present
655   // the results graphically or save the results to file.
656
657   fOutput = dynamic_cast<TList*> (GetOutputData(0));
658   if (!fOutput) {
659     Printf("ERROR: fOutput not available");
660     return;
661   }
662
663   fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
664   for (Int_t i = 0; i < 4; ++i)
665     fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
666   fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
667
668   if (!fMultiplicity)
669   {
670     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
671     return;
672   }
673
674   TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
675
676   fMultiplicity->SaveHistograms();
677   for (Int_t i = 0; i < 4; ++i)
678     if (fParticleCorrection[i])
679       fParticleCorrection[i]->SaveHistograms();
680   if (fParticleSpecies)
681     fParticleSpecies->Write();
682
683   TObjString option(fOption);
684   option.Write();
685
686   file->Close();
687
688   Printf("Writting result to multiplicityMC.root");
689 }