#include <TH2F.h>
#include <TH3F.h>
#include <TDirectory.h>
-#include <TVirtualFitter.h>
#include <TCanvas.h>
#include <TString.h>
#include <TF1.h>
ClassImp(AliMultiplicityCorrection)
-const Int_t AliMultiplicityCorrection::fgkMaxInput = 250; // bins in measured histogram
-const Int_t AliMultiplicityCorrection::fgkMaxParams = 250; // bins in unfolded histogram = number of fit params
-
-TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0;
-TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0;
-TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0;
-TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0;
-
-Bool_t AliMultiplicityCorrection::fgCreateBigBin = kTRUE;
-AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1;
-Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000;
-Int_t AliMultiplicityCorrection::fgNParamsMinuit = AliMultiplicityCorrection::fgkMaxParams;
-
-TF1* AliMultiplicityCorrection::fgNBD = 0;
-
-Float_t AliMultiplicityCorrection::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
-Int_t AliMultiplicityCorrection::fgBayesianIterations = 100; // number of iterations in Bayesian method
-
// These are the areas where the quality of the unfolding results are evaluated
// Default defined here, call SetQualityRegions to change them
// unit is in multiplicity (not in bin!)
// SPD: peak area - flat area - low stat area
-Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4, 60, 180};
-Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210};
+Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {1, 20, 70};
+Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {10, 65, 80};
//____________________________________________________________________
void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
if (SPDStudy)
{
// SPD: peak area - flat area - low stat area
- fgQualityRegionsB[0] = 4;
- fgQualityRegionsE[0] = 20;
+ fgQualityRegionsB[0] = 1;
+ fgQualityRegionsE[0] = 10;
- fgQualityRegionsB[1] = 60;
- fgQualityRegionsE[1] = 140;
+ fgQualityRegionsB[1] = 20;
+ fgQualityRegionsE[1] = 65;
- fgQualityRegionsB[2] = 180;
- fgQualityRegionsE[2] = 210;
+ fgQualityRegionsB[2] = 70;
+ fgQualityRegionsE[2] = 80;
Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
}
fMultiplicityVtx[i] = 0;
fMultiplicityMB[i] = 0;
fMultiplicityINEL[i] = 0;
+ fMultiplicityNSD[i] = 0;
}
for (Int_t i = 0; i < kCorrHists; ++i)
fQuality[i] = 0;
}
+//____________________________________________________________________
+AliMultiplicityCorrection* AliMultiplicityCorrection::Open(const char* fileName, const char* folderName)
+{
+ // opens the given file, reads the multiplicity from the given folder and returns the object
+
+ TFile* file = TFile::Open(fileName);
+ if (!file)
+ {
+ Printf("ERROR: Could not open %s", fileName);
+ return 0;
+ }
+
+ AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folderName, folderName);
+ mult->LoadHistograms();
+
+ // TODO closing the file does not work here, because the histograms cannot be read anymore. LoadHistograms need to be adapted
+
+ return mult;
+}
+
//____________________________________________________________________
AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
TNamed(name, title),
#define VTXBINNING 10, binLimitsVtx
#define NBINNING fgkMaxParams, binLimitsN*/
- #define NBINNING 501, -0.5, 500.5
- #define VTXBINNING 1, -10, 10
+ #define NBINNING 201, -0.5, 200.5
+
+ Double_t vtxRange[] = { 15, 6, 2 };
for (Int_t i = 0; i < kESDHists; ++i)
- fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
+ fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 1, -vtxRange[i], vtxRange[i], NBINNING);
for (Int_t i = 0; i < kMCHists; ++i)
{
- fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
+ fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[i%3]->Clone(Form("fMultiplicityVtx%d", i)));
fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
- fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
+ fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityMB%d", i)));
fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
- fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
+ fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityINEL%d", i)));
fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
+
+ fMultiplicityNSD[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityNSD%d", i)));
+ fMultiplicityNSD[i]->SetTitle("fMultiplicityNSD");
}
for (Int_t i = 0; i < kCorrHists; ++i)
{
- fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
+ fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 1, -vtxRange[i%3], vtxRange[i%3], NBINNING, NBINNING);
fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
}
TH1::AddDirectory(oldStatus);
+
+ AliUnfolding::SetNbins(120, 120);
+ AliUnfolding::SetSkipBinsBegin(1);
+ AliUnfolding::SetNormalizeInput(kTRUE);
}
//____________________________________________________________________
if (fMultiplicityINEL[i])
delete fMultiplicityINEL[i];
fMultiplicityINEL[i] = 0;
- }
+
+ if (fMultiplicityNSD[i])
+ delete fMultiplicityNSD[i];
+ fMultiplicityNSD[i] = 0;
+}
for (Int_t i = 0; i < kCorrHists; ++i)
{
TObject* obj;
// collections of all histograms
- TList collections[kESDHists+kMCHists*3+kCorrHists*2];
+ TList collections[kESDHists+kMCHists*4+kCorrHists*2];
Int_t count = 0;
while ((obj = iter->Next())) {
collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
+ collections[kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
}
for (Int_t i = 0; i < kCorrHists; ++i)
- collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
+ collections[kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
for (Int_t i = 0; i < kCorrHists; ++i)
- collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
+ collections[kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
count++;
}
fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
+ fMultiplicityNSD[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
}
for (Int_t i = 0; i < kCorrHists; ++i)
- fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
+ fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*4+i]);
for (Int_t i = 0; i < kCorrHists; ++i)
- fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
+ fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*4+kCorrHists+i]);
delete iter;
oldObjects.Add(fMultiplicityMB[i]);
if (fMultiplicityINEL[i])
oldObjects.Add(fMultiplicityINEL[i]);
+ if (fMultiplicityNSD[i])
+ oldObjects.Add(fMultiplicityNSD[i]);
}
for (Int_t i = 0; i < kCorrHists; ++i)
fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
+ fMultiplicityNSD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityNSD[i]->GetName()));
if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
success = kFALSE;
}
for (Int_t i = 0; i < kESDHists; ++i)
if (fMultiplicityESD[i])
+ {
fMultiplicityESD[i]->Write();
+ fMultiplicityESD[i]->ProjectionY(Form("%s_px", fMultiplicityESD[i]->GetName()), 1, fMultiplicityESD[i]->GetNbinsX())->Write();
+ }
for (Int_t i = 0; i < kMCHists; ++i)
{
if (fMultiplicityVtx[i])
+ {
fMultiplicityVtx[i]->Write();
+ fMultiplicityVtx[i]->ProjectionY(Form("%s_px", fMultiplicityVtx[i]->GetName()), 1, fMultiplicityVtx[i]->GetNbinsX())->Write();
+ }
if (fMultiplicityMB[i])
+ {
fMultiplicityMB[i]->Write();
+ fMultiplicityMB[i]->ProjectionY(Form("%s_px", fMultiplicityMB[i]->GetName()), 1, fMultiplicityMB[i]->GetNbinsX())->Write();
+ }
if (fMultiplicityINEL[i])
+ {
fMultiplicityINEL[i]->Write();
+ fMultiplicityINEL[i]->ProjectionY(Form("%s_px", fMultiplicityINEL[i]->GetName()), 1, fMultiplicityINEL[i]->GetNbinsX())->Write();
+ }
+ if (fMultiplicityNSD[i])
+ {
+ fMultiplicityNSD[i]->Write();
+ fMultiplicityNSD[i]->ProjectionY(Form("%s_px", fMultiplicityNSD[i]->GetName()), 1, fMultiplicityNSD[i]->GetNbinsX())->Write();
+ }
}
for (Int_t i = 0; i < kCorrHists; ++i)
}
//____________________________________________________________________
-void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
+void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, AliPWG0Helper::MCProcessType processType, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll)
{
//
// Fills an event from MC
{
fMultiplicityMB[0]->Fill(vtx, generated05);
fMultiplicityMB[1]->Fill(vtx, generated10);
- fMultiplicityMB[2]->Fill(vtx, generated15);
- fMultiplicityMB[3]->Fill(vtx, generated20);
- fMultiplicityMB[4]->Fill(vtx, generatedAll);
+ fMultiplicityMB[2]->Fill(vtx, generated14);
+ fMultiplicityMB[3]->Fill(vtx, generatedAll);
if (vertex)
{
fMultiplicityVtx[0]->Fill(vtx, generated05);
fMultiplicityVtx[1]->Fill(vtx, generated10);
- fMultiplicityVtx[2]->Fill(vtx, generated15);
- fMultiplicityVtx[3]->Fill(vtx, generated20);
- fMultiplicityVtx[4]->Fill(vtx, generatedAll);
+ fMultiplicityVtx[2]->Fill(vtx, generated14);
+ fMultiplicityVtx[3]->Fill(vtx, generatedAll);
}
}
fMultiplicityINEL[0]->Fill(vtx, generated05);
fMultiplicityINEL[1]->Fill(vtx, generated10);
- fMultiplicityINEL[2]->Fill(vtx, generated15);
- fMultiplicityINEL[3]->Fill(vtx, generated20);
- fMultiplicityINEL[4]->Fill(vtx, generatedAll);
+ fMultiplicityINEL[2]->Fill(vtx, generated14);
+ fMultiplicityINEL[3]->Fill(vtx, generatedAll);
+
+ if (processType != AliPWG0Helper::kSD)
+ {
+ fMultiplicityNSD[0]->Fill(vtx, generated05);
+ fMultiplicityNSD[1]->Fill(vtx, generated10);
+ fMultiplicityNSD[2]->Fill(vtx, generated14);
+ fMultiplicityNSD[3]->Fill(vtx, generatedAll);
+ }
}
//____________________________________________________________________
-void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
+void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured14)
{
//
// Fills an event from ESD
fMultiplicityESD[0]->Fill(vtx, measured05);
fMultiplicityESD[1]->Fill(vtx, measured10);
- fMultiplicityESD[2]->Fill(vtx, measured15);
- fMultiplicityESD[3]->Fill(vtx, measured20);
+ fMultiplicityESD[2]->Fill(vtx, measured14);
}
//____________________________________________________________________
-void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
+void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated14, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured14)
{
//
// Fills an event into the correlation map with the information from MC and ESD
fCorrelation[0]->Fill(vtx, generated05, measured05);
fCorrelation[1]->Fill(vtx, generated10, measured10);
- fCorrelation[2]->Fill(vtx, generated15, measured15);
- fCorrelation[3]->Fill(vtx, generated20, measured20);
-
- fCorrelation[4]->Fill(vtx, generatedAll, measured05);
- fCorrelation[5]->Fill(vtx, generatedAll, measured10);
- fCorrelation[6]->Fill(vtx, generatedAll, measured15);
- fCorrelation[7]->Fill(vtx, generatedAll, measured20);
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // prefers constant function (pol0)
-
- Double_t chi2 = 0;
-
- // ignore the first bin here. on purpose...
- for (Int_t i=2; i<fgNParamsMinuit; ++i)
- {
- Double_t right = params[i];
- Double_t left = params[i-1];
-
- if (right != 0)
- {
- Double_t diff = 1 - left / right;
- chi2 += diff * diff;
- }
- }
-
- return chi2 / 100.0;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // prefers linear function (pol1)
-
- Double_t chi2 = 0;
-
- // ignore the first bin here. on purpose...
- for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
- {
- if (params[i-1] == 0)
- continue;
-
- Double_t right = params[i];
- Double_t middle = params[i-1];
- Double_t left = params[i-2];
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2) / middle;
-
- chi2 += diff * diff;
- }
-
- return chi2;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // prefers linear function (pol1)
-
- Double_t chi2 = 0;
-
- /*Float_t der[fgkMaxParams];
-
- for (Int_t i=0; i<fgkMaxParams-1; ++i)
- der[i] = params[i+1] - params[i];
-
- for (Int_t i=0; i<fgkMaxParams-6; ++i)
- {
- for (Int_t j=1; j<=5; ++j)
- {
- Double_t diff = der[i+j] - der[i];
- chi2 += diff * diff;
- }
- }*/
-
- // ignore the first bin here. on purpose...
- for (Int_t i=2+1; i<fgNParamsMinuit; ++i)
- {
- if (params[i-1] == 0)
- continue;
-
- Double_t right = log(params[i]);
- Double_t middle = log(params[i-1]);
- Double_t left = log(params[i-2]);
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2) / middle;
+ fCorrelation[2]->Fill(vtx, generated14, measured14);
- chi2 += diff * diff;
- }
-
- return chi2 * 100;
+ fCorrelation[3]->Fill(vtx, generatedAll, measured05);
+ fCorrelation[4]->Fill(vtx, generatedAll, measured10);
+ fCorrelation[5]->Fill(vtx, generatedAll, measured14);
}
//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
- // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
-
- Double_t chi2 = 0;
-
- // ignore the first bin here. on purpose...
- for (Int_t i=3; i<fgNParamsMinuit; ++i)
- {
- Double_t right = params[i];
- Double_t middle = params[i-1];
- Double_t left = params[i-2];
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2);
-
- chi2 += diff * diff;
- }
-
- return chi2 * 1e4;
-}
-
-//____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
-{
- // homogenity term for minuit fitting
- // pure function of the parameters
- // calculates entropy, from
- // The method of reduced cross-entropy (M. Schmelling 1993)
-
- Double_t paramSum = 0;
- // ignore the first bin here. on purpose...
- for (Int_t i=1; i<fgNParamsMinuit; ++i)
- paramSum += params[i];
-
- Double_t chi2 = 0;
- for (Int_t i=1; i<fgNParamsMinuit; ++i)
- {
- //Double_t tmp = params[i] / paramSum;
- Double_t tmp = params[i];
- if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
- chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
- }
-
- return 100.0 + chi2;
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
-{
- //
- // fit function for minuit
- // does: nbd
- // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
- // func->SetParNames("scaling", "averagen", "k");
- // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
- // func->SetParLimits(1, 0.001, 1000);
- // func->SetParLimits(2, 0.001, 1000);
- //
-
- fgNBD->SetParameters(params[0], params[1], params[2]);
-
- Double_t params2[fgkMaxParams];
-
- for (Int_t i=0; i<fgkMaxParams; ++i)
- params2[i] = fgNBD->Eval(i);
-
- MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
-
- printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
-{
- //
- // fit function for minuit
- // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
- //
-
- // d
- static TVectorD paramsVector(fgNParamsMinuit);
- for (Int_t i=0; i<fgNParamsMinuit; ++i)
- paramsVector[i] = params[i] * params[i];
-
- // calculate penalty factor
- Double_t penaltyVal = 0;
- switch (fgRegularizationType)
- {
- case kNone: break;
- case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
- case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
- case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
- case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
- case kLog: penaltyVal = RegularizationLog(paramsVector); break;
- }
-
- //if (penaltyVal > 1e10)
- // paramsVector2.Print();
-
- penaltyVal *= fgRegularizationWeight;
-
- // Ad
- TVectorD measGuessVector(fgkMaxInput);
- measGuessVector = (*fgCorrelationMatrix) * paramsVector;
-
- // Ad - m
- measGuessVector -= (*fgCurrentESDVector);
-
- //measGuessVector.Print();
-
- TVectorD copy(measGuessVector);
-
- // (Ad - m) W
- // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
- // normal way is like this:
- // measGuessVector *= (*fgCorrelationCovarianceMatrix);
- // optimized way like this:
- for (Int_t i=0; i<fgkMaxInput; ++i)
- measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
-
- //measGuessVector.Print();
-
- // reduce influence of first bin
- copy(0) *= 0.01;
-
- // (Ad - m) W (Ad - m)
- // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
- // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
- Double_t chi2FromFit = measGuessVector * copy * 1e6;
-
- chi2 = chi2FromFit + penaltyVal;
-
- static Int_t callCount = 0;
- if ((callCount++ % 10000) == 0)
- printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::DrawGuess(Double_t *params)
-{
- //
- // draws residuals of solution suggested by params and effect of regularization
- //
-
- // d
- static TVectorD paramsVector(fgNParamsMinuit);
- for (Int_t i=0; i<fgNParamsMinuit; ++i)
- paramsVector[i] = params[i] * params[i];
-
- // Ad
- TVectorD measGuessVector(fgkMaxInput);
- measGuessVector = (*fgCorrelationMatrix) * paramsVector;
-
- // Ad - m
- measGuessVector -= (*fgCurrentESDVector);
-
- TVectorD copy(measGuessVector);
-
- // (Ad - m) W
- // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
- // normal way is like this:
- // measGuessVector *= (*fgCorrelationCovarianceMatrix);
- // optimized way like this:
- for (Int_t i=0; i<fgkMaxInput; ++i)
- measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
-
- // (Ad - m) W (Ad - m)
- // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
- // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
- //Double_t chi2FromFit = measGuessVector * copy * 1e6;
-
- // draw residuals
- TH1* residuals = new TH1F("residuals", "residuals", fgkMaxInput, -0.5, fgkMaxInput - 0.5);
- for (Int_t i=0; i<fgkMaxInput; ++i)
- residuals->SetBinContent(i+1, measGuessVector(i) * copy(i) * 1e6);
- new TCanvas; residuals->DrawCopy(); gPad->SetLogy();
-
- // draw penalty
- if (fgRegularizationType != kPol1) {
- Printf("Drawing not supported");
- return;
- }
-
- TH1* penalty = new TH1F("penalty", "penalty", fgNParamsMinuit, -0.5, fgNParamsMinuit - 0.5);
-
- for (Int_t i=2; i<fgNParamsMinuit; ++i)
- {
- if (params[i-1] == 0)
- continue;
-
- Double_t right = params[i];
- Double_t middle = params[i-1];
- Double_t left = params[i-2];
-
- Double_t der1 = (right - middle);
- Double_t der2 = (middle - left);
-
- Double_t diff = (der1 - der2) / middle;
-
- penalty->SetBinContent(i, diff * diff);
- }
- new TCanvas; penalty->DrawCopy(); gPad->SetLogy();
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
+void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
{
//
// fills fCurrentESD, fCurrentCorrelation
// resets fMultiplicityESDCorrected
//
- Bool_t display = kFALSE;
-
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
fMultiplicityESDCorrected[correlationID]->Reset();
}
}
- fCurrentCorrelation = hist->Project3D("zy");
- //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
- //fMultiplicityESDCorrected[correlationID]->Rebin(2);
+ fCurrentCorrelation = (TH2*) hist->Project3D("zy");
fCurrentCorrelation->Sumw2();
-
+
Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
- if (createBigBin)
- {
- fLastBinLimit = 0;
- for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
- {
- if (fCurrentESD->GetBinContent(i) <= 5)
- {
- fLastBinLimit = i;
- break;
- }
- }
-
- if (fLastBinLimit > 0)
- {
- TCanvas* canvas = 0;
-
- if (display)
- {
- canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
- canvas->Divide(2, 2);
-
- canvas->cd(1);
- fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
- fCurrentESD->SetStats(kFALSE);
- fCurrentESD->SetTitle(";measured multiplicity");
- fCurrentESD->DrawCopy();
- gPad->SetLogy();
-
- canvas->cd(2);
- fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
- fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
- fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
- fCurrentCorrelation->SetStats(kFALSE);
- fCurrentCorrelation->DrawCopy("COLZ");
- }
-
- printf("Bin limit in measured spectrum is %d.\n", fLastBinLimit);
- fCurrentESD->SetBinContent(fLastBinLimit, fCurrentESD->Integral(fLastBinLimit, fCurrentESD->GetNbinsX()));
- for (Int_t i=fLastBinLimit+1; i<=fCurrentESD->GetNbinsX(); ++i)
- {
- fCurrentESD->SetBinContent(i, 0);
- fCurrentESD->SetBinError(i, 0);
- }
- // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
- fCurrentESD->SetBinError(fLastBinLimit, TMath::Sqrt(fCurrentESD->GetBinContent(fLastBinLimit)));
-
- printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(fLastBinLimit), fCurrentESD->GetBinError(fLastBinLimit));
-
- for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
- {
- fCurrentCorrelation->SetBinContent(i, fLastBinLimit, fCurrentCorrelation->Integral(i, i, fLastBinLimit, fCurrentCorrelation->GetNbinsY()));
- // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
- fCurrentCorrelation->SetBinError(i, fLastBinLimit, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, fLastBinLimit)));
-
- for (Int_t j=fLastBinLimit+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
- {
- fCurrentCorrelation->SetBinContent(i, j, 0);
- fCurrentCorrelation->SetBinError(i, j, 0);
- }
- }
-
- printf("Adjusted correlation matrix!\n");
-
- if (canvas)
- {
- canvas->cd(3);
- fCurrentESD->DrawCopy();
- gPad->SetLogy();
-
- canvas->cd(4);
- fCurrentCorrelation->DrawCopy("COLZ");
- }
- }
- }
-
#if 0 // does not help
// null bins with one entry
Int_t nNulledBins = 0;
//
// calculates efficiency for given event type
//
-
+
TH1* divisor = 0;
switch (eventType)
{
- case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
- case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
- case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
+ case kTrVtx : break;
+ case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
+ case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
+ case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY("divisor", 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
}
- TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
- eff->Divide(eff, divisor, 1, 1, "B");
- return eff;
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight, Int_t minuitParams)
-{
- //
- // sets the parameters for chi2 minimization
- //
-
- fgRegularizationType = type;
- fgRegularizationWeight = weight;
-
- printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
-
- if (minuitParams > 0)
+ TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
+
+ if (eventType == kTrVtx)
{
- fgNParamsMinuit = minuitParams;
- printf("AliMultiplicityCorrection::SetRegularizationParameters --> Number of MINUIT iterations set to %d\n", minuitParams);
+ for (Int_t i=0; i<= eff->GetNbinsX()+1; i++)
+ eff->SetBinContent(i, 1);
}
+ else
+ eff->Divide(eff, divisor, 1, 1, "B");
+
+ return eff;
}
//____________________________________________________________________
-void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
+TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange)
{
//
- // sets the parameters for Bayesian unfolding
+ // calculates efficiency for given event type
//
- fgBayesianSmoothing = smoothing;
- fgBayesianIterations = nIterations;
-
- printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
+ TH1* divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e");
+ TH1* eff = fMultiplicityMB[inputRange]->ProjectionY("CurrentEfficiency", 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
+ eff->Divide(eff, divisor, 1, 1, "B");
+ return eff;
}
//____________________________________________________________________
-Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
+Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
{
//
- // implemenation of unfolding (internal function)
- //
- // unfolds <measured> using response from <correlation> and effiency <aEfficiency>
- // output is in <result>
- // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
- // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
+ // correct spectrum using minuit chi2 method
//
- // returns minuit status (0 = success), (-1 when check was set)
+ // for description of parameters, see AliUnfolding::Unfold
//
- //normalize measured
- measured->Scale(1.0 / measured->Integral());
-
- if (!fgCorrelationMatrix)
- fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgNParamsMinuit);
- if (!fgCorrelationCovarianceMatrix)
- fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
- if (!fgCurrentESDVector)
- fgCurrentESDVector = new TVectorD(fgkMaxInput);
- if (!fgEntropyAPriori)
- fgEntropyAPriori = new TVectorD(fgNParamsMinuit);
-
- fgCorrelationMatrix->Zero();
- fgCorrelationCovarianceMatrix->Zero();
- fgCurrentESDVector->Zero();
- fgEntropyAPriori->Zero();
-
- // normalize correction for given nPart
- for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
- {
- Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
- if (sum <= 0)
- continue;
- Float_t maxValue = 0;
- Int_t maxBin = -1;
- for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
- {
- // find most probably value
- if (maxValue < correlation->GetBinContent(i, j))
- {
- maxValue = correlation->GetBinContent(i, j);
- maxBin = j;
- }
-
- // npart sum to 1
- correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
- correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
-
- if (i <= fgNParamsMinuit && j <= fgkMaxInput)
- (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
- }
-
- //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
- }
-
- // Initialize TMinuit via generic fitter interface
- TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgNParamsMinuit);
- Double_t arglist[100];
-
- // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
- arglist[0] = 0;
- minuit->ExecuteCommand("SET PRINT", arglist, 1);
-
- // however, enable warnings
- //minuit->ExecuteCommand("SET WAR", arglist, 0);
-
- // set minimization function
- minuit->SetFCN(MinuitFitFunction);
-
- for (Int_t i=0; i<fgNParamsMinuit; i++)
- (*fgEntropyAPriori)[i] = 1;
-
- if (!initialConditions) {
- initialConditions = measured;
- } else {
- printf("Using different starting conditions...\n");
- //new TCanvas; initialConditions->DrawCopy();
+ Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
+
+ AliUnfolding::SetCreateOverflowBin(5);
+ AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+
+ if (!initialConditions)
+ {
+ initialConditions = (TH1*) fCurrentESD->Clone("initialConditions");
initialConditions->Scale(1.0 / initialConditions->Integral());
- }
-
- for (Int_t i=0; i<fgNParamsMinuit; i++)
- if (initialConditions->GetBinContent(i+1) > 0)
- (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
-
- Double_t* results = new Double_t[fgNParamsMinuit+1];
- for (Int_t i=0; i<fgNParamsMinuit; ++i)
- {
- results[i] = initialConditions->GetBinContent(i+1);
-
- // minimum value
if (!check)
- results[i] = TMath::Max(results[i], 1e-3);
-
- Float_t stepSize = 0.1;
-
- // minuit sees squared values to prevent it from going negative...
- results[i] = TMath::Sqrt(results[i]);
- //stepSize /= results[i]; // keep relative step size
-
- minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
- }
- //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0);
- // bin 0 is filled with value from bin 1 (otherwise it's 0)
- //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
- //results[0] = 0;
- //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
-
- for (Int_t i=0; i<fgkMaxInput; ++i)
- {
- (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
- if (measured->GetBinError(i+1) > 0)
- (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
-
- if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
- (*fgCorrelationCovarianceMatrix)(i, i) = 0;
- }
-
- Int_t dummy = 0;
- Double_t chi2 = 0;
- MinuitFitFunction(dummy, 0, chi2, results, 0);
- DrawGuess(results);
- printf("Chi2 of initial parameters is = %f\n", chi2);
-
- if (check)
- {
- delete[] results;
- return -1;
- }
-
- // first param is number of iterations, second is precision....
- arglist[0] = 1e6;
- //arglist[1] = 1e-5;
- //minuit->ExecuteCommand("SCAN", arglist, 0);
- Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
- printf("MINUIT status is %d\n", status);
- //minuit->ExecuteCommand("MIGRAD", arglist, 1);
- //minuit->ExecuteCommand("MIGRAD", arglist, 1);
- //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
- //minuit->ExecuteCommand("MINOS", arglist, 0);
-
- for (Int_t i=0; i<fgNParamsMinuit; ++i)
- {
- if (aEfficiency->GetBinContent(i+1) > 0)
- {
- results[i] = minuit->GetParameter(i);
- result->SetBinContent(i+1, results[i] * results[i] / aEfficiency->GetBinContent(i+1));
- // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
- result->SetBinError(i+1, minuit->GetParError(i) * results[i] / aEfficiency->GetBinContent(i+1));
- }
- else
{
- result->SetBinContent(i+1, 0);
- result->SetBinError(i+1, 0);
+ // set minimum value to prevent MINUIT just staying in the small value
+ for (Int_t i=1; i<=initialConditions->GetNbinsX(); i++)
+ initialConditions->SetBinContent(i, TMath::Max(initialConditions->GetBinContent(i), 1e-3));
}
}
- DrawGuess(results);
-
- delete[] results;
- return status;
+ return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
}
//____________________________________________________________________
-Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
+Int_t AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
{
//
- // correct spectrum using minuit chi2 method
+ // correct spectrum using minuit chi2 method with a NBD function
//
- // for description of parameters, see UnfoldWithMinuit
+ // for description of parameters, see AliUnfolding::Unfold
//
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType, fgCreateBigBin);
-
- return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
-}
-
-//____________________________________________________________________
-void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
-{
- //
- // correct spectrum using minuit chi2 method applying a NBD fit
- //
-
- Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
-
- SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
- //normalize ESD
- fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
-
- fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
- fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
- fgCurrentESDVector = new TVectorD(fgkMaxParams);
-
- // normalize correction for given nPart
- for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
- {
- Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
- if (sum <= 0)
- continue;
- for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
- {
- // npart sum to 1
- fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
- fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
-
- if (i <= fgkMaxParams && j <= fgkMaxParams)
- (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
- }
- }
-
- for (Int_t i=0; i<fgkMaxParams; ++i)
- {
- (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
- if (fCurrentESD->GetBinError(i+1) > 0)
- (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
- }
-
- // Create NBD function
- if (!fgNBD)
- {
- fgNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250);
- fgNBD->SetParNames("scaling", "averagen", "k");
- }
-
- // Initialize TMinuit via generic fitter interface
- TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
-
- minuit->SetFCN(MinuitNBD);
- minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
- minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
- minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
-
- Double_t arglist[100];
- arglist[0] = 0;
- minuit->ExecuteCommand("SET PRINT", arglist, 1);
- minuit->ExecuteCommand("MIGRAD", arglist, 0);
- //minuit->ExecuteCommand("MINOS", arglist, 0);
-
- fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
-
- for (Int_t i=0; i<fgkMaxParams; ++i)
- {
- printf("%d : %f\n", i, fgNBD->Eval(i));
- if (fgNBD->Eval(i) > 0)
- {
- fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
- fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
- }
- }
+
+ AliUnfolding::SetUnfoldingMethod(AliUnfolding::kFunction);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+
+ TF1* func = new TF1("nbd", "[0] * TMath::Gamma([2]+x) / TMath::Gamma([2]) / TMath::Gamma(x+1) * pow([1] / ([1]+[2]), x) * pow(1.0 + [1]/[2], -[2])");
+ func->SetParNames("scaling", "averagen", "k");
+ func->SetParLimits(0, 0, 1000);
+ func->SetParLimits(1, 1, 50);
+ func->SetParLimits(2, 1, 10);
+ func->SetParameters(1, 10, 2);
+ AliUnfolding::SetFunction(func);
+
+ return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, 0, fMultiplicityESDCorrected[correlationID]);
}
//____________________________________________________________________
}
//____________________________________________________________________
-void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple)
+void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple, EventType eventType)
{
// draw comparison plots
break;
}
Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
-
+
// scale to 1
mcHist->Sumw2();
if (mcHist->Integral() > 0)
// for that we convolute the response matrix with the gathered result
// if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
+
+ // undo trigger, vertex efficiency correction
+ fCurrentEfficiency = GetEfficiency(inputRange, eventType);
tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
+
TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
if (convolutedProj->Integral() > 0)
TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
+ // find bin limit
+ Int_t lastBin = 0;
+ for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
+ {
+ if (fCurrentESD->GetBinContent(i) <= 5)
+ {
+ lastBin = i;
+ break;
+ }
+ }
+
// TODO fix errors
Float_t chi2 = 0;
for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
if (fCurrentESD->GetBinError(i) > 0)
{
Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
- if (i > 1 && i <= fLastBinLimit)
+ if (i > 1 && i <= lastBin)
chi2 += value * value;
Printf("%d --> %f (%f)", i, value * value, chi2);
residual->SetBinContent(i, value);
//residualHist->Fit("gaus", "N");
//delete residualHist;
- printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit);
+ printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, lastBin);
TCanvas* canvas1 = 0;
if (simple)
{
- canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
+ canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 600);
canvas1->Divide(2, 1);
}
else
}
canvas1->cd(1);
+ canvas1->cd(1)->SetGridx();
+ canvas1->cd(1)->SetGridy();
+ canvas1->cd(1)->SetTopMargin(0.05);
+ canvas1->cd(1)->SetRightMargin(0.05);
+ canvas1->cd(1)->SetLeftMargin(0.12);
+ canvas1->cd(1)->SetBottomMargin(0.12);
TH1* proj = (TH1*) mcHist->Clone("proj");
// normalize without 0 bin
proj->Scale(1.0 / proj->Integral(2, proj->GetNbinsX()));
Printf("Normalized without 0 bin!");
proj->GetXaxis()->SetRangeUser(0, 200);
- proj->SetTitle(";true multiplicity;Entries");
+ proj->GetYaxis()->SetTitleOffset(1.4);
+ //proj->SetLabelSize(0.05, "xy");
+ //proj->SetTitleSize(0.05, "xy");
+ proj->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
proj->SetStats(kFALSE);
fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
+ fMultiplicityESDCorrected[esdCorrId]->SetMarkerColor(2);
+ //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(5);
// normalize without 0 bin
fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(2, fMultiplicityESDCorrected[esdCorrId]->GetNbinsX()));
Printf("Normalized without 0 bin!");
gPad->SetLogy();
TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
- legend->AddEntry(proj, "true distribution");
- legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
+ legend->AddEntry(proj, "True distribution");
+ legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "Unfolded distribution");
legend->SetFillColor(0);
+ legend->SetTextSize(0.04);
legend->Draw();
// unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
canvas1->cd(2);
+ canvas1->cd(2)->SetGridx();
+ canvas1->cd(2)->SetGridy();
+ canvas1->cd(2)->SetTopMargin(0.05);
+ canvas1->cd(2)->SetRightMargin(0.05);
+ canvas1->cd(2)->SetLeftMargin(0.12);
+ canvas1->cd(2)->SetBottomMargin(0.12);
gPad->SetLogy();
fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
//fCurrentESD->SetLineColor(2);
- fCurrentESD->SetTitle(";measured multiplicity;Entries");
+ fCurrentESD->SetTitle(Form(";Measured multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
fCurrentESD->SetStats(kFALSE);
+ fCurrentESD->GetYaxis()->SetTitleOffset(1.4);
+ //fCurrentESD->SetLabelSize(0.05, "xy");
+ //fCurrentESD->SetTitleSize(0.05, "xy");
fCurrentESD->DrawCopy("HIST E");
convolutedProj->SetLineColor(2);
+ convolutedProj->SetMarkerColor(2);
+ convolutedProj->SetMarkerStyle(5);
//proj2->SetMarkerColor(2);
//proj2->SetMarkerStyle(5);
- convolutedProj->DrawCopy("HIST SAME");
+ convolutedProj->DrawCopy("HIST SAME P");
legend = new TLegend(0.3, 0.8, 0.93, 0.93);
- legend->AddEntry(fCurrentESD, "measured distribution");
- legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
+ legend->AddEntry(fCurrentESD, "Measured distribution");
+ legend->AddEntry(convolutedProj, "R #otimes unfolded distribution", "P");
legend->SetFillColor(0);
+ legend->SetTextSize(0.04);
legend->Draw();
//TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
Double_t newChi2 = 0;
Double_t newChi2Limit150 = 0;
Int_t ndf = 0;
- for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
+ for (Int_t i=1; i<=diffMCUnfolded2->GetNbinsX(); ++i)
{
Double_t value = 0;
if (proj->GetBinError(i) > 0)
case kTrVtx : return fMultiplicityVtx[i]; break;
case kMB : return fMultiplicityMB[i]; break;
case kINEL : return fMultiplicityINEL[i]; break;
+ case kNSD : return fMultiplicityNSD[i]; break;
}
return 0;
}
+//____________________________________________________________________
+void AliMultiplicityCorrection::SetMultiplicityMC(Int_t i, EventType eventType, TH2F* hist)
+{
+ //
+ // returns the corresponding MC spectrum
+ //
+
+ switch (eventType)
+ {
+ case kTrVtx : fMultiplicityVtx[i] = hist; break;
+ case kMB : fMultiplicityMB[i] = hist; break;
+ case kINEL : fMultiplicityINEL[i] = hist; break;
+ case kNSD : fMultiplicityNSD[i] = hist; break;
+ }
+}
+
//____________________________________________________________________
TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
{
}
//____________________________________________________________________
-TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
+TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
{
//
// evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
// initialize seed with current time
gRandom->SetSeed(0);
+
+ if (methodType == AliUnfolding::kChi2Minimization)
+ AliUnfolding::SetCreateOverflowBin(5);
+ AliUnfolding::SetUnfoldingMethod(methodType);
const Int_t kErrorIterations = 150;
{
Printf("Iteration %d of %d...", n, kErrorIterations);
- // only done for chi2 minimization
- Bool_t createBigBin = kFALSE;
- if (methodType == kChi2Minimization)
- createBigBin = kTRUE;
-
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
TH1* measured = (TH1*) fCurrentESD->Clone("measured");
}
// only for bayesian method we have to do it before the call to Unfold...
- if (methodType == kBayesian)
+ if (methodType == AliUnfolding::kBayesian)
{
for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
{
{
result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
- Int_t returnCode = -1;
- if (methodType == kChi2Minimization)
- {
- returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
- }
- else if (methodType == kBayesian)
- returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
+ Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result);
if (returnCode != 0)
return 0;
for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
//new TCanvas; average->DrawClone();
-
+
// find cov. matrix
TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
}
- new TCanvas; covMatrix->DrawClone("COLZ");
+ TCanvas* c = new TCanvas; c->cd(); covMatrix->DrawCopy("COLZ");
// // fill 2D histogram that contains deviation from first
// TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
//new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
//new TCanvas; averageFirstRatio->DrawCopy();
+ static TH1* temp = 0;
+ if (!temp)
+ {
+ temp = (TH1*) standardDeviation->Clone("temp");
+ for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
+ temp->SetBinContent(x, temp->GetBinContent(x) * results[0]->GetBinContent(x));
+ }
+ else
+ {
+ // find difference from result[0] as TH2
+ TH2F* pulls = new TH2F("pulls", "pulls;multiplicity;difference", nBins, lowEdge, upEdge, 1000, -10, 10);
+ for (Int_t n=1; n<kErrorIterations; ++n)
+ for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
+ if (temp->GetBinContent(x) > 0)
+ pulls->Fill(results[n]->GetXaxis()->GetBinCenter(x), (results[0]->GetBinContent(x) - results[n]->GetBinContent(x)) / temp->GetBinContent(x));
+ new TCanvas("pulls", "pulls", 800, 600); pulls->DrawCopy(); pulls->FitSlicesY();
+ }
+
// clean up
for (Int_t n=0; n<kErrorIterations; ++n)
delete results[n];
// initialize seed with current time
gRandom->SetSeed(0);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
// normalize correction for given nPart
for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
- if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
+ AliUnfolding::SetBayesianParameters(regPar, nIterations);
+ AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
+ if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
return;
if (!determineError)
TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
result2->Reset();
- if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
+ if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0)
return;
result2->Scale(1.0 / result2->Integral());
delete error;
}
-//____________________________________________________________________
-Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
-{
- //
- // unfolds a spectrum
- //
-
- if (measured->Integral() <= 0)
- {
- Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
- return -1;
- }
-
- const Int_t kStartBin = 0;
-
- const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
- const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
-
- // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
- const Double_t kConvergenceLimit = kMaxT * 1e-6;
-
- // store information in arrays, to increase processing speed (~ factor 5)
- Double_t measuredCopy[kMaxM];
- Double_t measuredError[kMaxM];
- Double_t prior[kMaxT];
- Double_t result[kMaxT];
- Double_t efficiency[kMaxT];
-
- Double_t response[kMaxT][kMaxM];
- Double_t inverseResponse[kMaxT][kMaxM];
-
- // for normalization
- Float_t measuredIntegral = measured->Integral();
- for (Int_t m=0; m<kMaxM; m++)
- {
- measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
- measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
-
- for (Int_t t=0; t<kMaxT; t++)
- {
- response[t][m] = correlation->GetBinContent(t+1, m+1);
- inverseResponse[t][m] = 0;
- }
- }
-
- for (Int_t t=0; t<kMaxT; t++)
- {
- efficiency[t] = aEfficiency->GetBinContent(t+1);
- prior[t] = measuredCopy[t];
- result[t] = 0;
- }
-
- // pick prior distribution
- if (initialConditions)
- {
- printf("Using different starting conditions...\n");
- // for normalization
- Float_t inputDistIntegral = initialConditions->Integral();
- for (Int_t i=0; i<kMaxT; i++)
- prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
- }
-
- //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
-
- // unfold...
- for (Int_t i=0; i<nIterations || nIterations < 0; i++)
- {
- //printf(" iteration %i \n", i);
-
- // calculate IR from Beyes theorem
- // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
-
- Double_t chi2Measured = 0;
- for (Int_t m=0; m<kMaxM; m++)
- {
- Float_t norm = 0;
- for (Int_t t = kStartBin; t<kMaxT; t++)
- norm += response[t][m] * prior[t];
-
- // calc. chi2: (measured - response * prior) / error
- if (measuredError[m] > 0)
- {
- Double_t value = (measuredCopy[m] - norm) / measuredError[m];
- chi2Measured += value * value;
- }
-
- if (norm > 0)
- {
- for (Int_t t = kStartBin; t<kMaxT; t++)
- inverseResponse[t][m] = response[t][m] * prior[t] / norm;
- }
- else
- {
- for (Int_t t = kStartBin; t<kMaxT; t++)
- inverseResponse[t][m] = 0;
- }
- }
- //Printf("chi2Measured of the last prior is %e", chi2Measured);
-
- for (Int_t t = kStartBin; t<kMaxT; t++)
- {
- Float_t value = 0;
- for (Int_t m=0; m<kMaxM; m++)
- value += inverseResponse[t][m] * measuredCopy[m];
-
- if (efficiency[t] > 0)
- result[t] = value / efficiency[t];
- else
- result[t] = 0;
- }
-
- Double_t chi2LastIter = 0;
- // regularization (simple smoothing)
- for (Int_t t=kStartBin; t<kMaxT; t++)
- {
- Float_t newValue = 0;
- // 0 bin excluded from smoothing
- if (t > kStartBin+1 && t<kMaxT-1)
- {
- Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
-
- // weight the average with the regularization parameter
- newValue = (1 - regPar) * result[t] + regPar * average;
- }
- else
- newValue = result[t];
-
- // calculate chi2 (change from last iteration)
- if (prior[t] > 1e-5)
- {
- Double_t diff = (prior[t] - newValue) / prior[t];
- chi2LastIter += diff * diff;
- }
-
- prior[t] = newValue;
- }
- //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
- //convergence->Fill(i+1, chi2LastIter);
-
- if (nIterations < 0 && chi2LastIter < kConvergenceLimit)
- {
- Printf("Stopped Bayesian unfolding after %d iterations at chi2(change since last iteration) of %e; chi2Measured of the last prior is %e", i, chi2LastIter, chi2Measured);
- break;
- }
-
- //if (i % 5 == 0)
- // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
- } // end of iterations
-
- //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
- //delete convergence;
-
- for (Int_t t=0; t<kMaxT; t++)
- aResult->SetBinContent(t+1, result[t]);
-
- return 0;
-
- // ********
- // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
-
- /*printf("Calculating covariance matrix. This may take some time...\n");
-
- // TODO check if this is the right one...
- TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
-
- Int_t xBins = hInverseResponseBayes->GetNbinsX();
- Int_t yBins = hInverseResponseBayes->GetNbinsY();
-
- // calculate "unfolding matrix" Mij
- Float_t matrixM[251][251];
- for (Int_t i=1; i<=xBins; i++)
- {
- for (Int_t j=1; j<=yBins; j++)
- {
- if (fCurrentEfficiency->GetBinContent(i) > 0)
- matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
- else
- matrixM[i-1][j-1] = 0;
- }
- }
-
- Float_t* vectorn = new Float_t[yBins];
- for (Int_t j=1; j<=yBins; j++)
- vectorn[j-1] = fCurrentESD->GetBinContent(j);
-
- // first part of covariance matrix, depends on input distribution n(E)
- Float_t cov1[251][251];
-
- Float_t nEvents = fCurrentESD->Integral(); // N
-
- xBins = 20;
- yBins = 20;
-
- for (Int_t k=0; k<xBins; k++)
- {
- printf("In Cov1: %d\n", k);
- for (Int_t l=0; l<yBins; l++)
- {
- cov1[k][l] = 0;
-
- // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
- for (Int_t j=0; j<yBins; j++)
- cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
- * (1.0 - vectorn[j] / nEvents);
-
- // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
- for (Int_t i=0; i<yBins; i++)
- for (Int_t j=0; j<yBins; j++)
- {
- if (i == j)
- continue;
- cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
- * vectorn[j] / nEvents;
- }
- }
- }
-
- printf("Cov1 finished\n");
-
- TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
- cov->Reset();
-
- for (Int_t i=1; i<=xBins; i++)
- for (Int_t j=1; j<=yBins; j++)
- cov->SetBinContent(i, j, cov1[i-1][j-1]);
-
- new TCanvas;
- cov->Draw("COLZ");
-
- // second part of covariance matrix, depends on response matrix
- Float_t cov2[251][251];
-
- // Cov[P(Er|Cu), P(Es|Cu)] term
- Float_t covTerm[100][100][100];
- for (Int_t r=0; r<yBins; r++)
- for (Int_t u=0; u<xBins; u++)
- for (Int_t s=0; s<yBins; s++)
- {
- if (r == s)
- covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
- * (1.0 - hResponse->GetBinContent(u+1, r+1));
- else
- covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
- * hResponse->GetBinContent(u+1, s+1);
- }
-
- for (Int_t k=0; k<xBins; k++)
- for (Int_t l=0; l<yBins; l++)
- {
- cov2[k][l] = 0;
- printf("In Cov2: %d %d\n", k, l);
- for (Int_t i=0; i<yBins; i++)
- for (Int_t j=0; j<yBins; j++)
- {
- //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
- // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
- Float_t tmpCov = 0;
- for (Int_t r=0; r<yBins; r++)
- for (Int_t u=0; u<xBins; u++)
- for (Int_t s=0; s<yBins; s++)
- {
- if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
- || hResponse->GetBinContent(u+1, i+1) == 0)
- continue;
-
- tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
- * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
- * covTerm[r][u][s];
- }
-
- cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
- }
- }
-
- printf("Cov2 finished\n");
-
- for (Int_t i=1; i<=xBins; i++)
- for (Int_t j=1; j<=yBins; j++)
- cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
-
- new TCanvas;
- cov->Draw("COLZ");*/
-}
-
//____________________________________________________________________
Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
{
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+ SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
//normalize ESD
fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
- SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
+ SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
//normalize ESD
fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
if (id < 0 || id >= kESDHists)
return;
- TH2F* mc = fMultiplicityVtx[id];
-
+ // fill histogram used for random generation
+ TH1* tmp = fMultiplicityVtx[id]->ProjectionY("tmp");
+ tmp->Reset();
+
+ for (Int_t i=1; i<=tmp->GetNbinsX(); ++i)
+ tmp->SetBinContent(i, inputMC->Eval(tmp->GetXaxis()->GetBinCenter(i)) * tmp->GetXaxis()->GetBinWidth(i));
+
+ TH1* mcRnd = fMultiplicityVtx[id]->ProjectionY("mcRnd");
+ mcRnd->Reset();
+ mcRnd->FillRandom(tmp, tmp->Integral());
+
+ //new TCanvas; tmp->Draw();
+ //new TCanvas; mcRnd->Draw();
+
+ // and move into 2d histogram
+ TH1* mc = fMultiplicityVtx[id];
mc->Reset();
-
for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
{
- mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
- mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
+ mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, mcRnd->GetBinContent(i));
+ mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, TMath::Sqrt(mcRnd->GetBinContent(i)));
+ }
+
+ //new TCanvas; mc->Draw("COLZ");
+
+ // now randomize the measured histogram; funcMeasured is used as pilot function to generated the measured entries
+ TH1* funcMeasured = CalculateMultiplicityESD(tmp, id)->ProjectionY("funcMeasured");
+
+ //new TCanvas; funcMeasured->Draw();
+
+ fMultiplicityESD[id]->Reset();
+
+ TH1* measRnd = fMultiplicityESD[id]->ProjectionY("measRnd");
+ measRnd->FillRandom(funcMeasured, tmp->Integral());
+
+ //new TCanvas; measRnd->Draw();
+
+ fMultiplicityESD[id]->Reset();
+ for (Int_t i=1; i<=fMultiplicityESD[id]->GetNbinsY(); ++i)
+ {
+ fMultiplicityESD[id]->SetBinContent(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, measRnd->GetBinContent(i));
+ fMultiplicityESD[id]->SetBinError(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, TMath::Sqrt(measRnd->GetBinContent(i)));
}
-
- //new TCanvas;
- //mc->Draw("COLZ");
-
- TH1* proj = mc->ProjectionY();
-
- TString tmp(fMultiplicityESD[id]->GetName());
- delete fMultiplicityESD[id];
- fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
- fMultiplicityESD[id]->SetName(tmp);
}