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