18 #include "TDatabasePDG.h" 19 #include "TLorentzVector.h" 40 TLorentzVector totMom = *k4 + *
p4;
42 TVector3
beta = totMom.BoostVector();
63 TLorentzVector total_p4 = (*probe) + hit_nucleon;
64 TVector3 beta_COM_to_lab = total_p4.BoostVector();
65 TVector3 beta_lab_to_COM = -beta_COM_to_lab;
68 double gamma2 =
std::pow(total_p4.Gamma(), 2);
71 TLorentzVector leptonCOM = TLorentzVector( lepton );
72 leptonCOM.Boost( beta_lab_to_COM );
75 double theta0 = leptonCOM.Angle( beta_COM_to_lab );
80 TVector3 lepton_vel = ( lepton.Vect() * (1. / lepton.E()) );
81 TVector3 outNucleon_vel = ( outNucleon.Vect() * (1. / outNucleon.E()) );
82 double vel_diff = ( lepton_vel - outNucleon_vel ).
Mag();
86 double factor =
std::sqrt( 1. + (1. - cos_theta0_squared)*(gamma2 - 1.) )
87 *
std::pow(leptonCOM.P(), 2) / vel_diff;
96 double cos_theta_0,
double phi_0,
double& Eb,
102 if ( bind_nucleon ) {
104 hitNucleonBindingMode);
111 TDatabasePDG *tb = TDatabasePDG::Instance();
119 if (
std::sqrt(s) < lepMass + mNf )
return 0.;
121 double outLeptonEnergy = ( s - mNf*mNf + lepMass*lepMass ) / (2 *
std::sqrt(s));
123 if (outLeptonEnergy*outLeptonEnergy - lepMass*lepMass < 0.)
return 0.;
124 double outMomentum =
TMath::Sqrt(outLeptonEnergy*outLeptonEnergy - lepMass*lepMass);
136 TVector3 lepton3Mom(0., 0., outMomentum);
137 lepton3Mom.SetTheta( TMath::ACos(cos_theta_0) );
138 lepton3Mom.SetPhi( phi_0 );
142 TVector3 zvec(0., 0., 1.);
143 TVector3 rot = ( zvec.Cross(beta) ).
Unit();
144 double angle = beta.Angle( zvec );
148 if ( beta.Perp() == 0. && beta.Z() < 0. ) {
149 rot = TVector3(0., 1., 0.);
155 lepton3Mom.Rotate(angle, rot);
159 TLorentzVector
lepton(lepton3Mom, outLeptonEnergy);
163 TLorentzVector outNucleon(-1*lepton.Px(),-1*lepton.Py(),-1*lepton.Pz(),
168 outNucleon.Boost(beta);
176 TLorentzVector qP4 = *nuP4 -
lepton;
178 double Q2 = -1 * qP4.Mag2();
187 if (Q2 < Q2lim.min || Q2 > Q2lim.
max)
return 0.;
200 if ( binding_mode ==
"UseGroundStateRemnant" ) {
203 else if ( binding_mode ==
"UseNuclearModel" ) {
206 else if ( binding_mode ==
"OnShell" ) {
210 LOG(
"QELEvent",
pFATAL) <<
"Unrecognized setting \"" << binding_mode
211 <<
"\" requested in genie::utils::StringToQELBindingMode()";
228 double gamma = 1. /
std::sqrt(1. - beta.Mag2());
233 double lepton_E_COM = (sqrt_s*sqrt_s + ml*ml - mNf*mNf) / (2.*sqrt_s);
246 double ENi = p4Ni.E();
250 double ENi_on_shell =
std::sqrt( mNi*mNi + p4Ni.Vect().Mag2() );
252 double epsilon_B = ENi_on_shell - ENi;
254 double cos_theta0_max = ( probe_E_lab - gamma*lepton_E_COM - epsilon_B )
255 / ( gamma * lepton_p_COM * beta.Mag() );
256 return cos_theta0_max;
270 TDatabasePDG* tb = TDatabasePDG::Instance();
285 double Mi = tgt->
Mass();
307 int Af = tgt->
A() - 1;
318 ENi = Mi -
std::sqrt( Mf*Mf + p3Ni.Mag2() );
336 p4Ni->SetVect( p3Ni );
Cross Section Calculation Interface.
T max(const caf::Proxy< T > &a, T b)
const KPhaseSpace & PhaseSpace(void) const
double RemovalEnergy(void) const
double Q2(const Interaction *const i)
int HitNucPdg(void) const
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
int RecoilNucleonPdg(void) const
recoil nucleon pdg
A simple [min,max] interval for doubles.
double HitNucMass(void) const
bool IsNucleus(void) const
const TLorentzVector & HadSystP4(void) const
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
double Mass(Resonance_t res)
resonance mass (GeV)
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)
Range1D_t Q2Lim(void) const
Q2 limits.
void BindHitNucleon(Interaction &interaction, const NuclearModelI &nucl_model, double &Eb, QELEvGen_BindingMode_t hitNucleonBindingMode)
enum genie::EQELEvGenBindingMode QELEvGen_BindingMode_t
const TVector3 & Momentum3(void) const
Summary information for an interaction.
const TLorentzVector & HitNucP4(void) const
const TLorentzVector & FSLeptonP4(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void SetFSLeptonP4(const TLorentzVector &p4)
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
const Kinematics & Kine(void) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
static const double kASmallNum
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
TLorentzVector * HitNucP4Ptr(void) const
static PDGLibrary * Instance(void)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
double CosTheta0Max(const genie::Interaction &interaction)
const UInt_t kIAssumeFreeNucleon
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
double CMEnergy() const
centre-of-mass energy (sqrt s)
Target * TgtPtr(void) const
TParticlePDG * Find(int pdgc)
int IonPdgCode(int A, int Z)
void SetHadSystP4(const TLorentzVector &p4)
InitialState * InitStatePtr(void) const
const InitialState & InitState(void) const
const ProcessInfo & ProcInfo(void) const
const Target & Tgt(void) const
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
double ProbeE(RefFrame_t rf) const
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
Initial State information.