diff --git a/app/scale.cxx b/app/scale.cxx index 1b8557e..9ccb339 100644 --- a/app/scale.cxx +++ b/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 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(); diff --git a/docs/configuration.md b/docs/configuration.md index 92156e1..f22d754 100644 --- a/docs/configuration.md +++ b/docs/configuration.md @@ -25,16 +25,21 @@ There are other parameters that tune the functionality of the simulation, but the parameters listed here are directly related to how the physics of the simulation is done. -Parameter | Required | Location | Description ----------------------|----------|-----------------------|------------- -A' Mass | Yes | G4APrime | Mass of dark photon in MeV -MG/ME Data | Yes | g4db::G4DarkBreMModel | Library of dark brem events -Use Muons | Yes | g4db::G4DarkBreMModel | Choice to dark brem off muons or electrons -Minimum DB Threshold | No | g4db::G4DarkBreMModel | Minimum energy in GeV lepton needs to have to DB (default twice A' mass) -Mixing Strength | No | g4db::G4DarkBreMModel | Mixing strength between dark photon and SM photon (default 1.0) -Scaling Method | No | g4db::G4DarkBreMModel | Choice of how to scale the outgoing kinematics from the library (default ForwardOnly) -Xsec Method | No | g4db::G4DarkBreMModel | Choice of WW approximation variant to calculate the total cross section. Options include Full, Improved, HyperImproved, and Auto (default is Auto where Full is used if the A' to lepton mass ratio is below the max and Improved is used otherwise) -Maximum mass ratio | No | g4db::G4DarkBreMModel | If using Auto Xsec Method, the maximum ratio of A' to lepton mass to transition from Full WW to Improved WW (default 50.0) -A' LHE ID | No | g4db::G4DarkBreMModel | If MG/ME Data is in the form of LHE files, the PDG ID of the A' inside of those files (default 622) -Only One per Event | No | G4DarkBremsstrahlung | Limit the dark brem process to only one interaction per event (default false) -Global Bias | No | G4DarkBremsstrahlung | Bias the dark brem process everywhere by this factor (default 1.0) +Parameter | Required | Location | Description +----------------------|----------|-----------------------|------------- +A' Mass | Yes | G4APrime | Mass of dark photon in MeV +A' tau | No | G4APrime | Proper lifetime of dark photon in seconds, only relevant if decay mode is set to GeantDecay (default -1, meaning no decay) +Decay Mode | No | G4APrime | Choice of whether to allow the A' to decay and the decay distribution. The only decay channel is e+/e- (default NoDecay) +MG/ME Data | Yes | g4db::G4DarkBreMModel | Library of dark brem events +Use Muons | Yes | g4db::G4DarkBreMModel | Choice to dark brem off muons or electrons +Minimum DB Threshold | No | g4db::G4DarkBreMModel | Minimum energy in GeV lepton needs to have to DB (default twice A' mass) +Mixing Strength | No | g4db::G4DarkBreMModel | Mixing strength between dark photon and SM photon (default 1.0) +Scaling Method | No | g4db::G4DarkBreMModel | Choice of how to scale the outgoing kinematics from the library (default ForwardOnly) +Xsec Method | No | g4db::G4DarkBreMModel | Choice of WW approximation variant to calculate the total cross section. Options include Full, Improved, HyperImproved, and Auto (default is Auto where Full is used if the A' to lepton mass ratio is below the max and Improved is used otherwise) +Maximum mass ratio | No | g4db::G4DarkBreMModel | If using Auto Xsec Method, the maximum ratio of A' to lepton mass to transition from Full WW to Improved WW (default 50.0) +A' LHE ID | No | g4db::G4DarkBreMModel | If MG/ME Data is in the form of LHE files, the PDG ID of the A' inside of those files (default 622) +Scale A' | No | g4db::G4DarkBreMModel | Whether to scale the outgoing dark photon energy along with the recoil lepton (default false, meaning the dark photon will conserve total momentum while ignoring the nuclear recoil) +Min A' decay distance | No | g4db::G4DarkBreMModel | When A' decay mode is set to FlatDecay, this is the minimum lab-frame distance in mm the A' will travel before decaying (default 0) +Max A' decay distance | No | g4db::G4DarkBreMModel | When A' decay mode is set to FlatDecay, this is the maximum lab-frame distance in mm the A' will travel before decaying (default 1) +Only One per Event | No | G4DarkBremsstrahlung | Limit the dark brem process to only one interaction per event (default false) +Global Bias | No | G4DarkBremsstrahlung | Bias the dark brem process everywhere by this factor (default 1.0) diff --git a/include/G4DarkBreM/G4APrime.h b/include/G4DarkBreM/G4APrime.h index 840cf44..cccaed9 100644 --- a/include/G4DarkBreM/G4APrime.h +++ b/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,6 +75,9 @@ 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 @@ -79,7 +85,39 @@ class G4APrime : public G4ParticleDefinition { * 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_ diff --git a/include/G4DarkBreM/G4DarkBreMModel.h b/include/G4DarkBreM/G4DarkBreMModel.h index ef53a49..89339b8 100644 --- a/include/G4DarkBreM/G4DarkBreMModel.h +++ b/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 + scale(double target_Z, double incident_energy, double lepton_mass); /** * Simulates the emission of a dark photon + lepton @@ -403,6 +411,40 @@ 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 + * is designed to be preserved from MadGraph. + * The preservation of this angle and the overall scaling of the A' has only been validated for + * electrons with incident energies ranging from 1GeV to 10GeV and dark photon masses + * ranging from 1MeV to 100MeV. While the scaling is expected to be well behaved for muons + * and at other energy scales, users are encouraged to double-check this behavior in their + * situation and report issues to the repository. https://github.com/LDMX-Software/G4DarkBreM + * + * 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. + * This is the default because then the user is able to "reconstruct" the true incident lepton's + * momentum using the recoil lepton and produced dark photon's three-momenta. + * Scaling the A' momentum is only advisable if the A' momentum _itself_ is important + * to the search (for example, when the A' decays visibly and the decay is observed + * in which case the A' momentum deciding where the decay occurs is very important). + */ + 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 * diff --git a/include/G4DarkBreM/ParseLibrary.h b/include/G4DarkBreM/ParseLibrary.h index cec30d4..2b74add 100644 --- a/include/G4DarkBreM/ParseLibrary.h +++ b/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 diff --git a/src/G4DarkBreM/G4APrime.cxx b/src/G4DarkBreM/G4APrime.cxx index 7460853..2b79e15 100644 --- a/src/G4DarkBreM/G4APrime.cxx +++ b/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,19 @@ 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."); + if (decay_mode == G4APrime::DecayMode::GeantDecay && tau < 0.0) + throw std::runtime_error( + "Invalid configuration: DecayMode set to GeantDecay but tau is negative." + ); + + G4APrime::decay_mode_ = decay_mode; + /** * Here are the properties of the formal Geant4 dark photon we define. * @@ -46,10 +58,11 @@ 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 */, @@ -57,5 +70,19 @@ void G4APrime::Initialize(double mass, int id) { "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; } diff --git a/src/G4DarkBreM/G4DarkBreMModel.cxx b/src/G4DarkBreM/G4DarkBreMModel.cxx index 2a83d45..e910c8e 100644 --- a/src/G4DarkBreM/G4DarkBreMModel.cxx +++ b/src/G4DarkBreM/G4DarkBreMModel.cxx @@ -17,6 +17,7 @@ #include "G4RunManager.hh" //for VerboseLevel #include "G4SystemOfUnits.hh" #include "Randomize.hh" +#include "G4PhaseSpaceDecayChannel.hh" // Boost #include @@ -145,7 +146,9 @@ G4DarkBreMModel::G4DarkBreMModel(const std::string &library_path, bool muons, double threshold, double epsilon, ScalingMethod scaling_method, XsecMethod xsec_method, double max_R_for_full, - int aprime_lhe_id, bool load_library) + int aprime_lhe_id, bool load_library, + bool scale_APrime, double dist_decay_min, + double dist_decay_max) : PrototypeModel(muons), maxIterations_{10000}, threshold_{std::max(threshold, @@ -154,7 +157,10 @@ G4DarkBreMModel::G4DarkBreMModel(const std::string &library_path, bool muons, aprime_lhe_id_{aprime_lhe_id}, scaling_method_{scaling_method}, xsec_method_{xsec_method}, - library_path_{library_path} { + library_path_{library_path}, + scale_APrime_{scale_APrime}, + dist_decay_min_{dist_decay_min}, + dist_decay_max_{dist_decay_max} { if (xsec_method_ == XsecMethod::Auto) { static const double MA = G4APrime::APrime()->GetPDGMass() / GeV; const double lepton_mass{(muons_ ? G4MuonMinus::MuonMinus()->GetPDGMass() @@ -395,30 +401,55 @@ G4double G4DarkBreMModel::ComputeCrossSectionPerAtom(G4double lepton_ke, return cross; } -G4ThreeVector G4DarkBreMModel::scale(double target_Z, double incident_energy, - double lepton_mass) { +std::pair G4DarkBreMModel::scale( + double target_Z, double incident_energy, double lepton_mass) { // mass A' in GeV static const double MA = G4APrime::APrime()->GetPDGMass() / CLHEP::GeV; + OutgoingKinematics data = sample(target_Z, incident_energy); - double EAcc = - (data.lepton.e() - lepton_mass) * - ((incident_energy - lepton_mass - MA) / (data.E - lepton_mass - MA)) + - lepton_mass; + double energy_factor = ((incident_energy - lepton_mass - MA) / + (data.E - lepton_mass - MA)); + double EAcc = (data.lepton.e() - lepton_mass) * energy_factor + lepton_mass; double Pt = data.lepton.perp(); double P = sqrt(EAcc * EAcc - lepton_mass * lepton_mass); + + double APrimeE{0.0}; + double APrimePt{0.0}; + double APrimeP{0.0}; + + CLHEP::HepLorentzVector ap(data.centerMomentum.px() - data.lepton.px(), + data.centerMomentum.py() - data.lepton.py(), + data.centerMomentum.pz() - data.lepton.pz(), + data.centerMomentum.e() - data.lepton.e()); + if (scaling_method_ == ScalingMethod::ForwardOnly) { + if (scale_APrime_) { + APrimeE = (ap.e() - MA) * energy_factor + MA; + APrimePt = ap.perp(); + APrimeP = sqrt(APrimeE * APrimeE - MA * MA); + } + unsigned int i = 0; while (Pt * Pt + lepton_mass * lepton_mass > EAcc * EAcc) { // Skip events until the transverse energy is less than the total energy. i++; data = sample(target_Z, incident_energy); - EAcc = (data.lepton.e() - lepton_mass) * - ((incident_energy - lepton_mass - MA) / - (data.E - lepton_mass - MA)) + - lepton_mass; + energy_factor = ((incident_energy - lepton_mass - MA) / + (data.E - lepton_mass - MA)); + EAcc = (data.lepton.e() - lepton_mass) * energy_factor + lepton_mass; Pt = data.lepton.perp(); P = sqrt(EAcc * EAcc - lepton_mass * lepton_mass); + if (scale_APrime_) { + ap.setPx(data.centerMomentum.px() - data.lepton.px()); + ap.setPy(data.centerMomentum.py() - data.lepton.py()); + ap.setPz(data.centerMomentum.pz() - data.lepton.pz()); + ap.setE(data.centerMomentum.e() - data.lepton.e()); + APrimeE = (ap.e() - MA) * energy_factor + MA; + APrimePt = ap.perp(); + APrimeP = sqrt(APrimeE * APrimeE - MA * MA); + } + if (i > maxIterations_) { std::cerr << "Could not produce a realistic vertex with library energy " @@ -432,24 +463,40 @@ G4ThreeVector G4DarkBreMModel::scale(double target_Z, double incident_energy, } else if (scaling_method_ == ScalingMethod::CMScaling) { CLHEP::HepLorentzVector el(data.lepton.px(), data.lepton.py(), data.lepton.pz(), data.lepton.e()); + double ediff = data.E - incident_energy; CLHEP::HepLorentzVector newcm( data.centerMomentum.px(), data.centerMomentum.py(), data.centerMomentum.pz() - ediff, data.centerMomentum.e() - ediff); el.boost(-1. * data.centerMomentum.boostVector()); el.boost(newcm.boostVector()); - double newE = (data.lepton.e() - lepton_mass) * - ((incident_energy - lepton_mass - MA) / - (data.E - lepton_mass - MA)) + - lepton_mass; + + double newE = (data.lepton.e() - lepton_mass) * energy_factor + lepton_mass; el.setE(newE); EAcc = el.e(); Pt = el.perp(); P = el.vect().mag(); + + if (scale_APrime_) { + double oldE = ap.e(); + ap.boost(-1. * data.centerMomentum.boostVector()); + ap.boost(newcm.boostVector()); + + newE = (oldE - MA) * energy_factor + MA; + ap.setE(newE); + APrimePt = ap.perp(); + APrimeP = ap.vect().mag(); + } } else if (scaling_method_ == ScalingMethod::Undefined) { EAcc = data.lepton.e(); P = sqrt(EAcc * EAcc - lepton_mass * lepton_mass); Pt = data.lepton.perp(); + + if (scale_APrime_) { + APrimeE = ap.e(); + APrimeP = sqrt(APrimeE * APrimeE - MA * MA); + APrimePt = ap.perp(); + } } else { throw std::runtime_error( "Unrecognized ScalingMethod, should be set by using the enum class."); @@ -463,7 +510,26 @@ G4ThreeVector G4DarkBreMModel::scale(double target_Z, double incident_energy, recoil.set(std::sin(ThetaAcc) * std::cos(PhiAcc), std::sin(ThetaAcc) * std::sin(PhiAcc), std::cos(ThetaAcc)); recoil.setMag(recoilMag); - return recoil; + + // outgoing A' momentum + G4ThreeVector aprime; + if (scale_APrime_) { + // make the azimuthal (phi) angle difference the same as in the MG data + G4double aPrimePhi = PhiAcc + ap.phi() - data.lepton.phi(); + + G4double aPrimeMag = sqrt(APrimeE * APrimeE - MA * MA) * GeV; + double aPrimeTheta = std::asin(APrimePt / APrimeP); + aprime.set(std::sin(aPrimeTheta) * std::cos(aPrimePhi), + std::sin(aPrimeTheta) * std::sin(aPrimePhi), + std::cos(aPrimeTheta)); + aprime.setMag(aPrimeMag); + } else { + // define A' momentum to be 3-momentum conserving with the recoil + double incident_momentum_mag = sqrt(incident_energy*incident_energy + - lepton_mass*lepton_mass) * GeV; + aprime = G4ThreeVector(0, 0, incident_momentum_mag) - recoil; + } + return std::make_pair(recoil, aprime); } void G4DarkBreMModel::GenerateChange(G4ParticleChange &particleChange, @@ -476,16 +542,28 @@ void G4DarkBreMModel::GenerateChange(G4ParticleChange &particleChange, G4double incidentEnergy = step.GetPostStepPoint()->GetTotalEnergy() / CLHEP::GeV; - G4ThreeVector recoilMomentum = scale(element.GetZ(), incidentEnergy, Ml); - recoilMomentum.rotateUz(track.GetMomentumDirection()); + std::pair outgoingMomenta + = scale(element.GetZ(), incidentEnergy, Ml); + outgoingMomenta.first.rotateUz(track.GetMomentumDirection()); + outgoingMomenta.second.rotateUz(track.GetMomentumDirection()); + G4ThreeVector recoilMomentum = outgoingMomenta.first; + G4ThreeVector darkPhotonMomentum = outgoingMomenta.second; - // create g4dynamicparticle object for the dark photon. - // define its 3-momentum so we conserve 3-momentum with primary and recoil - // lepton NOTE: does _not_ take nucleus recoil into account - G4ThreeVector darkPhotonMomentum = track.GetMomentum() - recoilMomentum; G4DynamicParticle *dphoton = new G4DynamicParticle(G4APrime::APrime(), darkPhotonMomentum); + if (G4APrime::getDecayMode() == G4APrime::DecayMode::FlatDecay) { + double dist_decay = G4UniformRand() * (dist_decay_max_ - dist_decay_min_) + + dist_decay_min_; + CLHEP::HepLorentzVector p4 = dphoton->Get4Momentum(); + double tau_decay = dist_decay / p4.beta() / p4.gamma() / c_light; + dphoton->SetPreAssignedDecayProperTime(tau_decay); + + dphoton->SetPreAssignedDecayProducts( + (new G4PhaseSpaceDecayChannel("A^1", 1.0, 2, "e-", "e+"))-> + DecayIt(G4APrime::APrime()->G4APrime::APrime()->GetPDGMass())); + } + // stop tracking and create new secondary instead of primary if (alwaysCreateNewLepton_) { /* diff --git a/src/G4DarkBreM/ParseLibrary.cxx b/src/G4DarkBreM/ParseLibrary.cxx index cf72c36..736a910 100644 --- a/src/G4DarkBreM/ParseLibrary.cxx +++ b/src/G4DarkBreM/ParseLibrary.cxx @@ -70,7 +70,7 @@ namespace parse { * * 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 conjuction with the * recoil four-momentum to calculate the center of momentum vector. * * @param[in] reader input stream reading the file @@ -117,13 +117,10 @@ void lhe( a_py >> a_pz >> a_E >> a_M; if (ptype == aprime_lhe_id and state == 1) { OutgoingKinematics evnt; - double cmpx = a_px + e_px; - double cmpy = a_py + e_py; - double cmpz = a_pz + e_pz; - double cmE = a_E + e_E; evnt.lepton = CLHEP::HepLorentzVector(e_px, e_py, e_pz, e_E); evnt.centerMomentum = - CLHEP::HepLorentzVector(cmpx, cmpy, cmpz, cmE); + CLHEP::HepLorentzVector(e_px + a_px, e_py + a_py, e_pz + a_pz, + e_E + a_E); evnt.E = incident_energy; if (target_Z < 0) { throw std::runtime_error( @@ -152,10 +149,10 @@ void lhe( * 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 + * 8. The x-component of the CoM + * 9. The y-component of the CoM + * 10. The z-component of the CoM * * @note If developing this function, make sure to update dumpLibrary * so that they can be used in conjuction. @@ -182,14 +179,14 @@ void csv( if (not lss and cell.empty()) vals.push_back(-9999); if (vals.size() != 10) { throw std::runtime_error( - "Malformed row in CSV file: not exactly 9 columns"); + "Malformed row in CSV file: not exactly 10 columns"); } OutgoingKinematics ok; int target_Z = vals[0]; // implicit drop of any decimal points ok.E = vals[1]; ok.lepton = CLHEP::HepLorentzVector(vals[3], vals[4], vals[5], vals[2]); - ok.centerMomentum = - CLHEP::HepLorentzVector(vals[7], vals[8], vals[9], vals[6]); + ok.centerMomentum + = CLHEP::HepLorentzVector(vals[7], vals[8], vals[9], vals[6]); lib[target_Z][ok.E].push_back(ok); } } @@ -198,7 +195,7 @@ void csv( void parseLibrary( const std::string& path, int aprime_lhe_id, - std::map>>& lib) { + std::map>>& lib) { if (hasEnding(path, ".csv") or hasEnding(path, ".csv.gz") or hasEnding(path, ".lhe") or hasEnding(path, ".lhe.gz")) { /**