NNBarOscPrimaryVtxGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Jeremy Hewes, Georgia Karagiorgi
8  University of Manchester
9 
10  For documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14 */
15 //____________________________________________________________________________
16 
37 
38 using namespace genie;
39 
40 //____________________________________________________________________________
42 EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator")
43 {
44 
45 }
46 //____________________________________________________________________________
48  string config) :
49 EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator",config)
50 {
51 
52 }
53 //____________________________________________________________________________
55 {
56 
57 }
58 //____________________________________________________________________________
60 {
61  // spit out some output
62  Interaction * interaction = event->Summary();
63  fCurrInitStatePdg = interaction->InitState().Tgt().Pdg();
64  fCurrDecayMode = (NNBarOscMode_t) interaction->ExclTag().DecayMode();
65 
66  // spit out that info -j
67  LOG("NNBarOsc", pNOTICE)
68  << "Simulating decay " << genie::utils::nnbar_osc::AsString(fCurrDecayMode)
69  << " for an initial state with code: " << fCurrInitStatePdg;
70 
71  // check if nucleon is bound
73  // can take this out for nnbar, nucleon is always bound!
74 
75  // now call these four functions!
76  this->AddInitialState(event);
78  this->GenerateFermiMomentum(event);
79  this->GenerateDecayProducts(event);
80 }
81 //____________________________________________________________________________
83 {
84 
85 // add initial state to event record
86 
87 // event record looks like this:
88 // id: 0, nucleus (kIStInitialState)
89 // |
90 // |---> { id: 1, neutron (kIStDecayedState)
91 // { id: 2, nucleon (kIStDecayedState)
92 // { id: 3, remnant nucleus (kIStStableFinalState)
93 //
94 
95  TLorentzVector v4(0,0,0,0);
96 
100 
101  int ipdg = fCurrInitStatePdg;
102 
103  // add initial nucleus
104  double Mi = PDGLibrary::Instance()->Find(ipdg)->Mass();
105  TLorentzVector p4i(0,0,0,Mi);
106  event->AddParticle(ipdg,stis,-1,-1,-1,-1, p4i, v4);
107 
108  // add oscillating neutron
109  int neutpdg = kPdgNeutron;
110  double mneut = PDGLibrary::Instance()->Find(neutpdg)->Mass();
111  TLorentzVector p4neut(0,0,0,mneut);
112  event->AddParticle(neutpdg,stdc,0,-1,-1,-1, p4neut, v4);
113 
114  // add annihilation nucleon
116  double mn = PDGLibrary::Instance()->Find(dpdg)->Mass();
117  TLorentzVector p4n(0,0,0,mn);
118  event->AddParticle(dpdg,stdc, 0,-1,-1,-1, p4n, v4);
119 
120  // add nuclear remnant
121  int A = pdg::IonPdgCodeToA(ipdg);
122  int Z = pdg::IonPdgCodeToZ(ipdg);
123  A--; A--;
124  if(dpdg == kPdgProton) { Z--; }
125  int rpdg = pdg::IonPdgCode(A, Z);
126  double Mf = PDGLibrary::Instance()->Find(rpdg)->Mass();
127  TLorentzVector p4f(0,0,0,Mf);
128  event->AddParticle(rpdg,stfs,0,-1,-1,-1, p4f, v4);
129 }
130 //____________________________________________________________________________
132  GHepRecord * event) const
133 {
134  GHepParticle * initial_nucleus = event->Particle(0);
135  int A = initial_nucleus->A();
136  if(A <= 2) {
137  return;
138  }
139 
141 
142  double R0 = 1.3;
143  double dA = (double)A;
144  double R = R0 * TMath::Power(dA, 1./3.);
145 
146  LOG("NNBarOsc", pINFO)
147  << "Generating vertex according to a realistic nuclear density profile";
148 
149  // get inputs to the rejection method
150  double ymax = -1;
151  double rmax = 3*R;
152  double dr = R/40.;
153  for(double r = 0; r < rmax; r+=dr) {
154  ymax = TMath::Max(ymax, r*r * utils::nuclear::Density(r,A));
155  }
156  ymax *= 1.2;
157 
158  // select a vertex using the rejection method
159  TLorentzVector vtx(0,0,0,0);
160  unsigned int iter = 0;
161  while(1) {
162  iter++;
163 
164  // throw an exception if it hasn't find a solution after many attempts
165  if(iter > controls::kRjMaxIterations) {
166  LOG("NNBarOsc", pWARN)
167  << "Couldn't generate a vertex position after " << iter << " iterations";
169  exception.SetReason("Couldn't generate vertex");
170  exception.SwitchOnFastForward();
171  throw exception;
172  }
173 
174  double r = rmax * rnd->RndFsi().Rndm();
175  double t = ymax * rnd->RndFsi().Rndm();
176  double y = r*r * utils::nuclear::Density(r,A);
177  if(y > ymax) {
178  LOG("NNBarOsc", pERROR)
179  << "y = " << y << " > ymax = " << ymax << " for r = " << r << ", A = " << A;
180  }
181  bool accept = (t < y);
182  if(accept) {
183  double phi = 2*constants::kPi * rnd->RndFsi().Rndm();
184  double cosphi = TMath::Cos(phi);
185  double sinphi = TMath::Sin(phi);
186  double costheta = -1 + 2 * rnd->RndFsi().Rndm();
187  double sintheta = TMath::Sqrt(1-costheta*costheta);
188  vtx.SetX(r*sintheta*cosphi);
189  vtx.SetY(r*sintheta*sinphi);
190  vtx.SetZ(r*costheta);
191  vtx.SetT(0.);
192  break;
193  }
194  } // while 1
195 
196  // giving position to oscillating neutron
197  GHepParticle * oscillating_neutron = event->Particle(1);
198  assert(oscillating_neutron);
199  oscillating_neutron->SetPosition(vtx);
200 
201  // give same position to annihilation nucleon
202  GHepParticle * annihilation_nucleon = event->Particle(2);
203  assert(annihilation_nucleon);
204  annihilation_nucleon->SetPosition(vtx);
205 }
206 //____________________________________________________________________________
208  GHepRecord * event) const
209 {
210  GHepParticle * initial_nucleus = event->Particle(0);
211  int A = initial_nucleus->A();
212  if(A <= 2) {
213  return;
214  }
215 
216  GHepParticle * oscillating_neutron = event->Particle(1);
217  GHepParticle * annihilation_nucleon = event->Particle(2);
218  GHepParticle * remnant_nucleus = event->Particle(3);
219  assert(oscillating_neutron);
220  assert(annihilation_nucleon);
221  assert(remnant_nucleus);
222 
223  // generate a Fermi momentum & removal energy
224  Target tgt(initial_nucleus->Pdg());
225 
226  // start with oscillating neutron
228  // generate nuclear model & fermi momentum
230  TVector3 p3 = fNuclModel->Momentum3();
231  double w = fNuclModel->RemovalEnergy();
232 
233  // use mass & momentum to calculate energy
234  double mass = oscillating_neutron->Mass();
235  double energy = sqrt(pow(mass,2) + p3.Mag2()) - w;
236  // give new energy & momentum to particle
237  TLorentzVector p4(p3, energy);
238  oscillating_neutron->SetMomentum(p4);
239 
240  LOG("NNBarOsc", pINFO)
241  << "Generated neutron momentum: ("
242  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
243  << "|p| = " << p3.Mag();
244 
245  // then rinse repeat for the annihilation nucleon
246  tgt.SetHitNucPdg(annihilation_nucleon->Pdg());
247  // use nuclear model to generate fermi momentum
249  p3 = fNuclModel->Momentum3();
250  w = fNuclModel->RemovalEnergy();
251  // use mass & momentum to figure out energy
252  mass = annihilation_nucleon->Mass();
253  energy = sqrt(pow(mass,2) + p3.Mag2()) - w;
254  // give new energy & momentum to particle
255  p4 = TLorentzVector(p3, energy);
256  annihilation_nucleon->SetMomentum(p4);
257 
258  LOG("NNBarOsc", pINFO)
259  << "Generated nucleon momentum: ("
260  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
261  << "|p| = " << p3.Mag();
262 
263  // now figure out momentum for the nuclear remnant
264  p3 = -1 * (oscillating_neutron->GetP4()->Vect() + annihilation_nucleon->GetP4()->Vect());
265  // figure out energy from mass & momentum
266  mass = remnant_nucleus->Mass();
267  energy = sqrt(pow(mass,2) + p3.Mag2());
268  // give new energy & momentum to remnant
269  p4 = TLorentzVector(p3, energy);
270  remnant_nucleus->SetMomentum(p4);
271 }
272 //____________________________________________________________________________
274  GHepRecord * event) const
275 {
276  LOG("NNBarOsc", pINFO) << "Generating decay...";
277 
279  LOG("NNBarOsc", pINFO) << "Decay product IDs: " << pdgv;
280  assert ( pdgv.size() > 1);
281 
282  LOG("NNBarOsc", pINFO) << "Performing a phase space decay...";
283 
284  // Get the decay product masses
285 
286  vector<int>::const_iterator pdg_iter;
287  int idx = 0;
288  double * mass = new double[pdgv.size()];
289  double sum = 0;
290  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
291  int pdgc = *pdg_iter;
292  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
293  mass[idx++] = m;
294  sum += m;
295  }
296 
297  LOG("NNBarOsc", pINFO)
298  << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
299  int initial_nucleus_id = 0;
300  int oscillating_neutron_id = 1;
301  int annihilation_nucleon_id = 2;
302 
303  // get our annihilating nucleons
304  GHepParticle * initial_nucleus = event->Particle(initial_nucleus_id);
305  assert(initial_nucleus);
306  GHepParticle * oscillating_neutron = event->Particle(oscillating_neutron_id);
307  assert(oscillating_neutron);
308  GHepParticle * annihilation_nucleon = event->Particle(annihilation_nucleon_id);
309  assert(annihilation_nucleon);
310 
311  Target tgt(initial_nucleus->Pdg());
313 
314  // get their momentum 4-vectors and boost into rest frame
315  TLorentzVector * p4_1 = oscillating_neutron->GetP4();
316  TLorentzVector * p4_2 = annihilation_nucleon->GetP4();
317  TLorentzVector * p4d = new TLorentzVector(*p4_1 + *p4_2);
318  TVector3 boost = p4d->BoostVector();
319  p4d->Boost(-boost);
320 
321  // get decay position
322  TLorentzVector * v4d = annihilation_nucleon->GetX4();
323 
324  delete p4_1;
325  delete p4_2;
326 
327  LOG("NNBarOsc", pINFO)
328  << "Decaying system p4 = " << utils::print::P4AsString(p4d);
329 
330  // Set the decay
331  bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
332 
333  // If the decay is not energetically allowed, select a new final state
334  while(!permitted) {
335 
336  LOG("NNBarOsc", pINFO)
337  << "Not enough energy to generate decay products! Selecting a new final state.";
338 
339  int mode = 0;
340 
341  int initial_nucleus_pdg = initial_nucleus->Pdg();
342 
343  std::string nucleus_pdg = std::to_string(static_cast<long long>(initial_nucleus_pdg));
344  if (nucleus_pdg.size() != 10) {
345  LOG("NNBarOsc", pERROR)
346  << "Expecting the nuclear PDG code to be a 10-digit integer, but it is " << nucleus_pdg << ", which is "
347  << nucleus_pdg.size() << " digits long. Drop me an email at jezhewes@gmail.com ; exiting...";
348  exit(1);
349  }
350 
351  int n_nucleons = std::stoi(nucleus_pdg.substr(6,3)) - 1;
352  int n_protons = std::stoi(nucleus_pdg.substr(3,3));
353  double proton_frac = ((double)n_protons) / ((double) n_nucleons);
354  double neutron_frac = 1 - proton_frac;
355 
356  // set branching ratios, taken from bubble chamber data
357  const int n_modes = 16;
358  double br [n_modes] = { 0.010, 0.080, 0.100, 0.220,
359  0.360, 0.160, 0.070, 0.020,
360  0.015, 0.065, 0.110, 0.280,
361  0.070, 0.240, 0.100, 0.100 };
362 
363  for (int i = 0; i < n_modes; i++) {
364  // for first 7 branching ratios, multiply by relative proton density
365  if (i < 7)
366  br[i] *= proton_frac;
367  // for next 9, multiply by relative neutron density
368  else
369  br[i] *= neutron_frac;
370  }
371 
372  // randomly generate a number between 1 and 0
374  rnd->SetSeed(0);
375  double p = rnd->RndNum().Rndm();
376 
377  // loop through all modes, figure out which one our random number corresponds to
378  double threshold = 0;
379  for (int j = 0; j < n_modes; j++) {
380  threshold += br[j];
381  if (p < threshold) {
382  // once we've found our mode, stop looping
383  mode = j + 1;
384  break;
385  }
386  }
387 
388  // create new event record with new final state
389  // TODO - I don't think Jeremy meant to make a _new_ record here; it
390  // shadows the one passed in...
391  // EventRecord * event = new EventRecord;
392  Interaction * interaction = Interaction::NOsc((int)fCurrInitStatePdg, mode);
393  event->AttachSummary(interaction);
394 
395  fCurrDecayMode = (NNBarOscMode_t) interaction->ExclTag().DecayMode();
396 
398  LOG("NNBarOsc", pINFO) << "Decay product IDs: " << pdgv;
399  assert ( pdgv.size() > 1);
400 
401  // get the decay particles again
402  LOG("NNBarOsc", pINFO) << "Performing a phase space decay...";
403  idx = 0;
404  delete [] mass;
405  mass = new double[pdgv.size()];
406  sum = 0;
407  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
408  int pdgc = *pdg_iter;
409  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
410  mass[idx++] = m;
411  sum += m;
412  }
413 
414  LOG("NNBarOsc", pINFO)
415  << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
416  LOG("NNBarOsc", pINFO)
417  << "Decaying system p4 = " << utils::print::P4AsString(p4d);
418 
419  permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
420  }
421 
422  if(!permitted) {
423  LOG("NNBarOsc", pERROR)
424  << " *** Phase space decay is not permitted \n"
425  << " Total particle mass = " << sum << "\n"
426  << " Decaying system p4 = " << utils::print::P4AsString(p4d);
427  // clean-up
428  delete [] mass;
429  delete p4d;
430  delete v4d;
431  // throw exception
433  exception.SetReason("Decay not permitted kinematically");
434  exception.SwitchOnFastForward();
435  throw exception;
436  }
437 
438  // Get the maximum weight
439  //double wmax = fPhaseSpaceGenerator.GetWtMax();
440  double wmax = -1;
441  for(int i=0; i<200; i++) {
442  double w = fPhaseSpaceGenerator.Generate();
443  wmax = TMath::Max(wmax,w);
444  }
445  assert(wmax>0);
446  wmax *= 2;
447 
448  LOG("NNBarOsc", pNOTICE)
449  << "Max phase space gen. weight @ current hadronic system: " << wmax;
450 
451  // Generate an unweighted decay
453 
454  bool accept_decay=false;
455  unsigned int itry=0;
456  while(!accept_decay)
457  {
458  itry++;
459 
461  // report, clean-up and return
462  LOG("NNBarOsc", pWARN)
463  << "Couldn't generate an unweighted phase space decay after "
464  << itry << " attempts";
465  // clean up
466  delete [] mass;
467  delete p4d;
468  delete v4d;
469  // throw exception
471  exception.SetReason("Couldn't select decay after N attempts");
472  exception.SwitchOnFastForward();
473  throw exception;
474  }
475  double w = fPhaseSpaceGenerator.Generate();
476  if(w > wmax) {
477  LOG("NNBarOsc", pWARN)
478  << "Decay weight = " << w << " > max decay weight = " << wmax;
479  }
480  double gw = wmax * rnd->RndHadro().Rndm();
481  accept_decay = (gw<=w);
482 
483  LOG("NNBarOsc", pINFO)
484  << "Decay weight = " << w << " / R = " << gw
485  << " - accepted: " << accept_decay;
486 
487  } //!accept_decay
488 
489  // Insert final state products into a TClonesArray of TMCParticles
490  TLorentzVector v4(*v4d);
491  int idp = 0;
492  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
493  int pdgc = *pdg_iter;
494  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
495  GHepStatus_t ist =
497  p4fin->Boost(boost);
498  event->AddParticle(pdgc, ist, oscillating_neutron_id,-1,-1,-1, *p4fin, v4);
499  idp++;
500  }
501 
502  // Clean-up
503  delete [] mass;
504  delete p4d;
505  delete v4d;
506 }
507 //___________________________________________________________________________
509 {
510  Algorithm::Configure(config);
511  this->LoadConfig();
512 }
513 //___________________________________________________________________________
515 {
516  Algorithm::Configure(config);
517  this->LoadConfig();
518 }
519 //___________________________________________________________________________
521 {
522 // AlgConfigPool * confp = AlgConfigPool::Instance();
523 // const Registry * gc = confp->GlobalParameterList();
524 
525  fNuclModel = 0;
526 
527  RgKey nuclkey = "NuclearModel";
528  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
529  assert(fNuclModel);
530 }
531 //___________________________________________________________________________
532 
TLorentzVector * GetX4(void) const
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
double RemovalEnergy(void) const
Definition: NuclearModelI.h:57
int AnnihilatingNucleonPdgCode(NNBarOscMode_t ndm)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
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)
enum genie::EGHepStatus GHepStatus_t
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
int Pdg(void) const
Definition: Target.h:72
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
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:42
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void SetMomentum(const TLorentzVector &p4)
Definition: config.py:1
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:30
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
A list of PDG codes.
Definition: PDGCodeList.h:33
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
Double_t ymax
Definition: plot.C:25
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...
const double R0
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
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:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
double energy
Definition: plottest35.C:25
void AddInitialState(GHepRecord *event) const
int DecayMode(void) const
Definition: XclsTag.h:63
#define R(x)
string AsString(NNBarOscMode_t ndm)
const double j
Definition: BetheBloch.cxx:29
#define pINFO
Definition: Messenger.h:63
#define pWARN
Definition: Messenger.h:61
TRandom3 & RndNum(void) const
rnd number generator used by MC integrators & other numerical methods
Definition: RandomGen.h:78
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:54
static const double A
Definition: Units.h:82
PDGCodeList DecayProductList(NNBarOscMode_t ndm)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:188
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
void ProcessEventRecord(GHepRecord *event) const
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
exit(0)
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
virtual bool GenerateNucleon(const Target &) const =0
const InitialState & InitState(void) const
Definition: Interaction.h:69
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:53
const int kPdgProton
Definition: PDGCodes.h:65
void GenerateOscillatingNeutronPosition(GHepRecord *event) const
int A(void) const
#define pNOTICE
Definition: Messenger.h:62
const Target & Tgt(void) const
Definition: InitialState.h:67
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
Double_t sum
Definition: plot.C:31
Float_t w
Definition: plot.C:20
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
static const double kPi
Definition: Constants.h:38
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:40
void GenerateFermiMomentum(GHepRecord *event) const
static Interaction * NOsc(int tgt, int annihilation_mode=-1)
void GenerateDecayProducts(GHepRecord *event) const
void SetSeed(long int seed)
Definition: RandomGen.cxx:90
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353