Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
INukeNucleonCorr Class Reference

Correction to free NN xsec in nuclear matter. More...

#include <INukeNucleonCorr.h>

Public Member Functions

double getAvgCorrection (const double rho, const double A, const double Ek)
 get the correction for given four-momentum and density More...
 
void OutputFiles (int A, int Z)
 
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 More...
 

Static Public Member Functions

static INukeNucleonCorrgetInstance ()
 get single instance of INukeNucleonCorr; create if necessary More...
 

Private Member Functions

 INukeNucleonCorr ()
 private constructor (called only by getInstance()) More...
 
 INukeNucleonCorr (const INukeNucleonCorr &)
 block copy constructor More...
 
INukeNucleonCorroperator= (const INukeNucleonCorr &)
 block assignment operator More...
 
double beta (const double rho)
 potential component (Eq. 2.18) More...
 
double lambda (const double rho)
 potential component (Eq. 2.19) More...
 
void setFermiLevel (const double rho, const int A, const int Z)
 
double fermiMomentum (const int pdg)
 return proper Fermi momentum based on nucleon PDG More...
 
double mstar (const double rho, const double k2)
 m* calculated based on Eqs. 2.6 and 2.16 More...
 
double localFermiMom (const double rho, const int A, const int Z, const int pdg)
 calculate local Fermi momentum More...
 
TLorentzVector generateTargetNucleon (const double mass, const double fermiMomentum)
 generate target nucleon More...
 
double getCorrection (const double mass, const double rho, const TVector3 &k1, const TVector3 &k2, const TVector3 &k3, const TVector3 &k4)
 calculate xsec correction More...
 

Private Attributes

double fFermiMomProton
 
double fFermiMomNeutron
 

Static Private Attributes

static INukeNucleonCorrfInstance = NULL
 single instance of INukeNucleonCorr More...
 
static const unsigned int fRepeat = 1000
 number of repetition to get average correction More...
 
static const double fRho0 = 0.16
 equilibrium density More...
 
static const double fAlpha1
 alpha coefficient as defined by Eq. 2.17 More...
 
static const double fAlpha2
 alpha coefficient as defined by Eq. 2.17 More...
 
static const double fBeta1 = -116.00 / fRho0 / 1000.0
 beta coefficient as defined by Eq. 2.18 More...
 
static const double fLambda0 = 3.29 / (units::fermi)
 lambda coefficient as defined by Eq. 2.19 More...
 
static const double fLambda1 = -0.373 / (units::fermi) / fRho0
 lambda coefficient as defined by Eq. 2.19 More...
 
static const int fNDensityBins
 cache binning for density More...
 
static const double fDensityStep
 within this density step correction is assumed to be constant More...
 
static const int fNEnergyBins1
 cache binning for energy More...
 
static const double fMaxEnergy1
 above this energy correction is assumed to be constant More...
 
static const double fEnergyStep1
 within this energy step correction is assumed to be constant More...
 
static const int fNEnergyBins2
 cache binning for energy More...
 
static const double fMaxEnergy2
 above this energy correction is assumed to be constant More...
 
static const double fEnergyStep2
 within this energy step correction is assumed to be constant More...
 
static const int fNEnergyBins3
 cache binning for energy More...
 
static const double fMaxEnergy3
 above this energy correction is assumed to be constant More...
 
static const double fEnergyStep3
 within this energy step correction is assumed to be constant More...
 

Detailed Description

Correction to free NN xsec in nuclear matter.

Author
Kyle Bachinski, Tomasz Golan
Date
2015
Remarks
V.R. Pandharipande and S. C. Pieper, Phys. Rev. C45 (1992) 791

Definition at line 18 of file INukeNucleonCorr.h.

Constructor & Destructor Documentation

INukeNucleonCorr::INukeNucleonCorr ( )
inlineprivate

private constructor (called only by getInstance())

Definition at line 78 of file INukeNucleonCorr.h.

INukeNucleonCorr::INukeNucleonCorr ( const INukeNucleonCorr )
private

block copy constructor

Member Function Documentation

double INukeNucleonCorr::AvgCorrection ( const double  rho,
const int  A,
const int  Z,
const int  pdg,
const double  Ek 
)

generate kinematics fRepeat times to calculate average correction

Definition at line 147 of file INukeNucleonCorr.cxx.

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 }
void setFermiLevel(const double rho, const int A, const int Z)
double getCorrection(const double mass, const double rho, const TVector3 &k1, const TVector3 &k2, const TVector3 &k3, const TVector3 &k4)
calculate xsec correction
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
double fermiMomentum(const int pdg)
return proper Fermi momentum based on nucleon PDG
TLorentzVector generateTargetNucleon(const double mass, const double fermiMomentum)
generate target nucleon
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils.cxx:754
p
Definition: test.py:223
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
#define A
Definition: memgrp.cpp:38
const int kPdgProton
Definition: PDGCodes.h:81
const int kPdgNeutron
Definition: PDGCodes.h:83
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
double INukeNucleonCorr::beta ( const double  rho)
inlineprivate

potential component (Eq. 2.18)

Definition at line 84 of file INukeNucleonCorr.h.

double INukeNucleonCorr::fermiMomentum ( const int  pdg)
inlineprivate

return proper Fermi momentum based on nucleon PDG

Definition at line 94 of file INukeNucleonCorr.h.

TLorentzVector INukeNucleonCorr::generateTargetNucleon ( const double  mass,
const double  fermiMomentum 
)
private

generate target nucleon

generate random momentum direction and return 4-momentum of target nucleon

Definition at line 117 of file INukeNucleonCorr.cxx.

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 }
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
p
Definition: test.py:223
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
static const double kPi
Definition: Constants.h:37
double INukeNucleonCorr::getAvgCorrection ( const double  rho,
const double  A,
const double  Ek 
)

get the correction for given four-momentum and density

Definition at line 237 of file INukeNucleonCorr.cxx.

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 }
void read_file(string rfilename)
vector< vector< double > > UraniumValues
string dir
vector< vector< double > > infile_values
const int NColumns
vector< vector< double > > TinValues
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
vector< vector< double > > CarbonValues
vector< vector< double > > HeliumValues
vector< vector< double > > IronValues
#define A
Definition: memgrp.cpp:38
const int NRows
vector< vector< double > > clear
vector< vector< double > > CalciumValues
#define pNOTICE
Definition: Messenger.h:61
double INukeNucleonCorr::getCorrection ( const double  mass,
const double  rho,
const TVector3 &  k1,
const TVector3 &  k2,
const TVector3 &  k3,
const TVector3 &  k4 
)
private

calculate xsec correction

calculate correction given by eq. 2.9

Definition at line 136 of file INukeNucleonCorr.cxx.

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 }
double mstar(const double rho, const double k2)
m* calculated based on Eqs. 2.6 and 2.16
static INukeNucleonCorr* INukeNucleonCorr::getInstance ( )
inlinestatic

get single instance of INukeNucleonCorr; create if necessary

Definition at line 23 of file INukeNucleonCorr.h.

23 {return fInstance ? fInstance : (fInstance = new INukeNucleonCorr);}
static INukeNucleonCorr * fInstance
single instance of INukeNucleonCorr
INukeNucleonCorr()
private constructor (called only by getInstance())
double INukeNucleonCorr::lambda ( const double  rho)
inlineprivate

potential component (Eq. 2.19)

Definition at line 85 of file INukeNucleonCorr.h.

double INukeNucleonCorr::localFermiMom ( const double  rho,
const int  A,
const int  Z,
const int  pdg 
)
private

calculate local Fermi momentum

$k_F = (\frac{3}{2}\pi^2\rho)^{1/3}$

Definition at line 101 of file INukeNucleonCorr.cxx.

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 }
constexpr T pow(T x)
Definition: pow.h:72
#define A
Definition: memgrp.cpp:38
static constexpr double fermi
Definition: Units.h:55
const int kPdgProton
Definition: PDGCodes.h:81
static const double kPi
Definition: Constants.h:37
double INukeNucleonCorr::mstar ( const double  rho,
const double  k2 
)
private

m* calculated based on Eqs. 2.6 and 2.16

$m^* (k,\rho) = m \frac{(\Lambda^2 + k^2)^2}{\Lambda^2 + k^2)^2 - 2\Lambda^2\beta m}$

Definition at line 83 of file INukeNucleonCorr.cxx.

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 }
double beta(const double rho)
potential component (Eq. 2.18)
double lambda(const double rho)
potential component (Eq. 2.19)
const int kPdgProton
Definition: PDGCodes.h:81
const int kPdgNeutron
Definition: PDGCodes.h:83
INukeNucleonCorr& INukeNucleonCorr::operator= ( const INukeNucleonCorr )
private

block assignment operator

void INukeNucleonCorr::OutputFiles ( int  A,
int  Z 
)

Definition at line 334 of file INukeNucleonCorr.cxx.

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
std::string string
Definition: nybbler.cc:12
const double e
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
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
string genie_dir(std::getenv("GENIE"))
#define A
Definition: memgrp.cpp:38
QTextStream & endl(QTextStream &s)
void INukeNucleonCorr::setFermiLevel ( const double  rho,
const int  A,
const int  Z 
)
inlineprivate
Parameters
Zset up Fermi momenta

Definition at line 87 of file INukeNucleonCorr.h.

88  {
89  fFermiMomProton = localFermiMom (rho, A, Z, genie::kPdgProton); // local Fermi momentum for protons
90  fFermiMomNeutron = localFermiMom (rho, A, Z, genie::kPdgNeutron); // local Fermi momentum for neutrons
91  }
const int kPdgProton
Definition: PDGCodes.h:81
double localFermiMom(const double rho, const int A, const int Z, const int pdg)
calculate local Fermi momentum
const int kPdgNeutron
Definition: PDGCodes.h:83

Member Data Documentation

const double INukeNucleonCorr::fAlpha1
staticprivate

alpha coefficient as defined by Eq. 2.17

Definition at line 43 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fAlpha2
staticprivate

alpha coefficient as defined by Eq. 2.17

Definition at line 44 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fBeta1 = -116.00 / fRho0 / 1000.0
staticprivate

beta coefficient as defined by Eq. 2.18

Definition at line 45 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fDensityStep
staticprivate

within this density step correction is assumed to be constant

Definition at line 53 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fEnergyStep1
staticprivate

within this energy step correction is assumed to be constant

Definition at line 61 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fEnergyStep2
staticprivate

within this energy step correction is assumed to be constant

Definition at line 64 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fEnergyStep3
staticprivate

within this energy step correction is assumed to be constant

Definition at line 67 of file INukeNucleonCorr.h.

double INukeNucleonCorr::fFermiMomNeutron
private

Definition at line 74 of file INukeNucleonCorr.h.

double INukeNucleonCorr::fFermiMomProton
private

Definition at line 73 of file INukeNucleonCorr.h.

INukeNucleonCorr * INukeNucleonCorr::fInstance = NULL
staticprivate

single instance of INukeNucleonCorr

Definition at line 33 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fLambda0 = 3.29 / (units::fermi)
staticprivate

lambda coefficient as defined by Eq. 2.19

Definition at line 46 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fLambda1 = -0.373 / (units::fermi) / fRho0
staticprivate

lambda coefficient as defined by Eq. 2.19

Definition at line 47 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fMaxEnergy1
staticprivate

above this energy correction is assumed to be constant

Definition at line 60 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fMaxEnergy2
staticprivate

above this energy correction is assumed to be constant

Definition at line 63 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fMaxEnergy3
staticprivate

above this energy correction is assumed to be constant

Definition at line 66 of file INukeNucleonCorr.h.

const int INukeNucleonCorr::fNDensityBins
staticprivate

cache binning for density

Definition at line 49 of file INukeNucleonCorr.h.

const int INukeNucleonCorr::fNEnergyBins1
staticprivate

cache binning for energy

Definition at line 59 of file INukeNucleonCorr.h.

const int INukeNucleonCorr::fNEnergyBins2
staticprivate

cache binning for energy

Definition at line 62 of file INukeNucleonCorr.h.

const int INukeNucleonCorr::fNEnergyBins3
staticprivate

cache binning for energy

Definition at line 65 of file INukeNucleonCorr.h.

const unsigned int INukeNucleonCorr::fRepeat = 1000
staticprivate

number of repetition to get average correction

Definition at line 37 of file INukeNucleonCorr.h.

const double INukeNucleonCorr::fRho0 = 0.16
staticprivate

equilibrium density

Definition at line 41 of file INukeNucleonCorr.h.


The documentation for this class was generated from the following files: