--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+//
+// This class is used to store correlation maps, generated and reconstructed data of the jet spectrum
+// It also contains functions to correct the spectrum using the bayesian unfolding
+//
+
+#include "AliJetSpectrumUnfolding.h"
+
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TDirectory.h>
+#include <TVirtualFitter.h>
+#include <TCanvas.h>
+#include <TString.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TMath.h>
+#include <TCollection.h>
+#include <TLegend.h>
+#include <TLine.h>
+#include <TRandom.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+#include <TStyle.h>
+#include <TColor.h>
+
+#include <THnSparse.h>
+
+ClassImp(AliJetSpectrumUnfolding)
+
+const Int_t NBINSE = 50;
+const Int_t NBINSZ = 50;
+const Int_t NEVENTS = 500000;
+const Double_t axisLowerLimitE = 0.;
+const Double_t axisLowerLimitZ = 0.;
+const Double_t axisUpperLimitE = 250.;
+const Double_t axisUpperLimitZ = 1.;
+
+Float_t AliJetSpectrumUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
+Int_t AliJetSpectrumUnfolding::fgBayesianIterations = 100; // number of iterations in Bayesian method
+
+//____________________________________________________________________
+
+AliJetSpectrumUnfolding::AliJetSpectrumUnfolding() :
+ TNamed(), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
+ fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
+{
+ //
+ // default constructor
+ //
+
+ fGenSpectrum = 0;
+ fRecSpectrum = 0;
+ fUnfSpectrum = 0;
+ fCorrelation = 0;
+}
+
+//____________________________________________________________________
+AliJetSpectrumUnfolding::AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title) :
+ TNamed(name, title), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
+ fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
+{
+ //
+ // named constructor
+ //
+
+ // do not add this hists to the directory
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ #define ZBINNING NBINSZ, axisLowerLimitZ, axisUpperLimitZ // fragmentation of leading particle
+ #define EBINNING NBINSE, axisLowerLimitE, axisUpperLimitE // energy of the jet
+
+ fRecSpectrum = new TH2F("fRecSpectrum", "Reconstructed Spectrum;E^{jet}_{rec} [GeV];z^{lp}_{rec}", EBINNING, ZBINNING);
+ fGenSpectrum = new TH2F("fGenSpectrum", "Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}_{gen}", EBINNING, ZBINNING);
+ fUnfSpectrum = new TH2F("fUnfSpectrum", "Unfolded Spectrum;E^{jet} [GeV];z^{lp}", EBINNING, ZBINNING);
+
+ const Int_t nbin[4]={NBINSE, NBINSE, NBINSZ, NBINSZ};
+ //arrays for bin limits
+ Double_t LowEdge[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
+ Double_t UpEdge[4] = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
+
+ fCorrelation = new THnSparseF("fCorrelation", "Correlation Function", 4, nbin, LowEdge, UpEdge);
+
+ TH1::AddDirectory(oldStatus);
+}
+
+//____________________________________________________________________
+AliJetSpectrumUnfolding::~AliJetSpectrumUnfolding()
+{
+ //
+ // Destructor
+ //
+
+ if (fGenSpectrum)
+ delete fGenSpectrum;
+ fGenSpectrum = 0;
+
+ if (fRecSpectrum)
+ delete fRecSpectrum;
+ fRecSpectrum = 0;
+
+ if (fUnfSpectrum)
+ delete fUnfSpectrum;
+ fUnfSpectrum = 0;
+
+ if (fCorrelation)
+ delete fCorrelation;
+ fCorrelation = 0;
+
+}
+
+//____________________________________________________________________
+Long64_t AliJetSpectrumUnfolding::Merge(TCollection* list)
+{
+ // Merge a list of AliJetSpectrumUnfolding objects with this (needed for
+ // PROOF).
+ // Returns the number of merged objects (including this).
+
+ if (!list)
+ return 0;
+
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // collections of all histograms
+ TList collections[4];
+
+ Int_t count = 0;
+ while ((obj = iter->Next())) {
+
+ AliJetSpectrumUnfolding* entry = dynamic_cast<AliJetSpectrumUnfolding*> (obj);
+ if (entry == 0)
+ continue;
+
+ collections[0].Add(entry->fGenSpectrum);
+ collections[1].Add(entry->fRecSpectrum);
+ collections[2].Add(entry->fUnfSpectrum);
+ collections[3].Add(entry->fCorrelation);
+
+ count++;
+ }
+
+ fGenSpectrum->Merge(&collections[0]);
+ fRecSpectrum->Merge(&collections[1]);
+ fUnfSpectrum->Merge(&collections[2]);
+ fCorrelation->Merge(&collections[3]);
+
+ delete iter;
+
+ return count+1;
+}
+
+//____________________________________________________________________
+Bool_t AliJetSpectrumUnfolding::LoadHistograms(const Char_t* dir)
+{
+ //
+ // loads the histograms from a file
+ // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
+ //
+
+ if (!dir)
+ dir = GetName();
+
+ if (!gDirectory->cd(dir))
+ return kFALSE;
+
+ Bool_t success = kTRUE;
+
+ // store old histograms to delete them later
+ TList oldHistograms;
+ oldHistograms.SetOwner(1);
+
+ if (fGenSpectrum) oldHistograms.Add(fGenSpectrum);
+ if (fRecSpectrum) oldHistograms.Add(fRecSpectrum);
+ if (fUnfSpectrum) oldHistograms.Add(fUnfSpectrum);
+ if (fCorrelation) oldHistograms.Add(fCorrelation);
+
+ // load new histograms
+ fGenSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fGenSpectrum->GetName()));
+ if (!fGenSpectrum)
+ success = kFALSE;
+
+ fRecSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fRecSpectrum->GetName()));
+ if (!fRecSpectrum)
+ success = kFALSE;
+
+ fUnfSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fUnfSpectrum->GetName()));
+ if (!fUnfSpectrum)
+ success = kFALSE;
+
+ fCorrelation = dynamic_cast<THnSparseF*> (gDirectory->Get(fCorrelation->GetName()));
+ if (!fCorrelation)
+ success = kFALSE;
+
+ gDirectory->cd("..");
+
+ // delete old histograms
+ oldHistograms.Delete();
+
+ return success;
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::SaveHistograms()
+{
+ //
+ // saves the histograms
+ //
+
+ gDirectory->mkdir(GetName());
+ gDirectory->cd(GetName());
+
+ if (fGenSpectrum)
+ fGenSpectrum->Write();
+
+ if (fRecSpectrum)
+ fRecSpectrum->Write();
+
+ if (fUnfSpectrum)
+ fUnfSpectrum->Write();
+
+ if (fCorrelation)
+ fCorrelation->Write();
+
+ gDirectory->cd("..");
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::SetupCurrentHists(Bool_t createBigBin)
+{
+ //
+ // resets fUnfSpectrum
+ //
+
+ fUnfSpectrum->Reset();
+ fUnfSpectrum->Sumw2();
+
+ fCurrentRec = (TH2F*)fRecSpectrum->Clone("fCurrentRec");
+ fCurrentRec->Sumw2();
+
+ fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation");
+ fCurrentCorrelation->Sumw2();
+
+ if (createBigBin)
+ {
+ Int_t maxBinE = 0, maxBinZ = 0;
+ Float_t maxE = 0, maxZ = 0;
+ for (Int_t me=1; me<=fCurrentRec->GetNbinsX(); me++)
+ for (Int_t mz=1; mz<=fCurrentRec->GetNbinsY(); mz++)
+ {
+ if (fCurrentRec->GetBinContent(me,mz) <= 5 && me>NBINSE/2 && mz>NBINSZ/2)
+ {
+ maxBinE = me;
+ maxBinZ = mz;
+ maxE = fCurrentRec->GetXaxis()->GetBinCenter(me);
+ maxZ = fCurrentRec->GetYaxis()->GetBinCenter(mz);
+ break;
+ }
+ }
+
+ if (maxBinE > 0 || maxBinZ > 0)
+ {
+ printf("Bin limit in measured spectrum is e = %d and z = %d.\n", maxBinE, maxBinZ);
+ fCurrentRec->SetBinContent(maxBinE, maxBinZ, fCurrentRec->Integral(maxBinE, fCurrentRec->GetNbinsX(), maxBinZ, fCurrentRec->GetNbinsY()));
+ for (Int_t me=maxBinE+1; me<=fCurrentRec->GetNbinsX(); me++)
+ for (Int_t mz=maxBinZ+1; mz<=fCurrentRec->GetNbinsY(); mz++)
+ {
+ fCurrentRec->SetBinContent(me, mz, 0);
+ fCurrentRec->SetBinError(me, mz, 0);
+ }
+ // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+ fCurrentRec->SetBinError(maxBinE, maxBinZ, TMath::Sqrt(fCurrentRec->GetBinContent(maxBinE, maxBinZ)));
+
+ printf("This bin has now %f +- %f entries\n", fCurrentRec->GetBinContent(maxBinE, maxBinZ), fCurrentRec->GetBinError(maxBinE, maxBinZ));
+
+ /* for (Int_t te=1; te<=NBINSE; te++)
+ {
+ for (Int_t tz=1; tz<=NBINSZ; tz++)
+ {
+ Int_t binMin[4] = {te, maxBinE, tz, maxBinZ};
+ Int_t binMax[4] = {NBINSE, NBINSE, NBINSZ, NBINSZ};
+ Float_t sum=0;
+ for (Int_t ite=te; ite<=NBINSE; ite++)
+ for (Int_t itz=tz; itz<=NBINSZ; itz++)
+ for (Int_t ime=maxBinE; ime<=NBINSE; ime++)
+ for (Int_t imz=maxBinZ; imz<=NBINSZ; imz++)
+ {
+ Int_t bin[4] = {ite, ime, itz, imz};
+ sum += fCurrentCorrelation->GetBinContent(bin);
+ }
+ fCurrentCorrelation->SetBinContent(binMin, sum);
+ fCurrentCorrelation->SetBinError(binMin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
+ printf("create big bin1, nbins = %d, te = %d, tz = %d\n", NBINSE, te, tz);
+ for (Int_t me=maxBinE; me<=NBINSE; me++)
+ {
+ for (Int_t mz=maxBinZ; mz<=NBINSZ; mz++)
+ {
+ Int_t bin[4] = {te, me, tz, mz};
+ fCurrentCorrelation->SetBinContent(bin, 0.);
+ fCurrentCorrelation->SetBinError(bin, 0.);
+ printf("create big bin2\n");
+ }
+ }
+ }
+ }*/
+
+ for(Int_t idx = 0; idx<=fCurrentCorrelation->GetNbins(); idx++)
+ {
+ Int_t bin[4];
+ Float_t binContent = fCurrentCorrelation->GetBinContent(idx,bin);
+ Float_t binError = fCurrentCorrelation->GetBinError(idx);
+ Int_t binMin[4] = {bin[0], maxBinE, bin[2], maxBinZ};
+ if ( (bin[1]>maxBinE && bin[1]<=NBINSE) && (bin[3]>maxBinZ && bin[3]<=NBINSZ) )
+ {
+ fCurrentCorrelation->SetBinContent(binMin, binContent + fCurrentCorrelation->GetBinContent(binMin));
+ fCurrentCorrelation->SetBinError(binMin, binError + TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
+ fCurrentCorrelation->SetBinContent(bin, 0.);
+ fCurrentCorrelation->SetBinError(bin, 0.);
+ }
+ printf("create big bin1, nbins = %d, te = %d, tz = %d\n", NBINSE, bin[0], bin[1]);
+ }
+
+ printf("Adjusted correlation matrix!\n");
+ }
+ } // end Create Big Bin
+
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
+{
+ //
+ // sets the parameters for Bayesian unfolding
+ //
+
+ fgBayesianSmoothing = smoothing;
+ fgBayesianIterations = nIterations;
+
+ printf("AliJetSpectrumUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::NormalizeToBinWidth(TH2* hist)
+{
+ //
+ // normalizes a 2-d histogram to its bin width (x width * y width)
+ //
+
+ for (Int_t i=1; i<=hist->GetNbinsX(); i++)
+ for (Int_t j=1; j<=hist->GetNbinsY(); j++)
+ {
+ Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
+ hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
+ hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
+ }
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::DrawHistograms()
+{
+ //
+ // draws the histograms of this class
+ //
+
+ gStyle->SetPalette(1);
+
+ TCanvas* canvas1 = new TCanvas("fRecSpectrum", "fRecSpectrum", 900, 600);
+ gPad->SetLogz();
+ fRecSpectrum->DrawCopy("COLZ");
+
+ TCanvas* canvas2 = new TCanvas("fGenSpectrum", "fGenSpectrum", 900, 600);
+ gPad->SetLogz();
+ fGenSpectrum->DrawCopy("COLZ");
+
+ TCanvas* canvas3 = new TCanvas("fUnfSpectrum", "fUnfSpectrum", 900, 600);
+ gPad->SetLogz();
+ fUnfSpectrum->DrawCopy("COLZ");
+
+ TCanvas* canvas4 = new TCanvas("fCorrelation", "fCorrelation", 500, 500);
+ canvas1->Divide(2);
+
+ canvas4->cd(1);
+ gPad->SetLogz();
+ TH2D* h0 = fCorrelation->Projection(1,0);
+ h0->SetXTitle("E^{jet}_{gen} [GeV]");
+ h0->SetYTitle("E^{jet}_{rec} [GeV]");
+ h0->SetTitle("Projection: Jet Energy");
+ h0->DrawCopy("colz");
+
+ canvas1->cd(2);
+ gPad->SetLogz();
+ TH2D* h1 = fCorrelation->Projection(3,2);
+ h1->SetXTitle("z^{lp}_{gen}");
+ h1->SetYTitle("z^{lp}_{rec}");
+ h1->SetTitle("Projection: Leading Particle Fragmentation");
+ h1->DrawCopy("colz");
+
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::DrawComparison(const char* name, TH2* genHist)
+{
+
+ if (fUnfSpectrum->Integral() == 0)
+ {
+ printf("ERROR. Unfolded histogram is empty\n");
+ return;
+ }
+
+ //regain measured distribution used for unfolding, because the bins were modified in SetupCurrentHists
+ //in create big bin
+ fCurrentRec = (TH2F*)fRecSpectrum->Clone();
+ fCurrentRec->Sumw2();
+ fCurrentRec->Scale(1.0 / fCurrentRec->Integral());
+
+ // normalize unfolded result to 1
+ fUnfSpectrum->Scale(1.0 / fUnfSpectrum->Integral());
+
+ // find bin with <= 100 entries. this is used as limit for the chi2 calculation
+ Int_t mcBinLimitE = 0, mcBinLimitZ = 0;
+ for (Int_t i=0; i<genHist->GetNbinsX(); ++i)
+ for (Int_t j=0; j<genHist->GetNbinsY(); ++j)
+ {
+ if (genHist->GetBinContent(i,j) > 100)
+ {
+ mcBinLimitE = i;
+ mcBinLimitZ = j;
+ }
+ else
+ break;
+ }
+ Printf("AliJetSpectrumUnfolding::DrawComparison: Gen bin limit is (x,y) = (%d,%d)", mcBinLimitE,mcBinLimitZ);
+
+ // scale to 1 true spectrum
+ genHist->Sumw2();
+ genHist->Scale(1.0 / genHist->Integral());
+
+ // calculate residual
+ // for that we convolute the response matrix with the gathered result
+ TH2* tmpRecRecorrected = (TH2*) fUnfSpectrum->Clone("tmpRecRecorrected");
+ TH2* convoluted = CalculateRecSpectrum(tmpRecRecorrected);
+ if (convoluted->Integral() > 0)
+ convoluted->Scale(1.0 / convoluted->Integral());
+ else
+ printf("ERROR: convoluted is empty. Something went wrong calculating the convoluted histogram.\n");
+
+ TH2* residual = (TH2*) convoluted->Clone("residual");
+ residual->SetTitle("(R#otimesUnfolded - Reconstructed)/Reconstructed;E^{jet} [GeV]; z^{lp}");
+
+ fCurrentRec->Scale(1./fCurrentRec->Integral());
+ residual->Add(fCurrentRec, -1);
+ //residual->Divide(residual, fCurrentRec, 1, 1, "B");
+
+ // draw canvas
+ TCanvas* canvas1 = new TCanvas(name, name, 1000, 1000);
+ canvas1->Divide(2, 3);
+
+ Int_t style = 1;
+ const Int_t NRGBs = 5;
+ const Int_t NCont = 500;
+
+ Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+ Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+ Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+ Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+ TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+ gStyle->SetNumberContours(NCont);
+
+ canvas1->cd(1);
+ gStyle->SetPalette(style);
+ gPad->SetLogz();
+ genHist->SetTitle("Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}");
+ genHist->SetStats(0);
+ genHist->DrawCopy("colz");
+
+ canvas1->cd(2);
+ gStyle->SetPalette(style);
+ gPad->SetLogz();
+ fUnfSpectrum->SetStats(0);
+ fUnfSpectrum->DrawCopy("colz");
+
+ canvas1->cd(3);
+ gStyle->SetPalette(style);
+ gPad->SetLogz();
+ fCurrentRec->SetTitle(fRecSpectrum->GetTitle());
+ fCurrentRec->SetStats(0);
+ fCurrentRec->DrawCopy("colz");
+
+ canvas1->cd(4);
+ gStyle->SetPalette(style);
+ gPad->SetLogy();
+ TH1D* projGenX = genHist->ProjectionX();
+ projGenX->SetTitle("Projection: Jet Energy; E^{jet} [GeV]");
+ TH1D* projUnfX = fUnfSpectrum->ProjectionX();
+ TH1D* projRecX = fCurrentRec->ProjectionX();
+ projGenX->SetStats(0);
+ projRecX->SetStats(0);
+ projUnfX->SetStats(0);
+ projGenX->SetLineColor(8);
+ projRecX->SetLineColor(2);
+ projGenX->DrawCopy();
+ projUnfX->DrawCopy("same");
+ projRecX->DrawCopy("same");
+
+ TLegend* legend = new TLegend(0.6, 0.85, 0.98, 0.98);
+ legend->AddEntry(projGenX, "Generated Spectrum");
+ legend->AddEntry(projUnfX, "Unfolded Spectrum");
+ legend->AddEntry(projRecX, "Reconstructed Spectrum");
+ //legend->SetFillColor(0);
+ legend->Draw("same");
+
+ canvas1->cd(5);
+ gPad->SetLogy();
+ gStyle->SetPalette(style);
+ TH1D* projGenY = genHist->ProjectionY();
+ projGenY->SetTitle("Projection: Leading Particle Fragmentation; z^{lp}");
+ TH1D* projUnfY = fUnfSpectrum->ProjectionY();
+ TH1D* projRecY = fCurrentRec->ProjectionY();
+ projGenY->SetStats(0);
+ projRecY->SetStats(0);
+ projUnfY->SetStats(0);
+ projGenY->SetLineColor(8);
+ projRecY->SetLineColor(2);
+ projGenY->DrawCopy();
+ projUnfY->DrawCopy("same");
+ projRecY->DrawCopy("same");
+
+ TLegend* legend1 = new TLegend(0.6, 0.85, 0.98, 0.98);
+ legend1->AddEntry(projGenY, "Generated Spectrum");
+ legend1->AddEntry(projUnfY, "Unfolded Spectrum");
+ legend1->AddEntry(projRecY, "Recontructed Spectrum");
+ //legend1->SetFillColor(0);
+ legend1->Draw("same");
+
+ // Draw residuals
+ canvas1->cd(6);
+ gStyle->SetPalette(style);
+ gPad->SetLogz();
+ residual->SetStats(0);
+ residual->DrawCopy("colz");
+
+ canvas1->SaveAs(Form("%s.png", canvas1->GetName()));
+}
+
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::GetComparisonResults(Float_t* gen, Int_t* genLimit, Float_t* residuals, Float_t* ratioAverage) const
+{
+ // Returns the chi2 between the Generated and the unfolded Reconstructed spectrum as well as between the Reconstructed and the folded unfolded
+ // These values are computed during DrawComparison, thus this function picks up the
+ // last calculation
+
+ if (gen)
+ *gen = fLastChi2MC;
+ if (genLimit)
+ *genLimit = fLastChi2MCLimit;
+ if (residuals)
+ *residuals = fLastChi2Residuals;
+ if (ratioAverage)
+ *ratioAverage = fRatioAverage;
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::ApplyBayesianMethod(Float_t regPar, Int_t nIterations, TH2* initialConditions, Bool_t determineError)
+{
+ //
+ // correct spectrum using bayesian unfolding
+ //
+
+ // initialize seed with current time
+ gRandom->SetSeed(0);
+
+ printf("seting up current arrays and histograms...\n");
+ SetupCurrentHists(kFALSE); // kFALSE to not create big bin
+
+ // normalize Correlation Map to convert number of events into probabilities
+ /*for (Int_t te=1; te<=NBINSE; te++)
+ for (Int_t tz=1; tz<=NBINSZ; tz++)
+ {
+ Int_t bin[4];
+ Float_t sum=0.;
+ for (Int_t me = 1; me<=NBINSE; me++)
+ for (Int_t mz = 1; mz<=NBINSZ; mz++)
+ {
+ bin[0] = te; bin[1] = me;
+ bin[2] = tz; bin[3] = mz;
+ sum += fCurrentCorrelation->GetBinContent(bin);
+ }
+ if (sum > 0.)
+ for (Int_t me = 1; me<=NBINSE; me++)
+ for (Int_t mz = 1; mz<=NBINSZ; mz++)
+ {
+ bin[0] = te; bin[1] = me;
+ bin[2] = tz; bin[3] = mz;
+ fCurrentCorrelation->SetBinContent(bin, fCurrentCorrelation->GetBinContent(bin)/sum);
+ fCurrentCorrelation->SetBinError(bin, fCurrentCorrelation->GetBinError(bin)/sum);
+ }
+ }*/
+ Float_t sum[NBINSE+2][NBINSZ+2];
+ memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
+
+ for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
+ {
+ Int_t bin[4];
+ Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
+ if ( (bin[1]>0 && bin[1]<=NBINSE) && (bin[3]>0 && bin[3]<=NBINSZ) )
+ sum[bin[0]][bin[2]] += binContent;
+ }
+
+ for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
+ {
+ Int_t bin[4];
+ Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
+ Float_t binError = fCurrentCorrelation->GetBinError(bin);
+ if (sum[bin[0]][bin[2]]>0 && (bin[1]>0 && bin[1]<=NBINSE) &&
+ (bin[3]>0 && bin[3]<=NBINSZ) && (bin[0]>0 && bin[0]<=NBINSE) && (bin[2]>0 && bin[2]<=NBINSZ) )
+ {
+ fCurrentCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
+ fCurrentCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
+ }
+ }
+
+ printf("calling UnfoldWithBayesian\n");
+ Int_t success = UnfoldWithBayesian(fCurrentCorrelation, fCurrentRec, initialConditions, fUnfSpectrum, regPar, nIterations, kFALSE);
+
+ if ( success != 0)
+ return;
+
+ if (!determineError)
+ {
+ Printf("AliJetSpectrumUnfolding::ApplyBayesianMethod: WARNING: No errors calculated.");
+ return;
+ }
+
+ // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
+ // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
+ // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
+ // error of the unfolding method itself.
+
+ TH2* maxError = (TH2*) fUnfSpectrum->Clone("maxError");
+ maxError->Reset();
+
+ TH2* normalizedResult = (TH2*) fUnfSpectrum->Clone("normalizedResult");
+ normalizedResult->Scale(1.0 / normalizedResult->Integral());
+
+ const Int_t kErrorIterations = 20;
+
+ printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
+
+ TH2* randomized = (TH2*) fCurrentRec->Clone("randomized");
+ TH2* result2 = (TH2*) fUnfSpectrum->Clone("result2");
+ for (Int_t n=0; n<kErrorIterations; ++n)
+ {
+ // randomize the content of clone following a poisson with the mean = the value of that bin
+ for (Int_t x=1; x<=randomized->GetNbinsX(); x++)
+ for (Int_t y=1; y<=randomized->GetNbinsY(); y++)
+ {
+ Float_t randomValue = fCurrentRec->GetBinContent(x,y);
+ TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",randomValue*0.25, randomValue*1.25);
+ poisson->SetParameters(randomValue,0.);
+ randomValue = poisson->GetRandom();
+ //printf("%e --> %e\n", fCurrentRec->GetBinContent(x,y), (Double_t)randomValue);
+ randomized->SetBinContent(x, y, randomValue);
+ delete poisson;
+ }
+
+ result2->Reset();
+ if (UnfoldWithBayesian(fCurrentCorrelation, randomized, initialConditions, result2, regPar, nIterations) != 0)
+ return;
+
+ result2->Scale(1.0 / result2->Integral());
+
+ // calculate ratio
+ result2->Divide(normalizedResult);
+
+ //new TCanvas; result2->DrawCopy("HIST");
+
+ // find max. deviation
+ for (Int_t i=1; i<=result2->GetNbinsX(); i++)
+ for (Int_t j=1; j<=result2->GetNbinsY(); j++)
+ maxError->SetBinContent(i, j, TMath::Max(maxError->GetBinContent(i,j), TMath::Abs(1 - result2->GetBinContent(i,j))));
+ }
+ delete randomized;
+ delete result2;
+
+ for (Int_t i=1; i<=fUnfSpectrum->GetNbinsX(); i++)
+ for (Int_t j=1; j<=fUnfSpectrum->GetNbinsY(); j++)
+ fUnfSpectrum->SetBinError(i, j, fUnfSpectrum->GetBinError(i,j) + maxError->GetBinContent(i,j)*fUnfSpectrum->GetBinContent(i,j));
+
+ delete maxError;
+ delete normalizedResult;
+}
+
+//____________________________________________________________________
+Int_t AliJetSpectrumUnfolding::UnfoldWithBayesian(THnSparseF* correlation, TH2* measured, TH2* initialConditions, TH2* aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors)
+{
+ //
+ // unfolds a spectrum
+ //
+
+ if (measured->Integral() <= 0)
+ {
+ Printf("AliJetSpectrumUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
+ return 1;
+ }
+ const Int_t NFilledBins = correlation->GetNbins();
+ const Int_t kStartBin = 1;
+
+ const Int_t kMaxTZ = NBINSZ; // max true axis fragmentation function
+ const Int_t kMaxMZ = NBINSZ; // max measured axis fragmentation function
+ const Int_t kMaxTE = NBINSE; // max true axis energy
+ const Int_t kMaxME = NBINSE; // max measured axis energy
+
+ printf("NbinsE=%d - NbinsZ=%d\n", NBINSE, NBINSZ);
+
+ // store information in arrays, to increase processing speed
+ Double_t measuredCopy[kMaxME+1][kMaxMZ+1];
+ Double_t prior[kMaxTE+1][kMaxTZ+1];
+ Double_t errors[kMaxTE+1][kMaxTZ+1];
+ Double_t result[kMaxTE+1][kMaxTZ+1];
+
+ THnSparseF *inverseCorrelation;
+ inverseCorrelation = (THnSparseF*)correlation->Clone("inverseCorrelation");
+ inverseCorrelation->Reset();
+
+ Float_t inputDistIntegral = 1;
+ if (initialConditions)
+ {
+ printf("Using different starting conditions...\n");
+ inputDistIntegral = initialConditions->Integral();
+ }
+ Float_t measuredIntegral = measured->Integral();
+ for (Int_t me=1; me<=kMaxME; me++)
+ for (Int_t mz=1; mz<=kMaxMZ; mz++)
+ {
+ // normalization of the measured spectrum
+ measuredCopy[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
+ errors[me][mz] = measured->GetBinError(me, mz) / measuredIntegral;
+ // pick prior distribution and normalize it
+ if (initialConditions)
+ prior[me][mz] = initialConditions->GetBinContent(me,mz) / inputDistIntegral;
+ else
+ prior[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
+ }
+
+ // unfold...
+ for (Int_t i=0; i<nIterations; i++)
+ {
+ // calculate Inverse Correlation Map from Bayes theorem:
+ // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
+ /*Float_t norm = 0;
+ for (Int_t me=1; me<=kMaxME; me++)
+ for (Int_t mz=1; mz<=kMaxMZ; mz++)
+ {
+ norm = 0;
+ for (Int_t te=kStartBin; te<=kMaxTE; te++)
+ for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
+ {
+ Int_t bin[4] = {te, me, tz, mz};
+ norm += correlation->GetBinContent(bin)*prior[te][tz];
+ }
+ if (norm > 0)
+ for (Int_t te = kStartBin; te <= kMaxTE; te++)
+ for (Int_t tz = kStartBin; tz <= kMaxTZ; tz++)
+ {
+ Int_t bin[4] = {te, me, tz, mz};
+ inverseCorrelation->SetBinContent(bin, correlation->GetBinContent(bin)*prior[te][tz]/norm );
+ }
+ //else
+ // inverse response set to '0' wich has been already done in line 2069
+ }*/
+ inverseCorrelation->Reset();
+ Float_t norm[kMaxTE+2][kMaxTZ+2];
+ for (Int_t te=0; te<(kMaxTE+2); te++)
+ for (Int_t tz=0; tz<(kMaxTZ+2); tz++)
+ norm[te][tz]=0;
+ for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
+ {
+ Int_t bin[4];
+ Float_t binContent = correlation->GetBinContent(idx, bin);
+ if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ &&
+ bin[0]>0 && bin[0]<=NBINSE && bin[2]>0 && bin[2]<=NBINSZ)
+ norm[bin[1]][bin[3]] += binContent*prior[bin[0]][bin[2]];
+ }
+ Float_t chi2Measured=0, diff;
+ for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
+ {
+ Int_t bin[4];
+ Float_t binContent = correlation->GetBinContent(idx, bin);
+ if (norm[bin[1]][bin[3]]>0 && bin[1]>0 && bin[1]<=NBINSE &&
+ bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ)
+ {
+ inverseCorrelation->SetBinContent(bin, binContent*prior[bin[0]][bin[2]]/norm[bin[1]][bin[3]]);
+ if (errors[bin[1]][bin[3]]>0)
+ {
+ diff = ((measuredCopy[bin[1]][bin[3]]-norm[bin[1]][bin[3]])/(errors[bin[1]][bin[3]]));
+ chi2Measured += diff*diff;
+ }
+ }
+ }
+
+ // calculate "generated" spectrum
+ for (Int_t te = kStartBin; te<=kMaxTE; te++)
+ for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
+ {
+ Float_t value = 0;
+ for (Int_t me=1; me<=kMaxME; me++)
+ for (Int_t mz=1; mz<=kMaxMZ; mz++)
+ {
+ Int_t bin[4] = {te, me, tz, mz};
+ value += inverseCorrelation->GetBinContent(bin)*measuredCopy[me][mz];
+ }
+ result[te][tz] = value;
+ //printf("%e\n", result[te][tz]);
+ }
+
+ // regularization (simple smoothing)
+ Float_t chi2LastIter = 0;
+ for (Int_t te=kStartBin; te<=kMaxTE; te++)
+ for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
+ {
+ Float_t newValue = 0;
+ // 0 bin excluded from smoothing
+ if (( te >(kStartBin+1) && te<(kMaxTE-1) ) && ( tz > (kStartBin+1) && tz<(kMaxTZ-1) ))
+ {
+ 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.;
+
+ // weight the average with the regularization parameter
+ newValue = (1 - regPar) * result[te][tz] + regPar * average;
+ }
+ else
+ newValue = result[te][tz];
+ if (prior[te][tz]>1.e-5)
+ {
+ diff = ((prior[te][tz]-newValue)/prior[te][tz]);
+ chi2LastIter = diff*diff;
+ }
+ prior[te][tz] = newValue;
+ }
+ //printf(" iteration %d - chi2LastIter = %e - chi2Measured = %e \n", i, chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ), chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ));
+ if (chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-6 && chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-3)
+ break;
+ } // end of iterations
+
+ // propagate errors of the reconstructed distribution through the unfolding
+ for (Int_t te = kStartBin; te<=kMaxTE; te++)
+ for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
+ {
+ Float_t valueError = 0;
+ Float_t binError = 0;
+ for (Int_t me=1; me<=kMaxME; me++)
+ for (Int_t mz=1; mz<=kMaxMZ; mz++)
+ {
+ Int_t bin[4] = {te, me, tz, mz};
+ valueError += inverseCorrelation->GetBinContent(bin)*inverseCorrelation->GetBinContent(bin)*errors[me][mz]*errors[me][mz];
+ }
+ //if (errors[te][tz]!=0)printf("errors[%d][%d]=%e\n", te, tz, valueError);
+ aResult->SetBinContent(te, tz, prior[te][tz]);
+ aResult->SetBinError(te, tz, TMath::Sqrt(valueError));
+ }
+
+ // ***********************************************************************************************************
+ // Calculate the covariance matrix, all arguments are taken from G. D'Agostini (p.6-8)
+ if (calculateErrors)
+ {
+ printf("Covariance matrix will be calculated... this will take a lot of time (>1 day) ;)\n");
+
+ //Variables and Matrices that will be use along the calculation
+ const Int_t binsV[4] = {NBINSE,NBINSE, NBINSZ, NBINSZ};
+ const Double_t LowEdgeV[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
+ const Double_t UpEdgeV[4] = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
+
+ const Double_t Ntrue = (Double_t)measured->Integral();
+
+ THnSparseF *V = new THnSparseF("V","",4, binsV, LowEdgeV, UpEdgeV);
+ V->Reset();
+ Double_t invCorrContent1, Nt;
+ Double_t invCorrContent2, v11, v12, v2;
+ // calculate V1 and V2
+ for (Int_t idx1=0; idx1<=NFilledBins; idx1++)
+ {
+ printf("Covariance Matrix calculation: iteration idx1=%d of %d\n", idx1, NFilledBins);
+ for (Int_t idx2=0; idx2<=NFilledBins; idx2++)
+ {
+ Int_t bin1[4];
+ Int_t bin2[4];
+ invCorrContent1 = inverseCorrelation->GetBinContent(idx1, bin1);
+ invCorrContent2 = inverseCorrelation->GetBinContent(idx2, bin2);
+ v11=0; v12=0; v2=0;
+ if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE &&
+ bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ &&
+ bin2[0]>0 && bin2[0]<=NBINSE && bin2[1]>0 && bin2[1]<=NBINSE &&
+ bin2[2]>0 && bin2[2]<=NBINSZ && bin2[3]>0 && bin2[3]<=NBINSZ)
+ {
+ if (bin1[1]==bin2[1] && bin1[3]==bin2[3])
+ v11 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]
+ *(1. - measuredCopy[bin2[1]][bin2[3]]/Ntrue);
+ else
+ v12 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]*
+ measuredCopy[bin2[1]][bin2[3]]/Ntrue;
+ Nt = (Double_t)prior[bin2[0]][bin2[2]];
+ v2 = measuredCopy[bin1[1]][bin1[3]]*measuredCopy[bin2[1]][bin2[3]]*
+ invCorrContent1*invCorrContent2*
+ BayesUncertaintyTerms(inverseCorrelation, correlation, bin1, bin2, Nt);
+ Int_t binV[4] = {bin1[0],bin2[0],bin1[2],bin2[2]};
+ V->SetBinContent(binV,v11-v12 + v2);
+ }
+ }
+ }
+
+ for(Int_t te = 1; te<=NBINSE; te++)
+ for(Int_t tz = 1; tz<=NBINSZ; tz++)
+ {
+ Int_t binV[4] = {te,te,tz,tz};
+ aResult->SetBinError( te, tz, V->GetBinContent(binV) );
+ }
+
+ TFile* f = new TFile("Covariance_UnfSpectrum.root");
+ f->Open("RECREATE");
+ V->Write();
+ f->Close();
+ }
+
+ return 0;
+
+}
+
+//____________________________________________________________________
+Double_t AliJetSpectrumUnfolding::BayesUncertaintyTerms(THnSparseF *M, THnSparseF *C, Int_t* binTM, Int_t* binTM1, Double_t Nt)
+{
+ //
+ // helper function for the covariance matrix of the bayesian method
+ //
+
+ Double_t result = 0;
+ Float_t term[9];
+ Int_t tmpBin[4], tmpBin1[4];
+ const Int_t nFilledBins = C->GetNbins();
+ if (Nt==0)
+ return 0;
+
+ Float_t CorrContent;
+ Float_t InvCorrContent;
+
+ tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
+ tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
+ if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
+ {
+ if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
+ term[0] = BayesCov(M, C, tmpBin, tmpBin1)/
+ (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
+ term[2] = term[0]*M->GetBinContent(tmpBin1);
+ }
+ else
+ {
+ term[0] = 0;
+ term[2] = 0;
+ }
+
+ tmpBin[0]=binTM1[0]; tmpBin[1]=binTM[1]; tmpBin[2]=binTM1[2]; tmpBin[3]=binTM[3];
+ tmpBin1[0]=binTM1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=binTM1[3];
+ if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
+ term[6] = BayesCov(M, C, tmpBin, tmpBin1)*
+ M->GetBinContent(tmpBin)/
+ (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
+ else
+ term[6] = 0;
+
+ for(Int_t idx1=0; idx1<=nFilledBins; idx1++)
+ {
+ Int_t bin1[4];
+ CorrContent = C->GetBinContent(idx1, bin1);
+ InvCorrContent = M->GetBinContent(idx1, bin1);
+ if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE &&
+ bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ)
+ {
+ tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
+ tmpBin1[0]=binTM[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=bin1[3];
+ if (C->GetBinContent(tmpBin)!=0 &&
+ binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
+ term[1] = BayesCov(M, C, tmpBin, tmpBin1)/C->GetBinContent(tmpBin);
+ else
+ term[1] = 0;
+
+ tmpBin[0] =binTM[0]; tmpBin[1] =bin1[1]; tmpBin[2] =binTM[2]; tmpBin[3] =bin1[3];
+ tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
+ if (C->GetBinContent(tmpBin1)!=0)
+ {
+ if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
+ term[3] = BayesCov(M, C, tmpBin, tmpBin1)/
+ C->GetBinContent(tmpBin1);
+ term[5] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin1)/
+ C->GetBinContent(tmpBin1);
+ }
+ else
+ {
+ term[3] = 0;
+ term[5] = 0;
+ }
+
+ tmpBin[0] =binTM1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM1[2]; tmpBin[3] =binTM[3];
+ tmpBin1[0]=binTM1[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=bin1[3];
+ if (C->GetBinContent(tmpBin)!=0)
+ term[7] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin)/
+ C->GetBinContent(tmpBin);
+ else
+ term[7] = 0;
+
+ tmpBin[0] =bin1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =bin1[2]; tmpBin[3] =binTM[3];
+ tmpBin1[0]=bin1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=bin1[2]; tmpBin1[3]=binTM1[3];
+ if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
+ term[8] = BayesCov(M, C, tmpBin, tmpBin1)*
+ M->GetBinContent(tmpBin)*M->GetBinContent(tmpBin)/
+ (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
+ else
+ term[8] = 0;
+
+ for (Int_t i=0; i<9; i++)
+ result += term[i]/Nt;
+ }
+ }
+
+ return result;
+}
+
+//____________________________________________________________________
+Double_t AliJetSpectrumUnfolding::BayesCov(THnSparseF *M, THnSparseF *correlation, Int_t* binTM, Int_t* bin1)
+{
+ Double_t result, result1, result2, result3;
+
+ if (binTM[0]==bin1[0] && binTM[2]==bin1[2])
+ {
+ if (correlation->GetBinContent(bin1)!=0)
+ result1 = 1./correlation->GetBinContent(bin1);
+ else
+ result1 = 0;
+ result2 = 1.;
+ }
+ else
+ {
+ result1 = 0;
+ result2 = 0;
+ }
+
+ if (binTM[1]==bin1[1] && binTM[3]==bin1[3])
+ {
+ Int_t tmpbin[4] = {bin1[0], binTM[1], bin1[2], binTM[3]};
+ if(correlation->GetBinContent(tmpbin)!=0)
+ result3 = M->GetBinContent(tmpbin)/correlation->GetBinContent(tmpbin);
+ else
+ result3 = 0;
+ }
+ else
+ {
+ result1 = 0;
+ result3 = 0;
+ }
+
+ return result = result1 + result2 + result3;
+}
+
+//____________________________________________________________________
+TH2F* AliJetSpectrumUnfolding::CalculateRecSpectrum(TH2* inputGen)
+{
+ // runs the distribution given in inputGen through the correlation histogram identified by
+ // fCorrelation and produces a reconstructed spectrum
+
+ if (!inputGen)
+ return 0;
+
+ // normalize to convert number of events into probability
+ /*for (Int_t te=1; te<=NBINSE; te++)
+ for (Int_t tz=1; tz<=NBINSZ; tz++)
+ {
+ Int_t bin[4];
+ Float_t sum=0.;
+ for (Int_t me = 1; me<=NBINSE; me++)
+ for (Int_t mz = 1; mz<=NBINSZ; mz++)
+ {
+ bin[0] = te; bin[1] = me;
+ bin[2] = tz; bin[3] = mz;
+ sum += fCorrelation[correlationMap]->GetBinContent(bin);
+ }
+ if (sum > 0.)
+ for (Int_t me = 1; me<=NBINSE; me++)
+ for (Int_t mz = 1; mz<=NBINSZ; mz++)
+ {
+ bin[0] = te; bin[1] = me;
+ bin[2] = tz; bin[3] = mz;
+ fCorrelation[correlationMap]->SetBinContent(bin, fCorrelation[correlationMap]->GetBinContent(bin)/sum);
+ fCorrelation[correlationMap]->SetBinError(bin, fCorrelation[correlationMap]->GetBinError(bin)/sum);
+ }
+ }*/
+ // normalize to convert number of events into probability (the following loop is much faster)
+ Float_t sum[NBINSE+2][NBINSZ+2];
+ memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
+
+ for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
+ {
+ Int_t bin[4];
+ Float_t binContent = fCorrelation->GetBinContent(idx, bin);
+ if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ){
+ sum[bin[0]][bin[2]] += binContent;
+ }
+ }
+
+ for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
+ {
+ Int_t bin[4];
+ Float_t binContent = fCorrelation->GetBinContent(idx, bin);
+ Float_t binError = fCorrelation->GetBinError(bin);
+ if (sum[bin[0]][bin[2]]>0 && bin[1]>0 && bin[1]<=NBINSE &&
+ bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ)
+ {
+ fCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
+ fCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
+ }
+ }
+
+ TH2F* target = dynamic_cast<TH2F*> (fRecSpectrum->Clone(Form("reconstructed_%s", inputGen->GetName())));
+ target->Reset();
+
+ for (Int_t me=1; me<=NBINSE; ++me)
+ for (Int_t mz=1; mz<=NBINSZ; ++mz)
+ {
+ Float_t measured = 0;
+ Float_t error = 0;
+
+ for (Int_t te=1; te<=NBINSE; ++te)
+ for (Int_t tz=1; tz<=NBINSZ; ++tz)
+ {
+ Int_t bin[4] = {te, me, tz, mz};
+ measured += inputGen->GetBinContent(te,tz) * fCorrelation->GetBinContent(bin);
+ error += inputGen->GetBinError(te,tz) * fCorrelation->GetBinContent(bin);
+ }
+ target->SetBinContent(me, mz, measured);
+ target->SetBinError(me, mz, error);
+ }
+
+ return target;
+}
+
+//__________________________________________________________________________________________________
+void AliJetSpectrumUnfolding::SetGenRecFromFunc(TF2* inputGen)
+{
+ // uses the given function to fill the input Generated histogram and generates from that
+ // the reconstructed histogram by applying the response histogram
+ // this can be used to evaluate if the methods work indepedently of the input
+ // distribution
+
+ if (!inputGen)
+ return;
+
+ TH2F* histtmp = new TH2F("histtmp", "tmp", EBINNING, ZBINNING);
+ TH2F* gen = fGenSpectrum;
+
+ histtmp->Reset();
+ gen->Reset();
+
+ histtmp->FillRandom(inputGen->GetName(), NEVENTS);
+
+ for (Int_t i=1; i<=gen->GetNbinsX(); ++i)
+ for (Int_t j=1; j<=gen->GetNbinsY(); ++j)
+ {
+ gen->SetBinContent(i, j, histtmp->GetBinContent(i,j));
+ gen->SetBinError(i, j, histtmp->GetBinError(i,j));
+ }
+
+ delete histtmp;
+
+ //new TCanvas;
+ //gStyle->SetPalette(1);
+ //gPad->SetLogz();
+ //gen->Draw("COLZ");
+
+
+ TH2 *recsave = fRecSpectrum;
+
+ fRecSpectrum = CalculateRecSpectrum(gen);
+ fRecSpectrum->SetName(recsave->GetName());
+ delete recsave;
+
+ return;
+}
+//________________________________________________________________________________________
--- /dev/null
+//
+// script with functions to use AliJetSpectrumUnfolding class
+//
+
+
+void loadlibs(){
+ // load all the needed libs to run wit root only
+
+ gSystem->Load("libTree");
+ gSystem->Load("libPhysics");
+ gSystem->Load("libHist");
+ gSystem->Load("libVMC");
+ gSystem->Load("libSTEERBase");
+ gSystem->Load("libAOD");
+ gSystem->Load("libESD");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libJETAN");
+ gSystem->Load("libPWG4JetTasks");
+
+
+}
+
+
+
+void draw(const char* fileName = "unfolded_pwg4spec.root", const char* folder = "unfolding", Bool_t proj = kTRUE)
+{
+
+ loadlibs();
+
+ AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
+
+ TFile::Open(fileName);
+ jetSpec->LoadHistograms();
+
+
+ if (proj)
+ {
+ canvas1 = new TCanvas("Response Map Projection", "Response Map Projection", 500, 500);
+ canvas1->Divide(2);
+
+ Int_t style = 1;
+ const Int_t NRGBs = 5;
+ const Int_t NCont = 500;
+
+ Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+ Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+ Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+ Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+ TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+ gStyle->SetNumberContours(NCont);
+ gStyle->SetPalette(style);
+
+ canvas1->cd(1);
+ gPad->SetLogz();
+ h2 = (jetSpec->GetCorrelation())->Projection(2,3);
+ h2->SetXTitle("z^{lp}_{gen}");
+ h2->SetYTitle("z^{lp}_{rec}");
+ h2->DrawCopy("colz");
+
+ canvas1->cd(2);
+ gPad->SetLogz();
+ h3 = (jetSpec->GetCorrelation())->Projection(0,1);
+ h3->SetXTitle("E^{jet}_{gen} [GeV]");
+ h3->SetYTitle("E^{jet}_{rec} [GeV]");
+ h3->DrawCopy("colz");
+ }
+ jetSpec->DrawComparison("Draw_unfolded_pwg4spec", jetSpec->GetGenSpectrum());
+
+ return;
+}
+
+//________________________________________________________________________________________________________________
+void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameUnf = "unfolded_pwg4spec.root")
+{
+ // function to load jet spectra from the output file of the task AliAnalysisTaskJetSpectrum
+ // to do the unfolding
+
+ loadlibs();
+
+ AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
+
+ TFile::Open(fileNameRec);
+ jetSpec->LoadHistograms();
+
+ TFile::Open(fileNameGen);
+ TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
+ jetSpec->SetGenSpectrum(hist);
+
+ jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
+ // last parameter = calculateErrors <- this method to calculate the errors takes a lot of time
+
+ TFile* file = TFile::Open(fileNameUnf,"RECREATE");
+ jetSpec->SaveHistograms();
+ file->Close();
+}
+
+//___________________________________________________________________________
+void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile)
+{
+ // function to test 2D Bayesian unfolding with other spectra
+ // build from a function
+
+ loadlibs();
+
+
+
+ AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
+
+ TFile::Open(inFile);
+ jetSpec->LoadHistograms("unfolding");
+
+ TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3);
+ h2->DrawCopy("colz");
+ c1->Update();
+ if(getchar()=='q')return;
+
+
+ switch (caseNo)
+ {
+ case 1: func = new TF2("func", "501-x-y"); break;
+ case 2: func = new TF2("func", "1000 * 1/(y+x+1)"); break;
+ case 3: func = new TF2("func", "1./(y*pow(x,5.7))"); break;
+ case 4: func = new TF2("func", "exp(-0.1*x - 0.1*y)"); break;
+ case 5: func = new TF2("func", "x*x + (y**3)/100."); break;
+ case 6: func = new TF2("func", "1000*y*exp(-0.1*x)"); break;
+ case 7: func = new TF2("func", "exp(-((x-100.)/(0.3*x))**2 - ((y-0.6)/(0.8*y))**2)"); break;
+ default: return;
+ }
+
+ //new TCanvas; func->Draw();
+
+ jetSpec->SetGenRecFromFunc(func);
+
+ TFile* file = TFile::Open(outFile,"RECREATE");
+ jetSpec->SaveHistograms();
+
+ h2 = (jetSpec->GetCorrelation())->Projection(0,3);
+ h2->DrawCopy("colz");
+ c1->Update();
+ file->Close();
+
+ //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz");
+ //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz");
+}
+
+//_____________________________________________________________________________________________
+void buildResponseMap(const char* fileName = "responseMap.root")
+{
+ // function to build a Response Map with a gaussian distribution
+ loadlibs();
+
+ AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
+
+ TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)");
+
+ bool bPerfect = false;
+
+ Double_t var[4];
+ Float_t sigmax, sigmay;
+ for (Int_t tx=1; tx<=jetSpec->GetCorrelation()->GetAxis(0)->GetNbins(); tx++)
+ for (Int_t ty=1; ty<=jetSpec->GetCorrelation()->GetAxis(2)->GetNbins(); ty++)
+ {
+ var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx);
+ var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty);
+ sigmax = 0.2*var[0];
+ sigmay = 0.2*var[2];
+ func->SetParameters(var[0],sigmax,var[2],sigmay);
+ for (Int_t mx=1; mx<=jetSpec->GetCorrelation()->GetAxis(1)->GetNbins(); mx++)
+ for (Int_t my=1; my<=jetSpec->GetCorrelation()->GetAxis(3)->GetNbins(); my++)
+ {
+ var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx);
+ var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my);
+
+
+ if(bPerfect){
+ if (var[1]==var[0] && var[3]==var[2])
+ jetSpec->GetCorrelation()->Fill(var,1);
+ }
+ else {
+ // cut at sigma
+ if (TMath::Abs(var[1]-var[0]) < sigmax || TMath::Abs(var[3]-var[2]) < sigmay)
+ jetSpec->GetCorrelation()->Fill(var,func->Eval(var[1],var[3]));
+ }
+ }
+ }
+
+
+ TFile* file = TFile::Open(fileName,"RECREATE");
+ jetSpec->SaveHistograms();
+ file->Close();
+}
+
+//_____________________________________________________________________________
+void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameOut = "Uncertainties2DBayesUnf.root")
+{
+ // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding
+
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libJETAN");
+ gSystem->Load("libPWG4JetTasks");
+
+ AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
+
+ TFile::Open(fileNameRec);
+ jetSpec->LoadHistograms();
+
+ TFile::Open(fileNameGen);
+ TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
+ jetSpec->SetGenSpectrum(hist);
+
+ // create sigma histogram
+ TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma");
+ sigma->Reset();
+
+ THnSparseF* hcorr = (THnSparseF*)(jetSpec->GetCorrelation())->Clone();
+ TH2F* htrue = (TH2F*)(jetSpec->GetGenSpectrum())->Clone();
+ TH2F* hmeas = (TH2F*)(jetSpec->GetRecSpectrum())->Clone();
+ TH2F* hunfo = (TH2F*)(jetSpec->GetUnfSpectrum())->Clone();
+
+ Int_t nIterations = 1000;
+ Float_t binContent;
+ for(Int_t i=0; i<nIterations; i++)
+ {
+ printf("iteration = %d\n", i);
+ // reset histograms
+ jetSpec->GetRecSpectrum()->Reset();
+ jetSpec->GetGenSpectrum()->Reset();
+ jetSpec->GetCorrelation()->Reset();
+ jetSpec->GetUnfSpectrum()->Reset();
+
+ THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr");
+ TH2F* tmptrue = (TH2F*)htrue->Clone("tmptrue");
+
+ jetSpec->SetGenSpectrum(tmptrue);
+ jetSpec->SetCorrelation(tmpcorr);
+
+ // randomize reconstructed distribution
+ for (Int_t me=1; me<=hmeas->GetNbinsX(); me++)
+ for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++)
+ {
+ binContent = hmeas->GetBinContent(me,mz);
+ if (binContent>0)
+ {
+ TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",binContent*0.25, binContent*1.25);
+ poisson->SetParameters(binContent,0.);
+ binContent = poisson->GetRandom();
+ delete poisson;
+ }
+ jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent);
+ }
+
+ // unfold
+ jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0);
+
+ // calculate sigma^2
+ for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
+ for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
+ {
+ if (htrue->GetBinContent(te,tz)!=0)
+ {
+ binContent = (jetSpec->GetUnfSpectrum()->GetBinContent(te,tz) -
+ htrue->GetBinContent(te,tz))/htrue->GetBinContent(te,tz);
+ binContent *= binContent;
+ sigma->SetBinContent(te,tz, binContent + sigma->GetBinContent(te,tz));
+ }
+ }
+ }
+
+ // calculate sigma
+ for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
+ for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
+ {
+ binContent = sigma->GetBinContent(te,tz);
+ binContent = TMath::Sqrt(binContent/(Float_t)nIterations);
+ sigma->SetBinContent(te,tz, binContent);
+ }
+
+ const Int_t NRGBs = 5;
+ const Int_t NCont = 500;
+
+ Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+ Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+ Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+ Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+ TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+ gStyle->SetNumberContours(NCont);
+ gStyle->SetPalette(1);
+
+ new TCanvas();
+ gPad->SetLogz();
+ sigma->SetTitle("#sigma((U_{R} - U)/U)");
+ sigma->SetXTitle("E^{jet} [GeV]");
+ sigma->SetYTitle("z^{lp}");
+ sigma->DrawCopy();
+
+ TFile* file = TFile::Open(fileNameOut,"RECREATE");
+ sigma->Write();
+ file->Close();
+}
+
+//_______________________________________________________________________________________________________________
+void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
+{
+ // This functions avoids problems when the number of bins or the bin limits
+ // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
+ // are different.
+
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libJETAN");
+ gSystem->Load("libPWG4JetTasks");
+
+ file = new TFile(fileNameSpec,"read");
+ tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
+ THnSparseF *fhCorrelation = (THnSparseF*)(tlist->FindObject("fhCorrelation_less5tracks"));
+ THnSparseF *fhCorrelation2 = (THnSparseF*)(tlist->FindObject("fhCorrelation_5to10tracks"));
+ THnSparseF *fhCorrelation3 = (THnSparseF*)(tlist->FindObject("fhCorrelation_more10tracks"));
+ TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fhEGenZGen"));
+ TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fhERecZRec"));
+
+ fhCorrelation->Add(fhCorrelation2, 1);
+ fhCorrelation->Add(fhCorrelation3, 1);
+
+ delete fhCorrelation2;
+ delete fhCorrelation3;
+
+ AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
+
+ // generated jets (true distribution)
+ for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
+ for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
+ {
+ Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
+ Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
+ Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
+ Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
+ Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
+ Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
+ jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
+ jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
+ }
+ file = TFile::Open("gen_pwg4spec.root", "RECREATE");
+ jetSpec->SaveHistograms();
+ file->Close();
+ jetSpec->GetGenSpectrum()->Reset();
+ printf("true distribution has been set\n");
+
+ // reconstructed jets (measured distribution)
+ for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
+ for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
+ {
+ Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
+ Float_t zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
+ Int_t bine = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
+ Int_t binz = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
+ Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
+ Float_t err = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
+ jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
+ jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
+ }
+
+ // Response function
+ jetSpec->GetCorrelation()->Reset();
+ jetSpec->GetCorrelation()->Sumw2();
+
+ for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
+ {
+ //printf("idx: %d\n",idx);
+ Double_t var[4];
+ Int_t bin[4];
+ Float_t BinContent = fhCorrelation->GetBinContent(idx, bin);
+ var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]);
+ var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]);
+ var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]);
+ var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]);
+ bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]);
+ bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]);
+ bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]);
+ bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]);
+ Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin);
+ Float_t err = jetSpec->GetCorrelation()->GetBinError(bin);
+ jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
+ jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
+ }
+
+ file = TFile::Open("rec_pwg4spec.root", "RECREATE");
+ jetSpec->SaveHistograms();
+ file->Close();
+
+ printf("reconstructed distribution has been set\n");
+ printf("response map has been set\n");
+
+}
+
+
+
+
+void correct(){
+ // simple steering to correct a given distribution;
+ loadlibs();
+ FillSpecFromFile("/home/ydelgado/pcalice014.cern.ch/macros/jets/CAFdata/histos_pwg4spec.root");
+
+ char name[100];
+ sprintf(name, "unfolded_pwg4spec.root");
+
+ unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name);
+ //draw(name, "unfolding", 1);
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff