visible decay and A' kinematics scaling

app/scale.cxx

@@ -52,6 +52,7 @@ void usage() {
          "  -M,--ap-mass          : mass of dark photon in MeV\n"
          "  --muons               : pass to set lepton to muons (otherwise "
          "electrons)\n"
+         "  --scale-APrime        : pass to scale the APrime kinematics\n"
       << std::flush;
 }
 
@@ -71,6 +72,7 @@ int main(int argc, char* argv[]) try {
   std::string db_lib;
   double ap_mass{0.1};
   bool muons{false};
+  bool scale_APrime{false};
   for (int i_arg{1}; i_arg < argc; i_arg++) {
     std::string arg{argv[i_arg]};
     if (arg == "-h" or arg == "--help") {
@@ -108,6 +110,8 @@ int main(int argc, char* argv[]) try {
         return 1;
       }
       num_events = std::stoi(argv[++i_arg]);
+    } else if (arg == "--scale-APrime") {
+      scale_APrime = true;
     } else if (not arg.empty() and arg[0] == '-') {
       std::cerr << arg << " is not a recognized option" << std::endl;
       return 1;
@@ -129,32 +133,53 @@ int main(int argc, char* argv[]) try {
   }
 
   // the process accesses the A' mass from the G4 particle
-  G4APrime::Initialize(ap_mass / GeV);
+  G4APrime::Initialize(ap_mass);
+
   // create the model, this is where the LHE file is parsed
   //    into an in-memory library to sample and scale from
   g4db::G4DarkBreMModel db_model(
       db_lib, muons,
       0.0,  // threshold
       1.0,  // epsilon
-      g4db::G4DarkBreMModel::ScalingMethod::ForwardOnly);
+      g4db::G4DarkBreMModel::ScalingMethod::ForwardOnly,
+      g4db::G4DarkBreMModel::XsecMethod::Auto,
+      50.0, // max_R_for_full
+      622, // aprime_lhe_id
+      true, // load_library
+      scale_APrime);
   db_model.PrintInfo();
-  printf("   %-16s %f\n", "Lepton Mass [MeV]:", lepton_mass);
+  printf("   %-16s %f\n", "Lepton Mass [MeV]:", lepton_mass * GeV / MeV);
   printf("   %-16s %f\n", "A' Mass [MeV]:", ap_mass / MeV);
+
+  double lepton_mass_squared{lepton_mass * lepton_mass};
+  double ap_mass_squared{ap_mass * ap_mass};
 
   std::ofstream f{output_filename};
   if (not f.is_open()) {
     std::cerr << "Unable to open output file for writing." << std::endl;
     return -1;
   }
-  f << "recoil_energy,recoil_px,recoil_py,recoil_pz\n";
+  f << "target_Z,incident_energy,recoil_energy,recoil_px,recoil_py,recoil_pz,"
+    << "centerMomentum_energy,centerMomentum_px,centerMomentum_py,"
+    << "centerMomentum_pz\n";
 
   for (int i_event{0}; i_event < num_events; ++i_event) {
-    G4ThreeVector recoil =
+    std::pair<G4ThreeVector, G4ThreeVector> momenta = 
         db_model.scale(target_Z, incident_energy, lepton_mass);
-    double recoil_energy = sqrt(recoil.mag2() + lepton_mass * lepton_mass);
+    G4ThreeVector recoil = momenta.first;
+    double recoil_energy = sqrt(recoil.mag2() + lepton_mass_squared);
+    G4ThreeVector aprime = momenta.second;
+    double aprime_energy = sqrt(aprime.mag2() + ap_mass_squared);
 
-    f << recoil_energy << ',' << recoil.x() << ',' << recoil.y() << ','
-      << recoil.z() << '\n';
+    // convert to GeV to match MG output
+    f << target_Z << ',' << incident_energy << ',' 
+      << recoil_energy / GeV << ',' << recoil.x() / GeV << ',' 
+      << recoil.y() / GeV << ',' << recoil.z() / GeV << "," 
+      << (recoil_energy + aprime_energy) / GeV << "," 
+      << (recoil.x() + aprime.x()) / GeV << "," 
+      << (recoil.y() + aprime.y()) / GeV << "," 
+      << (recoil.z() + aprime.z() ) / GeV
+      << '\n';
   }
 
   f.close();

include/G4DarkBreM/G4APrime.h

@@ -27,34 +27,37 @@ class G4DecayTable;
  * ConstructParticle function of a physics constructor.
  */
 class G4APrime : public G4ParticleDefinition {
- private:
-  /** Reference to single particle definition of A' */
-  static G4APrime* theAPrime;
-
+ public:
   /**
-   * Constructor
+   * @enum DecayMode
    *
-   * Passes all parameters to the base class constructor
-   * to register this particle definition with Geant4.
+   * How to handle APrime decays.
    */
-  G4APrime(const G4String& Name, G4double mass, G4double width, G4double charge,
-           G4int iSpin, G4int iParity, G4int iConjugation, G4int iIsospin,
-           G4int iIsospin3, G4int gParity, const G4String& pType, G4int lepton,
-           G4int baryon, G4int encoding, G4bool stable, G4double lifetime,
-           G4DecayTable* decaytable)
-      : G4ParticleDefinition(Name, mass, width, charge, iSpin, iParity,
-                             iConjugation, iIsospin, iIsospin3, gParity, pType,
-                             lepton, baryon, encoding, stable, lifetime,
-                             decaytable) {}
+  enum class DecayMode {
+    /**
+     * No decay/stable -- this is the default.
+     */
+    NoDecay = 1,
 
-  /**
-   * Destructor
-   *
-   * Does nothing on purpose.
-   */
-  virtual ~G4APrime() {}
+    /**
+     * Flat decay -- the decay proper time will be randomly sampled
+     * from a uniform distribution whose maximum and minimum are set
+     * in the G4DarkBreMModel::GenerateChange() method based on parameters
+     * of that class.
+     *
+     * As with GeantDecay, the APrime will decay to e-/e+.
+     */
+     FlatDecay = 2,
+
+     /**
+      * Let Geant4 handle the decay. In this case, the lifetime
+      * must be included as a parameter of the constructor.
+      *
+      * e-/e+ is the only decay channel.
+      */
+      GeantDecay = 3
+  };
 
- public:
   /**
    * Accessor for APrime definition
    *
@@ -72,14 +75,49 @@ class G4APrime : public G4ParticleDefinition {
    *
    * @param[in] mass The mass of the APrime in MeV
    * @param[in] id The PDG ID number to use for the APrime particle
+   * @param[in] tau The proper lifetime 
+   *                (only used if decay_mode is set to GeantDecay)
+   * @param[in] decay_mode G4APrime::DecayMode on whether/how to decay the A'
    *
    * The default value for the PDG ID is set to 62 and has been arbitrarily
    * chosen from the range defined by 11(c) in the
    * [PDG ID numbering
    * scheme](https://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf), avoiding the
    * already-defined "one-of-a-kind" particles using 39, 41, and 42.
    */
-  static void Initialize(double mass, int id = 62);
+  static void Initialize(double mass, int id = 62, double tau = -1.0, 
+                         DecayMode decay_mode = DecayMode::NoDecay);
+
+  static DecayMode getDecayMode() { return decay_mode_; }
+
+  private:
+  /** Reference to single particle definition of A' */
+  static G4APrime* theAPrime;
+
+  /**
+   * Constructor
+   *
+   * Passes all parameters to the base class constructor
+   * to register this particle definition with Geant4.
+   */
+  G4APrime(const G4String& Name, G4double mass, G4double width, G4double charge,
+           G4int iSpin, G4int iParity, G4int iConjugation, G4int iIsospin,
+           G4int iIsospin3, G4int gParity, const G4String& pType, G4int lepton,
+           G4int baryon, G4int encoding, G4bool stable, G4double lifetime,
+           G4DecayTable* decaytable)
+      : G4ParticleDefinition(Name, mass, width, charge, iSpin, iParity,
+                             iConjugation, iIsospin, iIsospin3, gParity, pType,
+                             lepton, baryon, encoding, stable, lifetime,
+                             decaytable) {}
+
+  static DecayMode decay_mode_;
+
+  /**
+   * Destructor
+   *
+   * Does nothing on purpose.
+   */
+  virtual ~G4APrime() {}
 };
 
 #endif  // SIMCORE_DARKBREM_G4APRIME_H_

include/G4DarkBreM/G4DarkBreMModel.h

@@ -205,14 +205,21 @@ class G4DarkBreMModel : public PrototypeModel {
    * files.**
    * @param[in] load_library only used in cross section executable where it is
    * known that the library will not be used during program run
-   *
+   * @param[in] scale_APrime whether to scale the A' kinematics based on the
+   * library. If false, define A' momentum to conserve momentum without taking
+   * nuclear recoil into account.
+   * @param[in] dist_decay_min minimum flight distance [mm] at which to 
+   * decay the A' if G4APrime::DecayMode is set to FlatDecay
+   * @param[in] dist_decay_max maximum flight distance [mm] at which to 
+   * decay the A' if G4APrime::DecayMode is set to FlatDecay
    */
   G4DarkBreMModel(const std::string& library_path, bool muons,
                   double threshold = 0.0, double epsilon = 1.0,
                   ScalingMethod sm = ScalingMethod::ForwardOnly,
                   XsecMethod xm = XsecMethod::Auto,
                   double max_R_for_full = 50.0, int aprime_lhe_id = 622,
-                  bool load_library = true);
+                  bool load_library = true, bool scale_APrime = false,
+                  double dist_decay_min = 0.0, double dist_decay_max = 1.0);
 
   /**
    * Destructor
@@ -265,10 +272,11 @@ class G4DarkBreMModel : public PrototypeModel {
    * @param[in] target_Z atomic Z of target nucleus
    * @param[in] incident_energy incident total energy of the lepton [GeV]
    * @param[in] lepton_mass mass of incident lepton [GeV]
-   * @return G4ThreeVector representing the recoil lepton's outgoing momentum
+   * @return G4ThreeVectors representing the outgoing momenta of the recoil lepton 
+   *         and the A', respectively
    */
-  G4ThreeVector scale(double target_Z, double incident_energy,
-                      double lepton_mass);
+  std::pair<G4ThreeVector, G4ThreeVector> 
+	  scale(double target_Z, double incident_energy, double lepton_mass);
 
   /**
    * Simulates the emission of a dark photon + lepton
@@ -403,6 +411,30 @@ class G4DarkBreMModel : public PrototypeModel {
    */
   bool alwaysCreateNewLepton_{true};
 
+  /**
+   * whether to scale the outgoing A' momentum based on the MadGraph library.
+   * 
+   * The same scaling method is used as for the recoil lepton.
+   * The difference between the azimuthal angle of the A' and the recoil lepton
+   * will be preserved from MadGraph.
+   *
+   * If false, will define the 3-momentum of the A' to conserve 3-momentum
+   * with primary and recoil lepton, not taking into account the nuclear recoil
+   */  
+  bool scale_APrime_{false};
+
+  /**
+   * Minimum flight distance [mm] at which to decay the A' 
+   * if G4APrime::DecayMode is set to FlatDecay 
+   */
+  double dist_decay_min_{0.0};
+
+  /**
+   * Maximum flight distance [mm] at which to decay the A' 
+   * if G4APrime::DecayMode is set to FlatDecay 
+   */
+  double dist_decay_max_{1.0};
+
   /**
    * Storage of data from mad graph
    *

include/G4DarkBreM/ParseLibrary.h

@@ -19,10 +19,9 @@ namespace g4db {
  * Data frame to store necessary information from LHE files
  */
 struct OutgoingKinematics {
-  /// 4-momentum of lepton in center of momentum frame for electron-A'
-  /// system
+  /// 4-momentum of lepton in CoM frame
   CLHEP::HepLorentzVector lepton;
-  /// 4-vector pointing to center of momentum frame
+  /// 4-vector pointing to center of momentum of lepton and A'
   CLHEP::HepLorentzVector centerMomentum;
   /// energy of lepton before brem (used as key in mad graph data map)
   double E;
@@ -50,10 +49,10 @@ struct OutgoingKinematics {
  * 4. The x-component of the recoil momentum
  * 5. The y-component of the recoil momentum
  * 6. The z-component of the recoil momentum
- * 7. The total energy of the A'
- * 8. The x-component of the A' momentum
- * 9. The y-component of the A' momentum
- * 10. The z-component of the A' momentum
+ * 7. The total energy of the CoM frame
+ * 8. The x-component of the CoM
+ * 9. The y-component of the CoM
+ * 10. The z-component of the CoM
  *
  * ### LHE
  * The LHE files must have dark brem events in it where "dark brem event"
@@ -88,7 +87,7 @@ struct OutgoingKinematics {
  *
  * The `E` from the first line is used as the incident lepton energy.
  * The four-momentum from the middle line is the recoil lepton's four momentum,
- * and the four-momentum from the last line is used in conjuction with the
+ * and the four-momentum from the last line (A') is used in conjunction with the
  * recoil four-momentum to calculate the center of momentum vector.
  *
  * @param[in] path path to library to parse

src/G4DarkBreM/G4APrime.cxx

@@ -8,10 +8,14 @@
 
 #include "G4ParticleTable.hh"
 #include "G4PhysicalConstants.hh"
+#include "G4DecayTable.hh"
+#include "G4PhaseSpaceDecayChannel.hh"
 #include "globals.hh"
 
 G4APrime* G4APrime::theAPrime = 0;
 
+G4APrime::DecayMode G4APrime::decay_mode_ = G4APrime::DecayMode::NoDecay;
+
 G4APrime* G4APrime::APrime() {
   if (!theAPrime) {
     throw std::runtime_error(
@@ -22,11 +26,14 @@ G4APrime* G4APrime::APrime() {
   return theAPrime;
 }
 
-void G4APrime::Initialize(double mass, int id) {
+void G4APrime::Initialize(double mass, int id, double tau,
+                          G4APrime::DecayMode decay_mode) {
   if (theAPrime)
     throw std::runtime_error(
         "Attempting to initialize the APrime particle more than once.");
 
+  G4APrime::decay_mode_ = decay_mode;
+
   /**
    * Here are the properties of the formal Geant4 dark photon we define.
    *
@@ -46,16 +53,31 @@ void G4APrime::Initialize(double mass, int id) {
    * lepton number | 0
    * baryon number | 0
    * PDG ID encoding | **configured**
-   * is stable (no decay) | true
-   * lifetime | -1 (i.e. no decay)
-   * decay table | nullptr (i.e. no decay)
+   * is stable (no decay) | depends on DecayMode
+   * lifetime | depends on DecayMode
+   * decay table | depends on DecayMode
    */
+
   theAPrime = new G4APrime(
       "A^1" /* short name */, mass * MeV, 0. /* mass width */,
       0. /*electric charge */, 0 /* spin */, 0 /* parity */,
       0 /* conjugation */, 0 /* isospine */, 0 /* isospin3 */, 0 /* G parity */,
       "APrime" /* long name */, 0 /* lepton number */, 0 /* baryon number */,
       id, true /* stable? */, -1 /* lifetime */, nullptr /* decay table */
   );
+
+  if (decay_mode != G4APrime::DecayMode::NoDecay) {
+    G4DecayTable* table = new G4DecayTable();
+    G4VDecayChannel* mode = new G4PhaseSpaceDecayChannel("A^1", 1.0, 2, 
+                                                         "e-", "e+");
+    table->Insert(mode);
+
+    theAPrime->SetPDGStable(false);
+    theAPrime->SetPDGLifeTime(tau * second);
+    if (decay_mode == G4APrime::DecayMode::FlatDecay)
+      theAPrime->SetPDGLifeTime(0.0); // decay configured in G4DarkBreMModel
+    theAPrime->SetDecayTable(table);
+  }
+
   return;
 }