PWG4/JetTasks/AliJetSpectrumUnfolding.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 // This class is used to store correlation maps, generated and reconstructed data of the jet spectrum
  17 // It also contains functions to correct the spectrum using the bayesian unfolding
  18 //
  19
  20 #include "AliJetSpectrumUnfolding.h"
  21
  22 #include <TFile.h>
  23 #include <TH1F.h>
  24 #include <TH2F.h>
  25 #include <TH3F.h>
  26 #include <TDirectory.h>
  27 #include <TVirtualFitter.h>
  28 #include <TCanvas.h>
  29 #include <TString.h>
  30 #include <TF1.h>
  31 #include <TF2.h>
  32 #include <TMath.h>
  33 #include <TCollection.h>
  34 #include <TLegend.h>
  35 #include <TLine.h>
  36 #include <TRandom.h>
  37 #include <TProfile.h>
  38 #include <TProfile2D.h>
  39 #include <TStyle.h>
  40 #include <TColor.h>
  41
  42 #include <THnSparse.h>
  43
  44 ClassImp(AliJetSpectrumUnfolding)
  45
  46 const Int_t NBINSE  = 50;
  47 const Int_t NBINSZ  = 50;
  48 const Int_t NEVENTS = 500000;
  49 const Double_t axisLowerLimitE = 0.;
  50 const Double_t axisLowerLimitZ = 0.;
  51 const Double_t axisUpperLimitE = 250.;
  52 const Double_t axisUpperLimitZ = 1.;
  53
  54 Float_t AliJetSpectrumUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
  55 Int_t   AliJetSpectrumUnfolding::fgBayesianIterations = 100;         // number of iterations in Bayesian method
  56
  57 //____________________________________________________________________
  58
  59 AliJetSpectrumUnfolding::AliJetSpectrumUnfolding() :
  60   TNamed(), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
  61   fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
  62 {
  63   //
  64   // default constructor
  65   //
  66
  67   fGenSpectrum = 0;
  68   fRecSpectrum = 0;
  69   fUnfSpectrum = 0;
  70   fCorrelation = 0;
  71 }
  72
  73 //____________________________________________________________________
  74 AliJetSpectrumUnfolding::AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title) :
  75   TNamed(name, title), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
  76   fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
  77 {
  78   //
  79   // named constructor
  80   //
  81
  82   // do not add this hists to the directory
  83   Bool_t oldStatus = TH1::AddDirectoryStatus();
  84   TH1::AddDirectory(kFALSE);
  85
  86   #define ZBINNING NBINSZ, axisLowerLimitZ, axisUpperLimitZ       // fragmentation of leading particle
  87   #define EBINNING NBINSE, axisLowerLimitE, axisUpperLimitE       // energy of the jet
  88
  89   fRecSpectrum = new TH2F("fRecSpectrum", "Reconstructed Spectrum;E^{jet}_{rec} [GeV];z^{lp}_{rec}", EBINNING, ZBINNING);
  90   fGenSpectrum = new TH2F("fGenSpectrum", "Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}_{gen}", EBINNING, ZBINNING);
  91   fUnfSpectrum = new TH2F("fUnfSpectrum", "Unfolded Spectrum;E^{jet} [GeV];z^{lp}", EBINNING, ZBINNING);
  92
  93   const Int_t nbin[4]={NBINSE, NBINSE, NBINSZ, NBINSZ};
  94   //arrays for bin limits
  95   Double_t LowEdge[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
  96   Double_t UpEdge[4]  = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
  97
  98   fCorrelation = new THnSparseF("fCorrelation", "Correlation Function", 4, nbin, LowEdge, UpEdge);
  99
 100   TH1::AddDirectory(oldStatus);
 101 }
 102
 103 //____________________________________________________________________
 104 AliJetSpectrumUnfolding::~AliJetSpectrumUnfolding()
 105 {
 106   //
 107   // Destructor
 108   //
 109
 110   if (fGenSpectrum)
 111     delete fGenSpectrum;
 112   fGenSpectrum = 0;
 113
 114   if (fRecSpectrum)
 115     delete fRecSpectrum;
 116   fRecSpectrum = 0;
 117
 118   if (fUnfSpectrum)
 119     delete fUnfSpectrum;
 120   fUnfSpectrum = 0;
 121
 122   if (fCorrelation)
 123     delete fCorrelation;
 124   fCorrelation = 0;
 125
 126 }
 127
 128 //____________________________________________________________________
 129 Long64_t AliJetSpectrumUnfolding::Merge(TCollection* list)
 130 {
 131   // Merge a list of AliJetSpectrumUnfolding objects with this (needed for
 132   // PROOF).
 133   // Returns the number of merged objects (including this).
 134
 135   if (!list)
 136     return 0;
 137
 138   if (list->IsEmpty())
 139     return 1;
 140
 141   TIterator* iter = list->MakeIterator();
 142   TObject* obj;
 143
 144   // collections of all histograms
 145   TList collections[4];
 146
 147   Int_t count = 0;
 148   while ((obj = iter->Next())) {
 149
 150     AliJetSpectrumUnfolding* entry = dynamic_cast<AliJetSpectrumUnfolding*> (obj);
 151     if (entry == 0)
 152       continue;
 153
 154     collections[0].Add(entry->fGenSpectrum);
 155     collections[1].Add(entry->fRecSpectrum);
 156     collections[2].Add(entry->fUnfSpectrum);
 157     collections[3].Add(entry->fCorrelation);
 158
 159     count++;
 160   }
 161
 162   fGenSpectrum->Merge(&collections[0]);
 163   fRecSpectrum->Merge(&collections[1]);
 164   fUnfSpectrum->Merge(&collections[2]);
 165   fCorrelation->Merge(&collections[3]);
 166
 167   delete iter;
 168
 169   return count+1;
 170 }
 171
 172 //____________________________________________________________________
 173 Bool_t AliJetSpectrumUnfolding::LoadHistograms(const Char_t* dir)
 174 {
 175   //
 176   // loads the histograms from a file
 177   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 178   //
 179
 180   if (!dir)
 181     dir = GetName();
 182
 183   if (!gDirectory->cd(dir))
 184     return kFALSE;
 185
 186   Bool_t success = kTRUE;
 187
 188   // store old histograms to delete them later
 189   TList oldHistograms;
 190   oldHistograms.SetOwner(1);
 191
 192   if (fGenSpectrum)  oldHistograms.Add(fGenSpectrum);
 193   if (fRecSpectrum)  oldHistograms.Add(fRecSpectrum);
 194   if (fUnfSpectrum)  oldHistograms.Add(fUnfSpectrum);
 195   if (fCorrelation)  oldHistograms.Add(fCorrelation);
 196
 197   // load new histograms
 198   fGenSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fGenSpectrum->GetName()));
 199   if (!fGenSpectrum)
 200     success = kFALSE;
 201
 202   fRecSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fRecSpectrum->GetName()));
 203   if (!fRecSpectrum)
 204     success = kFALSE;
 205
 206   fUnfSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fUnfSpectrum->GetName()));
 207   if (!fUnfSpectrum)
 208     success = kFALSE;
 209
 210   fCorrelation = dynamic_cast<THnSparseF*> (gDirectory->Get(fCorrelation->GetName()));
 211   if (!fCorrelation)
 212     success = kFALSE;
 213
 214   gDirectory->cd("..");
 215
 216   // delete old histograms
 217   oldHistograms.Delete();
 218
 219   return success;
 220 }
 221
 222 //____________________________________________________________________
 223 void AliJetSpectrumUnfolding::SaveHistograms()
 224 {
 225   //
 226   // saves the histograms
 227   //
 228
 229   gDirectory->mkdir(GetName());
 230   gDirectory->cd(GetName());
 231
 232   if (fGenSpectrum)
 233     fGenSpectrum->Write();
 234
 235   if (fRecSpectrum)
 236     fRecSpectrum->Write();
 237
 238   if (fUnfSpectrum)
 239     fUnfSpectrum->Write();
 240
 241   if (fCorrelation)
 242     fCorrelation->Write();
 243
 244   gDirectory->cd("..");
 245 }
 246
 247 //____________________________________________________________________
 248 void AliJetSpectrumUnfolding::SetupCurrentHists(Bool_t createBigBin)
 249 {
 250   //
 251   // resets fUnfSpectrum
 252   //
 253
 254   fUnfSpectrum->Reset();
 255   fUnfSpectrum->Sumw2();
 256
 257   fCurrentRec = (TH2F*)fRecSpectrum->Clone("fCurrentRec");
 258   fCurrentRec->Sumw2();
 259
 260   fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation");
 261   fCurrentCorrelation->Sumw2();
 262
 263   Printf("Correlation Matrix has %.0E filled bins", fCurrentCorrelation->GetNbins());
 264
 265   if (createBigBin)
 266   {
 267     Int_t maxBinE = 0, maxBinZ = 0;
 268     Float_t maxE = 0, maxZ = 0;
 269     for (Int_t me=1; me<=fCurrentRec->GetNbinsX(); me++)
 270       for (Int_t mz=1; mz<=fCurrentRec->GetNbinsY(); mz++)
 271       {
 272         if (fCurrentRec->GetBinContent(me,mz) <= 5 && me>NBINSE/2 && mz>NBINSZ/2)
 273         {
 274           maxBinE = me;
 275           maxBinZ = mz;
 276           maxE = fCurrentRec->GetXaxis()->GetBinCenter(me);
 277           maxZ = fCurrentRec->GetYaxis()->GetBinCenter(mz);
 278           break;
 279         }
 280       }
 281
 282     if (maxBinE > 0 || maxBinZ > 0)
 283     {
 284       printf("Bin limit in measured spectrum is e = %d and z = %d.\n", maxBinE, maxBinZ);
 285       fCurrentRec->SetBinContent(maxBinE, maxBinZ, fCurrentRec->Integral(maxBinE, fCurrentRec->GetNbinsX(), maxBinZ, fCurrentRec->GetNbinsY()));
 286       for (Int_t me=maxBinE+1; me<=fCurrentRec->GetNbinsX(); me++)
 287         for (Int_t mz=maxBinZ+1; mz<=fCurrentRec->GetNbinsY(); mz++)
 288         {
 289           fCurrentRec->SetBinContent(me, mz, 0);
 290           fCurrentRec->SetBinError(me, mz, 0);
 291         }
 292       // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
 293       fCurrentRec->SetBinError(maxBinE, maxBinZ, TMath::Sqrt(fCurrentRec->GetBinContent(maxBinE, maxBinZ)));
 294
 295       printf("This bin has now %f +- %f entries\n", fCurrentRec->GetBinContent(maxBinE, maxBinZ), fCurrentRec->GetBinError(maxBinE, maxBinZ));
 296
 297     /*  for (Int_t te=1; te<=NBINSE; te++)
 298       {
 299         for (Int_t tz=1; tz<=NBINSZ; tz++)
 300         {
 301           Int_t binMin[4] = {te, maxBinE, tz, maxBinZ};
 302           Int_t binMax[4] = {NBINSE, NBINSE, NBINSZ, NBINSZ};
 303           Float_t sum=0;
 304           for (Int_t ite=te; ite<=NBINSE; ite++)
 305             for (Int_t itz=tz; itz<=NBINSZ; itz++)
 306               for (Int_t ime=maxBinE; ime<=NBINSE; ime++)
 307                 for (Int_t imz=maxBinZ; imz<=NBINSZ; imz++)
 308                 {
 309                   Int_t bin[4] = {ite, ime, itz, imz};
 310                   sum += fCurrentCorrelation->GetBinContent(bin);
 311                 }
 312           fCurrentCorrelation->SetBinContent(binMin, sum);
 313           fCurrentCorrelation->SetBinError(binMin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
 314           printf("create big bin1, nbins = %d, te  = %d, tz = %d\n", NBINSE, te, tz);
 315           for (Int_t me=maxBinE; me<=NBINSE; me++)
 316           {
 317             for (Int_t mz=maxBinZ; mz<=NBINSZ; mz++)
 318             {
 319               Int_t bin[4] = {te, me, tz, mz};
 320               fCurrentCorrelation->SetBinContent(bin, 0.);
 321               fCurrentCorrelation->SetBinError(bin, 0.);
 322               printf("create big bin2\n");
 323             }
 324           }
 325         }
 326       }*/
 327
 328       for(Int_t idx = 0; idx<=fCurrentCorrelation->GetNbins(); idx++)
 329       {
 330         Int_t bin[4];
 331         Float_t binContent = fCurrentCorrelation->GetBinContent(idx,bin);
 332         Float_t binError   = fCurrentCorrelation->GetBinError(idx);
 333         Int_t binMin[4] = {bin[0], maxBinE, bin[2], maxBinZ};
 334         if ( (bin[1]>maxBinE && bin[1]<=NBINSE) && (bin[3]>maxBinZ && bin[3]<=NBINSZ) )
 335         {
 336           fCurrentCorrelation->SetBinContent(binMin, binContent + fCurrentCorrelation->GetBinContent(binMin));
 337           fCurrentCorrelation->SetBinError(binMin, binError + TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
 338           fCurrentCorrelation->SetBinContent(bin, 0.);
 339           fCurrentCorrelation->SetBinError(bin, 0.);
 340         }
 341         printf("create big bin1, nbins = %d, te  = %d, tz = %d\n", NBINSE, bin[0], bin[1]);
 342       }
 343
 344       printf("Adjusted correlation matrix!\n");
 345     }
 346   } // end Create Big Bin
 347
 348 }
 349
 350 //____________________________________________________________________
 351 void AliJetSpectrumUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 352 {
 353   //
 354   // sets the parameters for Bayesian unfolding
 355   //
 356
 357   fgBayesianSmoothing = smoothing;
 358   fgBayesianIterations = nIterations;
 359
 360   printf("AliJetSpectrumUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
 361 }
 362
 363 //____________________________________________________________________
 364 void AliJetSpectrumUnfolding::NormalizeToBinWidth(TH2* hist)
 365 {
 366   //
 367   // normalizes a 2-d histogram to its bin width (x width * y width)
 368   //
 369
 370   for (Int_t i=1; i<=hist->GetNbinsX(); i++)
 371     for (Int_t j=1; j<=hist->GetNbinsY(); j++)
 372     {
 373       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
 374       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
 375       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
 376     }
 377 }
 378
 379 //____________________________________________________________________
 380 void AliJetSpectrumUnfolding::DrawHistograms()
 381 {
 382   //
 383   // draws the histograms of this class
 384   //
 385
 386   gStyle->SetPalette(1);
 387
 388   TCanvas* canvas1 = new TCanvas("fRecSpectrum", "fRecSpectrum", 900, 600);
 389   gPad->SetLogz();
 390   fRecSpectrum->DrawCopy("COLZ");
 391
 392   TCanvas* canvas2 = new TCanvas("fGenSpectrum", "fGenSpectrum", 900, 600);
 393   gPad->SetLogz();
 394   fGenSpectrum->DrawCopy("COLZ");
 395
 396   TCanvas* canvas3 = new TCanvas("fUnfSpectrum", "fUnfSpectrum", 900, 600);
 397   gPad->SetLogz();
 398   fUnfSpectrum->DrawCopy("COLZ");
 399
 400   TCanvas* canvas4 = new TCanvas("fCorrelation", "fCorrelation", 500, 500);
 401   canvas1->Divide(2);
 402
 403   canvas4->cd(1);
 404   gPad->SetLogz();
 405   TH2D* h0 = fCorrelation->Projection(1,0);
 406   h0->SetXTitle("E^{jet}_{gen} [GeV]");
 407   h0->SetYTitle("E^{jet}_{rec} [GeV]");
 408   h0->SetTitle("Projection: Jet Energy");
 409   h0->DrawCopy("colz");
 410
 411   canvas1->cd(2);
 412   gPad->SetLogz();
 413   TH2D* h1 = fCorrelation->Projection(3,2);
 414   h1->SetXTitle("z^{lp}_{gen}");
 415   h1->SetYTitle("z^{lp}_{rec}");
 416   h1->SetTitle("Projection: Leading Particle Fragmentation");
 417   h1->DrawCopy("colz");
 418
 419 }
 420
 421 //____________________________________________________________________
 422 void AliJetSpectrumUnfolding::DrawComparison(const char* name, TH2* genHist)
 423 {
 424
 425   if (fUnfSpectrum->Integral() == 0)
 426   {
 427     printf("ERROR. Unfolded histogram is empty\n");
 428     return;
 429   }
 430
 431   //regain measured distribution used for unfolding, because the bins were modified in SetupCurrentHists
 432   //in create big bin
 433   fCurrentRec = (TH2F*)fRecSpectrum->Clone();
 434   fCurrentRec->Sumw2();
 435   fCurrentRec->Scale(1.0 / fCurrentRec->Integral());
 436
 437   // normalize unfolded result to 1
 438   fUnfSpectrum->Scale(1.0 / fUnfSpectrum->Integral());
 439
 440   // find bin with <= 100 entries. this is used as limit for the chi2 calculation
 441   Int_t mcBinLimitE = 0, mcBinLimitZ = 0;
 442   for (Int_t i=0; i<genHist->GetNbinsX(); ++i)
 443     for (Int_t j=0; j<genHist->GetNbinsY(); ++j)
 444     {
 445       if (genHist->GetBinContent(i,j) > 100)
 446       {
 447         mcBinLimitE = i;
 448         mcBinLimitZ = j;
 449       }
 450       else
 451         break;
 452     }
 453   Printf("AliJetSpectrumUnfolding::DrawComparison: Gen bin limit is (x,y) = (%d,%d)", mcBinLimitE,mcBinLimitZ);
 454
 455   // scale to 1 true spectrum
 456   genHist->Sumw2();
 457   genHist->Scale(1.0 / genHist->Integral());
 458
 459   // calculate residual
 460   // for that we convolute the response matrix with the gathered result
 461   TH2* tmpRecRecorrected = (TH2*) fUnfSpectrum->Clone("tmpRecRecorrected");
 462   TH2* convoluted = CalculateRecSpectrum(tmpRecRecorrected);
 463   if (convoluted->Integral() > 0)
 464     convoluted->Scale(1.0 / convoluted->Integral());
 465   else
 466     printf("ERROR: convoluted is empty. Something went wrong calculating the convoluted histogram.\n");
 467
 468   TH2* residual = (TH2*) convoluted->Clone("residual");
 469   residual->SetTitle("(R#otimesUnfolded - Reconstructed)/Reconstructed;E^{jet} [GeV]; z^{lp}");
 470
 471   fCurrentRec->Scale(1./fCurrentRec->Integral());
 472   residual->Add(fCurrentRec, -1);
 473   //residual->Divide(residual, fCurrentRec, 1, 1, "B");
 474
 475   // draw canvas
 476   TCanvas* canvas1 = new TCanvas(name, name, 1000, 1000);
 477   canvas1->Divide(2, 3);
 478
 479   Int_t style = 1;
 480   const Int_t NRGBs = 5;
 481   const Int_t NCont = 500;
 482
 483   Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
 484   Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
 485   Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
 486   Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
 487   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
 488   gStyle->SetNumberContours(NCont);
 489
 490   canvas1->cd(1);
 491   gStyle->SetPalette(style);
 492   gPad->SetLogz();
 493   genHist->SetTitle("Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}");
 494   genHist->SetStats(0);
 495   genHist->DrawCopy("colz");
 496
 497   canvas1->cd(2);
 498   gStyle->SetPalette(style);
 499   gPad->SetLogz();
 500   fUnfSpectrum->SetStats(0);
 501   fUnfSpectrum->DrawCopy("colz");
 502
 503   canvas1->cd(3);
 504   gStyle->SetPalette(style);
 505   gPad->SetLogz();
 506   fCurrentRec->SetTitle(fRecSpectrum->GetTitle());
 507   fCurrentRec->SetStats(0);
 508   fCurrentRec->DrawCopy("colz");
 509
 510   canvas1->cd(4);
 511   gStyle->SetPalette(style);
 512   gPad->SetLogy();
 513   TH1D* projGenX = genHist->ProjectionX();
 514   projGenX->SetTitle("Projection: Jet Energy; E^{jet} [GeV]");
 515   TH1D* projUnfX = fUnfSpectrum->ProjectionX();
 516   TH1D* projRecX = fCurrentRec->ProjectionX();
 517   projGenX->SetStats(0);
 518   projRecX->SetStats(0);
 519   projUnfX->SetStats(0);
 520   projGenX->SetLineColor(8);
 521   projRecX->SetLineColor(2);
 522   projGenX->DrawCopy();
 523   projUnfX->DrawCopy("same");
 524   projRecX->DrawCopy("same");
 525
 526   TLegend* legend = new TLegend(0.6, 0.85, 0.98, 0.98);
 527   legend->AddEntry(projGenX, "Generated Spectrum");
 528   legend->AddEntry(projUnfX, "Unfolded Spectrum");
 529   legend->AddEntry(projRecX, "Reconstructed Spectrum");
 530   //legend->SetFillColor(0);
 531   legend->Draw("same");
 532
 533   canvas1->cd(5);
 534   gPad->SetLogy();
 535   gStyle->SetPalette(style);
 536   TH1D* projGenY = genHist->ProjectionY();
 537   projGenY->SetTitle("Projection: Leading Particle Fragmentation; z^{lp}");
 538   TH1D* projUnfY = fUnfSpectrum->ProjectionY();
 539   TH1D* projRecY = fCurrentRec->ProjectionY();
 540   projGenY->SetStats(0);
 541   projRecY->SetStats(0);
 542   projUnfY->SetStats(0);
 543   projGenY->SetLineColor(8);
 544   projRecY->SetLineColor(2);
 545   projGenY->DrawCopy();
 546   projUnfY->DrawCopy("same");
 547   projRecY->DrawCopy("same");
 548
 549   TLegend* legend1 = new TLegend(0.6, 0.85, 0.98, 0.98);
 550   legend1->AddEntry(projGenY, "Generated Spectrum");
 551   legend1->AddEntry(projUnfY, "Unfolded Spectrum");
 552   legend1->AddEntry(projRecY, "Recontructed Spectrum");
 553   //legend1->SetFillColor(0);
 554   legend1->Draw("same");
 555
 556   // Draw residuals
 557   canvas1->cd(6);
 558   gStyle->SetPalette(style);
 559   gPad->SetLogz();
 560   residual->SetStats(0);
 561   residual->DrawCopy("colz");
 562
 563   canvas1->SaveAs(Form("%s.png", canvas1->GetName()));
 564 }
 565
 566
 567 //____________________________________________________________________
 568 void AliJetSpectrumUnfolding::GetComparisonResults(Float_t* gen, Int_t* genLimit, Float_t* residuals, Float_t* ratioAverage) const
 569 {
 570   // Returns the chi2 between the Generated and the unfolded Reconstructed spectrum as well as between the Reconstructed and the folded unfolded
 571   // These values are computed during DrawComparison, thus this function picks up the
 572   // last calculation
 573
 574   if (gen)
 575     *gen = fLastChi2MC;
 576   if (genLimit)
 577     *genLimit = fLastChi2MCLimit;
 578   if (residuals)
 579     *residuals = fLastChi2Residuals;
 580   if (ratioAverage)
 581     *ratioAverage = fRatioAverage;
 582 }
 583
 584 //____________________________________________________________________
 585 void AliJetSpectrumUnfolding::ApplyBayesianMethod(Float_t regPar, Int_t nIterations, TH2* initialConditions, Bool_t determineError)
 586 {
 587   //
 588   // correct spectrum using bayesian unfolding
 589   //
 590
 591   // initialize seed with current time
 592   gRandom->SetSeed(0);
 593
 594   printf("seting up current arrays and histograms...\n");
 595   SetupCurrentHists(kFALSE); // kFALSE to not create big bin
 596
 597   // normalize Correlation Map to convert number of events into probabilities
 598   /*for (Int_t te=1; te<=NBINSE; te++)
 599     for (Int_t tz=1; tz<=NBINSZ; tz++)
 600     {
 601        Int_t bin[4];
 602        Float_t sum=0.;
 603        for (Int_t me = 1; me<=NBINSE; me++)
 604          for (Int_t mz = 1; mz<=NBINSZ; mz++)
 605          {
 606            bin[0] = te; bin[1] = me;
 607            bin[2] = tz; bin[3] = mz;
 608            sum += fCurrentCorrelation->GetBinContent(bin);
 609          }
 610        if (sum > 0.)
 611          for (Int_t me = 1; me<=NBINSE; me++)
 612            for (Int_t mz = 1; mz<=NBINSZ; mz++)
 613            {
 614              bin[0] = te; bin[1] = me;
 615              bin[2] = tz; bin[3] = mz;
 616              fCurrentCorrelation->SetBinContent(bin, fCurrentCorrelation->GetBinContent(bin)/sum);
 617              fCurrentCorrelation->SetBinError(bin, fCurrentCorrelation->GetBinError(bin)/sum);
 618            }
 619     }*/
 620   Float_t sum[NBINSE+2][NBINSZ+2];
 621   memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
 622
 623   for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
 624   {
 625     Int_t bin[4];
 626     Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
 627     if ( (bin[1]>0 && bin[1]<=NBINSE) && (bin[3]>0 && bin[3]<=NBINSZ) )
 628       sum[bin[0]][bin[2]] += binContent;
 629   }
 630
 631   for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
 632   {
 633     Int_t bin[4];
 634     Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
 635     Float_t binError   = fCurrentCorrelation->GetBinError(bin);
 636     if (sum[bin[0]][bin[2]]>0 && (bin[1]>0 && bin[1]<=NBINSE) &&
 637         (bin[3]>0 && bin[3]<=NBINSZ) && (bin[0]>0 && bin[0]<=NBINSE) && (bin[2]>0 && bin[2]<=NBINSZ) )
 638     {
 639       fCurrentCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
 640       fCurrentCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
 641     }
 642   }
 643
 644   printf("calling UnfoldWithBayesian\n");
 645   Int_t success = UnfoldWithBayesian(fCurrentCorrelation, fCurrentRec, initialConditions, fUnfSpectrum, regPar, nIterations, kFALSE);
 646
 647   if ( success != 0)
 648     return;
 649
 650   if (!determineError)
 651   {
 652     Printf("AliJetSpectrumUnfolding::ApplyBayesianMethod: WARNING: No errors calculated.");
 653     return;
 654   }
 655
 656   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
 657   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
 658   // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
 659   // error of the unfolding method itself.
 660
 661   TH2* maxError = (TH2*) fUnfSpectrum->Clone("maxError");
 662   maxError->Reset();
 663
 664   TH2* normalizedResult = (TH2*) fUnfSpectrum->Clone("normalizedResult");
 665   normalizedResult->Scale(1.0 / normalizedResult->Integral());
 666
 667   const Int_t kErrorIterations = 20;
 668
 669   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
 670
 671   TH2* randomized = (TH2*) fCurrentRec->Clone("randomized");
 672   TH2* result2 = (TH2*) fUnfSpectrum->Clone("result2");
 673   for (Int_t n=0; n<kErrorIterations; ++n)
 674   {
 675     // randomize the content of clone following a poisson with the mean = the value of that bin
 676     for (Int_t x=1; x<=randomized->GetNbinsX(); x++)
 677       for (Int_t y=1; y<=randomized->GetNbinsY(); y++)
 678       {
 679         Float_t randomValue = fCurrentRec->GetBinContent(x,y);
 680         TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",randomValue*0.25, randomValue*1.25);
 681         poisson->SetParameters(randomValue,0.);
 682         randomValue = poisson->GetRandom();
 683         //printf("%e --> %e\n", fCurrentRec->GetBinContent(x,y), (Double_t)randomValue);
 684         randomized->SetBinContent(x, y, randomValue);
 685         delete poisson;
 686       }
 687
 688     result2->Reset();
 689     if (UnfoldWithBayesian(fCurrentCorrelation, randomized, initialConditions, result2, regPar, nIterations) != 0)
 690       return;
 691
 692     result2->Scale(1.0 / result2->Integral());
 693
 694     // calculate ratio
 695     result2->Divide(normalizedResult);
 696
 697     //new TCanvas; result2->DrawCopy("HIST");
 698
 699     // find max. deviation
 700     for (Int_t i=1; i<=result2->GetNbinsX(); i++)
 701       for (Int_t j=1; j<=result2->GetNbinsY(); j++)
 702         maxError->SetBinContent(i, j, TMath::Max(maxError->GetBinContent(i,j), TMath::Abs(1 - result2->GetBinContent(i,j))));
 703   }
 704   delete randomized;
 705   delete result2;
 706
 707   for (Int_t i=1; i<=fUnfSpectrum->GetNbinsX(); i++)
 708     for (Int_t j=1; j<=fUnfSpectrum->GetNbinsY(); j++)
 709       fUnfSpectrum->SetBinError(i, j, fUnfSpectrum->GetBinError(i,j) + maxError->GetBinContent(i,j)*fUnfSpectrum->GetBinContent(i,j));
 710
 711   delete maxError;
 712   delete normalizedResult;
 713 }
 714
 715 //____________________________________________________________________
 716 Int_t AliJetSpectrumUnfolding::UnfoldWithBayesian(THnSparseF* correlation, TH2* measured, TH2* initialConditions, TH2* aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors)
 717 {
 718   //
 719   // unfolds a spectrum
 720   //
 721
 722   if (measured->Integral() <= 0)
 723   {
 724     Printf("AliJetSpectrumUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
 725     return 1;
 726   }
 727   const Int_t NFilledBins = correlation->GetNbins();
 728   const Int_t kStartBin = 1;
 729
 730   const Int_t kMaxTZ = NBINSZ; // max true axis fragmentation function
 731   const Int_t kMaxMZ = NBINSZ; // max measured axis fragmentation function
 732   const Int_t kMaxTE = NBINSE; // max true axis energy
 733   const Int_t kMaxME = NBINSE; // max measured axis energy
 734
 735   printf("NbinsE=%d - NbinsZ=%d\n", NBINSE, NBINSZ);
 736
 737   // store information in arrays, to increase processing speed
 738   Double_t measuredCopy[kMaxME+1][kMaxMZ+1];
 739   Double_t prior[kMaxTE+1][kMaxTZ+1];
 740   Double_t errors[kMaxTE+1][kMaxTZ+1];
 741   Double_t result[kMaxTE+1][kMaxTZ+1];
 742
 743   THnSparseF *inverseCorrelation;
 744   inverseCorrelation = (THnSparseF*)correlation->Clone("inverseCorrelation");
 745   inverseCorrelation->Reset();
 746
 747   Float_t inputDistIntegral = 1;
 748   if (initialConditions)
 749   {
 750     printf("Using different starting conditions...\n");
 751     inputDistIntegral = initialConditions->Integral();
 752   }
 753   Float_t measuredIntegral = measured->Integral();
 754   for (Int_t me=1; me<=kMaxME; me++)
 755     for (Int_t mz=1; mz<=kMaxMZ; mz++)
 756     {
 757       // normalization of the measured spectrum
 758       measuredCopy[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
 759       errors[me][mz] = measured->GetBinError(me, mz) / measuredIntegral;
 760       // pick prior distribution and normalize it
 761       if (initialConditions)
 762         prior[me][mz] = initialConditions->GetBinContent(me,mz) / inputDistIntegral;
 763       else
 764         prior[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
 765     }
 766
 767   // unfold...
 768   for (Int_t i=0; i<nIterations; i++)
 769   {
 770    // calculate Inverse Correlation Map from Bayes theorem:
 771    // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
 772    /*Float_t norm = 0;
 773    for (Int_t me=1; me<=kMaxME; me++)
 774       for (Int_t mz=1; mz<=kMaxMZ; mz++)
 775       {
 776         norm = 0;
 777         for (Int_t te=kStartBin; te<=kMaxTE; te++)
 778           for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
 779           {
 780             Int_t bin[4] = {te, me, tz, mz};
 781             norm += correlation->GetBinContent(bin)*prior[te][tz];
 782           }
 783         if (norm > 0)
 784           for (Int_t te = kStartBin; te <= kMaxTE; te++)
 785             for (Int_t tz = kStartBin; tz <= kMaxTZ; tz++)
 786             {
 787               Int_t bin[4] = {te, me, tz, mz};
 788               inverseCorrelation->SetBinContent(bin, correlation->GetBinContent(bin)*prior[te][tz]/norm );
 789             }
 790         //else
 791           // inverse response set to '0' wich has been already done in line 2069
 792       }*/
 793     inverseCorrelation->Reset();
 794     Float_t norm[kMaxTE+2][kMaxTZ+2];
 795     for (Int_t te=0; te<(kMaxTE+2); te++)
 796       for (Int_t tz=0; tz<(kMaxTZ+2); tz++)
 797         norm[te][tz]=0;
 798     for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
 799     {
 800       Int_t bin[4];
 801       Float_t binContent = correlation->GetBinContent(idx, bin);
 802       if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ &&
 803           bin[0]>0 && bin[0]<=NBINSE && bin[2]>0 && bin[2]<=NBINSZ)
 804         norm[bin[1]][bin[3]] += binContent*prior[bin[0]][bin[2]];
 805     }
 806     Float_t chi2Measured=0, diff;
 807     for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
 808     {
 809       Int_t bin[4];
 810       Float_t binContent = correlation->GetBinContent(idx, bin);
 811       if (norm[bin[1]][bin[3]]>0 && bin[1]>0 && bin[1]<=NBINSE &&
 812           bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ)
 813       {
 814         inverseCorrelation->SetBinContent(bin, binContent*prior[bin[0]][bin[2]]/norm[bin[1]][bin[3]]);
 815         if (errors[bin[1]][bin[3]]>0)
 816         {
 817           diff = ((measuredCopy[bin[1]][bin[3]]-norm[bin[1]][bin[3]])/(errors[bin[1]][bin[3]]));
 818           chi2Measured += diff*diff;
 819         }
 820       }
 821     }
 822
 823     // calculate "generated" spectrum
 824     for (Int_t te = kStartBin; te<=kMaxTE; te++)
 825       for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
 826       {
 827         Float_t value = 0;
 828         for (Int_t me=1; me<=kMaxME; me++)
 829           for (Int_t mz=1; mz<=kMaxMZ; mz++)
 830           {
 831             Int_t bin[4] = {te, me, tz, mz};
 832             value += inverseCorrelation->GetBinContent(bin)*measuredCopy[me][mz];
 833           }
 834         result[te][tz] = value;
 835         //printf("%e\n", result[te][tz]);
 836       }
 837
 838     // regularization (simple smoothing)
 839     Float_t chi2LastIter = 0;
 840     for (Int_t te=kStartBin; te<=kMaxTE; te++)
 841       for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
 842       {
 843         Float_t newValue = 0;
 844         // 0 bin excluded from smoothing
 845         if (( te >(kStartBin+1) && te<(kMaxTE-1) ) && ( tz > (kStartBin+1) && tz<(kMaxTZ-1) ))
 846         {
 847           Float_t average = ((result[te-1][tz-1] + result[te-1][tz] + result[te-1][tz+1])+(result[te][tz-1] + result[te][tz] + result[te][tz+1])+(result[te+1][tz-1] + result[te+1][tz] + result[te+1][tz+1]))/9.;
 848
 849           // weight the average with the regularization parameter
 850           newValue = (1 - regPar) * result[te][tz] + regPar * average;
 851         }
 852         else
 853           newValue = result[te][tz];
 854         if (prior[te][tz]>1.e-5)
 855         {
 856           diff = ((prior[te][tz]-newValue)/prior[te][tz]);
 857           chi2LastIter = diff*diff;
 858         }
 859         prior[te][tz] = newValue;
 860       }
 861     //printf(" iteration %d - chi2LastIter = %e - chi2Measured = %e \n", i, chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ), chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ));
 862     if (chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-6 && chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-3)
 863       break;
 864   } // end of iterations
 865
 866   // propagate errors of the reconstructed distribution through the unfolding
 867   for (Int_t te = kStartBin; te<=kMaxTE; te++)
 868       for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
 869       {
 870         Float_t valueError = 0;
 871         Float_t binError = 0;
 872         for (Int_t me=1; me<=kMaxME; me++)
 873           for (Int_t mz=1; mz<=kMaxMZ; mz++)
 874           {
 875             Int_t bin[4] = {te, me, tz, mz};
 876             valueError += inverseCorrelation->GetBinContent(bin)*inverseCorrelation->GetBinContent(bin)*errors[me][mz]*errors[me][mz];
 877           }
 878         //if (errors[te][tz]!=0)printf("errors[%d][%d]=%e\n", te, tz, valueError);
 879         aResult->SetBinContent(te, tz, prior[te][tz]);
 880         aResult->SetBinError(te, tz, TMath::Sqrt(valueError));
 881       }
 882
 883   // ***********************************************************************************************************
 884   // Calculate the covariance matrix, all arguments are taken from G. D'Agostini (p.6-8)
 885   if (calculateErrors)
 886   {
 887     printf("Covariance matrix will be calculated... this will take a lot of time (>1 day) ;)\n");
 888
 889     //Variables and Matrices that will be use along the calculation
 890     const Int_t binsV[4] = {NBINSE,NBINSE, NBINSZ, NBINSZ};
 891     const Double_t LowEdgeV[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
 892     const Double_t UpEdgeV[4] = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
 893
 894     const Double_t Ntrue = (Double_t)measured->Integral();
 895
 896     THnSparseF *V = new THnSparseF("V","",4, binsV, LowEdgeV, UpEdgeV);
 897     V->Reset();
 898     Double_t invCorrContent1, Nt;
 899     Double_t invCorrContent2, v11, v12, v2;
 900     // calculate V1 and V2
 901     for (Int_t idx1=0; idx1<=NFilledBins; idx1++)
 902     {
 903       printf("Covariance Matrix calculation: iteration idx1=%d of %d\n", idx1, NFilledBins);
 904       for (Int_t idx2=0; idx2<=NFilledBins; idx2++)
 905       {
 906         Int_t bin1[4];
 907         Int_t bin2[4];
 908         invCorrContent1 = inverseCorrelation->GetBinContent(idx1, bin1);
 909         invCorrContent2 = inverseCorrelation->GetBinContent(idx2, bin2);
 910         v11=0; v12=0; v2=0;
 911         if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE &&
 912            bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ &&
 913            bin2[0]>0 && bin2[0]<=NBINSE && bin2[1]>0 && bin2[1]<=NBINSE &&
 914            bin2[2]>0 && bin2[2]<=NBINSZ && bin2[3]>0 && bin2[3]<=NBINSZ)
 915         {
 916           if (bin1[1]==bin2[1] && bin1[3]==bin2[3])
 917             v11 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]
 918                   *(1. - measuredCopy[bin2[1]][bin2[3]]/Ntrue);
 919           else
 920             v12 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]*
 921                   measuredCopy[bin2[1]][bin2[3]]/Ntrue;
 922           Nt = (Double_t)prior[bin2[0]][bin2[2]];
 923           v2 = measuredCopy[bin1[1]][bin1[3]]*measuredCopy[bin2[1]][bin2[3]]*
 924                invCorrContent1*invCorrContent2*
 925                BayesUncertaintyTerms(inverseCorrelation, correlation, bin1, bin2, Nt);
 926           Int_t binV[4] = {bin1[0],bin2[0],bin1[2],bin2[2]};
 927           V->SetBinContent(binV,v11-v12 + v2);
 928         }
 929       }
 930     }
 931
 932     for(Int_t te = 1; te<=NBINSE; te++)
 933       for(Int_t tz = 1; tz<=NBINSZ; tz++)
 934       {
 935         Int_t binV[4] = {te,te,tz,tz};
 936         aResult->SetBinError( te, tz, V->GetBinContent(binV) );
 937       }
 938
 939     TFile* f = new TFile("Covariance_UnfSpectrum.root");
 940     f->Open("RECREATE");
 941     V->Write();
 942     f->Close();
 943   }
 944
 945   return 0;
 946
 947 }
 948
 949 //____________________________________________________________________
 950 Double_t AliJetSpectrumUnfolding::BayesUncertaintyTerms(THnSparseF *M, THnSparseF *C, Int_t* binTM, Int_t* binTM1, Double_t Nt)
 951 {
 952   //
 953   // helper function for the covariance matrix of the bayesian method
 954   //
 955
 956   Double_t result = 0;
 957   Float_t term[9];
 958   Int_t tmpBin[4], tmpBin1[4];
 959   const Int_t nFilledBins = C->GetNbins();
 960   if (Nt==0)
 961     return 0;
 962
 963   Float_t CorrContent;
 964   Float_t InvCorrContent;
 965
 966   tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1];  tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
 967   tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
 968   if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
 969   {
 970     if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
 971       term[0] = BayesCov(M, C, tmpBin, tmpBin1)/
 972                 (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
 973     term[2] = term[0]*M->GetBinContent(tmpBin1);
 974   }
 975   else
 976   {
 977     term[0] = 0;
 978     term[2] = 0;
 979   }
 980
 981   tmpBin[0]=binTM1[0]; tmpBin[1]=binTM[1]; tmpBin[2]=binTM1[2]; tmpBin[3]=binTM[3];
 982   tmpBin1[0]=binTM1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=binTM1[3];
 983   if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
 984     term[6] = BayesCov(M, C, tmpBin, tmpBin1)*
 985               M->GetBinContent(tmpBin)/
 986               (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
 987   else
 988     term[6] = 0;
 989
 990   for(Int_t idx1=0; idx1<=nFilledBins; idx1++)
 991   {
 992     Int_t bin1[4];
 993     CorrContent    = C->GetBinContent(idx1, bin1);
 994     InvCorrContent = M->GetBinContent(idx1, bin1);
 995     if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE &&
 996        bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ)
 997     {
 998       tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
 999       tmpBin1[0]=binTM[0]; tmpBin1[1]=bin1[1];  tmpBin1[2]=binTM[2]; tmpBin1[3]=bin1[3];
1000       if (C->GetBinContent(tmpBin)!=0 &&
1001           binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
1002         term[1] = BayesCov(M, C, tmpBin, tmpBin1)/C->GetBinContent(tmpBin);
1003       else
1004         term[1] = 0;
1005
1006       tmpBin[0] =binTM[0]; tmpBin[1] =bin1[1];   tmpBin[2] =binTM[2]; tmpBin[3] =bin1[3];
1007       tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
1008       if (C->GetBinContent(tmpBin1)!=0)
1009       {
1010         if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
1011           term[3] = BayesCov(M, C, tmpBin, tmpBin1)/
1012                     C->GetBinContent(tmpBin1);
1013         term[5] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin1)/
1014                   C->GetBinContent(tmpBin1);
1015       }
1016       else
1017       {
1018         term[3] = 0;
1019         term[5] = 0;
1020       }
1021
1022       tmpBin[0] =binTM1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM1[2]; tmpBin[3] =binTM[3];
1023       tmpBin1[0]=binTM1[0]; tmpBin1[1]=bin1[1];  tmpBin1[2]=binTM1[2]; tmpBin1[3]=bin1[3];
1024       if (C->GetBinContent(tmpBin)!=0)
1025         term[7] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin)/
1026                   C->GetBinContent(tmpBin);
1027       else
1028         term[7] = 0;
1029
1030       tmpBin[0] =bin1[0]; tmpBin[1] =binTM[1];  tmpBin[2] =bin1[2]; tmpBin[3] =binTM[3];
1031       tmpBin1[0]=bin1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=bin1[2]; tmpBin1[3]=binTM1[3];
1032       if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
1033         term[8] = BayesCov(M, C, tmpBin, tmpBin1)*
1034                   M->GetBinContent(tmpBin)*M->GetBinContent(tmpBin)/
1035                   (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
1036       else
1037         term[8] = 0;
1038
1039       for (Int_t i=0; i<9; i++)
1040         result += term[i]/Nt;
1041     }
1042   }
1043
1044   return result;
1045 }
1046
1047 //____________________________________________________________________
1048 Double_t AliJetSpectrumUnfolding::BayesCov(THnSparseF *M, THnSparseF *correlation, Int_t* binTM, Int_t* bin1)
1049 {
1050   Double_t result, result1, result2, result3;
1051
1052   if (binTM[0]==bin1[0] && binTM[2]==bin1[2])
1053   {
1054     if (correlation->GetBinContent(bin1)!=0)
1055       result1 = 1./correlation->GetBinContent(bin1);
1056     else
1057       result1 = 0;
1058     result2 = 1.;
1059   }
1060   else
1061   {
1062     result1 = 0;
1063     result2 = 0;
1064   }
1065
1066   if (binTM[1]==bin1[1] && binTM[3]==bin1[3])
1067   {
1068     Int_t tmpbin[4] = {bin1[0], binTM[1], bin1[2], binTM[3]};
1069     if(correlation->GetBinContent(tmpbin)!=0)
1070       result3 = M->GetBinContent(tmpbin)/correlation->GetBinContent(tmpbin);
1071     else
1072       result3 = 0;
1073   }
1074   else
1075   {
1076     result1 = 0;
1077     result3 = 0;
1078   }
1079
1080   return result = result1 + result2 + result3;
1081 }
1082
1083 //____________________________________________________________________
1084 TH2F* AliJetSpectrumUnfolding::CalculateRecSpectrum(TH2* inputGen)
1085 {
1086   // runs the distribution given in inputGen through the correlation histogram identified by
1087   // fCorrelation and produces a reconstructed spectrum
1088
1089   if (!inputGen)
1090     return 0;
1091
1092   // normalize to convert number of events into probability
1093   /*for (Int_t te=1; te<=NBINSE; te++)
1094     for (Int_t tz=1; tz<=NBINSZ; tz++)
1095     {
1096        Int_t bin[4];
1097        Float_t sum=0.;
1098        for (Int_t me = 1; me<=NBINSE; me++)
1099          for (Int_t mz = 1; mz<=NBINSZ; mz++)
1100          {
1101            bin[0] = te; bin[1] = me;
1102            bin[2] = tz; bin[3] = mz;
1103            sum += fCorrelation[correlationMap]->GetBinContent(bin);
1104          }
1105        if (sum > 0.)
1106          for (Int_t me = 1; me<=NBINSE; me++)
1107            for (Int_t mz = 1; mz<=NBINSZ; mz++)
1108            {
1109              bin[0] = te; bin[1] = me;
1110              bin[2] = tz; bin[3] = mz;
1111              fCorrelation[correlationMap]->SetBinContent(bin, fCorrelation[correlationMap]->GetBinContent(bin)/sum);
1112              fCorrelation[correlationMap]->SetBinError(bin, fCorrelation[correlationMap]->GetBinError(bin)/sum);
1113            }
1114     }*/
1115   // normalize to convert number of events into probability (the following loop is much faster)
1116   Float_t sum[NBINSE+2][NBINSZ+2];
1117   memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
1118
1119   for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
1120   {
1121     Int_t bin[4];
1122     Float_t binContent = fCorrelation->GetBinContent(idx, bin);
1123     if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ){
1124       sum[bin[0]][bin[2]] += binContent;
1125     }
1126   }
1127
1128   for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
1129   {
1130     Int_t bin[4];
1131     Float_t binContent = fCorrelation->GetBinContent(idx, bin);
1132     Float_t binError   = fCorrelation->GetBinError(bin);
1133     if (sum[bin[0]][bin[2]]>0 && bin[1]>0 && bin[1]<=NBINSE &&
1134         bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ)
1135     {
1136       fCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
1137       fCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
1138     }
1139   }
1140
1141   TH2F* target = dynamic_cast<TH2F*> (fRecSpectrum->Clone(Form("reconstructed_%s", inputGen->GetName())));
1142   target->Reset();
1143
1144   for (Int_t me=1; me<=NBINSE; ++me)
1145     for (Int_t mz=1; mz<=NBINSZ; ++mz)
1146     {
1147       Float_t measured = 0;
1148       Float_t error = 0;
1149
1150       for (Int_t te=1; te<=NBINSE; ++te)
1151         for (Int_t tz=1; tz<=NBINSZ; ++tz)
1152         {
1153           Int_t bin[4] = {te, me, tz, mz};
1154           measured += inputGen->GetBinContent(te,tz) * fCorrelation->GetBinContent(bin);
1155           error += inputGen->GetBinError(te,tz) * fCorrelation->GetBinContent(bin);
1156         }
1157       target->SetBinContent(me, mz, measured);
1158       target->SetBinError(me, mz, error);
1159     }
1160
1161   return target;
1162 }
1163
1164 //__________________________________________________________________________________________________
1165 void AliJetSpectrumUnfolding::SetGenRecFromFunc(TF2* inputGen)
1166 {
1167   // uses the given function to fill the input Generated histogram and generates from that
1168   // the reconstructed histogram by applying the response histogram
1169   // this can be used to evaluate if the methods work indepedently of the input
1170   // distribution
1171
1172   if (!inputGen)
1173     return;
1174
1175   TH2F* histtmp = new TH2F("histtmp", "tmp", EBINNING, ZBINNING);
1176   TH2F* gen  = fGenSpectrum;
1177
1178   histtmp->Reset();
1179   gen->Reset();
1180
1181   histtmp->FillRandom(inputGen->GetName(), NEVENTS);
1182
1183   for (Int_t i=1; i<=gen->GetNbinsX(); ++i)
1184     for (Int_t j=1; j<=gen->GetNbinsY(); ++j)
1185     {
1186       gen->SetBinContent(i, j, histtmp->GetBinContent(i,j));
1187       gen->SetBinError(i, j, histtmp->GetBinError(i,j));
1188     }
1189
1190   delete histtmp;
1191
1192   //new TCanvas;
1193   //gStyle->SetPalette(1);
1194   //gPad->SetLogz();
1195   //gen->Draw("COLZ");
1196
1197
1198   TH2 *recsave = fRecSpectrum;
1199
1200   fRecSpectrum = CalculateRecSpectrum(gen);
1201   fRecSpectrum->SetName(recsave->GetName());
1202   delete recsave;
1203
1204   return;
1205 }
1206 //________________________________________________________________________________________