13 #include <TLorentzVector.h> 23 using namespace genie;
28 double dq0,
double dq3,
double Enu,
double lmass,
29 double &tmu,
double &cost,
double &area)
31 tmu = Enu - dq0 - lmass;
39 double thisE = tmu + lmass;
40 double thisp =
sqrt( thisE * thisE - lmass * lmass);
41 double numerator = Enu * Enu + thisp * thisp - dq3 * dq3;
42 double denominator = 2.0 * thisp * Enu;
43 if(denominator <= 0.0 ){
45 if(denominator < 0.0){
49 else cost = numerator / denominator;
51 if(TMath::Abs(numerator) > TMath::Abs(denominator)){
63 double areaElement = 0.0;
67 insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - lmass * lmass;
68 numerator = dq3 / Enu;
69 if(insqrt < 0.0)areaElement=0.0;
77 double q0,
double q3,
double Enu,
double ml,
double& Tl,
double& costl)
87 double plsq = El * El - ml * ml;
94 double numerator = Enu * Enu + pl * pl - q3 * q3;
95 double denominator = 2.0 * pl * Enu;
96 if(denominator <= 0.0) {
98 if(denominator < 0.0) {
103 costl = numerator / denominator;
106 if(TMath::Abs(numerator) > TMath::Abs(denominator)){
116 double Tl,
double costl,
double Enu,
double ml,
double& q0,
double& q3)
130 double q3sq = pl * pl + Enu * Enu - 2.0 * pl * Enu * costl;
142 double dq0,
double dq3,
double Enu,
double ml)
147 insqrt = Enu * Enu - 2.0 * Enu * dq0 + dq0 * dq0 - ml * ml;
148 double numerator = dq3 / Enu;
169 int nearpdg = targetpdg;
173 nearpdg = targetpdg + 10000;
176 nearpdg = targetpdg - 10000;
182 if(NULL == partf || NULL == parti){
186 <<
"Can't get qvalue, nucleus " << targetpdg <<
" or " << nearpdg
187 <<
" is not in the table of nuclei in /data/evgen/pdg ";
191 double massf = partf->Mass();
192 double massi = parti->Mass();
199 return massf - massi;
215 nu_pdg, target_pdg, Enu, Ml, Tl, costhl, tensor_pdg, tensor_type);
219 int nupdg,
int targetpdg,
220 double Enu,
double Ml,
double Tl,
double costhl,
224 TLorentzVector v4lep;
225 TLorentzVector v4Nu(0,0,Enu,Enu);
228 double facconv = 0.0389391289e15;
233 double sinthl = 1. - costhl * costhl;
234 if(sinthl < 0.0) sinthl = 0.0;
237 double Cosh = TMath::Cos(TMath::ACos(costhl)/2.);
238 double Sinh = TMath::Sin(TMath::ACos(costhl)/2.);
241 v4lep.SetE( Tl + Ml );
243 double q0 = v4Nu.E() - v4lep.E();
245 q0nucleus = q0 - myQvalue;
250 double w1, w2, w3, w4, w5;
259 modkprime =
std::sqrt(TMath::Power(v4lep.E(),2)-TMath::Power(Ml,2));
260 v4lep.SetX(modkprime*sinthl);
262 v4lep.SetZ(modkprime*costhl);
266 v4q.SetX(v4Nu.X() - v4lep.X());
267 v4q.SetY(v4Nu.Y() - v4lep.Y());
268 v4q.SetZ(v4Nu.Z() - v4lep.Z());
271 const vector <genie::BLI2DNonUnifGrid *> &
272 tensor_table = hadtensor->
TensorTable(tensorpdg, tensor_type);
274 for (
int i=0 ;
i < 5;
i++){
275 wtotd[
i] = tensor_table[
i]->Evaluate(v4q.Vect().Mag(),v4q.E());
280 double W00 = wtotd[0];
281 double W0Z = wtotd[1];
282 double WXX = wtotd[2];
283 double WXY = wtotd[3];
284 double WZZ = wtotd[4];
287 w2=(W00+WXX+(q0*q0/(v4q.Vect().Mag()*v4q.Vect().Mag())
289 -(2.*q0/v4q.Vect().Mag()*W0Z))/2.;
290 w3=WXY/v4q.Vect().Mag();
291 w4=(WZZ-WXX)/(2.*v4q.Vect().Mag()*v4q.Vect().Mag());
292 w5=(W0Z-(q0/v4q.Vect().Mag()*(WZZ-WXX)))/v4q.Vect().Mag();
297 if (nupdg < 0) w3 = -1. * w3;
300 double xw1 = w1*costhl;
301 double xw2 = w2/2.*costhl;
302 double xw3 = w3/2.*(v4lep.E()+modkprime-(v4lep.E()+v4Nu.E())*costhl);
303 double xw4 = w4/2.*(Ml*Ml*costhl+2.*v4lep.E()*(v4lep.E()+modkprime)*Sinh*Sinh);
304 double xw5 = w5*(v4lep.E()+modkprime)/2.;
305 part1 = xw1 - xw2 + xw3 + xw4 - xw5;
307 double yw1 = 2.*w1*Sinh*Sinh;
308 double yw2 = w2*Cosh*Cosh;
309 double yw3 = w3*(v4lep.E()+v4Nu.E())*Sinh*Sinh;
310 double yw4 = Ml*Ml*part1/(v4lep.E()*(v4lep.E()+modkprime));
311 part2 = yw1 + yw2 - yw3 + yw4;
313 xsec = modkprime*v4lep.E()*
kGF2*2./
kPi*part2*facconv;
315 if( ! (xsec >= 0.0) ){
320 <<
"Got xsec = " << xsec
321 <<
"\n " << part1 <<
" " << part2
322 <<
"\n w[] : " << w1 <<
", " << w2 <<
", " << w3 <<
", " << w4 <<
", " << w5
323 <<
"\n wtotd[] : " << wtotd[0] <<
", " << wtotd[1] <<
", " << wtotd[2] <<
", " << wtotd[3] <<
", " << wtotd[4]
324 <<
"\n vec " << v4q.Vect().Mag() <<
", " << q0 <<
", " << v4q.Px() <<
", " << v4q.Py() <<
", " << v4q.Pz()
325 <<
"\n v4qX " << v4Nu.X() <<
", " << v4lep.X() <<
", " << costhl <<
", " << sinthl <<
", " << modkprime;
bool GetTlCostlFromq0q3(double q0, double q3, double Enu, double ml, double &Tl, double &costl)
double J(double q0, double q3, double Enu, double ml)
THE MAIN GENIE PROJECT NAMESPACE
static MECHadronTensor * Instance()
enum genie::MECHadronTensor::EMECHadronTensorType MECHadronTensorType_t
double Qvalue(int targetpdg, int nupdg)
Summary information for an interaction.
const vector< genie::BLI2DNonUnifGrid * > & TensorTable(int targetpdg, MECHadronTensor::MECHadronTensorType_t type)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
const Kinematics & Kine(void) const
Singleton class to load and store MEC hadron tensor tables, to aid in the implementation (and improve...
bool Getq0q3FromTlCostl(double Tl, double costl, double Enu, double ml, double &q0, double &q3)
double GetKV(KineVar_t kv) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static PDGLibrary * Instance(void)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
TParticlePDG * Find(int pdgc)
const InitialState & InitState(void) const
double GetTmuCostFromq0q3(double dq0, double dq3, double Enu, double lmass, double &tmu, double &cost, double &area)
const Target & Tgt(void) const
double TensorContraction(const Interaction *interaction, int tensor_pdg, MECHadronTensor::MECHadronTensorType_t tensor_type)
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...