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