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