GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SRCNuclearRecoil.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  or see $GENIE/LICENSE
6 
7  Author: Afroditi Papadopoulou <apapadop \at mit.edu>
8  Massachusetts Institute of Technology - October 04, 2019
9 
10  @ October 4, 2019 - Afroditi Papadopoulou (AP)
11  Created this new module that controls the addition of the recoil nucleon in the event record
12  and extracts its kinematics
13 */
14 //____________________________________________________________________________
15 
16 #include <cstdlib>
17 
18 #include <TLorentzVector.h>
19 #include <TVector3.h>
20 #include <TParticlePDG.h>
21 #include <TMath.h>
22 
28 
46 
47 using namespace genie;
48 using namespace genie::constants;
49 
50 //___________________________________________________________________________
52 SecondNucleonEmissionI("genie::SRCNuclearRecoil")
53 {
54 
55 }
56 //___________________________________________________________________________
58  SecondNucleonEmissionI("genie::SRCNuclearRecoil", config )
59 {
60 
61 }
62 
63 //___________________________________________________________________________
64 
66 {
67 
68 }
69 
70 //___________________________________________________________________________
71 
73 {
74 
75  const Interaction * interaction = evrec -> Summary();
76  const InitialState & init_state = interaction -> InitState();
77  const Target & tgt = init_state.Tgt();
78 
79  // do nothing for non-nuclear targets
80  if(! tgt.IsNucleus()) return;
81 
82  // access the hit nucleon and target nucleus at the GHEP record
83  GHepParticle * nucleon = evrec->HitNucleon();
84  GHepParticle * nucleus = evrec->TargetNucleus();
85  assert(nucleon);
86  assert(nucleus);
87 
88  // Set this to either a proton or neutron to eject a secondary particle
89  int eject_nucleon_pdg = this->SRCRecoilPDG( *nucleon, tgt );
90 
91  // Ejection of secondary particle
92  if (eject_nucleon_pdg != 0) { EmitSecondNucleon(evrec,eject_nucleon_pdg); }
93 
94 }
95 
96 //___________________________________________________________________________
97 
98 int SRCNuclearRecoil::SRCRecoilPDG( const GHepParticle & nucleon, const Target & tgt) const {
99 
100  int eject_nucleon_pdg = 0;
101 
102  // const int nucleus_pdgc = tgt->Pdg();
103  const int nucleon_pdgc = nucleon.Pdg();
104 
105  double pN2 = TMath::Power(nucleon.P4()->Rho(),2.); // (momentum of struck nucleon)^2
106 
107  double kF = fNuclModel -> LocalFermiMomentum( tgt,
108  nucleon_pdgc,
109  nucleon.X4()->Vect().Mag() );
110 
111  if (TMath::Sqrt(pN2) > kF) {
112  double Pp = (nucleon_pdgc == kPdgProton) ? fPPPairPercentage : fPNPairPercentage;
113  RandomGen * rnd = RandomGen::Instance();
114  double prob = rnd->RndGen().Rndm();
115  eject_nucleon_pdg = (prob > Pp) ? kPdgNeutron : kPdgProton;
116  }
117 
118  return eject_nucleon_pdg;
119 }
120 //___________________________________________________________________________
122 {
123  Algorithm::Configure(config);
124  this->LoadConfig();
125 }
126 //____________________________________________________________________________
127 void SRCNuclearRecoil::Configure(string config)
128 {
129  Algorithm::Configure(config);
130  this->LoadConfig();
131 }
132 //____________________________________________________________________________
134 {
135 
137 
138  this->GetParamDef("PNPairPercentage", fPNPairPercentage, 0.95);
139 
140  if (fPNPairPercentage < 0. || fPNPairPercentage > 1.) {
141 
142  LOG("SRCNuclearRecoil", pFATAL)
143  << "PNPairPercentage either less than 0 or greater than 1: Exiting" ;
144 
145  exit(78);
146  }
147 
149 
150 }
151 //____________________________________________________________________________
virtual bool EmitSecondNucleon(GHepRecord *evrec, const int eject_nucleon_pdg) const
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
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
void ProcessEventRecord(GHepRecord *event_rec) const
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
int SRCRecoilPDG(const GHepParticle &nucleon, const Target &tgt) const
Interface to drive the a second nucleon emission from a nucleus Specfic impelmentations will have dif...
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
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
const NuclearModelI * fNuclModel
nuclear model
const int kPdgProton
Definition: PDGCodes.h:81
void Configure(const Registry &config)
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:66
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
Initial State information.
Definition: InitialState.h:48