LPAIR+others: fix of single-dissociation symmetrisation

CepGen/Process/FactorisedProcess.cpp

@@ -33,13 +33,15 @@ namespace cepgen {
           psgen_(PhaseSpaceGeneratorFactory::get().build(
               steer<ParametersList>("kinematicsGenerator")
                   .set("ids", std::vector<int>(central.begin(), central.end())))),
+          symmetrise_(steer<bool>("symmetrise")),
           store_alphas_(steer<bool>("storeAlphas")) {
       event().map()[Particle::CentralSystem].resize(central.size());
     }
 
     FactorisedProcess::FactorisedProcess(const FactorisedProcess& proc)
         : Process(proc),
           psgen_(PhaseSpaceGeneratorFactory::get().build(proc.psgen_->parameters())),
+          symmetrise_(proc.symmetrise_),
           store_alphas_(proc.store_alphas_) {}
 
     void FactorisedProcess::addEventContent() {
@@ -73,13 +75,18 @@ namespace cepgen {
         defineVariable(mX2(), Mapping::square, kinematics().cuts().remnants.mx, "Positive-z beam remnant squared mass");
       if (!kinematics().incomingBeams().negative().elastic())
         defineVariable(mY2(), Mapping::square, kinematics().cuts().remnants.mx, "Negative-z beam remnant squared mass");
+      // symmetrisation of the phase space: factor 2 in Jacobian for single-dissociation
+      kin_prefactor_ = symmetrise_ && (kinematics().incomingBeams().mode() == mode::Kinematics::ElasticInelastic ||
+                                       kinematics().incomingBeams().mode() == mode::Kinematics::InelasticElastic)
+                           ? 2.
+                           : 1.;
     }
 
     double FactorisedProcess::computeWeight() {
       if (!psgen_->generate())
         return 0.;
       if (const auto cent_weight = computeFactorisedMatrixElement(); utils::positive(cent_weight))
-        return cent_weight * psgen_->weight();
+        return cent_weight * psgen_->weight() * kin_prefactor_;
       return 0.;
     }
 
@@ -96,6 +103,14 @@ namespace cepgen {
       part1.setMomentum(pA() - pX(), true);
       part2.setMomentum(pB() - pY(), true);
 
+      if (symmetrise_ && rnd_gen_->uniformInt(0, 1) == 1) {  // symmetrise the el-in and in-el cases
+        std::swap(pX(), pY());
+        std::swap(q1(), q2());
+        std::swap(pc(0), pc(1));
+        for (auto* mom : {&q1(), &q2(), &pX(), &pY(), &pc(0), &pc(1)})
+          mom->mirrorZ();
+      }
+
       // add couplings to metadata
       if (store_alphas_) {
         const auto two_part_mass = (part1.momentum() + part2.momentum()).mass();
@@ -122,6 +137,7 @@ namespace cepgen {
                            2, PartonFluxFactory::get().describeParameters("BudnevElastic").parameters()))
                    .setDescription("Parton fluxes modelling"));
       desc.add("kinematicsGenerator", PhaseSpaceGeneratorFactory::get().describeParameters("kt2to4"));
+      desc.add<bool>("symmetrise", false).setDescription("Symmetrise along z the central system?");
       desc.add<bool>("storeAlphas", false)
           .setDescription("store the electromagnetic and strong coupling constants to the event content?");
       return desc;

CepGen/Process/FactorisedProcess.h

@@ -56,7 +56,11 @@ namespace cepgen {
 
       /// Kinematic variables generator for the phase space coverage
       const std::unique_ptr<PhaseSpaceGenerator> psgen_;
+      const bool symmetrise_;
       const bool store_alphas_;
+
+    private:
+      double kin_prefactor_{1.};
     };
   }  // namespace proc
 }  // namespace cepgen

CepGenProcesses/LPAIR.cpp

@@ -120,6 +120,12 @@ class LPAIR final : public cepgen::proc::Process {
       defineVariable(mX2(), Mapping::power_law, mx_range(mA()), "MX2");
     if (beams_mode_ == mode::Kinematics::InelasticInelastic)  // second outgoing beam particle or remnant mass
       defineVariable(mY2(), Mapping::power_law, mx_range(mB()), "MY2");
+    if (symmetrise_ &&
+        (beams_mode_ == mode::Kinematics::InelasticElastic || beams_mode_ == mode::Kinematics::ElasticInelastic))
+      CG_INFO("LPAIR:prepareKinematics")
+          << "Single dissociation kinematics mode was enabled with symmetrisation of the outgoing system.\n\t"
+             "The generator-level cross section will be doubled, and beam particles, incoming partons, and central "
+             "system will be mirrored in z.";
   }
 
   void fillKinematics() override {
@@ -136,16 +142,14 @@ class LPAIR final : public cepgen::proc::Process {
     // randomly rotate all particles
     const short rany = rnd_gen_->uniformInt(0, 1) == 1 ? 1 : -1;
     const double ranphi = rnd_gen_->uniform(0., 2. * M_PI);
-    const bool mirror = rnd_gen_->uniformInt(0, 1) == 1;
-    for (auto* mom : {&q1(), &q2(), &pc(0), &pc(1), &pX(), &pY()}) {
+    for (auto* mom : {&q1(), &q2(), &pX(), &pY(), &pc(0), &pc(1)})
       mom->rotatePhi(ranphi, rany);
-      if (symmetrise_ && mirror)
-        mom->mirrorZ();
-    }
-    if (beams_mode_ == mode::Kinematics::ElasticInelastic) {  // mirror X/Y and dilepton systems if needed
+    if ((symmetrise_ && rnd_gen_->uniformInt(0, 1) == 1) ||
+        beams_mode_ == mode::Kinematics::ElasticInelastic) {  // mirror X/Y and dilepton systems if needed
       std::swap(pX(), pY());
+      std::swap(q1(), q2());
       std::swap(pc(0), pc(1));
-      for (auto* mom : {&pX(), &pY(), &pc(0), &pc(1)})
+      for (auto* mom : {&q1(), &q2(), &pX(), &pY(), &pc(0), &pc(1)})
         mom->mirrorZ();
     }
     // first outgoing beam

test/physics/process_sd_symmetrisation.cc

@@ -0,0 +1,75 @@
+/*
+ *  CepGen: a central exclusive processes event generator
+ *  Copyright (C) 2023-2024  Laurent Forthomme
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#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/Test.h"
+#include "CepGen/Utils/Timer.h"
+
+using namespace std;
+
+int main(int argc, char* argv[]) {
+  double num_sigma;
+  string str_fun, proc_name, integrator;
+
+  cepgen::ArgumentsParser(argc, argv)
+      .addOptionalArgument("process,p", "process to compute", &proc_name, "lpair")
+      .addOptionalArgument("num-sigma,n", "max. number of std.dev.", &num_sigma, 3.)
+      .addOptionalArgument("str-fun,s", "struct.functions modelling", &str_fun, "SuriYennie")
+      .addOptionalArgument("integrator,i", "type of integrator used", &integrator, "Vegas")
+      .parse();
+
+  cepgen::utils::Timer tmr;
+  cepgen::Generator gen;
+  gen.runParameters().integrator().setName(integrator);
+
+  cepgen::utils::AbortHandler ah;
+
+  auto pkin =
+      cepgen::ParametersList()
+          .set<double>("sqrtS", 13.e3)
+          .set<cepgen::ParametersList>(
+              "structureFunctions", cepgen::StructureFunctionsFactory::get().describeParameters(str_fun).parameters())
+          .set<double>("ptmin", 5.)
+          .set<int>("mode", 2)  // elastic-inelastic
+          .set<cepgen::Limits>("eta", {-2.5, 2.5})
+          .set<cepgen::Limits>("mx", {1.07, 1000.});
+
+  cepgen::Value cs_ei_no_symm, cs_ei_symm;
+  {
+    gen.runParameters().setProcess(
+        cepgen::ProcessFactory::get().build(proc_name, cepgen::ParametersList().set<int>("pair", 13)));
+    gen.runParameters().process().kinematics().setParameters(pkin);
+    cs_ei_no_symm = gen.computeXsection();
+  }
+  {  // inelastic-elastic
+    gen.runParameters().setProcess(cepgen::ProcessFactory::get().build(
+        proc_name, cepgen::ParametersList().set<int>("pair", 13).set<bool>("symmetrise", true)));
+    gen.runParameters().process().kinematics().setParameters(pkin);
+    cs_ei_symm = gen.computeXsection();
+  }
+  CG_TEST_VALUES(cs_ei_no_symm * 2., cs_ei_symm, num_sigma, "symmetrised SD == 2*non-symmetrised SD");
+
+  CG_TEST_SUMMARY;
+}