FermiMover.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 - October 08, 2004
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 07, 2009 - CA
14  Added option to simulate the recoil nucleon due to short range corelation.
15  The recoil nucleon is added in case that the selected hit nucleon has a
16  momentum selected from the NN corelation tail. The 2 nucleons are put
17  back-to-back. For the time-being using hard-coded relative fractions for
18  the nn, pp, np pairs.
19  The code for adding the recoil nuclear target at the GHEP record was moved
20  into this processing step.
21 
22  @ Mar 18, 2016- Joe Johnston (SD)
23  Call GenerateNucleon() with a target and a radius, so the local Fermi
24  gas model can access the radius.
25  Added a check to see if a local Fermi gas model is being used. If so,
26  use a local Fermi gas model when deciding whether to eject a recoil nucleon.
27 */
28 //____________________________________________________________________________
29 
30 #include <cstdlib>
31 
32 #include <TLorentzVector.h>
33 #include <TVector3.h>
34 #include <TParticlePDG.h>
35 #include <TMath.h>
36 
42 
60 
61 using namespace genie;
62 using namespace genie::constants;
63 
64 //___________________________________________________________________________
66 EventRecordVisitorI("genie::FermiMover")
67 {
68 
69 }
70 //___________________________________________________________________________
72 EventRecordVisitorI("genie::FermiMover", config)
73 {
74 
75 }
76 //___________________________________________________________________________
78 {
79 
80 }
81 //___________________________________________________________________________
83 {
84  // skip if not a nuclear target
85  if(! evrec->Summary()->InitState().Tgt().IsNucleus()) return;
86 
87  // skip if no hit nucleon is set
88  if(! evrec->HitNucleon()) return;
89 
90  // give hit nucleon a Fermi momentum
91  this->KickHitNucleon(evrec);
92 
93  // add a recoiled nucleus remnant
94  this->AddTargetNucleusRemnant(evrec);
95 }
96 //___________________________________________________________________________
98 {
99  Interaction * interaction = evrec -> Summary();
100  InitialState * init_state = interaction -> InitStatePtr();
101  Target * tgt = init_state -> TgtPtr();
102 
103  // do nothing for non-nuclear targets
104  if(!tgt->IsNucleus()) return;
105 
106  TLorentzVector * p4 = tgt->HitNucP4Ptr();
107 
108  // do nothing if the struct nucleon 4-momentum was set (eg as part of the
109  // initial state selection)
110  if(p4->Px()>0 || p4->Py()>0 || p4->Pz()>0) return;
111 
112  // access the hit nucleon and target nucleus at the GHEP record
113  GHepParticle * nucleon = evrec->HitNucleon();
114  GHepParticle * nucleus = evrec->TargetNucleus();
115  assert(nucleon);
116  assert(nucleus);
117 
118  // generate a Fermi momentum & removal energy
119  // call GenerateNucleon with a radius in case the model is LocalFGM
120  double rad = nucleon->X4()->Vect().Mag();
121  fNuclModel->GenerateNucleon(*tgt,rad);
122 
123  TVector3 p3 = fNuclModel->Momentum3();
124  double w = fNuclModel->RemovalEnergy();
125 
126  LOG("FermiMover", pINFO)
127  << "Generated nucleon momentum: ("
128  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
129  << "|p| = " << p3.Mag();
130  LOG("FermiMover", pINFO)
131  << "Generated nucleon removal energy: w = " << w;
132 
133  double pF2 = p3.Mag2(); // (fermi momentum)^2
134 
135  nucleon->SetRemovalEnergy(w);
136 
137  // struck nucleon energy:
138  // two possible prescriptions depending on whether you want to force
139  // the sruck nucleon to be on the mass-shell or not...
140 
141  double EN=0;
142  // Set this to either a proton or neutron to eject a secondary particle
143  int eject_nucleon_pdg = 0;
146  if (interaction_type == kFermiMoveEffectiveSF1p1h) {
147  EN = nucleon->Mass() - w -
148  pF2 / (2 * (nucleus->Mass() - nucleon->Mass()));
149  } else if (interaction_type == kFermiMoveEffectiveSF2p2h_eject ||
150  interaction_type == kFermiMoveEffectiveSF2p2h_noeject) {
151  TParticlePDG * other_nucleon;
152  if(nucleon->Pdg() == kPdgProton) {
153  other_nucleon = PDGLibrary::Instance()->Find(kPdgNeutron);
154  }
155  else {
156  other_nucleon = PDGLibrary::Instance()->Find(kPdgProton);
157  }
158  TParticlePDG * deuteron = PDGLibrary::Instance()->Find(1000010020);
159  EN = deuteron->Mass() - 2 * w -
160  TMath::Sqrt(pF2 + other_nucleon->Mass() * other_nucleon->Mass());
161  if (interaction_type == kFermiMoveEffectiveSF2p2h_eject &&
162  deuteron->PdgCode() != nucleus->Pdg()) {
163  // Eject a remnant nucleon for the 2p2h process. Don't eject other
164  // nucleon if target is deuterium- need to leave behind a
165  // valid remnant nucleus.
166  eject_nucleon_pdg = other_nucleon->PdgCode();
167  }
168  // Do default Fermi Moving
169  } else {
170  if (!fKeepNuclOnMassShell) {
171  //-- compute A,Z for final state nucleus & get its PDG code
172  int nucleon_pdgc = nucleon->Pdg();
173  bool is_p = pdg::IsProton(nucleon_pdgc);
174  int Z = (is_p) ? nucleus->Z()-1 : nucleus->Z();
175  int A = nucleus->A() - 1;
176 
177  TParticlePDG * fnucleus = 0;
178  int ipdgc = pdg::IonPdgCode(A, Z);
179  fnucleus = PDGLibrary::Instance()->Find(ipdgc);
180  if(!fnucleus) {
181  LOG("FermiMover", pFATAL)
182  << "No particle with [A = " << A << ", Z = " << Z
183  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
184  exit(1);
185  }
186  //-- compute the energy of the struck (off the mass-shell) nucleus
187 
188  double Mf = fnucleus -> Mass(); // remnant nucleus mass
189  double Mi = nucleus -> Mass(); // initial nucleus mass
190 
191  EN = Mi - TMath::Sqrt(pF2 + Mf*Mf);
192  } else {
193  double MN = nucleon->Mass();
194  double MN2 = TMath::Power(MN,2);
195  EN = TMath::Sqrt(MN2+pF2);
196  }
197  if (fSRCRecoilNucleon) {
198  const int nucleus_pdgc = nucleus->Pdg();
199  const int nucleon_pdgc = nucleon->Pdg();
200 
201  // Calculate the Fermi momentum, using a local Fermi gas if the
202  // nuclear model is LocalFGM, and RFG otherwise
203  double kF;
205  assert(pdg::IsProton(nucleon_pdgc) || pdg::IsNeutron(nucleon_pdgc));
206  int A = tgt->A();
207  bool is_p = pdg::IsProton(nucleon_pdgc);
208  double numNuc = (is_p) ? (double)tgt->Z():(double)tgt->N();
209  double radius = nucleon->X4()->Vect().Mag();
211  kF= TMath::Power(3*kPi2*numNuc*
212  genie::utils::nuclear::Density(radius,A),1.0/3.0) *hbarc;
213  }else{
215  const FermiMomentumTable * kft = kftp->GetTable("Default");
216  kF = kft->FindClosestKF(nucleus_pdgc, nucleon_pdgc);
217  }
218  if (TMath::Sqrt(pF2) > kF) {
219  double Pp = (nucleon->Pdg() == kPdgProton) ? 0.05 : 0.95;
221  double prob = rnd->RndGen().Rndm();
222  eject_nucleon_pdg = (prob > Pp) ? kPdgNeutron : kPdgProton;
223  }
224  }
225  }
226 
227  //-- update the struck nucleon 4p at the interaction summary and at
228  // the GHEP record
229  p4->SetPx( p3.Px() );
230  p4->SetPy( p3.Py() );
231  p4->SetPz( p3.Pz() );
232  p4->SetE ( EN );
233 
234  nucleon->SetMomentum(*p4); // update GHEP value
235 
236  // Sometimes, for interactions near threshold, Fermi momentum might bring
237  // the neutrino energy in the nucleon rest frame below threshold (for the
238  // selected interaction). In this case mark the event as unphysical and
239  // abort the current thread.
240  const KPhaseSpace & kps = interaction->PhaseSpace();
241  if(!kps.IsAboveThreshold()) {
242  LOG("FermiMover", pNOTICE)
243  << "Event below threshold after generating Fermi momentum";
244 
245  double Ethr = kps.Threshold();
246  double Ev = init_state->ProbeE(kRfHitNucRest);
247  LOG("FermiMover", pNOTICE)
248  << "Ev (@ nucleon rest frame) = " << Ev << ", Ethr = " << Ethr;
249 
250  evrec->EventFlags()->SetBitNumber(kBelowThrNRF, true);
252  exception.SetReason("E < Ethr after generating nucleon Fermi momentum");
253  exception.SwitchOnFastForward();
254  throw exception;
255  }
256  if(eject_nucleon_pdg != 0) {
257  this->Emit2ndNucleonFromSRC(evrec, eject_nucleon_pdg);
258  }
259 }
260 //___________________________________________________________________________
262  const int eject_pdg_code) const
263 {
264  LOG("FermiMover", pINFO) << "Adding a recoil nucleon "
265  "due to short range corelation";
266 
267  GHepParticle * nucleon = evrec->HitNucleon();
268 
270  int imom = evrec->TargetNucleusPosition();
271 
272  //-- Has opposite momentum from the struck nucleon
273  double vx = nucleon->Vx();
274  double vy = nucleon->Vy();
275  double vz = nucleon->Vz();
276  double px = -1.* nucleon->Px();
277  double py = -1.* nucleon->Py();
278  double pz = -1.* nucleon->Pz();
279  double M = PDGLibrary::Instance()->Find(eject_pdg_code)->Mass();
280  double E = TMath::Sqrt(px*px+py*py+pz*pz+M*M);
281 
282  evrec->AddParticle(
283  eject_pdg_code, status, imom, -1, -1, -1, px, py, pz, E, vx, vy, vz, 0);
284 }
285 //___________________________________________________________________________
287 {
288 // add the remnant nuclear target at the GHEP record
289 
290  LOG("FermiMover", pINFO) << "Adding final state nucleus";
291 
292  double Px = 0;
293  double Py = 0;
294  double Pz = 0;
295  double E = 0;
296 
297  GHepParticle * nucleus = evrec->TargetNucleus();
298  int A = nucleus->A();
299  int Z = nucleus->Z();
300 
301  int fd = nucleus->FirstDaughter();
302  int ld = nucleus->LastDaughter();
303 
304  for(int id = fd; id <= ld; id++) {
305 
306  // compute A,Z for final state nucleus & get its PDG code and its mass
307  GHepParticle * particle = evrec->Particle(id);
308  assert(particle);
309  int pdgc = particle->Pdg();
310  bool is_p = pdg::IsProton (pdgc);
311  bool is_n = pdg::IsNeutron(pdgc);
312 
313  if (is_p) Z--;
314  if (is_p || is_n) A--;
315 
316  Px += particle->Px();
317  Py += particle->Py();
318  Pz += particle->Pz();
319  E += particle->E();
320 
321  }//daughters
322 
323  TParticlePDG * remn = 0;
324  int ipdgc = pdg::IonPdgCode(A, Z);
325  remn = PDGLibrary::Instance()->Find(ipdgc);
326  if(!remn) {
327  LOG("HadronicVtx", pFATAL)
328  << "No particle with [A = " << A << ", Z = " << Z
329  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
330  assert(remn);
331  }
332 
333  double Mi = nucleus->Mass();
334  Px *= -1;
335  Py *= -1;
336  Pz *= -1;
337  E = Mi-E;
338 
339  // Add the nucleus to the event record
340  LOG("FermiMover", pINFO)
341  << "Adding nucleus [A = " << A << ", Z = " << Z
342  << ", pdgc = " << ipdgc << "]";
343 
344  int imom = evrec->TargetNucleusPosition();
345  evrec->AddParticle(
346  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
347 }
348 //___________________________________________________________________________
350 {
351  Algorithm::Configure(config);
352  this->LoadConfig();
353 }
354 //____________________________________________________________________________
356 {
357  Algorithm::Configure(config);
358  this->LoadConfig();
359 }
360 //____________________________________________________________________________
362 {
363  RgKey nuclkey = "NuclearModel";
364  fNuclModel = 0;
365  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
367 
368  this->GetParamDef("KeepHitNuclOnMassShell", fKeepNuclOnMassShell, false);
369  this->GetParamDef("SimRecoilNucleon", fSRCRecoilNucleon, false);
370 }
371 //____________________________________________________________________________
bool fSRCRecoilNucleon
simulate recoil nucleon due to short range corellation?
Definition: FermiMover.h:56
int Z(void) const
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
void Configure(const Registry &config)
Definition: FermiMover.cxx:349
int status
Definition: fabricate.py:1613
double RemovalEnergy(void) const
Definition: NuclearModelI.h:57
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double E(void) const
Get energy.
Definition: GHepParticle.h:92
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
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.)
int FirstDaughter(void) const
Definition: GHepParticle.h:69
static const double fermi
Definition: Units.h:63
int A(void) const
Definition: Target.h:71
double Threshold(void) const
Energy threshold.
Definition: KPhaseSpace.cxx:81
enum genie::EGHepStatus GHepStatus_t
#define pFATAL
Definition: Messenger.h:57
bool IsNucleus(void) const
Definition: Target.cxx:289
void ProcessEventRecord(GHepRecord *event_rec) const
Definition: FermiMover.cxx:82
static FermiMomentumTablePool * Instance(void)
void KickHitNucleon(GHepRecord *evrec) const
give hit nucleon a momentum
Definition: FermiMover.cxx:97
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)
A table of Fermi momentum constants.
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
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
FermiMoverInteractionType_t GetFermiMoverInteractionType(void) const
Definition: NuclearModelI.h:72
interaction_type
Definition: Constants.h:104
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
virtual NuclearModel_t ModelType(const Target &) const =0
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
Summary information for an interaction.
Definition: Interaction.h:56
int LastDaughter(void) const
Definition: GHepParticle.h:70
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
const FermiMomentumTable * GetTable(string name)
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
Kinematical phase space.
Definition: KPhaseSpace.h:34
const NuclearModelI * fNuclModel
nuclear model
Definition: FermiMover.h:57
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
int Z(void) const
Definition: Target.h:69
Double_t radius
#define pINFO
Definition: Messenger.h:63
static constexpr Double_t rad
Definition: Munits.h:162
void LoadConfig(void)
Definition: FermiMover.cxx:361
static const double kLightSpeed
Definition: Constants.h:32
enum genie::EFermiMoverInteractionType FermiMoverInteractionType_t
void SetRemovalEnergy(double Erm)
static const double A
Definition: Units.h:82
int N(void) const
Definition: Target.h:70
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:118
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
double Vz(void) const
Get production z.
Definition: GHepParticle.h:97
string RgKey
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:350
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:81
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
exit(0)
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
assert(nhit_max >=nhit_nbins)
virtual bool GenerateNucleon(const Target &) const =0
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
const InitialState & InitState(void) const
Definition: Interaction.h:69
const int kPdgProton
Definition: PDGCodes.h:65
void AddTargetNucleusRemnant(GHepRecord *evrec) const
^ emit a 2nd nucleon due to short range corellations
Definition: FermiMover.cxx:286
double Vy(void) const
Get production y.
Definition: GHepParticle.h:96
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
int A(void) const
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:67
bool fKeepNuclOnMassShell
keep hit bound nucleon on the mass shell?
Definition: FermiMover.h:55
Float_t w
Definition: plot.C:20
static const double kPlankConstant
Definition: Constants.h:33
void Emit2ndNucleonFromSRC(GHepRecord *evrec, const int eject_nucleon_pdg) const
Definition: FermiMover.cxx:261
double ProbeE(RefFrame_t rf) const
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
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
static const double kPi2
Definition: Constants.h:39
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:406
double Vx(void) const
Get production x.
Definition: GHepParticle.h:95
Initial State information.
Definition: InitialState.h:49
double Py(void) const
Get Py.
Definition: GHepParticle.h:90
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353