coding conventions
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Jul 2007 14:25:59 +0000 (14:25 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Jul 2007 14:25:59 +0000 (14:25 +0000)
plotsMultiplicity macro can now compile

PWG0/dNdEta/AliMultiplicityCorrection.cxx
PWG0/dNdEta/AliMultiplicityCorrection.h
PWG0/dNdEta/AliMultiplicityMCSelector.cxx
PWG0/dNdEta/AliMultiplicityMCSelector.h
PWG0/dNdEta/plotsMultiplicity.C
PWG0/dNdEta/runMultiplicitySelector.C

index e25996e..46aee66 100644 (file)
@@ -18,6 +18,7 @@
 // This class is used to store correction maps, raw input and results of the multiplicity
 // measurement with the ITS or TPC
 // It also contains functions to correct the spectrum using different methods.
+// e.g. chi2 minimization and bayesian unfolding
 //
 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
 
 
 ClassImp(AliMultiplicityCorrection)
 
-const Int_t AliMultiplicityCorrection::fgMaxInput  = 250;  // bins in measured histogram
-const Int_t AliMultiplicityCorrection::fgMaxParams = 250;  // bins in unfolded histogram = number of fit params
+const Int_t AliMultiplicityCorrection::fgkMaxInput  = 250;  // bins in measured histogram
+const Int_t AliMultiplicityCorrection::fgkMaxParams = 250;  // bins in unfolded histogram = number of fit params
 
-TH1* AliMultiplicityCorrection::fCurrentESD = 0;
-TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
-TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
-TMatrixD* AliMultiplicityCorrection::fCorrelationMatrix = 0;
-TMatrixD* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
-TVectorD* AliMultiplicityCorrection::fCurrentESDVector = 0;
-TVectorD* AliMultiplicityCorrection::fEntropyAPriori = 0;
-AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
-Float_t AliMultiplicityCorrection::fRegularizationWeight = 5000;
-TF1* AliMultiplicityCorrection::fNBD = 0;
+TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0;
+TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0;
+TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0;
+TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0;
+
+AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1;
+Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000;
+
+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
@@ -82,6 +85,8 @@ void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
 
     fgQualityRegionsB[2] = 180;
     fgQualityRegionsE[2] = 210;
+
+    Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
   }
   else
   {
@@ -94,12 +99,14 @@ void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
 
     fgQualityRegionsB[2] = 88;
     fgQualityRegionsE[2] = 108;
+
+    Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
   }
 }
 
 //____________________________________________________________________
 AliMultiplicityCorrection::AliMultiplicityCorrection() :
-  TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
+  TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
 {
   //
   // default constructor
@@ -120,6 +127,9 @@ AliMultiplicityCorrection::AliMultiplicityCorrection() :
     fCorrelation[i] = 0;
     fMultiplicityESDCorrected[i] = 0;
   }
+
+  for (Int_t i = 0; i < kQualityRegions; ++i)
+    fQuality[i] = 0;
 }
 
 //____________________________________________________________________
@@ -150,7 +160,7 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C
                           //1000.5 };
 
   #define VTXBINNING 10, binLimitsVtx
-  #define NBINNING fgMaxParams, binLimitsN*/
+  #define NBINNING fgkMaxParams, binLimitsN*/
 
   #define NBINNING 501, -0.5, 500.5
   #define VTXBINNING 1, -10, 10
@@ -185,6 +195,39 @@ AliMultiplicityCorrection::~AliMultiplicityCorrection()
   //
   // Destructor
   //
+
+  for (Int_t i = 0; i < kESDHists; ++i)
+  {
+    if (fMultiplicityESD[i])
+      delete fMultiplicityESD[i];
+    fMultiplicityESD[i] = 0;
+  }
+
+  for (Int_t i = 0; i < kMCHists; ++i)
+  {
+    if (fMultiplicityVtx[i])
+      delete fMultiplicityVtx[i];
+    fMultiplicityVtx[i] = 0;
+
+    if (fMultiplicityMB[i])
+      delete fMultiplicityMB[i];
+    fMultiplicityMB[i] = 0;
+
+    if (fMultiplicityINEL[i])
+      delete fMultiplicityINEL[i];
+    fMultiplicityINEL[i] = 0;
+  }
+
+  for (Int_t i = 0; i < kCorrHists; ++i)
+  {
+    if (fCorrelation[i])
+      delete fCorrelation[i];
+    fCorrelation[i] = 0;
+
+    if (fMultiplicityESDCorrected[i])
+      delete fMultiplicityESDCorrected[i];
+    fMultiplicityESDCorrected[i] = 0;
+  }
 }
 
 //____________________________________________________________________
@@ -410,7 +453,7 @@ Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
   Double_t chi2 = 0;
 
   // ignore the first bin here. on purpose...
-  for (Int_t i=2; i<fgMaxParams; ++i)
+  for (Int_t i=2; i<fgkMaxParams; ++i)
   {
     Double_t right  = params[i];
     Double_t left   = params[i-1];
@@ -435,7 +478,7 @@ Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
   Double_t chi2 = 0;
 
   // ignore the first bin here. on purpose...
-  for (Int_t i=2+1; i<fgMaxParams; ++i)
+  for (Int_t i=2+1; i<fgkMaxParams; ++i)
   {
     if (params[i-1] == 0)
       continue;
@@ -464,12 +507,12 @@ Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
 
   Double_t chi2 = 0;
 
-  /*Float_t der[fgMaxParams];
+  /*Float_t der[fgkMaxParams];
 
-  for (Int_t i=0; i<fgMaxParams-1; ++i)
+  for (Int_t i=0; i<fgkMaxParams-1; ++i)
     der[i] = params[i+1] - params[i];
 
-  for (Int_t i=0; i<fgMaxParams-6; ++i)
+  for (Int_t i=0; i<fgkMaxParams-6; ++i)
   {
     for (Int_t j=1; j<=5; ++j)
     {
@@ -479,7 +522,7 @@ Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
   }*/
 
   // ignore the first bin here. on purpose...
-  for (Int_t i=2+1; i<fgMaxParams; ++i)
+  for (Int_t i=2+1; i<fgkMaxParams; ++i)
   {
     if (params[i-1] == 0)
       continue;
@@ -510,7 +553,7 @@ Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& param
   Double_t chi2 = 0;
 
   // ignore the first bin here. on purpose...
-  for (Int_t i=3; i<fgMaxParams; ++i)
+  for (Int_t i=3; i<fgkMaxParams; ++i)
   {
     Double_t right  = params[i];
     Double_t middle = params[i-1];
@@ -537,16 +580,16 @@ Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
 
   Double_t paramSum = 0;
   // ignore the first bin here. on purpose...
-  for (Int_t i=1; i<fgMaxParams; ++i)
+  for (Int_t i=1; i<fgkMaxParams; ++i)
     paramSum += params[i];
 
   Double_t chi2 = 0;
-  for (Int_t i=1; i<fgMaxParams; ++i)
+  for (Int_t i=1; i<fgkMaxParams; ++i)
   {
     //Double_t tmp = params[i] / paramSum;
     Double_t tmp = params[i];
-    if (tmp > 0 && (*fEntropyAPriori)[i] > 0)
-      chi2 += tmp * TMath::Log(tmp / (*fEntropyAPriori)[i]);
+    if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
+      chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
   }
 
   return 10.0 + chi2;
@@ -565,12 +608,12 @@ void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Dou
   // func->SetParLimits(2, 0.001, 1000);
   //
 
-  fNBD->SetParameters(params[0], params[1], params[2]);
+  fgNBD->SetParameters(params[0], params[1], params[2]);
 
-  Double_t params2[fgMaxParams];
+  Double_t params2[fgkMaxParams];
 
-  for (Int_t i=0; i<fgMaxParams; ++i)
-    params2[i] = fNBD->Eval(i);
+  for (Int_t i=0; i<fgkMaxParams; ++i)
+    params2[i] = fgNBD->Eval(i);
 
   MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
 
@@ -586,13 +629,13 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c
   //
 
   // d
-  static TVectorD paramsVector(fgMaxParams);
-  for (Int_t i=0; i<fgMaxParams; ++i)
+  static TVectorD paramsVector(fgkMaxParams);
+  for (Int_t i=0; i<fgkMaxParams; ++i)
     paramsVector[i] = params[i] * params[i];
 
   // calculate penalty factor
   Double_t penaltyVal = 0;
-  switch (fRegularizationType)
+  switch (fgRegularizationType)
   {
     case kNone:       break;
     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
@@ -605,30 +648,30 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c
   //if (penaltyVal > 1e10)
   //  paramsVector2.Print();
 
-  penaltyVal *= fRegularizationWeight;
+  penaltyVal *= fgRegularizationWeight;
 
   // Ad
-  TVectorD measGuessVector(fgMaxInput);
-  measGuessVector = (*fCorrelationMatrix) * paramsVector;
+  TVectorD measGuessVector(fgkMaxInput);
+  measGuessVector = (*fgCorrelationMatrix) * paramsVector;
 
   // Ad - m
-  measGuessVector -= (*fCurrentESDVector);
+  measGuessVector -= (*fgCurrentESDVector);
 
   TVectorD copy(measGuessVector);
 
   // (Ad - m) W
-  // this step can be optimized because currently only the diagonal elements of fCorrelationCovarianceMatrix are used
+  // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
   // normal way is like this:
-  // measGuessVector *= (*fCorrelationCovarianceMatrix);
+  // measGuessVector *= (*fgCorrelationCovarianceMatrix);
   // optimized way like this:
-  for (Int_t i=0; i<fgMaxParams; ++i)
-    measGuessVector[i] *= (*fCorrelationCovarianceMatrix)(i, i);
+  for (Int_t i=0; i<fgkMaxParams; ++i)
+    measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
 
   //measGuessVector.Print();
 
   // (Ad - m) W (Ad - m)
   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
-  // big number ((*fCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
+  // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
   Double_t chi2FromFit = measGuessVector * copy * 1e6;
 
   chi2 = chi2FromFit + penaltyVal;
@@ -791,74 +834,93 @@ TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventT
 //____________________________________________________________________
 void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight)
 {
-  fRegularizationType = type;
-  fRegularizationWeight = weight;
+  //
+  // 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);
 }
 
 //____________________________________________________________________
-Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
+void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
 {
   //
-  // correct spectrum using minuit chi2 method
-  //
-  // if check is kTRUE the input MC solution (by definition the right solution) is used
-  // no fit is made and just the chi2 is calculated
+  // sets the parameters for Bayesian unfolding
   //
 
-  Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
-  Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
+  fgBayesianSmoothing = smoothing;
+  fgBayesianIterations = nIterations;
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kTRUE);
-  //normalize ESD
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+  printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
+}
 
-  fCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
-  fCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
-  fCurrentESDVector = new TVectorD(fgMaxInput);
-  fEntropyAPriori = new TVectorD(fgMaxParams);
+//____________________________________________________________________
+Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
+{
+  //
+  // 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
+  //
+  // returns minuit status (0 = success), (-1 when check was set)
+  //
 
-  /*new TCanvas; fCurrentESD->DrawCopy();
-  fCurrentESD = ((TH2*) fCurrentCorrelation)->ProjectionY("check-proj2");
-  fCurrentESD->Sumw2();
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
-  fCurrentESD->SetLineColor(2);
-  fCurrentESD->DrawCopy("SAME");*/
+  //normalize measured
+  measured->Scale(1.0 / measured->Integral());
+
+  if (!fgCorrelationMatrix)
+    fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgkMaxParams);
+  if (!fgCorrelationCovarianceMatrix)
+    fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
+  if (!fgCurrentESDVector)
+    fgCurrentESDVector = new TVectorD(fgkMaxInput);
+  if (!fgEntropyAPriori)
+    fgEntropyAPriori = new TVectorD(fgkMaxParams);
+
+  fgCorrelationMatrix->Zero();
+  fgCorrelationCovarianceMatrix->Zero();
+  fgCurrentESDVector->Zero();
+  fgEntropyAPriori->Zero();
 
   // normalize correction for given nPart
-  for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+  for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
   {
-    Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
+    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<=fCurrentCorrelation->GetNbinsY(); ++j)
+    for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
     {
       // find most probably value
-      if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
+      if (maxValue < correlation->GetBinContent(i, j))
       {
-        maxValue = fCurrentCorrelation->GetBinContent(i, j);
+        maxValue = correlation->GetBinContent(i, j);
         maxBin = j;
       }
 
       // npart sum to 1
-      fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
-      fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
+      correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
+      correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
 
-      if (i <= fgMaxParams && j <= fgMaxInput)
-        (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
+      if (i <= fgkMaxParams && 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, fgMaxParams);
+  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgkMaxParams);
   Double_t arglist[100];
 
-  // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find out why...
+  // 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);
 
@@ -868,38 +930,26 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha
   // set minimization function
   minuit->SetFCN(MinuitFitFunction);
 
-  for (Int_t i=0; i<fgMaxParams; i++)
-    (*fEntropyAPriori)[i] = 1;
+  for (Int_t i=0; i<fgkMaxParams; i++)
+    (*fgEntropyAPriori)[i] = 1;
 
-  if (inputDist)
+  if (initialConditions)
   {
     printf("Using different starting conditions...\n");
-    new TCanvas;
-    inputDist->DrawCopy();
+    //new TCanvas; initialConditions->DrawCopy();
 
-    inputDist->Scale(1.0 / inputDist->Integral());
-    for (Int_t i=0; i<fgMaxParams; i++)
-      if (inputDist->GetBinContent(i+1) > 0)
-        (*fEntropyAPriori)[i] = inputDist->GetBinContent(i+1);
+    initialConditions->Scale(1.0 / initialConditions->Integral());
+    for (Int_t i=0; i<fgkMaxParams; i++)
+      if (initialConditions->GetBinContent(i+1) > 0)
+        (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
   }
   else
-    inputDist = fCurrentESD;
-
-
-  //Float_t minStartValue = 0; //1e-3;
+    initialConditions = measured;
 
-  //new TCanvas; fMultiplicityVtx[mcTarget]->Draw("COLZ");
-  TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj");
-  //proj->Rebin(2);
-  proj->Scale(1.0 / proj->Integral());
-
-  Double_t results[fgMaxParams+1];
-  for (Int_t i=0; i<fgMaxParams; ++i)
+  Double_t results[fgkMaxParams+1];
+  for (Int_t i=0; i<fgkMaxParams; ++i)
   {
-    results[i] = inputDist->GetBinContent(i+1);
-
-    if (check)
-      results[i] = proj->GetBinContent(i+1);
+    results[i] = initialConditions->GetBinContent(i+1);
 
     // minimum value
     results[i] = TMath::Max(results[i], 1e-3);
@@ -912,24 +962,20 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha
 
     minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
   }
-  //minuit->SetParameter(fgMaxParams, "alpha", 1e4, 1, 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<fgMaxInput; ++i)
+  for (Int_t i=0; i<fgkMaxInput; ++i)
   {
-    //printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
-
-    (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
-    if (fCurrentESD->GetBinError(i+1) > 0)
-      (*fCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / fCurrentESD->GetBinError(i+1) / fCurrentESD->GetBinError(i+1);
+    (*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 ((*fCorrelationCovarianceMatrix)(i, i) > 1e7)
-      (*fCorrelationCovarianceMatrix)(i, i) = 0;
-
-    //printf("%d --> %e\n", i, (*fCorrelationCovarianceMatrix)(i, i));
+    if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
+      (*fgCorrelationCovarianceMatrix)(i, i) = 0;
   }
 
   Int_t dummy;
@@ -951,18 +997,18 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha
   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
   //minuit->ExecuteCommand("MINOS", arglist, 0);
 
-  for (Int_t i=0; i<fgMaxParams; ++i)
+  for (Int_t i=0; i<fgkMaxParams; ++i)
   {
-    if (fCurrentEfficiency->GetBinContent(i+1) > 0)
+    if (aEfficiency->GetBinContent(i+1) > 0)
     {
-      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
+      result->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / aEfficiency->GetBinContent(i+1));
       // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
-      fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) /  fCurrentEfficiency->GetBinContent(i+1));
+      result->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) /  aEfficiency->GetBinContent(i+1));
     }
     else
     {
-      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, 0);
-      fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
+      result->SetBinContent(i+1, 0);
+      result->SetBinError(i+1, 0);
     }
   }
 
@@ -970,6 +1016,21 @@ Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPha
 }
 
 //____________________________________________________________________
+Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
+{
+  //
+  // correct spectrum using minuit chi2 method
+  //
+  // for description of parameters, see UnfoldWithMinuit
+  //
+
+  Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kTRUE);
+
+  return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
+}
+
+//____________________________________________________________________
 void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
 {
   //
@@ -982,9 +1043,9 @@ void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSp
   //normalize ESD
   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 
-  fCorrelationMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
-  fCorrelationCovarianceMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
-  fCurrentESDVector = new TVectorD(fgMaxParams);
+  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)
@@ -998,23 +1059,23 @@ void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSp
       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
 
-      if (i <= fgMaxParams && j <= fgMaxParams)
-        (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
+      if (i <= fgkMaxParams && j <= fgkMaxParams)
+        (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
     }
   }
 
-  for (Int_t i=0; i<fgMaxParams; ++i)
+  for (Int_t i=0; i<fgkMaxParams; ++i)
   {
-    (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
+    (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
     if (fCurrentESD->GetBinError(i+1) > 0)
-      (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
+      (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
   }
 
   // Create NBD function
-  if (!fNBD)
+  if (!fgNBD)
   {
-    fNBD = 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);
-    fNBD->SetParNames("scaling", "averagen", "k");
+    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
@@ -1031,14 +1092,14 @@ void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSp
   minuit->ExecuteCommand("MIGRAD", arglist, 0);
   //minuit->ExecuteCommand("MINOS", arglist, 0);
 
-  fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
+  fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
 
-  for (Int_t i=0; i<fgMaxParams; ++i)
+  for (Int_t i=0; i<fgkMaxParams; ++i)
   {
-    printf("%d : %f\n", i, fNBD->Eval(i));
-    if (fNBD->Eval(i) > 0)
+    printf("%d : %f\n", i, fgNBD->Eval(i));
+    if (fgNBD->Eval(i) > 0)
     {
-      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fNBD->Eval(i));
+      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
     }
   }
@@ -1077,6 +1138,10 @@ void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
 //____________________________________________________________________
 void AliMultiplicityCorrection::DrawHistograms()
 {
+  //
+  // draws the histograms of this class
+  //
+
   printf("ESD:\n");
 
   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
@@ -1185,8 +1250,6 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
   else
     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
 
-  //NormalizeToBinWidth(proj2);
-
   TH1* residual = (TH1*) convolutedProj->Clone("residual");
   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
 
@@ -1240,10 +1303,6 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
 
   canvas1->cd(1);
   TH1* proj = (TH1*) mcHist->Clone("proj");
-  NormalizeToBinWidth(proj);
-
-  if (normalizeESD)
-    NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
 
   proj->GetXaxis()->SetRangeUser(0, 200);
   proj->SetTitle(";true multiplicity;Entries");
@@ -1510,9 +1569,9 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
 
 
     Double_t newChi2 = 0;
-    Double_t newChi2_150 = 0;
+    Double_t newChi2Limit150 = 0;
     Int_t ndf = 0;
-    for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgMaxParams+1); ++i)
+    for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
     {
       Double_t value = 0;
       if (proj->GetBinError(i) > 0)
@@ -1520,7 +1579,7 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
         newChi2 += value * value;
         if (i > 1 && i <= mcBinLimit)
-          newChi2_150 += value * value;
+          newChi2Limit150 += value * value;
         ++ndf;
 
         for (Int_t region=0; region<kQualityRegions; ++region)
@@ -1545,8 +1604,8 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     if (mcBinLimit > 1)
     {
       // TODO another FAKE
-      fLastChi2MC = newChi2_150 / (mcBinLimit - 1);
-      Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2_150, mcBinLimit - 1, fLastChi2MC);
+      fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
+      Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
     }
     else
       fLastChi2MC = -1;
@@ -1657,7 +1716,7 @@ void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage)
+void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
 {
   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
   // These values are computed during DrawComparison, thus this function picks up the
@@ -1691,7 +1750,7 @@ TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
 }
 
 //____________________________________________________________________
-TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType,  Bool_t randomizeMeasured, Bool_t randomizeResponse, Float_t regPar, Int_t nIterations, TH1* compareTo)
+TH1* AliMultiplicityCorrection::StatisticalUncertainty(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
@@ -1703,6 +1762,8 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
   // returns the error assigned to the measurement
   //
 
+  Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
+
   // initialize seed with current time
   gRandom->SetSeed(0);
 
@@ -1717,7 +1778,12 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
   {
     Printf("Iteration %d of %d...", n, kErrorIterations);
 
-    SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
+    // only done for chi2 minimization
+    Bool_t createBigBin = kFALSE;
+    if (methodType == kChi2Minimization)
+      createBigBin = kTRUE;
+
+    SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
 
     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
 
@@ -1743,28 +1809,32 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
       }
     }
 
-    for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+    // only for bayesian method we have to do it before the call to Unfold...
+    if (methodType == kBayesian)
     {
-      // with this it is normalized to 1
-      Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
-
-      // with this normalized to the given efficiency
-      if (fCurrentEfficiency->GetBinContent(i) > 0)
-        sum /= fCurrentEfficiency->GetBinContent(i);
-      else
-        sum = 0;
-
-      for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+      for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
       {
-        if (sum > 0)
-        {
-          fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
-          fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
-        }
+        // with this it is normalized to 1
+        Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
+
+        // with this normalized to the given efficiency
+        if (fCurrentEfficiency->GetBinContent(i) > 0)
+          sum /= fCurrentEfficiency->GetBinContent(i);
         else
+          sum = 0;
+
+        for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
         {
-          fCurrentCorrelation->SetBinContent(i, j, 0);
-          fCurrentCorrelation->SetBinError(i, j, 0);
+          if (sum > 0)
+          {
+            fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
+            fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
+          }
+          else
+          {
+            fCurrentCorrelation->SetBinContent(i, j, 0);
+            fCurrentCorrelation->SetBinError(i, j, 0);
+          }
         }
       }
     }
@@ -1777,10 +1847,20 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
       result->Sumw2();
     }
     else
-      result = UnfoldWithBayesian(measured, regPar, nIterations, 0);
+    {
+      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);
 
-    if (!result)
-      return 0;
+      if (returnCode != 0)
+        return 0;
+    }
 
     // normalize
     result->Scale(1.0 / result->Integral());
@@ -1836,24 +1916,49 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
       }
   new TCanvas; covMatrix->DrawClone("COLZ");
 
-  // fill 2D histogram that contains deviation from first
-  TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
-  for (Int_t n=1; n<kErrorIterations; ++n)
-    for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
-      deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
-  new TCanvas; deviations->DrawCopy("COLZ");
-
+//   // fill 2D histogram that contains deviation from first
+//   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
+//   for (Int_t n=1; n<kErrorIterations; ++n)
+//     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
+//       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
+//   //new TCanvas; deviations->DrawCopy("COLZ");
+// 
+//   // get standard deviation "by hand"
+//   for (Int_t x=1; x<=nBins; x++)
+//   {
+//     if (results[0]->GetBinContent(x) > 0)
+//     {
+//       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
+//       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
+//       //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
+//       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
+//     }
+//   }
+
+  // calculate standard deviation of (results[0] - results[n])
   TH1* standardDeviation = (TH1*) maxError->Clone("standardDeviation");
   standardDeviation->Reset();
 
-  // get standard deviation "by hand"
-  for (Int_t x=1; x<=nBins; x++)
+  for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
   {
     if (results[0]->GetBinContent(x) > 0)
     {
-      TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
-      Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
+      Double_t average = 0;
+      for (Int_t n=1; n<kErrorIterations; ++n)
+        average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
+      average /= kErrorIterations-1;
+
+      Double_t variance = 0;
+      for (Int_t n=1; n<kErrorIterations; ++n)
+      {
+        Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
+        variance += value * value;
+      }
+      variance /= kErrorIterations-1;
+
+      Double_t standardDev = TMath::Sqrt(variance);
       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
+      Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
     }
   }
 
@@ -1876,15 +1981,15 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
       if (average->GetBinContent(x) > 0)
         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
 
-  new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
+  //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
 
   // plot difference between average and MC/first result
   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
   averageFirstRatio->Reset();
   averageFirstRatio->Divide(results[0], average);
 
-  new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
-  new TCanvas; averageFirstRatio->DrawCopy();
+  //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
+  //new TCanvas; averageFirstRatio->DrawCopy();
 
   // clean up
   for (Int_t n=0; n<kErrorIterations; ++n)
@@ -1892,8 +1997,6 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
   delete[] results;
 
   // fill into result histogram
-  Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
-
   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
 
@@ -1904,7 +2007,7 @@ TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist, Bool_t determineError)
+void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
 {
   //
   // correct spectrum using bayesian method
@@ -1944,13 +2047,9 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 
-  TH1* result = UnfoldWithBayesian(fCurrentESD, regPar, nIterations, inputDist);
-  if (!result)
+  if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
     return;
 
-  for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
-    fMultiplicityESDCorrected[correlationID]->SetBinContent(i, result->GetBinContent(i));
-
   if (!determineError)
   {
     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
@@ -1962,10 +2061,10 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
   // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
   // error of the unfolding method itself.
 
-  TH1* maxError = (TH1*) result->Clone("maxError");
+  TH1* maxError = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("maxError");
   maxError->Reset();
 
-  TH1* normalizedResult = (TH1*) result->Clone("normalizedResult");
+  TH1* normalizedResult = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("normalizedResult");
   normalizedResult->Scale(1.0 / normalizedResult->Integral());
 
   const Int_t kErrorIterations = 20;
@@ -1973,6 +2072,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
 
   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
+  TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->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
@@ -1984,25 +2084,23 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
       randomized->SetBinError(x, TMath::Sqrt(randomValue));
     }
 
-    TH1* result2 = UnfoldWithBayesian(randomized, regPar, nIterations, inputDist);
-    if (!result2)
+    result2->Reset();
+    if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
       return;
 
     result2->Scale(1.0 / result2->Integral());
 
     // calculate ratio
-    TH1* ratio = (TH1*) result2->Clone("ratio");
-    ratio->Divide(normalizedResult);
+    result2->Divide(normalizedResult);
 
     //new TCanvas; ratio->DrawCopy("HIST");
 
     // find max. deviation
-    for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
-      maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
-
-    delete ratio;
-    delete result2;
+    for (Int_t x=1; x<=result2->GetNbinsX(); x++)
+      maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - result2->GetBinContent(x))));
   }
+  delete randomized;
+  delete result2;
 
   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
@@ -2012,7 +2110,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 }
 
 //____________________________________________________________________
-TH1* AliMultiplicityCorrection::UnfoldWithBayesian(TH1* measured, Float_t regPar, Int_t nIterations, TH1* inputDist)
+Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
 {
   //
   // unfolds a spectrum
@@ -2021,13 +2119,13 @@ TH1* AliMultiplicityCorrection::UnfoldWithBayesian(TH1* measured, Float_t regPar
   if (measured->Integral() <= 0)
   {
     Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
-    return 0;
+    return -1;
   }
 
   const Int_t kStartBin = 0;
 
-  const Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
-  const Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
+  const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
+  const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
 
   // store information in arrays, to increase processing speed (~ factor 5)
   Double_t measuredCopy[kMaxM];
@@ -2046,26 +2144,26 @@ TH1* AliMultiplicityCorrection::UnfoldWithBayesian(TH1* measured, Float_t regPar
 
     for (Int_t t=0; t<kMaxT; t++)
     {
-      response[t][m] = fCurrentCorrelation->GetBinContent(t+1, m+1);
+      response[t][m] = correlation->GetBinContent(t+1, m+1);
       inverseResponse[t][m] = 0;
     }
   }
 
   for (Int_t t=0; t<kMaxT; t++)
   {
-    efficiency[t] = fCurrentEfficiency->GetBinContent(t+1);
+    efficiency[t] = aEfficiency->GetBinContent(t+1);
     prior[t] = measuredCopy[t];
     result[t] = 0;
   }
 
   // pick prior distribution
-  if (inputDist)
+  if (initialConditions)
   {
     printf("Using different starting conditions...\n");
     // for normalization
-    Float_t inputDistIntegral = inputDist->Integral();
+    Float_t inputDistIntegral = initialConditions->Integral();
     for (Int_t i=0; i<kMaxT; i++)
-      prior[i] = inputDist->GetBinContent(i+1) / inputDistIntegral;
+      prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
   }
 
   //TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
@@ -2151,13 +2249,10 @@ TH1* AliMultiplicityCorrection::UnfoldWithBayesian(TH1* measured, Float_t regPar
 
   //delete convergence;
 
-  // TODO this does not work when the number of bins is differnt
-  TH1* resultHist = (TH1*) measured->Clone("resultHist");
-  resultHist->Reset();
   for (Int_t t=0; t<kMaxT; t++)
-    resultHist->SetBinContent(t+1, result[t]);
+    aResult->SetBinContent(t+1, result[t]);
 
-  return resultHist;
+  return 0;
 
   // ********
   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
@@ -2447,8 +2542,6 @@ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t ful
   //normalize ESD
   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
 
-  NormalizeToBinWidth((TH2*) fCurrentCorrelation);
-
   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
   correction->SetTitle("GaussianMean");
 
index 3da2d88..dc6f2f1 100644 (file)
@@ -8,6 +8,7 @@
 //
 // class that contains the correction matrix and the functions for
 // correction the multiplicity spectrum
+// implements a several unfolding methods: e.g. chi2 minimization and bayesian unfolding
 //
 
 class TH1;
@@ -18,6 +19,11 @@ class TH3F;
 class TF1;
 class TCollection;
 
+// defined here, because it does not seem possible to predeclare these (or i do not know how)
+// -->
+// $ROOTSYS/include/TVectorDfwd.h:21: conflicting types for `typedef struct TVectorT<Double_t> TVectorD'
+// PWG0/dNdEta/AliMultiplicityCorrection.h:21: previous declaration as `struct TVectorD'
+
 #include <TMatrixD.h>
 #include <TVectorD.h>
 
@@ -25,6 +31,7 @@ class AliMultiplicityCorrection : public TNamed {
   public:
     enum EventType { kTrVtx = 0, kMB, kINEL };
     enum RegularizationType { kNone = 0, kPol0, kPol1, kLog, kEntropy, kCurvature };
+    enum MethodType { kChi2Minimization = 0, kBayesian = 1 };
     enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8, kQualityRegions = 3 };
 
     AliMultiplicityCorrection();
@@ -43,13 +50,15 @@ class AliMultiplicityCorrection : public TNamed {
     void DrawHistograms();
     void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE);
 
-    Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0);
+    Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* initialConditions = 0);
     void SetRegularizationParameters(RegularizationType type, Float_t weight);
+    void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
 
     void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
 
-    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* inputDist = 0, Bool_t determineError = kTRUE);
-    TH1* BayesianStatisticsEffect(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, Float_t regPar = 1, Int_t nIterations = 100, TH1* compareTo = 0);
+    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 1, Int_t nIterations = 100, TH1* initialConditions = 0, Bool_t determineError = kTRUE);
+
+    TH1* StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo = 0);
 
     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
 
@@ -76,18 +85,18 @@ class AliMultiplicityCorrection : public TNamed {
     static void NormalizeToBinWidth(TH1* hist);
     static void NormalizeToBinWidth(TH2* hist);
 
-    void GetComparisonResults(Float_t* mc = 0, Int_t* mcLimit = 0, Float_t* residuals = 0, Float_t* ratioAverage = 0);
+    void GetComparisonResults(Float_t* mc = 0, Int_t* mcLimit = 0, Float_t* residuals = 0, Float_t* ratioAverage = 0) const;
 
     TH1* GetEfficiency(Int_t inputRange, EventType eventType);
 
     static void SetQualityRegions(Bool_t SPDStudy);
-    Float_t GetQuality(Int_t region) { return fQuality[region]; }
+    Float_t GetQuality(Int_t region) const { return fQuality[region]; }
 
     void FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y);
 
   protected:
-    static const Int_t fgMaxParams;  // bins in unfolded histogram = number of fit params
-    static const Int_t fgMaxInput;   // bins in measured histogram
+    static const Int_t fgkMaxParams;  //! bins in unfolded histogram = number of fit params
+    static const Int_t fgkMaxInput;   //! bins in measured histogram
 
     static Double_t RegularizationPol0(TVectorD& params);
     static Double_t RegularizationPol1(TVectorD& params);
@@ -101,44 +110,51 @@ class AliMultiplicityCorrection : public TNamed {
     void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
 
     Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u);
-    TH1* UnfoldWithBayesian(TH1* measured, Float_t regPar, Int_t nIterations, TH1* inputDist);
+    static Int_t UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations);
+    static Int_t UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check);
+
+    TH1* fCurrentESD;         //! static variable to be accessed by MINUIT
+    TH1* fCurrentCorrelation; //! static variable to be accessed by MINUIT
+    TH1* fCurrentEfficiency;  //! static variable to be accessed by MINUIT
 
-    static TH1* fCurrentESD;         // static variable to be accessed by MINUIT
-    static TH1* fCurrentCorrelation; // static variable to be accessed by MINUIT
-    static TH1* fCurrentEfficiency;  // static variable to be accessed by MINUIT
+    // static variable to be accessed by MINUIT
+    static TMatrixD* fgCorrelationMatrix;            //! contains fCurrentCorrelation in matrix form
+    static TMatrixD* fgCorrelationCovarianceMatrix;  //! contains the errors of fCurrentESD
+    static TVectorD* fgCurrentESDVector;             //! contains fCurrentESD
+    static TVectorD* fgEntropyAPriori;               //! a-priori distribution for entropy regularization
 
-    static TMatrixD* fCorrelationMatrix;            // contains fCurrentCorrelation in matrix form
-    static TMatrixD* fCorrelationCovarianceMatrix;  // contains the errors of fCurrentESD
-    static TVectorD* fCurrentESDVector;             // contains fCurrentESD
-    static TVectorD* fEntropyAPriori;               // a-priori distribution for entropy regularization
+    static TF1* fgNBD;   //! negative binomial distribution
 
-    static TF1* fNBD;   // negative binomial distribution
+    static RegularizationType fgRegularizationType; //! regularization that is used during Chi2 method
+    static Float_t fgRegularizationWeight;          //! factor for regularization term
 
-    static RegularizationType fRegularizationType; // regularization that is used during Chi2 method
-    static Float_t fRegularizationWeight;          // factor for regularization term
+    static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
+    static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
 
-    TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2 (0..3)
-    TH2F* fMultiplicityVtx[kMCHists];  // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2, inf (0..4)
-    TH2F* fMultiplicityMB[kMCHists];   // multiplicity histogram of triggered events                        : vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2, inf (0..4)
-    TH2F* fMultiplicityINEL[kMCHists]; // multiplicity histogram of all (inelastic) events                  : vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2, inf (0..4)
+    TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2 (0..3)
+
+    TH2F* fMultiplicityVtx[kMCHists];  // multiplicity histogram of events that have a reconstructed vertex : vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2, inf (0..4)
+    TH2F* fMultiplicityMB[kMCHists];   // multiplicity histogram of triggered events                        : vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2, inf (0..4)
+    TH2F* fMultiplicityINEL[kMCHists]; // multiplicity histogram of all (inelastic) events                  : vtx vs multiplicity; array: |eta| < 0.5, 0.9, 1.5, 2, inf (0..4)
 
     TH3F* fCorrelation[kCorrHists];              // vtx vs. (gene multiplicity (trig+vtx)) vs. (meas multiplicity); array: |eta| < 0.5, 1, 1.5, 2 (0..3 and 4..7), the first corrects to the eta range itself, the second to full phase space
+
     TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms
 
-    Float_t fLastChi2MC;        // last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
-    Int_t   fLastChi2MCLimit;   // bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
-    Float_t fLastChi2Residuals; // last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
-    Float_t fRatioAverage;      // last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
+    Float_t fLastChi2MC;        //! last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
+    Int_t   fLastChi2MCLimit;   //! bin where the last chi2 breached a certain threshold, used to evaluate the multiplicity reach (calc. in DrawComparison)
+    Float_t fLastChi2Residuals; //! last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
+    Float_t fRatioAverage;      //! last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
 
-    static Int_t   fgQualityRegionsB[kQualityRegions]; // begin, given in multiplicity units
-    static Int_t   fgQualityRegionsE[kQualityRegions]; // end
-    Float_t fQuality[kQualityRegions];        // stores the quality of the last comparison (calculated in DrawComparison). Contains 3 values that are averages of (MC - unfolded) / e(MC) in 3 regions, these are defined in fQualityRegionB,E
+    static Int_t   fgQualityRegionsB[kQualityRegions]; //! begin, given in multiplicity units
+    static Int_t   fgQualityRegionsE[kQualityRegions]; //! end
+    Float_t fQuality[kQualityRegions];                 //! stores the quality of the last comparison (calculated in DrawComparison). Contains 3 values that are averages of (MC - unfolded) / e(MC) in 3 regions, these are defined in fQualityRegionB,E
 
  private:
     AliMultiplicityCorrection(const AliMultiplicityCorrection&);
     AliMultiplicityCorrection& operator=(const AliMultiplicityCorrection&);
 
-  ClassDef(AliMultiplicityCorrection, 1);
+  ClassDef(AliMultiplicityCorrection, 2);
 };
 
 #endif
index 9a115aa..04226a0 100644 (file)
@@ -30,8 +30,8 @@
 #include "AliCorrection.h"
 #include "AliCorrectionMatrix3D.h"
 
-//#define TPCMEASUREMENT
-#define ITSMEASUREMENT
+#define TPCMEASUREMENT
+//#define ITSMEASUREMENT
 
 ClassImp(AliMultiplicityMCSelector)
 
@@ -83,6 +83,12 @@ void AliMultiplicityMCSelector::ReadUserObjects(TTree* tree)
 
   if (!fEsdTrackCuts)
      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
+
+  if (!fPtSpectrum && fInput)
+    fPtSpectrum = dynamic_cast<TH1*> (fInput->FindObject("pt-spectrum"));
+
+  if (!fPtSpectrum && tree)
+    fPtSpectrum = dynamic_cast<TH1*> (tree->GetUserInfo()->FindObject("pt-spectrum"));
 }
 
 void AliMultiplicityMCSelector::SlaveBegin(TTree* tree)
@@ -141,16 +147,24 @@ void AliMultiplicityMCSelector::SlaveBegin(TTree* tree)
 
   if (option.Contains("pt-spectrum-func"))
   {
-    fPtSpectrum = new TH1F("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
-    fPtSpectrum->SetBinContent(1, 1);  
+    if (fPtSpectrum)
+    {
+      Printf("Using function from input list for systematic p_t study");
+    }
+    else
+    {
+      fPtSpectrum = new TH1F("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
+      fPtSpectrum->SetBinContent(1, 1);
+    }
 
     if (fPtSpectrum)
       AliInfo("WARNING: Systematic study enabled. Pt spectrum will be modified");
   }
 
-
   if (option.Contains("particle-species"))
     fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "vtx:Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec:nolabel:doublePrim:doubleCount");
+
+  // TODO set seed for random generator
 }
 
 void AliMultiplicityMCSelector::Init(TTree* tree)
@@ -172,6 +186,7 @@ void AliMultiplicityMCSelector::Init(TTree* tree)
 
     #ifdef TPCMEASUREMENT
       AliESDtrackCuts::EnableNeededBranches(tree);
+      tree->SetBranchStatus("fTracks.fLabel", 1);
     #endif
 
     tree->SetCacheSize(0);
@@ -285,7 +300,7 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
 
   // tracks in different eta ranges
   Int_t nMCTracks05 = 0;
-  Int_t nMCTracks10 = 0;
+  Int_t nMCTracks09 = 0;
   Int_t nMCTracks15 = 0;
   Int_t nMCTracks20 = 0;
   Int_t nMCTracksAll = 0;
@@ -294,6 +309,12 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
   Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
   for (Int_t i = 0; i<4; ++i)
     nMCTracksSpecies[i] = 0;
+  // eta range for nMCTracksSpecies, nESDTracksSpecies
+#ifdef ITSMEASUREMENT
+  const Float_t etaRange = 2.0;
+#else
+  const Float_t etaRange = 0.9;
+#endif
 
   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
   {
@@ -327,24 +348,24 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
 
     // in case of systematic study, weight according to the change of the pt spectrum
     // it cannot be just multiplied because we cannot count "half of a particle"
-    // instead a random generator decides if the particle is counted twice (if value > 1) 
+    // instead a random generator decides if the particle is counted twice (if value > 1)
     // or not (if value < 0)
-    if (fPtSpectrum) 
+    if (fPtSpectrum)
     {
       Int_t bin = fPtSpectrum->FindBin(pt);
       if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
       {
-       Float_t factor = fPtSpectrum->GetBinContent(bin);
-       if (factor > 0)
-       {
-         Float_t random = gRandom->Uniform();
-         if (factor > 1 && random < factor - 1)
-         {
-           particleWeight = 2;
-         }
-         else if (factor < 1 && random < 1 - factor)
-           particleWeight = 0;
-       }
+        Float_t factor = fPtSpectrum->GetBinContent(bin);
+        if (factor > 0)
+        {
+          Float_t random = gRandom->Uniform();
+          if (factor > 1 && random < factor - 1)
+          {
+            particleWeight = 2;
+          }
+          else if (factor < 1 && random < 1 - factor)
+            particleWeight = 0;
+        }
       }
     }
 
@@ -353,8 +374,8 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
     if (TMath::Abs(particle->Eta()) < 0.5)
       nMCTracks05 += particleWeight;
 
-    if (TMath::Abs(particle->Eta()) < 1.0)
-      nMCTracks10 += particleWeight;
+    if (TMath::Abs(particle->Eta()) < 0.9)
+      nMCTracks09 += particleWeight;
 
     if (TMath::Abs(particle->Eta()) < 1.5)
       nMCTracks15 += particleWeight;
@@ -374,20 +395,20 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
         case 2212: id = 2; break;
         default: id = 3; break;
       }
-      if (TMath::Abs(particle->Eta()) < 2.0)
+      if (TMath::Abs(particle->Eta()) < etaRange)
         nMCTracksSpecies[id]++;
       if (fParticleCorrection[id])
         fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
     }
   }// end of mc particle
 
-  fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
+  fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, (Int_t) nMCTracks05, (Int_t) nMCTracks09, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll);
 
   if (!eventTriggered || !eventVertex)
     return kTRUE;
 
   Int_t nESDTracks05 = 0;
-  Int_t nESDTracks10 = 0;
+  Int_t nESDTracks09 = 0;
   Int_t nESDTracks15 = 0;
   Int_t nESDTracks20 = 0;
 
@@ -400,7 +421,7 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
   for (Int_t i=0; i<nPrim; i++)
     foundPrimaries[i] = kFALSE;
 
-  Bool_t* foundTracks = new Bool_t[nMCPart]; // to prevent double counting
+  Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
   for (Int_t i=0; i<nMCPart; i++)
     foundTracks[i] = kFALSE;
 
@@ -420,6 +441,9 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
 
     Float_t theta = mult->GetTheta(i);
     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
+
+    // TODO define needed, because this only works with new AliRoot
+    Int_t label = mult->GetLabel(i);
 #endif
 #ifdef TPCMEASUREMENT
   // get multiplicity from ESD tracks
@@ -436,13 +460,15 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
     }
 
     Double_t p[3];
-    esdTrack->GetConstrainedPxPyPz(p); // ### TODO should be okay because we have a vertex, however GetInnerPxPyPy / GetOuterPxPyPy also exist
+    esdTrack->GetConstrainedPxPyPz(p); 
     TVector3 vector(p);
 
     Float_t theta = vector.Theta();
     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
-    Float_t pt = vector.Pt();
 
+    Int_t label = esdTrack->GetLabel();
+
+    //Float_t pt = vector.Pt();
     //if (pt < kPtCut)
     //  continue;
 #endif
@@ -454,16 +480,16 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
     {
       TParticle* mother = 0;
 
-      // TODO define needed, because this only works with new AliRoot
-      Int_t label = mult->GetLabel(i);
-      if (label >= 0)
-       label = AliPWG0depHelper::FindPrimaryMotherLabel(stack, label);
-      if (label >= 0)
-       mother = stack->Particle(label);
+      // preserve label for later
+      Int_t labelCopy = label;
+      if (labelCopy >= 0)
+        labelCopy = AliPWG0depHelper::FindPrimaryMotherLabel(stack, labelCopy);
+      if (labelCopy >= 0)
+        mother = stack->Particle(labelCopy);
 
-      // in this case we do not count it!!!
+      // in case of pt study we do not count particles w/o label, because they cannot be scaled
       if (!mother)
-       continue;
+        continue;
 
       // it cannot be just multiplied because we cannot count "half of a particle"
       // instead a random generator decides if the particle is counted twice (if value > 1) 
@@ -471,17 +497,17 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
       Int_t bin = fPtSpectrum->FindBin(mother->Pt());
       if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
       {
-       Float_t factor = fPtSpectrum->GetBinContent(bin);
-       if (factor > 0)
-       {
-         Float_t random = gRandom->Uniform();
-         if (factor > 1 && random < factor - 1)
-         {
-           particleWeight = 2;
-         }
-         else if (factor < 1 && random < 1 - factor)
-           particleWeight = 0;
-       }
+        Float_t factor = fPtSpectrum->GetBinContent(bin);
+        if (factor > 0)
+        {
+          Float_t random = gRandom->Uniform();
+          if (factor > 1 && random < factor - 1)
+          {
+            particleWeight = 2;
+          }
+          else if (factor < 1 && random < 1 - factor)
+            particleWeight = 0;
+        }
       }
     }
 
@@ -490,8 +516,8 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
     if (TMath::Abs(eta) < 0.5)
       nESDTracks05 += particleWeight;
 
-    if (TMath::Abs(eta) < 1.0)
-      nESDTracks10 += particleWeight;
+    if (TMath::Abs(eta) < 0.9)
+      nESDTracks09 += particleWeight;
 
     if (TMath::Abs(eta) < 1.5)
       nESDTracks15 += particleWeight;
@@ -499,49 +525,27 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
     if (TMath::Abs(eta) < 2.0)
       nESDTracks20 += particleWeight;
 
+
     if (fParticleCorrection[0] || fParticleSpecies)
     {
-      // TODO define needed, because this only works with new AliRoot
-      Int_t label = mult->GetLabel(i);
-
-      // TODO double counting ok for particle ratio study!!!
-      // check if we already counted this particle, this way distinguishes double counted particles (bug) or double counted primaries due to secondaries (physics)
-      if (label >= 0)
-      {
-       if (foundTracks[label])
-       {
-         if (TMath::Abs(eta) < 2.0)
-           nESDTracksSpecies[6]++;
-         continue;
-       }
-       foundTracks[label] = kTRUE;
-      }
-
+      Int_t motherLabel = -1;
       TParticle* mother = 0;
 
       // find mother
       if (label >= 0)
-       label = AliPWG0depHelper::FindPrimaryMotherLabel(stack, label);
-      if (label >= 0)
-       mother = stack->Particle(label);
+        motherLabel = AliPWG0depHelper::FindPrimaryMotherLabel(stack, label);
+      if (motherLabel >= 0)
+        mother = stack->Particle(motherLabel);
 
       if (!mother)
       {
-       // count tracks that did not have a label
-       if (TMath::Abs(eta) < 2.0)
-         nESDTracksSpecies[4]++;
+        // count tracks that did not have a label
+        if (TMath::Abs(eta) < etaRange)
+          nESDTracksSpecies[4]++;
         continue;
-      } 
-
-      // particle (primary) already counted?
-      if (foundPrimaries[label])
-      {
-       if (TMath::Abs(eta) < 2.0)
-         nESDTracksSpecies[5]++;
-       continue;
       }
-      foundPrimaries[label] = kTRUE;
 
+      // get particle type (pion, proton, kaon, other)
       Int_t id = -1;
       switch (TMath::Abs(mother->GetPdgCode()))
       {
@@ -550,8 +554,31 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
         case 2212: id = 2; break;
         default: id = 3; break;
       }
-      if (TMath::Abs(eta) < 2.0)
+
+      // double counting is ok for particle ratio study
+      if (TMath::Abs(eta) < etaRange)
         nESDTracksSpecies[id]++;
+
+      // double counting is not ok for efficiency study
+
+      // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
+      if (foundTracks[label])
+      {
+        if (TMath::Abs(eta) < etaRange)
+          nESDTracksSpecies[6]++;
+        continue;
+      }
+      foundTracks[label] = kTRUE;
+
+      // particle (primary) already counted?
+      if (foundPrimaries[motherLabel])
+      {
+        if (TMath::Abs(eta) < etaRange)
+          nESDTracksSpecies[5]++;
+        continue;
+      }
+      foundPrimaries[motherLabel] = kTRUE;
+
       if (fParticleCorrection[id])
         fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
     }
@@ -561,18 +588,18 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
   delete[] foundPrimaries;
 
   if ((Int_t) nMCTracks20 == 0 && nESDTracks20 > 3)
-    printf("WARNING: Event %lld has %d generated and %d reconstructed...\n", entry, (Int_t) nMCTracks20, (Int_t) nESDTracks20);
+    printf("WARNING: Event %lld has %d generated and %d reconstructed...\n", entry, nMCTracks20, nESDTracks20);
 
   //printf("%.2f generated and %.2f reconstructed\n", nMCTracks20, nESDTracks20);
 
-  //Printf("%f %f %f %f %f %f %f %f %f %f", vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
+  //Printf("%f %f %f %f %f %f %f %f %f %f", vtxMC[2], nMCTracks05, nMCTracks09, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
 
-  //Printf("%f %f %f %f %f %f %f %f %f %f", vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
+  //Printf("%f %f %f %f %f %f %f %f %f %f", vtxMC[2], nMCTracks05, nMCTracks09, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
 
-  fMultiplicity->FillMeasured(vtx[2], (Int_t) nESDTracks05, (Int_t) nESDTracks10, (Int_t) nESDTracks15, (Int_t) nESDTracks20);
+  fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks09, nESDTracks15, nESDTracks20);
 
   // fill response matrix using vtxMC (best guess)
-  fMultiplicity->FillCorrection(vtxMC[2], (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks15, (Int_t) nMCTracks20, (Int_t) nMCTracksAll, (Int_t) nESDTracks05, (Int_t) nESDTracks10, (Int_t) nESDTracks15, (Int_t) nESDTracks20);
+  fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks09,  nMCTracks15,  nMCTracks20,  nMCTracksAll,  nESDTracks05,  nESDTracks09,  nESDTracks15,  nESDTracks20);
 
   if (fParticleSpecies)
     fParticleSpecies->Fill(vtxMC[2], nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
index 94b6828..79c552f 100644 (file)
@@ -9,7 +9,6 @@ class AliESDtrackCuts;
 class AliMultiplicityCorrection;
 class AliCorrection;
 class TNtuple;
-class TF1;
 class TH1;
 
 class AliMultiplicityMCSelector : public AliSelectorRL {
@@ -34,7 +33,9 @@ class AliMultiplicityMCSelector : public AliSelectorRL {
     AliCorrection* fParticleCorrection[4]; // correction from measured to generated particles for trigger, vertex sample in |eta| < 2;
                                            // for each of the species: pi, k, p, other; for systematic study of pt cut off
     Int_t fSelectProcessType;        // 0 = all (default), 1 = ND, 2 = SD, 3 = DD (for systematic study)
-    TNtuple *fParticleSpecies;       // per event: vtx_mc, (pi, k, p, rest (in |eta| < 2)) X (true, recon); (for systematic study)
+    TNtuple *fParticleSpecies;       // per event: vtx_mc, (pi, k, p, rest (in |eta| < 2)) X (true, recon) + (nolabel,
+                                     // doubleTracks, doublePrimaries) [doubleTracks + doublePrimaries are already part of
+                                     // rec. particles!)
 
     TH1* fPtSpectrum;                // function that modifies the pt spectrum (syst. study)
 
index 6d418e7..d2c7ac3 100644 (file)
@@ -2,9 +2,33 @@
 // plots for the note about multplicity measurements
 //
 
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TLine.h>
+#include <TF1.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TLegend.h>
+#include <TStopwatch.h>
+#include <TROOT.h>
+#include <TGraph.h>
+#include <TMath.h>
+
+#include "AliMultiplicityCorrection.h"
+#include "AliCorrection.h"
+#include "AliCorrectionMatrix3D.h"
+
+#endif
+
 const char* correctionFile = "multiplicityMC_2M.root";
 const char* measuredFile   = "multiplicityMC_1M_3.root";
 Int_t etaRange = 3;
+Int_t drawRatioRange = 150;
 
 const char* correctionFileTPC = "multiplicityMC_TPC_1.4M.root";
 const char* measuredFileTPC   = "multiplicityMC_TPC_0.6M.root";
@@ -15,6 +39,23 @@ void SetTPC()
   correctionFile = correctionFileTPC;
   measuredFile = measuredFileTPC;
   etaRange = etaRangeTPC;
+  drawRatioRange = 100;
+}
+
+void Smooth(TH1* hist, Int_t windowWidth = 20)
+{
+  TH1* clone = (TH1*) hist->Clone("clone");
+  for (Int_t bin=2; bin<=clone->GetNbinsX(); ++bin)
+  {
+    Int_t min = TMath::Max(2, bin-windowWidth);
+    Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
+    Float_t average = clone->Integral(min, max) / (max - min + 1);
+
+    hist->SetBinContent(bin, average);
+    hist->SetBinError(bin, 0);
+  }
+
+  delete clone;
 }
 
 void responseMatrixPlot()
@@ -52,10 +93,10 @@ TCanvas* DrawResultRatio(TH1* mcHist, TH1* result, TString epsName)
   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
   canvas->Range(0, 0, 1, 1);
 
-  TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
+  TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
   pad1->Draw();
 
-  TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
+  TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
   pad2->Draw();
 
   pad1->SetRightMargin(0.05);
@@ -147,10 +188,10 @@ TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsNa
   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 600);
   canvas->Range(0, 0, 1, 1);
 
-  TPad* pad1 = new TPad(Form("%s_pad1", epsName), "", 0, 0.5, 0.98, 0.98);
+  TPad* pad1 = new TPad(Form("%s_pad1", epsName.Data()), "", 0, 0.5, 0.98, 0.98);
   pad1->Draw();
 
-  TPad* pad2 = new TPad(Form("%s_pad2", epsName), "", 0, 0.02, 0.98, 0.5);
+  TPad* pad2 = new TPad(Form("%s_pad2", epsName.Data()), "", 0, 0.02, 0.98, 0.5);
   pad2->Draw();
 
   pad1->SetRightMargin(0.05);
@@ -168,7 +209,7 @@ TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsNa
   mcHist->GetYaxis()->SetTitleSize(0.06);
   mcHist->GetYaxis()->SetTitleOffset(0.6);
 
-  mcHist->GetXaxis()->SetRangeUser(0, 150);
+  mcHist->GetXaxis()->SetRangeUser(0, drawRatioRange);
 
   mcHist->SetTitle(";true multiplicity;Entries");
   mcHist->SetStats(kFALSE);
@@ -198,7 +239,7 @@ TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsNa
   result1->GetYaxis()->SetTitleSize(0.06);
   result1->GetYaxis()->SetTitleOffset(0.6);
 
-  result1->GetXaxis()->SetRangeUser(0, 150);
+  result1->GetXaxis()->SetRangeUser(0, drawRatioRange);
 
   result1->SetTitle(";true multiplicity;Entries");
   result1->SetStats(kFALSE);
@@ -215,23 +256,23 @@ TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsNa
 
   // get average of ratio
   Float_t sum = 0;
-  for (Int_t i=2; i<=150; ++i)
+  for (Int_t i=2; i<=drawRatioRange; ++i)
   {
     sum += TMath::Abs(ratio->GetBinContent(i) - 1);
   }
-  sum /= 149;
+  sum /= drawRatioRange-1;
 
-  printf("Average (2..150) of |ratio - 1| is %f\n", sum);
+  printf("Average (2..%d) of |ratio - 1| is %f\n", drawRatioRange, sum);
 
-  TLine* line = new TLine(0, 1, 150, 1);
+  TLine* line = new TLine(0, 1, drawRatioRange, 1);
   line->SetLineWidth(2);
   line->Draw();
 
-  line = new TLine(0, 1.1, 150, 1.1);
+  line = new TLine(0, 1.1, drawRatioRange, 1.1);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
-  line = new TLine(0, 0.9, 150, 0.9);
+  line = new TLine(0, 0.9, drawRatioRange, 0.9);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
@@ -243,7 +284,7 @@ TCanvas* Draw2ResultRatio(TH1* mcHist, TH1* result1, TH1* result2, TString epsNa
   return canvas;
 }
 
-TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE)
+TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString epsName, Bool_t firstMarker = kFALSE, const char** legendStrings = 0)
 {
   // compares n results with first results. E.g. one gained with the default response, another with a changed one to study
   // a systematic effect
@@ -253,9 +294,14 @@ TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString eps
 
   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
 
-  result->GetXaxis()->SetRangeUser(0, 150);
+  result->GetXaxis()->SetRangeUser(1, drawRatioRange);
   result->SetStats(kFALSE);
 
+  Int_t colors[] = {1, 2, 4, 6, 7, 8, 9, 10};
+
+  TLegend* legend = new TLegend(0.2, 0.75, 0.35, 0.95);
+  legend->SetFillColor(0);
+
   for (Int_t n=0; n<nResultSyst; ++n)
   {
     resultSyst[n]->Scale(1.0 / resultSyst[n]->Integral(2, 200));
@@ -269,28 +315,36 @@ TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString eps
     if (firstMarker)
       ratio->SetMarkerStyle(5);
 
-    ratio->SetLineColor(n+1);
+    ratio->SetLineColor(colors[n / 2]);
+    if ((n % 2))
+      ratio->SetLineStyle(2);
 
     ratio->DrawCopy((n == 0) ? ((firstMarker) ? "P" : "HIST") : "SAME HIST");
 
+    if (legendStrings && legendStrings[n])
+      legend->AddEntry(ratio, legendStrings[n]);
+
     // get average of ratio
     Float_t sum = 0;
-    for (Int_t i=2; i<=150; ++i)
+    for (Int_t i=2; i<=drawRatioRange; ++i)
       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
-    sum /= 149;
+    sum /= drawRatioRange-1;
 
-    printf("%d) Average (2..150) of |ratio - 1| is %f\n", n, sum);
+    printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, drawRatioRange, sum);
   }
 
-  TLine* line = new TLine(0, 1, 150, 1);
+  if (legendStrings)
+    legend->Draw();
+
+  TLine* line = new TLine(1, 1, drawRatioRange, 1);
   line->SetLineWidth(2);
   line->Draw();
 
-  line = new TLine(0, 1.1, 150, 1.1);
+  line = new TLine(1, 1.1, drawRatioRange, 1.1);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
-  line = new TLine(0, 0.9, 150, 0.9);
+  line = new TLine(1, 0.9, drawRatioRange, 0.9);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
@@ -303,11 +357,13 @@ TCanvas* DrawRatio(TH1* result, Int_t nResultSyst, TH1** resultSyst, TString eps
   return canvas;
 }
 
-TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
+TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName, Bool_t smooth = kFALSE, Bool_t dashed = kFALSE)
 {
   // draws the ratios of each mc to the corresponding result
 
   TCanvas* canvas = new TCanvas(epsName, epsName, 800, 400);
+  canvas->SetRightMargin(0.05);
+  canvas->SetTopMargin(0.05);
 
   for (Int_t n=0; n<nResultSyst; ++n)
   {
@@ -315,37 +371,50 @@ TCanvas* DrawRatio(Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 
-    result[n]->GetXaxis()->SetRangeUser(0, 150);
+    result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
     result[n]->SetStats(kFALSE);
 
     // calculate ratio
     TH1* ratio = (TH1*) result[n]->Clone("ratio");
     ratio->Divide(mc[n], ratio, 1, 1, "B");
-    ratio->SetTitle(";true multiplicity;Ratio (true / unfolded)");
+
+    // SetRangeUser(1, ...) would be the same, but the 0 should be still on the axis...
+    ratio->SetBinContent(1, 1); ratio->SetBinError(1, 0);
+
+    if (smooth)
+      Smooth(ratio);
+
+    ratio->SetTitle(Form(";true multiplicity;Ratio (true / unfolded)%s", ((smooth) ? " (smoothed)" : "")));
     ratio->GetYaxis()->SetRangeUser(0.55, 1.45);
 
-    ratio->SetLineColor(n+1);
+    if (dashed)
+    {
+      ratio->SetLineColor((n/2)+1);
+      ratio->SetLineStyle((n%2)+1);
+    }
+    else
+      ratio->SetLineColor(n+1);
 
     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
 
     // get average of ratio
     Float_t sum = 0;
-    for (Int_t i=2; i<=150; ++i)
+    for (Int_t i=2; i<=drawRatioRange; ++i)
       sum += TMath::Abs(ratio->GetBinContent(i) - 1);
-    sum /= 149;
+    sum /= drawRatioRange-1;
 
-    printf("%d) Average (2..150) of |ratio - 1| is %f\n", n, sum);
+    printf("%d) Average (2..%d) of |ratio - 1| is %f\n", n, drawRatioRange, sum);
   }
 
-  TLine* line = new TLine(0, 1, 150, 1);
+  TLine* line = new TLine(0, 1, drawRatioRange, 1);
   line->SetLineWidth(2);
   line->Draw();
 
-  line = new TLine(0, 1.1, 150, 1.1);
+  line = new TLine(0, 1.1, drawRatioRange, 1.1);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
-  line = new TLine(0, 0.9, 150, 0.9);
+  line = new TLine(0, 0.9, drawRatioRange, 0.9);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
@@ -381,7 +450,7 @@ TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1**
     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 
-    result[n]->GetXaxis()->SetRangeUser(0, 150);
+    result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
     result[n]->SetStats(kFALSE);
 
     // calculate ratio
@@ -396,22 +465,22 @@ TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1**
 
     // get average of ratio
     Float_t sum = 0;
-    for (Int_t i=2; i<=150; ++i)
+    for (Int_t i=2; i<=drawRatioRange; ++i)
       sum += TMath::Abs(ratio->GetBinContent(i));
-    sum /= 149;
+    sum /= drawRatioRange-1;
 
-    printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
+    printf("%d) Average (2..%d) of |ratio - ratioBase| is %f\n", n, drawRatioRange, sum);
   }
 
-  TLine* line = new TLine(0, 0, 150, 0);
+  TLine* line = new TLine(0, 0, drawRatioRange, 0);
   line->SetLineWidth(2);
   line->Draw();
 
-  line = new TLine(0, 0.1, 150, 0.1);
+  line = new TLine(0, 0.1, drawRatioRange, 0.1);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
-  line = new TLine(0, -0.1, 150, -0.1);
+  line = new TLine(0, -0.1, drawRatioRange, -0.1);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
@@ -424,22 +493,6 @@ TCanvas* DrawRatioDeduct(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1**
   return canvas;
 }
 
-void Smooth(TH1* hist, Int_t windowWidth = 20)
-{
-  TH1* clone = (TH1*) hist->Clone("clone");
-  for (Int_t bin=3; bin<=clone->GetNbinsX(); ++bin)
-  {
-    Int_t min = TMath::Max(3, bin-windowWidth);
-    Int_t max = TMath::Min(clone->GetNbinsX(), bin+windowWidth);
-    Float_t average = clone->Integral(min, max) / (max - min + 1);
-
-    hist->SetBinContent(bin, average);
-    hist->SetBinError(bin, 0);
-  }
-
-  delete clone;
-}
-
 TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst, TH1** mc, TH1** result, TString epsName)
 {
   // draws the ratios of each mc to the corresponding result
@@ -464,7 +517,7 @@ TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst,
     result[n]->Scale(1.0 / result[n]->Integral(2, 200));
     mc[n]->Scale(1.0 / mc[n]->Integral(2, 200));
 
-    result[n]->GetXaxis()->SetRangeUser(0, 150);
+    result[n]->GetXaxis()->SetRangeUser(0, drawRatioRange);
     result[n]->SetStats(kFALSE);
 
     // calculate ratio
@@ -472,15 +525,19 @@ TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst,
     ratio->Divide(mc[n], ratio, 1, 1, "B");
     ratio->Add(ratioBase, -1);
 
-    new TCanvas; ratio->DrawCopy();
+    //new TCanvas; ratio->DrawCopy();
+    // clear 0 bin
+    ratio->SetBinContent(1, 0); ratio->SetBinError(1, 0);
+
     Smooth(ratio);
-    ratio->SetLineColor(1);
-    ratio->DrawCopy("SAME");
+
+    //ratio->SetLineColor(1); ratio->DrawCopy("SAME");
 
     canvas->cd();
-    ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u)");
-    ratio->GetYaxis()->SetRangeUser(-1, 1);
-    ratio->SetLineColor(n+1);
+    ratio->SetTitle(";true multiplicity;Ratio_{syst} (t/u) - Ratio (t/u) (smoothed)");
+    ratio->GetYaxis()->SetRangeUser(-0.3, 0.3);
+    ratio->SetLineColor((n / 2)+1);
+    ratio->SetLineStyle((n % 2)+1);
     ratio->DrawCopy((n == 0) ? "HIST" : "SAME HIST");
 
     // get average of ratio
@@ -492,15 +549,15 @@ TCanvas* DrawRatioDeductSmooth(TH1* mcBase, TH1* resultBase, Int_t nResultSyst,
     printf("%d) Average (2..150) of |ratio - ratioBase| is %f\n", n, sum);
   }
 
-  TLine* line = new TLine(0, 0, 150, 0);
+  TLine* line = new TLine(0, 0, drawRatioRange, 0);
   line->SetLineWidth(2);
   line->Draw();
 
-  line = new TLine(0, 0.1, 150, 0.1);
+  line = new TLine(0, 0.1, drawRatioRange, 0.1);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
-  line = new TLine(0, -0.1, 150, -0.1);
+  line = new TLine(0, -0.1, drawRatioRange, -0.1);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
@@ -683,7 +740,7 @@ void chi2FluctuationResult()
   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 
   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+  //TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
 
   mult->DrawComparison("MinuitChi2", etaRange, kFALSE, kTRUE, mcHist, kTRUE);
 
@@ -815,9 +872,9 @@ void minimizationInfluenceAlpha()
 
   TFile::Open("eval-2M-1M/EvaluateChi2MethodDetail.root");
 
-  TH1* hist1 = gFile->Get("MinuitChi2_00_2_100.000000");
-  TH1* hist2 = gFile->Get("MinuitChi2_03_2_100000.000000");
-  TH1* hist3 = gFile->Get("MinuitChi2_06_2_100000000.000000");
+  TH1* hist1 = (TH1*) gFile->Get("MinuitChi2_00_2_100.000000");
+  TH1* hist2 = (TH1*) gFile->Get("MinuitChi2_03_2_100000.000000");
+  TH1* hist3 = (TH1*) gFile->Get("MinuitChi2_06_2_100000000.000000");
 
   mcHist->Rebin(2);  mcHist->Scale(0.5);
   hist1->Rebin(2);   hist1->Scale(0.5);
@@ -955,7 +1012,7 @@ void StartingConditions()
   mc->Sumw2();
   mc->Scale(1.0 / mc->Integral());
 
-  Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
+  //Int_t marker[] = {24, 25, 26, 27, 28, 2, 3, 4, 5};
 
   for (Int_t i=0; i<6; ++i)
   {
@@ -1008,7 +1065,7 @@ void StatisticsPlot()
 
   TCanvas* canvas = new TCanvas(name, name, 600, 400);
 
-  TGraph* fitResultsChi2 = file->Get("fitResultsChi2");
+  TGraph* fitResultsChi2 = (TGraph*) file->Get("fitResultsChi2");
   fitResultsChi2->SetTitle(";number of measured events;P_{1}");
   fitResultsChi2->GetYaxis()->SetRangeUser(0, 2);
   fitResultsChi2->Draw("AP");
@@ -1023,9 +1080,8 @@ void StatisticsPlot()
 
   const char* plotname = "chi2Result";
 
-  const char* name = "StatisticsPlotRatios";
-
-  TCanvas* canvas = new TCanvas(name, name, 600, 400);
+  name = "StatisticsPlotRatios";
+  canvas = new TCanvas(name, name, 600, 400);
 
   for (Int_t i=0; i<5; ++i)
   {
@@ -1055,7 +1111,7 @@ void SystematicLowEfficiency()
   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 
   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
-  TH1* result1 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
+  TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
 
   DrawResultRatio(mcHist, result1, "SystematicLowEfficiencyLow.eps");
 
@@ -1066,7 +1122,7 @@ void SystematicLowEfficiency()
   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
   mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 
-  TH1* result2 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
+  TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
 
   DrawResultRatio(mcHist, result2, "SystematicLowEfficiencyOK.eps");
 
@@ -1095,112 +1151,151 @@ void SystematicMisalignment()
   DrawResultRatio(mcHist, result, "SystematicMisalignment.eps");
 }
 
-void EfficiencySpecies(const char* fileName = "multiplicityMC_400k_syst.root")
+void SystematicMisalignmentTPC()
 {
   gSystem->Load("libPWG0base");
 
-  TFile::Open(fileName);
-  AliCorrection* correction[4];
+  SetTPC();
 
-  TCanvas* canvas = new TCanvas("EfficiencySpeciesSPD", "EfficiencySpeciesSPD", 600, 400);
-  canvas->SetGridx();
-  canvas->SetGridy();
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
+  TFile::Open(correctionFile);
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  mult->LoadHistograms("Multiplicity");
 
-  TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
-  legend->SetFillColor(0);
+  TFile::Open("multiplicityMC_TPC_100k_fullmis.root");
+  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+  mult2->LoadHistograms("Multiplicity");
+
+  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+  TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
+
+  DrawResultRatio(mcHist, result, "SystematicMisalignmentTPC.eps");
+}
+
+void EfficiencySpecies()
+{
+  gSystem->Load("libPWG0base");
 
   Int_t marker[] = {24, 25, 26};
   Int_t color[] = {1, 2, 4};
 
-  Float_t below = 0;
-  Float_t total = 0;
+  // SPD TPC
+  const char* fileName[] = { "multiplicityMC_400k_syst.root", "multiplicityMC_TPC_4kfiles_syst.root" };
+  Float_t etaRange[] = {0.49, 0.9};
+  const char* titles[] = { "SPD Tracklets", "TPC Tracks" };
+
+  TCanvas* canvas = new TCanvas("EfficiencySpecies", "EfficiencySpecies", 1000, 500);
+  canvas->Divide(2, 1);
 
-  for (Int_t i=0; i<3; ++i)
+  for (Int_t loop=0; loop<2; ++loop)
   {
-    Printf("correction %d", i);
+    Printf("%s", fileName[loop]);
 
-    TString name; name.Form("correction_%d", i);
-    correction[i] = new AliCorrection(name, name);
-    correction[i]->LoadHistograms();
+    TFile::Open(fileName[loop]);
+    AliCorrection* correction[4];
 
-    TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
-    TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
+    canvas->cd(loop+1);
 
-    // limit vtx axis
-    gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
-    meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
+    gPad->SetGridx();
+    gPad->SetGridy();
+    gPad->SetRightMargin(0.05);
+    //gPad->SetTopMargin(0.05);
 
-    // empty over/underflow bin in eta, setting range to +-2 is not enough because this is the maximum range, Project3D takes them into account then (might be a bug)
-    /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
-      for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
-      {
-        gene->SetBinContent(x, 0, z, 0);
-        gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
-        meas->SetBinContent(x, 0, z, 0);
-        meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
-      }*/
+    TLegend* legend = new TLegend(0.7, 0.4, 0.85, 0.6);
+    legend->SetFillColor(0);
+    legend->SetEntrySeparation(0.2);
+
+    Float_t below = 0;
+    Float_t total = 0;
 
-    // limit eta axis
-    gene->GetYaxis()->SetRangeUser(-0.49, 0.49);
-    meas->GetYaxis()->SetRangeUser(-0.49, 0.49);
+    for (Int_t i=0; i<3; ++i)
+    {
+      Printf("correction %d", i);
 
-    TH1* genePt = gene->Project3D(Form("z_%d", i));
-    TH1* measPt = meas->Project3D(Form("z_%d", i));
+      TString name; name.Form("correction_%d", i);
+      correction[i] = new AliCorrection(name, name);
+      correction[i]->LoadHistograms();
 
-    genePt->Sumw2();
-    measPt->Sumw2();
+      TH3* gene = correction[i]->GetTrackCorrection()->GetGeneratedHistogram();
+      TH3* meas = correction[i]->GetTrackCorrection()->GetMeasuredHistogram();
 
-    effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
-    effPt->Reset();
-    effPt->Divide(measPt, genePt, 1, 1, "B");
+      // limit vtx axis
+      gene->GetXaxis()->SetRangeUser(-3.9, 3.9);
+      meas->GetXaxis()->SetRangeUser(-3.9, 3.9);
 
-    Int_t bin = 0;
-    for (bin=20; bin>=1; bin--)
-    {
-      if (effPt->GetBinContent(bin) < 0.5)
-        break;
-    }
+      // empty over/underflow bin in eta, setting range to +-2 is not enough because this is the maximum range, Project3D takes them into account then (might be a bug)
+      /*for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
+        for (Int_t z = 1; z <= gene->GetNbinsZ(); z++)
+        {
+          gene->SetBinContent(x, 0, z, 0);
+          gene->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
+          meas->SetBinContent(x, 0, z, 0);
+          meas->SetBinContent(x, gene->GetNbinsY()+1, z, 0);
+        }*/
 
-    Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
+      // limit eta axis
+      gene->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
+      meas->GetYaxis()->SetRangeUser(-etaRange[loop], etaRange[loop]);
 
-    Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
-    Printf("%.4f of the particles are below that momentum", fraction);
+      TH1* genePt = gene->Project3D(Form("z_%d", i));
+      TH1* measPt = meas->Project3D(Form("z_%d", i));
 
-    below += genePt->Integral(1, bin);
-    total += genePt->Integral();
+      genePt->Sumw2();
+      measPt->Sumw2();
 
-    effPt->SetLineColor(color[i]);
-    effPt->SetMarkerColor(color[i]);
-    effPt->SetMarkerStyle(marker[i]);
+      TH1* effPt = (TH1*) genePt->Clone(Form("effPt_%d", i));
+      effPt->Reset();
+      effPt->Divide(measPt, genePt, 1, 1, "B");
 
-    effPt->GetXaxis()->SetRangeUser(0, 1);
-    effPt->GetYaxis()->SetRangeUser(0, 1);
+      Int_t bin = 0;
+      for (bin=20; bin>=1; bin--)
+      {
+        if (effPt->GetBinContent(bin) < 0.5)
+          break;
+      }
 
-    effPt->SetStats(kFALSE);
-    effPt->SetTitle("");
-    effPt->GetYaxis()->SetTitle("Efficiency");
+      Printf("Eff. below 50%% at bin %d, i.e. %.3f GeV/c", bin, effPt->GetXaxis()->GetBinUpEdge(bin));
 
-    effPt->DrawCopy((i == 0) ? "" : "SAME");
+      Float_t fraction = genePt->Integral(1, bin) / genePt->Integral();
+      Printf("%.4f of the particles are below that momentum", fraction);
 
-    legend->AddEntry(effPt, ((i == 0) ? "#pi" : ((i == 1) ? "K" : "p")));
-  }
+      below += genePt->Integral(1, bin);
+      total += genePt->Integral();
 
-  Printf("In total %.4f of the particles are below their pt cut off", (Float_t) below / total);
+      effPt->SetLineColor(color[i]);
+      effPt->SetMarkerColor(color[i]);
+      effPt->SetMarkerStyle(marker[i]);
 
-  legend->Draw();
+      effPt->GetXaxis()->SetRangeUser(0.06, 1);
+      effPt->GetYaxis()->SetRangeUser(0, 1);
+
+      effPt->GetYaxis()->SetTitleOffset(1.2);
+
+      effPt->SetStats(kFALSE);
+      effPt->SetTitle(titles[loop]);
+      effPt->GetYaxis()->SetTitle("Efficiency");
+
+      effPt->DrawCopy((i == 0) ? "" : "SAME");
+
+      legend->AddEntry(effPt, ((i == 0) ? "#pi^{#pm}" : ((i == 1) ? "K^{#pm}" : "p,#bar{p}")));
+    }
+
+    Printf("In total %.4f of the particles are below their effective pt cut off", (Float_t) below / total);
+
+    legend->Draw();
+  }
 
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 }
 
-void ParticleSpeciesComparison1()
+void ParticleSpeciesComparison1(const char* fileNameMC = "multiplicityMC_400k_syst_species.root", const char* fileNameESD = "multiplicityMC_100k_syst.root")
 {
   gSystem->Load("libPWG0base");
 
-  const char* fileNameMC = "multiplicityMC_400k_syst_species.root";
-  const char* fileNameESD = "multiplicityMC_100k_syst.root";
-  Bool_t chi2 = 0;
+  Bool_t chi2 = 1;
 
   TFile::Open(fileNameESD);
   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", etaRange));
@@ -1252,9 +1347,12 @@ void ParticleSpeciesComparison1()
     mc->SetBinError(i, 0);
   }
 
-  TCanvas* canvas = DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps");
+  const char* legendStrings[] = { "#pi^{#pm}", 0, "K^{#pm}", 0, "p,#bar{p}", 0 };
+
+  DrawRatio(results[0], nMax-1, results+1, "ParticleSpeciesComparison1_2.eps", kFALSE, legendStrings);
 
-  TFile::Open("bayesianUncertainty_400k_100k_syst.root");
+  //not valid: draw chi2 uncertainty on top!
+  /*TFile::Open("bayesianUncertainty_400k_100k_syst.root");
   TH1* errorHist = (TH1*) gFile->Get("errorBoth");
   errorHist->SetLineColor(1);
   errorHist->SetLineWidth(2);
@@ -1266,16 +1364,16 @@ void ParticleSpeciesComparison1()
   }
 
   errorHist->DrawCopy("SAME");
-  errorHist2->DrawCopy("SAME");
+  errorHist2->DrawCopy("SAME");*/
 
-  canvas->SaveAs(canvas->GetName());
+  //canvas->SaveAs(canvas->GetName());
 
-  TCanvas* canvas2 = DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE);
+  DrawRatio(mc, nMax, results, "ParticleSpeciesComparison1_3.eps", kTRUE, 0);
 
-  errorHist->DrawCopy("SAME");
-  errorHist2->DrawCopy("SAME");
+  //errorHist->DrawCopy("SAME");
+  //errorHist2->DrawCopy("SAME");
 
-  canvas2->SaveAs(canvas2->GetName());
+  //canvas2->SaveAs(canvas2->GetName());
 }
 
 void ParticleSpeciesComparison2()
@@ -1395,7 +1493,7 @@ void TriggerVertexCorrection()
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 }
 
-void bayesianUncertainty(Bool_t mc = kFALSE)
+void StatisticalUncertainty(Int_t methodType, Bool_t mc = kFALSE)
 {
   gSystem->Load("libPWG0base");
 
@@ -1411,17 +1509,18 @@ void bayesianUncertainty(Bool_t mc = kFALSE)
 
   TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
 
-  TH1* errorResponse = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorResponse");
-  TH1* errorMeasured = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorMeasured");
-  TH1* errorBoth = mult->BayesianStatisticsEffect(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, 1, 100, ((mc) ? mcHist : 0))->Clone("errorBoth");
+  TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
+
+  TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
+  TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliMultiplicityCorrection::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
 
   if (!mc)
   {
     TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
-    DrawResultRatio(mcHist, result, "bayesianUncertainty2.eps");
+    DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
   }
 
-  TCanvas* canvas = new TCanvas("bayesianUncertainty", "bayesianUncertainty", 600, 400);
+  TCanvas* canvas = new TCanvas("StatisticalUncertainty", "StatisticalUncertainty", 600, 400);
   canvas->SetGridx();
   canvas->SetGridy();
   canvas->SetRightMargin(0.05);
@@ -1438,7 +1537,7 @@ void bayesianUncertainty(Bool_t mc = kFALSE)
   errorMeasured->SetLineColor(2);
   errorMeasured->Draw("SAME");
 
-  errorBoth->SetLineColor(3);
+  errorBoth->SetLineColor(4);
   errorBoth->Draw("SAME");
 
   Printf("Average errorResponse: %f", errorResponse->Integral(2, 150) / 149);
@@ -1454,6 +1553,54 @@ void bayesianUncertainty(Bool_t mc = kFALSE)
   file->Close();
 }
 
+void StatisticalUncertaintyCompare(const char* det = "SPD")
+{
+  TFile* file1 = TFile::Open(Form("StatisticalUncertainty%sBayesian.root", det));
+  TH1* errorResponse = (TH1*) file1->Get("errorResponse");
+  TH1* errorMeasured = (TH1*) file1->Get("errorMeasured");
+  TH1* errorBoth = (TH1*) file1->Get("errorBoth");
+
+  TString str;
+  str.Form("StatisticalUncertaintyCompare%s", det);
+
+  TCanvas* canvas = new TCanvas(str, str, 600, 400);
+  canvas->SetGridx();
+  canvas->SetGridy();
+  canvas->SetRightMargin(0.05);
+  canvas->SetTopMargin(0.05);
+
+  errorResponse->SetLineColor(1);
+  errorResponse->GetXaxis()->SetRangeUser(1, (strcmp(det, "TPC") ? 200 : 100));
+  errorResponse->GetYaxis()->SetRangeUser(0, 0.3);
+  errorResponse->SetStats(kFALSE);
+  errorResponse->SetTitle(";true multiplicity;Uncertainty");
+
+  errorResponse->Draw();
+
+  errorMeasured->SetLineColor(2);
+  errorMeasured->Draw("SAME");
+
+  errorBoth->SetLineColor(4);
+  errorBoth->Draw("SAME");
+
+  TFile* file2 = TFile::Open(Form("StatisticalUncertainty%sChi2.root", det));
+  TH1* errorBoth2 = (TH1*) file2->Get("errorBoth");
+
+  errorBoth2->SetLineColor(4);
+  errorBoth2->SetLineStyle(2);
+  errorBoth2->Draw("SAME");
+
+  TLegend* legend = new TLegend(0.2, 0.6, 0.6, 0.9);
+  legend->SetFillColor(0);
+  legend->AddEntry(errorResponse, "response matrix (Bayesian)");
+  legend->AddEntry(errorMeasured, "measured (Bayesian)");
+  legend->AddEntry(errorBoth, "both (Bayesian)");
+  legend->AddEntry(errorBoth2, "both (#chi^{2} minimization)");
+  legend->Draw();
+
+  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+}
+
 void EfficiencyComparison(Int_t eventType = 2)
 {
  const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root", "multiplicityMC_400k_syst_xsection.root" };
@@ -1492,7 +1639,7 @@ void EfficiencyComparison(Int_t eventType = 2)
     else
       data[i]->LoadHistograms("Multiplicity_0");
 
-    TH1* eff = (TH1*) data[i]->GetEfficiency(3, eventType)->Clone(Form("eff_%d", i));
+    TH1* eff = (TH1*) data[i]->GetEfficiency(3, (AliMultiplicityCorrection::EventType) eventType)->Clone(Form("eff_%d", i));
     effArray[i] = eff;
 
     eff->GetXaxis()->SetRangeUser(0, 15);
@@ -1514,14 +1661,14 @@ void EfficiencyComparison(Int_t eventType = 2)
         AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multtmp", "Multtmp");
         mult->LoadHistograms(Form("Multiplicity_%d", j));
 
-        TH1* eff2 = mult->GetEfficiency(3, eventType);
+        TH1* eff2 = mult->GetEfficiency(3, (AliMultiplicityCorrection::EventType) eventType);
 
         for (Int_t bin=1; bin<=eff->GetNbinsX(); bin++)
         {
           // TODO we could also do assymetric errors here
           Float_t deviation = 10.0 * TMath::Abs(eff->GetBinContent(bin) - eff2->GetBinContent(bin));
 
-          eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), deviation));
+          eff->SetBinError(bin, TMath::Max(eff->GetBinError(bin), (Double_t) deviation));
         }
       }
 
@@ -1529,7 +1676,7 @@ void EfficiencyComparison(Int_t eventType = 2)
         if (eff->GetBinContent(bin) > 0)
           Printf("Bin %d: Error: %.2f", bin, 100.0 * eff->GetBinError(bin) / eff->GetBinContent(bin));
 
-      effError = (TH1*) eff2->Clone("effError");
+      effError = (TH1*) eff->Clone("effError");
       effError->Reset();
 
       for (Int_t bin=2; bin<=eff->GetNbinsX(); bin++)
@@ -1612,8 +1759,8 @@ void ModelDependencyPlot()
   canvas->cd(1);
   gPad->SetLogy();
 
-  TH1* full = proj->ProjectionY("full2");
-  TH1* selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
+  full = proj->ProjectionY("full2");
+  selected = proj->ProjectionX("selected2", proj->GetYaxis()->FindBin(selectedMult), proj->GetYaxis()->FindBin(selectedMult));
 
   full->SetStats(kFALSE);
   full->GetXaxis()->SetRangeUser(0, 200);
@@ -1629,7 +1776,7 @@ void ModelDependencyPlot()
 
   legend->Draw();
 
-  TLine* line = new TLine(selectedMult, 5, selectedMult, yMax);
+  line = new TLine(selectedMult, 5, selectedMult, yMax);
   line->SetLineWidth(2);
   line->Draw();
 
@@ -1748,7 +1895,7 @@ void FitPt(const char* fileName = "firstplots100k_truept.root")
   /// genePt->GetXaxis()->GetBinCenter(x));
 
   genePt->GetXaxis()->SetRangeUser(0, 7.9);
-  genePt->GetYaxis()->SetTitle("a.u.");
+  //genePt->GetYaxis()->SetTitle("a.u.");
 
   TF1* func = new TF1("func", "[0]*TMath::Exp([1]*x)+[2]/(1+(x*[4])**[3])", 0.001, 100);
   //func->SetLineColor(2);
@@ -1810,23 +1957,48 @@ void FitPt(const char* fileName = "firstplots100k_truept.root")
 
   new TCanvas;
   genePt->DrawCopy("P");
-  func->SetRange(0.02, 8);
+  func->SetRange(0.02, 5);
   func->DrawCopy("SAME");
   gPad->SetLogy();
 
   genePt->Fit(func, "0", "", 0.2, 4);
 
-  TCanvas* canvas = new TCanvas("FitPt", "FitPt", 600, 400);
-  canvas->SetGridx();
-  canvas->SetGridy();
-  canvas->SetRightMargin(0.05);
-  canvas->SetTopMargin(0.05);
+  TCanvas* canvas = new TCanvas("FitPt", "FitPt", 800, 400);
+  canvas->Divide(2, 1);
+  canvas->cd(1);
+
+  gPad->SetGridx();
+  gPad->SetGridy();
+  gPad->SetLeftMargin(0.13);
+  gPad->SetRightMargin(0.05);
+  gPad->SetTopMargin(0.05);
 
+  genePt->GetXaxis()->SetRangeUser(0, 4.9);
+  genePt->GetYaxis()->SetRangeUser(1e-2, 1e4);
+  genePt->GetYaxis()->SetTitleOffset(1.4);
+  genePt->GetXaxis()->SetTitleOffset(1.1);
   genePt->DrawCopy("P");
-  func->SetRange(0.02, 8);
+  func->SetRange(0.02, 5);
+  func->DrawCopy("SAME");
+  gPad->SetLogy();
+
+  canvas->cd(2);
+
+  TH1* genePtClone = (TH1*) genePt->Clone("genePtClone");
+  genePtClone->Reset();
+  genePtClone->DrawCopy("P");
+
+  gPad->SetGridx();
+  gPad->SetGridy();
+  gPad->SetLeftMargin(0.13);
+  gPad->SetRightMargin(0.05);
+  gPad->SetTopMargin(0.05);
+
   func->DrawCopy("SAME");
   gPad->SetLogy();
 
+  canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+
   TH1* first = (TH1*) func->GetHistogram()->Clone("first");
 
   TCanvas* canvas2 = new TCanvas("FitPt2", "FitPt2", 600, 400);
@@ -1845,9 +2017,9 @@ void FitPt(const char* fileName = "firstplots100k_truept.root")
       func2->SetParameter(param, func2->GetParameter(param) * factor);
       //func2->Print();
 
-      canvas->cd();
+      canvas->cd(2);
       func2->SetLineWidth(1);
-      func2->SetLineColor(param + 2);
+      func2->SetLineColor(2);
       func2->DrawCopy("SAME");
 
       canvas2->cd();
@@ -1861,15 +2033,46 @@ void FitPt(const char* fileName = "firstplots100k_truept.root")
   }
 
   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+  canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
 
   file->Close();
 }
 
-void SystematicpT()
+void DrawSystematicpT()
+{
+  TFile* file = TFile::Open("SystematicpT.root");
+
+  TH1* mcHist2 = (TH1*) file->Get("mymc_unity");
+  TH1* result2 = (TH1*) file->Get("result_unity");
+
+  TH1* mcHist[12];
+  TH1* result[12];
+
+  Int_t nParams = 5;
+
+  for (Int_t id=0; id<nParams*2; ++id)
+  {
+    mcHist[id] = (TH1*) file->Get(Form("mymc_%d_%d.root", id / 2, id % 2));
+    result[id] = (TH1*) file->Get(Form("result_%d_%d.root", id / 2, id % 2));
+  }
+
+  DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
+
+  //DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
+
+  DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps", kTRUE, kTRUE);
+
+  //DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
+
+  // does not make sense: mc is different
+  //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
+}
+
+void SystematicpT(Bool_t chi2 = 1)
 {
   gSystem->Load("libPWG0base");
 
-  TFile::Open("ptspectrum400.root");
+  TFile::Open("ptspectrum900.root");
   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
   mult->LoadHistograms("Multiplicity");
 
@@ -1878,21 +2081,25 @@ void SystematicpT()
   TH1* mcHist[12];
   TH1* result[12];
 
-  Int_t nParams = 1;
+  Int_t nParams = 5;
 
   for (Int_t param=0; param<nParams; param++)
   {
     for (Int_t sign=0; sign<2; sign++)
     {
       // calculate result with systematic effect
-      TFile* file = TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
+      TFile::Open(Form("ptspectrum100_%d_%d.root", param, sign));
       mult2->LoadHistograms("Multiplicity");
 
       mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 
-      mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
-      mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
-      //mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+      if (chi2)
+      {
+        mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+        mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+      }
+      else
+        mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0);
 
       Int_t id = param * 2 + sign;
 
@@ -1907,39 +2114,63 @@ void SystematicpT()
   // calculate normal result
   TFile::Open("ptspectrum100_1.root");
   mult2->LoadHistograms("Multiplicity");
-  TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc2");
+  TH1* mcHist2 = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc_unity");
 
   mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
 
-  mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
-  //mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  if (chi2)
+  {
+    mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  }
+  else
+    mult->ApplyBayesianMethod(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
 
-  TH1* result2 = mult->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
+  TH1* result2 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result_unity");
 
-  DrawResultRatio(mcHist2, result2, "SystematicpT_OK.eps");
+  TFile* file = TFile::Open("SystematicpT.root", "RECREATE");
+  mcHist2->Write();
+  result2->Write();
+  for (Int_t id=0; id<nParams*2; ++id)
+  {
+    mcHist[id]->Write();
+    result[id]->Write();
+  }
+  file->Close();
+
+  DrawSystematicpT();
+}
 
-  DrawRatioDeduct(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
+void SystematicpTCutOff()
+{
+  gSystem->Load("libPWG0base");
+  SetTPC();
 
-  TCanvas* canvas = DrawRatio(nParams*2, mcHist, result, "SystematicpT_Ratios.eps");
+  TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_unity.root");
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  mult->LoadHistograms("Multiplicity");
 
-  TFile::Open("bayesianUncertainty_pt_400_100.root");
-  TH1* errorHist = (TH1*) gFile->Get("errorBoth");
-  errorHist->SetLineColor(1);
-  errorHist->SetLineWidth(2);
-  TH1* errorHist2 = (TH1*) errorHist->Clone("errorHist2");
-  for (Int_t i=1; i<=errorHist->GetNbinsX(); i++)
-  {
-    errorHist->SetBinContent(i, errorHist->GetBinContent(i) + 1);
-    errorHist2->SetBinContent(i, 1 - errorHist2->GetBinContent(i));
-  }
+  TFile::Open(measuredFile);
+  AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity2", "Multiplicity2");
+  mult2->LoadHistograms("Multiplicity");
 
-  errorHist->DrawCopy("SAME");
-  errorHist2->DrawCopy("SAME");
+  // "normal" result
+  mult->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+  mult->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
 
-  canvas->SaveAs(canvas->GetName());
+  TH1* mcHist = mult2->GetMultiplicityVtx(etaRange)->ProjectionY("mymc");
+  TH1* result1 = (TH1*) mult->GetMultiplicityESDCorrected(etaRange)->Clone("result1");
 
-  DrawRatioDeductSmooth(mcHist2, result2, nParams*2, mcHist, result, "SystematicpT_Summary.eps");
+  // change of pt spectrum (down)
+  TFile::Open("multiplicityMC_TPC_1.3M_syst_pt_red.root");
+  AliMultiplicityCorrection* mult3 = new AliMultiplicityCorrection("Multiplicity3", "Multiplicity3");
+  mult3->LoadHistograms("Multiplicity");
 
-  // does not make sense: mc is different
-  //Draw2ResultRatio(mcHist, result1, result2, "SystematicpT.eps");
+  mult3->SetMultiplicityESD(etaRange, mult2->GetMultiplicityESD(etaRange));
+  mult3->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+  mult3->ApplyMinuitFit(etaRange, kFALSE, AliMultiplicityCorrection::kTrVtx);
+
+  TH1* result2 = (TH1*) mult3->GetMultiplicityESDCorrected(etaRange)->Clone("result2");
+
+  Draw2ResultRatio(mcHist, result1, result2, "SystematicLowEfficiency.eps");
 }
index 0a26fce..b3b6107 100644 (file)
 void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "")
 {
   if (aProof)
+  {
+    TProof::Mgr("lxb6046")->SetROOTVersion("vHEAD-07Jun2007b_dbg");
     connectProof("lxb6046");
+  }
 
   TString libraries("libEG;libGeom;libESD;libPWG0base");
   TString packages("PWG0base");
@@ -38,6 +41,16 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_
   TList inputList;
   inputList.Add(esdTrackCuts);
 
+  // pt study
+  if (TString(option).Contains("pt-spectrum-func"))
+  {
+    //TF1* func = new TF1("func", "0.7 + x", 0, 0.3);
+    //TF1* func = new TF1("func", "1.3 - x", 0, 0.3);
+    TF1* func = new TF1("func", "1", 0, 0.3);
+    new TCanvas; func->Draw();
+    inputList.Add(func->GetHistogram()->Clone("pt-spectrum"));
+  }
+
   TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE);
 
   TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector");
@@ -131,7 +144,7 @@ void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fol
   }
   else
   {
-    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.2, 100);
     mult->DrawComparison("Bayesian", histID, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
   }
 
@@ -1772,9 +1785,10 @@ void BuildResponseFromTree(const char* fileName, const char* target)
 
   Int_t tracks = 0; // control variables
   Int_t noLabel = 0;
+  Int_t secondaries = 0;
   Int_t doubleCount = 0;
 
-  for (Int_t num = 0; num < 7; num++)
+  for (Int_t num = 0; num < 1; num++)
   {
     AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(Form("Multiplicity_%d", num), Form("Multiplicity_%d", num));
 
@@ -1812,19 +1826,20 @@ void BuildResponseFromTree(const char* fileName, const char* target)
       tracks += f[9];
       noLabel += f[9];
 
-      // add the double counted
-      tracks += f[10];
-      doubleCount += f[10];
+      // secondaries are already part of meas!
+      secondaries += f[10];
+
+      // double counted are already part of meas!
+      doubleCount += f[11];
 
-      // TODO fake, add the ones w/o label and the double counted to check the method
+      // ones w/o label are added without weighting to allow comparison to default analysis. however this is only valid when their fraction is low enough!
       meas += f[9];
-      meas += f[10];
 
       //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
 
-      fMultiplicity->FillCorrection(f[0], 0, 0, 0, gene, 0, 0, 0, 0, meas);
-      fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 0, 0, 0, gene, 0);
-      fMultiplicity->FillMeasured(f[0], 0, 0, 0, meas);
+      fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, 0, meas, meas, meas, meas);
+      fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, gene, gene, gene, gene, 0);
+      fMultiplicity->FillMeasured(f[0], meas, meas, meas, meas);
     }
 
     //fMultiplicity->DrawHistograms();
@@ -1834,7 +1849,11 @@ void BuildResponseFromTree(const char* fileName, const char* target)
     file->Close();
 
     if (num == 0)
-      Printf("%d total tracks, %d w/o label = %.2f %%, %d double counted = %.2f %% ", tracks, noLabel, 100.0 * noLabel / tracks, doubleCount, 100.0 * doubleCount / tracks);
+    {
+      Printf("%d total tracks, %d w/o label = %.2f %%, %d double counted = %.2f %%, secondaries = %.2f %%", tracks, noLabel, 100.0 * noLabel / tracks, doubleCount, 100.0 * doubleCount / tracks, 100.0 * secondaries / tracks);
+      if ((Float_t) noLabel / tracks > 0.02)
+        Printf("WARNING: More than 2%% of tracks without label, this might bias the study!");
+    }
   }
 }