added init here

[u/mrichter/AliRoot.git] / ANALYSIS / AliUnfolding.cxx

diff --git a/ANALYSIS/AliUnfolding.cxx b/ANALYSIS/AliUnfolding.cxx
index a34974c8a3a05a9daf736f1d4f0b8fe8dd279b0d..82d0acc09248bf17782b640813e3fdf6b2e6a4e3 100644 (file)
--- a/ANALYSIS/AliUnfolding.cxx
+++ b/ANALYSIS/AliUnfolding.cxx
@@ -55,7 +55,8 @@ Float_t AliUnfolding::fgRegularizationWeight = 10000;
 Int_t AliUnfolding::fgSkipBinsBegin = 0;
 Float_t AliUnfolding::fgMinuitStepSize = 0.1;                 // (usually not needed to be changed) step size in minimization
 Float_t AliUnfolding::fgMinuitPrecision = 1e-6;               // minuit precision
-Int_t   AliUnfolding::fgMinuitMaxIterations = 1e6;           // minuit maximum number of iterations
+Int_t   AliUnfolding::fgMinuitMaxIterations = 1000000;            // minuit maximum number of iterations\r
+Double_t AliUnfolding::fgMinuitStrategy = 1.;                 // minuit strategy
 Bool_t AliUnfolding::fgMinimumInitialValue = kFALSE;          // set all initial values at least to the smallest value among the initial values
 Float_t AliUnfolding::fgMinimumInitialValueFix = -1;
 Bool_t AliUnfolding::fgNormalizeInput = kFALSE;               // normalize input spectrum
@@ -196,7 +197,7 @@ void AliUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
   fgBayesianSmoothing = smoothing;
   fgBayesianIterations = nIterations;
 
-  Printf("AliUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);
+  Printf("AliUnfolding::SetBayesianParameters --> Parameters set to %d iterations with smoothing %f", fgBayesianIterations, fgBayesianSmoothing);\r
 }
 
 //____________________________________________________________________
@@ -394,6 +395,8 @@ Int_t AliUnfolding::UnfoldWithMinuit(TH2* correlation, TH1* efficiency, TH1* mea
 
   minuit->SetMaxIterations(fgMinuitMaxIterations);
 
+  minuit->ExecuteCommand("SET STRATEGY",&fgMinuitStrategy,1);
+
   for (Int_t i=0; i<fgMaxParams; i++)
     (*fgEntropyAPriori)[i] = 1;
 
@@ -552,6 +555,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
   Double_t* result = new Double_t[kMaxT];
   Double_t* efficiency = new Double_t[kMaxT];
   Double_t* binWidths = new Double_t[kMaxT];
+  Double_t* normResponse = new Double_t[kMaxT]; \r
 
   Double_t** response = new Double_t*[kMaxT];
   Double_t** inverseResponse = new Double_t*[kMaxT];
@@ -559,6 +563,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
   {
     response[i] = new Double_t[kMaxM];
     inverseResponse[i] = new Double_t[kMaxM];
+    normResponse[i] = 0;\r
   }
 
   // for normalization
@@ -572,6 +577,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
     {
       response[t][m] = correlation->GetBinContent(t+1, m+1);
       inverseResponse[t][m] = 0;
+      normResponse[t] += correlation->GetBinContent(t+1, m+1);\r
     }
   }
 
@@ -587,6 +593,13 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
     prior[t] = measuredCopy[t];
     result[t] = 0;
     binWidths[t] = aResult->GetXaxis()->GetBinWidth(t+1);
+\r
+    for (Int_t m=0; m<kMaxM; m++) { // Normalise response matrix\r
+      if (normResponse[t] != 0) \r
+       response[t][m] /= normResponse[t];\r
+      else\r
+        Printf("AliUnfolding::UnfoldWithBayesian: Empty row,column in response matrix, for truth bin %d",t);\r
+    }\r
   }
 
   // pick prior distribution
@@ -616,7 +629,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
     {
       Float_t norm = 0;
       for (Int_t t = kStartBin; t<kMaxT; t++)
-        norm += response[t][m] * prior[t];
+        norm += response[t][m] * prior[t] * efficiency[t];\r
 
       // calc. chi2: (measured - response * prior) / error
       if (measuredError[m] > 0)
@@ -628,7 +641,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
       if (norm > 0)
       {
         for (Int_t t = kStartBin; t<kMaxT; t++)
-          inverseResponse[t][m] = response[t][m] * prior[t] / norm;
+          inverseResponse[t][m] = response[t][m] * prior[t] * efficiency[t] / norm;\r
       }
       else
       {
@@ -712,6 +725,7 @@ Int_t AliUnfolding::UnfoldWithBayesian(TH2* correlation, TH1* aEfficiency, TH1*
   delete[] result;
   delete[] efficiency;
   delete[] binWidths;
+  delete[] normResponse;\r
 
   for (Int_t i=0; i<kMaxT; i++)
   {