GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
INukeNucleonCorr.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  Author: Tomek Golan <tomasz.golan@uwr.edu.pl>, FNAL/Rochester
8  Steve Dytman <dytman+@pitt.edu>, Pittsburgh Univ.
9  Josh Kleckner <jok84@pitt.edu>, Pittsburgh Univ.
10 
11  Auguest 20, 2016
12 
13  Calculate the cross section attenuation factor due to medium
14 corrections for NA FSI from Pandharipande, Pieper (Phys Rev (2009).
15 Paper gives prescription for Pauli blocking and average nuclear
16 potential in (e,e'p). GENIE code was adapted from NuWro implementation.
17 
18  Important revisions after version 2.12.0 :
19  @ Aug, 2016 - TK
20  adapted to GENIE from NuWro. Use free NN xs.
21  @ Aug, 2017 - SD, JK
22  Original code stores values in rotating buffer. This won't be accurate
23  for many problems, esp. heavy targets and multiple nuclei. New code has
24  lookup tables in probe KE and nuclear density (rho) stored in text files
25  for He4, C12, Ca40, Fe56, Sn120, and U238. Use values from the text
26  files for KE and rho, interpolation in A.
27 */
28 //____________________________________________________________________________
37 using namespace genie;
38 using namespace genie::constants;
39 
40 #include <vector>
41 #include <string>
42 #include <fstream>
43 #include <sstream>
44 #include <iomanip>
45 #include <iostream>
46 #include <cmath>
47 #include <cstdlib>
48 
49 #include <TGraph.h>
50 using namespace std;
51 
52 
53 INukeNucleonCorr* INukeNucleonCorr::fInstance = NULL; // initialize instance with NULL
54 
55 
56 const int NRows = 200;
57 const int NColumns = 17;
58 
59 string genie_dir(std::getenv("GENIE"));
60 //string directory = genie_dir + string("/data/evgen/nncorr/");
61 string dir = genie_dir + string("/data/evgen/nncorr/");
62 string infile;
63 vector<vector<double> > infile_values;
64 vector<string> comments;
65 
66 
67 
68 // ----- STATIC CONSTANTS ----- //
69 
70 const unsigned int INukeNucleonCorr::fRepeat = 1000; // how many times kinematics should be generated to get avg corr
71 
72 const double INukeNucleonCorr::fRho0 = 0.16; // fm^-3
73 
74 const double INukeNucleonCorr::fBeta1 = -116.00 / fRho0 / 1000.0; // converted to GeV
75 const double INukeNucleonCorr::fLambda0 = 3.29 / (units::fermi); // converted to GeV
76 const double INukeNucleonCorr::fLambda1 = -0.373 / (units::fermi) / fRho0; // converted to GeV
77 
78 
79 
80 // ----- CALCULATIONS ----- //
81 
82 //! \f$m^* (k,\rho) = m \frac{(\Lambda^2 + k^2)^2}{\Lambda^2 + k^2)^2 - 2\Lambda^2\beta m}\f$
83 double INukeNucleonCorr::mstar (const double rho, const double k2)
84 {
85  // density [fm^-3], momentum square [GeV^2]
86 
87  static const double m = (PDGLibrary::Instance()->Find(kPdgProton)->Mass() +
88  PDGLibrary::Instance()->Find(kPdgNeutron)->Mass()) / 2.0;
89 
90  const double L = lambda (rho); // potential coefficient lambda
91  const double B = beta (rho); // potential coefficient beta
92 
93  const double L2 = L * L; // lambda^2 used twice in equation
94 
95  const double num = (L2 + k2) * (L2 + k2); // numerator
96 
97  return m * num / (num - 2.0 * L2 * B * m);
98 }
99 
100 //! \f$k_F = (\frac{3}{2}\pi^2\rho)^{1/3}\f$
101 double INukeNucleonCorr :: localFermiMom (const double rho, const int A, const int Z, const int pdg)
102 {
103  static double factor = 3.0 * kPi * kPi / 2.0;
104 
105  return pdg == kPdgProton ? pow (factor * rho * Z / A, 1.0 / 3.0) / (units::fermi) :
106  pow (factor * rho * (A - Z) / A, 1.0 / 3.0) / (units::fermi);
107 }
108  vector<vector<double> > HeliumValues;
109  vector<vector<double> > CarbonValues;
110  vector<vector<double> > CalciumValues;
111  vector<vector<double> > IronValues;
112  vector<vector<double> > TinValues;
113  vector<vector<double> > UraniumValues;
114  vector<vector<double> > clear;
115 
116 //! generate random momentum direction and return 4-momentum of target nucleon
117 TLorentzVector INukeNucleonCorr :: generateTargetNucleon (const double mass, const double fermiMom)
118 {
119  RandomGen * rnd = RandomGen::Instance();
120 
121  // get random momentum direction
122  const double costheta = 2.0 * rnd->RndGen().Rndm() - 1.0; // random cos (theta)
123  const double sintheta = sqrt (1.0 - costheta * costheta); // sin (theta)
124  const double phi = 2.0 * kPi * rnd->RndGen().Rndm(); // random phi
125 
126  // set nucleon 4-momentum
127  const double p = rnd->RndGen().Rndm() * fermiMom; // random nucleon momentum up to Fermi level
128 
129  const TVector3 p3 = TVector3 (p * sintheta * cos (phi), p * sintheta * sin (phi), p * costheta); // 3-momentum
130  const double energy = sqrt (p3.Mag2() + mass * mass); // energy
131 
132  return TLorentzVector (p3, energy);
133 }
134 
135 //! calculate correction given by eq. 2.9
136 double INukeNucleonCorr :: getCorrection (const double mass, const double rho,
137  const TVector3 &k1, const TVector3 &k2,
138  const TVector3 &k3, const TVector3 &k4)
139 {
140  const double num = (k1 - k2).Mag() * mstar (rho, (k3.Mag2() + k4.Mag2()) / 2.0) / mass / mass;
141  const double den = (k1 * (1.0 / mstar (rho, k1.Mag2())) - k2 * (1.0 / mstar (rho, k2.Mag2()))).Mag();
142 
143  return num / den;
144 }
145 
146 //! generate kinematics fRepeat times to calculate average correction
147 double INukeNucleonCorr :: AvgCorrection (const double rho, const int A, const int Z, const int pdg, const double Ek)
148 {
149 
150  RandomGen * rnd = RandomGen::Instance();
151 
152  setFermiLevel (rho, A, Z); // set Fermi momenta for protons and neutrons
153 
154  const double mass = PDGLibrary::Instance()->Find(pdg)->Mass(); // mass of incoming nucleon
155  const double energy = Ek + mass;
156 
157  TLorentzVector p (0.0, 0.0, sqrt (energy * energy - mass * mass), energy); // incoming particle 4-momentum
158  GHepParticle incomingParticle (pdg, kIStUndefined, -1,-1,-1,-1, p, TLorentzVector ()); // for IntBounce
159 
160  double corrPauliBlocking = 0.0; // correction coming from Pauli blocking
161  double corrPotential = 0.0; // correction coming from potential
162 
163  for (unsigned int i = 0; i < fRepeat; i++) // generate kinematics fRepeat times to get avg corrections
164  {
165  // get proton vs neutron randomly based on Z/A
166  const int targetPdg = rnd->RndGen().Rndm() < (double) Z / A ? kPdgProton : kPdgNeutron;
167 
168  const double targetMass = PDGLibrary::Instance()->Find(targetPdg)->Mass(); // set nucleon mass
169 
170  const TLorentzVector target = generateTargetNucleon (targetMass, fermiMomentum (targetPdg)); // generate target nucl
171 
172  TLorentzVector outNucl1, outNucl2, RemnP4; // final 4-momenta
173 
174  // random scattering angle
175  double C3CM = INukeHadroData2018::Instance()->IntBounce (&incomingParticle, targetPdg, pdg, kIHNFtElas);
176 
177  // generate kinematics
178  utils::intranuke2018::TwoBodyKinematics (mass, targetMass, p, target, outNucl1, outNucl2, C3CM, RemnP4);
179 
180  // update Pauli blocking correction
181  corrPauliBlocking += (outNucl1.Vect().Mag() > fermiMomentum (pdg) and outNucl2.Vect().Mag() > fermiMomentum (targetPdg));
182 
183  // update potential-based correction
184  corrPotential += getCorrection (mass, rho, p.Vect(), target.Vect(), outNucl1.Vect(), outNucl2.Vect());
185  }
186 
187  corrPauliBlocking /= fRepeat;
188  corrPotential /= fRepeat;
189 
190  return corrPotential * corrPauliBlocking;
191 
192 }
193 
194 
195 //This function reads the correction files that will be used to interpolate new correction values for some target//
196 void read_file(string rfilename)
197 {
198  ifstream file;
199  file.open((char*)rfilename.c_str(), ios::in);
200 
201  if (file.is_open())
202  {
203  string line;
204  int cur_line = 0;
205  while (getline(file,line))
206  {
207  if (line[0]=='#')
208  {
209  comments.push_back(line);
210  }
211  else {
212  vector<double> temp_vector;
213  istringstream iss(line);
214  string s;
215  for (int i=0; i<18; i++)
216  {
217  iss >> s;
218  temp_vector.push_back(atof(s.c_str()));
219  }
220  infile_values.push_back(temp_vector);
221  cur_line++;
222  }
223  }
224  LOG("INukeNucleonCorr",pNOTICE) << "Successful open file" << rfilename << "\n";
225 
226  }
227  else {
228  LOG("INukeNucleonCorr",pNOTICE) << "Could not open " << rfilename << "\n";
229 
230  }
231  file.close();
232 }
233 
234 
235 // This function interpolates and returns correction values
236 //
237 double INukeNucleonCorr :: getAvgCorrection(double rho, double A, double ke)
238 {
239  //Read in energy and density to determine the row and column of the correction table - adjust for variable binning - throws away some of the accuracy
240  int Column = round(rho*100);
241  if(rho<.01) Column = 1;
242  if (Column>=NColumns) Column = NColumns-1;
243  int Row = 0;
244  if(ke<=.002) Row = 1;
245  if(ke>.002&&ke<=.1) Row = round(ke*1000.);
246  if(ke>.1&&ke<=.5) Row = round(.1*1000.+(ke-.1)*200);
247  if(ke>.5&&ke<=1) Row = round(.1*1000.+(.5-.1)*200+(ke-.5)*40);
248  if(ke>1) Row = NRows-1;
249  //LOG ("INukeNucleonCorr",pNOTICE)
250  // << "row, column = " << Row << " " << Column;
251  //If the table of correction values has already been created
252  // return a value. Else, interpolate the needed correction table//
253 // static double cache[NRows][NColumns] = {{-1}};
254  static bool ReadFile = false;
255  if( ReadFile == true ) {
256  int Npoints = 6;
257  TGraph * Interp = new TGraph(Npoints);
258  //LOG("INukeNucleonCorr",pNOTICE)
259  // << HeliumValues[Row][Column];
260  Interp->SetPoint(0,4,HeliumValues[Row][Column]);
261  Interp->SetPoint(1,12,CarbonValues[Row][Column]);
262  Interp->SetPoint(2,40,CalciumValues[Row][Column]);
263  Interp->SetPoint(3,56,IronValues[Row][Column]);
264  Interp->SetPoint(4,120,TinValues[Row][Column]);
265  Interp->SetPoint(5,238,UraniumValues[Row][Column]);
266 
267  // Interpolated[e][r] = Interp->Eval(A);
268  // LOG("INukeNucleonCorr",pNOTICE)
269  // << "e,r,value= " << e << " " << r << " " << Interpolated[e][r];
270  double returnval = Interp->Eval(A);
271  delete Interp;
272  LOG("INukeNucleonCorr",pINFO)
273  << "Nucleon Corr interpolated correction factor = "
274  << returnval //cache[Row][Column]
275  << " for rho, KE, A= "<< rho << " " << ke << " " << A;
276  // return cache[Row][Column];
277  return returnval;
278  } else {
279  //Reading in correction files//
280  // string dir = genie_dir + string("/data/evgen/nncorr/");
281 
282  read_file(dir+"NNCorrection_2_4.txt");
285 
286  read_file(dir+"NNCorrection_6_12.txt");
289 
290  read_file(dir+"NNCorrection_20_40.txt");
293 
294  read_file(dir+"NNCorrection_26_56.txt");
297 
298  read_file(dir+"NNCorrection_50_120.txt");
301 
302  read_file(dir+"NNCorrection_92_238.txt");
305 
306  LOG("INukeNucleonCorr",pNOTICE)
307  << "Nucleon Corr interpolation files read in successfully";
308  //get interpolated value for first event.
309  int Npoints = 6;
310  TGraph * Interp = new TGraph(Npoints);
311  //LOG("INukeNucleonCorr",pNOTICE)
312  // << HeliumValues[Row][Column];
313  Interp->SetPoint(0,4,HeliumValues[Row][Column]);
314  Interp->SetPoint(1,12,CarbonValues[Row][Column]);
315  Interp->SetPoint(2,40,CalciumValues[Row][Column]);
316  Interp->SetPoint(3,56,IronValues[Row][Column]);
317  Interp->SetPoint(4,120,TinValues[Row][Column]);
318  Interp->SetPoint(5,238,UraniumValues[Row][Column]);
319 
320  // LOG("INukeNucleonCorr",pNOTICE)
321  // << "Row,Column,value= " << Row << " " << Column << " " << Interp->Eval(A);
322  double returnval = Interp->Eval(A);
323  delete Interp;
324  ReadFile = true;
325  LOG("INukeNucleonCorr",pINFO)
326  << "Nucleon Corr interpolated correction factor = "
327  << returnval
328  << " for rho, KE, A= "<< rho << " " << ke << " " << A;
329  return returnval;
330  }
331 }
332 
333 //This function outputs new correction files a new target if needed//
335 {
336  string outputdir = genie_dir + string("/data/evgen/nncorr/");
337  double pdgc;
338  string file;
339  string header;
340 
341 
342  file = outputdir+"NNCorrection.txt";
343  header = "##Correction values for protons (density(horizontal) and energy(vertical))";
344  pdgc = 2212;
345 
346 
347  double output[1002][18];
348  //label densities (columns)//
349  for(int r = 1; r < 18; r++)
350  { double rho = (r-1)*0.01;
351  output[0][r] = rho;}
352 
353  //label energies (rows) //
354  for(int e = 1; e < 1002; e++)
355  { double energy = (e-1)*0.001;
356  output[e][0] = energy;}
357 
358  //loop over each energy and density to get corrections and build the correction table//
359  for(int e = 1; e < 1002; e++){
360  for(int r = 1; r < 18; r++){
361  double energy = (e-1)*0.001;
362  double density = (r-1)*0.01;
363  double correction = INukeNucleonCorr::getInstance()-> AvgCorrection (density, A, Z, pdgc, energy);
364  output[e][r] = correction;
365  }
366  }
367  //output the new correction table //
368  ofstream outfile;
369  outfile.open((char*)file.c_str(), ios::trunc);
370  if (outfile.is_open()) {
371 
372  outfile << header << endl;
373  outfile << left << setw(12) << "##";
374  for (int e = 0; e < 1002; e++) {
375  for (int r = 0; r < 18; r++) {
376  if((e == 0) && (r == 0)) {}
377  else{outfile << left << setw(12) << output[e][r];}
378  }
379  outfile << endl;
380 
381  }
382  }
383  outfile.close();
384 }
static INukeNucleonCorr * getInstance()
get single instance of INukeNucleonCorr; create if necessary
double getCorrection(const double mass, const double rho, const TVector3 &k1, const TVector3 &k2, const TVector3 &k3, const TVector3 &k4)
calculate xsec correction
static INukeNucleonCorr * fInstance
single instance of INukeNucleonCorr
double mstar(const double rho, const double k2)
m* calculated based on Eqs. 2.6 and 2.16
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void read_file(string rfilename)
vector< vector< double > > UraniumValues
static constexpr double s
Definition: Units.h:95
string dir
vector< vector< double > > infile_values
const int NColumns
static const double fRho0
equilibrium density
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
vector< vector< double > > TinValues
const double e
TLorentzVector generateTargetNucleon(const double mass, const double fermiMomentum)
generate target nucleon
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string infile
double AvgCorrection(const double rho, const int A, const int Z, const int pdg, const double Ek)
generate kinematics fRepeat times to calculate average correction
static constexpr double A
Definition: Units.h:74
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
#define pINFO
Definition: Messenger.h:62
Correction to free NN xsec in nuclear matter.
vector< vector< double > > CarbonValues
static const double fLambda0
lambda coefficient as defined by Eq. 2.19
static const double fLambda1
lambda coefficient as defined by Eq. 2.19
void OutputFiles(int A, int Z)
vector< vector< double > > HeliumValues
double getAvgCorrection(const double rho, const double A, const double Ek)
get the correction for given four-momentum and density
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
static const double fBeta1
beta coefficient as defined by Eq. 2.18
static INukeHadroData2018 * Instance(void)
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
static const unsigned int fRepeat
number of repetition to get average correction
vector< vector< double > > IronValues
string genie_dir(std::getenv("GENIE"))
const int NRows
static constexpr double fermi
Definition: Units.h:55
vector< vector< double > > clear
const int kPdgProton
Definition: PDGCodes.h:81
vector< vector< double > > CalciumValues
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
double localFermiMom(const double rho, const int A, const int Z, const int pdg)
calculate local Fermi momentum
vector< string > comments
const int kPdgNeutron
Definition: PDGCodes.h:83
static constexpr double m
Definition: Units.h:71
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)