GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NNBarOscPrimaryVtxGenerator.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  Jeremy Hewes, Georgia Karagiorgi
7  University of Manchester
8 */
9 //____________________________________________________________________________
10 
31 
32 using namespace genie;
33 
34 //____________________________________________________________________________
36 EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator")
37 {
38 
39 }
40 //____________________________________________________________________________
42  string config) :
43 EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator",config)
44 {
45 
46 }
47 //____________________________________________________________________________
49 {
50 
51 }
52 //____________________________________________________________________________
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 }
75 //____________________________________________________________________________
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 }
124 //____________________________________________________________________________
126  GHepRecord * event) const
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 }
200 //____________________________________________________________________________
202  GHepRecord * event) const
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 }
266 //____________________________________________________________________________
268  GHepRecord * event) const
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 }
501 //___________________________________________________________________________
503 {
504  Algorithm::Configure(config);
505  this->LoadConfig();
506 }
507 //___________________________________________________________________________
509 {
510  Algorithm::Configure(config);
511  this->LoadConfig();
512 }
513 //___________________________________________________________________________
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 }
525 //___________________________________________________________________________
TLorentzVector * GetX4(void) const
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
double RemovalEnergy(void) const
Definition: NuclearModelI.h:65
int AnnihilatingNucleonPdgCode(NNBarOscMode_t ndm)
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double Density(double r, int A, double ring=0.)
GHepStatus_t DecayProductStatus(bool in_nucleus, int pdgc)
int Pdg(void) const
Definition: Target.h:71
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
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
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
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
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void AddInitialState(GHepRecord *event) const
int DecayMode(void) const
Definition: XclsTag.h:70
string AsString(NNBarOscMode_t ndm)
#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
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
void ProcessEventRecord(GHepRecord *event) const
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
const InitialState & InitState(void) const
Definition: Interaction.h:69
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
void GenerateOscillatingNeutronPosition(GHepRecord *event) const
int A(void) const
#define pNOTICE
Definition: Messenger.h:61
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
virtual bool GenerateNucleon(const Target &) const =0
const int kPdgNeutron
Definition: PDGCodes.h:83
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
static constexpr double m
Definition: Units.h:71
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
void GenerateFermiMomentum(GHepRecord *event) const
static Interaction * NOsc(int tgt, int annihilation_mode=-1)
enum genie::EGHepStatus GHepStatus_t
void GenerateDecayProducts(GHepRecord *event) const
void SetSeed(long int seed)
Definition: RandomGen.cxx:82
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345