Modifications to take into account new muon esd structure (P. Cortese)
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityCorrection.cxx
CommitLineData
0a173978 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//
22// Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
23
24#include "AliMultiplicityCorrection.h"
25
26#include <TFile.h>
27#include <TH1F.h>
28#include <TH2F.h>
29#include <TH3F.h>
30#include <TDirectory.h>
31#include <TVirtualFitter.h>
32#include <TCanvas.h>
33#include <TString.h>
9ca6be09 34#include <TF1.h>
d560b581 35#include <TMath.h>
3602328d 36#include <TCollection.h>
0a173978 37
38ClassImp(AliMultiplicityCorrection)
39
cfc19dd5 40const Int_t AliMultiplicityCorrection::fgMaxParams = 200;
9ca6be09 41TH1* AliMultiplicityCorrection::fCurrentESD = 0;
42TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
cfc19dd5 43TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
44TMatrixF* AliMultiplicityCorrection::fCorrelationMatrix = 0;
45TMatrixF* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
46TVectorF* AliMultiplicityCorrection::fCurrentESDVector = 0;
47AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
48Float_t AliMultiplicityCorrection::fRegularizationWeight = 1e4;
dd701109 49TF1* AliMultiplicityCorrection::fNBD = 0;
0a173978 50
51//____________________________________________________________________
52AliMultiplicityCorrection::AliMultiplicityCorrection() :
dd701109 53 TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
0a173978 54{
55 //
56 // default constructor
57 //
58
59 for (Int_t i = 0; i < kESDHists; ++i)
60 fMultiplicityESD[i] = 0;
61
62 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 63 {
64 fMultiplicityVtx[i] = 0;
65 fMultiplicityMB[i] = 0;
66 fMultiplicityINEL[i] = 0;
67 }
0a173978 68
69 for (Int_t i = 0; i < kCorrHists; ++i)
70 {
71 fCorrelation[i] = 0;
72 fMultiplicityESDCorrected[i] = 0;
73 }
74}
75
76//____________________________________________________________________
77AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
dd701109 78 TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0)
0a173978 79{
80 //
81 // named constructor
82 //
83
84 // do not add this hists to the directory
85 Bool_t oldStatus = TH1::AddDirectoryStatus();
86 TH1::AddDirectory(kFALSE);
87
b477f4f2 88 /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
0a173978 89 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,
90 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
91 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
92 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
93 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
94 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
95 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
96 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
97 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
9ca6be09 98 500.5 };
99 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
100 //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
101 //1000.5 };
0a173978 102
b477f4f2 103 #define VTXBINNING 10, binLimitsVtx
104 #define NBINNING fgMaxParams, binLimitsN*/
105
106 #define NBINNING 251, -0.5, 250.5
107 #define VTXBINNING 10, -10, 10
0a173978 108
109 for (Int_t i = 0; i < kESDHists; ++i)
b477f4f2 110 fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
0a173978 111
112 for (Int_t i = 0; i < kMCHists; ++i)
113 {
cfc19dd5 114 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
115 fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
116
117 fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
118 fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
119
120 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
121 fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
0a173978 122 }
123
124 for (Int_t i = 0; i < kCorrHists; ++i)
125 {
b477f4f2 126 fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
127 fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
0a173978 128 }
129
130 TH1::AddDirectory(oldStatus);
131}
132
133//____________________________________________________________________
134AliMultiplicityCorrection::~AliMultiplicityCorrection()
135{
136 //
137 // Destructor
138 //
139}
140
141//____________________________________________________________________
142Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
143{
144 // Merge a list of AliMultiplicityCorrection objects with this (needed for
145 // PROOF).
146 // Returns the number of merged objects (including this).
147
148 if (!list)
149 return 0;
150
151 if (list->IsEmpty())
152 return 1;
153
154 TIterator* iter = list->MakeIterator();
155 TObject* obj;
156
157 // collections of all histograms
cfc19dd5 158 TList collections[kESDHists+kMCHists*3+kCorrHists*2];
0a173978 159
160 Int_t count = 0;
161 while ((obj = iter->Next())) {
162
163 AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
164 if (entry == 0)
165 continue;
166
167 for (Int_t i = 0; i < kESDHists; ++i)
168 collections[i].Add(entry->fMultiplicityESD[i]);
169
170 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 171 {
172 collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
173 collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
174 collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
175 }
0a173978 176
177 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 178 collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
0a173978 179
180 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 181 collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
0a173978 182
183 count++;
184 }
185
186 for (Int_t i = 0; i < kESDHists; ++i)
187 fMultiplicityESD[i]->Merge(&collections[i]);
188
189 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 190 {
191 fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
192 fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
193 fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
194 }
0a173978 195
196 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 197 fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
0a173978 198
199 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 200 fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
0a173978 201
202 delete iter;
203
204 return count+1;
205}
206
207//____________________________________________________________________
208Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
209{
210 //
211 // loads the histograms from a file
212 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
213 //
214
215 if (!dir)
216 dir = GetName();
217
218 if (!gDirectory->cd(dir))
219 return kFALSE;
220
221 // TODO memory leak. old histograms needs to be deleted.
222
223 Bool_t success = kTRUE;
224
225 for (Int_t i = 0; i < kESDHists; ++i)
226 {
227 fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
228 if (!fMultiplicityESD[i])
229 success = kFALSE;
230 }
231
232 for (Int_t i = 0; i < kMCHists; ++i)
233 {
cfc19dd5 234 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
235 fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
236 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
237 if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
0a173978 238 success = kFALSE;
239 }
240
241 for (Int_t i = 0; i < kCorrHists; ++i)
242 {
243 fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
244 if (!fCorrelation[i])
245 success = kFALSE;
246 fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
247 if (!fMultiplicityESDCorrected[i])
248 success = kFALSE;
249 }
250
251 gDirectory->cd("..");
252
253 return success;
254}
255
256//____________________________________________________________________
257void AliMultiplicityCorrection::SaveHistograms()
258{
259 //
260 // saves the histograms
261 //
262
263 gDirectory->mkdir(GetName());
264 gDirectory->cd(GetName());
265
266 for (Int_t i = 0; i < kESDHists; ++i)
267 if (fMultiplicityESD[i])
268 fMultiplicityESD[i]->Write();
269
270 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 271 {
272 if (fMultiplicityVtx[i])
273 fMultiplicityVtx[i]->Write();
274 if (fMultiplicityMB[i])
275 fMultiplicityMB[i]->Write();
276 if (fMultiplicityINEL[i])
277 fMultiplicityINEL[i]->Write();
278 }
0a173978 279
280 for (Int_t i = 0; i < kCorrHists; ++i)
281 {
282 if (fCorrelation[i])
283 fCorrelation[i]->Write();
284 if (fMultiplicityESDCorrected[i])
285 fMultiplicityESDCorrected[i]->Write();
286 }
287
288 gDirectory->cd("..");
289}
290
291//____________________________________________________________________
cfc19dd5 292void 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)
0a173978 293{
294 //
295 // Fills an event from MC
296 //
297
cfc19dd5 298 if (triggered)
299 {
300 fMultiplicityMB[0]->Fill(vtx, generated05);
301 fMultiplicityMB[1]->Fill(vtx, generated10);
302 fMultiplicityMB[2]->Fill(vtx, generated15);
303 fMultiplicityMB[3]->Fill(vtx, generated20);
304 fMultiplicityMB[4]->Fill(vtx, generatedAll);
305
306 if (vertex)
307 {
308 fMultiplicityVtx[0]->Fill(vtx, generated05);
309 fMultiplicityVtx[1]->Fill(vtx, generated10);
310 fMultiplicityVtx[2]->Fill(vtx, generated15);
311 fMultiplicityVtx[3]->Fill(vtx, generated20);
312 fMultiplicityVtx[4]->Fill(vtx, generatedAll);
313 }
314 }
315
316 fMultiplicityINEL[0]->Fill(vtx, generated05);
317 fMultiplicityINEL[1]->Fill(vtx, generated10);
318 fMultiplicityINEL[2]->Fill(vtx, generated15);
319 fMultiplicityINEL[3]->Fill(vtx, generated20);
320 fMultiplicityINEL[4]->Fill(vtx, generatedAll);
0a173978 321}
322
323//____________________________________________________________________
324void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
325{
326 //
327 // Fills an event from ESD
328 //
329
330 fMultiplicityESD[0]->Fill(vtx, measured05);
331 fMultiplicityESD[1]->Fill(vtx, measured10);
332 fMultiplicityESD[2]->Fill(vtx, measured15);
333 fMultiplicityESD[3]->Fill(vtx, measured20);
334}
335
336//____________________________________________________________________
337void 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)
338{
339 //
340 // Fills an event into the correlation map with the information from MC and ESD
341 //
342
343 fCorrelation[0]->Fill(vtx, generated05, measured05);
344 fCorrelation[1]->Fill(vtx, generated10, measured10);
345 fCorrelation[2]->Fill(vtx, generated15, measured15);
346 fCorrelation[3]->Fill(vtx, generated20, measured20);
347
348 fCorrelation[4]->Fill(vtx, generatedAll, measured05);
349 fCorrelation[5]->Fill(vtx, generatedAll, measured10);
350 fCorrelation[6]->Fill(vtx, generatedAll, measured15);
351 fCorrelation[7]->Fill(vtx, generatedAll, measured20);
352}
353
354//____________________________________________________________________
9ca6be09 355Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
0a173978 356{
357 // homogenity term for minuit fitting
358 // pure function of the parameters
359 // prefers constant function (pol0)
360
361 Double_t chi2 = 0;
362
363 for (Int_t i=1; i<fgMaxParams; ++i)
364 {
dd701109 365 Double_t right = params[i];
366 Double_t left = params[i-1];
0a173978 367
dd701109 368 Double_t diff = (right - left);
0a173978 369 chi2 += diff * diff;
370 }
371
dd701109 372 return chi2 * 100;
0a173978 373}
374
375//____________________________________________________________________
9ca6be09 376Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
0a173978 377{
378 // homogenity term for minuit fitting
379 // pure function of the parameters
380 // prefers linear function (pol1)
381
382 Double_t chi2 = 0;
383
384 for (Int_t i=2; i<fgMaxParams; ++i)
385 {
dd701109 386 if (params[i-1] == 0)
0a173978 387 continue;
388
dd701109 389 Double_t right = params[i];
390 Double_t middle = params[i-1];
391 Double_t left = params[i-2];
392
393 Double_t der1 = (right - middle);
394 Double_t der2 = (middle - left);
395
396 Double_t diff = (der1 - der2) / middle;
397
398 chi2 += diff * diff;
399 }
0a173978 400
dd701109 401 return chi2;
402}
0a173978 403
dd701109 404//____________________________________________________________________
405Double_t AliMultiplicityCorrection::RegularizationTest(Double_t *params)
406{
407 // homogenity term for minuit fitting
408 // pure function of the parameters
409 // prefers linear function (pol1)
0a173978 410
dd701109 411 Double_t chi2 = 0;
0a173978 412
dd701109 413 Float_t der[fgMaxParams];
414
415 for (Int_t i=0; i<fgMaxParams-1; ++i)
416 der[i] = params[i+1] - params[i];
417
418 for (Int_t i=0; i<fgMaxParams-6; ++i)
419 {
420 for (Int_t j=1; j<=5; ++j)
421 {
422 Double_t diff = der[i+j] - der[i];
423 chi2 += diff * diff;
424 }
0a173978 425 }
426
dd701109 427 return chi2 * 100; // 10000
0a173978 428}
429
9ca6be09 430//____________________________________________________________________
431Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params)
432{
433 // homogenity term for minuit fitting
434 // pure function of the parameters
435 // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
436 // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
437
438 Double_t chi2 = 0;
439
440 for (Int_t i=2; i<fgMaxParams; ++i)
441 {
dd701109 442 Double_t right = params[i];
443 Double_t middle = params[i-1];
444 Double_t left = params[i-2];
9ca6be09 445
dd701109 446 Double_t der1 = (right - middle);
447 Double_t der2 = (middle - left);
9ca6be09 448
dd701109 449 Double_t secDer = (der1 - der2);
9ca6be09 450
451 // square and weight with the bin width
dd701109 452 chi2 += secDer * secDer;
9ca6be09 453 }
454
dd701109 455 return chi2 * 100;
3602328d 456}
457
458//____________________________________________________________________
459Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params)
460{
461 // homogenity term for minuit fitting
462 // pure function of the parameters
463 // calculates entropy, from
464 // The method of reduced cross-entropy (M. Schmelling 1993)
465
3602328d 466 Double_t paramSum = 0;
467 for (Int_t i=0; i<fgMaxParams; ++i)
dd701109 468 paramSum += params[i];
3602328d 469
470 Double_t chi2 = 0;
471 for (Int_t i=0; i<fgMaxParams; ++i)
472 {
dd701109 473 Double_t tmp = params[i] / paramSum;
3602328d 474 if (tmp > 0)
dd701109 475 chi2 += tmp * log(tmp);
3602328d 476 }
3602328d 477
dd701109 478 return chi2;
479}
480
481//____________________________________________________________________
482void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
483{
484 //
485 // fit function for minuit
486 // does: nbd
487 // 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);
488 // func->SetParNames("scaling", "averagen", "k");
489 //
490
491 fNBD->SetParameters(params[0], params[1], params[2]);
492
493 Double_t params2[fgMaxParams];
494
495 for (Int_t i=0; i<fgMaxParams; ++i)
496 params2[i] = fNBD->Eval(i);
497
498 MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
499
500 printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
9ca6be09 501}
502
0a173978 503//____________________________________________________________________
504void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
505{
506 //
507 // fit function for minuit
cfc19dd5 508 // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
0a173978 509 //
510
0a173978 511 static Int_t callCount = 0;
512
513 Double_t chi2FromFit = 0;
514
cfc19dd5 515 // d
516 TVectorF paramsVector(fgMaxParams);
517 for (Int_t i=0; i<fgMaxParams; ++i)
518 paramsVector[i] = params[i];
519
520 // Ad
521 paramsVector = (*fCorrelationMatrix) * paramsVector;
522
523 // Ad - m
524 paramsVector -= (*fCurrentESDVector);
525
526 TVectorF copy(paramsVector);
527
528 // (Ad - m) W
529 paramsVector *= (*fCorrelationCovarianceMatrix);
530
531 // (Ad - m) W (Ad - m)
532 chi2FromFit = paramsVector * copy;
533
cfc19dd5 534 Double_t penaltyVal = 0;
535
536 switch (fRegularizationType)
537 {
dd701109 538 case kNone: break;
539 case kPol0: penaltyVal = RegularizationPol0(params); break;
540 case kPol1: penaltyVal = RegularizationPol1(params); break;
541 case kCurvature: penaltyVal = RegularizationTotalCurvature(params); break;
542 case kEntropy: penaltyVal = RegularizationEntropy(params); break;
543 case kTest: penaltyVal = RegularizationTest(params); break;
cfc19dd5 544 }
545
546 penaltyVal *= fRegularizationWeight;
0a173978 547
3602328d 548 chi2 = chi2FromFit + penaltyVal;
0a173978 549
cfc19dd5 550 if ((callCount++ % 1000) == 0)
0a173978 551 printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
552}
553
554//____________________________________________________________________
cfc19dd5 555void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
0a173978 556{
557 //
9ca6be09 558 // fills fCurrentESD, fCurrentCorrelation
559 // resets fMultiplicityESDCorrected
0a173978 560 //
561
562 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
9ca6be09 563 fMultiplicityESDCorrected[correlationID]->Reset();
0a173978 564
9ca6be09 565 fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
566 fCurrentESD->Sumw2();
0a173978 567
568 // empty under/overflow bins in x, otherwise Project3D takes them into account
569 TH3* hist = fCorrelation[correlationID];
570 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
571 {
572 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
573 {
574 hist->SetBinContent(0, y, z, 0);
575 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
576 }
577 }
578
9ca6be09 579 fCurrentCorrelation = hist->Project3D("zy");
580 fCurrentCorrelation->Sumw2();
cfc19dd5 581
582 fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency");
583 TH2* divisor = 0;
584 switch (eventType)
585 {
586 case kTrVtx : divisor = fMultiplicityVtx[inputRange]; break;
587 case kMB: divisor = fMultiplicityMB[inputRange]; break;
588 case kINEL: divisor = fMultiplicityINEL[inputRange]; break;
589 }
590 fCurrentEfficiency->Divide(divisor->ProjectionY());
9ca6be09 591}
592
593//____________________________________________________________________
dd701109 594void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
9ca6be09 595{
596 //
597 // correct spectrum using minuit chi2 method
598 //
599 // if check is kTRUE the input MC solution (by definition the right solution) is used
600 // no fit is made and just the chi2 is calculated
601 //
602
603 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
604 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
605
dd701109 606 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
cfc19dd5 607
608 fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
609 fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
610 fCurrentESDVector = new TVectorF(fgMaxParams);
0a173978 611
612 // normalize correction for given nPart
9ca6be09 613 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0a173978 614 {
9ca6be09 615 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
0a173978 616 if (sum <= 0)
617 continue;
9ca6be09 618 Float_t maxValue = 0;
619 Int_t maxBin = -1;
620 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0a173978 621 {
9ca6be09 622 // find most probably value
623 if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
624 {
625 maxValue = fCurrentCorrelation->GetBinContent(i, j);
626 maxBin = j;
627 }
628
0a173978 629 // npart sum to 1
9ca6be09 630 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
631 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
cfc19dd5 632
633 if (i <= fgMaxParams && j <= fgMaxParams)
634 (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
0a173978 635 }
9ca6be09 636
cfc19dd5 637 //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
0a173978 638 }
639
640 // small hack to get around charge conservation for full phase space ;-)
641 /*if (fullPhaseSpace)
642 {
643 for (Int_t i=2; i<=50; i+=2)
9ca6be09 644 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
645 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
0a173978 646 }*/
647
0a173978 648 // Initialize TMinuit via generic fitter interface
649 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
650
651 minuit->SetFCN(MinuitFitFunction);
652
3602328d 653 static Double_t results[fgMaxParams];
0a173978 654
dd701109 655 if (inputDist)
656 {
657 printf("Using different starting conditions...\n");
658 new TCanvas;
659 inputDist->DrawCopy();
660 }
661 else
662 inputDist = fCurrentESD;
663
664 // normalize ESD
665 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
666 inputDist->Scale(1.0 / inputDist->Integral());
667
668 Float_t minStartValue = 1e-3;
669
670 TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj");
671 proj->Scale(1.0 / proj->Integral());
0a173978 672 for (Int_t i=0; i<fgMaxParams; ++i)
673 {
dd701109 674 results[i] = inputDist->GetBinContent(i+1);
9ca6be09 675 if (check)
676 results[i] = proj->GetBinContent(i+1);
dd701109 677
678 printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
cfc19dd5 679
680 (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
681 if (fCurrentESD->GetBinError(i+1) > 0)
dd701109 682 (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / fCurrentESD->GetBinError(i+1);
cfc19dd5 683
dd701109 684 minuit->SetParameter(i, Form("param%d", i), ((results[i] > minStartValue) ? results[i] : minStartValue), 0.1, 0, 1);
0a173978 685 }
dd701109 686 // bin 0 is filled with value from bin 1 (otherwise it's 0)
687 minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 1);
0a173978 688
689 Int_t dummy;
690 Double_t chi2;
691 MinuitFitFunction(dummy, 0, chi2, results, 0);
dd701109 692 printf("Chi2 of initial parameters is = %f\n", chi2);
0a173978 693
9ca6be09 694 if (check)
695 return;
0a173978 696
697 Double_t arglist[100];
3602328d 698 arglist[0] = 0;
699 minuit->ExecuteCommand("SET PRINT", arglist, 1);
0a173978 700 minuit->ExecuteCommand("MIGRAD", arglist, 1);
dd701109 701 //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
0a173978 702 //minuit->ExecuteCommand("MINOS", arglist, 0);
703
704 for (Int_t i=0; i<fgMaxParams; ++i)
705 {
dd701109 706 if (fCurrentEfficiency->GetBinContent(i+1) > 0)
707 {
708 fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
709 fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) / fCurrentEfficiency->GetBinContent(i+1));
710 }
711 else
712 {
713 fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, 0);
714 fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
715 }
716 }
717}
718
719//____________________________________________________________________
720void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
721{
722 //
723 // correct spectrum using minuit chi2 method applying a NBD fit
724 //
725
726 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
727
728 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
729
730 fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
731 fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
732 fCurrentESDVector = new TVectorF(fgMaxParams);
733
734 // normalize correction for given nPart
735 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
736 {
737 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
738 if (sum <= 0)
739 continue;
740 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
741 {
742 // npart sum to 1
743 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
744 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
745
746 if (i <= fgMaxParams && j <= fgMaxParams)
747 (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
748 }
749 }
750
751 for (Int_t i=0; i<fgMaxParams; ++i)
752 {
753 (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
754 if (fCurrentESD->GetBinError(i+1) > 0)
755 (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
756 }
757
758 // Create NBD function
759 if (!fNBD)
760 {
761 fNBD = 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);
762 fNBD->SetParNames("scaling", "averagen", "k");
0a173978 763 }
764
dd701109 765 // Initialize TMinuit via generic fitter interface
766 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
767
768 minuit->SetFCN(MinuitNBD);
769 minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
770 minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
771 minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
772
773 Double_t arglist[100];
774 arglist[0] = 0;
775 minuit->ExecuteCommand("SET PRINT", arglist, 1);
776 minuit->ExecuteCommand("MIGRAD", arglist, 1);
777 //minuit->ExecuteCommand("MINOS", arglist, 0);
778
779 fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
780
781 for (Int_t i=0; i<fgMaxParams; ++i)
782 {
783 printf("%d : %f\n", i, fNBD->Eval(i));
784 if (fNBD->Eval(i) > 0)
785 {
786 fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fNBD->Eval(i));
787 fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
788 }
789 }
0a173978 790}
791
792//____________________________________________________________________
793void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
794{
795 //
796 // normalizes a 1-d histogram to its bin width
797 //
798
799 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
800 {
801 hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
802 hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
803 }
804}
805
806//____________________________________________________________________
807void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
808{
809 //
810 // normalizes a 2-d histogram to its bin width (x width * y width)
811 //
812
813 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
814 for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
815 {
816 Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
817 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
818 hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
819 }
820}
821
0a173978 822//____________________________________________________________________
823void AliMultiplicityCorrection::DrawHistograms()
824{
825 printf("ESD:\n");
826
827 TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
828 canvas1->Divide(3, 2);
829 for (Int_t i = 0; i < kESDHists; ++i)
830 {
831 canvas1->cd(i+1);
832 fMultiplicityESD[i]->DrawCopy("COLZ");
833 printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
834 }
835
cfc19dd5 836 printf("Vtx:\n");
0a173978 837
838 TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
839 canvas2->Divide(3, 2);
840 for (Int_t i = 0; i < kMCHists; ++i)
841 {
842 canvas2->cd(i+1);
cfc19dd5 843 fMultiplicityVtx[i]->DrawCopy("COLZ");
844 printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
0a173978 845 }
846
847 TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
848 canvas3->Divide(3, 3);
849 for (Int_t i = 0; i < kCorrHists; ++i)
850 {
851 canvas3->cd(i+1);
b477f4f2 852 TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
853 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
854 {
855 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
856 {
857 hist->SetBinContent(0, y, z, 0);
858 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
859 }
860 }
861 TH1* proj = hist->Project3D("zy");
0a173978 862 proj->DrawCopy("COLZ");
863 }
864}
865
866//____________________________________________________________________
cfc19dd5 867void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist)
0a173978 868{
cfc19dd5 869 Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
870
9ca6be09 871 TString tmpStr;
cfc19dd5 872 tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
0a173978 873
b477f4f2 874 // calculate residual
875
876 // for that we convolute the response matrix with the gathered result
877 // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
dd701109 878 TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
879 tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
880 TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId, !normalizeESD);
b477f4f2 881 TH1* proj2 = convoluted->ProjectionY("proj2", -1, -1, "e");
882 NormalizeToBinWidth(proj2);
dd701109 883
884 TH1* residual = (TH1*) proj2->Clone("residual");
885 residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / measured");
b477f4f2 886
887 residual->Add(fCurrentESD, -1);
888 residual->Divide(residual, fCurrentESD, 1, 1, "B");
889
890 // TODO fix errors
891 for (Int_t i=1; i<residual->GetNbinsX(); ++i)
892 {
893 proj2->SetBinError(i, 0);
894 residual->SetBinError(i, 0);
895 }
896
dd701109 897 // normalize mc to 1
898 //mcHist->Scale(1.0 / mcHist->Integral(1, 200));
899
900 // normalize unfolded esd to 1
901 //fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(1, 200));
902
903 // normalize unfolded result to mc hist
904 fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
905 fMultiplicityESDCorrected[esdCorrId]->Scale(mcHist->Integral(1, 200));
906
cfc19dd5 907 TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 800);
b477f4f2 908 canvas1->Divide(2, 2);
0a173978 909
910 canvas1->cd(1);
cfc19dd5 911 TH1* proj = (TH1*) mcHist->Clone("proj");
0a173978 912 NormalizeToBinWidth(proj);
9ca6be09 913
914 if (normalizeESD)
915 NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
916
dd701109 917 proj->GetXaxis()->SetRangeUser(0, 200);
cfc19dd5 918 proj->DrawCopy("HIST");
0a173978 919 gPad->SetLogy();
920
cfc19dd5 921 //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
922 fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
923 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST");
924
925 Float_t chi2 = 0;
dd701109 926 fLastChi2MCLimit = 0;
927 for (Int_t i=2; i<=200; ++i)
cfc19dd5 928 {
dd701109 929 if (proj->GetBinContent(i) != 0)
cfc19dd5 930 {
dd701109 931 Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
cfc19dd5 932 chi2 += value * value;
dd701109 933
934 if (chi2 < 0.1)
935 fLastChi2MCLimit = i;
936
937 if (i == 100)
938 fLastChi2MC = chi2;
cfc19dd5 939 }
940 }
941
dd701109 942 printf("Difference (from MC) is %f for bin 2-100. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
cfc19dd5 943
944 canvas1->cd(2);
945
9ca6be09 946 NormalizeToBinWidth(fCurrentESD);
cfc19dd5 947 gPad->SetLogy();
dd701109 948 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
cfc19dd5 949 //fCurrentESD->SetLineColor(2);
950 fCurrentESD->DrawCopy("HIST");
b477f4f2 951
952 proj2->SetLineColor(2);
cfc19dd5 953 //proj2->SetMarkerColor(2);
954 //proj2->SetMarkerStyle(5);
955 proj2->DrawCopy("HIST SAME");
9ca6be09 956
cfc19dd5 957 chi2 = 0;
dd701109 958 for (Int_t i=2; i<=100; ++i)
cfc19dd5 959 {
960 if (fCurrentESD->GetBinContent(i) != 0)
961 {
962 Float_t value = (proj2->GetBinContent(i) - fCurrentESD->GetBinContent(i)) / fCurrentESD->GetBinContent(i);
963 chi2 += value * value;
964 }
965 }
0a173978 966
dd701109 967 printf("Difference (Residuals) is %f for bin 2-100\n", chi2);
cfc19dd5 968 fLastChi2Residuals = chi2;
0a173978 969
970 canvas1->cd(3);
971 TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
972 clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
dd701109 973 clone->SetTitle("Ratio;Npart;MC/Unfolded Measured");
0a173978 974 clone->GetYaxis()->SetRangeUser(0.8, 1.2);
dd701109 975 clone->GetXaxis()->SetRangeUser(0, 200);
0a173978 976 clone->DrawCopy("HIST");
977
978 /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
979 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
980 legend->AddEntry(fMultiplicityMC, "MC");
981 legend->AddEntry(fMultiplicityESD, "ESD");
982 legend->Draw();*/
3602328d 983
b477f4f2 984 canvas1->cd(4);
b477f4f2 985 residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
dd701109 986 residual->GetXaxis()->SetRangeUser(0, 200);
b477f4f2 987 residual->DrawCopy();
3602328d 988
b477f4f2 989 canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
0a173978 990}
991
992//____________________________________________________________________
dd701109 993void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals)
cfc19dd5 994{
995 // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
996 // These values are computed during DrawComparison, thus this function picks up the
997 // last calculation
998
dd701109 999 if (mc)
1000 *mc = fLastChi2MC;
1001 if (mcLimit)
1002 *mcLimit = fLastChi2MCLimit;
1003 if (residuals)
1004 *residuals = fLastChi2Residuals;
cfc19dd5 1005}
1006
1007
1008//____________________________________________________________________
1009TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1010{
1011 //
1012 // returns the corresponding MC spectrum
1013 //
1014
1015 switch (eventType)
1016 {
1017 case kTrVtx : return fMultiplicityVtx[i]; break;
1018 case kMB : return fMultiplicityMB[i]; break;
1019 case kINEL : return fMultiplicityINEL[i]; break;
1020 }
1021
1022 return 0;
1023}
1024
1025//____________________________________________________________________
dd701109 1026void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist)
0a173978 1027{
1028 //
1029 // correct spectrum using bayesian method
1030 //
1031
1032 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
0a173978 1033
cfc19dd5 1034 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1035
dd701109 1036 // normalize ESD
1037 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
0a173978 1038
1039 // normalize correction for given nPart
9ca6be09 1040 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0a173978 1041 {
dd701109 1042 // with this it is normalized to 1
1043 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1044
1045 // with this normalized to the given efficiency
1046 if (fCurrentEfficiency->GetBinContent(i) > 0)
1047 sum /= fCurrentEfficiency->GetBinContent(i);
1048 else
1049 sum = 0;
1050
9ca6be09 1051 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0a173978 1052 {
dd701109 1053 if (sum > 0)
1054 {
1055 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1056 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1057 }
1058 else
1059 {
1060 fCurrentCorrelation->SetBinContent(i, j, 0);
1061 fCurrentCorrelation->SetBinError(i, j, 0);
1062 }
0a173978 1063 }
1064 }
1065
dd701109 1066 // smooth input spectrum
1067 /*
1068 TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone");
cfc19dd5 1069
dd701109 1070 for (Int_t i=2; i<fCurrentESD->GetNbinsX(); ++i)
1071 if (esdClone->GetBinContent(i) < 1e-3)
1072 fCurrentESD->SetBinContent(i, (esdClone->GetBinContent(i-1) + esdClone->GetBinContent(i) + esdClone->GetBinContent(i+1)) / 3);
cfc19dd5 1073
dd701109 1074 delete esdClone;
1075
1076 // rescale to 1
1077 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1078 */
6127aca6 1079
6127aca6 1080 // pick prior distribution
dd701109 1081 TH1* hPrior = 0;
1082 if (inputDist)
1083 {
1084 printf("Using different starting conditions...\n");
1085 hPrior = (TH1*)inputDist->Clone("prior");
1086 }
1087 else
1088 hPrior = (TH1*)fCurrentESD->Clone("prior");
cfc19dd5 1089 Float_t norm = hPrior->Integral();
6127aca6 1090 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
6127aca6 1091 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1092
cfc19dd5 1093 // define temp hist
1094 TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1095 hTemp->Reset();
6127aca6 1096
cfc19dd5 1097 // just a shortcut
1098 TH2F* hResponse = (TH2F*) fCurrentCorrelation;
6127aca6 1099
cfc19dd5 1100 // histogram to store the inverse response calculated from bayes theorem (from prior and response) IR
6127aca6 1101 TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
6127aca6 1102 hInverseResponseBayes->Reset();
b477f4f2 1103
dd701109 1104 TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
1105
6127aca6 1106 // unfold...
dd701109 1107 for (Int_t i=0; i<nIterations; i++)
cfc19dd5 1108 {
1109 //printf(" iteration %i \n", i);
1110
6127aca6 1111 // calculate IR from Beyes theorem
1112 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
9ca6be09 1113
cfc19dd5 1114 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1115 {
1116 Float_t norm = 0;
1117 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1118 norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t);
1119 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1120 {
6127aca6 1121 if (norm==0)
cfc19dd5 1122 hInverseResponseBayes->SetBinContent(t,m,0);
6127aca6 1123 else
cfc19dd5 1124 hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
6127aca6 1125 }
1126 }
cfc19dd5 1127
1128 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1129 {
1130 Float_t value = 0;
6127aca6 1131 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
cfc19dd5 1132 value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m);
1133
dd701109 1134 if (fCurrentEfficiency->GetBinContent(t) > 0)
1135 hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t));
1136 else
1137 hTemp->SetBinContent(t, 0);
cfc19dd5 1138 }
1139
1140 // this is the last guess, fill before (!) smoothing
1141 for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1142 {
1143 fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1144 fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1145
1146 //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
6127aca6 1147 }
9ca6be09 1148
b477f4f2 1149 // regularization (simple smoothing)
cfc19dd5 1150 TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
9ca6be09 1151
cfc19dd5 1152 for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1153 {
dd701109 1154 Float_t average = (hTemp->GetBinContent(t-1)
1155 + hTemp->GetBinContent(t)
1156 + hTemp->GetBinContent(t+1)
1157 ) / 3.;
9ca6be09 1158
6127aca6 1159 // weight the average with the regularization parameter
cfc19dd5 1160 hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
6127aca6 1161 }
9ca6be09 1162
cfc19dd5 1163 // calculate chi2 (change from last iteration)
1164 Double_t chi2 = 0;
1165
9ca6be09 1166 // use smoothed true (normalized) as new prior
cfc19dd5 1167 Float_t norm = 1;
1168 //for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1169 // norm = norm + hTrueSmoothed->GetBinContent(t);
1170
1171 for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1172 {
1173 Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
1174 Float_t diff = hPrior->GetBinContent(t) - newValue;
1175 chi2 += (Double_t) diff * diff;
1176
1177 hPrior->SetBinContent(t, newValue);
6127aca6 1178 }
9ca6be09 1179
cfc19dd5 1180 //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1181
dd701109 1182 convergence->Fill(i+1, chi2);
1183
cfc19dd5 1184 //if (i % 5 == 0)
1185 // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
1186
9ca6be09 1187 delete hTrueSmoothed;
6127aca6 1188 } // end of iterations
1189
dd701109 1190 //new TCanvas;
1191 //convergence->DrawCopy();
1192 //gPad->SetLogy();
1193
1194 delete convergence;
6127aca6 1195
cfc19dd5 1196 return;
1197
1198 // ********
1199 // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
1200
dd701109 1201 /*printf("Calculating covariance matrix. This may take some time...\n");
1202
1203 // TODO check if this is the right one...
1204 TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
cfc19dd5 1205
1206 Int_t xBins = hInverseResponseBayes->GetNbinsX();
1207 Int_t yBins = hInverseResponseBayes->GetNbinsY();
1208
1209 // calculate "unfolding matrix" Mij
1210 Float_t matrixM[251][251];
1211 for (Int_t i=1; i<=xBins; i++)
1212 {
1213 for (Int_t j=1; j<=yBins; j++)
1214 {
1215 if (fCurrentEfficiency->GetBinContent(i) > 0)
1216 matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
1217 else
1218 matrixM[i-1][j-1] = 0;
1219 }
1220 }
1221
1222 Float_t* vectorn = new Float_t[yBins];
1223 for (Int_t j=1; j<=yBins; j++)
1224 vectorn[j-1] = fCurrentESD->GetBinContent(j);
1225
1226 // first part of covariance matrix, depends on input distribution n(E)
1227 Float_t cov1[251][251];
1228
1229 Float_t nEvents = fCurrentESD->Integral(); // N
1230
1231 xBins = 20;
1232 yBins = 20;
1233
1234 for (Int_t k=0; k<xBins; k++)
1235 {
1236 printf("In Cov1: %d\n", k);
1237 for (Int_t l=0; l<yBins; l++)
1238 {
1239 cov1[k][l] = 0;
1240
1241 // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
1242 for (Int_t j=0; j<yBins; j++)
1243 cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
1244 * (1.0 - vectorn[j] / nEvents);
1245
1246 // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
1247 for (Int_t i=0; i<yBins; i++)
1248 for (Int_t j=0; j<yBins; j++)
1249 {
1250 if (i == j)
1251 continue;
1252 cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
1253 * vectorn[j] / nEvents;
1254 }
1255 }
1256 }
1257
1258 printf("Cov1 finished\n");
1259
1260 TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
1261 cov->Reset();
1262
1263 for (Int_t i=1; i<=xBins; i++)
1264 for (Int_t j=1; j<=yBins; j++)
1265 cov->SetBinContent(i, j, cov1[i-1][j-1]);
1266
1267 new TCanvas;
1268 cov->Draw("COLZ");
1269
1270 // second part of covariance matrix, depends on response matrix
1271 Float_t cov2[251][251];
1272
1273 // Cov[P(Er|Cu), P(Es|Cu)] term
1274 Float_t covTerm[100][100][100];
1275 for (Int_t r=0; r<yBins; r++)
1276 for (Int_t u=0; u<xBins; u++)
1277 for (Int_t s=0; s<yBins; s++)
1278 {
1279 if (r == s)
1280 covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1281 * (1.0 - hResponse->GetBinContent(u+1, r+1));
1282 else
1283 covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1284 * hResponse->GetBinContent(u+1, s+1);
1285 }
1286
1287 for (Int_t k=0; k<xBins; k++)
1288 for (Int_t l=0; l<yBins; l++)
1289 {
1290 cov2[k][l] = 0;
1291 printf("In Cov2: %d %d\n", k, l);
1292 for (Int_t i=0; i<yBins; i++)
1293 for (Int_t j=0; j<yBins; j++)
1294 {
1295 //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
1296 // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
1297 Float_t tmpCov = 0;
1298 for (Int_t r=0; r<yBins; r++)
1299 for (Int_t u=0; u<xBins; u++)
1300 for (Int_t s=0; s<yBins; s++)
1301 {
1302 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
1303 || hResponse->GetBinContent(u+1, i+1) == 0)
1304 continue;
1305
1306 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
1307 * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
1308 * covTerm[r][u][s];
1309 }
1310
1311 cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
1312 }
1313 }
1314
1315 printf("Cov2 finished\n");
1316
1317 for (Int_t i=1; i<=xBins; i++)
1318 for (Int_t j=1; j<=yBins; j++)
1319 cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
1320
1321 new TCanvas;
dd701109 1322 cov->Draw("COLZ");*/
1323}
1324
1325//____________________________________________________________________
1326Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse,
1327 TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u)
1328{
1329 //
1330 // helper function for the covariance matrix of the bayesian method
1331 //
1332
1333 Float_t result = 0;
1334
1335 if (k == u && r == i)
1336 result += 1.0 / hResponse->GetBinContent(u+1, r+1);
1337
1338 if (k == u)
1339 result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
1340
1341 if (r == i)
1342 result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
1343
1344 result *= matrixM[k][i];
1345
1346 return result;
cfc19dd5 1347}
1348
1349//____________________________________________________________________
1350void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
1351{
1352 //
1353 // correct spectrum using bayesian method
1354 //
1355
dd701109 1356 Float_t regPar = 0;
cfc19dd5 1357
1358 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1359 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1360
1361 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1362
1363 // TODO should be taken from correlation map
1364 //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1365
1366 // normalize correction for given nPart
1367 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1368 {
1369 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1370 //Double_t sum = sumHist->GetBinContent(i);
1371 if (sum <= 0)
1372 continue;
1373 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1374 {
1375 // npart sum to 1
1376 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1377 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1378 }
6127aca6 1379 }
cfc19dd5 1380
1381 new TCanvas;
1382 fCurrentCorrelation->Draw("COLZ");
1383
1384 // FAKE
1385 fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1386
1387 // pick prior distribution
1388 TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1389 Float_t norm = 1; //hPrior->Integral();
1390 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1391 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1392
1393 // zero distribution
1394 TH1F* zero = (TH1F*)hPrior->Clone("zero");
1395
1396 // define temp hist
1397 TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1398 hTemp->Reset();
1399
1400 // just a shortcut
1401 TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1402
1403 // unfold...
dd701109 1404 Int_t iterations = 25;
cfc19dd5 1405 for (Int_t i=0; i<iterations; i++)
1406 {
1407 //printf(" iteration %i \n", i);
1408
1409 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1410 {
1411 Float_t value = 0;
1412 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1413 value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
1414 hTemp->SetBinContent(m, value);
1415 //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
1416 }
1417
1418 // regularization (simple smoothing)
1419 TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1420
1421 for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1422 {
1423 Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1424 + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1425 + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1426 average *= hTrueSmoothed->GetBinWidth(t);
1427
1428 // weight the average with the regularization parameter
1429 hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1430 }
1431
1432 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1433 hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
1434
1435 // fill guess
1436 for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1437 {
1438 fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1439 fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1440
1441 //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1442 }
1443
1444
1445 // calculate chi2 (change from last iteration)
1446 Double_t chi2 = 0;
1447
1448 // use smoothed true (normalized) as new prior
1449 Float_t norm = 1; //hTrueSmoothed->Integral();
1450
1451 for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1452 {
dd701109 1453 Float_t newValue = hTemp->GetBinContent(t)/norm;
cfc19dd5 1454 Float_t diff = hPrior->GetBinContent(t) - newValue;
1455 chi2 += (Double_t) diff * diff;
1456
1457 hPrior->SetBinContent(t, newValue);
1458 }
1459
1460 printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1461
dd701109 1462 //if (i % 5 == 0)
cfc19dd5 1463 DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1464
1465 delete hTrueSmoothed;
1466 } // end of iterations
1467
1468 DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1469}
1470
9ca6be09 1471//____________________________________________________________________
1472void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
1473{
1474 //
1475 // correct spectrum using a simple Gaussian approach, that is model-dependent
1476 //
1477
1478 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1479 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1480
cfc19dd5 1481 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
9ca6be09 1482
1483 NormalizeToBinWidth((TH2*) fCurrentCorrelation);
1484
1485 TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
1486 correction->SetTitle("GaussianMean");
1487
1488 TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
1489 correction->SetTitle("GaussianWidth");
1490
1491 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1492 {
1493 TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
1494 proj->Fit("gaus", "0Q");
1495 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1496 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1497
1498 continue;
1499
1500 // draw for debugging
1501 new TCanvas;
1502 proj->DrawCopy();
1503 proj->GetFunction("gaus")->DrawCopy("SAME");
1504 }
1505
1506 TH1* target = fMultiplicityESDCorrected[correlationID];
1507
1508 Int_t nBins = target->GetNbinsX()*10+1;
1509 Float_t* binning = new Float_t[nBins];
1510 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1511 for (Int_t j=0; j<10; ++j)
1512 binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
1513
1514 binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
1515
1516 TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
1517
1518 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1519 {
1520 Float_t mean = correction->GetBinContent(i);
1521 Float_t width = correctionWidth->GetBinContent(i);
1522
3602328d 1523 Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
1524 Int_t fillEnd = fineBinned->FindBin(mean + width * 5);
1525 //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
9ca6be09 1526
1527 for (Int_t j=fillBegin; j <= fillEnd; ++j)
1528 {
1529 fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
1530 }
1531 }
1532
1533 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1534 {
1535 Float_t sum = 0;
1536 for (Int_t j=1; j<=10; ++j)
1537 sum += fineBinned->GetBinContent((i-1)*10 + j);
1538 target->SetBinContent(i, sum / 10);
1539 }
1540
1541 delete[] binning;
1542
cfc19dd5 1543 DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
0a173978 1544}
3602328d 1545
1546//____________________________________________________________________
b477f4f2 1547TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized)
3602328d 1548{
1549 // runs the distribution given in inputMC through the response matrix identified by
1550 // correlationMap and produces a measured distribution
1551 // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
b477f4f2 1552 // if normalized is set, inputMC is expected to be normalized to the bin width
3602328d 1553
1554 if (!inputMC)
1555 return 0;
1556
1557 if (correlationMap < 0 || correlationMap >= kCorrHists)
1558 return 0;
1559
1560 // empty under/overflow bins in x, otherwise Project3D takes them into account
1561 TH3* hist = fCorrelation[correlationMap];
1562 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1563 {
1564 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1565 {
1566 hist->SetBinContent(0, y, z, 0);
1567 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1568 }
1569 }
1570
1571 TH1* corr = hist->Project3D("zy");
1572 corr->Sumw2();
1573
1574 // normalize correction for given nPart
1575 for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
1576 {
1577 Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
1578 if (sum <= 0)
1579 continue;
1580
1581 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1582 {
1583 // npart sum to 1
1584 corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
1585 corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
1586 }
1587 }
1588
1589 TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
1590 target->Reset();
1591
1592 for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
1593 {
1594 Float_t measured = 0;
1595 Float_t error = 0;
1596
1597 for (Int_t gen=1; gen<=target->GetNbinsY(); ++gen)
1598 {
1599 Int_t mcGenBin = inputMC->GetXaxis()->FindBin(target->GetYaxis()->GetBinCenter(gen));
1600
b477f4f2 1601 Float_t factor = 1;
1602 if (normalized)
1603 factor = inputMC->GetBinWidth(mcGenBin);
1604
1605 measured += inputMC->GetBinContent(mcGenBin) * factor * corr->GetBinContent(gen, meas);
1606 error += inputMC->GetBinError(mcGenBin) * factor * corr->GetBinContent(gen, meas);
3602328d 1607 }
1608
1609 // bin width of the <measured> axis has to be taken into account
b477f4f2 1610 //measured /= target->GetYaxis()->GetBinWidth(meas);
1611 //error /= target->GetYaxis()->GetBinWidth(meas);
3602328d 1612
cfc19dd5 1613 //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
3602328d 1614
1615 target->SetBinContent(target->GetNbinsX() / 2, meas, measured);
1616 target->SetBinError(target->GetNbinsX() / 2, meas, error);
1617 }
1618
1619 return target;
1620}
1621
1622//____________________________________________________________________
1623void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
1624{
1625 // uses the given function to fill the input MC histogram and generates from that
1626 // the measured histogram by applying the response matrix
1627 // this can be used to evaluate if the methods work indepedently of the input
1628 // distribution
1629 // WARNING does not respect the vertex distribution, just fills central vertex bin
1630
1631 if (!inputMC)
1632 return;
1633
1634 if (id < 0 || id >= kESDHists)
1635 return;
1636
dd701109 1637 TH2F* mc = fMultiplicityVtx[id];
3602328d 1638
1639 mc->Reset();
1640
1641 for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
b477f4f2 1642 {
3602328d 1643 mc->SetBinContent(mc->GetNbinsX() / 2, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
b477f4f2 1644 mc->SetBinError(mc->GetNbinsX() / 2, i, 0);
1645 }
3602328d 1646
1647 //new TCanvas;
1648 //mc->Draw("COLZ");
1649
1650 TH1* proj = mc->ProjectionY();
1651
dd701109 1652 TString tmp(fMultiplicityESD[id]->GetName());
1653 delete fMultiplicityESD[id];
3602328d 1654 fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
dd701109 1655 fMultiplicityESD[id]->SetName(tmp);
3602328d 1656}