adjustments in preparation for the implementation of ceres-solver

cards/fitPV19.yaml

@@ -10,10 +10,10 @@ qToQmax: 0.2
 # number and order of parameters matches those expected by the
 # particular parameterisation being used.
 Parameters:
-  - {name: "g2",     starting_value: 0.13,  step: 0.001, fix: false}#, lower_bound: 0, upper_bound: 10}
-  - {name: "Npt",    starting_value: 0.285, step: 0.003, fix: false}#, lower_bound: 0, upper_bound: 10}
-  - {name: "alpha",  starting_value: 2.98,  step: 0.030, fix: false}#, lower_bound: 0, upper_bound: 10}
-  - {name: "sigma",  starting_value: 0.173, step: 0.002, fix: false}#, lower_bound: 0, upper_bound: 10}
-  - {name: "lambda", starting_value: 0.39,  step: 0.005, fix: true }#, lower_bound: 0, upper_bound: 10}
-  - {name: "delta",  starting_value: 0,     step: 0.002, fix: false}#, lower_bound: 0, upper_bound: 10}
+  - {name: "g2",     starting_value: 0.03, step: 0.001, fix: false}#, lower_bound: 0, upper_bound: 10}
+  - {name: "Npt",    starting_value: 0.1,  step: 0.003, fix: false}#, lower_bound: 0, upper_bound: 10}
+  - {name: "alpha",  starting_value: 3.7,  step: 0.010, fix: false}#, lower_bound: 0, upper_bound: 10}
+  - {name: "sigma",  starting_value: 0.7,  step: 0.007, fix: false}#, lower_bound: 0, upper_bound: 10}
+  - {name: "lambda", starting_value: 0,    step: 0.005, fix: true }#, lower_bound: 0, upper_bound: 10}
+  - {name: "delta",  starting_value: 0,    step: 0.002, fix: true }#, lower_bound: 0, upper_bound: 10}
 

inc/NangaParbat/DWS.h

@@ -100,12 +100,13 @@ namespace NangaParbat
 
       // Correction to 'b' due to the bmin prescription. Notice that
       // the non-perturbative function does not tend to one anymore
-      // for 'b' tending to zero.
-      const double bmin  = 2 * exp( - apfel::emc) / sqrt(zeta);
-      const double power = 4;
-      const double bc    = b / pow( ( 1 - exp( - pow(b / bmin, power) ) ), 1 / power);
+      // for 'b' tending to zero. (Not sure this is correct).
+      //const double bmin  = 2 * exp( - apfel::emc) / sqrt(zeta);
+      //const double power = 4;
+      //const double bc    = b / pow( ( 1 - exp( - pow(b / bmin, power) ) ), 1 / power);
 
-      return exp( - ( g1 + g2 * log(zeta / Q02) / 2 ) * bc * bc ) * ( 1 - lambda * pow(g1 * bc / 2, 2) / ( 1 + lambda * g1 ) );
+      //return exp( - ( g1 + g2 * log(zeta / Q02) / 2 ) * bc * bc ) * ( 1 - lambda * pow(g1 * bc / 2, 2) / ( 1 + lambda * g1 ) );
+      return exp( - ( g1 + g2 * log(zeta / Q02) / 2 ) * b * b ) * ( 1 - lambda * pow(g1 * b / 2, 2) / ( 1 + lambda * g1 ) );
     };
   };
 }

inc/NangaParbat/chisquare.h

@@ -39,6 +39,14 @@ namespace NangaParbat
      */
     void AddBlock(std::pair<DataHandler,ConvolutionTable> const& DSBlock);
 
+    /**
+     * @brief Function that returns the residuals of the chi2 deriving
+     * from the Cholesky decomposition of the covariance matrix.
+     * @param ids: the dataset index
+     * @return the vector of residuals
+     */
+    std::vector<double> GetResiduals(int const& ids) const;
+
     /**
      * @brief Function that evaluates the &chi;<SUP>2</SUP>'s
      * @param ids: the dataset index (default: -1, the global &chi;<SUP>2</SUP> is computed)

inc/NangaParbat/minimisation.h

@@ -35,6 +35,7 @@ namespace NangaParbat
     double operator()(const std::vector<double>& pars) const;
 
     double Up() const { return 4; }
+
   private:
     mutable ChiSquare _chi2; //!< The "ChiSquare" object that returns the values of all chi2's
   };

run/ComputeChi2.cc

@@ -18,12 +18,14 @@ int main(int argc, char* argv[])
 
   // Allocate "Parameterisation" derived object
   //NangaParbat::DWS NPFunc{};
-  NangaParbat::PV17 NPFunc{};
+  //NangaParbat::PV17 NPFunc{};
+  NangaParbat::PV19 NPFunc{};
 
   // Define "ChiSquare" object with a given qT / Q cut
-  const double qToQmax = 0.3;
+  const double qToQmax = 0.2;
   NangaParbat::ChiSquare chi2{NPFunc, qToQmax};
   //chi2.SetParameters({0.065, 0.285, 2.98, 0.173, 0.39});
+  //chi2.SetParameters({0.01245723634707, 0.1180391594513, 3.418214687464, 1.050830257625, 0, 0});
 
   // Open datasets.yaml file that contains the list of tables to be
   // produced and push data sets into the a vector of DataHandler
@@ -48,11 +50,11 @@ int main(int argc, char* argv[])
   // Get number of data points for each experiment
   const std::vector<int> ndata = chi2.GetDataPointNumbers();
 
-  // // Compute total chi2 with the initial parameters
-  // std::cout << "\nTotal chi2 = " << chi2() << std::endl;
-  // for (int iexp = 0; iexp < (int) ndata.size(); iexp++)
-  //   std::cout << iexp << ") Partial chi2 / #d.p.= " << chi2(iexp) << " (#d.p = " << ndata[iexp] << ")" << std::endl;
-  // std::cout << "\n";
+  // Compute total chi2 with the initial parameters
+  std::cout << "\nTotal chi2 = " << chi2() << std::endl;
+  for (int iexp = 0; iexp < (int) ndata.size(); iexp++)
+    std::cout << iexp << ") Partial chi2 / #d.p.= " << chi2(iexp) << " (#d.p = " << ndata[iexp] << ")" << std::endl;
+  std::cout << "\n";
 
   // // Update parameters and recompute chi2's
   // chi2.SetParameters({1e-3, 1e-3, 1e-3, 1e-3});

run/RunFit.cc

@@ -99,8 +99,18 @@ int main(int argc, char* argv[])
   // Output of Minuit
   std::cout << "minimum: " << min << std::endl;
 
-  // Now print (This also produces plots in pdf with ROOT)
-  std::cout << "Total chi2 = " << chi2() << std::endl;
+  // Now print the total chi2
+  std::cout << "Total chi2 = " << chi2() << "\n" << std::endl;
+
+  // Get number of data points for each experiment
+  const std::vector<int> ndata = chi2.GetDataPointNumbers();
+
+  // Print individual chi2's
+  for (int iexp = 0; iexp < (int) ndata.size(); iexp++)
+    std::cout << iexp << ") Partial chi2 / #d.p.= " << chi2(iexp) << " (#d.p. = " << ndata[iexp] << ")" << std::endl;
+  std::cout << "\n";
+
+  // Finally, this also produces plots
   std::cout << chi2;
 
   // Delete parameterisation

src/chi2/chisquare.cc

@@ -59,6 +59,38 @@ namespace NangaParbat
     _ndata.push_back(idata - (kin.IntqT ? 1 : 0));
   };
 
+  //_________________________________________________________________________________
+  std::vector<double> ChiSquare::GetResiduals(int const& ids) const
+  {
+    if (ids < 0 || ids >= _DSVect.size())
+      throw std::runtime_error("[ChiSquare::GetResiduals]: index out of range");
+
+    // Get "DataHandler" and "ConvolutionTable" objects
+    const DataHandler      dh = _DSVect[ids].first;
+    const ConvolutionTable ct = _DSVect[ids].second;
+
+    // Get experimental central values
+    const std::vector<double> mean = dh.GetMeanValues();
+
+    // Get predictions
+    auto const fNP = [&] (double const& x, double const& b, double const& zeta, int const& ifun) -> double{ return _NPFunc.Evaluate(x, b, zeta, ifun); };
+    const std::vector<double> pred = ct.GetPredictions(fNP);
+
+    // Check that the number of points in the DataHandler and
+    // Convolution table objects is the same.
+    if (mean.size() != pred.size())
+      throw std::runtime_error("[ChiSquare::GetResiduals]: mismatch in the number of points");
+
+    // Compute residuals only for the points that pass the cut qT
+    // / Q, set the others to zero.
+    std::vector<double> res(_ndata[ids], 0.);
+    for (int j = 0; j < _ndata[ids]; j++)
+      res[j] = mean[j] - pred[j];
+
+    // Solve lower-diagonal system and return the result
+    return SolveLowerSystem(dh.GetCholeskyDecomposition(), res);
+  }
+
   //_________________________________________________________________________________
   double ChiSquare::Evaluate(int const& ids) const
   {
@@ -70,8 +102,6 @@ namespace NangaParbat
         istart = ids;
         iend   = ids + 1;
       }
-    else if (ids >= iend)
-      throw std::runtime_error("[ChiSquare::Evaluate]: index out of range");
 
     // Initialise chi2 and number of data points
     double chi2 = 0;
@@ -80,30 +110,8 @@ namespace NangaParbat
     // Loop over the the blocks
     for (int i = istart; i < iend; i++)
       {
-        // Get "DataHandler" and "ConvolutionTable" objects
-        const DataHandler      dh = _DSVect[i].first;
-        const ConvolutionTable ct = _DSVect[i].second;
-
-        // Get experimental central values
-        const std::vector<double> mean = dh.GetMeanValues();
-
-        // Get predictions
-        auto const fNP = [&] (double const& x, double const& b, double const& zeta, int const& ifun) -> double{ return _NPFunc.Evaluate(x, b, zeta, ifun); };
-        const std::vector<double> pred = ct.GetPredictions(fNP);
-
-        // Check that the number of points in the DataHandler and
-        // Convolution table objects is the same.
-        if (mean.size() != pred.size())
-          throw std::runtime_error("[ChiSquare::Evaluate]: mismatch in the number of points");
-
-        // Compute residuals only for the points that pass the cut qT
-        // / Q, set the others to zero.
-        std::vector<double> res(_ndata[i], 0.);
-        for (int j = 0; j < _ndata[i]; j++)
-          res[j] = mean[j] - pred[j];
-
-        // Solve lower-diagonal system
-        const std::vector<double> x = SolveLowerSystem(dh.GetCholeskyDecomposition(), res);
+        // Get residuals
+        const std::vector<double> x = GetResiduals(i);
 
         // Compute contribution to the chi2 as absolute value of "x"
         chi2 += std::inner_product(x.begin(), x.end(), x.begin(), 0.);
@@ -220,10 +228,8 @@ namespace NangaParbat
           }
         exp->SetLineColor(1);
         exp->SetMarkerStyle(20);
-        // theo->SetLineColor(2);
         theo->SetLineColor(kBlue-7);
         theo->SetLineWidth(3);
-        // theoshift->SetLineColor(3);
         theoshift->SetLineColor(kPink-6);
         theoshift->SetLineWidth(3);