1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
23 // Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
25 #include "AliMultiplicityCorrection.h"
31 #include <TDirectory.h>
32 #include <TVirtualFitter.h>
37 #include <TCollection.h>
42 #include <TProfile2D.h>
44 ClassImp(AliMultiplicityCorrection)
46 const Int_t AliMultiplicityCorrection::fgkMaxInput = 250; // bins in measured histogram
47 const Int_t AliMultiplicityCorrection::fgkMaxParams = 250; // bins in unfolded histogram = number of fit params
49 TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0;
50 TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0;
51 TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0;
52 TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0;
54 Bool_t AliMultiplicityCorrection::fgCreateBigBin = kTRUE;
55 AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1;
56 Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000;
57 Int_t AliMultiplicityCorrection::fgNParamsMinuit = AliMultiplicityCorrection::fgkMaxParams;
59 TF1* AliMultiplicityCorrection::fgNBD = 0;
61 Float_t AliMultiplicityCorrection::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
62 Int_t AliMultiplicityCorrection::fgBayesianIterations = 100; // number of iterations in Bayesian method
64 // These are the areas where the quality of the unfolding results are evaluated
65 // Default defined here, call SetQualityRegions to change them
66 // unit is in multiplicity (not in bin!)
68 // SPD: peak area - flat area - low stat area
69 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4, 60, 180};
70 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210};
72 //____________________________________________________________________
73 void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
76 // sets the quality region definition to TPC or SPD
81 // SPD: peak area - flat area - low stat area
82 fgQualityRegionsB[0] = 4;
83 fgQualityRegionsE[0] = 20;
85 fgQualityRegionsB[1] = 60;
86 fgQualityRegionsE[1] = 140;
88 fgQualityRegionsB[2] = 180;
89 fgQualityRegionsE[2] = 210;
91 Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
95 // TPC: peak area - flat area - low stat area
96 fgQualityRegionsB[0] = 4;
97 fgQualityRegionsE[0] = 12;
99 fgQualityRegionsB[1] = 25;
100 fgQualityRegionsE[1] = 55;
102 fgQualityRegionsB[2] = 88;
103 fgQualityRegionsE[2] = 108;
105 Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
109 //____________________________________________________________________
110 AliMultiplicityCorrection::AliMultiplicityCorrection() :
111 TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
114 // default constructor
117 for (Int_t i = 0; i < kESDHists; ++i)
118 fMultiplicityESD[i] = 0;
120 for (Int_t i = 0; i < kMCHists; ++i)
122 fMultiplicityVtx[i] = 0;
123 fMultiplicityMB[i] = 0;
124 fMultiplicityINEL[i] = 0;
127 for (Int_t i = 0; i < kCorrHists; ++i)
130 fMultiplicityESDCorrected[i] = 0;
133 for (Int_t i = 0; i < kQualityRegions; ++i)
137 //____________________________________________________________________
138 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
139 TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
145 // do not add this hists to the directory
146 Bool_t oldStatus = TH1::AddDirectoryStatus();
147 TH1::AddDirectory(kFALSE);
149 /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
150 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,
151 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
152 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
153 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
154 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
155 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
156 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
157 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
158 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
160 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
161 //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
164 #define VTXBINNING 10, binLimitsVtx
165 #define NBINNING fgkMaxParams, binLimitsN*/
167 #define NBINNING 501, -0.5, 500.5
168 #define VTXBINNING 1, -10, 10
170 for (Int_t i = 0; i < kESDHists; ++i)
171 fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
173 for (Int_t i = 0; i < kMCHists; ++i)
175 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
176 fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
178 fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
179 fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
181 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
182 fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
185 for (Int_t i = 0; i < kCorrHists; ++i)
187 fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
188 fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
191 TH1::AddDirectory(oldStatus);
194 //____________________________________________________________________
195 AliMultiplicityCorrection::~AliMultiplicityCorrection()
201 Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called");
203 for (Int_t i = 0; i < kESDHists; ++i)
205 if (fMultiplicityESD[i])
206 delete fMultiplicityESD[i];
207 fMultiplicityESD[i] = 0;
210 for (Int_t i = 0; i < kMCHists; ++i)
212 if (fMultiplicityVtx[i])
213 delete fMultiplicityVtx[i];
214 fMultiplicityVtx[i] = 0;
216 if (fMultiplicityMB[i])
217 delete fMultiplicityMB[i];
218 fMultiplicityMB[i] = 0;
220 if (fMultiplicityINEL[i])
221 delete fMultiplicityINEL[i];
222 fMultiplicityINEL[i] = 0;
225 for (Int_t i = 0; i < kCorrHists; ++i)
228 delete fCorrelation[i];
231 if (fMultiplicityESDCorrected[i])
232 delete fMultiplicityESDCorrected[i];
233 fMultiplicityESDCorrected[i] = 0;
237 //____________________________________________________________________
238 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
240 // Merge a list of AliMultiplicityCorrection objects with this (needed for
242 // Returns the number of merged objects (including this).
250 TIterator* iter = list->MakeIterator();
253 // collections of all histograms
254 TList collections[kESDHists+kMCHists*3+kCorrHists*2];
257 while ((obj = iter->Next())) {
259 AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
263 for (Int_t i = 0; i < kESDHists; ++i)
264 collections[i].Add(entry->fMultiplicityESD[i]);
266 for (Int_t i = 0; i < kMCHists; ++i)
268 collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
269 collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
270 collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
273 for (Int_t i = 0; i < kCorrHists; ++i)
274 collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
276 for (Int_t i = 0; i < kCorrHists; ++i)
277 collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
282 for (Int_t i = 0; i < kESDHists; ++i)
283 fMultiplicityESD[i]->Merge(&collections[i]);
285 for (Int_t i = 0; i < kMCHists; ++i)
287 fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
288 fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
289 fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
292 for (Int_t i = 0; i < kCorrHists; ++i)
293 fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
295 for (Int_t i = 0; i < kCorrHists; ++i)
296 fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
303 //____________________________________________________________________
304 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
307 // loads the histograms from a file
308 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
314 if (!gDirectory->cd(dir))
317 // TODO memory leak. old histograms needs to be deleted.
319 Bool_t success = kTRUE;
321 for (Int_t i = 0; i < kESDHists; ++i)
323 fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
324 if (!fMultiplicityESD[i])
328 for (Int_t i = 0; i < kMCHists; ++i)
330 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
331 fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
332 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
333 if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
337 for (Int_t i = 0; i < kCorrHists; ++i)
339 fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
340 if (!fCorrelation[i])
342 fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
343 if (!fMultiplicityESDCorrected[i])
347 gDirectory->cd("..");
352 //____________________________________________________________________
353 void AliMultiplicityCorrection::SaveHistograms()
356 // saves the histograms
359 gDirectory->mkdir(GetName());
360 gDirectory->cd(GetName());
362 for (Int_t i = 0; i < kESDHists; ++i)
363 if (fMultiplicityESD[i])
364 fMultiplicityESD[i]->Write();
366 for (Int_t i = 0; i < kMCHists; ++i)
368 if (fMultiplicityVtx[i])
369 fMultiplicityVtx[i]->Write();
370 if (fMultiplicityMB[i])
371 fMultiplicityMB[i]->Write();
372 if (fMultiplicityINEL[i])
373 fMultiplicityINEL[i]->Write();
376 for (Int_t i = 0; i < kCorrHists; ++i)
379 fCorrelation[i]->Write();
380 if (fMultiplicityESDCorrected[i])
381 fMultiplicityESDCorrected[i]->Write();
384 gDirectory->cd("..");
387 //____________________________________________________________________
388 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
391 // Fills an event from MC
396 fMultiplicityMB[0]->Fill(vtx, generated05);
397 fMultiplicityMB[1]->Fill(vtx, generated10);
398 fMultiplicityMB[2]->Fill(vtx, generated15);
399 fMultiplicityMB[3]->Fill(vtx, generated20);
400 fMultiplicityMB[4]->Fill(vtx, generatedAll);
404 fMultiplicityVtx[0]->Fill(vtx, generated05);
405 fMultiplicityVtx[1]->Fill(vtx, generated10);
406 fMultiplicityVtx[2]->Fill(vtx, generated15);
407 fMultiplicityVtx[3]->Fill(vtx, generated20);
408 fMultiplicityVtx[4]->Fill(vtx, generatedAll);
412 fMultiplicityINEL[0]->Fill(vtx, generated05);
413 fMultiplicityINEL[1]->Fill(vtx, generated10);
414 fMultiplicityINEL[2]->Fill(vtx, generated15);
415 fMultiplicityINEL[3]->Fill(vtx, generated20);
416 fMultiplicityINEL[4]->Fill(vtx, generatedAll);
419 //____________________________________________________________________
420 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
423 // Fills an event from ESD
426 fMultiplicityESD[0]->Fill(vtx, measured05);
427 fMultiplicityESD[1]->Fill(vtx, measured10);
428 fMultiplicityESD[2]->Fill(vtx, measured15);
429 fMultiplicityESD[3]->Fill(vtx, measured20);
432 //____________________________________________________________________
433 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)
436 // Fills an event into the correlation map with the information from MC and ESD
439 fCorrelation[0]->Fill(vtx, generated05, measured05);
440 fCorrelation[1]->Fill(vtx, generated10, measured10);
441 fCorrelation[2]->Fill(vtx, generated15, measured15);
442 fCorrelation[3]->Fill(vtx, generated20, measured20);
444 fCorrelation[4]->Fill(vtx, generatedAll, measured05);
445 fCorrelation[5]->Fill(vtx, generatedAll, measured10);
446 fCorrelation[6]->Fill(vtx, generatedAll, measured15);
447 fCorrelation[7]->Fill(vtx, generatedAll, measured20);
450 //____________________________________________________________________
451 Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
453 // homogenity term for minuit fitting
454 // pure function of the parameters
455 // prefers constant function (pol0)
459 // ignore the first bin here. on purpose...
460 for (Int_t i=2; i<fgNParamsMinuit; ++i)
462 Double_t right = params[i];
463 Double_t left = params[i-1];
467 Double_t diff = 1 - left / right;
475 //____________________________________________________________________
476 Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
478 // homogenity term for minuit fitting
479 // pure function of the parameters
480 // prefers linear function (pol1)
484 // ignore the first bin here. on purpose...
485 for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
487 if (params[i-1] == 0)
490 Double_t right = params[i];
491 Double_t middle = params[i-1];
492 Double_t left = params[i-2];
494 Double_t der1 = (right - middle);
495 Double_t der2 = (middle - left);
497 Double_t diff = (der1 - der2) / middle;
505 //____________________________________________________________________
506 Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
508 // homogenity term for minuit fitting
509 // pure function of the parameters
510 // prefers linear function (pol1)
514 /*Float_t der[fgkMaxParams];
516 for (Int_t i=0; i<fgkMaxParams-1; ++i)
517 der[i] = params[i+1] - params[i];
519 for (Int_t i=0; i<fgkMaxParams-6; ++i)
521 for (Int_t j=1; j<=5; ++j)
523 Double_t diff = der[i+j] - der[i];
528 // ignore the first bin here. on purpose...
529 for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
531 if (params[i-1] == 0)
534 Double_t right = log(params[i]);
535 Double_t middle = log(params[i-1]);
536 Double_t left = log(params[i-2]);
538 Double_t der1 = (right - middle);
539 Double_t der2 = (middle - left);
541 Double_t diff = (der1 - der2) / middle;
549 //____________________________________________________________________
550 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
552 // homogenity term for minuit fitting
553 // pure function of the parameters
554 // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
555 // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
559 // ignore the first bin here. on purpose...
560 for (Int_t i=3; i<fgNParamsMinuit; ++i)
562 Double_t right = params[i];
563 Double_t middle = params[i-1];
564 Double_t left = params[i-2];
566 Double_t der1 = (right - middle);
567 Double_t der2 = (middle - left);
569 Double_t diff = (der1 - der2);
577 //____________________________________________________________________
578 Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
580 // homogenity term for minuit fitting
581 // pure function of the parameters
582 // calculates entropy, from
583 // The method of reduced cross-entropy (M. Schmelling 1993)
585 Double_t paramSum = 0;
586 // ignore the first bin here. on purpose...
587 for (Int_t i=1; i<fgNParamsMinuit; ++i)
588 paramSum += params[i];
591 for (Int_t i=1; i<fgNParamsMinuit; ++i)
593 //Double_t tmp = params[i] / paramSum;
594 Double_t tmp = params[i];
595 if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
596 chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
602 //____________________________________________________________________
603 void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
606 // fit function for minuit
608 // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
609 // func->SetParNames("scaling", "averagen", "k");
610 // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
611 // func->SetParLimits(1, 0.001, 1000);
612 // func->SetParLimits(2, 0.001, 1000);
615 fgNBD->SetParameters(params[0], params[1], params[2]);
617 Double_t params2[fgkMaxParams];
619 for (Int_t i=0; i<fgkMaxParams; ++i)
620 params2[i] = fgNBD->Eval(i);
622 MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
624 printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
627 //____________________________________________________________________
628 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
631 // fit function for minuit
632 // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
636 static TVectorD paramsVector(fgNParamsMinuit);
637 for (Int_t i=0; i<fgNParamsMinuit; ++i)
638 paramsVector[i] = params[i] * params[i];
640 // calculate penalty factor
641 Double_t penaltyVal = 0;
642 switch (fgRegularizationType)
645 case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
646 case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
647 case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
648 case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
649 case kLog: penaltyVal = RegularizationLog(paramsVector); break;
652 //if (penaltyVal > 1e10)
653 // paramsVector2.Print();
655 penaltyVal *= fgRegularizationWeight;
658 TVectorD measGuessVector(fgkMaxInput);
659 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
662 measGuessVector -= (*fgCurrentESDVector);
664 //measGuessVector.Print();
666 TVectorD copy(measGuessVector);
669 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
670 // normal way is like this:
671 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
672 // optimized way like this:
673 for (Int_t i=0; i<fgkMaxInput; ++i)
674 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
676 //measGuessVector.Print();
678 // reduce influence of first bin
681 // (Ad - m) W (Ad - m)
682 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
683 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
684 Double_t chi2FromFit = measGuessVector * copy * 1e6;
686 chi2 = chi2FromFit + penaltyVal;
688 static Int_t callCount = 0;
689 if ((callCount++ % 10000) == 0)
690 printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
693 //____________________________________________________________________
694 void AliMultiplicityCorrection::DrawGuess(Double_t *params)
697 // draws residuals of solution suggested by params and effect of regularization
701 static TVectorD paramsVector(fgNParamsMinuit);
702 for (Int_t i=0; i<fgNParamsMinuit; ++i)
703 paramsVector[i] = params[i] * params[i];
706 TVectorD measGuessVector(fgkMaxInput);
707 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
710 measGuessVector -= (*fgCurrentESDVector);
712 TVectorD copy(measGuessVector);
715 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
716 // normal way is like this:
717 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
718 // optimized way like this:
719 for (Int_t i=0; i<fgkMaxInput; ++i)
720 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
722 // (Ad - m) W (Ad - m)
723 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
724 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
725 //Double_t chi2FromFit = measGuessVector * copy * 1e6;
728 TH1* residuals = new TH1F("residuals", "residuals", fgkMaxInput, -0.5, fgkMaxInput - 0.5);
729 for (Int_t i=0; i<fgkMaxInput; ++i)
730 residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
731 new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
734 if (fgRegularizationType != kPol1) {
735 Printf("Drawing not supported");
739 TH1* penalty = new TH1F("penalty", "penalty", fgNParamsMinuit, -0.5, fgNParamsMinuit - 0.5);
741 for (Int_t i=2; i<fgNParamsMinuit; ++i)
743 if (params[i-1] == 0)
746 Double_t right = params[i];
747 Double_t middle = params[i-1];
748 Double_t left = params[i-2];
750 Double_t der1 = (right - middle);
751 Double_t der2 = (middle - left);
753 Double_t diff = (der1 - der2) / middle;
755 penalty->SetBinContent(i, diff * diff);
757 new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
760 //____________________________________________________________________
761 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
764 // fills fCurrentESD, fCurrentCorrelation
765 // resets fMultiplicityESDCorrected
768 Bool_t display = kFALSE;
770 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
772 fMultiplicityESDCorrected[correlationID]->Reset();
773 fMultiplicityESDCorrected[correlationID]->Sumw2();
775 fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD"); //"fCurrentESD", -1, -1, "e");
776 fCurrentESD->Sumw2();
778 // empty under/overflow bins in x, otherwise Project3D takes them into account
779 TH3* hist = fCorrelation[correlationID];
780 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
782 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
784 hist->SetBinContent(0, y, z, 0);
785 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
789 fCurrentCorrelation = hist->Project3D("zy");
790 //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
791 //fMultiplicityESDCorrected[correlationID]->Rebin(2);
792 fCurrentCorrelation->Sumw2();
794 Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
799 for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
801 if (fCurrentESD->GetBinContent(i) <= 5)
808 if (fLastBinLimit > 0)
814 canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
815 canvas->Divide(2, 2);
818 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
819 fCurrentESD->SetStats(kFALSE);
820 fCurrentESD->SetTitle(";measured multiplicity");
821 fCurrentESD->DrawCopy();
825 fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
826 fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
827 fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
828 fCurrentCorrelation->SetStats(kFALSE);
829 fCurrentCorrelation->DrawCopy("COLZ");
832 printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit);
833 fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX()));
834 for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i)
836 fCurrentESD->SetBinContent(i, 0);
837 fCurrentESD->SetBinError(i, 0);
839 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
840 fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit)));
842 printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit));
844 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
846 fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY()));
847 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
848 fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit)));
850 for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
852 fCurrentCorrelation->SetBinContent(i, j, 0);
853 fCurrentCorrelation->SetBinError(i, j, 0);
857 printf("Adjusted correlation matrix!\n");
862 fCurrentESD->DrawCopy();
866 fCurrentCorrelation->DrawCopy("COLZ");
871 #if 0 // does not help
872 // null bins with one entry
873 Int_t nNulledBins = 0;
874 for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
875 for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
877 if (fCurrentCorrelation->GetBinContent(x, y) == 1)
879 fCurrentCorrelation->SetBinContent(x, y, 0);
880 fCurrentCorrelation->SetBinError(x, y, 0);
885 Printf("Nulled %d bins", nNulledBins);
888 fCurrentEfficiency = GetEfficiency(inputRange, eventType);
889 //fCurrentEfficiency->Rebin(2);
890 //fCurrentEfficiency->Scale(0.5);
893 //____________________________________________________________________
894 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
897 // calculates efficiency for given event type
903 case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
904 case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
905 case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
907 TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
908 eff->Divide(eff, divisor, 1, 1, "B");
912 //____________________________________________________________________
913 void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams)
916 // sets the parameters for chi2 minimization
919 fgRegularizationType = type;
920 fgRegularizationWeight = weight;
922 printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
924 if (minuitParams > 0)
926 fgNParamsMinuit = minuitParams;
927 printf("AliMultiplicityCorrection::SetRegularizationParameters --> Number of MINUIT iterations set to %d\n", minuitParams);
931 //____________________________________________________________________
932 void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
935 // sets the parameters for Bayesian unfolding
938 fgBayesianSmoothing = smoothing;
939 fgBayesianIterations = nIterations;
941 printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
944 //____________________________________________________________________
945 Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
948 // implemenation of unfolding (internal function)
950 // unfolds <measured> using response from <correlation> and effiency <aEfficiency>
951 // output is in <result>
952 // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
953 // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
955 // returns minuit status (0 = success), (-1 when check was set)
959 measured->Scale(1.0 / measured->Integral());
961 if (!fgCorrelationMatrix)
962 fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgNParamsMinuit);
963 if (!fgCorrelationCovarianceMatrix)
964 fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
965 if (!fgCurrentESDVector)
966 fgCurrentESDVector = new TVectorD(fgkMaxInput);
967 if (!fgEntropyAPriori)
968 fgEntropyAPriori = new TVectorD(fgNParamsMinuit);
970 fgCorrelationMatrix->Zero();
971 fgCorrelationCovarianceMatrix->Zero();
972 fgCurrentESDVector->Zero();
973 fgEntropyAPriori->Zero();
975 // normalize correction for given nPart
976 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
978 Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
981 Float_t maxValue = 0;
983 for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
985 // find most probably value
986 if (maxValue < correlation->GetBinContent(i, j))
988 maxValue = correlation->GetBinContent(i, j);
993 correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
994 correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
996 if (i <= fgNParamsMinuit && j <= fgkMaxInput)
997 (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
1000 //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
1003 // Initialize TMinuit via generic fitter interface
1004 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgNParamsMinuit);
1005 Double_t arglist[100];
1007 // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
1009 minuit->ExecuteCommand("SET PRINT", arglist, 1);
1011 // however, enable warnings
1012 //minuit->ExecuteCommand("SET WAR", arglist, 0);
1014 // set minimization function
1015 minuit->SetFCN(MinuitFitFunction);
1017 for (Int_t i=0; i<fgNParamsMinuit; i++)
1018 (*fgEntropyAPriori)[i] = 1;
1020 if (!initialConditions) {
1021 initialConditions = measured;
1023 printf("Using different starting conditions...\n");
1024 //new TCanvas; initialConditions->DrawCopy();
1025 initialConditions->Scale(1.0 / initialConditions->Integral());
1028 for (Int_t i=0; i<fgNParamsMinuit; i++)
1029 if (initialConditions->GetBinContent(i+1) > 0)
1030 (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
1032 Double_t* results = new Double_t[fgNParamsMinuit+1];
1033 for (Int_t i=0; i<fgNParamsMinuit; ++i)
1035 results[i] = initialConditions->GetBinContent(i+1);
1039 results[i] = TMath::Max(results[i], 1e-3);
1041 Float_t stepSize = 0.1;
1043 // minuit sees squared values to prevent it from going negative...
1044 results[i] = TMath::Sqrt(results[i]);
1045 //stepSize /= results[i]; // keep relative step size
1047 minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
1049 //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0);
1050 // bin 0 is filled with value from bin 1 (otherwise it's 0)
1051 //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
1053 //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
1055 for (Int_t i=0; i<fgkMaxInput; ++i)
1057 (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
1058 if (measured->GetBinError(i+1) > 0)
1059 (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
1061 if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
1062 (*fgCorrelationCovarianceMatrix)(i, i) = 0;
1067 MinuitFitFunction(dummy, 0, chi2, results, 0);
1069 printf("Chi2 of initial parameters is = %f\n", chi2);
1077 // first param is number of iterations, second is precision....
1079 //arglist[1] = 1e-5;
1080 //minuit->ExecuteCommand("SCAN", arglist, 0);
1081 Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
1082 printf("MINUIT status is %d\n", status);
1083 //minuit->ExecuteCommand("MIGRAD", arglist, 1);
1084 //minuit->ExecuteCommand("MIGRAD", arglist, 1);
1085 //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
1086 //minuit->ExecuteCommand("MINOS", arglist, 0);
1088 for (Int_t i=0; i<fgNParamsMinuit; ++i)
1090 if (aEfficiency->GetBinContent(i+1) > 0)
1092 results[i] = minuit->GetParameter(i);
1093 result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1));
1094 // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
1095 result->SetBinError(i+1, minuit->GetParError(i) * results[i] / aEfficiency->GetBinContent(i+1));
1099 result->SetBinContent(i+1, 0);
1100 result->SetBinError(i+1, 0);
1110 //____________________________________________________________________
1111 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
1114 // correct spectrum using minuit chi2 method
1116 // for description of parameters, see UnfoldWithMinuit
1119 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1120 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, fgCreateBigBin);
1122 return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
1125 //____________________________________________________________________
1126 void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
1129 // correct spectrum using minuit chi2 method applying a NBD fit
1132 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1134 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
1136 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1138 fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1139 fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1140 fgCurrentESDVector = new TVectorD(fgkMaxParams);
1142 // normalize correction for given nPart
1143 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1145 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1148 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1151 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1152 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1154 if (i <= fgkMaxParams && j <= fgkMaxParams)
1155 (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
1159 for (Int_t i=0; i<fgkMaxParams; ++i)
1161 (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
1162 if (fCurrentESD->GetBinError(i+1) > 0)
1163 (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
1166 // Create NBD function
1169 fgNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250);
1170 fgNBD->SetParNames("scaling", "averagen", "k");
1173 // Initialize TMinuit via generic fitter interface
1174 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
1176 minuit->SetFCN(MinuitNBD);
1177 minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
1178 minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
1179 minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
1181 Double_t arglist[100];
1183 minuit->ExecuteCommand("SET PRINT", arglist, 1);
1184 minuit->ExecuteCommand("MIGRAD", arglist, 0);
1185 //minuit->ExecuteCommand("MINOS", arglist, 0);
1187 fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
1189 for (Int_t i=0; i<fgkMaxParams; ++i)
1191 printf("%d : %f\n", i, fgNBD->Eval(i));
1192 if (fgNBD->Eval(i) > 0)
1194 fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
1195 fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
1200 //____________________________________________________________________
1201 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
1204 // normalizes a 1-d histogram to its bin width
1207 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1209 hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
1210 hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
1214 //____________________________________________________________________
1215 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
1218 // normalizes a 2-d histogram to its bin width (x width * y width)
1221 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1222 for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
1224 Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
1225 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
1226 hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
1230 //____________________________________________________________________
1231 void AliMultiplicityCorrection::DrawHistograms()
1234 // draws the histograms of this class
1239 TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
1240 canvas1->Divide(3, 2);
1241 for (Int_t i = 0; i < kESDHists; ++i)
1244 fMultiplicityESD[i]->DrawCopy("COLZ");
1245 printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1250 TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1251 canvas2->Divide(3, 2);
1252 for (Int_t i = 0; i < kMCHists; ++i)
1255 fMultiplicityVtx[i]->DrawCopy("COLZ");
1256 printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
1259 TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1260 canvas3->Divide(3, 3);
1261 for (Int_t i = 0; i < kCorrHists; ++i)
1264 TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
1265 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1267 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1269 hist->SetBinContent(0, y, z, 0);
1270 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1273 TH1* proj = hist->Project3D("zy");
1274 proj->DrawCopy("COLZ");
1278 //____________________________________________________________________
1279 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
1283 Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1286 tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1288 if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1290 printf("ERROR. Unfolded histogram is empty\n");
1294 //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1295 fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
1296 fCurrentESD->Sumw2();
1297 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1299 // normalize unfolded result to 1
1300 fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1302 //fCurrentESD->Scale(mcHist->Integral(2, 200));
1305 /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1306 ratio->Divide(mcHist);
1307 ratio->Draw("HIST");
1308 ratio->Fit("pol0", "W0", "", 20, 120);
1309 Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1311 mcHist->Scale(scalingFactor);*/
1313 // find bin with <= 50 entries. this is used as limit for the chi2 calculation
1314 Int_t mcBinLimit = 0;
1315 for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
1317 if (mcHist->GetBinContent(i) > 50)
1324 Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
1328 mcHist->Scale(1.0 / mcHist->Integral());
1330 // calculate residual
1332 // for that we convolute the response matrix with the gathered result
1333 // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
1334 TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1335 tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
1336 TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1337 TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1338 if (convolutedProj->Integral() > 0)
1340 convolutedProj->Scale(1.0 / convolutedProj->Integral());
1343 printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1345 TH1* residual = (TH1*) convolutedProj->Clone("residual");
1346 residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
1348 residual->Add(fCurrentESD, -1);
1349 //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1351 TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
1355 for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
1357 if (fCurrentESD->GetBinError(i) > 0)
1359 Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1360 if (i > 1 && i <= fLastBinLimit)
1361 chi2 += value * value;
1362 Printf("%d --> %f (%f)", i, value * value, chi2);
1363 residual->SetBinContent(i, value);
1364 residualHist->Fill(value);
1368 //printf("Residual bin %d set to 0\n", i);
1369 residual->SetBinContent(i, 0);
1371 convolutedProj->SetBinError(i, 0);
1372 residual->SetBinError(i, 0);
1374 fLastChi2Residuals = chi2;
1376 new TCanvas; residualHist->DrawCopy();
1378 //residualHist->Fit("gaus", "N");
1379 //delete residualHist;
1381 printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit);
1383 TCanvas* canvas1 = 0;
1386 canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1387 canvas1->Divide(2, 1);
1391 canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1392 canvas1->Divide(2, 3);
1396 TH1* proj = (TH1*) mcHist->Clone("proj");
1398 proj->GetXaxis()->SetRangeUser(0, 200);
1399 proj->SetTitle(";true multiplicity;Entries");
1400 proj->SetStats(kFALSE);
1402 fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1404 if (proj->GetEntries() > 0) {
1405 proj->DrawCopy("HIST");
1406 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
1409 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
1413 TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1414 legend->AddEntry(proj, "true distribution");
1415 legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1416 legend->SetFillColor(0);
1418 // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1423 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1424 //fCurrentESD->SetLineColor(2);
1425 fCurrentESD->SetTitle(";measured multiplicity;Entries");
1426 fCurrentESD->SetStats(kFALSE);
1427 fCurrentESD->DrawCopy("HIST E");
1429 convolutedProj->SetLineColor(2);
1430 //proj2->SetMarkerColor(2);
1431 //proj2->SetMarkerStyle(5);
1432 convolutedProj->DrawCopy("HIST SAME");
1434 legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1435 legend->AddEntry(fCurrentESD, "measured distribution");
1436 legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1437 legend->SetFillColor(0);
1440 //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1441 //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1445 fLastChi2MCLimit = 0;
1447 for (Int_t i=2; i<=200; ++i)
1449 if (proj->GetBinContent(i) != 0)
1451 Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
1452 chi2 += value * value;
1453 chi += TMath::Abs(value);
1455 //printf("%d %f\n", i, chi);
1458 fLastChi2MCLimit = i;
1469 for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
1471 if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
1473 Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1475 value = 1e8; //prevent arithmetic exception
1476 else if (value < -1e8)
1480 chi2 += value * value;
1481 chi += TMath::Abs(value);
1483 diffMCUnfolded->SetBinContent(i, value);
1487 //printf("diffMCUnfolded bin %d set to 0\n", i);
1488 diffMCUnfolded->SetBinContent(i, 0);
1491 fLastChi2MCLimit = i;
1498 printf("limits %d %d\n", fLastChi2MCLimit, limit);
1499 fLastChi2MCLimit = limit;
1501 printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
1507 diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1508 //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1509 diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1510 diffMCUnfolded->DrawCopy("HIST");
1512 TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1513 for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1514 fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1516 //new TCanvas; fluctuation->DrawCopy();
1517 delete fluctuation;*/
1519 /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1520 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1521 legend->AddEntry(fMultiplicityMC, "MC");
1522 legend->AddEntry(fMultiplicityESD, "ESD");
1526 //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
1527 residual->GetXaxis()->SetRangeUser(0, 200);
1528 residual->DrawCopy();
1532 TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1533 ratio->Divide(mcHist);
1534 ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1535 ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1536 ratio->GetXaxis()->SetRangeUser(0, 200);
1537 ratio->SetStats(kFALSE);
1538 ratio->Draw("HIST");
1540 Double_t ratioChi2 = 0;
1542 fLastChi2MCLimit = 0;
1543 for (Int_t i=2; i<=150; ++i)
1545 Float_t value = ratio->GetBinContent(i) - 1;
1547 value = 1e8; //prevent arithmetic exception
1548 else if (value < -1e8)
1551 ratioChi2 += value * value;
1552 fRatioAverage += TMath::Abs(value);
1554 if (ratioChi2 < 0.1)
1555 fLastChi2MCLimit = i;
1557 fRatioAverage /= 149;
1559 printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1561 fLastChi2MC = ratioChi2;
1565 const Int_t kFFT = 128;
1566 Double_t fftReal[kFFT];
1567 Double_t fftImag[kFFT];
1569 for (Int_t i=0; i<kFFT; ++i)
1571 fftReal[i] = ratio->GetBinContent(i+1+10);
1573 //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1577 FFT(-1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1579 TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1580 fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1581 TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1582 TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1584 fftResult2->Reset();
1585 fftResult3->Reset();
1586 fftResult->GetYaxis()->UnZoom();
1587 fftResult2->GetYaxis()->UnZoom();
1588 fftResult3->GetYaxis()->UnZoom();
1590 //Printf("AFTER FFT");
1591 for (Int_t i=0; i<kFFT; ++i)
1593 //Printf("%d: %f", i, fftReal[i]);
1594 fftResult->SetBinContent(i+1, fftReal[i]);
1595 /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1597 Printf("Nulled %d", i);
1602 fftResult->SetLineColor(1);
1603 fftResult->DrawCopy();
1607 for (Int_t i=0; i<kFFT; ++i)
1609 //Printf("%d: %f", i, fftImag[i]);
1610 fftResult2->SetBinContent(i+1, fftImag[i]);
1612 fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1615 fftResult2->SetLineColor(2);
1616 fftResult2->DrawCopy("SAME");
1618 fftResult3->SetLineColor(4);
1619 fftResult3->DrawCopy("SAME");
1621 for (Int_t i=1; i<kFFT - 1; ++i)
1623 if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1631 //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1632 //fftImag[i] = (fftImag[i-1] + fftImag[i+1]) / 2;
1633 //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1637 FFT(1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1639 TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1640 fftResult4->Reset();
1642 //Printf("AFTER BACK-TRAFO");
1643 for (Int_t i=0; i<kFFT; ++i)
1645 //Printf("%d: %f", i, fftReal[i]);
1646 fftResult4->SetBinContent(i+1+10, fftReal[i]);
1650 fftResult4->SetLineColor(4);
1651 fftResult4->DrawCopy("SAME");
1653 // plot (MC - Unfolded) / error (MC)
1656 TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1657 diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1659 Int_t ndfQual[kQualityRegions];
1660 for (Int_t region=0; region<kQualityRegions; ++region)
1662 fQuality[region] = 0;
1663 ndfQual[region] = 0;
1667 Double_t newChi2 = 0;
1668 Double_t newChi2Limit150 = 0;
1670 for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
1673 if (proj->GetBinError(i) > 0)
1675 value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1676 newChi2 += value * value;
1677 if (i > 1 && i <= mcBinLimit)
1678 newChi2Limit150 += value * value;
1681 for (Int_t region=0; region<kQualityRegions; ++region)
1682 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
1684 fQuality[region] += TMath::Abs(value);
1689 diffMCUnfolded2->SetBinContent(i, value);
1692 // normalize region to the number of entries
1693 for (Int_t region=0; region<kQualityRegions; ++region)
1695 if (ndfQual[region] > 0)
1696 fQuality[region] /= ndfQual[region];
1697 Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1702 // TODO another FAKE
1703 fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1704 Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
1709 Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, newChi2 / ndf);
1711 diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1712 //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1713 diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1714 diffMCUnfolded2->DrawCopy("HIST");
1717 canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1720 //____________________________________________________________________
1721 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1723 /*-------------------------------------------------------------------------
1724 This computes an in-place complex-to-complex FFT
1725 x and y are the real and imaginary arrays of 2^m points.
1726 dir = 1 gives forward transform
1727 dir = -1 gives reverse transform
1732 1 \ - j k 2 pi n / N
1733 X(n) = --- > x(k) e = forward transform
1742 X(n) = > x(k) e = forward transform
1748 Long_t nn, i, i1, j, k, i2, l, l1, l2;
1749 Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1751 /* Calculate the number of points */
1753 for (i = 0; i < m; i++)
1756 /* Do the bit reversal */
1759 for (i= 0; i < nn - 1; i++) {
1776 /* Compute the FFT */
1780 for (l = 0; l < m; l++) {
1785 for (j = 0;j < l1; j++) {
1786 for (i = j; i < nn; i += l2) {
1788 t1 = u1 * x[i1] - u2 * y[i1];
1789 t2 = u1 * y[i1] + u2 * x[i1];
1795 z = u1 * c1 - u2 * c2;
1796 u2 = u1 * c2 + u2 * c1;
1799 c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1802 c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1805 /* Scaling for forward transform */
1807 for (i=0;i<nn;i++) {
1808 x[i] /= (Double_t)nn;
1809 y[i] /= (Double_t)nn;
1814 //____________________________________________________________________
1815 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
1817 // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1818 // These values are computed during DrawComparison, thus this function picks up the
1824 *mcLimit = fLastChi2MCLimit;
1826 *residuals = fLastChi2Residuals;
1828 *ratioAverage = fRatioAverage;
1831 //____________________________________________________________________
1832 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1835 // returns the corresponding MC spectrum
1840 case kTrVtx : return fMultiplicityVtx[i]; break;
1841 case kMB : return fMultiplicityMB[i]; break;
1842 case kINEL : return fMultiplicityINEL[i]; break;
1848 //____________________________________________________________________
1849 TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
1851 // calculate standard deviation of (results[0] - results[k]) k=1...max-1
1852 // per bin one gets: sigma(r[0] - r[n]) / r[0]
1854 TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
1855 standardDeviation->Reset();
1857 for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1859 if (results[0]->GetBinContent(x) > 0)
1861 Double_t average = 0;
1862 for (Int_t n=1; n<max; ++n)
1863 average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1866 Double_t variance = 0;
1867 for (Int_t n=1; n<max; ++n)
1869 Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1870 variance += value * value;
1874 Double_t standardDev = TMath::Sqrt(variance);
1875 standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1876 //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1880 return standardDeviation;
1883 //____________________________________________________________________
1884 TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
1887 // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1888 // the function unfolds the spectrum using the default response matrix and several modified ones
1889 // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1890 // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1891 // in <compareTo> (optional)
1893 // returns the error assigned to the measurement
1896 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1898 // initialize seed with current time
1899 gRandom->SetSeed(0);
1901 const Int_t kErrorIterations = 150;
1904 TH1* firstResult = 0;
1906 TH1** results = new TH1*[kErrorIterations];
1908 for (Int_t n=0; n<kErrorIterations; ++n)
1910 Printf("Iteration %d of %d...", n, kErrorIterations);
1912 // only done for chi2 minimization
1913 Bool_t createBigBin = kFALSE;
1914 if (methodType == kChi2Minimization)
1915 createBigBin = kTRUE;
1917 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
1919 TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1923 if (randomizeResponse)
1925 // randomize response matrix
1926 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1927 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1928 fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1931 if (randomizeMeasured)
1933 // randomize measured spectrum
1934 for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1936 Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1937 measured->SetBinContent(x, randomValue);
1938 measured->SetBinError(x, TMath::Sqrt(randomValue));
1943 // only for bayesian method we have to do it before the call to Unfold...
1944 if (methodType == kBayesian)
1946 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1948 // with this it is normalized to 1
1949 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1951 // with this normalized to the given efficiency
1952 if (fCurrentEfficiency->GetBinContent(i) > 0)
1953 sum /= fCurrentEfficiency->GetBinContent(i);
1957 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1961 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1962 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1966 fCurrentCorrelation->SetBinContent(i, j, 0);
1967 fCurrentCorrelation->SetBinError(i, j, 0);
1974 if (n == 0 && compareTo)
1976 // in this case we just store the histogram we want to compare to
1977 result = (TH1*) compareTo->Clone("compareTo");
1982 result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
1984 Int_t returnCode = -1;
1985 if (methodType == kChi2Minimization)
1987 returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
1989 else if (methodType == kBayesian)
1990 returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
1992 if (returnCode != 0)
1997 result->Scale(1.0 / result->Integral());
2001 firstResult = (TH1*) result->Clone("firstResult");
2003 maxError = (TH1*) result->Clone("maxError");
2009 TH1* ratio = (TH1*) firstResult->Clone("ratio");
2010 ratio->Divide(result);
2012 // find max. deviation
2013 for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
2014 maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
2019 results[n] = result;
2022 // find covariance matrix
2023 // results[n] is X_x
2024 // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
2026 Int_t nBins = results[0]->GetNbinsX();
2027 Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
2028 Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
2030 // find average, E(X)
2031 TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
2032 for (Int_t n=1; n<kErrorIterations; ++n)
2033 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2034 average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
2035 //new TCanvas; average->DrawClone();
2038 TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
2040 for (Int_t n=1; n<kErrorIterations; ++n)
2041 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2042 for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
2044 // (X_x - E(X_x)) * (X_y - E(X_y)
2045 Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
2046 covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
2048 new TCanvas; covMatrix->DrawClone("COLZ");
2050 // // fill 2D histogram that contains deviation from first
2051 // TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
2052 // for (Int_t n=1; n<kErrorIterations; ++n)
2053 // for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2054 // deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
2055 // //new TCanvas; deviations->DrawCopy("COLZ");
2057 // // get standard deviation "by hand"
2058 // for (Int_t x=1; x<=nBins; x++)
2060 // if (results[0]->GetBinContent(x) > 0)
2062 // TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
2063 // Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
2064 // //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
2065 // Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
2069 TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
2071 // compare maxError to RMS of profile (variable name: average)
2072 // first: calculate rms in percent of value
2073 TH1* rmsError = (TH1*) maxError->Clone("rmsError");
2076 // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
2077 average->SetErrorOption("s");
2078 for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
2079 if (average->GetBinContent(x) > 0)
2080 rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
2082 // find maxError deviation from average (not from "first result"/mc as above)
2083 TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
2085 for (Int_t n=1; n<kErrorIterations; ++n)
2086 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2087 if (average->GetBinContent(x) > 0)
2088 maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
2090 //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
2092 // plot difference between average and MC/first result
2093 TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
2094 averageFirstRatio->Reset();
2095 averageFirstRatio->Divide(results[0], average);
2097 //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
2098 //new TCanvas; averageFirstRatio->DrawCopy();
2101 for (Int_t n=0; n<kErrorIterations; ++n)
2105 // fill into result histogram
2106 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2107 fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
2109 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2110 fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2112 return standardDeviation;
2115 //____________________________________________________________________
2116 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
2119 // correct spectrum using bayesian method
2122 // initialize seed with current time
2123 gRandom->SetSeed(0);
2125 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2127 // normalize correction for given nPart
2128 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2130 // with this it is normalized to 1
2131 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2133 // with this normalized to the given efficiency
2134 if (fCurrentEfficiency->GetBinContent(i) > 0)
2135 sum /= fCurrentEfficiency->GetBinContent(i);
2139 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2143 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2144 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2148 fCurrentCorrelation->SetBinContent(i, j, 0);
2149 fCurrentCorrelation->SetBinError(i, j, 0);
2154 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2156 if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
2159 if (!determineError)
2161 Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
2165 // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
2166 // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
2167 // spectrum calculated. This is performed N times and the sigma is taken as the statistical
2168 // error of the unfolding method itself.
2170 const Int_t kErrorIterations = 20;
2172 printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
2174 TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
2175 TH1* resultArray[kErrorIterations+1];
2176 for (Int_t n=0; n<kErrorIterations; ++n)
2178 // randomize the content of clone following a poisson with the mean = the value of that bin
2179 for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
2181 Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
2182 //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
2183 randomized->SetBinContent(x, randomValue);
2184 randomized->SetBinError(x, TMath::Sqrt(randomValue));
2187 TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
2189 if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
2192 result2->Scale(1.0 / result2->Integral());
2194 resultArray[n+1] = result2;
2198 resultArray[0] = fMultiplicityESDCorrected[correlationID];
2199 TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
2201 for (Int_t n=0; n<kErrorIterations; ++n)
2202 delete resultArray[n+1];
2204 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2205 fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2210 //____________________________________________________________________
2211 Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
2214 // unfolds a spectrum
2217 if (measured->Integral() <= 0)
2219 Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
2223 const Int_t kStartBin = 0;
2225 const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
2226 const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
2228 // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
2229 const Double_t kConvergenceLimit = kMaxT * 1e-6;
2231 // store information in arrays, to increase processing speed (~ factor 5)
2232 Double_t measuredCopy[kMaxM];
2233 Double_t measuredError[kMaxM];
2234 Double_t prior[kMaxT];
2235 Double_t result[kMaxT];
2236 Double_t efficiency[kMaxT];
2238 Double_t response[kMaxT][kMaxM];
2239 Double_t inverseResponse[kMaxT][kMaxM];
2241 // for normalization
2242 Float_t measuredIntegral = measured->Integral();
2243 for (Int_t m=0; m<kMaxM; m++)
2245 measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
2246 measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
2248 for (Int_t t=0; t<kMaxT; t++)
2250 response[t][m] = correlation->GetBinContent(t+1, m+1);
2251 inverseResponse[t][m] = 0;
2255 for (Int_t t=0; t<kMaxT; t++)
2257 efficiency[t] = aEfficiency->GetBinContent(t+1);
2258 prior[t] = measuredCopy[t];
2262 // pick prior distribution
2263 if (initialConditions)
2265 printf("Using different starting conditions...\n");
2266 // for normalization
2267 Float_t inputDistIntegral = initialConditions->Integral();
2268 for (Int_t i=0; i<kMaxT; i++)
2269 prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
2272 //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
2275 for (Int_t i=0; i<nIterations || nIterations < 0; i++)
2277 //printf(" iteration %i \n", i);
2279 // calculate IR from Beyes theorem
2280 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
2282 Double_t chi2Measured = 0;
2283 for (Int_t m=0; m<kMaxM; m++)
2286 for (Int_t t = kStartBin; t<kMaxT; t++)
2287 norm += response[t][m] * prior[t];
2289 // calc. chi2: (measured - response * prior) / error
2290 if (measuredError[m] > 0)
2292 Double_t value = (measuredCopy[m] - norm) / measuredError[m];
2293 chi2Measured += value * value;
2298 for (Int_t t = kStartBin; t<kMaxT; t++)
2299 inverseResponse[t][m] = response[t][m] * prior[t] / norm;
2303 for (Int_t t = kStartBin; t<kMaxT; t++)
2304 inverseResponse[t][m] = 0;
2307 //Printf("chi2Measured of the last prior is %e", chi2Measured);
2309 for (Int_t t = kStartBin; t<kMaxT; t++)
2312 for (Int_t m=0; m<kMaxM; m++)
2313 value += inverseResponse[t][m] * measuredCopy[m];
2315 if (efficiency[t] > 0)
2316 result[t] = value / efficiency[t];
2321 Double_t chi2LastIter = 0;
2322 // regularization (simple smoothing)
2323 for (Int_t t=kStartBin; t<kMaxT; t++)
2325 Float_t newValue = 0;
2326 // 0 bin excluded from smoothing
2327 if (t > kStartBin+1 && t<kMaxT-1)
2329 Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
2331 // weight the average with the regularization parameter
2332 newValue = (1 - regPar) * result[t] + regPar * average;
2335 newValue = result[t];
2337 // calculate chi2 (change from last iteration)
2338 if (prior[t] > 1e-5)
2340 Double_t diff = (prior[t] - newValue) / prior[t];
2341 chi2LastIter += diff * diff;
2344 prior[t] = newValue;
2346 //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
2347 //convergence->Fill(i+1, chi2LastIter);
2349 if (nIterations < 0 && chi2LastIter < kConvergenceLimit)
2351 Printf("Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
2356 // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
2357 } // end of iterations
2359 //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
2360 //delete convergence;
2362 for (Int_t t=0; t<kMaxT; t++)
2363 aResult->SetBinContent(t+1, result[t]);
2368 // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
2370 /*printf("Calculating covariance matrix. This may take some time...\n");
2372 // TODO check if this is the right one...
2373 TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2375 Int_t xBins = hInverseResponseBayes->GetNbinsX();
2376 Int_t yBins = hInverseResponseBayes->GetNbinsY();
2378 // calculate "unfolding matrix" Mij
2379 Float_t matrixM[251][251];
2380 for (Int_t i=1; i<=xBins; i++)
2382 for (Int_t j=1; j<=yBins; j++)
2384 if (fCurrentEfficiency->GetBinContent(i) > 0)
2385 matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
2387 matrixM[i-1][j-1] = 0;
2391 Float_t* vectorn = new Float_t[yBins];
2392 for (Int_t j=1; j<=yBins; j++)
2393 vectorn[j-1] = fCurrentESD->GetBinContent(j);
2395 // first part of covariance matrix, depends on input distribution n(E)
2396 Float_t cov1[251][251];
2398 Float_t nEvents = fCurrentESD->Integral(); // N
2403 for (Int_t k=0; k<xBins; k++)
2405 printf("In Cov1: %d\n", k);
2406 for (Int_t l=0; l<yBins; l++)
2410 // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
2411 for (Int_t j=0; j<yBins; j++)
2412 cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
2413 * (1.0 - vectorn[j] / nEvents);
2415 // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
2416 for (Int_t i=0; i<yBins; i++)
2417 for (Int_t j=0; j<yBins; j++)
2421 cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
2422 * vectorn[j] / nEvents;
2427 printf("Cov1 finished\n");
2429 TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
2432 for (Int_t i=1; i<=xBins; i++)
2433 for (Int_t j=1; j<=yBins; j++)
2434 cov->SetBinContent(i, j, cov1[i-1][j-1]);
2439 // second part of covariance matrix, depends on response matrix
2440 Float_t cov2[251][251];
2442 // Cov[P(Er|Cu), P(Es|Cu)] term
2443 Float_t covTerm[100][100][100];
2444 for (Int_t r=0; r<yBins; r++)
2445 for (Int_t u=0; u<xBins; u++)
2446 for (Int_t s=0; s<yBins; s++)
2449 covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2450 * (1.0 - hResponse->GetBinContent(u+1, r+1));
2452 covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2453 * hResponse->GetBinContent(u+1, s+1);
2456 for (Int_t k=0; k<xBins; k++)
2457 for (Int_t l=0; l<yBins; l++)
2460 printf("In Cov2: %d %d\n", k, l);
2461 for (Int_t i=0; i<yBins; i++)
2462 for (Int_t j=0; j<yBins; j++)
2464 //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
2465 // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
2467 for (Int_t r=0; r<yBins; r++)
2468 for (Int_t u=0; u<xBins; u++)
2469 for (Int_t s=0; s<yBins; s++)
2471 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
2472 || hResponse->GetBinContent(u+1, i+1) == 0)
2475 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
2476 * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
2480 cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
2484 printf("Cov2 finished\n");
2486 for (Int_t i=1; i<=xBins; i++)
2487 for (Int_t j=1; j<=yBins; j++)
2488 cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
2491 cov->Draw("COLZ");*/
2494 //____________________________________________________________________
2495 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
2498 // helper function for the covariance matrix of the bayesian method
2503 if (k == u && r == i)
2504 result += 1.0 / hResponse->GetBinContent(u+1, r+1);
2507 result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
2510 result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
2512 result *= matrixM[k][i];
2517 //____________________________________________________________________
2518 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
2521 // correct spectrum using bayesian method
2526 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2527 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2529 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2531 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2533 // TODO should be taken from correlation map
2534 //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2536 // normalize correction for given nPart
2537 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2539 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2540 //Double_t sum = sumHist->GetBinContent(i);
2543 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2546 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2547 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2552 fCurrentCorrelation->Draw("COLZ");
2555 fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
2557 // pick prior distribution
2558 TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
2559 Float_t norm = 1; //hPrior->Integral();
2560 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
2561 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
2563 // zero distribution
2564 TH1F* zero = (TH1F*)hPrior->Clone("zero");
2567 TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
2571 TH2F* hResponse = (TH2F*) fCurrentCorrelation;
2574 Int_t iterations = 25;
2575 for (Int_t i=0; i<iterations; i++)
2577 //printf(" iteration %i \n", i);
2579 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2582 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
2583 value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
2584 hTemp->SetBinContent(m, value);
2585 //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
2588 // regularization (simple smoothing)
2589 TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
2591 for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
2593 Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
2594 + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
2595 + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
2596 average *= hTrueSmoothed->GetBinWidth(t);
2598 // weight the average with the regularization parameter
2599 hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
2602 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2603 hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
2606 for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
2608 fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
2609 fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
2611 //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
2615 // calculate chi2 (change from last iteration)
2618 // use smoothed true (normalized) as new prior
2619 Float_t norm = 1; //hTrueSmoothed->Integral();
2621 for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
2623 Float_t newValue = hTemp->GetBinContent(t)/norm;
2624 Float_t diff = hPrior->GetBinContent(t) - newValue;
2625 chi2 += (Double_t) diff * diff;
2627 hPrior->SetBinContent(t, newValue);
2630 printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2633 DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2635 delete hTrueSmoothed;
2636 } // end of iterations
2638 DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2641 //____________________________________________________________________
2642 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
2645 // correct spectrum using a simple Gaussian approach, that is model-dependent
2648 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2649 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2651 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
2653 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2655 TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
2656 correction->SetTitle("GaussianMean");
2658 TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
2659 correction->SetTitle("GaussianWidth");
2661 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2663 TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
2664 proj->Fit("gaus", "0Q");
2665 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
2666 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
2670 // draw for debugging
2673 proj->GetFunction("gaus")->DrawCopy("SAME");
2676 TH1* target = fMultiplicityESDCorrected[correlationID];
2678 Int_t nBins = target->GetNbinsX()*10+1;
2679 Float_t* binning = new Float_t[nBins];
2680 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2681 for (Int_t j=0; j<10; ++j)
2682 binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2684 binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2686 TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2688 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2690 Float_t mean = correction->GetBinContent(i);
2691 Float_t width = correctionWidth->GetBinContent(i);
2693 Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2694 Int_t fillEnd = fineBinned->FindBin(mean + width * 5);
2695 //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2697 for (Int_t j=fillBegin; j <= fillEnd; ++j)
2699 fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2703 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2706 for (Int_t j=1; j<=10; ++j)
2707 sum += fineBinned->GetBinContent((i-1)*10 + j);
2708 target->SetBinContent(i, sum / 10);
2713 DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
2716 //____________________________________________________________________
2717 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
2719 // runs the distribution given in inputMC through the response matrix identified by
2720 // correlationMap and produces a measured distribution
2721 // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
2722 // if normalized is set, inputMC is expected to be normalized to the bin width
2727 if (correlationMap < 0 || correlationMap >= kCorrHists)
2730 // empty under/overflow bins in x, otherwise Project3D takes them into account
2731 TH3* hist = fCorrelation[correlationMap];
2732 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
2734 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
2736 hist->SetBinContent(0, y, z, 0);
2737 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2741 TH2* corr = (TH2*) hist->Project3D("zy");
2742 //corr->Rebin2D(2, 1);
2745 // normalize correction for given nPart
2746 for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2748 Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2752 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2755 corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2756 corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2760 TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2763 for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2765 Float_t measured = 0;
2768 for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2770 Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2772 measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2773 error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2776 //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2778 target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2779 target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2785 //____________________________________________________________________
2786 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2788 // uses the given function to fill the input MC histogram and generates from that
2789 // the measured histogram by applying the response matrix
2790 // this can be used to evaluate if the methods work indepedently of the input
2792 // WARNING does not respect the vertex distribution, just fills central vertex bin
2797 if (id < 0 || id >= kESDHists)
2800 TH2F* mc = fMultiplicityVtx[id];
2804 for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2806 mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
2807 mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
2813 TH1* proj = mc->ProjectionY();
2815 TString tmp(fMultiplicityESD[id]->GetName());
2816 delete fMultiplicityESD[id];
2817 fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
2818 fMultiplicityESD[id]->SetName(tmp);