Implementation fit avec minuit pour le fit des baselines et decodage des...

Implementation fit avec minuit pour le fit des baselines et decodage des options -nocxf -phifreq ds main, Reza 26/02/2019

gcxfitbaseline.h

@@ -12,6 +12,8 @@
 #include "cxbeam.h"
 #include "trkfit.h"
 
+#include "tkfit_minuit.h"
+
 class My6CxSignalsB {
 public:
   My6CxSignalsB(vector< AcxDataSet> & acxd, vector< TrackSet > & trks, bool fgno_aphi=false)
@@ -24,18 +26,21 @@ public:
     n_cxcor=acxd_[0].getNbCrossCor();
     v_xi2.resize(acxd_[0].getNbCrossCor());
     v_npts_xi2.resize(acxd_[0].getNbCrossCor());
+    ndatapoints_=0;
     for(size_t kcx=0; kcx<n_cxcor; kcx++) {
       v_npts_xi2[kcx]=0;
       for(size_t i=0; i<acxd_.size(); i++) 
-      v_npts_xi2[kcx]+=acxd_[i].v_time_data.size();
+	v_npts_xi2[kcx]+=acxd_[i].v_time_data.size();
+      ndatapoints_+=v_npts_xi2[kcx];
     }
+    ndatapoints_*=2;
     //DBG    cout << " My6CxSignalsB() "<< endl;
   }
 
   inline void SetPrintLevel(int lev=0) { prtlevel_=lev; }
 
   // get the the TimeVec in dataset i , corresponding to the j-th track  
-  virtual TVector<double> getTimeVec(size_t i, size_t j)
+  virtual TVector<double> getTimeVec(size_t i, size_t j)  const 
   {
     vector<double> & tvec = acxd_[i].v_time_data[j];
     TVector<double>  tmv(tvec);
@@ -43,7 +48,7 @@ public:
   }
 
   // get the the signal in dataset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5) 
-  virtual int getDataSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig) 
+  virtual int getDataSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig)  const 
   {
     //DEL    sig.resize(v_time_data[j].size());
     vector< complex<double> > & vcxdata =  acxd_[i].vv_cxdata[j][kcx];
@@ -53,7 +58,7 @@ public:
   }
 
   // get the j-th signal for cross-correlation kcx (0<=kcx<=5) 
-  virtual TVector< complex<double> > getDataSignal(size_t i, size_t j, size_t kcx)
+  virtual TVector< complex<double> > getDataSignal(size_t i, size_t j, size_t kcx)  const 
   {
     vector< complex<double> >sig;
     getDataSignal(i, j, kcx, sig);
@@ -63,21 +68,21 @@ public:
   }
 
   // get the beam for the i-th data set and the cross-correlation kcx (0<=kcx<=5)
-  inline CxBeam & getBeam(size_t i, size_t kcx)
+  inline CxBeam const & getBeam(size_t i, size_t kcx)   const 
   {
     return acxd_[i].v_cxbeams[kcx];
   }
 
   //---- This method should be identical to the one in MyCxSignal (gcxfit.h) 
   // Model for linearly varying phase with frequency : Phase = Phi0 + a_phi * (freq-1250)/250 
-  inline double getPhase4Freq(double phi0, double a_phi, double freq)
+  inline double getPhase4Freq(double phi0, double a_phi, double freq)  const 
   {
     return (phi0+a_phi*(freq-1250.)*phlinfac_);
   }
 
   // get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5) 
   virtual int getExpectedSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig, 
-				double phi0, double aphi, Vector3d & baseline_shift)
+				double phi0, double aphi, Vector3d & baseline_shift)  const 
   {
     //DBG   cout<<"*DBG*A*My6CxSignalsB::getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<" kcx="<<kcx<<endl;
     vector<double> & vtime =  acxd_[i].v_time_data[j];
@@ -133,7 +138,7 @@ public:
      parm[14] = Z_shift_4  
   */
   // determine la phase et le baseline_shift a partir du tableau parm, pour la ligne de base kcx (0<=kcx<=5)
-  void getFromFitParam(size_t kcx, const double* parm, double & phi0, double & aphi, Vector3d & baseline_shift)
+  void getFromFitParam(size_t kcx, const double* parm, double & phi0, double & aphi, Vector3d & baseline_shift)  const 
   {
     size_t na[6]={0,0,0,1,1,2};
     size_t nb[6]={1,2,3,2,3,3};
@@ -175,7 +180,7 @@ public:
     return;
   }
   // get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5) 
-  virtual int getExpectedSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig, const double* parm)
+  virtual int getExpectedSignal(size_t i, size_t j, size_t kcx, vector< complex<double> > & sig, const double* parm)  const 
   {
     double phi0,aphi;
     Vector3d baseline_shift;
@@ -184,7 +189,7 @@ public:
   }
 
   // get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5) 
-  virtual TVector< complex<double> > getExpectedSignal(size_t i, size_t j, size_t kcx, double phi0, double aphi, Vector3d & baseline_shift)
+  virtual TVector< complex<double> > getExpectedSignal(size_t i, size_t j, size_t kcx, double phi0, double aphi, Vector3d & baseline_shift)  const 
   {
     vector< complex<double> > sig;
     //DBG    cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
@@ -195,7 +200,7 @@ public:
   }
 
   // get the expected signal for trackset i , corresponding to the j-th track for cross-correlation kcx (0<=kcx<=5) 
-  virtual TVector< complex<double> > getExpectedSignal(size_t i, size_t j, size_t kcx, const double* parm)
+  virtual TVector< complex<double> > getExpectedSignal(size_t i, size_t j, size_t kcx, const double* parm)  const 
   {
     vector< complex<double> > sig;
     //DBG    cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
@@ -205,7 +210,7 @@ public:
     return rvs;
   }
 
-  virtual double getXi2(const double* parm, int & npts) 
+  virtual double getXi2(const double* parm, int & npts)   const 
   {
     double xi2=0.;
     npts=0;
@@ -245,11 +250,12 @@ public:
   size_t n_cxcor;     // Number of Cx correlations ( 6 for PAON4 )
   vector< AcxDataSet> & acxd_;
   vector< TrackSet > & trk_;
-  vector< double > v_xi2;
+  mutable vector< double > v_xi2;
   vector< size_t > v_npts_xi2;
   int prtlevel_;
   double phlinfac_;
-  bool fgno_aphi_;   // if true -> pas de terme a_phi (variation lineaire de phase avec al frequence ds le tableau param de getXi2 
+  bool fgno_aphi_;   // if true -> pas de terme a_phi (variation lineaire de phase avec al frequence ds le tableau param de getXi2
+  size_t ndatapoints_; 
 }; 
 
 //-----------------------------------------------------------------
@@ -261,6 +267,9 @@ public:
   {
   }
 
+  size_t getNbParams() const { return NPar(); }
+  size_t getNbDataPoints() const { return ndatapoints_; }
+
   /*    definition du tableau des parametres de fit 
      parm[0] = Phi0_2          Phi0 Antenne 2 
      parm[1] = a_Phi_2         aPhi (slope)  Antenne 2 
@@ -302,4 +311,44 @@ public:
   }
 
 };
+
+#ifndef TKF_AVEC_MINUIT
+typedef My6CxGenXi2B  TkF_6CxXi2B ;
+#else 
+// Ajustement avec Minuit 
+class My6Cx_Xi2B_Minuit : public Trk_FCNBase , public My6CxSignalsB {
+public:
+  My6Cx_Xi2B_Minuit(vector< AcxDataSet> & acxd, vector< TrackSet > & trks)
+    : My6CxSignalsB(acxd, trks)
+  {
+    nbparam_=15;   // nb_param= 15 = 3 * 2 (phi0,aphi) + 3*3 position shifts
+    myparam_ = new double[nbparam_];
+  }
+  virtual ~My6Cx_Xi2B_Minuit() 
+  {
+    delete[] myparam_;
+  }
+
+  virtual size_t getNbParams() const { return nbparam_; }
+  virtual size_t getNbDataPoints() const { return ndatapoints_; }
+
+  virtual double operator() (const std::vector< double > & p) const 
+  {
+    if (p.size() !=  nbparam_)  throw ParmError("My6Cx_Xi2B_Minuit::operator()(p)  p.size() != nbparam_");
+    for(size_t i=0; i<nbparam_; i++) myparam_[i]=p[i];
+    int ndataused=0;
+    return getXi2(myparam_, ndataused);
+  }
+
+  virtual r_8 Up() const { return errDef_;}
+  
+  size_t nbparam_;
+  double * myparam_;
+  double errDef_;
+};
+
+typedef My6Cx_Xi2B_Minuit  TkF_6CxXi2B ;
+#endif   /*  TKF_AVEC_MINUIT */
+
+
 #endif

tkfit_minuit.h

@@ -82,6 +82,12 @@ public:
     return;
   }
 
+  void SetFix(size_t num)
+  {
+    upar_.Fix((unsigned int)num);
+    return;
+  }
+
   void SetMaxStep(int maxstep) { max_step_=maxstep; }
 
   static string Num2String(double v) 

trkacxfit.cc

@@ -26,7 +26,8 @@ csh> ../AnaPaon4/Objs/trkacxfit -D 4.5 -out A21_fac.txt  trk_1.d
 @zenang ZenithAngle 
 ### Input data and track files - multiple @trk cards 
 # filenames with .ppf extension , Tstart,Tend in minutes 
-@trk AutoCrossVisiDataTableFile  Tstart,Tend  Freq  TrackFile 
+@trk AutoCrossVisiDataTableFile  Tstart,Tend  Freq  TrackFile [NOAC NOCX  NOACX]
+### Optional last argument [ NOAC NOCX  NOACX ] not used currently 
  ----------------------------------------------------------------------------*/
 
 #include <stdlib.h>
@@ -88,6 +89,9 @@ static int prtlevel=0;  // print level
 static string default_input_path;
 static vector<string> trksetfiles;   // datacard files defining the track/data sets 
 
+static bool do_cxfit = true ;       // if false, perform AutoCor fit only
+static bool fg_phi0only = true;     // if false, linear varying phase with frequency, in AutoCor and Cross-Cor fits
+
 // ---- fit simultane sur les 6 cross-correlations 
 static bool do_baselinefit = false;   // if true , perform baseline fit 
 static bool fg_fixbaseline = false;   // if true , do phases fit with fixed baselines 
@@ -130,32 +134,35 @@ int main (int narg, char* arg[])
       acxfit.doACfit(acofile);  
       string acckfile = fnbuff+checkfilename;
       acxfit.saveExpectedAC(acckfile);
-      sprintf(fnbuff,"cx_%d_",(int)(i+1));
-      string cxofile = fnbuff+outfilename;
-      string cxckfile = fnbuff+checkfilename;
-      acxfit.doCxfit(cxofile, fguseAac4Cx); 
-      acxfit.saveExpectedCx(cxckfile);
-
+      if (do_cxfit)  {  //  Fit des cross-cor aussi 
+	sprintf(fnbuff,"cx_%d_",(int)(i+1));
+	string cxofile = fnbuff+outfilename;
+	string cxckfile = fnbuff+checkfilename;
+	acxfit.doCxfit(cxofile, fguseAac4Cx); 
+	acxfit.saveExpectedCx(cxckfile);
+      }
       v_acxd.push_back(acxd);
       v_trk.push_back(tks);
     }
 
-    string cxbofile = "cxb6_"+outfilename;
-    string cxbckfile = "cxb6_"+checkfilename;
-    CxBaselineFitter cxbfit(v_acxd, v_trk);
-    //cxbfit.doSimplexMinimize();
-    cxbfit.initFitParams();
-    cxbfit.saveExpectedCx(cxbckfile);
-    if (do_baselinesimplex) {
-      cxbfit.doSimplexMinimize();
-      cxbckfile="cxb6s_"+checkfilename;
-      cxbfit.saveExpectedCx(cxbckfile);
-      //      cxbfit.doCheck();
-    }
-    if (do_baselinefit) {
-      cxbfit.dofit(cxbofile,fg_fixbaseline);
-      cxbckfile = "cxb6f_"+checkfilename;
+    if (do_baselinefit || do_baselinesimplex)  {   // fit simultane des 6 cross-cor 
+      string cxbofile = "cxb6_"+outfilename;
+      string cxbckfile = "cxb6_"+checkfilename;
+      CxBaselineFitter cxbfit(v_acxd, v_trk);
+      //cxbfit.doSimplexMinimize();
+      cxbfit.initFitParams();
       cxbfit.saveExpectedCx(cxbckfile);
+      if (do_baselinesimplex) {
+	cxbfit.doSimplexMinimize();
+	cxbckfile="cxb6s_"+checkfilename;
+	cxbfit.saveExpectedCx(cxbckfile);
+	//      cxbfit.doCheck();
+      }
+      if (do_baselinefit) {
+	cxbfit.dofit(cxbofile,fg_fixbaseline);
+	cxbckfile = "cxb6f_"+checkfilename;
+	cxbfit.saveExpectedCx(cxbckfile);
+      }
     }
     Timer tm("trkacxfit");
   }  // End of try bloc 
@@ -192,7 +199,8 @@ int decode_args(int narg, char* arg[])
     cout << "  trkacxfit : fit array geometry and   \n"
 	 << "  Usage: trkacxfit [-options]  track_Set_1 [track_Set_2]  [track_Set_3 ...] \n"
 	 << "  options: -inp def_input_path -out OutFileName -ckf CheckFileName \n" 
-	 << "           -docx6f -fixb -docx6s -ngb  -D dish_diameter -prt PrintLevel  \n"<<endl;
+	 << "           -nocxf -phifreq -docx6f -fixb -docx6s  \n"
+	 <<"            -D dish_diameter -ngb  -prt PrintLevel  \n"<<endl;
     if (!fglonghelp) {
       cout << " trkacxfit -h  to get option description " << endl;
       return 2;
@@ -203,6 +211,9 @@ int decode_args(int narg, char* arg[])
       	 << "  -ckf CheckFileName : Output file (PPF) to save expected signals (def chktrkfit.ppf)  \n"
 	 << "    Output file names are constructed from OutFileName and CheckFileName prepending \n" 
 	 << "  ac_JJ_ or cx_JJ_  where JJ=1...n is the track_set number \n"
+      	 << "  -nocxf : Don't perform fit on cross-cor (AutoCor only fit) - default YES -> DO Cx fits) \n"
+	 << "  -phifreq : Use linearly varying phase with frequency model - default No fixed Phi=Phi0  \n"
+	 <<"             Linear phi model : Phi(freq)=Phi_0 + a_Phi*(freq-1250)/250 \n"
 	 << "  -docx6f : Perform simultaneous fit over 6 cross-corr (baseline and phases determination) - default NO  \n"
 	 << "  -fixb : perform the previous fit (docx6f) with fixed baselines - default NO  \n"
 	 << "  -docx6s : try simplex minimisation over 6 cross-corr (baseline and phases determination) - default NO  \n"
@@ -235,6 +246,12 @@ int decode_args(int narg, char* arg[])
     else if (fbo=="-ngb") {  // Use Non gaussian beam 
       fggaussbeam=false;   arg++;  narg--;  lastargs.clear(); 
     }
+    else if (fbo=="-nocxf") {  // DON'T FIT Cross-Correlations 
+      do_cxfit=false;   arg++;  narg--;  lastargs.clear(); 
+    }
+    else if (fbo=="-phifreq") {  // Use linear model for phase(frequency)  - If not Phi(freq) = Phi0  
+      fg_phi0only=false;   arg++;  narg--;  lastargs.clear(); 
+    }
     else if (fbo=="-docx6f") {  // Perform  simultaneous fit over 6 cross-corr 
       do_baselinefit=true;   arg++;  narg--;  lastargs.clear(); 
     }
@@ -268,9 +285,13 @@ int decode_args(int narg, char* arg[])
     return 9;
   }
 
+  if (!do_cxfit) do_baselinefit=do_baselinesimplex=false;
+  
   cout << " ------------------- trkacxfit/run parameters:"<<endl; 
   cout << " Beam: D_dish (initial guess)="<<D_dish<<" GaussianBeam ? "<<(fggaussbeam?"Yes":"No")<<endl;
   cout << " OutFileName= "<<outfilename<<" CheckFileName= "<<checkfilename<<endl;
+  cout << " Perform fits on cross-correlations ? " << (do_cxfit?"Yes":"No")
+       << "   Phase model : Phase(freq)= " << (fg_phi0only?"Phi0":"Phi0+aPhi*(freq-1250)/250") << endl;
   cout << " Perform baseline fit ? " << (do_baselinefit?"Yes":"No")
        <<(fg_fixbaseline?"  With FIXED baselines":" (Phases & baselines)")
        <<" Simplex-Minimisation ? "<<(do_baselinesimplex?"Yes":"No")<<endl;

trkfit.cc

@@ -903,10 +903,12 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline, bool fgphi0o
       }
     }
 
-  My6CxGenXi2B gxi2(v_acxd, v_trks); 
-  GeneralFit mFit(&gxi2);
+  TkF_6CxXi2B gxi2(v_acxd, v_trks);   // My6CxGenXi2B
+  //  GeneralFit mFit(&gxi2);
+  TkF_Fitter mFit(gxi2);
+
   mFit.SetData(&gdata);        // connect data to the fitter , here the data is unused - gxi2 includes its data 
-  mFit.SetMaxStep(1000);
+  mFit.SetMaxStep(3000);
   
   // SetParam(int n,double value, double step,double min=1., double max=-1.);
   for(size_t i=0; i<(NB_ANTENNES-1); i++) {
@@ -929,7 +931,7 @@ int CxBaselineFitter::dofit(string outfilename, bool fgfixbaseline, bool fgphi0o
     }
   }
   cout << " Performing the fit (tot_npoints_fit= "<<tot_npoints_fit<<" ?= (npoints2="<<npoints2<<") ..."<< endl;
-  rcfit = mFit.Fit();  xi2red=-99999.;
+  rcfit = mFit.doFit();  xi2red=-99999.;
   cout<< "------ Fit result Reduce_Chisquare = " << mFit.GetChi2Red()<< " nstep="<<mFit.GetNStep() << " rc="<<rcfit<<endl;
   mFit.PrintFit();