GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NucleonDecayPrimaryVtxGenerator.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
30 
31 using namespace genie;
32 
33 //____________________________________________________________________________
35 EventRecordVisitorI("genie::NucleonDecayPrimaryVtxGenerator")
36 {
37 
38 }
39 //____________________________________________________________________________
41  string config) :
42 EventRecordVisitorI("genie::NucleonDecayPrimaryVtxGenerator",config)
43 {
44 
45 }
46 //____________________________________________________________________________
48 {
49 
50 }
51 //____________________________________________________________________________
53  GHepRecord * event) const
54 {
55  Interaction * interaction = event->Summary();
56  fCurrInitStatePdg = interaction->InitState().Tgt().Pdg();
58  fCurrDecayedNucleon = interaction->InitState().Tgt().HitNucPdg();
59 
60  LOG("NucleonDecay", pNOTICE)
62  << " for an initial state with code: " << fCurrInitStatePdg;
63 
65 
66  this->AddInitialState(event);
67  this->GenerateDecayedNucleonPosition(event);
68  this->GenerateFermiMomentum(event);
69  this->GenerateDecayProducts(event);
70 }
71 //____________________________________________________________________________
73  GHepRecord * event) const
74 {
75 //
76 // Add initial state in the event record.
77 //
78 // If the decayed nucleon is one bound in a nucleus, the event record is
79 // initialized as follows:
80 // id: 0, nucleus (kIStInitialState)
81 // |
82 // |---> { id: 1, nucleon (kIStDecayedState)
83 // { id: 2, remnant nucleus (kIStStableFinalState)
84 //
85 // If the decayed nucleon is a free one, the event record is initialized as
86 // follows:
87 // id: 0, nucleon (kIStInitialState)
88 // |
89 // |---> id: 1, nucleon (kIStDecayedState)
90 //
91 
92  TLorentzVector v4(0,0,0,0);
93 
97 
98  int ipdg = fCurrInitStatePdg;
99 
100  // Decayed nucleon is a bound one.
101  if(fNucleonIsBound)
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 decayed nucleon
109  int dpdg = fCurrDecayedNucleon;
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--;
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 
125  // Decayed nucleon is a free one
126  else
127  {
128  // Initial state is either a neutron or a proton.
129  // Convert the initial state PDG code from the ion convention (10LZZZAAAI)
130  // to the usual code for neutrons or protons
131  int ipdg_short = -1;
132  if(ipdg == kPdgTgtFreeP) ipdg_short = kPdgProton;
133  if(ipdg == kPdgTgtFreeN) ipdg_short = kPdgNeutron;
134 
135  // Decayed nucleon code
136  int dpdg = fCurrDecayedNucleon;
137 
138  if(dpdg != ipdg_short) {
139  LOG("NucleonDecay", pWARN)
140  << "Couldn't generate given decay ("
142  << " for given initial state (PDG = " << ipdg_short << ")";
144  exception.SetReason("Decay-mode / Initial-state mismatch!");
145  exception.SwitchOnFastForward();
146  throw exception;
147  }
148  // add initial nucleon
149  double mn = PDGLibrary::Instance()->Find(ipdg)->Mass();
150  TLorentzVector p4i(0,0,0,mn);
151  event->AddParticle(dpdg,stis,-1,-1,-1,-1, p4i, v4);
152  // add decayed nucleon
153  event->AddParticle(dpdg,stdc,0,-1,-1,-1, p4i, v4);
154  }
155 }
156 //____________________________________________________________________________
158  GHepRecord * event) const
159 {
160  GHepParticle * initial_nucleus = event->Particle(0);
161  int A = initial_nucleus->A();
162  if(A <= 2) {
163  return;
164  }
165 
166  RandomGen * rnd = RandomGen::Instance();
167 
168  double R0 = 1.3;
169  double dA = (double)A;
170  double R = R0 * TMath::Power(dA, 1./3.);
171 
172  LOG("NucleonDecay", pINFO)
173  << "Generating vertex according to a realistic nuclear density profile";
174 
175  // get inputs to the rejection method
176  double ymax = -1;
177  double rmax = 3*R;
178  double dr = R/40.;
179  for(double r = 0; r < rmax; r+=dr) {
180  ymax = TMath::Max(ymax, r*r * utils::nuclear::Density(r,A));
181  }
182  ymax *= 1.2;
183 
184  // select a vertex using the rejection method
185  TLorentzVector vtx(0,0,0,0);
186  unsigned int iter = 0;
187  while(1) {
188  iter++;
189 
190  // throw an exception if it hasn't find a solution after many attempts
191  if(iter > controls::kRjMaxIterations) {
192  LOG("NucleonDecay", pWARN)
193  << "Couldn't generate a vertex position after " << iter << " iterations";
195  exception.SetReason("Couldn't generate vertex");
196  exception.SwitchOnFastForward();
197  throw exception;
198  }
199 
200  double r = rmax * rnd->RndFsi().Rndm();
201  double t = ymax * rnd->RndFsi().Rndm();
202  double y = r*r * utils::nuclear::Density(r,A);
203  if(y > ymax) {
204  LOG("NucleonDecay", pERROR)
205  << "y = " << y << " > ymax = " << ymax << " for r = " << r << ", A = " << A;
206  }
207  bool accept = (t < y);
208  if(accept) {
209  double phi = 2*constants::kPi * rnd->RndFsi().Rndm();
210  double cosphi = TMath::Cos(phi);
211  double sinphi = TMath::Sin(phi);
212  double costheta = -1 + 2 * rnd->RndFsi().Rndm();
213  double sintheta = TMath::Sqrt(1-costheta*costheta);
214  vtx.SetX(r*sintheta*cosphi);
215  vtx.SetY(r*sintheta*sinphi);
216  vtx.SetZ(r*costheta);
217  vtx.SetT(0.);
218  break;
219  }
220  } // while 1
221 
222  GHepParticle * decayed_nucleon = event->Particle(1);
223  assert(decayed_nucleon);
224  decayed_nucleon->SetPosition(vtx);
225 }
226 //____________________________________________________________________________
228  GHepRecord * event) const
229 {
230  GHepParticle * initial_nucleus = event->Particle(0);
231  int A = initial_nucleus->A();
232  if(A <= 2) {
233  return;
234  }
235 
236  GHepParticle * decayed_nucleon = event->Particle(1);
237  GHepParticle * remnant_nucleus = event->Particle(2);
238  assert(decayed_nucleon);
239  assert(remnant_nucleus);
240 
241  // generate a Fermi momentum & removal energy
242  Target tgt(initial_nucleus->Pdg());
243  tgt.SetHitNucPdg(decayed_nucleon->Pdg());
245  TVector3 p3 = fNuclModel->Momentum3();
246  double w = fNuclModel->RemovalEnergy();
247 
248  LOG("FermiMover", pINFO)
249  << "Generated nucleon momentum: ("
250  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
251  << "|p| = " << p3.Mag();
252  LOG("NucleonDecay", pINFO)
253  << "Generated nucleon removal energy: w = " << w;
254 
255  double pF2 = p3.Mag2(); // (fermi momentum)^2
256 
257  double Mi = PDGLibrary::Instance()->Find(initial_nucleus->Pdg())-> Mass(); // initial nucleus mass
258  double Mf = PDGLibrary::Instance()->Find(remnant_nucleus->Pdg())-> Mass(); // remnant nucleus mass
259 
260  double EN = Mi - TMath::Sqrt(pF2 + Mf*Mf);
261 
262  TLorentzVector p4nucleon( p3.Px(), p3.Py(), p3.Pz(), EN);
263  TLorentzVector p4remnant(-1*p3.Px(), -1*p3.Py(), -1*p3.Pz(), Mi-EN);
264 
265  decayed_nucleon->SetMomentum(p4nucleon);
266  remnant_nucleus->SetMomentum(p4remnant);
267 }
268 //____________________________________________________________________________
270  GHepRecord * event) const
271 {
272  LOG("NucleonDecay", pINFO) << "Generating decay...";
273 
275  LOG("NucleonDecay", pINFO) << "Decay product IDs: " << pdgv;
276  assert ( pdgv.size() > 1);
277 
278  LOG("NucleonDecay", pINFO) << "Performing a phase space decay...";
279 
280  // Get the decay product masses
281 
282  vector<int>::const_iterator pdg_iter;
283  int i = 0;
284  double * mass = new double[pdgv.size()];
285  double sum = 0;
286  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
287  int pdgc = *pdg_iter;
288  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
289  mass[i++] = m;
290  sum += m;
291  }
292 
293  LOG("NucleonDecay", pINFO)
294  << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
295 
296  int decayed_nucleon_id = 1;
297  GHepParticle * decayed_nucleon = event->Particle(decayed_nucleon_id);
298  assert(decayed_nucleon);
299  TLorentzVector * p4d = decayed_nucleon->GetP4();
300  TLorentzVector * v4d = decayed_nucleon->GetX4();
301 
302  LOG("NucleonDecay", pINFO)
303  << "Decaying system p4 = " << utils::print::P4AsString(p4d);
304 
305  // Set the decay
306  bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
307  if(!permitted) {
308  LOG("NucleonDecay", pERROR)
309  << " *** Phase space decay is not permitted \n"
310  << " Total particle mass = " << sum << "\n"
311  << " Decaying system p4 = " << utils::print::P4AsString(p4d);
312  // clean-up
313  delete [] mass;
314  delete p4d;
315  delete v4d;
316  // throw exception
318  exception.SetReason("Decay not permitted kinematically");
319  exception.SwitchOnFastForward();
320  throw exception;
321  }
322 
323  // Get the maximum weight
324  //double wmax = fPhaseSpaceGenerator.GetWtMax();
325  double wmax = -1;
326  for(int idec=0; idec<200; idec++) {
327  double w = fPhaseSpaceGenerator.Generate();
328  wmax = TMath::Max(wmax,w);
329  }
330  assert(wmax>0);
331  wmax *= 2;
332 
333  LOG("NucleonDecay", pNOTICE)
334  << "Max phase space gen. weight @ current hadronic system: " << wmax;
335 
336  // Generate an unweighted decay
337  RandomGen * rnd = RandomGen::Instance();
338 
339  bool accept_decay=false;
340  unsigned int itry=0;
341  while(!accept_decay)
342  {
343  itry++;
344 
346  // report, clean-up and return
347  LOG("NucleonDecay", pWARN)
348  << "Couldn't generate an unweighted phase space decay after "
349  << itry << " attempts";
350  // clean up
351  delete [] mass;
352  delete p4d;
353  delete v4d;
354  // throw exception
356  exception.SetReason("Couldn't select decay after N attempts");
357  exception.SwitchOnFastForward();
358  throw exception;
359  }
360  double w = fPhaseSpaceGenerator.Generate();
361  if(w > wmax) {
362  LOG("NucleonDecay", pWARN)
363  << "Decay weight = " << w << " > max decay weight = " << wmax;
364  }
365  double gw = wmax * rnd->RndHadro().Rndm();
366  accept_decay = (gw<=w);
367 
368  LOG("NucleonDecay", pINFO)
369  << "Decay weight = " << w << " / R = " << gw
370  << " - accepted: " << accept_decay;
371 
372  } //!accept_decay
373 
374  // Insert final state products into a TClonesArray of GHepParticle's
375  TLorentzVector v4(*v4d);
376  int idp = 0;
377  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
378  int pdgc = *pdg_iter;
379  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
380  GHepStatus_t ist =
382  event->AddParticle(pdgc, ist, decayed_nucleon_id,-1,-1,-1, *p4fin, v4);
383  idp++;
384  }
385 
386  // Clean-up
387  delete [] mass;
388  delete p4d;
389  delete v4d;
390 }
391 //____________________________________________________________________________
393 {
394  Algorithm::Configure(config);
395  this->LoadConfig();
396 }
397 //___________________________________________________________________________
399 {
400  Algorithm::Configure(config);
401  this->LoadConfig();
402 }
403 //___________________________________________________________________________
405 {
406 // AlgConfigPool * confp = AlgConfigPool::Instance();
407 // const Registry * gc = confp->GlobalParameterList();
408 
409  fNuclModel = 0;
410 
411  RgKey nuclkey = "NuclearModel";
412  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
413  assert(fNuclModel);
414 }
415 //___________________________________________________________________________
TLorentzVector * GetX4(void) const
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
void GenerateDecayedNucleonPosition(GHepRecord *event) const
double RemovalEnergy(void) const
Definition: NuclearModelI.h:65
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
int HitNucPdg(void) const
Definition: Target.cxx:304
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double Density(double r, int A, double ring=0.)
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
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
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
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
const int kPdgTgtFreeN
Definition: PDGCodes.h:200
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
int DecayMode(void) const
Definition: XclsTag.h:70
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
PDGCodeList DecayProductList(NucleonDecayMode_t ndm, int npdg=0)
GHepStatus_t DecayProductStatus(bool in_nucleus, int pdgc)
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
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
enum genie::ENucleonDecayMode NucleonDecayMode_t
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
string AsString(NucleonDecayMode_t ndm, int npdg=0)
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
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
enum genie::EGHepStatus GHepStatus_t
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345