GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::NNBarOscPrimaryVtxGenerator Class Reference

#include <NNBarOscPrimaryVtxGenerator.h>

Inheritance diagram for genie::NNBarOscPrimaryVtxGenerator:
Inheritance graph
[legend]
Collaboration diagram for genie::NNBarOscPrimaryVtxGenerator:
Collaboration graph
[legend]

Public Member Functions

 NNBarOscPrimaryVtxGenerator ()
 
 NNBarOscPrimaryVtxGenerator (string config)
 
 ~NNBarOscPrimaryVtxGenerator ()
 
void ProcessEventRecord (GHepRecord *event) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
- 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...
 

Private Member Functions

void LoadConfig (void)
 
void AddInitialState (GHepRecord *event) const
 
void GenerateOscillatingNeutronPosition (GHepRecord *event) const
 
void GenerateFermiMomentum (GHepRecord *event) const
 
void GenerateDecayProducts (GHepRecord *event) const
 

Private Attributes

int fCurrInitStatePdg
 
NNBarOscMode_t fCurrDecayMode
 
bool fNucleonIsBound
 
TGenPhaseSpace fPhaseSpaceGenerator
 
const NuclearModelIfNuclModel
 

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)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (const std::string &comm_name)
 
- 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
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (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...
 
- 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< bool > fOwnerships
 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

Definition at line 35 of file NNBarOscPrimaryVtxGenerator.h.

Constructor & Destructor Documentation

NNBarOscPrimaryVtxGenerator::NNBarOscPrimaryVtxGenerator ( )

Definition at line 35 of file NNBarOscPrimaryVtxGenerator.cxx.

35  :
36 EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator")
37 {
38 
39 }
NNBarOscPrimaryVtxGenerator::NNBarOscPrimaryVtxGenerator ( string  config)

Definition at line 41 of file NNBarOscPrimaryVtxGenerator.cxx.

42  :
43 EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator",config)
44 {
45 
46 }
NNBarOscPrimaryVtxGenerator::~NNBarOscPrimaryVtxGenerator ( )

Definition at line 48 of file NNBarOscPrimaryVtxGenerator.cxx.

49 {
50 
51 }

Member Function Documentation

void NNBarOscPrimaryVtxGenerator::AddInitialState ( GHepRecord event) const
private

Definition at line 76 of file NNBarOscPrimaryVtxGenerator.cxx.

References genie::units::A, genie::utils::nnbar_osc::AnnihilatingNucleonPdgCode(), fCurrDecayMode, fCurrInitStatePdg, genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCode(), genie::pdg::IonPdgCodeToA(), genie::pdg::IonPdgCodeToZ(), genie::kIStDecayedState, genie::kIStInitialState, genie::kIStStableFinalState, genie::kPdgNeutron, and genie::kPdgProton.

Referenced by ProcessEventRecord().

77 {
78 
79 // add initial state to event record
80 
81 // event record looks like this:
82 // id: 0, nucleus (kIStInitialState)
83 // |
84 // |---> { id: 1, neutron (kIStDecayedState)
85 // { id: 2, nucleon (kIStDecayedState)
86 // { id: 3, remnant nucleus (kIStStableFinalState)
87 //
88 
89  TLorentzVector v4(0,0,0,0);
90 
94 
95  int ipdg = fCurrInitStatePdg;
96 
97  // add initial nucleus
98  double Mi = PDGLibrary::Instance()->Find(ipdg)->Mass();
99  TLorentzVector p4i(0,0,0,Mi);
100  event->AddParticle(ipdg,stis,-1,-1,-1,-1, p4i, v4);
101 
102  // add oscillating neutron
103  int neutpdg = kPdgNeutron;
104  double mneut = PDGLibrary::Instance()->Find(neutpdg)->Mass();
105  TLorentzVector p4neut(0,0,0,mneut);
106  event->AddParticle(neutpdg,stdc,0,-1,-1,-1, p4neut, v4);
107 
108  // add annihilation nucleon
110  double mn = PDGLibrary::Instance()->Find(dpdg)->Mass();
111  TLorentzVector p4n(0,0,0,mn);
112  event->AddParticle(dpdg,stdc, 0,-1,-1,-1, p4n, v4);
113 
114  // add nuclear remnant
115  int A = pdg::IonPdgCodeToA(ipdg);
116  int Z = pdg::IonPdgCodeToZ(ipdg);
117  A--; A--;
118  if(dpdg == kPdgProton) { Z--; }
119  int rpdg = pdg::IonPdgCode(A, Z);
120  double Mf = PDGLibrary::Instance()->Find(rpdg)->Mass();
121  TLorentzVector p4f(0,0,0,Mf);
122  event->AddParticle(rpdg,stfs,0,-1,-1,-1, p4f, v4);
123 }
int AnnihilatingNucleonPdgCode(NNBarOscMode_t ndm)
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
static constexpr double A
Definition: Units.h:74
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:55
const int kPdgProton
Definition: PDGCodes.h:81
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
const int kPdgNeutron
Definition: PDGCodes.h:83
enum genie::EGHepStatus GHepStatus_t
void NNBarOscPrimaryVtxGenerator::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 502 of file NNBarOscPrimaryVtxGenerator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

503 {
504  Algorithm::Configure(config);
505  this->LoadConfig();
506 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void NNBarOscPrimaryVtxGenerator::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 508 of file NNBarOscPrimaryVtxGenerator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

509 {
510  Algorithm::Configure(config);
511  this->LoadConfig();
512 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void NNBarOscPrimaryVtxGenerator::GenerateDecayProducts ( GHepRecord event) const
private

accept_decay

Definition at line 267 of file NNBarOscPrimaryVtxGenerator.cxx.

References genie::XclsTag::DecayMode(), genie::utils::nnbar_osc::DecayProductList(), genie::utils::nnbar_osc::DecayProductStatus(), genie::Interaction::ExclTag(), fCurrDecayMode, fCurrInitStatePdg, genie::PDGLibrary::Find(), fNucleonIsBound, fPhaseSpaceGenerator, genie::GHepParticle::GetP4(), genie::GHepParticle::GetX4(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::controls::kMaxUnweightDecayIterations, genie::kPdgNeutron, LOG, genie::units::m, genie::Interaction::NOsc(), genie::utils::print::P4AsString(), genie::GHepParticle::Pdg(), pERROR, pINFO, pNOTICE, pWARN, genie::RandomGen::RndHadro(), genie::RandomGen::RndNum(), genie::Target::SetHitNucPdg(), genie::exceptions::EVGThreadException::SetReason(), genie::RandomGen::SetSeed(), and genie::exceptions::EVGThreadException::SwitchOnFastForward().

Referenced by ProcessEventRecord().

269 {
270  LOG("NNBarOsc", pINFO) << "Generating decay...";
271 
273  LOG("NNBarOsc", pINFO) << "Decay product IDs: " << pdgv;
274  assert ( pdgv.size() > 1);
275 
276  LOG("NNBarOsc", pINFO) << "Performing a phase space decay...";
277 
278  // Get the decay product masses
279 
280  vector<int>::const_iterator pdg_iter;
281  int idx = 0;
282  double * mass = new double[pdgv.size()];
283  double sum = 0;
284  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
285  int pdgc = *pdg_iter;
286  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
287  mass[idx++] = m;
288  sum += m;
289  }
290 
291  LOG("NNBarOsc", pINFO)
292  << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
293  int initial_nucleus_id = 0;
294  int oscillating_neutron_id = 1;
295  int annihilation_nucleon_id = 2;
296 
297  // get our annihilating nucleons
298  GHepParticle * initial_nucleus = event->Particle(initial_nucleus_id);
299  assert(initial_nucleus);
300  GHepParticle * oscillating_neutron = event->Particle(oscillating_neutron_id);
301  assert(oscillating_neutron);
302  GHepParticle * annihilation_nucleon = event->Particle(annihilation_nucleon_id);
303  assert(annihilation_nucleon);
304 
305  Target tgt(initial_nucleus->Pdg());
307 
308  // get their momentum 4-vectors and boost into rest frame
309  TLorentzVector * p4_1 = oscillating_neutron->GetP4();
310  TLorentzVector * p4_2 = annihilation_nucleon->GetP4();
311  TLorentzVector * p4d = new TLorentzVector(*p4_1 + *p4_2);
312  TVector3 boost = p4d->BoostVector();
313  p4d->Boost(-boost);
314 
315  // get decay position
316  TLorentzVector * v4d = annihilation_nucleon->GetX4();
317 
318  delete p4_1;
319  delete p4_2;
320 
321  LOG("NNBarOsc", pINFO)
322  << "Decaying system p4 = " << utils::print::P4AsString(p4d);
323 
324  // Set the decay
325  bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
326 
327  // If the decay is not energetically allowed, select a new final state
328  while(!permitted) {
329 
330  LOG("NNBarOsc", pINFO)
331  << "Not enough energy to generate decay products! Selecting a new final state.";
332 
333  int mode = 0;
334 
335  int initial_nucleus_pdg = initial_nucleus->Pdg();
336 
337  std::string nucleus_pdg = std::to_string(static_cast<long long>(initial_nucleus_pdg));
338  if (nucleus_pdg.size() != 10) {
339  LOG("NNBarOsc", pERROR)
340  << "Expecting the nuclear PDG code to be a 10-digit integer, but it is " << nucleus_pdg << ", which is "
341  << nucleus_pdg.size() << " digits long. Drop me an email at jezhewes@gmail.com ; exiting...";
342  exit(1);
343  }
344 
345  int n_nucleons = std::stoi(nucleus_pdg.substr(6,3)) - 1;
346  int n_protons = std::stoi(nucleus_pdg.substr(3,3));
347  double proton_frac = ((double)n_protons) / ((double) n_nucleons);
348  double neutron_frac = 1 - proton_frac;
349 
350  // set branching ratios, taken from bubble chamber data
351  const int n_modes = 16;
352  double br [n_modes] = { 0.010, 0.080, 0.100, 0.220,
353  0.360, 0.160, 0.070, 0.020,
354  0.015, 0.065, 0.110, 0.280,
355  0.070, 0.240, 0.100, 0.100 };
356 
357  for (int i = 0; i < n_modes; i++) {
358  // for first 7 branching ratios, multiply by relative proton density
359  if (i < 7)
360  br[i] *= proton_frac;
361  // for next 9, multiply by relative neutron density
362  else
363  br[i] *= neutron_frac;
364  }
365 
366  // randomly generate a number between 1 and 0
367  RandomGen * rnd = RandomGen::Instance();
368  rnd->SetSeed(0);
369  double p = rnd->RndNum().Rndm();
370 
371  // loop through all modes, figure out which one our random number corresponds to
372  double threshold = 0;
373  for (int j = 0; j < n_modes; j++) {
374  threshold += br[j];
375  if (p < threshold) {
376  // once we've found our mode, stop looping
377  mode = j + 1;
378  break;
379  }
380  }
381 
382  // create new event record with new final state
383  // TODO - I don't think Jeremy meant to make a _new_ record here; it
384  // shadows the one passed in...
385  // EventRecord * event = new EventRecord;
386  Interaction * interaction = Interaction::NOsc((int)fCurrInitStatePdg, mode);
387  event->AttachSummary(interaction);
388 
389  fCurrDecayMode = (NNBarOscMode_t) interaction->ExclTag().DecayMode();
390 
392  LOG("NNBarOsc", pINFO) << "Decay product IDs: " << pdgv;
393  assert ( pdgv.size() > 1);
394 
395  // get the decay particles again
396  LOG("NNBarOsc", pINFO) << "Performing a phase space decay...";
397  idx = 0;
398  delete [] mass;
399  mass = new double[pdgv.size()];
400  sum = 0;
401  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
402  int pdgc = *pdg_iter;
403  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
404  mass[idx++] = m;
405  sum += m;
406  }
407 
408  LOG("NNBarOsc", pINFO)
409  << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
410  LOG("NNBarOsc", pINFO)
411  << "Decaying system p4 = " << utils::print::P4AsString(p4d);
412 
413  permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
414  }
415 
416  if(!permitted) {
417  LOG("NNBarOsc", pERROR)
418  << " *** Phase space decay is not permitted \n"
419  << " Total particle mass = " << sum << "\n"
420  << " Decaying system p4 = " << utils::print::P4AsString(p4d);
421  // clean-up
422  delete [] mass;
423  delete p4d;
424  delete v4d;
425  // throw exception
427  exception.SetReason("Decay not permitted kinematically");
428  exception.SwitchOnFastForward();
429  throw exception;
430  }
431 
432  // Get the maximum weight
433  //double wmax = fPhaseSpaceGenerator.GetWtMax();
434  double wmax = -1;
435  for(int i=0; i<200; i++) {
436  double w = fPhaseSpaceGenerator.Generate();
437  wmax = TMath::Max(wmax,w);
438  }
439  assert(wmax>0);
440  wmax *= 2;
441 
442  LOG("NNBarOsc", pNOTICE)
443  << "Max phase space gen. weight @ current hadronic system: " << wmax;
444 
445  // Generate an unweighted decay
446  RandomGen * rnd = RandomGen::Instance();
447 
448  bool accept_decay=false;
449  unsigned int itry=0;
450  while(!accept_decay)
451  {
452  itry++;
453 
455  // report, clean-up and return
456  LOG("NNBarOsc", pWARN)
457  << "Couldn't generate an unweighted phase space decay after "
458  << itry << " attempts";
459  // clean up
460  delete [] mass;
461  delete p4d;
462  delete v4d;
463  // throw exception
465  exception.SetReason("Couldn't select decay after N attempts");
466  exception.SwitchOnFastForward();
467  throw exception;
468  }
469  double w = fPhaseSpaceGenerator.Generate();
470  if(w > wmax) {
471  LOG("NNBarOsc", pWARN)
472  << "Decay weight = " << w << " > max decay weight = " << wmax;
473  }
474  double gw = wmax * rnd->RndHadro().Rndm();
475  accept_decay = (gw<=w);
476 
477  LOG("NNBarOsc", pINFO)
478  << "Decay weight = " << w << " / R = " << gw
479  << " - accepted: " << accept_decay;
480 
481  } //!accept_decay
482 
483  // Insert final state products into a TClonesArray of GHepParticle's
484  TLorentzVector v4(*v4d);
485  int idp = 0;
486  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
487  int pdgc = *pdg_iter;
488  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
489  GHepStatus_t ist =
491  p4fin->Boost(boost);
492  event->AddParticle(pdgc, ist, oscillating_neutron_id,-1,-1,-1, *p4fin, v4);
493  idp++;
494  }
495 
496  // Clean-up
497  delete [] mass;
498  delete p4d;
499  delete v4d;
500 }
TLorentzVector * GetX4(void) const
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
GHepStatus_t DecayProductStatus(bool in_nucleus, int pdgc)
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
A list of PDG codes.
Definition: PDGCodeList.h:32
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
enum genie::ENNBarOscMode NNBarOscMode_t
TLorentzVector * GetP4(void) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int DecayMode(void) const
Definition: XclsTag.h:70
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
TRandom3 & RndNum(void) const
rnd number generator used by MC integrators &amp; other numerical methods
Definition: RandomGen.h:77
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
PDGCodeList DecayProductList(NNBarOscMode_t ndm)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
const int kPdgNeutron
Definition: PDGCodes.h:83
static constexpr double m
Definition: Units.h:71
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
static Interaction * NOsc(int tgt, int annihilation_mode=-1)
enum genie::EGHepStatus GHepStatus_t
void SetSeed(long int seed)
Definition: RandomGen.cxx:82
void NNBarOscPrimaryVtxGenerator::GenerateFermiMomentum ( GHepRecord event) const
private

Definition at line 201 of file NNBarOscPrimaryVtxGenerator.cxx.

References genie::units::A, genie::GHepParticle::A(), fNuclModel, genie::NuclearModelI::GenerateNucleon(), genie::GHepParticle::GetP4(), genie::kPdgNeutron, LOG, genie::GHepParticle::Mass(), genie::NuclearModelI::Momentum3(), genie::GHepParticle::Pdg(), pINFO, genie::NuclearModelI::RemovalEnergy(), genie::Target::SetHitNucPdg(), and genie::GHepParticle::SetMomentum().

Referenced by ProcessEventRecord().

203 {
204  GHepParticle * initial_nucleus = event->Particle(0);
205  int A = initial_nucleus->A();
206  if(A <= 2) {
207  return;
208  }
209 
210  GHepParticle * oscillating_neutron = event->Particle(1);
211  GHepParticle * annihilation_nucleon = event->Particle(2);
212  GHepParticle * remnant_nucleus = event->Particle(3);
213  assert(oscillating_neutron);
214  assert(annihilation_nucleon);
215  assert(remnant_nucleus);
216 
217  // generate a Fermi momentum & removal energy
218  Target tgt(initial_nucleus->Pdg());
219 
220  // start with oscillating neutron
222  // generate nuclear model & fermi momentum
224  TVector3 p3 = fNuclModel->Momentum3();
225  double w = fNuclModel->RemovalEnergy();
226 
227  // use mass & momentum to calculate energy
228  double mass = oscillating_neutron->Mass();
229  double energy = sqrt(pow(mass,2) + p3.Mag2()) - w;
230  // give new energy & momentum to particle
231  TLorentzVector p4(p3, energy);
232  oscillating_neutron->SetMomentum(p4);
233 
234  LOG("NNBarOsc", pINFO)
235  << "Generated neutron momentum: ("
236  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
237  << "|p| = " << p3.Mag();
238 
239  // then rinse repeat for the annihilation nucleon
240  tgt.SetHitNucPdg(annihilation_nucleon->Pdg());
241  // use nuclear model to generate fermi momentum
243  p3 = fNuclModel->Momentum3();
244  w = fNuclModel->RemovalEnergy();
245  // use mass & momentum to figure out energy
246  mass = annihilation_nucleon->Mass();
247  energy = sqrt(pow(mass,2) + p3.Mag2()) - w;
248  // give new energy & momentum to particle
249  p4 = TLorentzVector(p3, energy);
250  annihilation_nucleon->SetMomentum(p4);
251 
252  LOG("NNBarOsc", pINFO)
253  << "Generated nucleon momentum: ("
254  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
255  << "|p| = " << p3.Mag();
256 
257  // now figure out momentum for the nuclear remnant
258  p3 = -1 * (oscillating_neutron->GetP4()->Vect() + annihilation_nucleon->GetP4()->Vect());
259  // figure out energy from mass & momentum
260  mass = remnant_nucleus->Mass();
261  energy = sqrt(pow(mass,2) + p3.Mag2());
262  // give new energy & momentum to remnant
263  p4 = TLorentzVector(p3, energy);
264  remnant_nucleus->SetMomentum(p4);
265 }
double RemovalEnergy(void) const
Definition: NuclearModelI.h:65
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
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
static constexpr double A
Definition: Units.h:74
TLorentzVector * GetP4(void) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
#define pINFO
Definition: Messenger.h:62
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
int A(void) const
virtual bool GenerateNucleon(const Target &) const =0
const int kPdgNeutron
Definition: PDGCodes.h:83
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
void NNBarOscPrimaryVtxGenerator::GenerateOscillatingNeutronPosition ( GHepRecord event) const
private

Definition at line 125 of file NNBarOscPrimaryVtxGenerator.cxx.

References genie::units::A, genie::GHepParticle::A(), genie::utils::nuclear::Density(), genie::RandomGen::Instance(), genie::constants::kPi, genie::controls::kRjMaxIterations, LOG, pERROR, pINFO, pWARN, genie::RandomGen::RndFsi(), genie::GHepParticle::SetPosition(), genie::exceptions::EVGThreadException::SetReason(), and genie::exceptions::EVGThreadException::SwitchOnFastForward().

Referenced by ProcessEventRecord().

127 {
128  GHepParticle * initial_nucleus = event->Particle(0);
129  int A = initial_nucleus->A();
130  if(A <= 2) {
131  return;
132  }
133 
134  RandomGen * rnd = RandomGen::Instance();
135 
136  double R0 = 1.3;
137  double dA = (double)A;
138  double R = R0 * TMath::Power(dA, 1./3.);
139 
140  LOG("NNBarOsc", pINFO)
141  << "Generating vertex according to a realistic nuclear density profile";
142 
143  // get inputs to the rejection method
144  double ymax = -1;
145  double rmax = 3*R;
146  double dr = R/40.;
147  for(double r = 0; r < rmax; r+=dr) {
148  ymax = TMath::Max(ymax, r*r * utils::nuclear::Density(r,A));
149  }
150  ymax *= 1.2;
151 
152  // select a vertex using the rejection method
153  TLorentzVector vtx(0,0,0,0);
154  unsigned int iter = 0;
155  while(1) {
156  iter++;
157 
158  // throw an exception if it hasn't find a solution after many attempts
159  if(iter > controls::kRjMaxIterations) {
160  LOG("NNBarOsc", pWARN)
161  << "Couldn't generate a vertex position after " << iter << " iterations";
163  exception.SetReason("Couldn't generate vertex");
164  exception.SwitchOnFastForward();
165  throw exception;
166  }
167 
168  double r = rmax * rnd->RndFsi().Rndm();
169  double t = ymax * rnd->RndFsi().Rndm();
170  double y = r*r * utils::nuclear::Density(r,A);
171  if(y > ymax) {
172  LOG("NNBarOsc", pERROR)
173  << "y = " << y << " > ymax = " << ymax << " for r = " << r << ", A = " << A;
174  }
175  bool accept = (t < y);
176  if(accept) {
177  double phi = 2*constants::kPi * rnd->RndFsi().Rndm();
178  double cosphi = TMath::Cos(phi);
179  double sinphi = TMath::Sin(phi);
180  double costheta = -1 + 2 * rnd->RndFsi().Rndm();
181  double sintheta = TMath::Sqrt(1-costheta*costheta);
182  vtx.SetX(r*sintheta*cosphi);
183  vtx.SetY(r*sintheta*sinphi);
184  vtx.SetZ(r*costheta);
185  vtx.SetT(0.);
186  break;
187  }
188  } // while 1
189 
190  // giving position to oscillating neutron
191  GHepParticle * oscillating_neutron = event->Particle(1);
192  assert(oscillating_neutron);
193  oscillating_neutron->SetPosition(vtx);
194 
195  // give same position to annihilation nucleon
196  GHepParticle * annihilation_nucleon = event->Particle(2);
197  assert(annihilation_nucleon);
198  annihilation_nucleon->SetPosition(vtx);
199 }
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
double Density(double r, int A, double ring=0.)
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
void SetPosition(const TLorentzVector &v4)
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
int A(void) const
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
void NNBarOscPrimaryVtxGenerator::LoadConfig ( void  )
private

Definition at line 514 of file NNBarOscPrimaryVtxGenerator.cxx.

References fNuclModel, and genie::Algorithm::SubAlg().

Referenced by Configure().

515 {
516 // AlgConfigPool * confp = AlgConfigPool::Instance();
517 // const Registry * gc = confp->GlobalParameterList();
518 
519  fNuclModel = 0;
520 
521  RgKey nuclkey = "NuclearModel";
522  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
523  assert(fNuclModel);
524 }
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
string RgKey
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
void NNBarOscPrimaryVtxGenerator::ProcessEventRecord ( GHepRecord event) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 53 of file NNBarOscPrimaryVtxGenerator.cxx.

References AddInitialState(), genie::utils::nnbar_osc::AsString(), genie::XclsTag::DecayMode(), genie::Interaction::ExclTag(), fCurrDecayMode, fCurrInitStatePdg, fNucleonIsBound, GenerateDecayProducts(), GenerateFermiMomentum(), GenerateOscillatingNeutronPosition(), genie::Interaction::InitState(), genie::pdg::IonPdgCodeToA(), LOG, genie::Target::Pdg(), pNOTICE, and genie::InitialState::Tgt().

54 {
55  // spit out some output
56  Interaction * interaction = event->Summary();
57  fCurrInitStatePdg = interaction->InitState().Tgt().Pdg();
58  fCurrDecayMode = (NNBarOscMode_t) interaction->ExclTag().DecayMode();
59 
60  // spit out that info -j
61  LOG("NNBarOsc", pNOTICE)
62  << "Simulating decay " << genie::utils::nnbar_osc::AsString(fCurrDecayMode)
63  << " for an initial state with code: " << fCurrInitStatePdg;
64 
65  // check if nucleon is bound
67  // can take this out for nnbar, nucleon is always bound!
68 
69  // now call these four functions!
70  this->AddInitialState(event);
72  this->GenerateFermiMomentum(event);
73  this->GenerateDecayProducts(event);
74 }
int Pdg(void) const
Definition: Target.h:71
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
enum genie::ENNBarOscMode NNBarOscMode_t
void AddInitialState(GHepRecord *event) const
int DecayMode(void) const
Definition: XclsTag.h:70
string AsString(NNBarOscMode_t ndm)
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
const InitialState & InitState(void) const
Definition: Interaction.h:69
void GenerateOscillatingNeutronPosition(GHepRecord *event) const
#define pNOTICE
Definition: Messenger.h:61
const Target & Tgt(void) const
Definition: InitialState.h:66
void GenerateFermiMomentum(GHepRecord *event) const
void GenerateDecayProducts(GHepRecord *event) const

Member Data Documentation

NNBarOscMode_t genie::NNBarOscPrimaryVtxGenerator::fCurrDecayMode
mutableprivate
int genie::NNBarOscPrimaryVtxGenerator::fCurrInitStatePdg
mutableprivate
bool genie::NNBarOscPrimaryVtxGenerator::fNucleonIsBound
mutableprivate

Definition at line 60 of file NNBarOscPrimaryVtxGenerator.h.

Referenced by GenerateDecayProducts(), and ProcessEventRecord().

const NuclearModelI* genie::NNBarOscPrimaryVtxGenerator::fNuclModel
private

Definition at line 63 of file NNBarOscPrimaryVtxGenerator.h.

Referenced by GenerateFermiMomentum(), and LoadConfig().

TGenPhaseSpace genie::NNBarOscPrimaryVtxGenerator::fPhaseSpaceGenerator
mutableprivate

Definition at line 61 of file NNBarOscPrimaryVtxGenerator.h.

Referenced by GenerateDecayProducts().


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