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