Fixing compatibility problem with NangaParbat-Beta

MapCollaboration · Dec 7, 2022 · fc23e37 · fc23e37
1 parent e7f86a0
commit fc23e37
Show file tree

Hide file tree

Showing 4 changed files with 8 additions and 37 deletions.
diff --git a/inc/MontBlanc/predictionshandler.h b/inc/MontBlanc/predictionshandler.h
@@ -63,6 +63,5 @@ namespace MontBlanc
     std::vector<double>                                  _ChargeMap;
     std::vector<apfel::Set<apfel::Operator>>             _FKt;
     apfel::Set<apfel::Distribution>                      _D;
-    apfel::Set<apfel::Distribution>                      _Do;
   };
 }
diff --git a/run/ComputeChi2s.cc b/run/ComputeChi2s.cc
@@ -163,7 +163,7 @@ void compute_chi2s(std::string ResultFolder, int member_index, std::string LHAPD
       std::vector<std::shared_ptr<NangaParbat::Cut>> cuts;
       for (auto const& c :ds["cuts"])
         cuts.push_back(NangaParbat::CutFactory::GetInstance(*DH, c["name"].as<std::string>(), c["min"].as<double>(), c["max"].as<double>(),
-							    (c["pars"] ? c["pars"].as<std::vector<double>>() : std::vector<double>{})));
+                                                            (c["pars"] ? c["pars"].as<std::vector<double>>() : std::vector<double> {})));
 
       // Predictions
       NangaParbat::ConvolutionTable *PH = new MontBlanc::PredictionsHandler{config["Predictions"], *DH, g, cuts};

diff --git a/run/Optimize.cc b/run/Optimize.cc
@@ -149,16 +149,16 @@ int main(int argc, char *argv[])
       // unfluctuated data set.
       NangaParbat::DataHandler *DHc = new NangaParbat::DataHandler{ds["name"].as<std::string>(), YAML::LoadFile(DataFolder + ds["file"].as<std::string>())};
 
-      // Check whethe the hadronic species set in the configuration card
-      // matches with that of the data set.
+      // Check whether the hadronic species set in the configuration
+      // card matches with that of the data set.
       if (DH->GetHadron() != hadron)
         throw std::runtime_error("This data set corresponds to a hadronic species different from that set in the input configuration card");
 
       // Accumulate kinematic cuts
       std::vector<std::shared_ptr<NangaParbat::Cut> > cuts;
       for (auto const &c : ds["cuts"])
         cuts.push_back(NangaParbat::CutFactory::GetInstance(*DH, c["name"].as<std::string>(), c["min"].as<double>(), c["max"].as<double>(),
-							    (c["pars"] ? c["pars"].as<std::vector<double>>() : std::vector<double>{})));
+                                                            (c["pars"] ? c["pars"].as<std::vector<double>>() : std::vector<double> {})));
 
       // Compute predictions within kinematic cuts
       MontBlanc::PredictionsHandler PH{config["Predictions"], *DH, g, cuts};

diff --git a/src/predictions/predictionshandler.cc b/src/predictions/predictionshandler.cc
@@ -76,7 +76,6 @@ namespace MontBlanc
       {*(BuildDglap(InitializeDglapObjectsQCDT(*_g, _Thresholds, true), _mu0, PerturbativeOrder, Alphas)), 100, 1, 100, 3}};
 
     // Zero operator
-    const apfel::Operator Id  {*_g, apfel::Identity{}};
     const apfel::Operator Zero{*_g, apfel::Null{}};
 
     // Set cuts in the mother class
@@ -126,7 +125,7 @@ namespace MontBlanc
             if (PerturbativeOrder > 1)
               Ki += pow(as / apfel::FourPi, 2) * F2Obj.C2.at(comp);
 
-            // Convolution coefficient functions with the evolution
+            // Convolute coefficient functions with the evolution
             // operators
             for (int j = 0; j < 13; j++)
               {
@@ -341,8 +340,7 @@ namespace MontBlanc
                 Qmin = std::max(sqrt(_bins[i].xmin * _bins[i].ymin) * Vs, DH.GetKinematics().var1b.first);
                 Qmax = sqrt(_bins[i].xmax * _bins[i].ymax) * Vs;
               }
-            else if (_obs == NangaParbat::DataHandler::Observable::dsigma_dxdQdz ||
-                     _obs == NangaParbat::DataHandler::Observable::opposite_sign_ratio)
+            else if (_obs == NangaParbat::DataHandler::Observable::dsigma_dxdQdz)
               {
                 Qmin = _bins[i].Qmin;
                 Qmax = _bins[i].Qmax;
@@ -379,8 +377,7 @@ namespace MontBlanc
                   if (DH.GetKinematics().PSRed)
                     xbmax = std::min(xbmax, 1 / ( 1 + pow(DH.GetKinematics().pTMin / Q, 2) ));
                 }
-              else if (_obs == NangaParbat::DataHandler::Observable::dsigma_dxdQdz ||
-                       _obs == NangaParbat::DataHandler::Observable::opposite_sign_ratio)
+              else if (_obs == NangaParbat::DataHandler::Observable::dsigma_dxdQdz)
                 {
                   xbmin = _bins[i].xmin;
                   xbmax = _bins[i].xmax;
@@ -408,8 +405,7 @@ namespace MontBlanc
                   if (DH.GetKinematics().PSRed)
                     xbmax = std::min(xbmax, 1 / ( 1 + pow(DH.GetKinematics().pTMin / Q, 2) ));
                 }
-              else if (_obs == NangaParbat::DataHandler::Observable::dsigma_dxdQdz ||
-                       _obs == NangaParbat::DataHandler::Observable::opposite_sign_ratio)
+              else if (_obs == NangaParbat::DataHandler::Observable::dsigma_dxdQdz)
                 {
                   xbmin = _bins[i].xmin;
                   xbmax = _bins[i].xmax;
@@ -509,18 +505,6 @@ namespace MontBlanc
   {
     // Construct set of distributions
     _D = apfel::Set<apfel::Distribution> {_ChargeMap * InDistFunc(_mu0)};
-
-    // If we need to compute also predictions with opposite sign, also
-    // allocate _Do.
-    if (_obs == NangaParbat::DataHandler::Observable::opposite_sign_ratio)
-      {
-        // Reverse sign
-        std::vector<double> ChargeMapo = _ChargeMap;
-        for (int i = 1; i <= 6; i++)
-          ChargeMapo[2 * i] *= - 1;
-
-        _Do = apfel::Set<apfel::Distribution> {ChargeMapo * InDistFunc(_mu0)};
-      }
   }
 
   //_________________________________________________________________________
@@ -539,18 +523,6 @@ namespace MontBlanc
       else
         preds[id] = (_cutmask[id] ? (_FKt[id] * _D).Combine().Evaluate(_bins[id].zav) / _bins[id].zav * _qTfact[id] : 0);
 
-    // If needed, compute also predictions with opposite sign using
-    // _Do and take the ratio. If a point does not pass the cut divide
-    // by 1. This is just to avoid dividing zero by another zero.
-    if (_obs == NangaParbat::DataHandler::Observable::opposite_sign_ratio)
-      {
-        for (int id = 0; id < (int) _bins.size(); id++)
-          if (_bins[id].Intz)
-            preds[id] /= (_cutmask[id] ? ((_FKt[id] * _Do).Combine() * [] (double const& z) -> double{ return 1 / z; }).Integrate(_bins[id].zmin, _bins[id].zmax)
-                          / ( _bins[id].zmax - _bins[id].zmin ) * _qTfact[id] : 1);
-          else
-            preds[id] /= (_cutmask[id] ? (_FKt[id] * _Do).Combine().Evaluate(_bins[id].zav) / _bins[id].zav * _qTfact[id] : 1);
-      }
     return preds;
   }