Ajout du fichier pour interface avec le fit minuit, Reza 25/02/2019

tkfit_minuit.h

+/*  PAON4 analysis software 
+    class to encapsulate Minuit for trkfit.h .cc and trkacxfit.cc 
+    R. Ansari, Fevrier 2019                                             */
+
+#ifndef TKFIT_MINUIT_H_SEEN
+#define TKFIT_MINUIT_H_SEEN
+/Minuit
+
+#include "Math/MinimizerOptions.h"
+#include "Minuit2/FCNBase.h"
+#include "Minuit2/FunctionMinimum.h"
+#include "Minuit2/MnUserParameterState.h"
+#include "Minuit2/MnPrint.h"
+#include "Minuit2/MnMigrad.h"
+#include "Minuit2/MnMinimize.h"
+#include "Minuit2/MnMinos.h"
+#include "Minuit2/MnContours.h"
+#include "Minuit2/MnPlot.h"
+#include "Minuit2/MinosError.h"
+#include "Minuit2/ContoursError.h"
+
+using namespace ROOT::Minuit2;
+
+
+void DumpMinuitParam(const MnUserParameters& upar) {
+  int pr = cout.precision();
+#define PRECISION 13
+#define WIDTH     12
+
+  for(std::vector<MinuitParameter>::const_iterator ipar = upar.Parameters().begin(); 
+      ipar != upar.Parameters().end(); ipar++) {
+    
+    cout << std::setw(4) << (*ipar).Number() << std::setw(5) << "||"; 
+    cout << std::setw(10) << (*ipar).Name()   << std::setw(3) << "||";
+    if((*ipar).IsConst()) {
+      cout << "  const  ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
+    } else if((*ipar).IsFixed()) {
+      cout << "  fixed  ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
+    } else if((*ipar).HasLimits()) {
+      if((*ipar).Error() > 0.) {
+    cout << " limited ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value();
+    cout << " || " << std::setw(12) << (*ipar).Error() << " ||";
+    cout << " [" << ((*ipar).HasLowerLimit()?Num2String((*ipar).LowerLimit()):" - ") 
+         << ", " << ( (*ipar).HasUpperLimit()?Num2String((*ipar).UpperLimit()):" - ") << "]" << std::endl;
+      } else
+    cout << "  free   ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() 
+         << " || " << std::setw(12) << "no" << std::endl;
+    } else {
+      if((*ipar).Error() > 0.)
+    cout << "  free   ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() 
+         << " || " << std::setw(12) << (*ipar).Error() << std::endl;
+      else
+    cout << "  free   ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() 
+         << " || " << std::setw(12) << "no" << std::endl;
+    }
+      
+  }//loop on parameter
+  cout << std::endl;
+  cout.precision(pr);
+#undef PRECISION
+#undef WIDTH
+
+  
+}
+4) Classe de minimisation qui herrite de FCNBase. Par exemple 
+class AutoCorrChi2 : public FCNBase {
+public:
+  
+  AutoCorrChi2(GeneralFunction* func, GeneralFitData* data): func_(func), errDef_(1.) {
+    npts_ = data->NData();
+    for(size_t i=0;i<npts_;i++){
+      x_.push_back(data->X1(i));
+      y_.push_back(data->Val(i));
+      erry_.push_back(data->EVal(i));
+      bool flag = data->IsValid(i) == 1 ? true: false;
+      valid_.push_back(flag);
+    }
+  }
+     
+  //Dtor
+  virtual ~AutoCorrChi2() {}
+  
+  //Specialization of the baseclass
+  virtual r_8 operator() (const vector< r_8 > &par) const { 
+
+    r_8 chi2 =0.;
+
+    for (size_t i=0;i<npts_;i++){
+      if(valid_[i]) {
+
+    r_8 xp[1]; xp[0] = x_[i];
+
+    r_8 val = (func_->Value(xp, par.data()) - y_[i])/erry_[i];
+    chi2 += val*val;
+                      
+      }
+    }
+    
+    return chi2;
+  }
+
+  virtual r_8 Up() const { return errDef_;}
+  
+private:
+  
+  size_t npts_; //nbre of points in the fit
+  GeneralFunction* func_;  //function to fit
+  vector<r_8> x_;          //position values
+  vector<r_8> y_;          //measurements
+  vector<r_8> erry_;       //error on measurements
+  vector<bool> valid_;     // if the measurement is valid or not
+  r_8 errDef_;             // 1: 1-sigma Chi2, n*n for n-sigma errors
+}; //AutoCorrChi2
+
+5) Ensuite il faut instancier la fonction qui sera donner a Minuit
+AutoCorrChi2 fitf(&funcAuto,&dataFit);
+
+funcAuto est la GeneralFunction au sens de Sophya
+
+dataFit sont les donnes au sens de GeneralfitData de Sophya
+
+6) La declartion des parametres
+
+ MnUserParameters upar;
+
+upar.Add(nom, initial, err, min, max)
+exemple upar.Add("Deff",4.5,0.1,3.,6.);
+
+si on veut fixer un parametre > upar.Fix(nom)
+
+au contraire upar.Release(nom) pour le reliberer
+7) Migrad instantiation et minimisation
+      MnMigrad migrad(fitf, upar);
+      unsigned int nparFree = migrad.VariableParameters();
+      cout << "Nbre of free parameters = " << nparFree << endl;
+
+      unsigned int maxfcn = 200+100*nparFree + 100*nparFree*nparFree; 
+      FunctionMinimum min = migrad(maxfcn,0.001);
+      cout << "Minimizer..." << endl;
+       
+      cout << min << endl;
+
+8) recuperation des valeurs fittees
+
+Test si le fit a converge
+      if(!min.IsValid()) {
+    rc = -1; //TODO be more explicit
+    cout << "Minuit FAIL: Freq " <<  freqCentre << endl;
+      }
+ensuite
+min.UserParameters().Value("Deff");
+
+et si on veut l'erreur
+min.UserParameters().Error("Deff");
+
+
+#endif
+