INukeNucleonCorr.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: 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 //____________________________________________________________________________
36 using namespace genie;
37 
38 #include <vector>
39 #include <string>
40 #include <fstream>
41 #include <sstream>
42 #include <iomanip>
43 #include <iostream>
44 #include <cmath>
45 #include <cstdlib>
46 
47 #include <TGraph.h>
48 using namespace std;
49 
50 
51 INukeNucleonCorr* INukeNucleonCorr::fInstance = NULL; // initialize instance with NULL
52 
53 
54 const int NRows = 200;
55 const int NColumns = 17;
56 
57 string genie_dir(std::getenv("GENIE"));
58 //string directory = genie_dir + string("/data/evgen/nncorr/");
59 string dir = genie_dir + string("/data/evgen/nncorr/");
60 string infile;
61 vector<vector<double> > infile_values;
62 vector<string> comments;
63 
64 
65 
66 // ----- STATIC CONSTANTS ----- //
67 
68 const unsigned int INukeNucleonCorr::fRepeat = 1000; // how many times kinematics should be generated to get avg corr
69 
70 const double INukeNucleonCorr::fRho0 = 0.16; // fm^-3
71 
72 const double INukeNucleonCorr::fBeta1 = -116.00 / fRho0 / 1000.0; // converted to GeV
73 const double INukeNucleonCorr::fLambda0 = 3.29 / (units::fermi); // converted to GeV
74 const double INukeNucleonCorr::fLambda1 = -0.373 / (units::fermi) / fRho0; // converted to GeV
75 
76 
77 
78 // ----- CALCULATIONS ----- //
79 
80 //! \f$m^* (k,\rho) = m \frac{(\Lambda^2 + k^2)^2}{\Lambda^2 + k^2)^2 - 2\Lambda^2\beta m}\f$
81 double INukeNucleonCorr::mstar (const double rho, const double k2)
82 {
83  // density [fm^-3], momentum square [GeV^2]
84 
85  static const double m = (PDGLibrary::Instance()->Find(kPdgProton)->Mass() +
86  PDGLibrary::Instance()->Find(kPdgNeutron)->Mass()) / 2.0;
87 
88  const double L = lambda (rho); // potential coefficient lambda
89  const double B = beta (rho); // potential coefficient beta
90 
91  const double L2 = L * L; // lambda^2 used twice in equation
92 
93  const double num = (L2 + k2) * (L2 + k2); // numerator
94 
95  return m * num / (num - 2.0 * L2 * B * m);
96 }
97 
98 //! \f$k_F = (\frac{3}{2}\pi^2\rho)^{1/3}\f$
99 double INukeNucleonCorr :: localFermiMom (const double rho, const int A, const int Z, const int pdg)
100 {
101  static double factor = 3.0 * M_PI * M_PI / 2.0;
102 
103  return pdg == kPdgProton ? pow (factor * rho * Z / A, 1.0 / 3.0) / (units::fermi) :
104  pow (factor * rho * (A - Z) / A, 1.0 / 3.0) / (units::fermi);
105 }
106  vector<vector<double> > HeliumValues;
107  vector<vector<double> > CarbonValues;
108  vector<vector<double> > CalciumValues;
109  vector<vector<double> > IronValues;
110  vector<vector<double> > TinValues;
111  vector<vector<double> > UraniumValues;
112  vector<vector<double> > clear;
113 
114 //! generate random momentum direction and return 4-momentum of target nucleon
115 TLorentzVector INukeNucleonCorr :: generateTargetNucleon (const double mass, const double fermiMom)
116 {
118 
119  // get random momentum direction
120  const double costheta = 2.0 * rnd->RndGen().Rndm() - 1.0; // random cos (theta)
121  const double sintheta = sqrt (1.0 - costheta * costheta); // sin (theta)
122  const double phi = 2.0 * M_PI * rnd->RndGen().Rndm(); // random phi
123 
124  // set nucleon 4-momentum
125  const double p = rnd->RndGen().Rndm() * fermiMom; // random nucleon momentum up to Fermi level
126 
127  const TVector3 p3 = TVector3 (p * sintheta * cos (phi), p * sintheta * sin (phi), p * costheta); // 3-momentum
128  const double energy = sqrt (p3.Mag2() + mass * mass); // energy
129 
130  return TLorentzVector (p3, energy);
131 }
132 
133 //! calculate correction given by eq. 2.9
134 double INukeNucleonCorr :: getCorrection (const double mass, const double rho,
135  const TVector3 &k1, const TVector3 &k2,
136  const TVector3 &k3, const TVector3 &k4)
137 {
138  const double num = (k1 - k2).Mag() * mstar (rho, (k3.Mag2() + k4.Mag2()) / 2.0) / mass / mass;
139  const double den = (k1 * (1.0 / mstar (rho, k1.Mag2())) - k2 * (1.0 / mstar (rho, k2.Mag2()))).Mag();
140 
141  return num / den;
142 }
143 
144 //! generate kinematics fRepeat times to calculate average correction
145 double INukeNucleonCorr :: AvgCorrection (const double rho, const int A, const int Z, const int pdg, const double Ek)
146 {
147 
149 
150  setFermiLevel (rho, A, Z); // set Fermi momenta for protons and neutrons
151 
152  const double mass = PDGLibrary::Instance()->Find(pdg)->Mass(); // mass of incoming nucleon
153  const double energy = Ek + mass;
154 
155  TLorentzVector p (0.0, 0.0, sqrt (energy * energy - mass * mass), energy); // incoming particle 4-momentum
156  GHepParticle incomingParticle (pdg, kIStUndefined, -1,-1,-1,-1, p, TLorentzVector ()); // for IntBounce
157 
158  double corrPauliBlocking = 0.0; // correction coming from Pauli blocking
159  double corrPotential = 0.0; // correction coming from potential
160 
161  for (unsigned int i = 0; i < fRepeat; i++) // generate kinematics fRepeat times to get avg corrections
162  {
163  // get proton vs neutron randomly based on Z/A
164  const int targetPdg = rnd->RndGen().Rndm() < (double) Z / A ? kPdgProton : kPdgNeutron;
165 
166  const double targetMass = PDGLibrary::Instance()->Find(targetPdg)->Mass(); // set nucleon mass
167 
168  const TLorentzVector target = generateTargetNucleon (targetMass, fermiMomentum (targetPdg)); // generate target nucl
169 
170  TLorentzVector outNucl1, outNucl2, RemnP4; // final 4-momenta
171 
172  // random scattering angle
173  double C3CM = INukeHadroData2018::Instance()->IntBounce (&incomingParticle, targetPdg, pdg, kIHNFtElas);
174 
175  // generate kinematics
176  utils::intranuke2018::TwoBodyKinematics (mass, targetMass, p, target, outNucl1, outNucl2, C3CM, RemnP4);
177 
178  // update Pauli blocking correction
179  corrPauliBlocking += (outNucl1.Vect().Mag() > fermiMomentum (pdg) and outNucl2.Vect().Mag() > fermiMomentum (targetPdg));
180 
181  // update potential-based correction
182  corrPotential += getCorrection (mass, rho, p.Vect(), target.Vect(), outNucl1.Vect(), outNucl2.Vect());
183  }
184 
185  corrPauliBlocking /= fRepeat;
186  corrPotential /= fRepeat;
187 
188  return corrPotential * corrPauliBlocking;
189 
190 }
191 
192 
193 //This function reads the correction files that will be used to interpolate new correction values for some target//
194 void read_file(string rfilename)
195 {
196  ifstream file;
197  file.open((char*)rfilename.c_str(), ios::in);
198 
199  if (file.is_open())
200  {
201  string line;
202  int cur_line = 0;
203  while (getline(file,line))
204  {
205  if (line[0]=='#')
206  {
207  comments.push_back(line);
208  }
209  else {
210  vector<double> temp_vector;
211  istringstream iss(line);
212  string s;
213  for (int i=0; i<18; i++)
214  {
215  iss >> s;
216  temp_vector.push_back(atof(s.c_str()));
217  }
218  infile_values.push_back(temp_vector);
219  cur_line++;
220  }
221  }
222  LOG("INukeNucleonCorr",pNOTICE) << "Successful open file" << rfilename << "\n";
223 
224  }
225  else {
226  LOG("INukeNucleonCorr",pNOTICE) << "Could not open " << rfilename << "\n";
227 
228  }
229  file.close();
230 }
231 
232 
233 // This function interpolates and returns correction values
234 //
235 double INukeNucleonCorr :: getAvgCorrection(double rho, double A, double ke)
236 {
237  //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
238  int Column = round(rho*100);
239  if(rho<.01) Column = 1;
240  if (Column>=NColumns) Column = NColumns-1;
241  int Row = 0;
242  if(ke<=.002) Row = 1;
243  if(ke>.002&&ke<=.1) Row = round(ke*1000.);
244  if(ke>.1&&ke<=.5) Row = round(.1*1000.+(ke-.1)*200);
245  if(ke>.5&&ke<=1) Row = round(.1*1000.+(.5-.1)*200+(ke-.5)*40);
246  if(ke>1) Row = NRows-1;
247  //LOG ("INukeNucleonCorr",pNOTICE)
248  // << "row, column = " << Row << " " << Column;
249  //If the table of correction values has already been created
250  // return a value. Else, interpolate the needed correction table//
251 // static double cache[NRows][NColumns] = {{-1}};
252  static bool ReadFile = false;
253  if( ReadFile == true ) {
254  int Npoints = 6;
255  TGraph * Interp = new TGraph(Npoints);
256  //LOG("INukeNucleonCorr",pNOTICE)
257  // << HeliumValues[Row][Column];
258  Interp->SetPoint(0,4,HeliumValues[Row][Column]);
259  Interp->SetPoint(1,12,CarbonValues[Row][Column]);
260  Interp->SetPoint(2,40,CalciumValues[Row][Column]);
261  Interp->SetPoint(3,56,IronValues[Row][Column]);
262  Interp->SetPoint(4,120,TinValues[Row][Column]);
263  Interp->SetPoint(5,238,UraniumValues[Row][Column]);
264 
265  // Interpolated[e][r] = Interp->Eval(A);
266  // LOG("INukeNucleonCorr",pNOTICE)
267  // << "e,r,value= " << e << " " << r << " " << Interpolated[e][r];
268  double returnval = Interp->Eval(A);
269  delete Interp;
270  LOG("INukeNucleonCorr",pINFO)
271  << "Nucleon Corr interpolated correction factor = "
272  << returnval //cache[Row][Column]
273  << " for rho, KE, A= "<< rho << " " << ke << " " << A;
274  // return cache[Row][Column];
275  return returnval;
276  } else {
277  //Reading in correction files//
278  // string dir = genie_dir + string("/data/evgen/nncorr/");
279 
280  read_file(dir+"NNCorrection_2_4.txt");
281  HeliumValues = infile_values;
282  infile_values = clear;
283 
284  read_file(dir+"NNCorrection_6_12.txt");
285  CarbonValues = infile_values;
286  infile_values = clear;
287 
288  read_file(dir+"NNCorrection_20_40.txt");
289  CalciumValues = infile_values;
290  infile_values = clear;
291 
292  read_file(dir+"NNCorrection_26_56.txt");
293  IronValues = infile_values;
294  infile_values = clear;
295 
296  read_file(dir+"NNCorrection_50_120.txt");
297  TinValues = infile_values;
298  infile_values = clear;
299 
300  read_file(dir+"NNCorrection_92_238.txt");
301  UraniumValues = infile_values;
302  infile_values = clear;
303 
304  LOG("INukeNucleonCorr",pNOTICE)
305  << "Nucleon Corr interpolation files read in successfully";
306  //get interpolated value for first event.
307  int Npoints = 6;
308  TGraph * Interp = new TGraph(Npoints);
309  //LOG("INukeNucleonCorr",pNOTICE)
310  // << HeliumValues[Row][Column];
311  Interp->SetPoint(0,4,HeliumValues[Row][Column]);
312  Interp->SetPoint(1,12,CarbonValues[Row][Column]);
313  Interp->SetPoint(2,40,CalciumValues[Row][Column]);
314  Interp->SetPoint(3,56,IronValues[Row][Column]);
315  Interp->SetPoint(4,120,TinValues[Row][Column]);
316  Interp->SetPoint(5,238,UraniumValues[Row][Column]);
317 
318  // LOG("INukeNucleonCorr",pNOTICE)
319  // << "Row,Column,value= " << Row << " " << Column << " " << Interp->Eval(A);
320  double returnval = Interp->Eval(A);
321  delete Interp;
322  ReadFile = true;
323  LOG("INukeNucleonCorr",pINFO)
324  << "Nucleon Corr interpolated correction factor = "
325  << returnval
326  << " for rho, KE, A= "<< rho << " " << ke << " " << A;
327  return returnval;
328  }
329 }
330 
331 //This function outputs new correction files a new target if needed//
333 {
334  string outputdir = genie_dir + string("/data/evgen/nncorr/");
335  double pdgc;
336  string file;
337  string header;
338 
339 
340  file = outputdir+"NNCorrection.txt";
341  header = "##Correction values for protons (density(horizontal) and energy(vertical))";
342  pdgc = 2212;
343 
344 
345  double output[1002][18];
346  //label densities (columns)//
347  for(int r = 1; r < 18; r++)
348  { double rho = (r-1)*0.01;
349  output[0][r] = rho;}
350 
351  //label energies (rows) //
352  for(int e = 1; e < 1002; e++)
353  { double energy = (e-1)*0.001;
354  output[e][0] = energy;}
355 
356  //loop over each energy and density to get corrections and build the correction table//
357  for(int e = 1; e < 1002; e++){
358  for(int r = 1; r < 18; r++){
359  double energy = (e-1)*0.001;
360  double density = (r-1)*0.01;
361  double correction = INukeNucleonCorr::getInstance()-> AvgCorrection (density, A, Z, pdgc, energy);
362  output[e][r] = correction;
363  }
364  }
365  //output the new correction table //
366  ofstream outfile;
367  outfile.open((char*)file.c_str(), ios::trunc);
368  if (outfile.is_open()) {
369 
370  outfile << header << endl;
371  outfile << left << setw(12) << "##";
372  for (int e = 0; e < 1002; e++) {
373  for (int r = 0; r < 18; r++) {
374  if((e == 0) && (r == 0)) {}
375  else{outfile << left << setw(12) << output[e][r];}
376  }
377  outfile << endl;
378 
379  }
380  }
381  outfile.close();
382 }
ofstream output
const XML_Char * target
Definition: expat.h:268
static INukeNucleonCorr * getInstance()
get single instance of INukeNucleonCorr; create if necessary
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
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:79
void read_file(string rfilename)
static const double fermi
Definition: Units.h:63
Float_t den
Definition: plot.C:36
vector< vector< double > > UraniumValues
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
vector< vector< double > > clear
Double_t beta
string dir
vector< vector< double > > infile_values
const int NColumns
static const double fRho0
equilibrium density
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
fvar< T > trunc(const fvar< T > &x)
Definition: trunc.hpp:19
#define M_PI
Definition: SbMath.h:34
Float_t Z
Definition: plot.C:38
const XML_Char * s
Definition: expat.h:262
vector< vector< double > > TinValues
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:97
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 L
std::string getenv(std::string const &name)
double energy
Definition: plottest35.C:25
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:63
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)
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
vector< vector< double > > HeliumValues
static const double A
Definition: Units.h:82
double getAvgCorrection(const double rho, const double A, const double Ek)
get the correction for given four-momentum and density
ifstream in
Definition: comparison.C:7
T sin(T number)
Definition: d0nt_math.hpp:132
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static const double fBeta1
beta coefficient as defined by Eq. 2.18
static INukeHadroData2018 * Instance(void)
int num
Definition: f2_nu.C:119
cosmicTree ReadFile("CosmicList.txt","isSA/I:isA17/I:HadEfrac/F:CCE/F:CCESA/F")
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:81
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
static const unsigned int fRepeat
number of repetition to get average correction
float Mag() const
vector< vector< double > > IronValues
string genie_dir(std::getenv("GENIE"))
TFile * file
Definition: cellShifts.C:17
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
const int NRows
const int kPdgProton
Definition: PDGCodes.h:65
vector< vector< double > > CalciumValues
tuple density
Definition: barite.py:33
#define pNOTICE
Definition: Messenger.h:62
double localFermiMom(const double rho, const int A, const int Z, const int pdg)
calculate local Fermi momentum
Float_t e
Definition: plot.C:35
vector< string > comments
const int kPdgNeutron
Definition: PDGCodes.h:67
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
FILE * outfile
Definition: dump_event.C:13