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