PWG0/dNdEta/AliMultiplicityCorrection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // This class is used to store correction maps, raw input and results of the multiplicity
  19 // measurement with the ITS or TPC
  20 // It also contains functions to correct the spectrum using different methods.
  21 //
  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>
  34 #include <TF1.h>
  35 #include <TMath.h>
  36
  37 ClassImp(AliMultiplicityCorrection)
  38
  39 const Int_t AliMultiplicityCorrection::fgMaxParams = 90;
  40 TH1* AliMultiplicityCorrection::fCurrentESD = 0;
  41 TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
  42
  43 //____________________________________________________________________
  44 AliMultiplicityCorrection::AliMultiplicityCorrection() :
  45   TNamed()
  46 {
  47   //
  48   // default constructor
  49   //
  50
  51   for (Int_t i = 0; i < kESDHists; ++i)
  52     fMultiplicityESD[i] = 0;
  53
  54   for (Int_t i = 0; i < kMCHists; ++i)
  55     fMultiplicityMC[i] = 0;
  56
  57   for (Int_t i = 0; i < kCorrHists; ++i)
  58   {
  59     fCorrelation[i] = 0;
  60     fMultiplicityESDCorrected[i] = 0;
  61   }
  62 }
  63
  64 //____________________________________________________________________
  65 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
  66   TNamed(name, title)
  67 {
  68   //
  69   // named constructor
  70   //
  71
  72   // do not add this hists to the directory
  73   Bool_t oldStatus = TH1::AddDirectoryStatus();
  74   TH1::AddDirectory(kFALSE);
  75
  76   Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
  77   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,
  78                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
  79                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
  80                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
  81                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
  82                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
  83                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
  84                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
  85                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
  86                           500.5 };
  87                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
  88                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
  89                           //1000.5 };
  90
  91   #define BINNING fgMaxParams, binLimitsN
  92
  93   for (Int_t i = 0; i < kESDHists; ++i)
  94     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 10, binLimitsVtx, BINNING);
  95
  96   for (Int_t i = 0; i < kMCHists; ++i)
  97   {
  98     fMultiplicityMC[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityMC%d", i)));
  99     fMultiplicityMC[i]->SetTitle("fMultiplicityMC");
 100   }
 101
 102   for (Int_t i = 0; i < kCorrHists; ++i)
 103   {
 104     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 10, binLimitsVtx, BINNING, BINNING);
 105     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", BINNING);
 106   }
 107
 108   TH1::AddDirectory(oldStatus);
 109 }
 110
 111 //____________________________________________________________________
 112 AliMultiplicityCorrection::~AliMultiplicityCorrection()
 113 {
 114   //
 115   // Destructor
 116   //
 117 }
 118
 119 //____________________________________________________________________
 120 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
 121 {
 122   // Merge a list of AliMultiplicityCorrection objects with this (needed for
 123   // PROOF).
 124   // Returns the number of merged objects (including this).
 125
 126   if (!list)
 127     return 0;
 128
 129   if (list->IsEmpty())
 130     return 1;
 131
 132   TIterator* iter = list->MakeIterator();
 133   TObject* obj;
 134
 135   // collections of all histograms
 136   TList collections[kESDHists+kMCHists+kCorrHists*2];
 137
 138   Int_t count = 0;
 139   while ((obj = iter->Next())) {
 140
 141     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
 142     if (entry == 0)
 143       continue;
 144
 145     for (Int_t i = 0; i < kESDHists; ++i)
 146       collections[i].Add(entry->fMultiplicityESD[i]);
 147
 148     for (Int_t i = 0; i < kMCHists; ++i)
 149       collections[kESDHists+i].Add(entry->fMultiplicityMC[i]);
 150
 151     for (Int_t i = 0; i < kCorrHists; ++i)
 152       collections[kESDHists+kMCHists+i].Add(entry->fCorrelation[i]);
 153
 154     for (Int_t i = 0; i < kCorrHists; ++i)
 155       collections[kESDHists+kMCHists+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
 156
 157     count++;
 158   }
 159
 160   for (Int_t i = 0; i < kESDHists; ++i)
 161     fMultiplicityESD[i]->Merge(&collections[i]);
 162
 163   for (Int_t i = 0; i < kMCHists; ++i)
 164     fMultiplicityMC[i]->Merge(&collections[kESDHists+i]);
 165
 166   for (Int_t i = 0; i < kCorrHists; ++i)
 167     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists+i]);
 168
 169   for (Int_t i = 0; i < kCorrHists; ++i)
 170     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists+kCorrHists+i]);
 171
 172   delete iter;
 173
 174   return count+1;
 175 }
 176
 177 //____________________________________________________________________
 178 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
 179 {
 180   //
 181   // loads the histograms from a file
 182   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 183   //
 184
 185   if (!dir)
 186     dir = GetName();
 187
 188   if (!gDirectory->cd(dir))
 189     return kFALSE;
 190
 191   // TODO memory leak. old histograms needs to be deleted.
 192
 193   Bool_t success = kTRUE;
 194
 195   for (Int_t i = 0; i < kESDHists; ++i)
 196   {
 197     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
 198     if (!fMultiplicityESD[i])
 199       success = kFALSE;
 200   }
 201
 202   for (Int_t i = 0; i < kMCHists; ++i)
 203   {
 204     fMultiplicityMC[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMC[i]->GetName()));
 205     if (!fMultiplicityMC[i])
 206       success = kFALSE;
 207   }
 208
 209   for (Int_t i = 0; i < kCorrHists; ++i)
 210   {
 211     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
 212     if (!fCorrelation[i])
 213       success = kFALSE;
 214     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
 215     if (!fMultiplicityESDCorrected[i])
 216       success = kFALSE;
 217   }
 218
 219   gDirectory->cd("..");
 220
 221   return success;
 222 }
 223
 224 //____________________________________________________________________
 225 void AliMultiplicityCorrection::SaveHistograms()
 226 {
 227   //
 228   // saves the histograms
 229   //
 230
 231   gDirectory->mkdir(GetName());
 232   gDirectory->cd(GetName());
 233
 234   for (Int_t i = 0; i < kESDHists; ++i)
 235     if (fMultiplicityESD[i])
 236       fMultiplicityESD[i]->Write();
 237
 238   for (Int_t i = 0; i < kMCHists; ++i)
 239     if (fMultiplicityMC[i])
 240       fMultiplicityMC[i]->Write();
 241
 242   for (Int_t i = 0; i < kCorrHists; ++i)
 243   {
 244     if (fCorrelation[i])
 245       fCorrelation[i]->Write();
 246     if (fMultiplicityESDCorrected[i])
 247       fMultiplicityESDCorrected[i]->Write();
 248   }
 249
 250   gDirectory->cd("..");
 251 }
 252
 253 //____________________________________________________________________
 254 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
 255 {
 256   //
 257   // Fills an event from MC
 258   //
 259
 260   fMultiplicityMC[0]->Fill(vtx, generated05);
 261   fMultiplicityMC[1]->Fill(vtx, generated10);
 262   fMultiplicityMC[2]->Fill(vtx, generated15);
 263   fMultiplicityMC[3]->Fill(vtx, generated20);
 264   fMultiplicityMC[4]->Fill(vtx, generatedAll);
 265 }
 266
 267 //____________________________________________________________________
 268 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
 269 {
 270   //
 271   // Fills an event from ESD
 272   //
 273
 274   fMultiplicityESD[0]->Fill(vtx, measured05);
 275   fMultiplicityESD[1]->Fill(vtx, measured10);
 276   fMultiplicityESD[2]->Fill(vtx, measured15);
 277   fMultiplicityESD[3]->Fill(vtx, measured20);
 278 }
 279
 280 //____________________________________________________________________
 281 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)
 282 {
 283   //
 284   // Fills an event into the correlation map with the information from MC and ESD
 285   //
 286
 287   fCorrelation[0]->Fill(vtx, generated05, measured05);
 288   fCorrelation[1]->Fill(vtx, generated10, measured10);
 289   fCorrelation[2]->Fill(vtx, generated15, measured15);
 290   fCorrelation[3]->Fill(vtx, generated20, measured20);
 291
 292   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
 293   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
 294   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
 295   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
 296 }
 297
 298 //____________________________________________________________________
 299 Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
 300 {
 301   // homogenity term for minuit fitting
 302   // pure function of the parameters
 303   // prefers constant function (pol0)
 304
 305   Double_t chi2 = 0;
 306
 307   for (Int_t i=1; i<fgMaxParams; ++i)
 308   {
 309     if (params[i] == 0)
 310       continue;
 311
 312     Double_t right  = params[i] / fCurrentESD->GetBinWidth(i+1);
 313     Double_t left   = params[i-1] / fCurrentESD->GetBinWidth(i);
 314
 315     Double_t diff = (right - left) / right;
 316     chi2 += diff * diff;
 317   }
 318
 319   return chi2;
 320 }
 321
 322 //____________________________________________________________________
 323 Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
 324 {
 325   // homogenity term for minuit fitting
 326   // pure function of the parameters
 327   // prefers linear function (pol1)
 328
 329   Double_t chi2 = 0;
 330
 331   for (Int_t i=2; i<fgMaxParams; ++i)
 332   {
 333     if (params[i] == 0 || params[i-1] == 0)
 334       continue;
 335
 336     Double_t right  = params[i] / fCurrentESD->GetBinWidth(i+1);
 337     Double_t middle = params[i-1] / fCurrentESD->GetBinWidth(i);
 338     Double_t left   = params[i-2] / fCurrentESD->GetBinWidth(i-1);
 339
 340     Double_t der1 = (right - middle) / fCurrentESD->GetBinWidth(i);
 341     Double_t der2 = (middle - left)  / fCurrentESD->GetBinWidth(i-1);
 342
 343     Double_t diff = der1 - der2;
 344
 345     // TODO give additonal weight to big bins
 346     chi2 += diff * diff * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
 347
 348     //printf("%d --> %f\n", i, diff);
 349   }
 350
 351   return chi2 / 1e5 / 2;
 352 }
 353
 354 //____________________________________________________________________
 355 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params)
 356 {
 357   // homogenity term for minuit fitting
 358   // pure function of the parameters
 359   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
 360   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
 361
 362   Double_t chi2 = 0;
 363
 364   for (Int_t i=2; i<fgMaxParams; ++i)
 365   {
 366     if (params[i] == 0 || params[i-1] == 0)
 367       continue;
 368
 369     Double_t right  = params[i] / fCurrentESD->GetBinWidth(i+1);
 370     Double_t middle = params[i-1] / fCurrentESD->GetBinWidth(i);
 371     Double_t left   = params[i-2] / fCurrentESD->GetBinWidth(i-1);
 372
 373     Double_t der1 = (right - middle) / fCurrentESD->GetBinWidth(i);
 374     Double_t der2 = (middle - left)  / fCurrentESD->GetBinWidth(i-1);
 375
 376     Double_t secDer = (der1 - der2) / fCurrentESD->GetBinWidth(i);
 377
 378     // square and weight with the bin width
 379     chi2 += secDer * secDer * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
 380
 381     //printf("%d --> %f\n", i, secDer);
 382   }
 383
 384   return chi2 / 1e5;
 385 }
 386
 387 //____________________________________________________________________
 388 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
 389 {
 390   //
 391   // fit function for minuit
 392   //
 393
 394   // TODO take errors into account
 395
 396   static Int_t callCount = 0;
 397
 398   Double_t chi2FromFit = 0;
 399
 400   // loop over multiplicity (x axis of fMultiplicityESD)
 401   for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
 402   {
 403     if (i > fCurrentCorrelation->GetNbinsY())
 404       break;
 405
 406     Double_t sum = 0;
 407     //Double_t error = 0;
 408     // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks
 409     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j)
 410     {
 411       if (j > fgMaxParams)
 412         break;
 413
 414       sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1];
 415
 416       //if (params[j-1] > 0)
 417       //  error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1];
 418       //printf("%f  ", sum);
 419     }
 420
 421     Double_t diff = fCurrentESD->GetBinContent(i) - sum;
 422     if (fCurrentESD->GetBinContent(i) > 0)
 423       diff /= fCurrentESD->GetBinContent(i);
 424     else
 425       diff /= fCurrentESD->Integral();
 426
 427     // weight with bin width
 428     //diff *= fCurrentESD->GetBinWidth(i);
 429
 430     //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i);
 431     //if (error <= 0)
 432     // error = 1;
 433
 434     //Double_t tmp = diff / error;
 435     //chi2 += tmp * tmp;
 436     chi2FromFit += diff * diff;
 437
 438     //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff);
 439     //printf("Diff for bin %d is %f\n", i, diff);
 440   }
 441
 442   Double_t penaltyVal = RegularizationTotalCurvature(params);
 443
 444   chi2 = chi2FromFit * chi2FromFit + penaltyVal * penaltyVal;
 445
 446   if ((callCount++ % 100) == 0)
 447     printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
 448 }
 449
 450 //____________________________________________________________________
 451 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace)
 452 {
 453   //
 454   // fills fCurrentESD, fCurrentCorrelation
 455   // resets fMultiplicityESDCorrected
 456   //
 457
 458   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 459   fMultiplicityESDCorrected[correlationID]->Reset();
 460
 461   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
 462   fCurrentESD->Sumw2();
 463
 464   // empty under/overflow bins in x, otherwise Project3D takes them into account
 465   TH3* hist = fCorrelation[correlationID];
 466   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
 467   {
 468     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
 469     {
 470       hist->SetBinContent(0, y, z, 0);
 471       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
 472     }
 473   }
 474
 475   fCurrentCorrelation = hist->Project3D("zy");
 476   fCurrentCorrelation->Sumw2();
 477 }
 478
 479 //____________________________________________________________________
 480 void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check)
 481 {
 482   //
 483   // correct spectrum using minuit chi2 method
 484   //
 485   // if check is kTRUE the input MC solution (by definition the right solution) is used
 486   // no fit is made and just the chi2 is calculated
 487   //
 488
 489   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 490   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 491
 492   SetupCurrentHists(inputRange, fullPhaseSpace);
 493
 494   // normalize correction for given nPart
 495   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 496   {
 497     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
 498     if (sum <= 0)
 499       continue;
 500     Float_t maxValue = 0;
 501     Int_t maxBin = -1;
 502     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 503     {
 504       // find most probably value
 505       if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
 506       {
 507         maxValue = fCurrentCorrelation->GetBinContent(i, j);
 508         maxBin = j;
 509       }
 510
 511       // npart sum to 1
 512       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
 513       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
 514     }
 515
 516     printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
 517   }
 518
 519   // small hack to get around charge conservation for full phase space ;-)
 520   /*if (fullPhaseSpace)
 521   {
 522     for (Int_t i=2; i<=50; i+=2)
 523       for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 524         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
 525   }*/
 526
 527   TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
 528   canvas1->Divide(2, 1);
 529   canvas1->cd(1);
 530   fCurrentESD->DrawCopy();
 531   canvas1->cd(2);
 532   fCurrentCorrelation->DrawCopy("COLZ");
 533
 534   // Initialize TMinuit via generic fitter interface
 535   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
 536
 537   minuit->SetFCN(MinuitFitFunction);
 538
 539   Double_t results[fgMaxParams];
 540
 541   TH1* proj = fMultiplicityMC[mcTarget]->ProjectionY();
 542   for (Int_t i=0; i<fgMaxParams; ++i)
 543   {
 544     //results[i] = 1;
 545     results[i] = fCurrentESD->GetBinContent(i+1);
 546     if (check)
 547       results[i] = proj->GetBinContent(i+1);
 548     minuit->SetParameter(i, Form("param%d", i), results[i], ((results[i] > 1) ? (results[i] / 10) : 0), 0, fCurrentESD->GetMaximum() * 100);
 549   }
 550
 551   Int_t dummy;
 552   Double_t chi2;
 553   MinuitFitFunction(dummy, 0, chi2, results, 0);
 554   printf("Chi2 of right solution is = %f\n", chi2);
 555
 556   if (check)
 557     return;
 558
 559   Double_t arglist[100];
 560   arglist[0] = 100000;
 561   //minuit->ExecuteCommand("SET PRINT", arglist, 1);
 562   minuit->ExecuteCommand("MIGRAD", arglist, 1);
 563   //minuit->ExecuteCommand("MINOS", arglist, 0);
 564
 565   for (Int_t i=0; i<fgMaxParams; ++i)
 566   {
 567     results[i] = minuit->GetParameter(i);
 568     fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]);
 569     fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
 570   }
 571
 572   printf("Penalty is %f\n", RegularizationTotalCurvature(results));
 573
 574   DrawComparison("MinuitChi2", mcTarget, correlationID);
 575 }
 576
 577 //____________________________________________________________________
 578 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
 579 {
 580   //
 581   // normalizes a 1-d histogram to its bin width
 582   //
 583
 584   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
 585   {
 586     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
 587     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
 588   }
 589 }
 590
 591 //____________________________________________________________________
 592 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
 593 {
 594   //
 595   // normalizes a 2-d histogram to its bin width (x width * y width)
 596   //
 597
 598   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
 599     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
 600     {
 601       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
 602       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
 603       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
 604     }
 605 }
 606
 607 //____________________________________________________________________
 608 void AliMultiplicityCorrection::ApplyMinuitFitAll()
 609 {
 610   //
 611   // fit all eta ranges
 612   //
 613
 614   for (Int_t i=0; i<kESDHists; ++i)
 615   {
 616     ApplyMinuitFit(i, kFALSE);
 617     ApplyMinuitFit(i, kTRUE);
 618   }
 619 }
 620
 621 //____________________________________________________________________
 622 void AliMultiplicityCorrection::DrawHistograms()
 623 {
 624   printf("ESD:\n");
 625
 626   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
 627   canvas1->Divide(3, 2);
 628   for (Int_t i = 0; i < kESDHists; ++i)
 629   {
 630     canvas1->cd(i+1);
 631     fMultiplicityESD[i]->DrawCopy("COLZ");
 632     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
 633   }
 634
 635   printf("MC:\n");
 636
 637   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
 638   canvas2->Divide(3, 2);
 639   for (Int_t i = 0; i < kMCHists; ++i)
 640   {
 641     canvas2->cd(i+1);
 642     fMultiplicityMC[i]->DrawCopy("COLZ");
 643     printf("%d --> %f\n", i, fMultiplicityMC[i]->ProjectionY()->GetMean());
 644   }
 645
 646   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
 647   canvas3->Divide(3, 3);
 648   for (Int_t i = 0; i < kCorrHists; ++i)
 649   {
 650     canvas3->cd(i+1);
 651     TH1* proj = fCorrelation[i]->Project3D("zy");
 652     proj->DrawCopy("COLZ");
 653   }
 654 }
 655
 656 //____________________________________________________________________
 657 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t mcID, Int_t esdCorrId, Bool_t normalizeESD)
 658 {
 659   TString tmpStr;
 660   tmpStr.Form("%s_DrawComparison_%d_%d", name, mcID, esdCorrId);
 661
 662   TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 900, 300);
 663   canvas1->Divide(3, 1);
 664
 665   canvas1->cd(1);
 666   TH1* proj = fMultiplicityMC[mcID]->ProjectionY();
 667   NormalizeToBinWidth(proj);
 668
 669   if (normalizeESD)
 670     NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
 671
 672   proj->GetXaxis()->SetRangeUser(0, 200);
 673   proj->DrawCopy("HIST");
 674   gPad->SetLogy();
 675
 676   NormalizeToBinWidth(fCurrentESD);
 677   fCurrentESD->SetLineColor(2);
 678   fCurrentESD->DrawCopy("HISTSAME");
 679
 680   fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
 681   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAMEP");
 682
 683   canvas1->cd(2);
 684   fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, 200);
 685   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST");
 686
 687   canvas1->cd(3);
 688   TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
 689   clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
 690   clone->SetTitle("Ratio;Npart;MC/ESD");
 691   clone->GetYaxis()->SetRangeUser(0.8, 1.2);
 692   clone->GetXaxis()->SetRangeUser(0, 200);
 693   clone->DrawCopy("HIST");
 694
 695   /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
 696   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
 697   legend->AddEntry(fMultiplicityMC, "MC");
 698   legend->AddEntry(fMultiplicityESD, "ESD");
 699   legend->Draw();*/
 700 }
 701
 702 //____________________________________________________________________
 703 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
 704 {
 705   //
 706   // correct spectrum using bayesian method
 707   //
 708
 709   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 710   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 711
 712   SetupCurrentHists(inputRange, fullPhaseSpace);
 713
 714   // normalize correction for given nPart
 715   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 716   {
 717     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
 718     if (sum <= 0)
 719       continue;
 720     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
 721     {
 722       // npart sum to 1
 723       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
 724       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
 725     }
 726   }
 727
 728   Float_t regPar = 0.01;
 729
 730   Float_t norm = 0;
 731   // pick prior distribution
 732   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
 733   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
 734     norm = norm + hPrior->GetBinContent(t);
 735   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) {
 736     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
 737
 738     printf(" bin %d content %f \n", t, hPrior->GetBinContent(t));
 739
 740   }
 741
 742   // define hist to store guess of true
 743   TH1F* hTrue  = (TH1F*)fCurrentESD->Clone("prior");
 744   //  hTrue->Reset();
 745
 746   // clone response R
 747   TH2F* hResponse = (TH2F*)fCurrentCorrelation->Clone("pij");
 748
 749   // histogram to store the inverse response calculated from bayes theorem (from prior and response)
 750   // IR
 751   //TAxis* xAxis = hResponse->GetYaxis();
 752   //TAxis* yAxis = hResponse->GetXaxis();
 753
 754   TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
 755   //new TH2F("pji","pji",
 756   //                                     xAxis->GetNbins(),xAxis->GetXbins()->GetArray(),
 757   //                                     yAxis->GetNbins(),yAxis->GetXbins()->GetArray());
 758   hInverseResponseBayes->Reset();
 759
 760   // unfold...
 761   Int_t iterations = 50;
 762   for (Int_t i=0; i<iterations; i++) {
 763     printf(" iteration %i \n", i);
 764
 765     // calculate IR from Beyes theorem
 766     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
 767     for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
 768       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) {
 769
 770         norm = 0;
 771         for(Int_t t2 = 1; t2<=hResponse->GetNbinsX(); t2++)
 772           norm += (hResponse->GetBinContent(t2,m) * hPrior->GetBinContent(t2));
 773
 774         if (norm==0)
 775           hInverseResponseBayes->SetBinContent(t,m,0);
 776         else
 777           hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
 778       }
 779     }
 780     // calculate true assuming IR is the correct inverse response
 781     for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
 782       hTrue ->SetBinContent(t, 0);
 783       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
 784         hTrue->SetBinContent(t, hTrue->GetBinContent(t) + fCurrentESD->GetBinContent(m)*hInverseResponseBayes->GetBinContent(t,m));
 785     }
 786
 787     // regularization (simple smoothing) NB : not good bin width!!!
 788     TH1F* hTrueSmoothed = (TH1F*)hTrue->Clone("truesmoothed");
 789
 790     for (Int_t t=2; t<hTrue->GetNbinsX()-1; t++) {
 791       Float_t average = (hTrue->GetBinContent(t-1) / hTrue->GetBinWidth(t-1)
 792                          + hTrue->GetBinContent(t) / hTrue->GetBinWidth(t)
 793                          + hTrue->GetBinContent(t+1) / hTrue->GetBinWidth(t+1)) / 3.;
 794       average *= hTrueSmoothed->GetBinWidth(t);
 795
 796       // weight the average with the regularization parameter
 797       hTrueSmoothed->SetBinContent(t, (1-regPar)*hTrue->GetBinContent(t) + (regPar*average));
 798     }
 799
 800     // use smoothed true (normalized) as new prior
 801     norm = 0;
 802     for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
 803       norm = norm + hTrueSmoothed->GetBinContent(t);
 804     for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) {
 805       hPrior->SetBinContent(t, hTrueSmoothed->GetBinContent(t)/norm);
 806       hTrue->SetBinContent(t, hTrueSmoothed->GetBinContent(t));
 807     }
 808
 809     delete hTrueSmoothed;
 810   } // end of iterations
 811
 812   for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++) {
 813     fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTrue->GetBinContent(t));
 814     fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05*hTrue->GetBinContent(t));
 815
 816     printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
 817   }
 818
 819   DrawComparison("Bayesian", mcTarget, correlationID);
 820 }
 821
 822 //____________________________________________________________________
 823 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
 824 {
 825   //
 826   // correct spectrum using a simple Gaussian approach, that is model-dependent
 827   //
 828
 829   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 830   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 831
 832   SetupCurrentHists(inputRange, fullPhaseSpace);
 833
 834   NormalizeToBinWidth((TH2*) fCurrentCorrelation);
 835
 836   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
 837   correction->SetTitle("GaussianMean");
 838
 839   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
 840   correction->SetTitle("GaussianWidth");
 841
 842   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 843   {
 844     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
 845     proj->Fit("gaus", "0Q");
 846     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
 847     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
 848
 849     continue;
 850
 851     // draw for debugging
 852     new TCanvas;
 853     proj->DrawCopy();
 854     proj->GetFunction("gaus")->DrawCopy("SAME");
 855   }
 856
 857   TH1* target = fMultiplicityESDCorrected[correlationID];
 858
 859   Int_t nBins = target->GetNbinsX()*10+1;
 860   Float_t* binning = new Float_t[nBins];
 861   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
 862     for (Int_t j=0; j<10; ++j)
 863       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
 864
 865   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
 866
 867   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
 868
 869   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
 870   {
 871     Float_t mean = correction->GetBinContent(i);
 872     Float_t width = correctionWidth->GetBinContent(i);
 873
 874     Int_t fillBegin = fineBinned->FindBin(mean - width * 3);
 875     Int_t fillEnd   = fineBinned->FindBin(mean + width * 3);
 876     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
 877
 878     for (Int_t j=fillBegin; j <= fillEnd; ++j)
 879     {
 880       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
 881     }
 882   }
 883
 884   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
 885   {
 886     Float_t sum = 0;
 887     for (Int_t j=1; j<=10; ++j)
 888       sum += fineBinned->GetBinContent((i-1)*10 + j);
 889     target->SetBinContent(i, sum / 10);
 890   }
 891
 892   delete[] binning;
 893
 894   DrawComparison("Gaussian", mcTarget, correlationID, kFALSE);
 895 }