MEC utilities. More...

Namespaces
	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

Function Documentation

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 = ( 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;
 }

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;
 }

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;
 }

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;
 }

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;
 }

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);
 }

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.
 }

Variable Documentation

const double genie::utils::mec::Q0LimitMaxXSec = 2.

Definition at line 83 of file MECUtils.h.

const double genie::utils::mec::QMagLimitMaxXSec = 2.

Definition at line 86 of file MECUtils.h.

Namespaces

Functions

Variables

Detailed Description

Function Documentation

Variable Documentation