36 using namespace genie;
88 const double L = lambda (rho);
89 const double B =
beta (rho);
91 const double L2 = L *
L;
93 const double num = (L2 + k2) * (L2 + k2);
95 return m * num / (num - 2.0 * L2 * B *
m);
101 static double factor = 3.0 *
M_PI *
M_PI / 2.0;
120 const double costheta = 2.0 * rnd->
RndGen().Rndm() - 1.0;
121 const double sintheta =
sqrt (1.0 - costheta * costheta);
122 const double phi = 2.0 *
M_PI * rnd->
RndGen().Rndm();
125 const double p = rnd->
RndGen().Rndm() * fermiMom;
127 const TVector3
p3 = TVector3 (p * sintheta *
cos (phi), p * sintheta *
sin (phi), p * costheta);
128 const double energy =
sqrt (p3.Mag2() + mass * mass);
130 return TLorentzVector (p3, energy);
135 const TVector3 &k1,
const TVector3 &k2,
136 const TVector3 &k3,
const TVector3 &k4)
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();
150 setFermiLevel (rho, A, Z);
153 const double energy = Ek + mass;
155 TLorentzVector
p (0.0, 0.0,
sqrt (energy * energy - mass * mass), energy);
158 double corrPauliBlocking = 0.0;
159 double corrPotential = 0.0;
161 for (
unsigned int i = 0;
i < fRepeat;
i++)
168 const TLorentzVector
target = generateTargetNucleon (targetMass, fermiMomentum (targetPdg));
170 TLorentzVector outNucl1, outNucl2, RemnP4;
179 corrPauliBlocking += (outNucl1.Vect().Mag() > fermiMomentum (pdg) and outNucl2.Vect().Mag() > fermiMomentum (targetPdg));
182 corrPotential += getCorrection (mass, rho, p.Vect(), target.Vect(), outNucl1.Vect(), outNucl2.Vect());
185 corrPauliBlocking /= fRepeat;
186 corrPotential /= fRepeat;
188 return corrPotential * corrPauliBlocking;
197 file.open((
char*)rfilename.c_str(),
ios::in);
203 while (getline(file,line))
207 comments.push_back(line);
210 vector<double> temp_vector;
211 istringstream iss(line);
213 for (
int i=0;
i<18;
i++)
216 temp_vector.push_back(atof(s.c_str()));
218 infile_values.push_back(temp_vector);
222 LOG(
"INukeNucleonCorr",
pNOTICE) <<
"Successful open file" << rfilename <<
"\n";
226 LOG(
"INukeNucleonCorr",
pNOTICE) <<
"Could not open " << rfilename <<
"\n";
238 int Column = round(rho*100);
239 if(rho<.01) Column = 1;
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;
253 if( ReadFile ==
true ) {
255 TGraph * Interp =
new TGraph(Npoints);
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]);
268 double returnval = Interp->Eval(A);
271 <<
"Nucleon Corr interpolated correction factor = " 273 <<
" for rho, KE, A= "<< rho <<
" " << ke <<
" " <<
A;
282 infile_values =
clear;
286 infile_values =
clear;
290 infile_values =
clear;
294 infile_values =
clear;
296 read_file(dir+
"NNCorrection_50_120.txt");
298 infile_values =
clear;
300 read_file(dir+
"NNCorrection_92_238.txt");
302 infile_values =
clear;
305 <<
"Nucleon Corr interpolation files read in successfully";
308 TGraph * Interp =
new TGraph(Npoints);
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]);
320 double returnval = Interp->Eval(A);
324 <<
"Nucleon Corr interpolated correction factor = " 326 <<
" for rho, KE, A= "<< rho <<
" " << ke <<
" " <<
A;
340 file = outputdir+
"NNCorrection.txt";
341 header =
"##Correction values for protons (density(horizontal) and energy(vertical))";
347 for(
int r = 1;
r < 18;
r++)
348 {
double rho = (
r-1)*0.01;
352 for(
int e = 1;
e < 1002;
e++)
357 for(
int e = 1;
e < 1002;
e++){
358 for(
int r = 1;
r < 18;
r++){
362 output[
e][
r] = correction;
367 outfile.open((
char*)file.c_str(), ios::trunc);
368 if (outfile.is_open()) {
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];}
static INukeNucleonCorr * getInstance()
get single instance of INukeNucleonCorr; create if necessary
THE MAIN GENIE PROJECT NAMESPACE
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.
void read_file(string rfilename)
static const double fermi
vector< vector< double > > UraniumValues
vector< vector< double > > clear
vector< vector< double > > infile_values
static const double fRho0
equilibrium density
A singleton holding random number generator classes. All random number generation in GENIE should tak...
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...
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)
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
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)
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)
static const double fBeta1
beta coefficient as defined by Eq. 2.18
static INukeHadroData2018 * Instance(void)
cosmicTree ReadFile("CosmicList.txt","isSA/I:isA17/I:HadEfrac/F:CCE/F:CCESA/F")
TRandom3 & RndGen(void) const
rnd number generator for generic usage
TParticlePDG * Find(int pdgc)
static const unsigned int fRepeat
number of repetition to get average correction
vector< vector< double > > IronValues
string genie_dir(std::getenv("GENIE"))
vector< vector< double > > CalciumValues
double localFermiMom(const double rho, const int A, const int Z, const int pdg)
calculate local Fermi momentum
vector< string > comments
STDHEP-like event record entry that can fit a particle or a nucleus.
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)