AliPWG0depHelper: function to find mother among the primaries
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 11 Jun 2007 16:52:23 +0000 (16:52 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 11 Jun 2007 16:52:23 +0000 (16:52 +0000)
Multiplicity study: extended in many ways + adding plots for mult. note

12 files changed:
PWG0/AliCorrection.cxx
PWG0/AliPWG0depHelper.cxx
PWG0/AliPWG0depHelper.h
PWG0/PWG0selectorsLinkDef.h
PWG0/dNdEta/AliMultiplicityCorrection.cxx
PWG0/dNdEta/AliMultiplicityCorrection.h
PWG0/dNdEta/AliMultiplicityMCSelector.cxx
PWG0/dNdEta/AliMultiplicityMCSelector.h
PWG0/dNdEta/AlidNdEtaSystematicsSelector.cxx
PWG0/dNdEta/runMultiplicitySelector.C
PWG0/dNdEta/testAnalysis2.C
PWG0/libPWG0selectors.pkg

index a29dd0d..319d17e 100644 (file)
@@ -41,8 +41,8 @@ AliCorrection::AliCorrection(const Char_t* name, const Char_t* title) : TNamed(n
   //Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
   //Float_t binLimitsVtx[] = {-20,-19,-18,-17,-16,-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};
   Float_t binLimitsVtx[] = {-20,-15,-10,-8,-6,-4,-2,0,2,4,6,8,10,15,20};
-  Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
-                           0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
+  Float_t binLimitsEta[] = {-2,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
+//  Float_t binLimitsEta[] = { -7.0, -6.9, -6.8, -6.7, -6.6, -6.5, -6.4, -6.3, -6.2, -6.1, -6.0, -5.9, -5.8, -5.7, -5.6, -5.5, -5.4, -5.3, -5.2, -5.1, -5.0, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4.0, -3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3.0, -2.9, -2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, -0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0 };
 
   TH3F* dummyBinning = new TH3F("","",14, binLimitsVtx, 40, binLimitsEta , 28, binLimitsPt);
 
index 8a5dd88..d721aba 100644 (file)
 /* $Id$ */
 
 #include <TList.h>
+#include <TParticle.h>
 
 #include <AliPWG0depHelper.h>
 
 #include <AliHeader.h>
+#include <AliStack.h>
+#include <AliLog.h>
 
 #include <AliGenEventHeader.h>
 #include <AliGenPythiaEventHeader.h>
@@ -29,7 +32,7 @@
 ClassImp(AliPWG0depHelper)
 
 //____________________________________________________________________
-const Int_t AliPWG0depHelper::GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug) {
+Int_t AliPWG0depHelper::GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug) {
   //
   // get the process type of the event.
   //
@@ -68,3 +71,44 @@ const Int_t AliPWG0depHelper::GetPythiaEventProcessType(AliHeader* aHeader, Bool
 
   return pythiaGenHeader->ProcessType();
 }
+
+//____________________________________________________________________
+TParticle* AliPWG0depHelper::FindPrimaryMother(AliStack* stack, Int_t label)
+{
+  //
+  // Finds the first mother among the primary particles of the particle identified by <label>,
+  // i.e. the primary that "caused" this particle
+  //
+
+  Int_t nPrim  = stack->GetNprimary();
+  TParticle* mother = stack->Particle(label);
+
+  if (!mother)
+  {
+    AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not retrieve particle with label %d.", label));
+    return 0;
+  }
+
+  while (label >= nPrim)
+  {
+    //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
+
+    if (mother->GetMother(0) == -1)
+    {
+      AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
+      mother = 0;
+      break;
+    }
+
+    label = mother->GetMother(0);
+
+    mother = stack->Particle(label);
+    if (!mother)
+    {
+      AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
+      break;
+    }
+  }
+
+  return mother;
+}
index 62f7b24..bb07f68 100644 (file)
@@ -8,11 +8,14 @@
 // static helper functions that depend on more than ESD
 
 class AliHeader;
+class TParticle;
+class AliStack;
 
 class AliPWG0depHelper : public TObject
 {
   public:
-    static const Int_t GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
+    static Int_t GetPythiaEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE);
+    static TParticle* FindPrimaryMother(AliStack* stack, Int_t label);
 
   protected:
     ClassDef(AliPWG0depHelper, 0)
index dad9831..3f55c79 100644 (file)
@@ -18,5 +18,6 @@
 #pragma link C++ class AliROCESDAnalysisSelector+;
 #pragma link C++ class AliROCRawAnalysisSelector+;
 #pragma link C++ class AliROCClusterAnalysisSelector+;
+#pragma link C++ class AliHighMultiplicitySelector+;
 
 #endif
index bdb827f..45bc241 100644 (file)
 #include <TF1.h>
 #include <TMath.h>
 #include <TCollection.h>
+#include <TLegend.h>
+#include <TLine.h>
 
 ClassImp(AliMultiplicityCorrection)
 
-const Int_t AliMultiplicityCorrection::fgMaxParams = 200;
+const Int_t AliMultiplicityCorrection::fgMaxInput  = 250;  // bins in measured histogram
+const Int_t AliMultiplicityCorrection::fgMaxParams = 250;  // bins in unfolded histogram = number of fit params
+
 TH1* AliMultiplicityCorrection::fCurrentESD = 0;
 TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
 TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
-TMatrixF* AliMultiplicityCorrection::fCorrelationMatrix = 0;
-TMatrixF* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
-TVectorF* AliMultiplicityCorrection::fCurrentESDVector = 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 = 1e4;
+Float_t AliMultiplicityCorrection::fRegularizationWeight = 5000;
 TF1* AliMultiplicityCorrection::fNBD = 0;
 
 //____________________________________________________________________
@@ -103,8 +108,8 @@ AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const C
   #define VTXBINNING 10, binLimitsVtx
   #define NBINNING fgMaxParams, binLimitsN*/
 
-  #define NBINNING 251, -0.5, 250.5
-  #define VTXBINNING 10, -10, 10
+  #define NBINNING 501, -0.5, 500.5
+  #define VTXBINNING 1, -10, 10
 
   for (Int_t i = 0; i < kESDHists; ++i)
     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
@@ -352,7 +357,7 @@ void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, I
 }
 
 //____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
 {
   // homogenity term for minuit fitting
   // pure function of the parameters
@@ -360,20 +365,24 @@ Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
 
   Double_t chi2 = 0;
 
-  for (Int_t i=1; i<fgMaxParams; ++i)
+  // ignore the first bin here. on purpose...
+  for (Int_t i=2; i<fgMaxParams; ++i)
   {
     Double_t right  = params[i];
     Double_t left   = params[i-1];
 
-    Double_t diff = (right - left);
-    chi2 += diff * diff;
+    if (right != 0)
+    {
+      Double_t diff = 1 - left / right;
+      chi2 += diff * diff;
+    }
   }
 
-  return chi2 * 100;
+  return chi2 / 100.0;
 }
 
 //____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
 {
   // homogenity term for minuit fitting
   // pure function of the parameters
@@ -381,7 +390,8 @@ Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
 
   Double_t chi2 = 0;
 
-  for (Int_t i=2; i<fgMaxParams; ++i)
+  // ignore the first bin here. on purpose...
+  for (Int_t i=2+1; i<fgMaxParams; ++i)
   {
     if (params[i-1] == 0)
       continue;
@@ -402,7 +412,7 @@ Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
 }
 
 //____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTest(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationTest(TVectorD& params)
 {
   // homogenity term for minuit fitting
   // pure function of the parameters
@@ -428,7 +438,7 @@ Double_t AliMultiplicityCorrection::RegularizationTest(Double_t *params)
 }
 
 //____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
 {
   // homogenity term for minuit fitting
   // pure function of the parameters
@@ -437,7 +447,8 @@ Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *param
 
   Double_t chi2 = 0;
 
-  for (Int_t i=2; i<fgMaxParams; ++i)
+  // ignore the first bin here. on purpose...
+  for (Int_t i=3; i<fgMaxParams; ++i)
   {
     Double_t right  = params[i];
     Double_t middle = params[i-1];
@@ -446,17 +457,16 @@ Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *param
     Double_t der1 = (right - middle);
     Double_t der2 = (middle - left);
 
-    Double_t secDer = (der1 - der2);
+    Double_t diff = (der1 - der2);
 
-    // square and weight with the bin width
-    chi2 += secDer * secDer;
+    chi2 += diff * diff;
   }
 
-  return chi2 * 100;
+  return chi2 * 1e4;
 }
 
 //____________________________________________________________________
-Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params)
+Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
 {
   // homogenity term for minuit fitting
   // pure function of the parameters
@@ -464,18 +474,19 @@ Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params)
   // The method of reduced cross-entropy (M. Schmelling 1993)
 
   Double_t paramSum = 0;
-  for (Int_t i=0; i<fgMaxParams; ++i)
+  // ignore the first bin here. on purpose...
+  for (Int_t i=1; i<fgMaxParams; ++i)
     paramSum += params[i];
 
   Double_t chi2 = 0;
-  for (Int_t i=0; i<fgMaxParams; ++i)
+  for (Int_t i=1; i<fgMaxParams; ++i)
   {
     Double_t tmp = params[i] / paramSum;
-    if (tmp > 0)
-      chi2 += tmp * log(tmp);
+    if (tmp > 0 && (*fEntropyAPriori)[i] > 0)
+      chi2 += tmp * log(tmp / (*fEntropyAPriori)[i]);
   }
 
-  return chi2;
+  return 10.0 + chi2;
 }
 
 //____________________________________________________________________
@@ -486,6 +497,9 @@ void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Dou
   // does: nbd
   // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
   // func->SetParNames("scaling", "averagen", "k");
+  // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
+  // func->SetParLimits(1, 0.001, 1000);
+  // func->SetParLimits(2, 0.001, 1000);
   //
 
   fNBD->SetParameters(params[0], params[1], params[2]);
@@ -508,51 +522,54 @@ void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& c
   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
   //
 
-  static Int_t callCount = 0;
-
-  Double_t chi2FromFit = 0;
-
   // d
-  TVectorF paramsVector(fgMaxParams);
+  TVectorD paramsVector(fgMaxParams);
   for (Int_t i=0; i<fgMaxParams; ++i)
-    paramsVector[i] = params[i];
+    paramsVector[i] = params[i] * params[i];
+
+  // calculate penalty factor
+  Double_t penaltyVal = 0;
+  switch (fRegularizationType)
+  {
+    case kNone:       break;
+    case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
+    case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
+    case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
+    case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
+    case kTest:       penaltyVal = RegularizationTest(paramsVector); break;
+  }
+
+  //if (penaltyVal > 1e10)
+  //  paramsVector2.Print();
+
+  penaltyVal *= fRegularizationWeight;
 
   // Ad
-  paramsVector = (*fCorrelationMatrix) * paramsVector;
+  TVectorD measGuessVector(fgMaxInput);
+  measGuessVector = (*fCorrelationMatrix) * paramsVector;
 
   // Ad - m
-  paramsVector -= (*fCurrentESDVector);
+  measGuessVector -= (*fCurrentESDVector);
 
-  TVectorF copy(paramsVector);
+  TVectorD copy(measGuessVector);
 
   // (Ad - m) W
-  paramsVector *= (*fCorrelationCovarianceMatrix);
+  measGuessVector *= (*fCorrelationCovarianceMatrix);
 
-  // (Ad - m) W (Ad - m)
-  chi2FromFit = paramsVector * copy;
-
-  Double_t penaltyVal = 0;
+  //measGuessVector.Print();
 
-  switch (fRegularizationType)
-  {
-    case kNone:       break;
-    case kPol0:       penaltyVal = RegularizationPol0(params); break;
-    case kPol1:       penaltyVal = RegularizationPol1(params); break;
-    case kCurvature:  penaltyVal = RegularizationTotalCurvature(params); break;
-    case kEntropy:    penaltyVal = RegularizationEntropy(params); break;
-    case kTest:    penaltyVal = RegularizationTest(params); break;
-  }
-
-  penaltyVal *= fRegularizationWeight;
+  // (Ad - m) W (Ad - m)
+  Double_t chi2FromFit = measGuessVector * copy * 1e6;
 
   chi2 = chi2FromFit + penaltyVal;
 
-  if ((callCount++ % 1000) == 0)
-    printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
+  static Int_t callCount = 0;
+  if ((callCount++ % 10000) == 0)
+    printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
+void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
 {
   //
   // fills fCurrentESD, fCurrentCorrelation
@@ -577,8 +594,73 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
   }
 
   fCurrentCorrelation = hist->Project3D("zy");
+  //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
+  //fMultiplicityESDCorrected[correlationID]->Rebin(2);
   fCurrentCorrelation->Sumw2();
 
+  if (createBigBin)
+  {
+    Int_t maxBin = 0;
+    for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
+    {
+      if (fCurrentESD->GetBinContent(i) <= 5)
+      {
+        maxBin = i;
+        break;
+      }
+    }
+
+    if (maxBin > 0)
+    {
+      TCanvas* canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
+      canvas->Divide(2, 2);
+
+      canvas->cd(1);
+      fCurrentESD->DrawCopy();
+      gPad->SetLogy();
+
+      canvas->cd(2);
+      fCurrentCorrelation->DrawCopy("COLZ");
+
+      printf("Bin limit in measured spectrum is %d.\n", maxBin);
+      fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
+      for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i)
+      {
+        fCurrentESD->SetBinContent(i, 0);
+        fCurrentESD->SetBinError(i, 0);
+      }
+      // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+      fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin)));
+
+      printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin));
+
+      for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
+      {
+        fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY()));
+        // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+        fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin)));
+
+        for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+        {
+          fCurrentCorrelation->SetBinContent(i, j, 0);
+          fCurrentCorrelation->SetBinError(i, j, 0);
+        }
+      }
+
+      printf("Adjusted correlation matrix!\n");
+
+      canvas->cd(3);
+      fCurrentESD->DrawCopy();
+      gPad->SetLogy();
+
+      canvas->cd(4);
+      fCurrentCorrelation->DrawCopy("COLZ");
+    }
+  }
+
+  //normalize ESD
+  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+
   fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency");
   TH2* divisor = 0;
   switch (eventType)
@@ -588,10 +670,21 @@ void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullP
     case kINEL: divisor = fMultiplicityINEL[inputRange]; break;
   }
   fCurrentEfficiency->Divide(divisor->ProjectionY());
+  //fCurrentEfficiency->Rebin(2);
+  //fCurrentEfficiency->Scale(0.5);
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
+void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight)
+{
+  fRegularizationType = type;
+  fRegularizationWeight = 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)
 {
   //
   // correct spectrum using minuit chi2 method
@@ -603,11 +696,19 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE); // TODO FAKE kTRUE
 
-  fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
-  fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
-  fCurrentESDVector = new TVectorF(fgMaxParams);
+  fCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
+  fCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
+  fCurrentESDVector = new TVectorD(fgMaxInput);
+  fEntropyAPriori = new TVectorD(fgMaxParams);
+
+  /*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 correction for given nPart
   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
@@ -630,61 +731,83 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
 
-      if (i <= fgMaxParams && j <= fgMaxParams)
+      if (i <= fgMaxParams && j <= fgMaxInput)
         (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
     }
 
     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
   }
 
-  // small hack to get around charge conservation for full phase space ;-)
-  /*if (fullPhaseSpace)
-  {
-    for (Int_t i=2; i<=50; i+=2)
-      for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
-        fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
-  }*/
-
   // Initialize TMinuit via generic fitter interface
   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
+  Double_t arglist[100];
+  arglist[0] = 0;
+  minuit->ExecuteCommand("SET PRINT", arglist, 1);
 
   minuit->SetFCN(MinuitFitFunction);
 
-  static Double_t results[fgMaxParams];
+  for (Int_t i=0; i<fgMaxParams; i++)
+    (*fEntropyAPriori)[i] = 1;
 
   if (inputDist)
   {
     printf("Using different starting conditions...\n");
     new TCanvas;
     inputDist->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);
   }
   else
     inputDist = fCurrentESD;
 
-  // normalize ESD
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
-  inputDist->Scale(1.0 / inputDist->Integral());
 
-  Float_t minStartValue = 1e-3;
+  //Float_t minStartValue = 0; //1e-3;
 
+  //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];
   for (Int_t i=0; i<fgMaxParams; ++i)
   {
     results[i] = inputDist->GetBinContent(i+1);
+
     if (check)
       results[i] = proj->GetBinContent(i+1);
 
-    printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
+    // minimum value
+    results[i] = TMath::Max(results[i], 1e-3);
+
+    Float_t stepSize = 0.1;
+
+    // minuit sees squared values to prevent it from going negative...
+    results[i] = TMath::Sqrt(results[i]);
+    //stepSize /= results[i]; // keep relative step size
+
+    minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
+  }
+  // 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)
+  {
+    //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) = 1.0 / fCurrentESD->GetBinError(i+1);
+      (*fCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / fCurrentESD->GetBinError(i+1) / fCurrentESD->GetBinError(i+1);
+
+    if ((*fCorrelationCovarianceMatrix)(i, i) > 1e7)
+      (*fCorrelationCovarianceMatrix)(i, i) = 0;
 
-    minuit->SetParameter(i, Form("param%d", i), ((results[i] > minStartValue) ? results[i] : minStartValue), 0.1, 0, 1);
+    //printf("%d --> %e\n", i, (*fCorrelationCovarianceMatrix)(i, i));
   }
-  // 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, 1);
 
   Int_t dummy;
   Double_t chi2;
@@ -692,12 +815,16 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
   printf("Chi2 of initial parameters is = %f\n", chi2);
 
   if (check)
-    return;
-
-  Double_t arglist[100];
-  arglist[0] = 0;
-  minuit->ExecuteCommand("SET PRINT", arglist, 1);
-  minuit->ExecuteCommand("MIGRAD", arglist, 1);
+    return -1;
+
+  // first param is number of iterations, second is precision....
+  arglist[0] = 1e6;
+  //arglist[1] = 1e-5;
+  //minuit->ExecuteCommand("SCAN", arglist, 0);
+  Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
+  printf("MINUIT status is %d\n", status);
+  //minuit->ExecuteCommand("MIGRAD", arglist, 1);
+  //minuit->ExecuteCommand("MIGRAD", arglist, 1);
   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
   //minuit->ExecuteCommand("MINOS", arglist, 0);
 
@@ -705,8 +832,9 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
   {
     if (fCurrentEfficiency->GetBinContent(i+1) > 0)
     {
-      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
-      fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) / fCurrentEfficiency->GetBinContent(i+1));
+      fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / fCurrentEfficiency->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));
     }
     else
     {
@@ -714,6 +842,8 @@ void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhas
       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
     }
   }
+
+  return status;
 }
 
 //____________________________________________________________________
@@ -725,11 +855,11 @@ void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSp
 
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
+  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
 
-  fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
-  fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
-  fCurrentESDVector = new TVectorF(fgMaxParams);
+  fCorrelationMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
+  fCorrelationCovarianceMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
+  fCurrentESDVector = new TVectorD(fgMaxParams);
 
   // normalize correction for given nPart
   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
@@ -773,7 +903,7 @@ void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSp
   Double_t arglist[100];
   arglist[0] = 0;
   minuit->ExecuteCommand("SET PRINT", arglist, 1);
-  minuit->ExecuteCommand("MIGRAD", arglist, 1);
+  minuit->ExecuteCommand("MIGRAD", arglist, 0);
   //minuit->ExecuteCommand("MINOS", arglist, 0);
 
   fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
@@ -864,48 +994,107 @@ void AliMultiplicityCorrection::DrawHistograms()
 }
 
 //____________________________________________________________________
-void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist)
+void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
 {
+  //mcHist->Rebin(2);
+
   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 
   TString tmpStr;
   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
 
+  if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
+  {
+    printf("ERROR. Unfolded histogram is empty\n");
+    return;
+  }
+
+  //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
+  fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
+  fCurrentESD->Sumw2();
+  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+
+  // normalize unfolded result to 1
+  fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
+
+  //fCurrentESD->Scale(mcHist->Integral(2, 200));
+
+  //new TCanvas;
+  /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
+  ratio->Divide(mcHist);
+  ratio->Draw("HIST");
+  ratio->Fit("pol0", "W0", "", 20, 120);
+  Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
+  delete ratio;
+  mcHist->Scale(scalingFactor);*/
+
+  mcHist->Scale(1.0 / mcHist->Integral());
+
   // calculate residual
 
   // for that we convolute the response matrix with the gathered result
   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
-  TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId, !normalizeESD);
-  TH1* proj2 = convoluted->ProjectionY("proj2", -1, -1, "e");
-  NormalizeToBinWidth(proj2);
+  TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
+  TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
+  if (convolutedProj->Integral() > 0)
+  {
+    convolutedProj->Scale(1.0 / convolutedProj->Integral());
+  }
+  else
+    printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
+
+  //NormalizeToBinWidth(proj2);
 
-  TH1* residual = (TH1*) proj2->Clone("residual");
-  residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / measured");
+  TH1* residual = (TH1*) convolutedProj->Clone("residual");
+  residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
 
   residual->Add(fCurrentESD, -1);
-  residual->Divide(residual, fCurrentESD, 1, 1, "B");
+  //residual->Divide(residual, fCurrentESD, 1, 1, "B");
+
+  TH1* residualHist = new TH1F("residualHist", "residualHist", 50, -5, 5);
 
   // TODO fix errors
-  for (Int_t i=1; i<residual->GetNbinsX(); ++i)
+  Float_t chi2 = 0;
+  for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
   {
-    proj2->SetBinError(i, 0);
+    if (fCurrentESD->GetBinError(i) > 0)
+    {
+      Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
+      if (i > 1)
+        chi2 += value * value;
+      residual->SetBinContent(i, value);
+      residualHist->Fill(value);
+    }
+    else
+    {
+      //printf("Residual bin %d set to 0\n", i);
+      residual->SetBinContent(i, 0);
+    }
+    convolutedProj->SetBinError(i, 0);
     residual->SetBinError(i, 0);
-  }
 
-  // normalize mc to 1
-  //mcHist->Scale(1.0 / mcHist->Integral(1, 200));
+    if (i == 200)
+      fLastChi2Residuals = chi2;
+  }
 
-  // normalize unfolded esd to 1
-  //fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral(1, 200));
+  residualHist->Fit("gaus", "N");
+  delete residualHist;
 
-  // normalize unfolded result to mc hist
-  fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
-  fMultiplicityESDCorrected[esdCorrId]->Scale(mcHist->Integral(1, 200));
+  printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
 
-  TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 800);
-  canvas1->Divide(2, 2);
+  TCanvas* canvas1 = 0;
+  if (simple)
+  {
+    canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
+    canvas1->Divide(2, 1);
+  }
+  else
+  {
+    canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
+    canvas1->Divide(2, 3);
+  }
 
   canvas1->cd(1);
   TH1* proj = (TH1*) mcHist->Clone("proj");
@@ -915,76 +1104,162 @@ void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRang
     NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
 
   proj->GetXaxis()->SetRangeUser(0, 200);
+  proj->SetTitle(";true multiplicity;Entries");
+  proj->SetStats(kFALSE);
   proj->DrawCopy("HIST");
   gPad->SetLogy();
 
   //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
-  fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST");
+  fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
 
-  Float_t chi2 = 0;
+  TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
+  legend->AddEntry(proj, "true distribution");
+  legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
+  legend->SetFillColor(0);
+  legend->Draw();
+  // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
+
+  canvas1->cd(2);
+
+  gPad->SetLogy();
+  fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
+  //fCurrentESD->SetLineColor(2);
+  fCurrentESD->SetTitle(";measured multiplicity;Entries");
+  fCurrentESD->SetStats(kFALSE);
+  fCurrentESD->DrawCopy("HIST E");
+
+  convolutedProj->SetLineColor(2);
+  //proj2->SetMarkerColor(2);
+  //proj2->SetMarkerStyle(5);
+  convolutedProj->DrawCopy("HIST SAME");
+
+  legend = new TLegend(0.3, 0.8, 0.93, 0.93);
+  legend->AddEntry(fCurrentESD, "measured distribution");
+  legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
+  legend->SetFillColor(0);
+  legend->Draw();
+
+  TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
+  diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
+
+  /*Float_t chi2 = 0;
+  Float_t chi = 0;
   fLastChi2MCLimit = 0;
+  Int_t limit = 0;
   for (Int_t i=2; i<=200; ++i)
   {
     if (proj->GetBinContent(i) != 0)
     {
       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
       chi2 += value * value;
+      chi += TMath::Abs(value);
 
-      if (chi2 < 0.1)
-        fLastChi2MCLimit = i;
-
-      if (i == 100)
-        fLastChi2MC = chi2;
-    }
-  }
+      //printf("%d %f\n", i, chi);
 
-  printf("Difference (from MC) is %f for bin 2-100. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
+      if (chi2 < 0.2)
+        fLastChi2MCLimit = i;
 
-  canvas1->cd(2);
+      if (chi < 3)
+        limit = i;
 
-  NormalizeToBinWidth(fCurrentESD);
-  gPad->SetLogy();
-  fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
-  //fCurrentESD->SetLineColor(2);
-  fCurrentESD->DrawCopy("HIST");
-
-  proj2->SetLineColor(2);
-  //proj2->SetMarkerColor(2);
-  //proj2->SetMarkerStyle(5);
-  proj2->DrawCopy("HIST SAME");
+    }
+  }*/
 
   chi2 = 0;
-  for (Int_t i=2; i<=100; ++i)
+  Float_t chi = 0;
+  Int_t limit = 0;
+  for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
   {
-    if (fCurrentESD->GetBinContent(i) != 0)
+    if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
     {
-      Float_t value = (proj2->GetBinContent(i) - fCurrentESD->GetBinContent(i)) / fCurrentESD->GetBinContent(i);
-      chi2 += value * value;
+      Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
+      if (value > 1e8)
+        value = 1e8; //prevent arithmetic exception
+      else if (value < -1e8)
+        value = -1e8;
+      if (i > 1)
+      {
+        chi2 += value * value;
+        chi += TMath::Abs(value);
+      }
+      diffMCUnfolded->SetBinContent(i, value);
     }
+    else
+    {
+      //printf("diffMCUnfolded bin %d set to 0\n", i);
+      diffMCUnfolded->SetBinContent(i, 0);
+    }
+    if (chi2 < 1000)
+      fLastChi2MCLimit = i;
+    if (chi < 1000)
+      limit = i;
+    if (i == 150)
+      fLastChi2MC = chi2;
   }
 
-  printf("Difference (Residuals) is %f for bin 2-100\n", chi2);
-  fLastChi2Residuals = chi2;
+  printf("limits %d %d\n", fLastChi2MCLimit, limit);
+  fLastChi2MCLimit = limit;
+
+  printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);
 
-  canvas1->cd(3);
-  TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
-  clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
-  clone->SetTitle("Ratio;Npart;MC/Unfolded Measured");
-  clone->GetYaxis()->SetRangeUser(0.8, 1.2);
-  clone->GetXaxis()->SetRangeUser(0, 200);
-  clone->DrawCopy("HIST");
+  if (!simple)
+  {
+    canvas1->cd(3);
+
+    diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e");
+    //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
+    diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
+    diffMCUnfolded->DrawCopy("HIST");
 
-  /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
-  legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
-  legend->AddEntry(fMultiplicityMC, "MC");
-  legend->AddEntry(fMultiplicityESD, "ESD");
-  legend->Draw();*/
+    TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
+    for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
+      fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
 
-  canvas1->cd(4);
-  residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
-  residual->GetXaxis()->SetRangeUser(0, 200);
-  residual->DrawCopy();
+    new TCanvas;
+    fluctuation->Draw();
+
+    /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
+    legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
+    legend->AddEntry(fMultiplicityMC, "MC");
+    legend->AddEntry(fMultiplicityESD, "ESD");
+    legend->Draw();*/
+
+    canvas1->cd(4);
+    //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
+    residual->GetXaxis()->SetRangeUser(0, 200);
+    residual->DrawCopy();
+
+    canvas1->cd(5);
+
+    TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
+    ratio->Divide(mcHist);
+    ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
+    ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
+    ratio->GetXaxis()->SetRangeUser(0, 200);
+    ratio->SetStats(kFALSE);
+    ratio->Draw("HIST");
+
+    Double_t ratioChi2 = 0;
+    fLastChi2MCLimit = 0;
+    for (Int_t i=2; i<=150; ++i)
+    {
+      Float_t value = ratio->GetBinContent(i) - 1;
+      if (value > 1e8)
+        value = 1e8; //prevent arithmetic exception
+      else if (value < -1e8)
+        value = -1e8;
+
+      ratioChi2 += value * value;
+
+      if (ratioChi2 < 0.1)
+        fLastChi2MCLimit = i;
+    }
+
+    printf("Sum over (ratio-1)^2 (2..150) is %f\n", ratioChi2);
+    // TODO FAKE
+    fLastChi2MC = ratioChi2;
+  }
 
   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
 }
@@ -1031,10 +1306,17 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
 
-  // normalize ESD
-  fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
+  // smooth efficiency
+  //TH1* tmp = (TH1*) fCurrentEfficiency->Clone("eff_clone");
+  //for (Int_t i=2; i<fCurrentEfficiency->GetNbinsX(); ++i)
+  //  fCurrentEfficiency->SetBinContent(i, (tmp->GetBinContent(i-1) + tmp->GetBinContent(i) + tmp->GetBinContent(i+1)) / 3);
+
+  // set efficiency to 1 above 150
+  // FAKE TEST
+  //for (Int_t i=150; i<=fCurrentEfficiency->GetNbinsX(); ++i)
+  //  fCurrentEfficiency->SetBinContent(i, 1);
 
   // normalize correction for given nPart
   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
@@ -1063,6 +1345,27 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
     }
   }
 
+  // normalize nTrack
+  /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
+  {
+    // with this it is normalized to 1
+    Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, 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);
+      }
+      else
+      {
+        fCurrentCorrelation->SetBinContent(i, j, 0);
+        fCurrentCorrelation->SetBinError(i, j, 0);
+      }
+    }
+  }*/
+
   // smooth input spectrum
   /*
   TH1* esdClone = (TH1*) fCurrentESD->Clone("esdClone");
@@ -1077,6 +1380,18 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
   */
 
+  /*new TCanvas;
+  fCurrentEfficiency->Draw();
+
+  new TCanvas;
+  fCurrentCorrelation->DrawCopy("COLZ");
+
+  new TCanvas;
+  ((TH2*) fCurrentCorrelation)->ProjectionX()->DrawCopy();
+
+  new TCanvas;
+  ((TH2*) fCurrentCorrelation)->ProjectionY()->DrawCopy();*/
+
   // pick prior distribution
   TH1* hPrior = 0;
   if (inputDist)
@@ -1103,6 +1418,8 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
 
   TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
 
+  const Int_t kStartBin = 1;
+
   // unfold...
   for (Int_t i=0; i<nIterations; i++)
   {
@@ -1114,9 +1431,9 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
     {
       Float_t norm = 0;
-      for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
+      for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
         norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t);
-      for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
+      for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
       {
         if (norm==0)
           hInverseResponseBayes->SetBinContent(t,m,0);
@@ -1125,7 +1442,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
       }
     }
 
-    for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
+    for (Int_t t = kStartBin; t<=hResponse->GetNbinsX(); t++)
     {
       Float_t value = 0;
       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
@@ -1138,18 +1455,24 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
     }
 
     // this is the last guess, fill before (!) smoothing
-    for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
+    for (Int_t t=kStartBin; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
     {
+      //as bin error put the difference to the last iteration
+      //fMultiplicityESDCorrected[correlationID]->SetBinError(t, hTemp->GetBinContent(t) - fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
-      fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
+      fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t));
 
       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
     }
 
+    /*new TCanvas;
+    fMultiplicityESDCorrected[correlationID]->DrawCopy();
+    gPad->SetLogy();*/
+
     // regularization (simple smoothing)
     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
 
-    for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
+    for (Int_t t=kStartBin+2; t<hTrueSmoothed->GetNbinsX(); t++)
     {
       Float_t average = (hTemp->GetBinContent(t-1)
                          + hTemp->GetBinContent(t)
@@ -1168,7 +1491,7 @@ void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t ful
     //for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
     //  norm = norm + hTrueSmoothed->GetBinContent(t);
 
-    for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
+    for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
     {
       Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
       Float_t diff = hPrior->GetBinContent(t) - newValue;
@@ -1358,7 +1681,7 @@ void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullP
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
+  SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
 
   // TODO should be taken from correlation map
   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
@@ -1478,7 +1801,7 @@ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t ful
   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
 
-  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
+  SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
 
   NormalizeToBinWidth((TH2*) fCurrentCorrelation);
 
@@ -1544,7 +1867,7 @@ void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t ful
 }
 
 //____________________________________________________________________
-TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized)
+TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
 {
   // runs the distribution given in inputMC through the response matrix identified by
   // correlationMap and produces a measured distribution
@@ -1568,7 +1891,8 @@ TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t co
     }
   }
 
-  TH1* corr = hist->Project3D("zy");
+  TH2* corr = (TH2*) hist->Project3D("zy");
+  //corr->Rebin2D(2, 1);
   corr->Sumw2();
 
   // normalize correction for given nPart
@@ -1594,26 +1918,18 @@ TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t co
     Float_t measured = 0;
     Float_t error = 0;
 
-    for (Int_t gen=1; gen<=target->GetNbinsY(); ++gen)
+    for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
     {
-      Int_t mcGenBin = inputMC->GetXaxis()->FindBin(target->GetYaxis()->GetBinCenter(gen));
-
-      Float_t factor = 1;
-      if (normalized)
-        factor = inputMC->GetBinWidth(mcGenBin);
+      Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
 
-      measured += inputMC->GetBinContent(mcGenBin) * factor * corr->GetBinContent(gen, meas);
-      error += inputMC->GetBinError(mcGenBin) * factor * corr->GetBinContent(gen, meas);
+      measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
+      error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
     }
 
-    // bin width of the <measured> axis has to be taken into account
-    //measured /= target->GetYaxis()->GetBinWidth(meas);
-    //error /= target->GetYaxis()->GetBinWidth(meas);
-
     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
 
-    target->SetBinContent(target->GetNbinsX() / 2, meas, measured);
-    target->SetBinError(target->GetNbinsX() / 2, meas, error);
+    target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
+    target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
   }
 
   return target;
index 45c94f5..a04f580 100644 (file)
@@ -18,13 +18,13 @@ class TH3F;
 class TF1;
 class TCollection;
 
-#include <TMatrixF.h>
-#include <TVectorF.h>
+#include <TMatrixD.h>
+#include <TVectorD.h>
 
 class AliMultiplicityCorrection : public TNamed {
   public:
     enum EventType { kTrVtx = 0, kMB, kINEL };
-    enum RegularizationType { kNone = 0, kPol0, kPol1, kCurvature, kEntropy, kTest };
+    enum RegularizationType { kNone = 0, kPol0, kPol1, kEntropy, kCurvature, kTest };
 
     AliMultiplicityCorrection();
     AliMultiplicityCorrection(const Char_t* name, const Char_t* title);
@@ -40,14 +40,14 @@ class AliMultiplicityCorrection : public TNamed {
     Bool_t LoadHistograms(const Char_t* dir);
     void SaveHistograms();
     void DrawHistograms();
-    void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist);
+    void DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple = kFALSE);
 
-    void ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0);
-    void SetRegularizationParameters(RegularizationType type, Float_t weight) { fRegularizationType = type; fRegularizationWeight = weight; };
+    Int_t ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check = kFALSE, TH1* inputDist = 0);
+    void SetRegularizationParameters(RegularizationType type, Float_t weight);
 
     void ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace);
 
-    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.07, Int_t nIterations = 30, TH1* inputDist = 0);
+    void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar = 0.1, Int_t nIterations = 15, TH1* inputDist = 0);
 
     void ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace);
 
@@ -69,7 +69,7 @@ class AliMultiplicityCorrection : public TNamed {
     void SetMultiplicityESDCorrected(Int_t i, TH1F* hist) { fMultiplicityESDCorrected[i] = hist; }
 
     void SetGenMeasFromFunc(TF1* inputMC, Int_t id);
-    TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized = kFALSE);
+    TH2F* CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap);
 
     static void NormalizeToBinWidth(TH1* hist);
     static void NormalizeToBinWidth(TH2* hist);
@@ -79,18 +79,19 @@ class AliMultiplicityCorrection : public TNamed {
   protected:
     enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8 };
 
-    static const Int_t fgMaxParams; // number of fit params
+    static const Int_t fgMaxParams;  // bins in unfolded histogram = number of fit params
+    static const Int_t fgMaxInput;   // bins in measured histogram
 
-    static Double_t RegularizationPol0(Double_t *params);
-    static Double_t RegularizationPol1(Double_t *params);
-    static Double_t RegularizationTotalCurvature(Double_t *params);
-    static Double_t RegularizationEntropy(Double_t *params);
-    static Double_t RegularizationTest(Double_t *params);
+    static Double_t RegularizationPol0(TVectorD& params);
+    static Double_t RegularizationPol1(TVectorD& params);
+    static Double_t RegularizationTotalCurvature(TVectorD& params);
+    static Double_t RegularizationEntropy(TVectorD& params);
+    static Double_t RegularizationTest(TVectorD& params);
 
     static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t);
     static void MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3);
 
-    void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType);
+    void SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin);
 
     Float_t BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u);
 
@@ -98,9 +99,10 @@ class AliMultiplicityCorrection : public TNamed {
     static TH1* fCurrentCorrelation; // static variable to be accessed by MINUIT
     static TH1* fCurrentEfficiency;  // static variable to be accessed by MINUIT
 
-    static TMatrixF* fCorrelationMatrix;            // contains fCurrentCorrelation in matrix form
-    static TMatrixF* fCorrelationCovarianceMatrix;  // contains the errors of fCurrentESD
-    static TVectorF* fCurrentESDVector;             // contains fCurrentESD
+    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* fNBD;   // negative binomial distribution
 
index dbecc09..fdc6d08 100644 (file)
@@ -11,6 +11,8 @@
 #include <TH2F.h>
 #include <TH3F.h>
 #include <TParticle.h>
+#include <TRandom.h>
+#include <TNtuple.h>
 
 #include <AliLog.h>
 #include <AliESD.h>
 
 #include "esdTrackCuts/AliESDtrackCuts.h"
 #include "AliPWG0Helper.h"
+#include "AliPWG0depHelper.h"
 #include "dNdEta/AliMultiplicityCorrection.h"
+#include "AliCorrection.h"
+#include "AliCorrectionMatrix3D.h"
 
 //#define TPCMEASUREMENT
 #define ITSMEASUREMENT
@@ -31,11 +36,17 @@ ClassImp(AliMultiplicityMCSelector)
 AliMultiplicityMCSelector::AliMultiplicityMCSelector() :
   AliSelectorRL(),
   fMultiplicity(0),
-  fEsdTrackCuts(0)
+  fEsdTrackCuts(0),
+  fSystSkipParticles(kFALSE),
+  fSelectProcessType(0),
+  fParticleSpecies(0)
 {
   //
   // Constructor. Initialization of pointers
   //
+
+  for (Int_t i = 0; i<4; i++)
+    fParticleCorrection[i] = 0;
 }
 
 AliMultiplicityMCSelector::~AliMultiplicityMCSelector()
@@ -82,6 +93,32 @@ void AliMultiplicityMCSelector::SlaveBegin(TTree* tree)
   ReadUserObjects(tree);
 
   fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TString option(GetOption());
+  if (option.Contains("skip-particles"))
+  {
+    fSystSkipParticles = kTRUE;
+    AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
+  }
+
+  if (option.Contains("particle-efficiency"))
+    for (Int_t i = 0; i<4; i++)
+      fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
+
+  if (option.Contains("only-process-type-nd"))
+    fSelectProcessType = 1;
+
+  if (option.Contains("only-process-type-sd"))
+    fSelectProcessType = 2;
+
+  if (option.Contains("only-process-type-dd"))
+    fSelectProcessType = 3;
+
+  if (fSelectProcessType != 0)
+    AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
+
+  if (option.Contains("particle-species"))
+    fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec");
 }
 
 void AliMultiplicityMCSelector::Init(TTree* tree)
@@ -157,7 +194,34 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
     return kFALSE;
   }
 
+  if (fSelectProcessType > 0)
+  {
+    // getting process information; NB: this only works for Pythia
+    Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
+
+    Bool_t processEvent = kFALSE;
+
+    // non diffractive
+    if (fSelectProcessType == 1 && processtype !=92 && processtype !=93 && processtype != 94)
+      processEvent = kTRUE;
+
+    // single diffractive
+    if (fSelectProcessType == 2 && (processtype == 92 || processtype == 93))
+      processEvent = kTRUE;
+
+    // double diffractive
+    if (fSelectProcessType == 3 && processtype == 94)
+      processEvent = kTRUE;
+
+    if (!processEvent)
+    {
+      AliDebug(AliLog::kDebug, Form("Skipping this event, because it is not of the requested process type (%d)", processtype));
+      return kTRUE;
+    }
+  }
+
   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
+  eventTriggered = kTRUE; // no generated trigger for the simulated events
   Bool_t eventVertex = AliPWG0Helper::IsVertexReconstructed(fESD);
 
   // get the ESD vertex
@@ -191,6 +255,11 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
   Int_t nMCTracks20 = 0;
   Int_t nMCTracksAll = 0;
 
+  // tracks per particle species, in |eta| < 2 (systematic study)
+  Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
+  for (Int_t i = 0; i<4; ++i)
+    nMCTracksSpecies[i] = 0;
+
   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
   {
     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
@@ -203,8 +272,16 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
       continue;
     }
 
-    if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+    Bool_t debug = kFALSE;
+    if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
+    {
+      //printf("%d) DROPPED ", iMC);
+      //particle->Print();
       continue;
+    }
+
+    //printf("%d) OK      ", iMC);
+    //particle->Print();
 
     //if (particle->Pt() < kPtCut)
     //  continue;
@@ -222,9 +299,21 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
       nMCTracks20++;
 
     nMCTracksAll++;
+
+    Int_t id = -1;
+    switch (TMath::Abs(particle->GetPdgCode()))
+    {
+      case 211: id = 0; break;
+      case 321: id = 1; break;
+      case 2212: id = 2; break;
+      default: id = 3; break;
+    }
+    if (TMath::Abs(particle->Eta()) < 2.0)
+      nMCTracksSpecies[id]++;
+    if (fParticleCorrection[id])
+      fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
   }// end of mc particle
 
-  // FAKE: put back vtxMC[2]
   fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll);
 
   if (!eventTriggered || !eventVertex)
@@ -235,6 +324,11 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
   Int_t nESDTracks15 = 0;
   Int_t nESDTracks20 = 0;
 
+  // tracks per particle species, in |eta| < 2 (systematic study)
+  Int_t nESDTracksSpecies[4]; // (pi, K, p, other)
+  for (Int_t i = 0; i<4; ++i)
+    nESDTracksSpecies[i] = 0;
+
 #ifdef ITSMEASUREMENT
   // get multiplicity from ITS tracklets
   for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
@@ -245,6 +339,10 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
     if (mult->GetDeltaPhi(i) < -1000)
       continue;
 
+    // systematic study: 10% lower efficiency
+    if (fSystSkipParticles && gRandom->Uniform() < 0.1)
+      continue;
+
     Float_t theta = mult->GetTheta(i);
     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
 
@@ -259,6 +357,29 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
 
     if (TMath::Abs(eta) < 2.0)
       nESDTracks20++;
+
+    // TODO define needed, because this only works with new AliRoot
+    Int_t label = mult->GetLabel(i);
+    if (label < 0)
+      continue;
+
+    TParticle* mother = AliPWG0depHelper::FindPrimaryMother(stack, label);
+
+    if (!mother)
+      continue;
+
+    Int_t id = -1;
+    switch (TMath::Abs(mother->GetPdgCode()))
+    {
+      case 211: id = 0; break;
+      case 321: id = 1; break;
+      case 2212: id = 2; break;
+      default: id = 3; break;
+    }
+    if (TMath::Abs(eta) < 2.0)
+      nESDTracksSpecies[id]++;
+    if (fParticleCorrection[id])
+      fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
   }
 #endif
 
@@ -298,14 +419,42 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
 
     if (TMath::Abs(eta) < 2.0)
       nESDTracks20++;
+
+    Int_t label = esdTrack->GetLabel();
+    if (label < 0)
+      continue;
+
+    TParticle* mother = AliPWG0depHelper::FindPrimaryMother(stack, label);
+
+    if (!mother)
+      continue;
+
+    Int_t id = -1;
+    switch (TMath::Abs(mother->GetPdgCode()))
+    {
+      case 211: id = 0; break;
+      case 321: id = 1; break;
+      case 2212: id = 2; break;
+      default: id = 3; break;
+    }
+    if (TMath::Abs(eta) < 2.0)
+      nESDTracksSpecies[id]++;
+    if (fParticleCorrection[id])
+      fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], mother->Eta(), mother->Pt());
   }
 #endif
 
+  if (nMCTracks20 == 0 && nESDTracks20 > 3)
+    printf("WARNING: Event %lld has %d generated and %d reconstructed...\n", entry, nMCTracks05, nESDTracks05);
+
   fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
 
   // TODO should this be vtx[2] or vtxMC[2] ?
   fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
 
+  if (fParticleSpecies)
+    fParticleSpecies->Fill(nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3]);
+
   return kTRUE;
 }
 
@@ -325,6 +474,10 @@ void AliMultiplicityMCSelector::SlaveTerminate()
   }
 
   fOutput->Add(fMultiplicity);
+  for (Int_t i = 0; i < 4; ++i)
+    fOutput->Add(fParticleCorrection[i]);
+
+  fOutput->Add(fParticleSpecies);
 }
 
 void AliMultiplicityMCSelector::Terminate()
@@ -336,6 +489,9 @@ void AliMultiplicityMCSelector::Terminate()
   AliSelector::Terminate();
 
   fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
+  for (Int_t i = 0; i < 4; ++i)
+    fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
+  fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
 
   if (!fMultiplicity)
   {
@@ -346,6 +502,10 @@ void AliMultiplicityMCSelector::Terminate()
   TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
 
   fMultiplicity->SaveHistograms();
+  for (Int_t i = 0; i < 4; ++i)
+    if (fParticleCorrection[i])
+      fParticleCorrection[i]->SaveHistograms();
+  fParticleSpecies->Write();
 
   file->Close();
 
index 3dfec89..18496b5 100644 (file)
@@ -7,6 +7,8 @@
 
 class AliESDtrackCuts;
 class AliMultiplicityCorrection;
+class AliCorrection;
+class TNtuple;
 
 class AliMultiplicityMCSelector : public AliSelectorRL {
   public:
@@ -24,7 +26,13 @@ class AliMultiplicityMCSelector : public AliSelectorRL {
     void ReadUserObjects(TTree* tree);
 
     AliMultiplicityCorrection* fMultiplicity; // object containing the extracted data
-    AliESDtrackCuts* fEsdTrackCuts;     // Object containing the parameters of the esd track cuts
+    AliESDtrackCuts* fEsdTrackCuts;           // Object containing the parameters of the esd track cuts
+
+    Bool_t fSystSkipParticles;     // if true skips particles (systematic study)
+    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: (pi, k, p, rest (in |eta| < 2)) X (true, recon); (for systematic study)
 
  private:
     AliMultiplicityMCSelector(const AliMultiplicityMCSelector&);
index 4499a02..4e73859 100644 (file)
@@ -224,7 +224,7 @@ Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
   
   // getting process information NB: this only works for Pythia !!!
   Int_t processtype = AliPWG0depHelper::GetPythiaEventProcessType(header);
-    
+
   // can only read pythia headers, either directly or from cocktalil header
   AliGenEventHeader* genHeader = (AliGenEventHeader*)(header->GenEventHeader());
   
@@ -244,7 +244,7 @@ Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
   
   // non diffractive
-  if (processtype!=92 && processtype!=93 && processtype!=94) { 
+  if (processtype!=92 && processtype!=93 && processtype!=94) {
     // NB: passing the wrong process type here (1), since the process type is defined by the index in the array (here non-diffractive)
     if (triggerBiasStudy) fdNdEtaCorrectionTriggerBias[0]->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);      
     if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[0] ->FillEvent(vtxMC[2], nGoodTracks, eventTriggered, vertexReconstructed, 1);
@@ -391,30 +391,8 @@ Bool_t AlidNdEtaSystematicsSelector::Process(Long64_t entry)
       if (vertexRecoStudy)  fdNdEtaCorrectionVertexReco[2] ->FillTrackedParticle(vtxMC[2], eta, pt);
     }
 
-
-    TParticle* mother = particle;
     // find primary particle that created this particle
-    while (label >= nPrim)
-    {
-      //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
-
-      if (mother->GetMother(0) == -1)
-      {
-        AliDebug(AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
-        mother = 0;
-        break;
-      }
-
-      label = mother->GetMother(0);
-
-      mother = stack->Particle(label);
-      if (!mother)
-      {
-        AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (find mother loop).", label));
-        break;
-      }
-    }
-
+    TParticle* mother = AliPWG0depHelper::FindPrimaryMother(stack, label);
     if (!mother)
       continue;
 
index 9dd5e60..28364d6 100644 (file)
@@ -10,7 +10,7 @@
 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)
-    connectProof("jgrosseo@lxb6046");
+    connectProof("proof40@lxb6046");
 
   TString libraries("libEG;libGeom;libESD;libPWG0base");
   TString packages("PWG0base");
@@ -48,6 +48,7 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_
   if (aDebug != kFALSE)
     selectorName += "g";
 
+  //Int_t result = chain->Process(selectorName, option);
   Int_t result = executeQuery(chain, &inputList, selectorName, option);
 
   if (result != 0)
@@ -67,6 +68,17 @@ void draw(const char* fileName = "multiplicityMC.root")
   mult->LoadHistograms("Multiplicity");
 
   mult->DrawHistograms();
+
+  TH2* hist = (TH2*) gROOT->FindObject("fCorrelation3_zy");
+  canvas = new TCanvas("c1", "c1", 600, 500);
+  hist->SetStats(kFALSE);
+  hist->Draw("COLZ");
+  hist->SetTitle(";true multiplicity in |#eta| < 2;measured multiplicity in |#eta| < 2");
+  hist->GetYaxis()->SetTitleOffset(1.1);
+  gPad->SetRightMargin(0.15);
+  gPad->SetLogz();
+
+  canvas->SaveAs("Plot_Correlation.pdf");
 }
 
 void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
@@ -83,17 +95,16 @@ void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t ev
   //return;
 
 
-  mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  /*mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
   mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
-
-  return;
+  return;*/
 
   TStopwatch timer;
 
   timer.Start();
 
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1);
-  mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.2);
+  mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
   mult->DrawComparison("MinuitChi2", hist, kFALSE, kFALSE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
 
   timer.Stop();
@@ -110,7 +121,7 @@ void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t ev
   return mult;
 }
 
-void* fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Int_t histID = 2)
+void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
 {
   gSystem->Load("libPWG0base");
 
@@ -121,33 +132,84 @@ void* fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fi
 
   TFile::Open(fileNameESD);
   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
-  TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
-  //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityMB%d", histID));
+  TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID)));
+  //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID));
 
   mult->SetMultiplicityESD(histID, hist);
 
-  mult->SetMultiplicityVtx(histID, hist2);
-  mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
-  return;
-
-  mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1.1);
-  mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
-  mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
+  // small hack to get around charge conservation for full phase space ;-)
+  if (fullPhaseSpace)
+  {
+    TH1* corr = mult->GetCorrelation(histID + 4);
 
-  return;
+    for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
+      for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
+      {
+        corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
+        corr->SetBinError(i, j, corr->GetBinError(i-1, j));
+      }
+  }
 
-  //mult->ApplyGaussianMethod(histID, kFALSE);
+  /*mult->SetMultiplicityVtx(histID, hist2);
+  mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  return;*/
 
-  for (Float_t f=0.1; f<=0.11; f+=0.05)
+  if (chi2)
+  {
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
+    mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); //, hist2->ProjectionY("mymchist"));
+    mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
+  }
+  else
   {
-    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, f);
-    mult->DrawComparison(Form("Bayesian_%f", f), histID, kFALSE, kTRUE, hist2->ProjectionY());
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+    mult->DrawComparison("Bayesian", histID, kFALSE, kTRUE, hist2->ProjectionY("mymchist2"));
   }
 
   //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7);
   //mult->ApplyMinuitFit(histID, kFALSE);
   //mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
 
+}
+
+void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
+{
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  TFile::Open(fileNameESD);
+  TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
+  TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID)));
+  //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID));
+
+  mult->SetMultiplicityESD(histID, hist);
+
+  // small hack to get around charge conservation for full phase space ;-)
+  if (fullPhaseSpace)
+  {
+    TH1* corr = mult->GetCorrelation(histID + 4);
+
+    for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
+      for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
+      {
+        corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
+        corr->SetBinError(i, j, corr->GetBinError(i-1, j));
+      }
+  }
+
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+  mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
+  mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
+
+  TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult");
+
+  mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000);
+  mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result);
+  mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
 
   return mult;
 }
@@ -177,7 +239,7 @@ const char* GetEventTypeName(Int_t type)
   return 0;
 }
 
-void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
 {
   gSystem->mkdir(targetDir);
 
@@ -219,7 +281,7 @@ void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", cons
     legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
     legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
 
-    for (Float_t weight = 0; weight < 0.301; weight += 0.02)
+    for (Float_t weight = 0; weight < 1.01; weight += 0.1)
     {
       Float_t chi2MC = 0;
       Float_t residuals = 0;
@@ -339,7 +401,201 @@ void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.r
   canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
 }
 
-void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
+{
+  gSystem->mkdir(targetDir);
+
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
+  TFile::Open(fileNameESD);
+  multESD->LoadHistograms("Multiplicity");
+  mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+
+  Int_t count = 0; // just to order the saved images...
+
+  TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE");
+
+  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+  {
+    TString tmp;
+    tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
+
+    TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
+
+    for (Int_t i = 1; i <= 7; i++)
+    {
+      Int_t iter = i * 20;
+      TGraph* fitResultsMC = new TGraph;
+      fitResultsMC->SetTitle(Form("%d iterations", iter));
+      fitResultsMC->GetXaxis()->SetTitle("smoothing parameter");
+      fitResultsMC->GetYaxis()->SetTitle("P_{1} (2 <= t <= 150)");
+      fitResultsMC->SetName(Form("%s_MC_%d", tmp.Data(), i));
+      TGraph* fitResultsRes = new TGraph;
+      fitResultsRes->SetTitle(Form("%d iterations", iter));
+      fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i));
+      fitResultsRes->GetXaxis()->SetTitle("smoothing parameter");
+      fitResultsRes->GetYaxis()->SetTitle("P_{2}");
+
+      fitResultsMC->SetFillColor(0);
+      fitResultsRes->SetFillColor(0);
+      fitResultsMC->SetMarkerStyle(19+i);
+      fitResultsRes->SetMarkerStyle(19+i);
+      fitResultsRes->SetMarkerColor(kRed);
+      fitResultsRes->SetLineColor(kRed);
+
+      for (Float_t weight = 0.0; weight < 1.01; weight += 0.2)
+      {
+        Float_t chi2MC = 0;
+        Float_t residuals = 0;
+
+        mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter);
+        mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
+        mult->GetComparisonResults(&chi2MC, 0, &residuals);
+
+        fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
+        fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
+      }
+
+      fitResultsMC->Write();
+      fitResultsRes->Write();
+    }
+  }
+}
+
+void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE)
+{
+  gSystem->Load("libPWG0base");
+
+  TString name;
+  TFile* graphFile = 0;
+  if (type == -1)
+  {
+    name = "EvaluateChi2Method";
+    graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
+  }
+  else
+  {
+    name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
+    graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
+  }
+
+  TCanvas* canvas = new TCanvas(name, name, 800, 500);
+  if (type == -1)
+    canvas->SetLogx();
+  canvas->SetLogy();
+
+  TLegend* legend = new TLegend(0.6, 0.75, 0.88, 0.88);
+  legend->SetFillColor(0);
+
+  Int_t count = 1;
+
+  Float_t xMin = 1e20;
+  Float_t xMax = 0;
+
+  Float_t yMin = 1e20;
+  Float_t yMax = 0;
+
+  TString xaxis, yaxis;
+
+  while (1)
+  {
+    TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
+    TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
+
+    if (!mc || !res)
+      break;
+
+    xaxis = mc->GetXaxis()->GetTitle();
+    yaxis = mc->GetYaxis()->GetTitle();
+
+    mc->Print();
+    res->Print();
+
+    count++;
+
+    if (plotRes)
+    {
+      xMin = TMath::Min(TMath::Min(xMin, mc->GetXaxis()->GetXmin()), res->GetXaxis()->GetXmin());
+      yMin = TMath::Min(TMath::Min(yMin, mc->GetYaxis()->GetXmin()), res->GetYaxis()->GetXmin());
+
+      xMax = TMath::Max(TMath::Max(xMax, mc->GetXaxis()->GetXmax()), res->GetXaxis()->GetXmax());
+      yMax = TMath::Max(TMath::Max(yMax, mc->GetYaxis()->GetXmax()), res->GetYaxis()->GetXmax());
+    }
+    else
+    {
+      xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
+      yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
+
+      xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
+      yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
+    }
+  }
+
+  if (type >= 0)
+  {
+    xaxis = "smoothing parameter";
+  }
+  else if (type == -1)
+  {
+    xaxis = "weight parameter";
+    //yaxis = "P_{3} (2 <= t <= 150)";
+  }
+  yaxis = "P_{1} (2 <= t <= 150)";
+
+  printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
+
+  TGraph* dummy = new TGraph;
+  dummy->SetPoint(0, xMin, yMin);
+  dummy->SetPoint(1, xMax, yMax);
+  dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
+
+  dummy->SetMarkerColor(0);
+  dummy->Draw("AP");
+
+  count = 1;
+
+  while (1)
+  {
+    TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
+    TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
+
+    printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
+
+    if (!mc || !res)
+      break;
+
+    printf("Loaded %d sucessful.\n", count);
+
+    if (type == -1)
+    {
+      legend->AddEntry(mc, Form("Eq. (%d)", count+9));
+    }
+    else
+      legend->AddEntry(mc);
+
+    mc->Draw("SAME PC");
+
+    if (plotRes)
+    {
+      legend->AddEntry(res);
+      res->Draw("SAME PC");
+    }
+
+    count++;
+  }
+
+  legend->Draw();
+
+  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
+  canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
+}
+
+void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
 {
   gSystem->mkdir(targetDir);
 
@@ -356,25 +612,24 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch
 
   mult->SetMultiplicityESD(histID, hist);
 
-  TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600);
-  TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
-  legend->SetFillColor(0);
-
-  Float_t min = 1e10;
-  Float_t max = 0;
-
-  TGraph* first = 0;
   Int_t count = 0; // just to order the saved images...
 
-  Bool_t firstLoop = kTRUE;
+  TGraph* fitResultsMC = 0;
+  TGraph* fitResultsRes = 0;
+
+  TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
 
   for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
-  //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
+//  for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
   {
-    TGraph* fitResultsMC = new TGraph;
-    fitResultsMC->SetTitle(";Weight Parameter");
-    TGraph* fitResultsRes = new TGraph;
-    fitResultsRes->SetTitle(";Weight Parameter");
+    fitResultsMC = new TGraph;
+    fitResultsMC->SetTitle(Form("%s", GetRegName(type)));
+    fitResultsMC->GetXaxis()->SetTitle("Weight Parameter");
+    fitResultsMC->SetName(Form("EvaluateChi2Method_MC_%d", type));
+    fitResultsRes = new TGraph;
+    fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type)));
+    fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type));
+    fitResultsRes->GetXaxis()->SetTitle("Weight Parameter");
 
     fitResultsMC->SetFillColor(0);
     fitResultsRes->SetFillColor(0);
@@ -383,51 +638,44 @@ void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const ch
     fitResultsRes->SetMarkerColor(kRed);
     fitResultsRes->SetLineColor(kRed);
 
-    legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetRegName(type)));
-    legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetRegName(type)));
+    for (Int_t i=0; i<5; ++i)
+    {
+      Float_t weight = TMath::Power(TMath::Sqrt(10), i+6);
+      //Float_t weight = TMath::Power(10, i+2);
 
-    if (first == 0)
-      first = fitResultsMC;
+      //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5;
 
-    for (Float_t weight = 1e-4; weight < 2e4; weight *= 10)
-    //for (Float_t weight = 0.1; weight < 10; weight *= TMath::Sqrt(TMath::Sqrt(10)))
-    {
       Float_t chi2MC = 0;
       Float_t residuals = 0;
+      Float_t chi2Limit = 0;
+
+      TString runName;
+      runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight);
 
       mult->SetRegularizationParameters(type, weight);
-      mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
-      mult->DrawComparison(Form("%s/MinuitChi2_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
+      Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
+      mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY());
+      if (status != 0)
+      {
+        printf("MINUIT did not succeed. Skipping...\n");
+        continue;
+      }
+
       mult->GetComparisonResults(&chi2MC, 0, &residuals);
+      TH1* result = mult->GetMultiplicityESDCorrected(histID);
+      result->SetName(runName);
+      result->Write();
 
       fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
       fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
-
-      min = TMath::Min(min, TMath::Min(chi2MC, residuals));
-      max = TMath::Max(max, TMath::Max(chi2MC, residuals));
     }
 
-    fitResultsMC->Print();
-    fitResultsRes->Print();
-
-    canvas->cd();
-    fitResultsMC->Draw(Form("%s CP", (firstLoop) ? "A" : "SAME"));
-    fitResultsRes->Draw("SAME CP");
-
-    firstLoop = kFALSE;
+    graphFile->cd();
+    fitResultsMC->Write();
+    fitResultsRes->Write();
   }
 
-  gPad->SetLogx();
-  gPad->SetLogy();
-  printf("min = %f, max = %f\n", min, max);
-  if (min <= 0)
-    min = 1e-5;
-  first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
-
-  legend->Draw();
-
-  canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
-  canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
+  graphFile->Close();
 }
 
 void EvaluateChi2MethodAll()
@@ -448,7 +696,7 @@ void EvaluateBayesianMethodAll()
   EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
 }
 
-void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
+void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
 {
   gSystem->mkdir(targetDir);
 
@@ -465,16 +713,17 @@ void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char*
 
   mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
 
-  TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 800);
-  canvas->Divide(3, 2);
+  TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200);
+  canvas->Divide(3, 3);
 
   Int_t count = 0;
 
-  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
+  for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
   {
-    TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY();
+    TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc");
+    mc->Sumw2();
 
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
     mult->ApplyMinuitFit(histID, kFALSE, type);
     mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc);
     TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
@@ -486,12 +735,12 @@ void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char*
     mc->GetXaxis()->SetRangeUser(0, 150);
     chi2Result->GetXaxis()->SetRangeUser(0, 150);
 
-    // skip errors for now
+/*    // skip errors for now
     for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
     {
       chi2Result->SetBinError(i, 0);
       bayesResult->SetBinError(i, 0);
-    }
+    }*/
 
     canvas->cd(++count);
     mc->SetFillColor(kYellow);
@@ -503,27 +752,27 @@ void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char*
     gPad->SetLogy();
 
     canvas->cd(count + 3);
-    chi2Result->Divide(chi2Result, mc, 1, 1, "B");
-    bayesResult->Divide(bayesResult, mc, 1, 1, "B");
+    chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio");
+    bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio");
+    chi2ResultRatio->Divide(chi2Result, mc, 1, 1, "");
+    bayesResultRatio->Divide(bayesResult, mc, 1, 1, "");
 
-    // skip errors for now
-    for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
-    {
-      chi2Result->SetBinError(i, 0);
-      bayesResult->SetBinError(i, 0);
-    }
+    chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5);
 
-    chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+    chi2ResultRatio->DrawCopy("HIST");
+    bayesResultRatio->DrawCopy("SAME HIST");
 
-    chi2Result->DrawCopy("");
-    bayesResult->DrawCopy("SAME");
+    canvas->cd(count + 6);
+    chi2Result->Divide(chi2Result, bayesResult, 1, 1, "");
+    chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+    chi2Result->DrawCopy("HIST");
   }
 
   canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
   canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
 }
 
-void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
+void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
 {
   gSystem->Load("libPWG0base");
 
@@ -532,16 +781,17 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID
   TFile::Open(fileNameMC);
   mult->LoadHistograms("Multiplicity");
 
-  const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
+  const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" };
 
-  TGraph* fitResultsChi2 = new TGraph;
-  fitResultsChi2->SetTitle(";Nevents;Chi2");
-  TGraph* fitResultsBayes = new TGraph;
-  fitResultsBayes->SetTitle(";Nevents;Chi2");
-  TGraph* fitResultsChi2Limit = new TGraph;
-  fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
-  TGraph* fitResultsBayesLimit = new TGraph;
-  fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
+  TGraph* fitResultsChi2 = new TGraph;       fitResultsChi2->SetTitle(";Nevents;Chi2");
+  TGraph* fitResultsBayes = new TGraph;      fitResultsBayes->SetTitle(";Nevents;Chi2");
+  TGraph* fitResultsChi2Limit = new TGraph;  fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
+  TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
+
+  fitResultsChi2->SetName("fitResultsChi2");
+  fitResultsBayes->SetName("fitResultsBayes");
+  fitResultsChi2Limit->SetName("fitResultsChi2Limit");
+  fitResultsBayesLimit->SetName("fitResultsBayesLimit");
 
   TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
   canvas->Divide(5, 2);
@@ -549,6 +799,9 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID
   Float_t min = 1e10;
   Float_t max = 0;
 
+  TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE");
+  file->Close();
+
   for (Int_t i=0; i<5; ++i)
   {
     TFile::Open(files[i]);
@@ -557,9 +810,9 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID
 
     mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
     Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries();
-    TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
+    TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
 
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
     mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
     mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
 
@@ -571,11 +824,11 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID
     min = TMath::Min(min, chi2MC);
     max = TMath::Max(max, chi2MC);
 
-    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
+    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
 
-    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1);
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.07, 10);
     mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
-    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
+    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
     mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
     fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC);
     fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
@@ -617,6 +870,12 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID
 
     chi2Result->DrawCopy("");
     bayesResult->DrawCopy("SAME");
+
+    TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
+    mc->Write();
+    chi2Result->Write();
+    bayesResult->Write();
+    file->Close();
   }
 
   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
@@ -647,9 +906,16 @@ void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID
 
   canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
   canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
+
+  TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
+  fitResultsChi2->Write();
+  fitResultsBayes->Write();
+  fitResultsChi2Limit->Write();
+  fitResultsBayesLimit->Write();
+  file->Close();
 }
 
-void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
+void StartingConditions(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
 {
   gSystem->Load("libPWG0base");
 
@@ -658,14 +924,14 @@ void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t hi
   TFile::Open(fileNameMC);
   mult->LoadHistograms("Multiplicity");
 
-  const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
+  const char* files[] = { "multiplicityMC_1M_3.root", "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root" }
 
   // this one we try to unfold
   TFile::Open(files[0]);
   AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
   multESD->LoadHistograms("Multiplicity");
   mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
-  TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
+  TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc");
 
   TGraph* fitResultsChi2 = new TGraph;
   fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
@@ -691,6 +957,10 @@ void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t hi
 
   TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
 
+  TFile* file = TFile::Open("StartingConditions.root", "RECREATE");
+  mc->Write();
+  file->Close();
+
   for (Int_t i=0; i<8; ++i)
   {
     if (i == 0)
@@ -722,7 +992,7 @@ void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t hi
         startCond->SetBinContent(j, 1);
     }
 
-    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
     mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
     mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
 
@@ -734,13 +1004,13 @@ void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t hi
     min = TMath::Min(min, chi2MC);
     max = TMath::Max(max, chi2MC);
 
-    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
+    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
     if (!firstChi)
       firstChi = (TH1*) chi2Result->Clone("firstChi");
 
-    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1, 30, startCond);
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, startCond);
     mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
-    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
+    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
     if (!firstBayesian)
       firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
 
@@ -748,6 +1018,11 @@ void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t hi
     fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
     fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
 
+    TFile* file = TFile::Open("StartingConditions.root", "UPDATE");
+    chi2Result->Write();
+    bayesResult->Write();
+    file->Close();
+
     min = TMath::Min(min, chi2MC);
     max = TMath::Max(max, chi2MC);
     mc->GetXaxis()->SetRangeUser(0, 150);
@@ -838,8 +1113,217 @@ void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t hi
   canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
 }
 
+void DifferentSamples(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
+{
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TFile::Open(fileNameMC);
+  mult->LoadHistograms("Multiplicity");
+
+  const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" };
+
+  TGraph* fitResultsChi2 = new TGraph;
+  fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
+  TGraph* fitResultsBayes = new TGraph;
+  fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
+  TGraph* fitResultsChi2Limit = new TGraph;
+  fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
+  TGraph* fitResultsBayesLimit = new TGraph;
+  fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
+
+  TCanvas* canvasA = new TCanvas("DifferentSamplesA", "DifferentSamplesA", 1200, 600);
+  canvasA->Divide(4, 2);
+
+  TCanvas* canvasB = new TCanvas("DifferentSamplesB", "DifferentSamplesB", 1200, 600);
+  canvasB->Divide(4, 2);
+
+  TCanvas* canvas4 = new TCanvas("DifferentSamples4", "DifferentSamples4", 1000, 400);
+  canvas4->Divide(2, 1);
+
+  TCanvas* canvas3 = new TCanvas("DifferentSamples3", "DifferentSamples3", 1000, 400);
+  canvas3->Divide(2, 1);
+
+  Float_t min = 1e10;
+  Float_t max = 0;
+
+  TH1* firstChi = 0;
+  TH1* firstBayesian = 0;
+
+  TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
+
+  TFile* file = TFile::Open("DifferentSamples.root", "RECREATE");
+  file->Close();
+
+  for (Int_t i=0; i<8; ++i)
+  {
+    TFile::Open(files[i]);
+    AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
+    multESD->LoadHistograms("Multiplicity");
+    mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
+    TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
+    mc->Sumw2();
+
+    mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
+    mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE);
+    mult->DrawComparison(Form("DifferentSamples_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
+
+    Float_t chi2MC = 0;
+    Int_t chi2MCLimit = 0;
+    mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+    fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
+    fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
+    min = TMath::Min(min, chi2MC);
+    max = TMath::Max(max, chi2MC);
+
+    TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
+    if (!firstChi)
+      firstChi = (TH1*) chi2Result->Clone("firstChi");
+
+    mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100);
+    mult->DrawComparison(Form("DifferentSamples_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
+    TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
+    if (!firstBayesian)
+      firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
+
+    TFile* file = TFile::Open("DifferentSamples.root", "UPDATE");
+    mc->Write();
+    chi2Result->Write();
+    bayesResult->Write();
+    file->Close();
+
+    mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
+    fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
+    fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
+
+    min = TMath::Min(min, chi2MC);
+    max = TMath::Max(max, chi2MC);
+    mc->GetXaxis()->SetRangeUser(0, 150);
+    chi2Result->GetXaxis()->SetRangeUser(0, 150);
+
+    // skip errors for now
+    for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+    {
+      chi2Result->SetBinError(j, 0);
+      bayesResult->SetBinError(j, 0);
+    }
+
+    canvas4->cd(1);
+    TH1* tmp = (TH1*) chi2Result->Clone("tmp");
+    tmp->SetTitle("Unfolded/MC;Npart;Ratio");
+    tmp->Divide(mc);
+    tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+    tmp->GetXaxis()->SetRangeUser(0, 200);
+    tmp->SetLineColor(i+1);
+    tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+    canvas4->cd(2);
+    tmp = (TH1*) bayesResult->Clone("tmp");
+    tmp->SetTitle("Unfolded/MC;Npart;Ratio");
+    tmp->Divide(mc);
+    tmp->SetLineColor(i+1);
+    tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+    tmp->GetXaxis()->SetRangeUser(0, 200);
+    tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+    canvas3->cd(1);
+    TH1* tmp = (TH1*) chi2Result->Clone("tmp");
+    tmp->SetTitle("Ratio to first result;Npart;Ratio");
+    tmp->Divide(firstChi);
+    tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+    tmp->GetXaxis()->SetRangeUser(0, 200);
+    tmp->SetLineColor(i+1);
+    legend->AddEntry(tmp, Form("%d", i));
+    tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+    canvas3->cd(2);
+    tmp = (TH1*) bayesResult->Clone("tmp");
+    tmp->SetTitle("Ratio to first result;Npart;Ratio");
+    tmp->Divide(firstBayesian);
+    tmp->SetLineColor(i+1);
+    tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
+    tmp->GetXaxis()->SetRangeUser(0, 200);
+    tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
+
+    if (i < 4)
+    {
+      canvasA->cd(i+1);
+    }
+    else
+      canvasB->cd(i+1-4);
+
+    mc->SetFillColor(kYellow);
+    mc->DrawCopy();
+    chi2Result->SetLineColor(kRed);
+    chi2Result->DrawCopy("SAME");
+    bayesResult->SetLineColor(kBlue);
+    bayesResult->DrawCopy("SAME");
+    gPad->SetLogy();
+
+    if (i < 4)
+    {
+      canvasA->cd(i+5);
+    }
+    else
+      canvasB->cd(i+5-4);
+
+    chi2Result->Divide(chi2Result, mc, 1, 1, "B");
+    bayesResult->Divide(bayesResult, mc, 1, 1, "B");
+
+    // skip errors for now
+    for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
+    {
+      chi2Result->SetBinError(j, 0);
+      bayesResult->SetBinError(j, 0);
+    }
+
+    chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
+    chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
+
+    chi2Result->DrawCopy("");
+    bayesResult->DrawCopy("SAME");
+  }
+
+  canvas3->cd(1);
+  legend->Draw();
+
+  canvasA->SaveAs(Form("%s.gif", canvasA->GetName()));
+  canvasB->SaveAs(Form("%s.gif", canvasB->GetName()));
+
+  TCanvas* canvas2 = new TCanvas("DifferentSamples2", "DifferentSamples2", 800, 400);
+  canvas2->Divide(2, 1);
+
+  canvas2->cd(1);
+  fitResultsChi2->SetMarkerStyle(20);
+  fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
+  fitResultsChi2->Draw("AP");
+
+  fitResultsBayes->SetMarkerStyle(3);
+  fitResultsBayes->SetMarkerColor(2);
+  fitResultsBayes->Draw("P SAME");
+
+  gPad->SetLogy();
+
+  canvas2->cd(2);
+  fitResultsChi2Limit->SetMarkerStyle(20);
+  fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
+  fitResultsChi2Limit->Draw("AP");
+
+  fitResultsBayesLimit->SetMarkerStyle(3);
+  fitResultsBayesLimit->SetMarkerColor(2);
+  fitResultsBayesLimit->Draw("P SAME");
+
+  canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
+  canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
+  canvas4->SaveAs(Form("%s.gif", canvas4->GetName()));
+}
+
 void Merge(Int_t n, const char** files, const char* output)
 {
+  // const char* files[] = { "multiplicityMC_100k_1.root",  "multiplicityMC_100k_2.root",  "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root",  "multiplicityMC_100k_5.root",  "multiplicityMC_100k_6.root",  "multiplicityMC_100k_7.root",  "multiplicityMC_100k_8.root" };
+
+
   gSystem->Load("libPWG0base");
 
   AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
@@ -859,7 +1343,7 @@ void Merge(Int_t n, const char** files, const char* output)
 
   data[0]->Merge(&list);
 
-  data[0]->DrawHistograms();
+  //data[0]->DrawHistograms();
 
   TFile::Open(output, "RECREATE");
   data[0]->SaveHistograms();
@@ -879,7 +1363,7 @@ void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
 
   if (caseNo >= 4)
   {
-    func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
+    func = 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);
     func->SetParNames("scaling", "averagen", "k");
   }
 
@@ -890,15 +1374,18 @@ void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
     case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
     case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
     case 4: func->SetParameters(1e7, 10, 2); break;
-    case 5: func->SetParameters(1e7, 20, 3); break;
+    case 5: func->SetParameters(1, 13, 7); break;
     case 6: func->SetParameters(1e7, 30, 4); break;
-    case 7: func->SetParameters(1e7, 70, 2); break;
+    case 7: func->SetParameters(1e7, 30, 2); break;
     case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
 
     default: return;
   }
 
-  mult->SetGenMeasFromFunc(func, 2);
+  new TCanvas;
+  func->Draw();
+
+  mult->SetGenMeasFromFunc(func, 3);
 
   TFile::Open("out.root", "RECREATE");
   mult->SaveHistograms();
@@ -906,7 +1393,7 @@ void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
   //mult->ApplyBayesianMethod(2, kFALSE);
   //mult->ApplyMinuitFit(2, kFALSE);
   //mult->ApplyGaussianMethod(2, kFALSE);
-  mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
+  //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
 }
 
 void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2)
@@ -947,7 +1434,7 @@ void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t co
   }
 
   new TCanvas;
-  proj->Draw("COLZ");
+  proj->DrawCopy("COLZ");
 
   TH1* scaling = proj->ProjectionY("scaling", 1, 1);
   scaling->Reset();
@@ -997,18 +1484,19 @@ void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t co
     printf("Fitting %d...\n", i);
 
     TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e");
+
     //hist->GetXaxis()->SetRangeUser(0, 50);
     //lognormal->SetParameter(0, hist->GetMaximum());
     hist->Fit(fitWith, "0 M", "");
 
     TF1* func = hist->GetListOfFunctions()->FindObject(fitWith);
 
-    if (((i-1) % 15 == 0) || ((i % 5 == 0) && i < 30))
+    if (0 && (i % 5 == 0))
     {
-      new TCanvas;
+      pad = new TCanvas;
       hist->Draw();
       func->Clone()->Draw("SAME");
-      gPad->SetLogy();
+      pad->SetLogy();
     }
 
     scaling->Fill(i, func->GetParameter(0));
@@ -1035,21 +1523,28 @@ void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t co
 
   //TF1* scalingFit = new TF1("mypol0", "[0]");
   TF1* scalingFit = over;
-  scaling->Fit(scalingFit, "", "", 3, 100);
+  scaling->Fit(scalingFit, "", "", 3, 140);
+  scalingFit->SetRange(0, 200);
+  scalingFit->Draw("SAME");
 
   c1->cd(2);
   mean->Draw("P");
 
   //TF1* meanFit = log;
   TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
-  mean->Fit(meanFit, "", "", 3, 100);
+  mean->Fit(meanFit, "", "", 3, 140);
+  meanFit->SetRange(0, 200);
+  meanFit->Draw("SAME");
 
   c1->cd(3);
   width->Draw("P");
 
   //TF1* widthFit = over;
-  TF1* widthFit = new TF1("mypol2", "[0]+[1]*x+[2]*x*x");
-  width->Fit(widthFit, "", "", 5, 100);
+  TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)");
+  widthFit->SetParLimits(2, 1e-5, 1e5);
+  width->Fit(widthFit, "", "", 5, 140);
+  widthFit->SetRange(0, 200);
+  widthFit->Draw("SAME");
 
   // build new correction matrix
   TH2* new = (TH2*) proj->Clone("new");
@@ -1067,7 +1562,7 @@ void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t co
 
     for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
     {
-      if (i < 21)
+      if (i < 11)
       {
         // leave bins 1..20 untouched
         new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j));
@@ -1137,3 +1632,26 @@ void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t co
   proj1->Draw();
   proj2->Draw("SAME");
 }
+
+void GetCrossSections(const char* fileName)
+{
+  gSystem->Load("libPWG0base");
+
+  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
+
+  TFile::Open(fileName);
+  mult->LoadHistograms("Multiplicity");
+
+  TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2");
+  xSection2->Sumw2();
+  xSection2->Scale(1.0 / xSection2->Integral());
+
+  TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15");
+  xSection15->Sumw2();
+  xSection15->Scale(1.0 / xSection15->Integral());
+
+  TFile::Open("crosssection.root", "RECREATE");
+  xSection2->Write();
+  xSection15->Write();
+  gFile->Close();
+}
index b2d03a2..c03055e 100644 (file)
 void testAnalysis2(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* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", const char* option = "", const char* proofServer = "jgrosseo@lxb6046")
 {
   if (aProof)
-  {
     connectProof(proofServer);
-    gProof->AddInput(new TParameter<long>("PROOF_MaxSlavesPerNode", (long)2));
-    gProof->AddInput(new TNamed("PROOF_Packetizer", "TAdaptivePacketizer"));
-  }
 
   TString libraries("libEG;libGeom;libESD;libPWG0base");
-  TString packages("adaptivepacketizer;PWG0base");
+  TString packages("PWG0base");
 
   if (!prepareQuery(libraries, packages, kTRUE))
     return;
 
-  //TODO somehow prevent loading several times
   gROOT->ProcessLine(".L CreateCuts.C");
   gROOT->ProcessLine(".L drawPlots.C");
 
index a2bfbba..b817c45 100644 (file)
@@ -12,10 +12,11 @@ HDRS = dNdEta/AlidNdEtaCorrectionSelector.h \
        esdTrackCuts/AliTestESDtrackCutsSelector.h \
        TPC/AliROCESDAnalysisSelector.h \
        TPC/AliROCRawAnalysisSelector.h \
-       TPC/AliROCClusterAnalysisSelector.h 
+       TPC/AliROCClusterAnalysisSelector.h \
+       highMultiplicity/AliHighMultiplicitySelector.h
 
 SRCS = $(HDRS:.h=.cxx)
 
 DHDR= PWG0selectorsLinkDef.h
 
-EINCLUDE=PYTHIA6 EVGEN TPC RAW
+EINCLUDE=PYTHIA6 EVGEN TPC RAW ITS