MEC utilities. More...

Namespaces
namespace	gsl

Functions
double	GetTmuCostFromq0q3 (double dq0, double dq3, double Enu, double lmass, double &tmu, double &cost, double &area)
bool	GetTlCostlFromq0q3 (double q0, double q3, double Enu, double ml, double &Tl, double &costl)
bool	Getq0q3FromTlCostl (double Tl, double costl, double Enu, double ml, double &q0, double &q3)
double	J (double q0, double q3, double Enu, double ml)
double	Qvalue (int targetpdg, int nupdg)
double	OldTensorContraction (int nupdg, int targetpdg, double Enu, double Ml, double Tl, double costhl, int tensorpdg, genie::HadronTensorType_t tensor_type, char *tensor_model)
double	GetMaxXSecTlctl (const XSecAlgorithmI &xsec_model, const Interaction &inter, const double tolerance=0.01, const double safety_factor=1.2, const int max_n_layers=100)

Variables
const double	Q0LimitMaxXSec = 2.
const double	QMagLimitMaxXSec = 2.

Detailed Description

MEC utilities.

Author

Created:\n

License:\n Copyright (c) 2003-2025, The GENIE Collaboration: For the full text of the license visit http://copyright.genie-mc.org

Function Documentation

◆ GetMaxXSecTlctl()

double genie::utils::mec::GetMaxXSecTlctl	(	const XSecAlgorithmI &	xsec_model,
		const Interaction &	inter,
		const double	tolerance = 0.01,
		const double	safety_factor = 1.2,
		const int	max_n_layers = 100 )

Definition at line 334 of file MECUtils.cxx.

{
  // Clone the input interaction so that we can modify the Tl and ctl values
  // without messing anything up
  Interaction* interaction = new Interaction( inter );
 
  // Choose the appropriate minimum Q^2 value based on the interaction
  // mode (this is important for EM interactions since the differential
  // cross section blows up as Q^2 --> 0)
  double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit
  if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics
    ::electromagnetic::kMinQ2Limit; // EM limit
 
  const double Enu = interaction->InitState().ProbeE( kRfLab );
  const double ProbeMass = interaction->InitState().Probe()->Mass();
  const double LepMass = interaction->FSPrimLepton()->Mass();
 
  // Determine the bounds for the current layer being scanned.
 
  // Maximum energy transfer to consider
  const double Q0Max = std::min( utils::mec::Q0LimitMaxXSec, Enu );
 
  // Maximum momentum transfer to consider
  const double QMagMax = std::min( utils::mec::QMagLimitMaxXSec, 2.*Enu );
 
  // Translate these bounds into limits for the lepton
  // kinetic energy and Q^2
  const double TMax = std::max( 0., Enu - LepMass );
  double CurTMin = std::max( 0., TMax - Q0Max );
  double CurTMax = TMax;
  double CurQ2Min = Q2min;
  // Overshoots the true maximum Q^2, but not so severely that it's expected to
  // be a problem.
  double CurQ2Max = QMagMax * QMagMax;
 
  // This is based on the technique used to compute the maximum differential
  // cross section for rejection sampling in QELEventGenerator. A brute-force
  // scan over the two-dimensional kPSTlctl phase space is used to find the
  // maximum cross section. Multiple layers are used to "zoom in" on the
  // vicinity of the maximum. Tighter and tighter layers are made until
  // the maximum number of iterations is reached or the xsec stabilizes
  // to within a given tolerance.
  const int N_T = 10; // number of kinetic energy steps per layer
  const int N_Q2 = 10; // number of scattering cosine steps per layer
  double T_at_xsec_max = CurTMax;
  double Q2_at_xsec_max = CurQ2Max;
  double cur_max_xsec = -1.;
 
  for ( int ilayer = 0; ilayer < max_n_layers; ++ilayer ) {
 
    double last_layer_xsec_max = cur_max_xsec;
 
    double T_increment = ( CurTMax - CurTMin ) / ( N_T - 1 );
    double Q2_increment   = ( CurQ2Max - CurQ2Min ) / ( N_Q2 - 1 );
 
    // Scan through the 2D phase space using the bounds of the current layer
    for ( int iT = 0; iT < N_T; ++iT ) {
      double T = CurTMin + iT * T_increment;
      for ( int iQ2 = 0; iQ2 < N_Q2; ++iQ2 ) {
        double Q2 = CurQ2Min + iQ2 * Q2_increment;
 
        // Compute the scattering cosine at fixed Tl given a value of Q^2.
        double Elep = T + LepMass;
        double pnu = std::sqrt( std::max(0., Enu*Enu - ProbeMass*ProbeMass) );
        double plep = std::sqrt( std::max(0., std::pow(T + LepMass, 2)
          - LepMass*LepMass) );
 
        double Costh = plep==0?0:( 2.*Enu*Elep - ProbeMass*ProbeMass - LepMass*LepMass
          - Q2 ) / ( 2. * pnu * plep );
        // Respect the bounds of the cosine function.
        Costh = std::min( std::max(-1., Costh), 1. );
 
        // Update the values of Tl and ctl in the interaction, then
        // compute the differential cross section
        interaction->KinePtr()->SetKV( kKVTl, T );
        interaction->KinePtr()->SetKV( kKVctl, Costh );
        double xsec = xsec_model.XSec( interaction, kPSTlctl );
 
        // If the current differential cross section exceeds the previous
        // maximum value, update the stored value and information about where it
        // lies in the kPSTlctl phase space.
        if ( xsec > cur_max_xsec ) {
          T_at_xsec_max = T;
          Q2_at_xsec_max = Q2;
          cur_max_xsec = xsec;
        }
 
      } // Done with cos(theta) scan
    } // Done with lepton kinetic energy scan
 
    LOG("MEC", pDEBUG) << "Layer " << ilayer << " has max xsec = "
      << cur_max_xsec << " for T = " << T_at_xsec_max << ", Q2 = "
      << Q2_at_xsec_max;
 
    // ** Calculate the range for the next layer **
 
    // Don't let the minimum kinetic energy fall below zero
    CurTMin = std::max( 0., T_at_xsec_max - T_increment );
 
    // We know a priori that the maximum cross section should occur near the
    // maximum lepton kinetic energy, so keep the old maximum value when
    // recalculating the new range. It was originally set to something near the
    // maximum allowed value, so this should zero in on the appropriate region
    // in future scans. The updated maximum from the old code is shown below,
    // but due to problems with round-off in the kinematic limits, it can
    // sometimes miss the edge of the phase space if the T_increment value is
    // coarse enough. - S. Gardiner, 1 July 2020
    CurTMax = TMax; // T_at_xsec_max + T_increment;
 
    // We use the same idea here. A priori, we know that the maximum cross section
    // should lie near the minimum Q^2 value. Round-off gives us the same problems
    // in finding the low Q^2 edge of phase space, so just keep the minimum equal
    // to the cut for the next scan. The old assignment (which mostly worked
    // but had problems for a coarse Q2_increment value) is commented out below.
    // - S. Gardiner, 1 July 2020
    CurQ2Min = Q2min; // std::max( Q2min, Q2_at_xsec_max - Q2_increment );
    CurQ2Max = Q2_at_xsec_max + Q2_increment;
 
    // If the xsec has stabilized within the requested tolerance, then
    // stop adding more layers
    double improvement_factor = cur_max_xsec / last_layer_xsec_max;
    if ( ilayer && std::abs(improvement_factor - 1.) < tolerance ) break;
 
  } // Done with iterations over layers
 
 
  // We're done. Delete the cloned interaction object before returning the
  // estimate of the maximum differential cross section.
  delete interaction;
 
  double XSecMax = cur_max_xsec * safety_factor;
  return XSecMax;
}

References genie::Interaction::FSPrimLepton(), genie::Interaction::InitState(), genie::ProcessInfo::IsEM(), genie::Interaction::KinePtr(), genie::kKVctl, genie::kKVTl, genie::controls::kMinQ2Limit, genie::kPSTlctl, genie::kRfLab, LOG, pDEBUG, genie::InitialState::Probe(), genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), Q0LimitMaxXSec, QMagLimitMaxXSec, genie::Kinematics::SetKV(), and genie::XSecAlgorithmI::XSec().

Referenced by genie::QELEventGeneratorSuSA::ComputeMaxXSec(), and genie::MECGenerator::SelectSuSALeptonKinematics().

◆ Getq0q3FromTlCostl()

bool genie::utils::mec::Getq0q3FromTlCostl	(	double	Tl,
		double	costl,
		double	Enu,
		double	ml,
		double &	q0,
		double &	q3 )

Definition at line 121 of file MECUtils.cxx.

{
  q0 = -999;
  q3 = -999;
 
  double El = Tl + ml;
  q0 = Enu - El;
  if(q0 < 0) {
     q0 = -999;
     q3 = -999;
     return false;
  }
 
  double pl = TMath::Sqrt(El * El - ml * ml);
  double q3sq = pl * pl + Enu * Enu - 2.0 * pl * Enu * costl;
  if(q3sq < 0) {
     q0 = -999;
     q3 = -999;
     return false;
  }
  q3 = TMath::Sqrt(q3sq);
 
  return true;
}

Referenced by genie::utils::mec::gsl::d2Xsec_dTCosth::DoEval(), genie::SuSAv2MECPXSec::PairRatio(), genie::QELEventGeneratorSuSA::SelectLeptonKinematics(), genie::MECGenerator::SelectNSVLeptonKinematics(), genie::SuSAv2MECPXSec::XSec(), and genie::SuSAv2QELPXSec::XSec().

◆ GetTlCostlFromq0q3()

bool genie::utils::mec::GetTlCostlFromq0q3	(	double	q0,
		double	q3,
		double	Enu,
		double	ml,
		double &	Tl,
		double &	costl )

Definition at line 82 of file MECUtils.cxx.

{
  Tl = Enu - q0 - ml;
  if(Tl < 0.) {
    costl = -999;
    Tl    = -999;
    return false;
  }
 
  double El = Tl + ml;
  double plsq = El * El - ml * ml;
  if(plsq < 0.) {
    costl = -999;
    Tl    = -999;
    return false;
  }
  double pl = TMath::Sqrt(plsq);
  double numerator =  Enu * Enu + pl * pl - q3 * q3;
  double denominator = 2.0 * pl * Enu;
  if(denominator <= 0.0) {
    costl = 0.0;
    if(denominator < 0.0) {
      return false;
    }
  }
  else {
    costl = numerator / denominator;
  }
 
  if(TMath::Abs(numerator) > TMath::Abs(denominator)){
    costl = -999;
    Tl    = -999;
    return false;
  }
 
  return true;
}

Referenced by genie::QELEventGeneratorSuSA::SelectLeptonKinematics().

◆ GetTmuCostFromq0q3()

double genie::utils::mec::GetTmuCostFromq0q3	(	double	dq0,
		double	dq3,
		double	Enu,
		double	lmass,
		double &	tmu,
		double &	cost,
		double &	area )

Definition at line 33 of file MECUtils.cxx.

{
 tmu = Enu - dq0 - lmass;
  if(tmu < 0.0){
    cost = -999;
    tmu = -999;
    area = 0.0;
    return -999;
  }
 
  double thisE = tmu + lmass;
  double thisp = sqrt( thisE * thisE - lmass * lmass);
  double numerator =  Enu * Enu + thisp * thisp - dq3 * dq3;
  double denominator = 2.0 * thisp * Enu;
  if(denominator <= 0.0 ){
    cost = 0.0;
    if(denominator < 0.0){
      return -999;
    }
  }
  else cost = numerator / denominator;
 
  if(TMath::Abs(numerator) > TMath::Abs(denominator)){
    cost = -999;
    tmu = -999;
    area = 0.0;
    return -999;
  }
 
  // xCrossSect is not yet in the right units for this particular case.
  // need areaElement to go dsigma/dTmudcost to dsigma/dq0dq3
  // Recompute the area element jacobian
 
  // dT/dq0 x dC/dq3 - dT/dq3 x dC/dq0
  double areaElement = 0.0;
  //double veryCloseToZero = 0.000001;  // in GeV, this should be safe.
  double insqrt = 0.0;
  numerator = 0.0;
  insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - lmass * lmass;
  numerator = dq3 / Enu;
  if(insqrt < 0.0)areaElement=0.0;
  else areaElement = numerator / TMath::Sqrt(insqrt);
  area = areaElement;
 
  return 0;
}

◆ J()

double genie::utils::mec::J	(	double	q0,
		double	q3,
		double	Enu,
		double	ml )

Definition at line 147 of file MECUtils.cxx.

{
  // dT/dq0 x dC/dq3 - dT/dq3 x dC/dq0
  double area = 0.0;
  double insqrt = 0.0;
  insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - ml * ml;
  double numerator = dq3 / Enu;
  if(insqrt < 0.0) {
     area=0.0;
  }
  else {
     area = numerator / TMath::Sqrt(insqrt);
  }
  return area;
}

◆ OldTensorContraction()

double genie::utils::mec::OldTensorContraction	(	int	nupdg,
		int	targetpdg,
		double	Enu,
		double	Ml,
		double	Tl,
		double	costhl,
		int	tensorpdg,
		genie::HadronTensorType_t	tensor_type,
		char *	tensor_model )

Definition at line 208 of file MECUtils.cxx.

{
  TLorentzVector v4lep;
  TLorentzVector v4Nu(0,0,Enu,Enu); // assuming traveling along z:
  TLorentzVector v4q;
  double q0nucleus;
  double facconv = 0.0389391289e15; // std::pow(0.19733,2)*1e15;
 
  double myQvalue = 0.0;
 
  // Angles
  double sinthl = 1. - costhl * costhl;
  if(sinthl < 0.0) sinthl = 0.0;
  else sinthl = TMath::Sqrt(sinthl);
 
  double Cosh = TMath::Cos(TMath::ACos(costhl)/2.);
  double Sinh = TMath::Sin(TMath::ACos(costhl)/2.);
 
  // Lepton
  v4lep.SetE( Tl + Ml );
  // energy transfer from the lepton
  double q0 = v4Nu.E() - v4lep.E();
  // energy transfer that actually gets to the nucleons
  q0nucleus = q0 - myQvalue;
 
  //Get q3
  double pl = TMath::Sqrt(v4lep.E() * v4lep.E() - Ml * Ml);
  double q3sq = pl * pl + Enu * Enu - 2.0 * pl * Enu * costhl;
  double q3 = sqrt(q3sq);
 
  // Define some calculation placeholders
  double part1, part2;
  double modkprime ;
  double w1, w2, w3, w4, w5;
  double wtotd[5];
 
  double xsec = 0;
  if (q0nucleus <= 0){
    //nothing transfered to nucleus thus no 2p2h
    xsec = 0.;
  } else {
    // energy was transfered to the nucleus
    modkprime = std::sqrt(TMath::Power(v4lep.E(),2)-TMath::Power(Ml,2));
    v4lep.SetX(modkprime*sinthl);
    v4lep.SetY(0);
    v4lep.SetZ(modkprime*costhl);
 
    //q: v4q = v4Nu - v4lep;
    v4q.SetE(q0nucleus);
    v4q.SetX(v4Nu.X() - v4lep.X());
    v4q.SetY(v4Nu.Y() - v4lep.Y());
    v4q.SetZ(v4Nu.Z() - v4lep.Z());
 
 
    // Get the appropriate hadron tensor model object
    const genie::HadronTensorModelI* ht_model
      = dynamic_cast<const genie::HadronTensorModelI*>(
      genie::AlgFactory::Instance()->GetAlgorithm( tensor_model, "Default" ));
 
    const LabFrameHadronTensorI* tensor_table
      = dynamic_cast<const LabFrameHadronTensorI*>( ht_model->GetTensor(targetpdg,
      tensor_type) );
 
    double W00 = (tensor_table->tt(q0nucleus,q3)).real();
    double W0Z = (tensor_table->tz(q0nucleus,q3)).real();
    double WXX = (tensor_table->xx(q0nucleus,q3)).real();
    double WXY = -1.0*(tensor_table->xy(q0nucleus,q3)).imag();
    double WZZ = (tensor_table->zz(q0nucleus,q3)).real();
 
    w1=WXX/2.;
    w2=(W00+WXX+(q0*q0/(v4q.Vect().Mag()*v4q.Vect().Mag())
                           *(WZZ-WXX))
          -(2.*q0/v4q.Vect().Mag()*W0Z))/2.;
    w3=WXY/v4q.Vect().Mag();
    w4=(WZZ-WXX)/(2.*v4q.Vect().Mag()*v4q.Vect().Mag());
    w5=(W0Z-(q0/v4q.Vect().Mag()*(WZZ-WXX)))/v4q.Vect().Mag();
    //w6 we have no need for w6, noted at the end of IIA.
 
    // adjust for anti neutrinos
 
    if (nupdg < 0) w3 = -1. * w3;
 
    // calculate cross section, in parts
    double xw1 = w1*costhl;
    double xw2 = w2/2.*costhl;
    double xw3 = w3/2.*(v4lep.E()+modkprime-(v4lep.E()+v4Nu.E())*costhl);
    double xw4 = w4/2.*(Ml*Ml*costhl+2.*v4lep.E()*(v4lep.E()+modkprime)*Sinh*Sinh);
    double xw5 = w5*(v4lep.E()+modkprime)/2.;
    part1 = xw1 - xw2 + xw3 + xw4 - xw5;
 
    double yw1 = 2.*w1*Sinh*Sinh;
    double yw2 = w2*Cosh*Cosh;
    double yw3 = w3*(v4lep.E()+v4Nu.E())*Sinh*Sinh;
    double yw4 = Ml*Ml*part1/(v4lep.E()*(v4lep.E()+modkprime));
    part2 = yw1 + yw2 - yw3 + yw4;
 
    xsec = modkprime*v4lep.E()*kGF2*2./kPi*part2*facconv;
 
    if( ! (xsec >= 0.0) ){
      // Traps negative numbers and nan's.
      // Sometimes costhl is just over 1.0 due to fluke or numerical
      // precision, so sinthl would be undefined.
      xsec = 0;
    }
  }
  if((Enu>1.15 && Enu<1.25) && (v4q.Vect().Mag()>0.615 && v4q.Vect().Mag()<0.625) &&  (q0<0.500 && q0>0.490))
    LOG("MECUtils", pFATAL)//SB
         << "Got xsec = " << xsec
         << "\n " << part1 << " " << part2
         << "\n w[] : " << w1 << ", " << w2 << ", " << w3 << ", " << w4 << ", " << w5
         << "\n wtotd[] : " << wtotd[0] << ", " << wtotd[1] << ", " << wtotd[2] << ", " << wtotd[3] << ", " << wtotd[4]
         << "\n vec " << v4q.Vect().Mag() << ", " << q0  << ", " << v4q.Px() << ", " << v4q.Py() << ", " << v4q.Pz()
         << "\n v4qX " << v4Nu.X() << ", " << v4lep.X() << ", " << costhl << ", " << sinthl << ", " << modkprime
   << "\n input " << Enu << ", " << Ml << ", " << Tl << ", " <<costhl << ", " << v4q.Vect().Mag() << ", " << q0;
 
  //LOG("MECUtils", pFATAL) << "xsec(Enu = " << Enu << " GeV, Ml = " << Ml << " GeV; " << "Tl = " << Tl << " GeV, costhl = " << costhl << ") = " << xsec << " x 1E-41 cm^2";
 
  //return (xsec * (1.0E-39 * units::cm2));
  return (xsec);
}

References genie::AlgFactory::GetAlgorithm(), genie::HadronTensorModelI::GetTensor(), genie::AlgFactory::Instance(), genie::constants::kGF2, genie::constants::kPi, LOG, pFATAL, genie::HadronTensorI::tt(), genie::HadronTensorI::tz(), genie::HadronTensorI::xx(), genie::HadronTensorI::xy(), and genie::HadronTensorI::zz().

◆ Qvalue()

double genie::utils::mec::Qvalue	(	int	targetpdg,
		int	nupdg )

Definition at line 164 of file MECUtils.cxx.

{
// The had tensor calculation requires parameters that describe the
// nuclear density. For the Valencia model they are taken by
// C. Garcia-Recio, Nieves, Oset Nucl.Phys.A 547 (1992) 473-487 Table I
// and simplifies that the p and n density are not necessarily the same?
// Standard tables such as deVries et al are similar ~5%
// This is the only nuclear input on the hadron side of the Hadron Tensor.
// and is what makes PseudoPb and PseudoFe different than Ni56 and Rf208
//
 
  int nearpdg = targetpdg;
 
  if(nupdg > 0){
    // get Z+1 for CC interaction e.g. 1000210400 from 1000200400
    nearpdg = targetpdg + 10000;
  } else {
    // get Z-1 for CC interaction e.g. 1000210400 from 1000200400
    nearpdg = targetpdg - 10000;
  }
 
  TParticlePDG *partf = PDGLibrary::Instance()->Find(nearpdg);
  TParticlePDG *parti = PDGLibrary::Instance()->Find(targetpdg);
 
  if(NULL == partf || NULL == parti){
    // maybe not every valid nucleus in the table has Z+1 or Z-1
    // for example, Ca40 did not.
    LOG("MECUtils", pFATAL)
      << "Can't get qvalue, nucleus " << targetpdg << " or " << nearpdg
      << " is not in the table of nuclei in /data/evgen/pdg ";
    exit(-1);
  }
 
  double massf = partf->Mass();
  double massi = parti->Mass();
 
  // keep this here as a reminder of what not to do
  // +/- emass;  // subtract electron from Z+1, add it to Z-1
  // the lookup table gives the nucleus mass, not the atomic
  // mass, so don't need this.
 
  return massf - massi;  // in GeV.
}