Public Member Functions | Protected Member Functions | Private Types | Static Private Member Functions | Private Attributes | List of all members
genie::DarkSectorDecayer Class Reference

Dark Sector decayer module. More...

#include <DarkSectorDecayer.h>

Inheritance diagram for genie::DarkSectorDecayer:
genie::EventRecordVisitorI genie::Algorithm

Public Member Functions

 DarkSectorDecayer ()
 
 DarkSectorDecayer (string config)
 
virtual ~DarkSectorDecayer ()
 
virtual void Configure (const Registry &config)
 
virtual void Configure (string config)
 
void ProcessEventRecord (GHepRecord *event) const
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Protected Member Functions

virtual void LoadConfig (void)
 
bool ToBeDecayed (const GHepParticle &p) const
 
std::vector< GHepParticleDecay (const GHepParticle &mother, const std::vector< int > &pdg_daughters) const
 
int SelectDecayChannel (const std::vector< DecayChannel > &dcs, double total_amplitude) const
 
std::vector< DecayChannelDarkMediatorDecayChannels (void) const
 
std::vector< DecayChannelDarkNeutrinoDecayChannels (int mother_pdg) const
 
void SetSpaceTime (std::vector< GHepParticle > &pp, const GHepParticle &mother, double total_amplitude) const
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 

Private Types

using DecayChannel = std::pair< std::vector< int >, double >
 

Static Private Member Functions

static string ParticleGunKineAsString (const TLorentzVector &vec4)
 

Private Attributes

TGenPhaseSpace fPhaseSpaceGenerator
 
double fEps2
 
std::array< double, 4 > fMixing2s
 
double fAlpha_D
 
double fDNuMass
 
double fDNuMass2
 
double fDMediatorMass
 
double fDMediatorMass2
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Dark Sector decayer module.

A simple decay simulation... .... Is a concerete implementation of the EventRecordVisitorI interface.

Author
Iker de Icaza <i.de-icaza-astiz sussex.ac.uk> University of Sussex

Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

July XX, 2020

Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 37 of file DarkSectorDecayer.h.

Member Typedef Documentation

using genie::DarkSectorDecayer::DecayChannel = std::pair<std::vector<int>, double>
private

Definition at line 39 of file DarkSectorDecayer.h.

Constructor & Destructor Documentation

DarkSectorDecayer::DarkSectorDecayer ( )

Definition at line 43 of file DarkSectorDecayer.cxx.

43  :
44  EventRecordVisitorI("genie::DarkSectorDecayer")
45 {
46 
47 }
DarkSectorDecayer::DarkSectorDecayer ( string  config)

Definition at line 49 of file DarkSectorDecayer.cxx.

49  :
50  EventRecordVisitorI("genie::DarkSectorDecayer", config)
51 {
52 
53 }
static Config * config
Definition: config.cpp:1054
DarkSectorDecayer::~DarkSectorDecayer ( )
virtual

Definition at line 55 of file DarkSectorDecayer.cxx.

56 {
57 
58 }

Member Function Documentation

void DarkSectorDecayer::Configure ( const Registry config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 342 of file DarkSectorDecayer.cxx.

343 {
344  Algorithm::Configure(config);
345  this->LoadConfig();
346 }
virtual void LoadConfig(void)
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void DarkSectorDecayer::Configure ( string  config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 348 of file DarkSectorDecayer.cxx.

349 {
351  this->LoadConfig();
352 }
virtual void LoadConfig(void)
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
std::vector< DarkSectorDecayer::DecayChannel > DarkSectorDecayer::DarkMediatorDecayChannels ( void  ) const
protected

Definition at line 216 of file DarkSectorDecayer.cxx.

217 {
218  // eq (4) and (5) and maybe some other higher order variations
219 
220  static constexpr std::array<int, 3> neutrinos = {kPdgNuE, kPdgNuMu, kPdgNuTau};
221  static constexpr std::array<int, 3> antineutrinos = {kPdgAntiNuE, kPdgAntiNuMu, kPdgAntiNuTau};
222  std::vector<DarkSectorDecayer::DecayChannel> dcs;
223 
224  for(size_t i=0; i<neutrinos.size(); ++i){
225  for(size_t j=0; j<antineutrinos.size(); ++j){// for antineutrinos
226  const double decay_width = fAlpha_D/3. * fMixing2s[i]*fMixing2s[j] * fDMediatorMass;
227  dcs.push_back(DecayChannel{{neutrinos[i], antineutrinos[j]}, decay_width});
228  }
229  }
230 
231  static const double electron_threshold = 2.*PDGLibrary::Instance()->Find(kPdgElectron)->Mass() ;
232  if(fDMediatorMass > electron_threshold ){
233  double ratio = electron_threshold / fDMediatorMass ;
234  double phase_space_correction = sqrt(1. - ratio*ratio ) ;
235  const double decay_width = kAem*fEps2/3. * fDMediatorMass * phase_space_correction ;
236  dcs.push_back(DecayChannel{{kPdgElectron, kPdgPositron}, decay_width});
237  }
238 
239  static const double muon_threshold = 2.*PDGLibrary::Instance()->Find(kPdgMuon)->Mass() ;
240  if(fDMediatorMass > muon_threshold ){
241  double ratio = muon_threshold / fDMediatorMass ;
242  double phase_space_correction = sqrt(1. - ratio*ratio ) ;
243  const double decay_width = kAem*fEps2/3. * fDMediatorMass * phase_space_correction ;
244  dcs.push_back(DecayChannel{{kPdgMuon, kPdgAntiMuon}, decay_width});
245  }
246  return dcs;
247 }
const int kPdgNuE
Definition: PDGCodes.h:28
const int kPdgAntiNuE
Definition: PDGCodes.h:29
std::array< double, 4 > fMixing2s
const int kPdgNuMu
Definition: PDGCodes.h:30
const int kPdgAntiMuon
Definition: PDGCodes.h:38
const int kPdgElectron
Definition: PDGCodes.h:35
std::pair< std::vector< int >, double > DecayChannel
static const double kAem
Definition: Constants.h:56
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
const int kPdgNuTau
Definition: PDGCodes.h:32
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
const int kPdgMuon
Definition: PDGCodes.h:37
const int kPdgPositron
Definition: PDGCodes.h:36
std::vector< DarkSectorDecayer::DecayChannel > DarkSectorDecayer::DarkNeutrinoDecayChannels ( int  mother_pdg) const
protected

Definition at line 249 of file DarkSectorDecayer.cxx.

250 {
251  // eq (3) and higher order variations
252 
253  static constexpr std::array<int, 3> neutrinos = {kPdgNuE, kPdgNuMu, kPdgNuTau};
254  static constexpr std::array<int, 3> antineutrinos = {kPdgAntiNuE, kPdgAntiNuMu, kPdgAntiNuTau};
255  std::vector<DarkSectorDecayer::DecayChannel> dcs;
256 
257  if(fDNuMass > fDMediatorMass){
258  for(size_t i=0; i<neutrinos.size(); ++i){
259  const double mass2ratio = fDMediatorMass2/fDNuMass2;
260  const double p0 = 0.5 * fAlpha_D * fMixing2s[3] * fMixing2s[i];
261  const double p1 = fDNuMass*fDNuMass2/fDMediatorMass2;
262  const double p2 = 1 - mass2ratio;
263  const double p3 = 1 + mass2ratio - 2*mass2ratio*mass2ratio;
264  const double decay_width = p0 * p1 * p2 * p3;
265 
266  if(mother_pdg == kPdgDarkNeutrino){
267  dcs.push_back(DecayChannel{{neutrinos[i], kPdgDNuMediator}, decay_width});
268  }
269  else if(mother_pdg == kPdgAntiDarkNeutrino){
270  dcs.push_back(DecayChannel{{antineutrinos[i], kPdgDNuMediator}, decay_width});
271  }
272  }
273  }
274  return dcs;
275 }
const int kPdgNuE
Definition: PDGCodes.h:28
const int kPdgAntiNuE
Definition: PDGCodes.h:29
std::array< double, 4 > fMixing2s
const int kPdgNuMu
Definition: PDGCodes.h:30
std::pair< std::vector< int >, double > DecayChannel
const int kPdgAntiDarkNeutrino
Definition: PDGCodes.h:223
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
const int kPdgNuTau
Definition: PDGCodes.h:32
const int kPdgDNuMediator
Definition: PDGCodes.h:224
const int kPdgDarkNeutrino
Definition: PDGCodes.h:222
std::vector< GHepParticle > DarkSectorDecayer::Decay ( const GHepParticle mother,
const std::vector< int > &  pdg_daughters 
) const
protected

Definition at line 118 of file DarkSectorDecayer.cxx.

121 {
122  TLorentzVector mother_p4 = *(mother.P4());
123  LOG("DarkSectorDecayer", pINFO)
124  << "Decaying a " << mother.Name()
125  << " with P4 = " << utils::print::P4AsString(&mother_p4);
126 
127  unsigned int nd = pdg_daughters.size();
128  std::vector<double> mass(nd, 0.);
129 
130  for(unsigned int iparticle = 0; iparticle < nd; iparticle++) {
131  TParticlePDG * daughter = PDGLibrary::Instance()->Find(pdg_daughters[iparticle]);
132  assert(daughter);
133 
134  mass[iparticle] = daughter->Mass();
135 
136  SLOG("DarkSectorDecayer", pINFO)
137  << "+ daughter[" << iparticle << "]: "
138  << daughter->GetName() << " (pdg-code = "
139  << pdg_daughters[iparticle] << ", mass = " << mass[iparticle] << ")";
140  }
141 
142  bool is_permitted = fPhaseSpaceGenerator.SetDecay( mother_p4, nd, mass.data() );
143  assert(is_permitted);
144 
145  // Find the maximum phase space decay weight
146  double wmax = -1;
147  for(int i=0; i<50; i++) {
148  double w = fPhaseSpaceGenerator.Generate();
149  wmax = TMath::Max(wmax,w);
150  }
151  assert(wmax>0);
152  LOG("DarkSectorDecayer", pINFO)
153  << "Max phase space gen. weight for current decay: " << wmax;
154 
155  // Generating un-weighted decays
156  RandomGen * rnd = RandomGen::Instance();
157  wmax *= 2;
158  bool accept_decay=false;
159  unsigned int itry=0;
160 
161  while(!accept_decay){
162  itry++;
163  assert(itry<kMaxUnweightDecayIterations);
164 
165  double w = fPhaseSpaceGenerator.Generate();
166  double gw = wmax * rnd->RndDec().Rndm();
167 
168  if(w>wmax) {
169  LOG("DarkSectorDecayer", pWARN)
170  << "Current decay weight = " << w << " > wmax = " << wmax;
171  }
172  LOG("DarkSectorDecayer", pINFO)
173  << "Current decay weight = " << w << " / R = " << gw;
174 
175  accept_decay = (gw<=w);
176  }
177 
178  // A decay was generated - Copy to the event record
179  std::vector<GHepParticle> particles;
180  // Loop over daughter list and add corresponding GHepParticles
181  for(unsigned int id = 0; id < nd; id++) {
182  const TLorentzVector * daughter_p4 = fPhaseSpaceGenerator.GetDecay(id);
183  SLOG("DarkSectorDecayer", pINFO)
184  << "Adding daughter particle with PDG code = " << pdg_daughters[id]
185  << " with P4 = " << utils::print::P4AsShortString( daughter_p4 );
186  SLOG("DarkSectorDecayer", pDEBUG)
187  << "Particle Gun Kinematics: "
188  << "PDG : " << pdg_daughters[id] << ", "
190  GHepStatus_t daughter_status_code = (pdg_daughters[id]==kPdgDNuMediator)
192  particles.push_back(GHepParticle(pdg_daughters[id], daughter_status_code,
193  -1, -1, -1, -1,
194  *daughter_p4, TLorentzVector()));
195  }
196 
197  return particles;
198 }
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:45
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
string Name(void) const
Name that corresponds to the PDG code.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
TGenPhaseSpace fPhaseSpaceGenerator
#define pWARN
Definition: Messenger.h:60
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
const int kPdgDNuMediator
Definition: PDGCodes.h:224
static string ParticleGunKineAsString(const TLorentzVector &vec4)
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:56
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
void DarkSectorDecayer::LoadConfig ( void  )
protectedvirtual

Definition at line 354 of file DarkSectorDecayer.cxx.

355 {
356 
357  // Check particles are in the PDG library, quit if they don't exist
358  constexpr std::array<int, 3> pdgc_mothers = {kPdgDNuMediator,
361  for (auto & pdg_code : pdgc_mothers){
362  TParticlePDG * mother = PDGLibrary::Instance()->Find(pdg_code);
363  if(!mother) {
364  LOG("DarkSectorDecayer", pERROR)
365  << "\n *** The particle with PDG code = " << pdg_code
366  << " was not found in PDGLibrary";
367  exit(78);
368  }
369  }
370 
371  bool good_configuration = true ;
372 
373  double DKineticMixing = 0.; // \varepsilon
374  this->GetParam("Dark-KineticMixing", DKineticMixing);
375  fEps2 = DKineticMixing * DKineticMixing;
376 
377  bool force_unitarity = false ;
378  GetParam("Dark-Mixing-ForceUnitarity", force_unitarity ) ;
379 
380  unsigned int n_min_mixing = force_unitarity ? 3 : 4 ;
381 
382  std::vector<double> DMixing2s; // |U_{\alpha 4}|^2
383  this->GetParamVect("Dark-Mixings2", DMixing2s);
384 
385  // check whether we have enough mixing elements
386  if ( DMixing2s.size () < n_min_mixing ) {
387 
388  good_configuration = false ;
389  LOG("DarkSectorDecayer", pERROR )
390  << "Not enough mixing elements specified, only specified "
391  << DMixing2s.size() << " / " << n_min_mixing ;
392  }
393 
394  double tot_mix = 0. ;
395  for( unsigned int i = 0; i < n_min_mixing ; ++i ) {
396  if ( DMixing2s[i] < 0. ) {
397  good_configuration = false ;
398  LOG("DarkSectorDecayer", pERROR )
399  << "Mixing " << i << " non positive: " << DMixing2s[i] ;
400  continue ;
401  }
402  tot_mix += fMixing2s[i] = DMixing2s[i] ;
403  }
404 
405  if ( force_unitarity ) {
406  fMixing2s[3] = 1. - tot_mix ;
407  }
408 
409  this->GetParam("Dark-Alpha", fAlpha_D);
410 
411  fDNuMass = 0.;
412  this->GetParam("Dark-NeutrinoMass", fDNuMass);
414 
415  fDMediatorMass = 0.;
416  this->GetParam("Dark-MediatorMass", fDMediatorMass);
418 
419  // The model is build on the assumption that the mass of the
420  // mediator is smaller than the mass of the dark neutrino. For the
421  // cross section, this is not a problem.
422  // The decayer though is sensitive to this as the only known decay
423  // amplitude of the dark neutrino requires the mediator in the final
424  // state.
425  // Until the decay amplitude in neutrino is not available
426  // we need to check that the mass hierarchy is respected.
427  if ( fDMediatorMass >= fDNuMass ) {
428  good_configuration = false ;
429  LOG("DarkSectorDecayer", pERROR )
430  << "Dark mediator mass (" << fDMediatorMass
431  << " GeV) too heavy for the dark neutrino ("
432  << fDNuMass << " GeV) to decay" ;
433  }
434 
435  // The other check we need is that the mass of the mediator
436  // has to be smaller than twice the pion mass
437  // Because again in that case we would not be able to
438  // have a proper decay rate since the Mediator would decay in
439  // pion but since we don't have the decay amplitude the
440  // decay rate would be wrong
441  double pion_threshold = PDGLibrary::Instance()->Find( kPdgPiP )->Mass() ;
442  if ( fDMediatorMass >= pion_threshold ) {
443  good_configuration = false ;
444  LOG("DarkSectorDecayer", pERROR )
445  << "Dark mediator mass (" << fDMediatorMass
446  << " GeV) too heavy with respect to the pion decay threshold ("
447  << pion_threshold << " GeV)" ;
448  }
449 
450  if ( ! good_configuration ) {
451  LOG("DarkSectorDecayer", pFATAL ) << "Wrong configuration. Exiting" ;
452  exit ( 78 ) ;
453  }
454 
455 }
#define pERROR
Definition: Messenger.h:59
#define pFATAL
Definition: Messenger.h:56
std::array< double, 4 > fMixing2s
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
const int kPdgAntiDarkNeutrino
Definition: PDGCodes.h:223
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgPiP
Definition: PDGCodes.h:158
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
const int kPdgDNuMediator
Definition: PDGCodes.h:224
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const int kPdgDarkNeutrino
Definition: PDGCodes.h:222
string DarkSectorDecayer::ParticleGunKineAsString ( const TLorentzVector &  vec4)
staticprivate

Definition at line 326 of file DarkSectorDecayer.cxx.

327 {
328  std::ostringstream fmt;
329 
330  double P0 = vec4.Vect().Mag();
331  double thetaYZ = TMath::ASin(vec4.Py()/P0);
332  double thetaXZ = TMath::ASin(vec4.Px()/(P0 * TMath::Cos(thetaYZ)));
333  double rad_to_degrees = 180./kPi;
334 
335  fmt << "P0 = " << P0
336  << ", ThetaXZ = " << thetaXZ*rad_to_degrees
337  << ", ThetaYZ = " << thetaYZ*rad_to_degrees;
338 
339  return fmt.str();
340 }
static const double kPi
Definition: Constants.h:37
void DarkSectorDecayer::ProcessEventRecord ( GHepRecord event) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 60 of file DarkSectorDecayer.cxx.

61 {
62  LOG("DarkSectorDecayer", pINFO)
63  << "Running dark sector decayer ";
64 
65  // Loop over particles, find unstable ones and decay them
66  TObjArrayIter piter(event);
67  GHepParticle * p = 0;
68  int ipos = -1;
69 
70  while( (p = (GHepParticle *) piter.Next()) ) {
71  ipos++;
72  LOG("DarkSectorDecayer", pDEBUG) << "Checking: " << p->Name();
73 
74  if(!this->ToBeDecayed(*p)) continue;
75 
76  GHepParticle& mother = *p; // rename p now we know it will decay
77  std::vector<DarkSectorDecayer::DecayChannel> dcs;
78  int pdg_code = mother.Pdg();
79  if(pdg_code == kPdgDNuMediator){
81  }
82  else if(pdg_code == kPdgDarkNeutrino ||
83  pdg_code == kPdgAntiDarkNeutrino){
84  dcs = DarkNeutrinoDecayChannels(pdg_code);
85  }
86 
87  for ( const auto & dc : dcs ) {
88  std::stringstream amplitudes_msg;
89  amplitudes_msg << "Decay amplitude: " << dc.second << " GeV "
90  << " -> " << 1. / dc.second / units::second << "s for channel [" ;
91  for ( const auto & pdgc : dc.first ) {
92  amplitudes_msg << pdgc << " " ;
93  }
94  amplitudes_msg << "]";
95  LOG("DarkSectorDecayer", pDEBUG) << amplitudes_msg.str();
96  }
97 
98  double total_amplitude = std::accumulate(dcs.begin(), dcs.end(), 0.,
99  [](double total,
101  {return total + dc.second;});
102 
103  int dcid = SelectDecayChannel(dcs, total_amplitude);
104  std::vector<GHepParticle> daughters = Decay(mother, dcs[dcid].first);
105  SetSpaceTime(daughters, mother, total_amplitude);
106 
107  for(auto & daughter: daughters){
108  daughter.SetFirstMother(ipos);
109  event->AddParticle(daughter);
110  }
111  }
112 
113  LOG("DarkSectorDecayer", pNOTICE)
114  << "Done finding & decaying dark sector particles";
115 
116 }
std::pair< std::vector< int >, double > DecayChannel
const int kPdgAntiDarkNeutrino
Definition: PDGCodes.h:223
std::vector< GHepParticle > Decay(const GHepParticle &mother, const std::vector< int > &pdg_daughters) const
void SetSpaceTime(std::vector< GHepParticle > &pp, const GHepParticle &mother, double total_amplitude) const
int Pdg(void) const
Definition: GHepParticle.h:63
string Name(void) const
Name that corresponds to the PDG code.
std::vector< DecayChannel > DarkMediatorDecayChannels(void) const
static constexpr double second
Definition: Units.h:37
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
p
Definition: test.py:223
#define pINFO
Definition: Messenger.h:62
bool ToBeDecayed(const GHepParticle &p) const
const int kPdgDNuMediator
Definition: PDGCodes.h:224
int SelectDecayChannel(const std::vector< DecayChannel > &dcs, double total_amplitude) const
#define pNOTICE
Definition: Messenger.h:61
const int kPdgDarkNeutrino
Definition: PDGCodes.h:222
std::vector< DecayChannel > DarkNeutrinoDecayChannels(int mother_pdg) const
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
#define pDEBUG
Definition: Messenger.h:63
int DarkSectorDecayer::SelectDecayChannel ( const std::vector< DecayChannel > &  dcs,
double  total_amplitude 
) const
protected

Definition at line 200 of file DarkSectorDecayer.cxx.

203 {
204  // Select a decay based on the amplitudes
205  unsigned int ich = 0, sel_ich; // id of selected decay channel
206  RandomGen * rnd = RandomGen::Instance();
207  double x = total_amplitude * rnd->RndDec().Rndm();
208  double partial_sum = 0. ;
209  do {
210  sel_ich = ich;
211  partial_sum += dcs.at(ich++).second;
212  } while (x > partial_sum );
213  return sel_ich;
214 }
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
list x
Definition: train.py:276
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:56
void DarkSectorDecayer::SetSpaceTime ( std::vector< GHepParticle > &  pp,
const GHepParticle mother,
double  total_amplitude 
) const
protected

Definition at line 277 of file DarkSectorDecayer.cxx.

280  {
281 
282  // convert decay amplitude into time
283  const double lifetime = 1e24/units::second/total_amplitude; // units of 10ˆ-24 s
284 
285  RandomGen * rnd = RandomGen::Instance();
286  double t = rnd->RndDec().Exp(lifetime);
287 
288  // t is the decay time in the mother reference frame
289  // it needs to be boosted by a factor gamma
290  t *= mother.P4() -> Gamma() ;
291 
292  // get beta of decaying particle
293  const TLorentzVector & mother_X4 = *(mother.X4());
294  TVector3 mother_boost = mother.P4()->BoostVector();
295 
296  // transport decay_particle with respect to their mother
297  constexpr double speed_of_light = units::second/units::meter; // this gives us the speed of light in m/s
298  TVector3 daughter_position = mother_X4.Vect() + mother_boost * (speed_of_light * t * 1e-9);// in fm
299  TLorentzVector daughter_X4 = TLorentzVector(daughter_position, (mother_X4.T() + t));
300 
301  for(auto & p : pp){
302  p.SetPosition(daughter_X4);
303  }
304 }
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
static constexpr double second
Definition: Units.h:37
const double e
p
Definition: test.py:223
static constexpr double meter
Definition: Units.h:35
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:56
bool DarkSectorDecayer::ToBeDecayed ( const GHepParticle p) const
protected

Definition at line 306 of file DarkSectorDecayer.cxx.

307 {
308  GHepStatus_t status_code = p.Status();
309  if(status_code != kIStDecayedState) return false;
310 
311  int pdg_code = p.Pdg();
312  bool is_handled = false;
313  if(pdg_code == kPdgDNuMediator ||
314  pdg_code == kPdgDarkNeutrino ||
315  pdg_code == kPdgAntiDarkNeutrino){
316  is_handled = true;
317  }
318 
319  LOG("DarkSectorDecayer", pDEBUG)
320  << "Can decay particle with PDG code = " << pdg_code
321  << "? " << ((is_handled)? "Yes" : "No");
322 
323  return is_handled;
324 }
const int kPdgAntiDarkNeutrino
Definition: PDGCodes.h:223
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
int Pdg(void) const
Definition: GHepParticle.h:63
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgDNuMediator
Definition: PDGCodes.h:224
const int kPdgDarkNeutrino
Definition: PDGCodes.h:222
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

double genie::DarkSectorDecayer::fAlpha_D
private

Definition at line 78 of file DarkSectorDecayer.h.

double genie::DarkSectorDecayer::fDMediatorMass
private

Definition at line 81 of file DarkSectorDecayer.h.

double genie::DarkSectorDecayer::fDMediatorMass2
private

Definition at line 81 of file DarkSectorDecayer.h.

double genie::DarkSectorDecayer::fDNuMass
private

Definition at line 80 of file DarkSectorDecayer.h.

double genie::DarkSectorDecayer::fDNuMass2
private

Definition at line 80 of file DarkSectorDecayer.h.

double genie::DarkSectorDecayer::fEps2
private

Definition at line 76 of file DarkSectorDecayer.h.

std::array<double, 4> genie::DarkSectorDecayer::fMixing2s
private

Definition at line 77 of file DarkSectorDecayer.h.

TGenPhaseSpace genie::DarkSectorDecayer::fPhaseSpaceGenerator
mutableprivate

Definition at line 74 of file DarkSectorDecayer.h.


The documentation for this class was generated from the following files: