Fix the caching of the parameters in TF1NormSum to avoid recomputing integrals...

Fix the caching of the parameters in TF1NormSum to avoid recomputing integrals for same parameter values

hist/hist/src/TF1NormSum.cxx

@@ -286,7 +286,7 @@ void TF1NormSum::SetParameters(const double* params)//params should have the siz
 {
    // Initialize array of all parameters.
    // double *params must contains first an array of the coefficients, then an array of the parameters.
-   
+
    for (unsigned int n=0; n<fNOfFunctions; n++)                         //initialization of the coefficients
    {
       fCoeffs[n] = params[n];
@@ -312,18 +312,24 @@ void TF1NormSum::SetParameters(const double* params)//params should have the siz
          }
          else
          {
-            noncstparams[n][k] = params[k+fNOfFunctions+offset];
+            noncstparams[n][k] = params[k+fNOfFunctions+offset]; 
             totalparams[n][i]  = params[k+fNOfFunctions+offset];
             //std::cout << " params " << k+fNOfFunctions+offset << " = " << params[k+fNOfFunctions+offset] << std::endl;
             k++;
          }
       }
       fParams[n]    =  noncstparams[n].data();                          // fParams doesn't take the coefficients, and neither the cst parameters
-      fFunctions[n] -> SetParameters(totalparams[n].data());
+      // check if function has really changed the parameters
+      bool equalParams = true;
+      for (int l = 0; l < fFunctions[n]->GetNpar(); ++l)
+         equalParams &= (fFunctions[n]->GetParameter(l) == totalparams[n][l]);
+      
+      if (!equalParams) fFunctions[n] -> SetParameters(totalparams[n].data());
       
       fixedvalue = 1.;
       if (fFunctions[n] -> GetNumber() == 200)  fixedvalue = 0.;        // if this is an exponential, the fixed value is zero
-      fFunctions[n] -> FixParameter(fCstIndexes[n], fixedvalue);
+      if (!equalParams && fFunctions[n]->GetParameter(fCstIndexes[n]) != fixedvalue) 
+         fFunctions[n] -> FixParameter(fCstIndexes[n], fixedvalue);
       // fFunctions[n]->Print();
       //std::cout << "coeff " << n << " : " << fCoeffs[n] << std::endl;
       

tutorials/fit/fitNormSum.C

-#include <stdio.h>
+
 #include <TMath.h>
 #include <TCanvas.h>
-#include <TROOT.h>
-#include <TChain.h>
-#include <TObject.h>
-#include <TRandom.h>
 #include <TF1NormSum.h>
 #include <TF1.h>
-#include <TH1F.h>
-#include <Math/PdfFuncMathCore.h>  //for the crystalball function
+#include <TH1.h>
 
 
 using namespace std;
@@ -35,19 +30,23 @@ void fitNormSum()
    
    //***************************************************************************************************
 
-   Int_t NEvents = 1e6;
+
+   const int nsig = 5.E4;
+   const int nbkg = 1.e6;
+   Int_t NEvents = nsig+nbkg;
    Int_t NBins   = 1e3;
-   
+
+   double signal_mean = 3;
    TF1 *f_cb    = new TF1("MyCrystalBall","crystalball",-5.,5.);
    TF1 *f_exp   = new TF1("MyExponential","expo",-5.,5.);
    
    // I.:
    f_exp-> SetParameters(1.,-0.3);
-   f_cb -> SetParameters(1,3,0.3,2,1.5);
+   f_cb -> SetParameters(1,signal_mean,0.3,2,1.5);
    
    // CONSTRUCTION OF THE TF1NORMSUM OBJECT ........................................
    // 1) :
-   TF1NormSum *fnorm_exp_cb = new TF1NormSum(f_exp,f_cb,1.e6,5e4);
+   TF1NormSum *fnorm_exp_cb = new TF1NormSum(f_cb,f_exp,nsig,nbkg);
    // 4) :
    
    TF1   * f_sum = new TF1("fsum", *fnorm_exp_cb, -5., 5., fnorm_exp_cb->GetNpar());
@@ -55,27 +54,42 @@ void fitNormSum()
    
    // III.:
    f_sum->SetParameters( fnorm_exp_cb->GetParameters().data() );
-   f_sum->SetParName(0,"NBackground");
-   f_sum->SetParName(1,"NSignal");
+   f_sum->SetParName(1,"NBackground");
+   f_sum->SetParName(0,"NSignal");
    for (int i = 2; i < f_sum->GetNpar(); ++i) 
       f_sum->SetParName(i,fnorm_exp_cb->GetParName(i) );
   
-   //HISTOGRAM TO FIT ..............................................................
-   TH1F *h_sum = new TH1F("h_ExpCB", "Exponential Bkg + CrystalBall function", NBins, -5., 5.);
+   //GENERATE HISTOGRAM TO FIT ..............................................................
+   TStopwatch w;
+   w.Start(); 
+   TH1D *h_sum = new TH1D("h_ExpCB", "Exponential Bkg + CrystalBall function", NBins, -5., 5.);
    for (int i=0; i<NEvents; i++)
    {
       h_sum -> Fill(f_sum -> GetRandom());
    }
+   printf("Time to generate %d events:  ",NEvents);
+   w.Print();
+   //TH1F *h_orig = new TH1F(*h_sum);
+   
+   // need to scale histogram with width since we are fitting a density
    h_sum -> Sumw2();
-   h_sum -> Scale(1., "width");
-   h_sum->Draw();
-   TH1F *h_copy = new TH1F(*h_sum);
+   h_sum -> Scale(1., "width");   
 
-   //fit
+   //fit - use Minuit2 if available
+   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
    new TCanvas("Fit","Fit",800,1000);
-   h_copy -> Fit("fsum");
-   h_copy -> Draw();
-   
-   
+   // do a least-square fit of the spectrum
+   auto result = h_sum -> Fit("fsum","SQ");
+   result->Print();
+   h_sum -> Draw();
+   printf("Time to fit using ROOT TF1Normsum: ");
+   w.Print();
+
+   // test if parameters are fine
+   std::vector<double>  pref = {nsig, nbkg, signal_mean};
+   for (unsigned int i = 0; i< pref.size(); ++i)  {
+      if (!TMath::AreEqualAbs(pref[i], f_sum->GetParameter(i), f_sum->GetParError(i)*10.) )
+         Error("testFitNormSum","Difference found in fitted %s - difference is %g sigma",f_sum->GetParName(i), (f_sum->GetParameter(i)-pref[i])/f_sum->GetParError(i));
+   }
    
 }