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