refactoring multiplicity analysis
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / AliMultiplicityCorrection.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // This class is used to store correction maps, raw input and results of the multiplicity
19 // measurement with the ITS or TPC
20 // It also contains functions to correct the spectrum using different methods.
21 // e.g. chi2 minimization and bayesian unfolding
22 //
23 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
24
25 #include "AliMultiplicityCorrection.h"
26
27 #include <TFile.h>
28 #include <TH1F.h>
29 #include <TH2F.h>
30 #include <TH3F.h>
31 #include <TDirectory.h>
32 #include <TCanvas.h>
33 #include <TString.h>
34 #include <TF1.h>
35 #include <TMath.h>
36 #include <TCollection.h>
37 #include <TLegend.h>
38 #include <TLine.h>
39 #include <TRandom.h>
40 #include <TProfile.h>
41 #include <TProfile2D.h>
42
43 ClassImp(AliMultiplicityCorrection)
44
45 // These are the areas where the quality of the unfolding results are evaluated
46 // Default defined here, call SetQualityRegions to change them
47 // unit is in multiplicity (not in bin!)
48
49 // SPD:   peak area - flat area - low stat area
50 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {1,  20, 70};
51 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {10, 65, 80};
52
53 //____________________________________________________________________
54 void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
55 {
56   //
57   // sets the quality region definition to TPC or SPD
58   //
59
60   if (SPDStudy)
61   {
62     // SPD:   peak area - flat area - low stat area
63     fgQualityRegionsB[0] = 1;
64     fgQualityRegionsE[0] = 10;
65
66     fgQualityRegionsB[1] = 20;
67     fgQualityRegionsE[1] = 65;
68
69     fgQualityRegionsB[2] = 70;
70     fgQualityRegionsE[2] = 80;
71
72     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
73   }
74   else
75   {
76     // TPC:   peak area - flat area - low stat area
77     fgQualityRegionsB[0] = 4;
78     fgQualityRegionsE[0] = 12;
79
80     fgQualityRegionsB[1] = 25;
81     fgQualityRegionsE[1] = 55;
82
83     fgQualityRegionsB[2] = 88;
84     fgQualityRegionsE[2] = 108;
85
86     Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
87   }
88 }
89
90 //____________________________________________________________________
91 AliMultiplicityCorrection::AliMultiplicityCorrection() :
92   TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
93 {
94   //
95   // default constructor
96   //
97
98   for (Int_t i = 0; i < kESDHists; ++i)
99     fMultiplicityESD[i] = 0;
100
101   for (Int_t i = 0; i < kMCHists; ++i)
102   {
103     fMultiplicityVtx[i] = 0;
104     fMultiplicityMB[i] = 0;
105     fMultiplicityINEL[i] = 0;
106     fMultiplicityNSD[i] = 0;
107   }
108
109   for (Int_t i = 0; i < kCorrHists; ++i)
110   {
111     fCorrelation[i] = 0;
112     fMultiplicityESDCorrected[i] = 0;
113   }
114
115   for (Int_t i = 0; i < kQualityRegions; ++i)
116     fQuality[i] = 0;
117 }
118
119 //____________________________________________________________________
120 AliMultiplicityCorrection* AliMultiplicityCorrection::Open(const char* fileName, const char* folderName)
121 {
122   // opens the given file, reads the multiplicity from the given folder and returns the object
123   
124   TFile* file = TFile::Open(fileName);
125   if (!file)
126   {
127     Printf("ERROR: Could not open %s", fileName);
128     return 0;
129   }
130   
131   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folderName, folderName);
132   mult->LoadHistograms();
133   
134   // TODO closing the file does not work here, because the histograms cannot be read anymore. LoadHistograms need to be adapted
135   
136   return mult;
137 }
138
139 //____________________________________________________________________
140 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
141   TNamed(name, title),
142   fCurrentESD(0),
143   fCurrentCorrelation(0),
144   fCurrentEfficiency(0),
145   fLastBinLimit(0),
146   fLastChi2MC(0),
147   fLastChi2MCLimit(0),
148   fLastChi2Residuals(0),
149   fRatioAverage(0)
150 {
151   //
152   // named constructor
153   //
154
155   // do not add this hists to the directory
156   Bool_t oldStatus = TH1::AddDirectoryStatus();
157   TH1::AddDirectory(kFALSE);
158
159   /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
160   Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
161                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
162                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
163                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
164                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
165                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
166                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
167                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
168                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
169                           500.5 };
170                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
171                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
172                           //1000.5 };
173
174   #define VTXBINNING 10, binLimitsVtx
175   #define NBINNING fgkMaxParams, binLimitsN*/
176
177   #define NBINNING 201, -0.5, 200.5
178   #define VTXBINNING 1, -6, 6
179
180   for (Int_t i = 0; i < kESDHists; ++i)
181     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
182
183   for (Int_t i = 0; i < kMCHists; ++i)
184   {
185     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
186     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
187
188     fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
189     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
190
191     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
192     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
193     
194     fMultiplicityNSD[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityNSD%d", i)));
195     fMultiplicityNSD[i]->SetTitle("fMultiplicityNSD");
196   }
197
198   for (Int_t i = 0; i < kCorrHists; ++i)
199   {
200     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
201     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
202   }
203
204   TH1::AddDirectory(oldStatus);
205
206   AliUnfolding::SetNbins(120, 120);
207   AliUnfolding::SetSkipBinsBegin(1);
208   AliUnfolding::SetNormalizeInput(kTRUE);
209 }
210
211 //____________________________________________________________________
212 AliMultiplicityCorrection::~AliMultiplicityCorrection()
213 {
214   //
215   // Destructor
216   //
217
218   Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called");
219
220   for (Int_t i = 0; i < kESDHists; ++i)
221   {
222     if (fMultiplicityESD[i])
223       delete fMultiplicityESD[i];
224     fMultiplicityESD[i] = 0;
225   }
226
227   for (Int_t i = 0; i < kMCHists; ++i)
228   {
229     if (fMultiplicityVtx[i])
230       delete fMultiplicityVtx[i];
231     fMultiplicityVtx[i] = 0;
232
233     if (fMultiplicityMB[i])
234       delete fMultiplicityMB[i];
235     fMultiplicityMB[i] = 0;
236
237     if (fMultiplicityINEL[i])
238       delete fMultiplicityINEL[i];
239     fMultiplicityINEL[i] = 0;
240   
241     if (fMultiplicityNSD[i])
242       delete fMultiplicityNSD[i];
243     fMultiplicityNSD[i] = 0;
244 }
245
246   for (Int_t i = 0; i < kCorrHists; ++i)
247   {
248     if (fCorrelation[i])
249       delete fCorrelation[i];
250     fCorrelation[i] = 0;
251
252     if (fMultiplicityESDCorrected[i])
253       delete fMultiplicityESDCorrected[i];
254     fMultiplicityESDCorrected[i] = 0;
255   }
256 }
257
258 //____________________________________________________________________
259 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
260 {
261   // Merge a list of AliMultiplicityCorrection objects with this (needed for
262   // PROOF).
263   // Returns the number of merged objects (including this).
264
265   if (!list)
266     return 0;
267
268   if (list->IsEmpty())
269     return 1;
270
271   TIterator* iter = list->MakeIterator();
272   TObject* obj;
273
274   // collections of all histograms
275   TList collections[kESDHists+kMCHists*4+kCorrHists*2];
276
277   Int_t count = 0;
278   while ((obj = iter->Next())) {
279
280     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
281     if (entry == 0) 
282       continue;
283
284     for (Int_t i = 0; i < kESDHists; ++i)
285       collections[i].Add(entry->fMultiplicityESD[i]);
286
287     for (Int_t i = 0; i < kMCHists; ++i)
288     {
289       collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
290       collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
291       collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
292       collections[kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
293     }
294
295     for (Int_t i = 0; i < kCorrHists; ++i)
296       collections[kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
297
298     for (Int_t i = 0; i < kCorrHists; ++i)
299       collections[kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
300
301     count++;
302   }
303
304   for (Int_t i = 0; i < kESDHists; ++i)
305     fMultiplicityESD[i]->Merge(&collections[i]);
306
307   for (Int_t i = 0; i < kMCHists; ++i)
308   {
309     fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
310     fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
311     fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
312     fMultiplicityNSD[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
313   }
314
315   for (Int_t i = 0; i < kCorrHists; ++i)
316     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*4+i]);
317
318   for (Int_t i = 0; i < kCorrHists; ++i)
319     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*4+kCorrHists+i]);
320
321   delete iter;
322
323   return count+1;
324 }
325
326 //____________________________________________________________________
327 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
328 {
329   //
330   // loads the histograms from a file
331   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
332   //
333
334   if (!dir)
335     dir = GetName();
336
337   if (!gDirectory->cd(dir))
338     return kFALSE;
339
340   // store old hists to delete them later
341   TList oldObjects;
342   oldObjects.SetOwner(1);
343   for (Int_t i = 0; i < kESDHists; ++i)
344     if (fMultiplicityESD[i])
345       oldObjects.Add(fMultiplicityESD[i]);
346
347   for (Int_t i = 0; i < kMCHists; ++i)
348   {
349     if (fMultiplicityVtx[i])
350       oldObjects.Add(fMultiplicityVtx[i]);
351     if (fMultiplicityMB[i])
352       oldObjects.Add(fMultiplicityMB[i]);
353     if (fMultiplicityINEL[i])
354       oldObjects.Add(fMultiplicityINEL[i]);
355     if (fMultiplicityNSD[i])
356       oldObjects.Add(fMultiplicityNSD[i]);
357   }
358
359   for (Int_t i = 0; i < kCorrHists; ++i)
360     if (fCorrelation[i])
361       oldObjects.Add(fCorrelation[i]);
362
363   // load histograms
364
365   Bool_t success = kTRUE;
366
367   for (Int_t i = 0; i < kESDHists; ++i)
368   {
369     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
370     if (!fMultiplicityESD[i])
371       success = kFALSE;
372   }
373
374   for (Int_t i = 0; i < kMCHists; ++i)
375   {
376     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
377     fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
378     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
379     fMultiplicityNSD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityNSD[i]->GetName()));
380     if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
381       success = kFALSE;
382   }
383
384   for (Int_t i = 0; i < kCorrHists; ++i)
385   {
386     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
387     if (!fCorrelation[i])
388       success = kFALSE;
389     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
390     if (!fMultiplicityESDCorrected[i])
391       success = kFALSE;
392   }
393
394   gDirectory->cd("..");
395
396   // delete old hists
397   oldObjects.Delete();
398
399   return success;
400 }
401
402 //____________________________________________________________________
403 void AliMultiplicityCorrection::SaveHistograms(const char* dir)
404 {
405   //
406   // saves the histograms
407   //
408
409   if (!dir)
410     dir = GetName();
411
412   gDirectory->mkdir(dir);
413   gDirectory->cd(dir);
414
415   for (Int_t i = 0; i < kESDHists; ++i)
416     if (fMultiplicityESD[i])
417     {
418       fMultiplicityESD[i]->Write();
419       fMultiplicityESD[i]->ProjectionY(Form("%s_px", fMultiplicityESD[i]->GetName()), 1, fMultiplicityESD[i]->GetNbinsX())->Write();
420     }
421
422   for (Int_t i = 0; i < kMCHists; ++i)
423   {
424     if (fMultiplicityVtx[i])
425     {
426       fMultiplicityVtx[i]->Write();
427       fMultiplicityVtx[i]->ProjectionY(Form("%s_px", fMultiplicityVtx[i]->GetName()), 1, fMultiplicityVtx[i]->GetNbinsX())->Write();
428     }
429     if (fMultiplicityMB[i])
430     {
431       fMultiplicityMB[i]->Write();
432       fMultiplicityMB[i]->ProjectionY(Form("%s_px", fMultiplicityMB[i]->GetName()), 1, fMultiplicityMB[i]->GetNbinsX())->Write();
433     }
434     if (fMultiplicityINEL[i])
435     {
436       fMultiplicityINEL[i]->Write();
437       fMultiplicityINEL[i]->ProjectionY(Form("%s_px", fMultiplicityINEL[i]->GetName()), 1, fMultiplicityINEL[i]->GetNbinsX())->Write();
438     }
439     if (fMultiplicityNSD[i])
440     {
441       fMultiplicityNSD[i]->Write();
442       fMultiplicityNSD[i]->ProjectionY(Form("%s_px", fMultiplicityNSD[i]->GetName()), 1, fMultiplicityNSD[i]->GetNbinsX())->Write();
443     }
444   }
445
446   for (Int_t i = 0; i < kCorrHists; ++i)
447   {
448     if (fCorrelation[i])
449       fCorrelation[i]->Write();
450     if (fMultiplicityESDCorrected[i])
451       fMultiplicityESDCorrected[i]->Write();
452   }
453
454   gDirectory->cd("..");
455 }
456
457 //____________________________________________________________________
458 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, AliPWG0Helper::MCProcessType processType, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
459 {
460   //
461   // Fills an event from MC
462   //
463
464   if (triggered)
465   {
466     fMultiplicityMB[0]->Fill(vtx, generated05);
467     fMultiplicityMB[1]->Fill(vtx, generated10);
468     fMultiplicityMB[2]->Fill(vtx, generated15);
469     fMultiplicityMB[3]->Fill(vtx, generated20);
470     fMultiplicityMB[4]->Fill(vtx, generatedAll);
471
472     if (vertex)
473     {
474       fMultiplicityVtx[0]->Fill(vtx, generated05);
475       fMultiplicityVtx[1]->Fill(vtx, generated10);
476       fMultiplicityVtx[2]->Fill(vtx, generated15);
477       fMultiplicityVtx[3]->Fill(vtx, generated20);
478       fMultiplicityVtx[4]->Fill(vtx, generatedAll);
479     }
480   }
481
482   fMultiplicityINEL[0]->Fill(vtx, generated05);
483   fMultiplicityINEL[1]->Fill(vtx, generated10);
484   fMultiplicityINEL[2]->Fill(vtx, generated15);
485   fMultiplicityINEL[3]->Fill(vtx, generated20);
486   fMultiplicityINEL[4]->Fill(vtx, generatedAll);
487   
488   if (processType != AliPWG0Helper::kSD)
489   {
490     fMultiplicityNSD[0]->Fill(vtx, generated05);
491     fMultiplicityNSD[1]->Fill(vtx, generated10);
492     fMultiplicityNSD[2]->Fill(vtx, generated15);
493     fMultiplicityNSD[3]->Fill(vtx, generated20);
494     fMultiplicityNSD[4]->Fill(vtx, generatedAll);
495   }
496 }
497
498 //____________________________________________________________________
499 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
500 {
501   //
502   // Fills an event from ESD
503   //
504
505   fMultiplicityESD[0]->Fill(vtx, measured05);
506   fMultiplicityESD[1]->Fill(vtx, measured10);
507   fMultiplicityESD[2]->Fill(vtx, measured15);
508   fMultiplicityESD[3]->Fill(vtx, measured20);
509 }
510
511 //____________________________________________________________________
512 void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
513 {
514   //
515   // Fills an event into the correlation map with the information from MC and ESD
516   //
517
518   fCorrelation[0]->Fill(vtx, generated05, measured05);
519   fCorrelation[1]->Fill(vtx, generated10, measured10);
520   fCorrelation[2]->Fill(vtx, generated15, measured15);
521   fCorrelation[3]->Fill(vtx, generated20, measured20);
522
523   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
524   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
525   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
526   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
527 }
528
529 //____________________________________________________________________
530 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
531 {
532   //
533   // fills fCurrentESD, fCurrentCorrelation
534   // resets fMultiplicityESDCorrected
535   //
536
537   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
538
539   fMultiplicityESDCorrected[correlationID]->Reset();
540   fMultiplicityESDCorrected[correlationID]->Sumw2();
541
542   // project without under/overflow bins
543   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
544   fCurrentESD->Sumw2();
545
546   // empty under/overflow bins in x, otherwise Project3D takes them into account
547   TH3* hist = fCorrelation[correlationID];
548   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
549   {
550     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
551     {
552       hist->SetBinContent(0, y, z, 0);
553       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
554     }
555   }
556
557   fCurrentCorrelation = (TH2*) hist->Project3D("zy");
558   fCurrentCorrelation->Sumw2();
559   
560   Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
561
562 #if 0 // does not help
563   // null bins with one entry
564   Int_t nNulledBins = 0;
565   for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
566     for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
567     {
568       if (fCurrentCorrelation->GetBinContent(x, y) == 1)
569       {
570         fCurrentCorrelation->SetBinContent(x, y, 0);
571         fCurrentCorrelation->SetBinError(x, y, 0);
572
573         ++nNulledBins;
574       }
575     }
576   Printf("Nulled %d bins", nNulledBins);
577 #endif
578
579   fCurrentEfficiency = GetEfficiency(inputRange, eventType);
580   //fCurrentEfficiency->Rebin(2);
581   //fCurrentEfficiency->Scale(0.5);
582 }
583
584 //____________________________________________________________________
585 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
586 {
587   //
588   // calculates efficiency for given event type
589   //
590   
591   TH1* divisor = 0;
592   switch (eventType)
593   {
594     case kTrVtx : break;
595     case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
596     case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
597     case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY("divisor", 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
598   }
599   TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
600   
601   if (eventType == kTrVtx)
602   {
603     for (Int_t i=0; i<= eff->GetNbinsX()+1; i++)
604       eff->SetBinContent(i, 1);
605   }
606   else
607     eff->Divide(eff, divisor, 1, 1, "B");
608     
609   return eff;
610 }
611
612 //____________________________________________________________________
613 TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange)
614 {
615   //
616   // calculates efficiency for given event type
617   //
618
619   TH1* divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e");
620   TH1* eff = fMultiplicityMB[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
621   eff->Divide(eff, divisor, 1, 1, "B");
622   return eff;
623 }
624
625 //____________________________________________________________________
626 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
627 {
628   //
629   // correct spectrum using minuit chi2 method
630   //
631   // for description of parameters, see AliUnfolding::Unfold
632   //
633
634   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
635   
636   AliUnfolding::SetCreateOverflowBin(5);
637   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
638   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
639   
640   if (!initialConditions)
641   {
642     initialConditions = (TH1*) fCurrentESD->Clone("initialConditions");
643     initialConditions->Scale(1.0 / initialConditions->Integral());
644     if (!check)
645     {
646       // set minimum value to prevent MINUIT just staying in the small value
647       for (Int_t i=1; i<=initialConditions->GetNbinsX(); i++)
648         initialConditions->SetBinContent(i, TMath::Max(initialConditions->GetBinContent(i), 1e-3));
649     }
650   }
651
652   return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
653 }
654
655 //____________________________________________________________________
656 Int_t AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
657 {
658   //
659   // correct spectrum using minuit chi2 method with a NBD function
660   //
661   // for description of parameters, see AliUnfolding::Unfold
662   //
663
664   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
665   
666   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kFunction);
667   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
668   
669         TF1* func = new TF1("nbd", "[0] * TMath::Gamma([2]+x) / TMath::Gamma([2]) / TMath::Gamma(x+1) * pow([1] / ([1]+[2]), x) * pow(1.0 + [1]/[2], -[2])");
670         func->SetParNames("scaling", "averagen", "k");
671         func->SetParLimits(0, 0, 1000);
672         func->SetParLimits(1, 1, 50);
673         func->SetParLimits(2, 1, 10);
674         func->SetParameters(1, 10, 2);
675   AliUnfolding::SetFunction(func);
676   
677   return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, 0, fMultiplicityESDCorrected[correlationID]);
678 }
679
680 //____________________________________________________________________
681 void AliMultiplicityCorrection::DrawHistograms()
682 {
683   //
684   // draws the histograms of this class
685   //
686
687   printf("ESD:\n");
688
689   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
690   canvas1->Divide(3, 2);
691   for (Int_t i = 0; i < kESDHists; ++i)
692   {
693     canvas1->cd(i+1);
694     fMultiplicityESD[i]->DrawCopy("COLZ");
695     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
696   }
697
698   printf("Vtx:\n");
699
700   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
701   canvas2->Divide(3, 2);
702   for (Int_t i = 0; i < kMCHists; ++i)
703   {
704     canvas2->cd(i+1);
705     fMultiplicityVtx[i]->DrawCopy("COLZ");
706     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
707   }
708
709   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
710   canvas3->Divide(3, 3);
711   for (Int_t i = 0; i < kCorrHists; ++i)
712   {
713     canvas3->cd(i+1);
714     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
715     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
716     {
717       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
718       {
719         hist->SetBinContent(0, y, z, 0);
720         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
721       }
722     }
723     TH1* proj = hist->Project3D("zy");
724     proj->DrawCopy("COLZ");
725   }
726 }
727
728 //____________________________________________________________________
729 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple)
730 {
731   // draw comparison plots
732
733
734   //mcHist->Rebin(2);
735
736   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
737
738   TString tmpStr;
739   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
740
741   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
742   {
743     printf("ERROR. Unfolded histogram is empty\n");
744     return;
745   }
746
747   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
748   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", 1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins());
749   fCurrentESD->Sumw2();
750   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
751
752   // normalize unfolded result to 1
753   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
754
755   //fCurrentESD->Scale(mcHist->Integral(2, 200));
756
757   //new TCanvas;
758   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
759   ratio->Divide(mcHist);
760   ratio->Draw("HIST");
761   ratio->Fit("pol0", "W0", "", 20, 120);
762   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
763   delete ratio;
764   mcHist->Scale(scalingFactor);*/
765
766   // find bin with <= 50 entries. this is used as limit for the chi2 calculation
767   Int_t mcBinLimit = 0;
768   for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
769   {
770     if (mcHist->GetBinContent(i) > 50)
771     {
772       mcBinLimit = i;
773     }
774     else
775       break;
776   }
777   Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
778   
779   // scale to 1
780   mcHist->Sumw2();
781   if (mcHist->Integral() > 0)
782     mcHist->Scale(1.0 / mcHist->Integral());
783
784   // calculate residual
785
786   // for that we convolute the response matrix with the gathered result
787   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
788   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
789   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
790   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
791   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
792   if (convolutedProj->Integral() > 0)
793   {
794     convolutedProj->Scale(1.0 / convolutedProj->Integral());
795   }
796   else
797     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
798
799   TH1* residual = (TH1*) convolutedProj->Clone("residual");
800   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
801
802   residual->Add(fCurrentESD, -1);
803   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
804
805   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
806
807   // find bin limit
808   Int_t lastBin = 0;
809   for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
810   {
811     if (fCurrentESD->GetBinContent(i) <= 5)
812     {
813       lastBin = i;
814       break;
815     }
816   }
817   
818   // TODO fix errors
819   Float_t chi2 = 0;
820   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
821   {
822     if (fCurrentESD->GetBinError(i) > 0)
823     {
824       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
825       if (i > 1 && i <= lastBin)
826         chi2 += value * value;
827       Printf("%d --> %f (%f)", i, value * value, chi2);
828       residual->SetBinContent(i, value);
829       residualHist->Fill(value);
830     }
831     else
832     {
833       //printf("Residual bin %d set to 0\n", i);
834       residual->SetBinContent(i, 0);
835     }
836     convolutedProj->SetBinError(i, 0);
837     residual->SetBinError(i, 0);
838   }
839   fLastChi2Residuals = chi2;
840
841   new TCanvas; residualHist->DrawCopy();
842
843   //residualHist->Fit("gaus", "N");
844   //delete residualHist;
845
846   printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, lastBin);
847
848   TCanvas* canvas1 = 0;
849   if (simple)
850   {
851     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 600);
852     canvas1->Divide(2, 1);
853   }
854   else
855   {
856     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
857     canvas1->Divide(2, 3);
858   }
859
860   canvas1->cd(1);
861   canvas1->cd(1)->SetGridx();
862   canvas1->cd(1)->SetGridy();
863   canvas1->cd(1)->SetTopMargin(0.05);
864   canvas1->cd(1)->SetRightMargin(0.05);
865   canvas1->cd(1)->SetLeftMargin(0.12);
866   canvas1->cd(1)->SetBottomMargin(0.12);
867   TH1* proj = (TH1*) mcHist->Clone("proj");
868
869   // normalize without 0 bin
870   proj->Scale(1.0 / proj->Integral(2, proj->GetNbinsX()));
871   Printf("Normalized without 0 bin!");
872   proj->GetXaxis()->SetRangeUser(0, 200);
873   proj->GetYaxis()->SetTitleOffset(1.4);
874   //proj->SetLabelSize(0.05, "xy");
875   //proj->SetTitleSize(0.05, "xy");
876   proj->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
877   proj->SetStats(kFALSE);
878
879   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
880   fMultiplicityESDCorrected[esdCorrId]->SetMarkerColor(2);
881   //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(5);
882   // normalize without 0 bin
883   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(2, fMultiplicityESDCorrected[esdCorrId]->GetNbinsX()));
884   Printf("Normalized without 0 bin!");
885
886   if (proj->GetEntries() > 0) {
887     proj->DrawCopy("HIST");
888     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
889   }
890   else
891     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
892
893   gPad->SetLogy();
894
895   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
896   legend->AddEntry(proj, "True distribution");
897   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "Unfolded distribution");
898   legend->SetFillColor(0);
899   legend->SetTextSize(0.04);
900   legend->Draw();
901   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
902
903   canvas1->cd(2);
904   canvas1->cd(2)->SetGridx();
905   canvas1->cd(2)->SetGridy();
906   canvas1->cd(2)->SetTopMargin(0.05);
907   canvas1->cd(2)->SetRightMargin(0.05);
908   canvas1->cd(2)->SetLeftMargin(0.12);
909   canvas1->cd(2)->SetBottomMargin(0.12);
910
911   gPad->SetLogy();
912   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
913   //fCurrentESD->SetLineColor(2);
914   fCurrentESD->SetTitle(Form(";Measured multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
915   fCurrentESD->SetStats(kFALSE);
916   fCurrentESD->GetYaxis()->SetTitleOffset(1.4);
917   //fCurrentESD->SetLabelSize(0.05, "xy");
918   //fCurrentESD->SetTitleSize(0.05, "xy");
919   fCurrentESD->DrawCopy("HIST E");
920
921   convolutedProj->SetLineColor(2);
922   convolutedProj->SetMarkerColor(2);
923   convolutedProj->SetMarkerStyle(5);
924   //proj2->SetMarkerColor(2);
925   //proj2->SetMarkerStyle(5);
926   convolutedProj->DrawCopy("HIST SAME P");
927
928   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
929   legend->AddEntry(fCurrentESD, "Measured distribution");
930   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution", "P");
931   legend->SetFillColor(0);
932   legend->SetTextSize(0.04);
933   legend->Draw();
934
935   //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
936   //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
937
938   /*Float_t chi2 = 0;
939   Float_t chi = 0;
940   fLastChi2MCLimit = 0;
941   Int_t limit = 0;
942   for (Int_t i=2; i<=200; ++i)
943   {
944     if (proj->GetBinContent(i) != 0)
945     {
946       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
947       chi2 += value * value;
948       chi += TMath::Abs(value);
949
950       //printf("%d %f\n", i, chi);
951
952       if (chi2 < 0.2)
953         fLastChi2MCLimit = i;
954
955       if (chi < 3)
956         limit = i;
957
958     }
959   }*/
960
961   /*chi2 = 0;
962   Float_t chi = 0;
963   Int_t limit = 0;
964   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
965   {
966     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
967     {
968       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
969       if (value > 1e8)
970         value = 1e8; //prevent arithmetic exception
971       else if (value < -1e8)
972         value = -1e8;
973       if (i > 1)
974       {
975         chi2 += value * value;
976         chi += TMath::Abs(value);
977       }
978       diffMCUnfolded->SetBinContent(i, value);
979     }
980     else
981     {
982       //printf("diffMCUnfolded bin %d set to 0\n", i);
983       diffMCUnfolded->SetBinContent(i, 0);
984     }
985     if (chi2 < 1000)
986       fLastChi2MCLimit = i;
987     if (chi < 1000)
988       limit = i;
989     if (i == 150)
990       fLastChi2MC = chi2;
991   }
992
993   printf("limits %d %d\n", fLastChi2MCLimit, limit);
994   fLastChi2MCLimit = limit;
995
996   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
997
998   if (!simple)
999   {
1000     /*canvas1->cd(3);
1001
1002     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1003     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1004     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1005     diffMCUnfolded->DrawCopy("HIST");
1006
1007     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1008     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1009       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1010
1011     //new TCanvas; fluctuation->DrawCopy();
1012     delete fluctuation;*/
1013
1014     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1015     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1016     legend->AddEntry(fMultiplicityMC, "MC");
1017     legend->AddEntry(fMultiplicityESD, "ESD");
1018     legend->Draw();*/
1019
1020     canvas1->cd(4);
1021     residual->GetYaxis()->SetRangeUser(-5, 5);
1022     residual->GetXaxis()->SetRangeUser(0, 200);
1023     residual->DrawCopy();
1024
1025     canvas1->cd(5);
1026
1027     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1028     ratio->Divide(mcHist);
1029     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1030     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1031     ratio->GetXaxis()->SetRangeUser(0, 200);
1032     ratio->SetStats(kFALSE);
1033     ratio->Draw("HIST");
1034
1035     Double_t ratioChi2 = 0;
1036     fRatioAverage = 0;
1037     fLastChi2MCLimit = 0;
1038     for (Int_t i=2; i<=150; ++i)
1039     {
1040       Float_t value = ratio->GetBinContent(i) - 1;
1041       if (value > 1e8)
1042         value = 1e8; //prevent arithmetic exception
1043       else if (value < -1e8)
1044         value = -1e8;
1045
1046       ratioChi2 += value * value;
1047       fRatioAverage += TMath::Abs(value);
1048
1049       if (ratioChi2 < 0.1)
1050         fLastChi2MCLimit = i;
1051     }
1052     fRatioAverage /= 149;
1053
1054     printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1055     // TODO FAKE
1056     fLastChi2MC = ratioChi2;
1057
1058     // FFT of ratio
1059     canvas1->cd(6);
1060     const Int_t kFFT = 128;
1061     Double_t fftReal[kFFT];
1062     Double_t fftImag[kFFT];
1063
1064     for (Int_t i=0; i<kFFT; ++i)
1065     {
1066       fftReal[i] = ratio->GetBinContent(i+1+10);
1067       // test: ;-)
1068       //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1069       fftImag[i] = 0;
1070     }
1071
1072     FFT(-1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
1073
1074     TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1075     fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1076     TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1077     TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1078     fftResult->Reset();
1079     fftResult2->Reset();
1080     fftResult3->Reset();
1081     fftResult->GetYaxis()->UnZoom();
1082     fftResult2->GetYaxis()->UnZoom();
1083     fftResult3->GetYaxis()->UnZoom();
1084
1085     //Printf("AFTER FFT");
1086     for (Int_t i=0; i<kFFT; ++i)
1087     {
1088       //Printf("%d: %f", i, fftReal[i]);
1089       fftResult->SetBinContent(i+1, fftReal[i]);
1090       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1091       {
1092         Printf("Nulled %d", i);
1093         fftReal[i] = 0;
1094       }*/
1095     }
1096
1097     fftResult->SetLineColor(1);
1098     fftResult->DrawCopy();
1099
1100
1101     //Printf("IMAG");
1102     for (Int_t i=0; i<kFFT; ++i)
1103     {
1104       //Printf("%d: %f", i, fftImag[i]);
1105       fftResult2->SetBinContent(i+1, fftImag[i]);
1106
1107       fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1108     }
1109
1110     fftResult2->SetLineColor(2);
1111     fftResult2->DrawCopy("SAME");
1112
1113     fftResult3->SetLineColor(4);
1114     fftResult3->DrawCopy("SAME");
1115
1116     for (Int_t i=1; i<kFFT - 1; ++i)
1117     {
1118       if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1119       {
1120         fftReal[i-1] = 0;
1121         fftReal[i] = 0;
1122         fftReal[i+1] = 0;
1123         fftImag[i-1] = 0;
1124         fftImag[i] = 0;
1125         fftImag[i+1] = 0;
1126         //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1127         //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
1128         //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1129       }
1130     }
1131
1132     FFT(1, TMath::Nint(TMath::Log(kFFT) / TMath::Log(2)), fftReal, fftImag);
1133
1134     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1135     fftResult4->Reset();
1136
1137     //Printf("AFTER BACK-TRAFO");
1138     for (Int_t i=0; i<kFFT; ++i)
1139     {
1140       //Printf("%d: %f", i, fftReal[i]);
1141       fftResult4->SetBinContent(i+1+10, fftReal[i]);
1142     }
1143
1144     canvas1->cd(5);
1145     fftResult4->SetLineColor(4);
1146     fftResult4->DrawCopy("SAME");
1147
1148     // plot (MC - Unfolded) / error (MC)
1149     canvas1->cd(3);
1150
1151     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1152     diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1153
1154     Int_t ndfQual[kQualityRegions];
1155     for (Int_t region=0; region<kQualityRegions; ++region)
1156     {
1157       fQuality[region] = 0;
1158       ndfQual[region] = 0;
1159     }
1160
1161
1162     Double_t newChi2 = 0;
1163     Double_t newChi2Limit150 = 0;
1164     Int_t ndf = 0;
1165     for (Int_t i=1; i<=diffMCUnfolded2->GetNbinsX(); ++i)
1166     {
1167       Double_t value = 0;
1168       if (proj->GetBinError(i) > 0)
1169       {
1170         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1171         newChi2 += value * value;
1172         if (i > 1 && i <= mcBinLimit)
1173           newChi2Limit150 += value * value;
1174         ++ndf;
1175
1176         for (Int_t region=0; region<kQualityRegions; ++region)
1177           if (diffMCUnfolded2->GetXaxis()->GetBinCenter(i) >= fgQualityRegionsB[region] - 0.1 && diffMCUnfolded2->GetXaxis()->GetBinCenter(i) <= fgQualityRegionsE[region] + 0.1) // 0.1 to avoid e.g. 3.9999 < 4 problem
1178           {
1179             fQuality[region] += TMath::Abs(value);
1180             ++ndfQual[region];
1181           }
1182       }
1183
1184       diffMCUnfolded2->SetBinContent(i, value);
1185     }
1186
1187     // normalize region to the number of entries
1188     for (Int_t region=0; region<kQualityRegions; ++region)
1189     {
1190       if (ndfQual[region] > 0)
1191         fQuality[region] /= ndfQual[region];
1192       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1193     }
1194
1195     if (mcBinLimit > 1)
1196     {
1197       // TODO another FAKE
1198       fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1199       Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
1200     }
1201     else
1202       fLastChi2MC = -1;
1203
1204     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
1205
1206     diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1207     diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5);
1208     diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1209     diffMCUnfolded2->DrawCopy("HIST");
1210   }
1211
1212   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1213 }
1214
1215 //____________________________________________________________________
1216 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1217 {
1218 /*-------------------------------------------------------------------------
1219    This computes an in-place complex-to-complex FFT
1220    x and y are the real and imaginary arrays of 2^m points.
1221    dir =  1 gives forward transform
1222    dir = -1 gives reverse transform
1223
1224      Formula: forward
1225                   N-1
1226                   ---
1227               1   \          - j k 2 pi n / N
1228       X(n) = ---   >   x(k) e                    = forward transform
1229               N   /                                n=0..N-1
1230                   ---
1231                   k=0
1232
1233       Formula: reverse
1234                   N-1
1235                   ---
1236                   \          j k 2 pi n / N
1237       X(n) =       >   x(k) e                    = forward transform
1238                   /                                n=0..N-1
1239                   ---
1240                   k=0
1241 */
1242
1243    Long_t   nn, i, i1, j, k, i2, l, l1, l2;
1244    Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1245
1246    /* Calculate the number of points */
1247    nn = 1;
1248    for (i = 0; i < m; i++)
1249        nn *= 2;
1250
1251    /* Do the bit reversal */
1252    i2 = nn >> 1;
1253    j = 0;
1254    for (i= 0; i < nn - 1; i++) {
1255        if (i < j) {
1256            tx = x[i];
1257            ty = y[i];
1258            x[i] = x[j];
1259            y[i] = y[j];
1260            x[j] = tx;
1261            y[j] = ty;
1262        }
1263        k = i2;
1264        while (k <= j) {
1265            j -= k;
1266            k >>= 1;
1267        }
1268        j += k;
1269    }
1270
1271    /* Compute the FFT */
1272    c1 = -1.0;
1273    c2 = 0.0;
1274    l2 = 1;
1275    for (l = 0; l < m; l++) {
1276        l1 = l2;
1277        l2 <<= 1;
1278        u1 = 1.0;
1279        u2 = 0.0;
1280        for (j = 0;j < l1; j++) {
1281            for (i = j; i < nn; i += l2) {
1282                i1 = i + l1;
1283                t1 = u1 * x[i1] - u2 * y[i1];
1284                t2 = u1 * y[i1] + u2 * x[i1];
1285                x[i1] = x[i] - t1;
1286                y[i1] = y[i] - t2;
1287                x[i] += t1;
1288                y[i] += t2;
1289            }
1290            z =  u1 * c1 - u2 * c2;
1291            u2 = u1 * c2 + u2 * c1;
1292            u1 = z;
1293        }
1294        c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1295        if (dir == 1)
1296            c2 = -c2;
1297        c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1298    }
1299
1300    /* Scaling for forward transform */
1301    if (dir == 1) {
1302        for (i=0;i<nn;i++) {
1303            x[i] /= (Double_t)nn;
1304            y[i] /= (Double_t)nn;
1305        }
1306    }
1307 }
1308
1309 //____________________________________________________________________
1310 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
1311 {
1312   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1313   // These values are computed during DrawComparison, thus this function picks up the
1314   // last calculation
1315
1316   if (mc)
1317     *mc = fLastChi2MC;
1318   if (mcLimit)
1319     *mcLimit = fLastChi2MCLimit;
1320   if (residuals)
1321     *residuals = fLastChi2Residuals;
1322   if (ratioAverage)
1323     *ratioAverage = fRatioAverage;
1324 }
1325
1326 //____________________________________________________________________
1327 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1328 {
1329   //
1330   // returns the corresponding MC spectrum
1331   //
1332
1333   switch (eventType)
1334   {
1335     case kTrVtx : return fMultiplicityVtx[i]; break;
1336     case kMB : return fMultiplicityMB[i]; break;
1337     case kINEL : return fMultiplicityINEL[i]; break;
1338     case kNSD : return fMultiplicityNSD[i]; break;
1339   }
1340
1341   return 0;
1342 }
1343
1344 //____________________________________________________________________
1345 void AliMultiplicityCorrection::SetMultiplicityMC(Int_t i, EventType eventType, TH2F* hist)
1346 {
1347   //
1348   // returns the corresponding MC spectrum
1349   //
1350
1351   switch (eventType)
1352   {
1353     case kTrVtx : fMultiplicityVtx[i] = hist; break;
1354     case kMB : fMultiplicityMB[i] = hist; break;
1355     case kINEL : fMultiplicityINEL[i] = hist; break;
1356     case kNSD : fMultiplicityNSD[i] = hist; break;
1357   }
1358 }
1359
1360 //____________________________________________________________________
1361 TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
1362 {
1363   // calculate standard deviation of (results[0] - results[k]) k=1...max-1
1364   // per bin one gets: sigma(r[0] - r[n]) / r[0]
1365
1366   TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
1367   standardDeviation->Reset();
1368
1369   for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1370   {
1371     if (results[0]->GetBinContent(x) > 0)
1372     {
1373       Double_t average = 0;
1374       for (Int_t n=1; n<max; ++n)
1375         average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1376       average /= max-1;
1377
1378       Double_t variance = 0;
1379       for (Int_t n=1; n<max; ++n)
1380       {
1381         Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1382         variance += value * value;
1383       }
1384       variance /= max-1;
1385
1386       Double_t standardDev = TMath::Sqrt(variance);
1387       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1388       //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1389     }
1390   }
1391
1392   return standardDeviation;
1393 }
1394
1395 //____________________________________________________________________
1396 TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
1397 {
1398   //
1399   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1400   // the function unfolds the spectrum using the default response matrix and several modified ones
1401   // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1402   // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1403   // in <compareTo> (optional)
1404   //
1405   // returns the error assigned to the measurement
1406   //
1407
1408   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1409
1410   // initialize seed with current time
1411   gRandom->SetSeed(0);
1412   
1413   if (methodType == AliUnfolding::kChi2Minimization)
1414     AliUnfolding::SetCreateOverflowBin(5);
1415   AliUnfolding::SetUnfoldingMethod(methodType);
1416
1417   const Int_t kErrorIterations = 150;
1418
1419   TH1* maxError = 0;
1420   TH1* firstResult = 0;
1421
1422   TH1** results = new TH1*[kErrorIterations];
1423
1424   for (Int_t n=0; n<kErrorIterations; ++n)
1425   {
1426     Printf("Iteration %d of %d...", n, kErrorIterations);
1427
1428     SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1429
1430     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1431
1432     if (n > 0)
1433     {
1434       if (randomizeResponse)
1435       {
1436         // randomize response matrix
1437         for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1438           for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1439             fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1440       }
1441
1442       if (randomizeMeasured)
1443       {
1444         // randomize measured spectrum
1445         for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1446         {
1447           Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1448           measured->SetBinContent(x, randomValue);
1449           measured->SetBinError(x, TMath::Sqrt(randomValue));
1450         }
1451       }
1452     }
1453
1454     // only for bayesian method we have to do it before the call to Unfold...
1455     if (methodType == AliUnfolding::kBayesian)
1456     {
1457       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1458       {
1459         // with this it is normalized to 1
1460         Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1461
1462         // with this normalized to the given efficiency
1463         if (fCurrentEfficiency->GetBinContent(i) > 0)
1464           sum /= fCurrentEfficiency->GetBinContent(i);
1465         else
1466           sum = 0;
1467
1468         for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1469         {
1470           if (sum > 0)
1471           {
1472             fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1473             fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1474           }
1475           else
1476           {
1477             fCurrentCorrelation->SetBinContent(i, j, 0);
1478             fCurrentCorrelation->SetBinError(i, j, 0);
1479           }
1480         }
1481       }
1482     }
1483
1484     TH1* result = 0;
1485     if (n == 0 && compareTo)
1486     {
1487       // in this case we just store the histogram we want to compare to
1488       result = (TH1*) compareTo->Clone("compareTo");
1489       result->Sumw2();
1490     }
1491     else
1492     {
1493       result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
1494
1495       Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result);
1496
1497       if (returnCode != 0)
1498         return 0;
1499     }
1500
1501     // normalize
1502     result->Scale(1.0 / result->Integral());
1503
1504     if (n == 0)
1505     {
1506       firstResult = (TH1*) result->Clone("firstResult");
1507
1508       maxError = (TH1*) result->Clone("maxError");
1509       maxError->Reset();
1510     }
1511     else
1512     {
1513       // calculate ratio
1514       TH1* ratio = (TH1*) firstResult->Clone("ratio");
1515       ratio->Divide(result);
1516
1517       // find max. deviation
1518       for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
1519         maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
1520
1521       delete ratio;
1522     }
1523
1524     results[n] = result;
1525   }
1526
1527   // find covariance matrix
1528   // results[n] is X_x
1529   // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
1530
1531   Int_t nBins = results[0]->GetNbinsX();
1532   Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
1533   Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
1534
1535   // find average, E(X)
1536   TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
1537   for (Int_t n=1; n<kErrorIterations; ++n)
1538     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1539       average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
1540   //new TCanvas; average->DrawClone();
1541   
1542   // find cov. matrix
1543   TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
1544
1545   for (Int_t n=1; n<kErrorIterations; ++n)
1546     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1547       for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
1548       {
1549         // (X_x - E(X_x)) * (X_y - E(X_y)
1550         Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
1551         covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
1552       }
1553   TCanvas* c = new TCanvas; c->cd(); covMatrix->DrawCopy("COLZ");
1554
1555 //   // fill 2D histogram that contains deviation from first
1556 //   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
1557 //   for (Int_t n=1; n<kErrorIterations; ++n)
1558 //     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1559 //       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
1560 //   //new TCanvas; deviations->DrawCopy("COLZ");
1561 // 
1562 //   // get standard deviation "by hand"
1563 //   for (Int_t x=1; x<=nBins; x++)
1564 //   {
1565 //     if (results[0]->GetBinContent(x) > 0)
1566 //     {
1567 //       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
1568 //       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
1569 //       //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1570 //       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1571 //     }
1572 //   }
1573
1574   TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
1575
1576   // compare maxError to RMS of profile (variable name: average)
1577   // first: calculate rms in percent of value
1578   TH1* rmsError = (TH1*) maxError->Clone("rmsError");
1579   rmsError->Reset();
1580
1581   // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
1582   average->SetErrorOption("s");
1583   for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
1584     if (average->GetBinContent(x) > 0)
1585       rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
1586
1587   // find maxError deviation from average (not from "first result"/mc as above)
1588   TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
1589   maxError2->Reset();
1590   for (Int_t n=1; n<kErrorIterations; ++n)
1591     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1592       if (average->GetBinContent(x) > 0)
1593         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
1594
1595   //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
1596
1597   // plot difference between average and MC/first result
1598   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
1599   averageFirstRatio->Reset();
1600   averageFirstRatio->Divide(results[0], average);
1601
1602   //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
1603   //new TCanvas; averageFirstRatio->DrawCopy();
1604
1605   static TH1* temp = 0;
1606   if (!temp)
1607   {
1608     temp = (TH1*) standardDeviation->Clone("temp");
1609     for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1610       temp->SetBinContent(x, temp->GetBinContent(x) * results[0]->GetBinContent(x));
1611   }
1612   else
1613   {
1614     // find difference from result[0] as TH2
1615     TH2F* pulls = new TH2F("pulls", "pulls;multiplicity;difference", nBins, lowEdge, upEdge, 1000, -10, 10);
1616     for (Int_t n=1; n<kErrorIterations; ++n)
1617       for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1618         if (temp->GetBinContent(x) > 0)
1619           pulls->Fill(results[n]->GetXaxis()->GetBinCenter(x), (results[0]->GetBinContent(x) - results[n]->GetBinContent(x)) / temp->GetBinContent(x));
1620     new TCanvas("pulls", "pulls", 800, 600); pulls->DrawCopy(); pulls->FitSlicesY();
1621   }
1622
1623   // clean up
1624   for (Int_t n=0; n<kErrorIterations; ++n)
1625     delete results[n];
1626   delete[] results;
1627
1628   // fill into result histogram
1629   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1630     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
1631
1632   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1633     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
1634
1635   return standardDeviation;
1636 }
1637
1638 //____________________________________________________________________
1639 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
1640 {
1641   //
1642   // correct spectrum using bayesian method
1643   //
1644
1645   // initialize seed with current time
1646   gRandom->SetSeed(0);
1647
1648   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1649
1650   // normalize correction for given nPart
1651   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1652   {
1653     // with this it is normalized to 1
1654     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1655
1656     // with this normalized to the given efficiency
1657     if (fCurrentEfficiency->GetBinContent(i) > 0)
1658       sum /= fCurrentEfficiency->GetBinContent(i);
1659     else
1660       sum = 0;
1661
1662     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1663     {
1664       if (sum > 0)
1665       {
1666         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1667         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1668       }
1669       else
1670       {
1671         fCurrentCorrelation->SetBinContent(i, j, 0);
1672         fCurrentCorrelation->SetBinError(i, j, 0);
1673       }
1674     }
1675   }
1676
1677   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1678
1679   AliUnfolding::SetBayesianParameters(regPar, nIterations);
1680   AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
1681   if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
1682     return;
1683
1684   if (!determineError)
1685   {
1686     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
1687     return;
1688   }
1689
1690   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
1691   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
1692   // spectrum calculated. This is performed N times and the sigma is taken as the statistical
1693   // error of the unfolding method itself.
1694
1695   const Int_t kErrorIterations = 20;
1696
1697   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
1698
1699   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
1700   TH1* resultArray[kErrorIterations+1];
1701   for (Int_t n=0; n<kErrorIterations; ++n)
1702   {
1703     // randomize the content of clone following a poisson with the mean = the value of that bin
1704     for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
1705     {
1706       Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1707       //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
1708       randomized->SetBinContent(x, randomValue);
1709       randomized->SetBinError(x, TMath::Sqrt(randomValue));
1710     }
1711
1712     TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
1713     result2->Reset();
1714     if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0)
1715       return;
1716
1717     result2->Scale(1.0 / result2->Integral());
1718
1719     resultArray[n+1] = result2;
1720   }
1721   delete randomized;
1722
1723   resultArray[0] = fMultiplicityESDCorrected[correlationID];
1724   TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
1725
1726   for (Int_t n=0; n<kErrorIterations; ++n)
1727     delete resultArray[n+1];
1728
1729   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1730     fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
1731
1732   delete error;
1733 }
1734
1735 //____________________________________________________________________
1736 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
1737 {
1738   //
1739   // helper function for the covariance matrix of the bayesian method
1740   //
1741
1742   Float_t result = 0;
1743
1744   if (k == u && r == i)
1745     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
1746
1747   if (k == u)
1748     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
1749
1750   if (r == i)
1751     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
1752
1753   result *= matrixM[k][i];
1754
1755   return result;
1756 }
1757
1758 //____________________________________________________________________
1759 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
1760 {
1761   //
1762   // correct spectrum using bayesian method
1763   //
1764
1765   Float_t regPar = 0;
1766
1767   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1768   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1769
1770   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1771   //normalize ESD
1772   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1773
1774   // TODO should be taken from correlation map
1775   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1776
1777   // normalize correction for given nPart
1778   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1779   {
1780     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1781     //Double_t sum = sumHist->GetBinContent(i);
1782     if (sum <= 0)
1783       continue;
1784     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1785     {
1786       // npart sum to 1
1787       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1788       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1789     }
1790   }
1791
1792   new TCanvas;
1793   fCurrentCorrelation->Draw("COLZ");
1794
1795   // FAKE
1796   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1797
1798   // pick prior distribution
1799   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1800   Float_t norm = 1; //hPrior->Integral();
1801   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1802     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1803
1804   // zero distribution
1805   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
1806
1807   // define temp hist
1808   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1809   hTemp->Reset();
1810
1811   // just a shortcut
1812   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1813
1814   // unfold...
1815   Int_t iterations = 25;
1816   for (Int_t i=0; i<iterations; i++)
1817   {
1818     //printf(" iteration %i \n", i);
1819
1820     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1821     {
1822       Float_t value = 0;
1823       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1824         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
1825       hTemp->SetBinContent(m, value);
1826       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
1827     }
1828
1829     // regularization (simple smoothing)
1830     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1831
1832     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1833     {
1834       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1835                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1836                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1837       average *= hTrueSmoothed->GetBinWidth(t);
1838
1839       // weight the average with the regularization parameter
1840       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1841     }
1842
1843     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1844       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
1845
1846     // fill guess
1847     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1848     {
1849       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1850       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1851
1852       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1853     }
1854
1855
1856     // calculate chi2 (change from last iteration)
1857     Double_t chi2 = 0;
1858
1859     // use smoothed true (normalized) as new prior
1860     norm = 1; //hTrueSmoothed->Integral();
1861
1862     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1863     {
1864       Float_t newValue = hTemp->GetBinContent(t)/norm;
1865       Float_t diff = hPrior->GetBinContent(t) - newValue;
1866       chi2 += (Double_t) diff * diff;
1867
1868       hPrior->SetBinContent(t, newValue);
1869     }
1870
1871     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1872
1873     //if (i % 5 == 0)
1874       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1875
1876     delete hTrueSmoothed;
1877   } // end of iterations
1878
1879   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1880 }
1881
1882 //____________________________________________________________________
1883 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
1884 {
1885   //
1886   // correct spectrum using a simple Gaussian approach, that is model-dependent
1887   //
1888
1889   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1890   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1891
1892   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
1893   //normalize ESD
1894   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1895
1896   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
1897   correction->SetTitle("GaussianMean");
1898
1899   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
1900   correction->SetTitle("GaussianWidth");
1901
1902   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1903   {
1904     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
1905     proj->Fit("gaus", "0Q");
1906     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1907     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1908
1909     continue;
1910
1911     // draw for debugging
1912     new TCanvas;
1913     proj->DrawCopy();
1914     proj->GetFunction("gaus")->DrawCopy("SAME");
1915   }
1916
1917   TH1* target = fMultiplicityESDCorrected[correlationID];
1918
1919   Int_t nBins = target->GetNbinsX()*10+1;
1920   Float_t* binning = new Float_t[nBins];
1921   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1922     for (Int_t j=0; j<10; ++j)
1923       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
1924
1925   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
1926
1927   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
1928
1929   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1930   {
1931     Float_t mean = correction->GetBinContent(i);
1932     Float_t width = correctionWidth->GetBinContent(i);
1933
1934     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
1935     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
1936     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1937
1938     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1939     {
1940       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
1941     }
1942   }
1943
1944   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1945   {
1946     Float_t sum = 0;
1947     for (Int_t j=1; j<=10; ++j)
1948       sum += fineBinned->GetBinContent((i-1)*10 + j);
1949     target->SetBinContent(i, sum / 10);
1950   }
1951
1952   delete[] binning;
1953
1954   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
1955 }
1956
1957 //____________________________________________________________________
1958 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
1959 {
1960   // runs the distribution given in inputMC through the response matrix identified by
1961   // correlationMap and produces a measured distribution
1962   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
1963   // if normalized is set, inputMC is expected to be normalized to the bin width
1964
1965   if (!inputMC)
1966     return 0;
1967
1968   if (correlationMap < 0 || correlationMap >= kCorrHists)
1969     return 0;
1970
1971   // empty under/overflow bins in x, otherwise Project3D takes them into account
1972   TH3* hist = fCorrelation[correlationMap];
1973   for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
1974   {
1975     for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
1976     {
1977       hist->SetBinContent(0, y, z, 0);
1978       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1979     }
1980   }
1981
1982   TH2* corr = (TH2*) hist->Project3D("zy");
1983   //corr->Rebin2D(2, 1);
1984   corr->Sumw2();
1985   Printf("Correction histogram used for convolution has %f entries", corr->Integral());
1986
1987   // normalize correction for given nPart
1988   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
1989   {
1990     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
1991     if (sum <= 0)
1992       continue;
1993
1994     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1995     {
1996       // npart sum to 1
1997       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
1998       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
1999     }
2000   }
2001
2002   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2003   target->Reset();
2004
2005   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2006   {
2007     Float_t measured = 0;
2008     Float_t error = 0;
2009
2010     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2011     {
2012       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2013
2014       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2015       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2016     }
2017
2018     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2019
2020     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2021     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2022   }
2023
2024   return target;
2025 }
2026
2027 //____________________________________________________________________
2028 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2029 {
2030   // uses the given function to fill the input MC histogram and generates from that
2031   // the measured histogram by applying the response matrix
2032   // this can be used to evaluate if the methods work indepedently of the input
2033   // distribution
2034   // WARNING does not respect the vertex distribution, just fills central vertex bin
2035
2036   if (!inputMC)
2037     return;
2038
2039   if (id < 0 || id >= kESDHists)
2040     return;
2041
2042   // fill histogram used for random generation
2043   TH1* tmp = fMultiplicityVtx[id]->ProjectionY("tmp");
2044   tmp->Reset();
2045
2046   for (Int_t i=1; i<=tmp->GetNbinsX(); ++i)
2047     tmp->SetBinContent(i, inputMC->Eval(tmp->GetXaxis()->GetBinCenter(i)) * tmp->GetXaxis()->GetBinWidth(i));
2048     
2049   TH1* mcRnd = fMultiplicityVtx[id]->ProjectionY("mcRnd");
2050   mcRnd->Reset();
2051   mcRnd->FillRandom(tmp, tmp->Integral());
2052   
2053   //new TCanvas; tmp->Draw();
2054   //new TCanvas; mcRnd->Draw();
2055   
2056   // and move into 2d histogram
2057   TH1* mc = fMultiplicityVtx[id];
2058   mc->Reset();
2059   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2060   {
2061     mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, mcRnd->GetBinContent(i));
2062     mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, TMath::Sqrt(mcRnd->GetBinContent(i)));
2063   }
2064   
2065   //new TCanvas; mc->Draw("COLZ");
2066
2067   // now randomize the measured histogram; funcMeasured is used as pilot function to generated the measured entries
2068   TH1* funcMeasured = CalculateMultiplicityESD(tmp, id)->ProjectionY("funcMeasured");
2069   
2070   //new TCanvas; funcMeasured->Draw();
2071   
2072   fMultiplicityESD[id]->Reset();
2073   
2074   TH1* measRnd = fMultiplicityESD[id]->ProjectionY("measRnd");
2075   measRnd->FillRandom(funcMeasured, tmp->Integral());
2076   
2077   //new TCanvas; measRnd->Draw();
2078   
2079   fMultiplicityESD[id]->Reset();
2080   for (Int_t i=1; i<=fMultiplicityESD[id]->GetNbinsY(); ++i)
2081   {
2082     fMultiplicityESD[id]->SetBinContent(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, measRnd->GetBinContent(i));
2083     fMultiplicityESD[id]->SetBinError(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, TMath::Sqrt(measRnd->GetBinContent(i)));
2084   }
2085 }