Fix of LPAIR z-dependent distributions in SD mode (#58)

cepgen · Mar 8, 2024 · b6606e2 · b6606e2
1 parent 3a6be7e
commit b6606e2
Show file tree

Hide file tree

Showing 6 changed files with 94 additions and 14 deletions.
diff --git a/CepGen/Core/Generator.cpp b/CepGen/Core/Generator.cpp
@@ -65,6 +65,7 @@ namespace cepgen {
     worker_->setIntegrator(integrator_.get());
     xsect_ = Value{-1., -1.};
     parameters_->prepareRun();
+    initialised_ = false;
   }
 
   const RunParameters& Generator::runParameters() const {
@@ -166,8 +167,6 @@ namespace cepgen {
   }
 
   void Generator::initialise() {
-    if (initialised_)
-      return;
     if (!parameters_)
       throw CG_FATAL("Generator:generate") << "No steering parameters specified!";
 

diff --git a/CepGen/Generator.h b/CepGen/Generator.h
@@ -95,12 +95,8 @@ namespace cepgen {
     RunParameters& runParameters();              ///< Run parameters block
     void setRunParameters(RunParameters* ip);    ///< Feed the generator with a RunParameters object
 
-    void resetIntegrator();  ///< Reset integrator algorithm from the user-specified configuration
     void setIntegrator(std::unique_ptr<Integrator>);  ///< Specify an integrator algorithm configuration
-
-    void clearRun();  ///< Remove all references to a previous generation/run
-
-    void integrate();  ///< Integrate the functional over the phase space of interest
+    void integrate();                                 ///< Integrate the functional over the phase space of interest
 
     /// Compute the cross section for the run parameters
     /// \return The computed cross-section and uncertainty, in pb
@@ -123,7 +119,10 @@ namespace cepgen {
     double computePoint(const std::vector<double>& x);
 
   private:
-    void initialise();                           ///< Initialise event generation
+    void initialise();       ///< Initialise event generation
+    void clearRun();         ///< Remove all references to a previous generation/run
+    void resetIntegrator();  ///< Reset integrator algorithm from the user-specified configuration
+
     std::unique_ptr<RunParameters> parameters_;  ///< Run parameters for event generation and cross-section computation
     std::unique_ptr<GeneratorWorker> worker_;    ///< Generator worker instance
     std::unique_ptr<Integrator> integrator_;     ///< Integration algorithm

diff --git a/CepGen/Utils/Hist1D.cpp b/CepGen/Utils/Hist1D.cpp
@@ -240,5 +240,21 @@ namespace cepgen {
       }
       return gsl_histogram_pdf_sample(pdf_.get(), rng.uniform());
     }
+
+    double Hist1D::chi2test(const Hist1D& oth, size_t& ndfval) const {
+      if (nbins() != oth.nbins())
+        return 0.;
+      double chi2val = 0.;
+      ndfval = nbins();
+      for (size_t i = 0; i < nbins(); ++i) {
+        const auto bin_val1 = value(i), bin_val2 = oth.value(i);
+        if (bin_val1 == 0. && bin_val2 == 0.) {
+          --ndfval;
+          continue;
+        }
+        chi2val += std::pow((double)bin_val1 - (double)bin_val2, 2) / ((double)bin_val1 + (double)bin_val2);
+      }
+      return chi2val;
+    }
   }  // namespace utils
 }  // namespace cepgen
diff --git a/CepGen/Utils/Histogram.h b/CepGen/Utils/Histogram.h
@@ -89,6 +89,11 @@ namespace cepgen {
       /// Sample individual "events" from a distribution
       double sample(RandomGenerator&) const;
 
+      /// Perform a chi^2 test between two histograms
+      /// \param[out] ndf number of degrees of freedom (non-empty bins)
+      /// \return chi^2-value of the equivalence test
+      double chi2test(const Hist1D&, size_t& ndf) const;
+
       /// Retrieve the value + uncertainty for all bins
       std::vector<Value> values() const;
       /// Retrieve the value + uncertainty for one bin

diff --git a/CepGenProcesses/LPAIR.cpp b/CepGenProcesses/LPAIR.cpp
@@ -132,10 +132,6 @@ class LPAIR final : public cepgen::proc::Process {
     // boost of the outgoing beams
     pX().setMass(mX());
     pY().setMass(mY());
-    if (beams_mode_ == mode::Kinematics::ElasticInelastic) {  // mirror X/Y and dilepton systems if needed
-      std::swap(pX(), pY());
-      std::swap(pc(0), pc(1));
-    }
     pX().betaGammaBoost(gamma_cm_, beta_gamma_cm_);
     pY().betaGammaBoost(gamma_cm_, beta_gamma_cm_);
     // incoming partons
@@ -151,6 +147,14 @@ class LPAIR final : public cepgen::proc::Process {
       if (symmetrise_ && mirror)
         mom->mirrorZ();
     }
+    if (beams_mode_ == mode::Kinematics::ElasticInelastic) {  // mirror X/Y and dilepton systems if needed
+      std::swap(pX(), pY());
+      std::swap(pc(0), pc(1));
+      pX().mirrorZ();
+      pY().mirrorZ();
+      pc(0).mirrorZ();
+      pc(1).mirrorZ();
+    }
     // first outgoing beam
     event()
         .oneWithRole(Particle::OutgoingBeam1)

diff --git a/test/physics/symm_singlediss.cc b/test/physics/symm_singlediss.cc
@@ -18,26 +18,30 @@
 
 #include "CepGen/Core/RunParameters.h"
 #include "CepGen/Generator.h"
+#include "CepGen/Modules/DrawerFactory.h"
 #include "CepGen/Modules/ProcessFactory.h"
 #include "CepGen/Modules/StructureFunctionsFactory.h"
 #include "CepGen/Process/Process.h"
 #include "CepGen/Utils/AbortHandler.h"
 #include "CepGen/Utils/ArgumentsParser.h"
+#include "CepGen/Utils/Drawer.h"
 #include "CepGen/Utils/Test.h"
 #include "CepGen/Utils/Timer.h"
 
 using namespace std;
 
 int main(int argc, char* argv[]) {
   double num_sigma;
-  int str_fun;
-  string proc_name, integrator;
+  int num_gen, str_fun;
+  string proc_name, integrator, plotter;
 
   cepgen::ArgumentsParser(argc, argv)
       .addOptionalArgument("process,p", "process to compute", &proc_name, "lpair")
+      .addOptionalArgument("num-gen,g", "number of events to generate", &num_gen, 25'000)
       .addOptionalArgument("num-sigma,n", "max. number of std.dev.", &num_sigma, 3.)
       .addOptionalArgument("str-fun,s", "struct.functions modelling", &str_fun, 11)
       .addOptionalArgument("integrator,i", "type of integrator used", &integrator, "Vegas")
+      .addOptionalArgument("plotter,t", "type of plotter to use", &plotter, "")
       .parse();
 
   cepgen::utils::Timer tmr;
@@ -58,16 +62,69 @@ int main(int argc, char* argv[]) {
   gen.runParameters().setProcess(
       cepgen::ProcessFactory::get().build(proc_name, cepgen::ParametersList().set<int>("pair", 13)));
   cepgen::Value cs_ei, cs_ie;
+
+  auto h_eta_lead_ei = cepgen::utils::Hist1D(50, {-2.5, 2.5}, "eta_lead_ei", "el-inel"),
+       h_eta_lead_ie = cepgen::utils::Hist1D(50, {-2.5, 2.5}, "eta_lead_ie", "inel-el"),
+       h_eta_sublead_ei = cepgen::utils::Hist1D(50, {-2.5, 2.5}, "eta_sublead_ei", "el-inel"),
+       h_eta_sublead_ie = cepgen::utils::Hist1D(50, {-2.5, 2.5}, "eta_sublead_ie", "inel-el"),
+       h_mdiff_ei = cepgen::utils::Hist1D(50, {0., 100.}, "mdiff_ei", "el-inel"),
+       h_mdiff_ie = cepgen::utils::Hist1D(50, {0., 100.}, "mdiff_ie", "inel-el");
+
   {  // elastic-inelastic
     pkin.set<int>("mode", 2);
     gen.runParameters().process().kinematics().setParameters(pkin);
     cs_ei = gen.computeXsection();
+    if (num_gen > 0)
+      gen.generate(num_gen, [&](const cepgen::Event& evt, size_t) {
+        const auto &mom1 = evt(cepgen::Particle::Role::CentralSystem).at(0).momentum(),
+                   &mom2 = evt(cepgen::Particle::Role::CentralSystem).at(1).momentum();
+        if (mom1.pt() > mom2.pt()) {
+          h_eta_lead_ei.fill(mom1.eta());
+          h_eta_sublead_ei.fill(mom2.eta());
+        } else {
+          h_eta_lead_ei.fill(mom2.eta());
+          h_eta_sublead_ei.fill(mom1.eta());
+        }
+        h_mdiff_ei.fill(evt(cepgen::Particle::Role::OutgoingBeam2).at(0).momentum().mass());
+      });
   }
   {  // inelastic-elastic
     pkin.set<int>("mode", 3);
     gen.runParameters().process().kinematics().setParameters(pkin);
     cs_ie = gen.computeXsection();
+    if (num_gen > 0)
+      gen.generate(num_gen, [&](const cepgen::Event& evt, size_t) {
+        const auto &mom1 = evt(cepgen::Particle::Role::CentralSystem).at(0).momentum(),
+                   &mom2 = evt(cepgen::Particle::Role::CentralSystem).at(1).momentum();
+        if (mom1.pt() > mom2.pt()) {
+          h_eta_lead_ie.fill(mom1.eta());
+          h_eta_sublead_ie.fill(mom2.eta());
+        } else {
+          h_eta_lead_ie.fill(mom2.eta());
+          h_eta_sublead_ie.fill(mom1.eta());
+        }
+        h_mdiff_ie.fill(evt(cepgen::Particle::Role::OutgoingBeam1).at(0).momentum().mass());
+      });
   }
   CG_TEST_VALUES(cs_ei, cs_ie, num_sigma, "el-inel == inel-el");
+
+  size_t ndf;
+  CG_TEST(h_eta_lead_ei.chi2test(h_eta_lead_ie, ndf) / ndf > 1., "leading lepton eta");
+  CG_TEST(h_eta_sublead_ei.chi2test(h_eta_sublead_ie, ndf) / ndf > 1., "subleading lepton eta");
+  CG_TEST(h_mdiff_ei.chi2test(h_mdiff_ie, ndf) / ndf < 1.1, "diffractive system mass");
+
+  if (!plotter.empty()) {
+    auto plt = cepgen::DrawerFactory::get().build(plotter);
+    plt->draw({&h_eta_lead_ie, &h_eta_lead_ei},
+              "leading_eta",
+              "leading lepton $\\eta$",
+              cepgen::utils::Drawer::Mode::nostack);
+    plt->draw({&h_eta_sublead_ie, &h_eta_sublead_ei},
+              "subleading_eta",
+              "subleading lepton $\\eta$",
+              cepgen::utils::Drawer::Mode::nostack);
+    plt->draw({&h_mdiff_ie, &h_mdiff_ei}, "mdiff", "diffractive system mass", cepgen::utils::Drawer::Mode::nostack);
+  }
+
   CG_TEST_SUMMARY;
 }