GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FermiMover.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 
7  Costas Andreopoulos <c.andreopoulos \at cern.ch>
8  University of Liverpool - 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 //___________________________________________________________________________
71 FermiMover::FermiMover(string config) :
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  // handle the addition of the recoil nucleon
95 
96  // add a recoiled nucleus remnant
97  this->AddTargetNucleusRemnant(evrec);
98 }
99 //___________________________________________________________________________
101 {
102  Interaction * interaction = evrec -> Summary();
103  InitialState * init_state = interaction -> InitStatePtr();
104  Target * tgt = init_state -> TgtPtr();
105 
106  // do nothing for non-nuclear targets
107  if(!tgt->IsNucleus()) return;
108 
109  TLorentzVector * p4 = tgt->HitNucP4Ptr();
110 
111  // do nothing if the struct nucleon 4-momentum was set (eg as part of the
112  // initial state selection)
113  if(p4->Px()>0 || p4->Py()>0 || p4->Pz()>0) return;
114 
115  // access the hit nucleon and target nucleus at the GHEP record
116  GHepParticle * nucleon = evrec->HitNucleon();
117  GHepParticle * nucleus = evrec->TargetNucleus();
118  assert(nucleon);
119  assert(nucleus);
120 
121  // generate a Fermi momentum & removal energy
122  // call GenerateNucleon with a radius in case the model is LocalFGM
123  double rad = nucleon->X4()->Vect().Mag();
124  fNuclModel->GenerateNucleon(*tgt,rad);
125 
126  TVector3 p3 = fNuclModel->Momentum3();
127  double w = fNuclModel->RemovalEnergy();
128 
129  LOG("FermiMover", pINFO)
130  << "Generated nucleon momentum: ("
131  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
132  << "|p| = " << p3.Mag();
133  LOG("FermiMover", pINFO)
134  << "Generated nucleon removal energy: w = " << w;
135 
136  double pF2 = p3.Mag2(); // (fermi momentum)^2
137 
138  nucleon->SetRemovalEnergy(w);
139 
140  // struck nucleon energy:
141  // two possible prescriptions depending on whether you want to force
142  // the sruck nucleon to be on the mass-shell or not...
143 
144  double EN=0;
146 
147  // EffectiveSF treatment or momentum-dependent removal energy
148  if (interaction_type == kFermiMoveEffectiveSF1p1h || fMomDepErmv ) {
149  EN = nucleon->Mass() - w - pF2 / (2 * (nucleus->Mass() - nucleon->Mass()));
150  } else if (interaction_type == kFermiMoveEffectiveSF2p2h_eject ||
151  interaction_type == kFermiMoveEffectiveSF2p2h_noeject) {
152 
153  int other_nucleon_pdg = nucleon->Pdg() == kPdgProton ? kPdgNeutron : kPdgProton ;
154  TParticlePDG * other_nucleon = PDGLibrary::Instance()->Find( other_nucleon_pdg );
155 
156  TParticlePDG * deuteron = PDGLibrary::Instance()->Find(1000010020);
157  EN = deuteron->Mass() - 2 * w - TMath::Sqrt(pF2 + other_nucleon->Mass() * other_nucleon->Mass());
158 
159  // Do default Fermi Moving
160  } else {
161  if (!fKeepNuclOnMassShell) {
162  //-- compute A,Z for final state nucleus & get its PDG code
163  int nucleon_pdgc = nucleon->Pdg();
164  bool is_p = pdg::IsProton(nucleon_pdgc);
165  int Z = (is_p) ? nucleus->Z()-1 : nucleus->Z();
166  int A = nucleus->A() - 1;
167 
168  TParticlePDG * fnucleus = 0;
169  int ipdgc = pdg::IonPdgCode(A, Z);
170  fnucleus = PDGLibrary::Instance()->Find(ipdgc);
171  if(!fnucleus) {
172  LOG("FermiMover", pFATAL)
173  << "No particle with [A = " << A << ", Z = " << Z
174  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
175  exit(1);
176  }
177  //-- compute the energy of the struck (off the mass-shell) nucleus
178 
179  double Mf = fnucleus -> Mass(); // remnant nucleus mass
180  double Mi = nucleus -> Mass(); // initial nucleus mass
181 
182  EN = Mi - TMath::Sqrt(pF2 + Mf*Mf);
183  } else {
184  double MN = nucleon->Mass();
185  double MN2 = TMath::Power(MN,2);
186  EN = TMath::Sqrt(MN2+pF2);
187  }
188 
189  }
190 
191  //-- update the struck nucleon 4p at the interaction summary and at
192  // the GHEP record
193  p4->SetPx( p3.Px() );
194  p4->SetPy( p3.Py() );
195  p4->SetPz( p3.Pz() );
196  p4->SetE ( EN );
197 
198  nucleon->SetMomentum(*p4); // update GHEP value
199 
200  // Sometimes, for interactions near threshold, Fermi momentum might bring
201  // the neutrino energy in the nucleon rest frame below threshold (for the
202  // selected interaction). In this case mark the event as unphysical and
203  // abort the current thread.
204  const KPhaseSpace & kps = interaction->PhaseSpace();
205  if(!kps.IsAboveThreshold()) {
206  LOG("FermiMover", pNOTICE)
207  << "Event below threshold after generating Fermi momentum";
208 
209  double Ethr = kps.Threshold();
210  double Ev = init_state->ProbeE(kRfHitNucRest);
211  LOG("FermiMover", pNOTICE)
212  << "Ev (@ nucleon rest frame) = " << Ev << ", Ethr = " << Ethr;
213 
214  evrec->EventFlags()->SetBitNumber(kBelowThrNRF, true);
216  exception.SetReason("E < Ethr after generating nucleon Fermi momentum");
217  exception.SwitchOnFastForward();
218  throw exception;
219  }
220 
221 
222 }
223 //___________________________________________________________________________
225 {
226 // add the remnant nuclear target at the GHEP record
227 
228  LOG("FermiMover", pINFO) << "Adding final state nucleus";
229 
230  double Px = 0;
231  double Py = 0;
232  double Pz = 0;
233  double E = 0;
234 
235  GHepParticle * nucleus = evrec->TargetNucleus();
236  int A = nucleus->A();
237  int Z = nucleus->Z();
238 
239  int fd = nucleus->FirstDaughter();
240  int ld = nucleus->LastDaughter();
241 
242  for(int id = fd; id <= ld; id++) {
243 
244  // compute A,Z for final state nucleus & get its PDG code and its mass
245  GHepParticle * particle = evrec->Particle(id);
246  assert(particle);
247  int pdgc = particle->Pdg();
248  bool is_p = pdg::IsProton (pdgc);
249  bool is_n = pdg::IsNeutron(pdgc);
250 
251  if (is_p) Z--;
252  if (is_p || is_n) A--;
253 
254  Px += particle->Px();
255  Py += particle->Py();
256  Pz += particle->Pz();
257  E += particle->E();
258 
259  }//daughters
260 
261  TParticlePDG * remn = 0;
262  int ipdgc = pdg::IonPdgCode(A, Z);
263  remn = PDGLibrary::Instance()->Find(ipdgc);
264  if(!remn) {
265  LOG("HadronicVtx", pFATAL)
266  << "No particle with [A = " << A << ", Z = " << Z
267  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
268  assert(remn);
269  }
270 
271  double Mi = nucleus->Mass();
272  Px *= -1;
273  Py *= -1;
274  Pz *= -1;
275  E = Mi-E;
276 
277  // Add the nucleus to the event record
278  LOG("FermiMover", pINFO)
279  << "Adding nucleus [A = " << A << ", Z = " << Z
280  << ", pdgc = " << ipdgc << "]";
281 
282  int imom = evrec->TargetNucleusPosition();
283  evrec->AddParticle(
284  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
285 }
286 //___________________________________________________________________________
287 void FermiMover::Configure(const Registry & config)
288 {
289  Algorithm::Configure(config);
290  this->LoadConfig();
291 }
292 //____________________________________________________________________________
293 void FermiMover::Configure(string config)
294 {
295  Algorithm::Configure(config);
296  this->LoadConfig();
297 }
298 //____________________________________________________________________________
300 {
301  RgKey nuclkey = "NuclearModel";
302  fNuclModel = 0;
303  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
304  assert(fNuclModel);
305 
306  this->GetParamDef("KeepHitNuclOnMassShell", fKeepNuclOnMassShell, false);
307 
308  bool mom_dep_energy_removal_def = false;
309  this->GetParamDef("LFG-MomentumDependentErmv", mom_dep_energy_removal_def, false ) ;
310  // it defaults to whatever the nuclear model sets. Since only the LFG has this option
311  // this simple search is enough.
312 
313  this->GetParamDef("MomentumDependentErmv", fMomDepErmv, mom_dep_energy_removal_def);
314 
315  RgKey nuclearrecoilkey = "SecondNucleonEmitter" ;
316  fSecondEmitter = dynamic_cast<const SecondNucleonEmissionI *> (this->SubAlg(nuclearrecoilkey));
317 
318  assert(fSecondEmitter);
319 
320 }
321 //____________________________________________________________________________
int Z(void) const
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
static constexpr double rad
Definition: Units.h:164
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
void Configure(const Registry &config)
Definition: FermiMover.cxx:287
double RemovalEnergy(void) const
Definition: NuclearModelI.h:65
double E(void) const
Get energy.
Definition: GHepParticle.h:91
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
int FirstDaughter(void) const
Definition: GHepParticle.h:68
double Threshold(void) const
Energy threshold.
Definition: KPhaseSpace.cxx:82
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
void ProcessEventRecord(GHepRecord *event_rec) const
Definition: FermiMover.cxx:82
void KickHitNucleon(GHepRecord *evrec) const
give hit nucleon a momentum
Definition: FermiMover.cxx:100
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)
double Mass(void) const
Mass that corresponds to the PDG code.
bool fMomDepErmv
use momentum dependent calculation of Ermv
Definition: FermiMover.h:59
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
FermiMoverInteractionType_t GetFermiMoverInteractionType(void) const
Definition: NuclearModelI.h:80
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
int Pdg(void) const
Definition: GHepParticle.h:63
const SecondNucleonEmissionI * fSecondEmitter
Definition: FermiMover.h:62
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
Summary information for an interaction.
Definition: Interaction.h:56
int LastDaughter(void) const
Definition: GHepParticle.h:69
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#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
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
Kinematical phase space.
Definition: KPhaseSpace.h:33
const NuclearModelI * fNuclModel
nuclear model
Definition: FermiMover.h:60
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
void LoadConfig(void)
Definition: FermiMover.cxx:299
enum genie::EFermiMoverInteractionType FermiMoverInteractionType_t
void SetRemovalEnergy(double Erm)
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
Interface to drive the a second nucleon emission from a nucleus Specfic impelmentations will have dif...
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:117
string RgKey
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:313
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
const InitialState & InitState(void) const
Definition: Interaction.h:69
const int kPdgProton
Definition: PDGCodes.h:81
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
void AddTargetNucleusRemnant(GHepRecord *evrec) const
add a recoiled nucleus remnant
Definition: FermiMover.cxx:224
int A(void) const
#define pNOTICE
Definition: Messenger.h:61
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:66
bool fKeepNuclOnMassShell
keep hit bound nucleon on the mass shell?
Definition: FermiMover.h:58
virtual bool GenerateNucleon(const Target &) const =0
double ProbeE(RefFrame_t rf) const
const int kPdgNeutron
Definition: PDGCodes.h:83
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
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
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:370
Initial State information.
Definition: InitialState.h:48
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345