NucBindEnergyAggregator.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 - November 19, 2004
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14 */
15 //____________________________________________________________________________
16 
17 #include <TLorentzVector.h>
18 
31 
32 using namespace genie;
33 using namespace genie::constants;
34 
35 //___________________________________________________________________________
37 EventRecordVisitorI("genie::NucBindEnergyAggregator")
38 {
39 
40 }
41 //___________________________________________________________________________
43 EventRecordVisitorI("genie::NucBindEnergyAggregator", config)
44 {
45 
46 }
47 //___________________________________________________________________________
49 {
50 
51 }
52 //___________________________________________________________________________
54 {
55  // Return if the neutrino was not scatterred off a nuclear target
56  GHepParticle * nucltgt = evrec->TargetNucleus();
57  if (!nucltgt) {
58  LOG("Nuclear", pINFO)
59  << "No nuclear target found - Not subtracting any binding energy";
60  return;
61  }
62 
63  // Loop over particles, find final state nucleons for which a removal
64  // energy has been set and subtract it from their energy
65  TIter stdhep_iter(evrec);
66  GHepParticle * p = 0;
67 
68  while( (p = (GHepParticle * ) stdhep_iter.Next()) ) {
69 
70  bool is_nucleon = pdg::IsNeutronOrProton(p->Pdg());
71  bool in_fin_state = (p->Status() == kIStStableFinalState);
72  bool had_bind_e = (p->RemovalEnergy() > 0.);
73 
74  bool handle = is_nucleon && in_fin_state && had_bind_e;
75  if(!handle) continue;
76 
77  //-- ask for the binding energy set by the nuclear model
78  double bindE = p->RemovalEnergy();
79  LOG("Nuclear", pINFO) << "Binding energy = " << bindE;
80 
81  //-- subtract this energy from the final state nucleon
82  LOG("Nuclear", pINFO)
83  << "Subtracting the binding energy from the escaped nucleon";
84 
85  double M = p->Mass();
86  double En = p->Energy();
87  double KE = En-M;
88 
89  LOG("Nuclear", pINFO) << "Kinetic energy before subtraction = " << KE;
90  KE -= bindE;
91  KE = TMath::Max(0.,KE);
92 
93  LOG("Nuclear", pINFO) << "Kinetic energy after subtraction = " << KE;
94 
95  En = KE+M;
96 
97  if(En>M || !fAllowRecombination) {
98  double pmag_old = p->P4()->P();
99  double pmag_new = TMath::Sqrt(utils::math::NonNegative(En*En-M*M));
100  double scale = pmag_new / pmag_old;
101  LOG("Nuclear", pINFO)
102  << "|pnew| = " << pmag_new << ", |pold| = " << pmag_old
103  << ", scale = " << scale;
104 
105  double pxn = scale * p->Px();
106  double pyn = scale * p->Py();
107  double pzn = scale * p->Pz();
108 
109  double pxb = (1-scale) * p->Px();
110  double pyb = (1-scale) * p->Py();
111  double pzb = (1-scale) * p->Pz();
112 
113  p->SetEnergy ( En );
114  p->SetPx ( pxn );
115  p->SetPy ( pyn );
116  p->SetPz ( pzn );
117 
118  //-- and add a GHEP entry to record this in the event record and
119  // conserve energy/momentum
120  LOG("Nuclear", pINFO)
121  << "Adding a [BindingE] to account for nuclear binding energy";
122 
124  -1,-1,-1,-1, pxb,pyb,pzb,bindE, 0,0,0,0);
125  } else {
126  LOG("Nuclear", pNOTICE)
127  << "Nucleon is above the Fermi sea but can't escape the nucleus";
128  LOG("Nuclear", pNOTICE)
129  << "Recombining remnant nucleus + f/s nucleon";
130 
131  LOG("Nuclear", pERROR)
132  << "*** This functionality is temporarily disabled";
133 /*
134  LOG("Nuclear", pNOTICE) << *evrec;
135 
136  // find the remnant nucleus
137  int rnucpos = evrec->RemnantNucleusPosition();
138  assert(rnucpos);
139 
140  GHepParticle * rnucl = evrec->Particle(rnucpos);
141 
142  // mark both the remnant nucleus and the final state nucleon as
143  // intermediate states
144  rnucl -> SetStatus(kIStIntermediateState);
145  p -> SetStatus(kIStIntermediateState);
146 
147  // figure out the recombined nucleus PDG code
148  int Z = rnucl->Z();
149  int A = rnucl->A();
150  if(pdg::IsProton(p->Pdg())) Z++;
151  A++;
152  int ipdgc = pdg::IonPdgCode(A,Z);
153 
154  // add-up their 4-momenta
155  double pxnuc = rnucl->Px() + p->Px();
156  double pynuc = rnucl->Py() + p->Py();
157  double pznuc = rnucl->Pz() + p->Pz();
158  double Enuc = rnucl->E() + p->E();
159 
160  evrec->AddParticle(ipdgc, kIStStableFinalState,
161  rnucpos,-1,-1,-1, pxnuc,pynuc,pznuc,Enuc, 0,0,0,0);
162 */
163  }
164  }
165 }
166 //___________________________________________________________________________
167 /*
168 GHepParticle * NucBindEnergyAggregator::FindMotherNucleus(
169  int ipos, GHepRecord * evrec) const
170 {
171  GHepParticle * p = evrec->Particle(ipos);
172 
173  //-- get its mothet
174  int mother_pos = p->FirstMother();
175 
176  //-- if mother is set
177  if(mother_pos != -1) {
178  GHepParticle * mother = evrec->Particle(mother_pos);
179 
180  //-- check its status
181  if( mother->Status() == kIStNucleonTarget ) {
182 
183  //-- get the mother's mother
184  int grandmother_pos = mother->FirstMother();
185 
186  //-- if grandmother is set get its PDG code a check if it is an ion
187  if(grandmother_pos != -1) {
188  GHepParticle * grandmother =
189  evrec->Particle(grandmother_pos);
190 
191  int grandmother_pdgc = grandmother->Pdg();
192  if( pdg::IsIon(grandmother_pdgc) ) return grandmother;
193 
194  } // gmpos != -1
195  } // mother-status
196  } //mpos != -1
197 
198  return 0;
199 } */
200 //___________________________________________________________________________
202 {
203  Algorithm::Configure(config);
204  this->LoadConfig();
205 }
206 //____________________________________________________________________________
208 {
209  Algorithm::Configure(config);
210  this->LoadConfig();
211 }
212 //____________________________________________________________________________
214 {
215  GetParamDef("AllowNuclRecombination",fAllowRecombination, true ) ;
216 
217 }
218 //____________________________________________________________________________
219 
Basic constants.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
void SetPz(double pz)
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
const char * p
Definition: xmltok.h:285
const int kPdgBindino
Definition: PDGCodes.h:189
double RemovalEnergy(void) const
Get removal energy.
Definition: GHepParticle.h:101
Definition: config.py:1
double Mass(void) const
Mass that corresponds to the PDG code.
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
void SetPx(double px)
double Energy(void) const
Get energy.
Definition: GHepParticle.h:93
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
int Pdg(void) const
Definition: GHepParticle.h:64
Double_t scale
Definition: plot.C:25
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
void Configure(const Registry &config)
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
#define pINFO
Definition: Messenger.h:63
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:314
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
double NonNegative(double x)
Definition: MathUtils.cxx:280
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
void SetEnergy(double E)
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
void SetPy(double py)
void ProcessEventRecord(GHepRecord *event_rec) const
double Py(void) const
Get Py.
Definition: GHepParticle.h:90