GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DarkSectorDecayer.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Author: Iker de Icaza <i.de-icaza-astiz \at sussex.ac.uk>
7  University of Sussex
8 
9  Costas Andreopoulos <c.andreopoulos \at cern.ch>
10  University of Liverpool
11 */
12 //____________________________________________________________________________
13 
14 #include <cmath>
15 #include <numeric>
16 #include <sstream>
17 
18 #include <TMath.h>
19 
35 #include "Math/GSLMinimizer.h"
36 #include <Math/WrappedParamFunction.h>
37 
38 using namespace genie;
39 using namespace genie::controls;
40 using namespace genie::constants;
41 
42 //____________________________________________________________________________
44  EventRecordVisitorI("genie::DarkSectorDecayer")
45 {
46 
47 }
48 //____________________________________________________________________________
50  EventRecordVisitorI("genie::DarkSectorDecayer", config)
51 {
52 
53 }
54 //____________________________________________________________________________
56 {
57 
58 }
59 //____________________________________________________________________________
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 }
117 //____________________________________________________________________________
118 std::vector<GHepParticle> DarkSectorDecayer::Decay(
119  const GHepParticle & mother,
120  const std::vector<int> & pdg_daughters) const
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 }
199 //____________________________________________________________________________
201  const std::vector<DarkSectorDecayer::DecayChannel> & dcs,
202  const double total_amplitude) const
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 }
215 //____________________________________________________________________________
216 std::vector<DarkSectorDecayer::DecayChannel> DarkSectorDecayer::DarkMediatorDecayChannels(void) const
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 }
248 //____________________________________________________________________________
249 std::vector<DarkSectorDecayer::DecayChannel> DarkSectorDecayer::DarkNeutrinoDecayChannels(int mother_pdg) const
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 }
276 //____________________________________________________________________________
278  std::vector<GHepParticle> & pp,
279  const GHepParticle & mother,
280  double total_amplitude) const {
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 }
305 //____________________________________________________________________________
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 }
325 //____________________________________________________________________________
326 string DarkSectorDecayer::ParticleGunKineAsString(const TLorentzVector & vec4)
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 }
341 //____________________________________________________________________________
343 {
344  Algorithm::Configure(config);
345  this->LoadConfig();
346 }
347 //____________________________________________________________________________
348 void DarkSectorDecayer::Configure(string config)
349 {
350  Algorithm::Configure(config);
351  this->LoadConfig();
352 }
353 //____________________________________________________________________________
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 }
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:45
const int kPdgNuE
Definition: PDGCodes.h:28
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
const int kPdgAntiNuE
Definition: PDGCodes.h:29
#define pFATAL
Definition: Messenger.h:56
virtual void Configure(const Registry &config)
std::array< double, 4 > fMixing2s
const int kPdgNuMu
Definition: PDGCodes.h:30
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 kPdgAntiMuon
Definition: PDGCodes.h:38
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
const int kPdgElectron
Definition: PDGCodes.h:35
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
std::pair< std::vector< int >, double > DecayChannel
const int kPdgAntiDarkNeutrino
Definition: PDGCodes.h:222
std::vector< GHepParticle > Decay(const GHepParticle &mother, const std::vector< int > &pdg_daughters) const
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
void SetSpaceTime(std::vector< GHepParticle > &pp, const GHepParticle &mother, double total_amplitude) const
int Pdg(void) const
Definition: GHepParticle.h:63
void ProcessEventRecord(GHepRecord *event) const
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
const double e
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual void LoadConfig(void)
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
const int kPdgPiP
Definition: PDGCodes.h:158
#define pINFO
Definition: Messenger.h:62
TGenPhaseSpace fPhaseSpaceGenerator
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:60
bool ToBeDecayed(const GHepParticle &p) const
const int kPdgNuTau
Definition: PDGCodes.h:32
static constexpr double meter
Definition: Units.h:35
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const int kPdgDNuMediator
Definition: PDGCodes.h:223
static string ParticleGunKineAsString(const TLorentzVector &vec4)
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
int SelectDecayChannel(const std::vector< DecayChannel > &dcs, double total_amplitude) const
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const int kPdgMuon
Definition: PDGCodes.h:37
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
const int kPdgPositron
Definition: PDGCodes.h:36
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
const int kPdgDarkNeutrino
Definition: PDGCodes.h:221
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:56
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
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63