Skip to content

Commit

Permalink
Turn RPA effects off in the Nieves CCQE model if the kIAssumeFreeNucleon
Browse files Browse the repository at this point in the history 
bit is set in the input interaction.
  • Loading branch information
@sjgardiner
sjgardiner committed Apr 5, 2019
1 parent faa11d9 commit e55a227
Show file tree
Hide file tree
Showing 5 changed files with 43 additions and 69 deletions.
1 change: 1 addition & 0 deletions src/Physics/QuasiElastic/EventGen/QELEventGenerator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -221,6 +221,7 @@ void QELEventGenerator::ProcessEventRecord(GHepRecord * evrec) const
// reset bits
interaction->ResetBit(kISkipProcessChk);
interaction->ResetBit(kISkipKinematicChk);
interaction->ResetBit(kIAssumeFreeNucleon);

// get neutrino energy at struck nucleon rest frame and the
// struck nucleon mass (can be off the mass shell)
Expand Down
73 changes: 31 additions & 42 deletions src/Physics/QuasiElastic/XSection/NievesQELCCPXSec.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,16 +157,8 @@ double NievesQELCCPXSec::XSec(const Interaction * interaction,
double coulombFactor = 1.0;
double plLocal = leptonMom.P();

// Store the extra parameters needed to compute the contraction of the
// leptonic and hadronic tensors LmunuAnumu
bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
double r = target.HitNucPosition();
bool tgtIsNucleus = target.IsNucleus();
int tgt_pdgc = target.HitNucPdg();
int A = target.A();
int Z = target.Z();
int N = target.N();
bool hitNucIsProton = pdg::IsProton(target.HitNucPdg());

if ( fCoulomb ) {
// Coulomb potential
Expand Down Expand Up @@ -299,8 +291,8 @@ double NievesQELCCPXSec::XSec(const Interaction * interaction,
// the effective 4-momentum transfer q tilde (corrected for the nucleon
// binding energy and Coulomb effects)
double LmunuAnumuResult = LmunuAnumu(neutrinoMom, inNucleonMomOnShell,
leptonMom, qTildeP4, mNucleon, r, is_neutrino, tgtIsNucleus, tgt_pdgc,
A, Z, N, hitNucIsProton);
leptonMom, qTildeP4, mNucleon, is_neutrino, target,
interaction->TestBit( kIAssumeFreeNucleon ));

// Calculate xsec
double xsec = Gfactor*coulombFactor*LmunuAnumuResult;
Expand Down Expand Up @@ -493,14 +485,11 @@ void NievesQELCCPXSec::LoadConfig(void)
}
//___________________________________________________________________________
void NievesQELCCPXSec::CNCTCLimUcalc(TLorentzVector qTildeP4,
double M, double r, bool is_neutrino,
bool tgtIsNucleus, int tgt_pdgc,
int A, int Z, int N, bool hitNucIsProton,
double & CN, double & CT, double & CL,
double & imaginaryU,
double & t0, double & r00) const
double M, double r, bool is_neutrino, bool tgtIsNucleus, int tgt_pdgc,
int A, int Z, int N, bool hitNucIsProton, double & CN, double & CT, double & CL,
double & imaginaryU, double & t0, double & r00, bool assumeFreeNucleon) const
{
if ( tgtIsNucleus ) {
if ( tgtIsNucleus && !assumeFreeNucleon ) {
double dq = qTildeP4.Vect().Mag();
double dq2 = TMath::Power(dq,2);
double q2 = 1 * qTildeP4.Mag2();
Expand Down Expand Up @@ -569,7 +558,8 @@ void NievesQELCCPXSec::CNCTCLimUcalc(TLorentzVector qTildeP4,

CT = 1.0/TMath::Power(abs(1.0-relLinTot*Vt),2);
CL = 1.0/TMath::Power(abs(1.0-relLinTot*Vl),2);
}else{
}
else {
//Polarization Coefficients: all equal to 1.0 for free nucleon
CN = 1.0;
CT = 1.0;
Expand Down Expand Up @@ -900,14 +890,18 @@ int NievesQELCCPXSec::leviCivita(int input[]) const{
// expressions used here are valid in a frame in which the
// initial nucleus is at rest, and qTilde must be in the z direction.
double NievesQELCCPXSec::LmunuAnumu(const TLorentzVector neutrinoMom,
const TLorentzVector inNucleonMomOnShell,
const TLorentzVector leptonMom,
const TLorentzVector qTildeP4,
double M, double r, bool is_neutrino,
bool tgtIsNucleus,
int tgt_pdgc, int A, int Z, int N,
bool hitNucIsProton) const
const TLorentzVector inNucleonMomOnShell, const TLorentzVector leptonMom,
const TLorentzVector qTildeP4, double M, bool is_neutrino,
const Target& target, bool assumeFreeNucleon) const
{
double r = target.HitNucPosition();
bool tgtIsNucleus = target.IsNucleus();
int tgt_pdgc = target.Pdg();
int A = target.A();
int Z = target.Z();
int N = target.N();
bool hitNucIsProton = pdg::IsProton( target.HitNucPdg() );

const double k[4] = {neutrinoMom.E(),neutrinoMom.Px(),neutrinoMom.Py(),neutrinoMom.Pz()};
const double kPrime[4] = {leptonMom.E(),leptonMom.Px(),
leptonMom.Py(),leptonMom.Pz()};
Expand Down Expand Up @@ -946,16 +940,17 @@ double NievesQELCCPXSec::LmunuAnumu(const TLorentzVector neutrinoMom,

double t0,r00;
double CN=1.,CT=1.,CL=1.,imU=0;
CNCTCLimUcalc(qTildeP4,M,r,is_neutrino,tgtIsNucleus,
tgt_pdgc,A,Z,N,hitNucIsProton,
CN,CT,CL,imU,t0,r00);
if(imU > kASmallNum)
CNCTCLimUcalc(qTildeP4, M, r, is_neutrino, tgtIsNucleus,
tgt_pdgc, A, Z, N, hitNucIsProton, CN, CT, CL, imU,
t0, r00, assumeFreeNucleon);

if ( imU > kASmallNum )
return 0.;

if(! fRPA){
CN=1.0;
CT=1.0;
CL=1.0;
if ( !fRPA || assumeFreeNucleon ) {
CN = 1.0;
CT = 1.0;
CL = 1.0;
}

double tulin[4] = {0.,0.,0.,0.};
Expand Down Expand Up @@ -1134,7 +1129,7 @@ double NievesQELCCPXSec::LmunuAnumu(const TLorentzVector neutrinoMom,

// Since the real parts of A and L are both symmetric and the imaginary
// parts are antisymmetric, the contraction should be real
if(imag(sum) > kASmallNum)
if ( imag(sum) > kASmallNum )
LOG("Nieves",pWARN) << "Imaginary part of tensor contraction is nonzero "
<< "in QEL XSec, real(sum) = " << real(sum)
<< "imag(sum) = " << imag(sum);
Expand Down Expand Up @@ -1241,12 +1236,6 @@ void NievesQELCCPXSec::CompareNievesTensors(const Interaction* in)
double M = target.HitNucMass();
double ml = interaction->FSPrimLepton()->Mass();
bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
bool tgtIsNucleus = target.IsNucleus();
int tgt_pdgc = target.HitNucPdg();
int A = target.A();
int Z = target.Z();
int N = target.N();
bool hitNucIsProton = pdg::IsProton(target.HitNucPdg());

// Iterate over lepton energy (which then affects q, which is passed to
// LmunuAnumu using in and out NucleonMom
Expand Down Expand Up @@ -1299,8 +1288,8 @@ void NievesQELCCPXSec::CompareNievesTensors(const Interaction* in)
// TODO: apply Coulomb correction to 3-momentum transfer dq

fFormFactors.Calculate(interaction);
LmunuAnumu(neutrinoMom,inNucleonMomOnShell,leptonMom,qTildeP4,
M,r,is_neutrino,tgtIsNucleus,tgt_pdgc,A,Z,N,hitNucIsProton);
LmunuAnumu(neutrinoMom, inNucleonMomOnShell, leptonMom, qTildeP4,
M, is_neutrino, target, false);
}
return;
} // END TESTING CODE
Expand Down
17 changes: 6 additions & 11 deletions src/Physics/QuasiElastic/XSection/NievesQELCCPXSec.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,10 +113,9 @@ class NievesQELCCPXSec : public XSecAlgorithmI {
// variables. If target is not a nucleus, then CN, CN, and CL are all 1.0.
// r must be in units of fm.
void CNCTCLimUcalc(TLorentzVector qTildeP4, double M, double r,
bool is_neutrino, bool tgtIsNucleus, int tgt_pdgc,
int A, int Z, int N, bool hitNucIsProton,
double & CN, double & CT, double & CL,
double & imU, double & t0, double & r00) const;
bool is_neutrino, bool tgtIsNucleus, int tgt_pdgc, int A, int Z, int N,
bool hitNucIsProton, double & CN, double & CT, double & CL, double & imU,
double & t0, double & r00, bool assumeFreeNucleon) const;

//Equations to calculate the relativistic Lindhard function for Amunu
std::complex<double> relLindhardIm(double q0gev, double dqgev,
Expand All @@ -141,13 +140,9 @@ class NievesQELCCPXSec : public XSecAlgorithmI {
int leviCivita(int input[]) const;

double LmunuAnumu(const TLorentzVector neutrinoMom,
const TLorentzVector inNucleonMom,
const TLorentzVector leptonMom,
const TLorentzVector outNucleonMom,
double M, double r, bool is_neutrino,
bool tgtIsNucleus,
int tgt_pdgc, int A, int Z, int N,
bool hitNucIsProton) const;
const TLorentzVector inNucleonMom, const TLorentzVector leptonMom,
const TLorentzVector outNucleonMom, double M, bool is_neutrino,
const Target& target, bool assumeFreeNucleon) const;

// NOTE: THE FOLLOWING CODE IS FOR TESTING PURPOSES ONLY
// Used to print tensor elements and various inputs for comparison to Nieves'
Expand Down
14 changes: 0 additions & 14 deletions src/Physics/QuasiElastic/XSection/QELUtils.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,20 +248,6 @@ double genie::utils::CosTheta0Max(const genie::Interaction& interaction) {
return cos_theta0_max;
}

std::string genie::utils::VecToString(const TVector3& vec)
{
std::stringstream ss;
ss << '(' << vec.X() << ", " << vec.Y() << ", " << vec.Z() << ')';
return ss.str();
}

std::string genie::utils::LVecToString(const TLorentzVector& lvec)
{
std::stringstream ss;
ss << "E = " << lvec.E() << ", p = ";
return ss.str() + genie::utils::VecToString( lvec.Vect() );
}

void genie::utils::BindHitNucleon(genie::Interaction& interaction,
const genie::NuclearModelI& nucl_model, double& Eb,
genie::QELEvGen_BindingMode_t hitNucleonBindingMode)
Expand Down
7 changes: 5 additions & 2 deletions src/Physics/QuasiElastic/XSection/QELUtils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,19 @@ namespace genie {
} QELEvGen_BindingMode_t;

namespace utils {

double EnergyDeltaFunctionSolutionQEL(const Interaction& inter);

QELEvGen_BindingMode_t StringToQELBindingMode( const std::string& mode_str );

double ComputeFullQELPXSec(Interaction* interaction,
const NuclearModelI* nucl_model, const XSecAlgorithmI* xsec_model,
double cos_theta_0, double phi_0, double& Eb,
QELEvGen_BindingMode_t hitNucleonBindingMode, double min_angle_EM = 0.,
bool bind_nucleon = true);

double CosTheta0Max(const genie::Interaction& interaction);
std::string VecToString(const TVector3& vec);
std::string LVecToString(const TLorentzVector& lvec);

void BindHitNucleon(Interaction& interaction, const NuclearModelI& nucl_model,
double& Eb, QELEvGen_BindingMode_t hitNucleonBindingMode);
}
Expand Down

0 comments on commit e55a227

Please sign in to comment.