NievesQELCCPXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Joe Johnston, University of Pittsburgh
8  Steven Dytman, University of Pittsburgh
9 
10  For the class documentation see the corresponding header file.
11 
12 */
13 //____________________________________________________________________________
14 
15 #include <TMath.h>
16 #include <TVector3.h>
17 #include <TLorentzVector.h>
18 #include <Math/IFunction.h>
19 #include <Math/Integrator.h>
20 #include <complex>
21 
46 
47 #include <iostream> // Used for testing code
48 #include <fstream> // Used for testing code
50 
51 using namespace genie;
52 using namespace genie::constants;
53 using namespace genie::controls;
54 using namespace genie::utils;
55 
56 //____________________________________________________________________________
58 XSecAlgorithmI("genie::NievesQELCCPXSec")
59 {
60 
61 }
62 //____________________________________________________________________________
64 XSecAlgorithmI("genie::NievesQELCCPXSec", config)
65 {
66 
67 }
68 //____________________________________________________________________________
70 {
71 
72  }
73 //____________________________________________________________________________
74 double NievesQELCCPXSec::XSec(const Interaction * interaction,
75  KinePhaseSpace_t kps) const
76 {
77  /*// TESTING CODE:
78  // The first time this method is called, output tensor elements and other
79  // kinmeatics variables for various kinematics. This can the be compared
80  // to Nieves' fortran code for validation purposes
81  if(fCompareNievesTensors){
82  LOG("Nieves",pNOTICE) << "Printing tensor elements for specific "
83  << "kinematics for testing purposes";
84  CompareNievesTensors(interaction);
85  fCompareNievesTensors = false;
86  exit(0);
87  }
88  // END TESTING CODE*/
89 
90 
91  if ( !this->ValidProcess (interaction) ) return 0.;
92  if ( !this->ValidKinematics(interaction) ) return 0.;
93 
94  // Get kinematics and init-state parameters
95  const Kinematics & kinematics = interaction -> Kine();
96  const InitialState & init_state = interaction -> InitState();
97  const Target & target = init_state.Tgt();
98 
99  // HitNucMass() looks up the PDGLibrary (on-shell) value for the initial
100  // struck nucleon
101  double mNi = target.HitNucMass();
102 
103  // Hadronic matrix element for CC neutrino interactions should really use
104  // the "nucleon mass," i.e., the mean of the proton and neutrino masses.
105  // This expression would also work for NC and EM scattering (since the
106  // initial and final on-shell nucleon masses would be the same)
107  double mNucleon = ( mNi + interaction->RecoilNucleon()->Mass() ) / 2.;
108 
109  // Create a copy of the struck nucleon 4-momentum that is forced
110  // to be on-shell (this will be needed later for the tensor contraction,
111  // in which the nucleon is treated in this way)
112  double inNucleonOnShellEnergy = std::sqrt( std::pow(mNi, 2)
113  + std::pow(target.HitNucP4().P(), 2) );
114 
115  // The Nieves CCQE model follows the de Forest prescription: free nucleon
116  // (i.e., on-shell) form factors and spinors are used, but an effective
117  // value of the 4-momentum transfer "qTilde" is used when computing the
118  // contraction of the hadronic tensor. See comments in the
119  // FullDifferentialXSec() method of LwlynSmithQELCCPXSec for more details.
120  TLorentzVector inNucleonMomOnShell( target.HitNucP4().Vect(),
121  inNucleonOnShellEnergy );
122 
123  // Get the four kinematic vectors and caluclate GFactor
124  // Create copies of all kinematics, so they can be rotated
125  // and boosted to the nucleon rest frame (because the tensor
126  // constraction below only applies for the initial nucleon
127  // at rest and q in the z direction)
128 
129  // All 4-momenta should already be stored, with the hit nucleon off-shell
130  // as appropriate
131  TLorentzVector* tempNeutrino = init_state.GetProbeP4(kRfLab);
132  TLorentzVector neutrinoMom = *tempNeutrino;
133  delete tempNeutrino;
134  TLorentzVector inNucleonMom = target.HitNucP4();
135  TLorentzVector leptonMom = kinematics.FSLeptonP4();
136  TLorentzVector outNucleonMom = kinematics.HadSystP4();
137 
138  // Apply Pauli blocking if enabled
139  if ( fDoPauliBlocking && target.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) {
140  int final_nucleon_pdg = interaction->RecoilNucleonPdg();
141  double kF = fPauliBlocker->GetFermiMomentum(target, final_nucleon_pdg,
142  target.HitNucPosition());
143  double pNf = outNucleonMom.P();
144  if ( pNf < kF ) return 0.;
145  }
146 
147  // Use the lab kinematics to calculate the Gfactor, in order to make
148  // the XSec differential in initial nucleon momentum and energy
149  // Divide by 4.0 because Nieves' conventions for the leptonic and hadronic
150  // tensor contraction differ from LwlynSmith by a factor of 4
151  double Gfactor = kGF2*fCos8c2 / (8.0*kPi*kPi*inNucleonOnShellEnergy
152  *neutrinoMom.E()*outNucleonMom.E()*leptonMom.E()) / 4.0;
153 
154  // Calculate Coulomb corrections
155  double ml = interaction->FSPrimLepton()->Mass();
156  double ml2 = TMath::Power(ml, 2);
157  double coulombFactor = 1.0;
158  double plLocal = leptonMom.P();
159 
160  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
161  double r = target.HitNucPosition();
162 
163  if ( fCoulomb ) {
164  // Coulomb potential
165  double Vc = vcr(& target, r);
166 
167  // Outgoing lepton energy and momentum including Coulomb potential
168  int sign = is_neutrino ? 1 : -1;
169  double El = leptonMom.E();
170  double pl = leptonMom.P();
171  double ElLocal = El - sign*Vc;
172 
173  if ( ElLocal - ml <= 0. ) {
174  LOG("Nieves", pDEBUG) << "Event should be rejected. Coulomb effects"
175  << " push kinematics below threshold. Returning xsec = 0.0";
176  return 0.0;
177  }
178 
179  // The Coulomb correction factor blows up as pl -> 0. To guard against
180  // unphysically huge corrections here, require that the lepton kinetic energy
181  // (at infinity) is larger than the magnitude of the Coulomb potential
182  // (should be around a few MeV)
183  double KEl = El - ml;
184  if ( KEl <= std::abs(Vc) ) {
185  LOG("Nieves", pDEBUG) << "Outgoing lepton has a very small kinetic energy."
186  << " Protecting against near-singularities in the Coulomb correction"
187  << " factor by returning xsec = 0.0";
188  return 0.0;
189  }
190 
191  // Local value of the lepton 3-momentum magnitude for the Coulomb
192  // correction
193  plLocal = TMath::Sqrt( ElLocal*ElLocal - ml2 );
194 
195  // Correction factor
196  coulombFactor= (plLocal * ElLocal) / (pl * El);
197 
198  }
199 
200  // When computing the contraction of the leptonic and hadronic tensors,
201  // we need to use an effective value of the 4-momentum transfer q.
202  // The energy transfer (q0) needs to be modified to account for the binding
203  // energy of the struck nucleon, while the 3-momentum transfer needs to
204  // be corrected for Coulomb effects.
205  //
206  // See the original Valencia model paper:
207  // https://journals.aps.org/prc/abstract/10.1103/PhysRevC.70.055503
208 
209  double q0Tilde = outNucleonMom.E() - inNucleonMomOnShell.E();
210 
211  // If binding energy effects pull us into an unphysical region, return
212  // zero for the differential cross section
213  if ( q0Tilde <= 0. && target.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) return 0.;
214 
215  // Note that we're working in the lab frame (i.e., the rest frame
216  // of the target nucleus). We can therefore use Nieves' explicit
217  // form of the Amunu tensor if we rotate the 3-momenta so that
218  // qTilde is in the +z direction
219  TVector3 neutrinoMom3 = neutrinoMom.Vect();
220  TVector3 leptonMom3 = leptonMom.Vect();
221 
222  TVector3 inNucleonMom3 = inNucleonMom.Vect();
223  TVector3 outNucleonMom3 = outNucleonMom.Vect();
224 
225  // If Coulomb corrections are being used, adjust the lepton 3-momentum used
226  // to get q3VecTilde so that its magnitude matches the local
227  // Coulomb-corrected value calculated earlier. Note that, although the
228  // treatment of Coulomb corrections by Nieves et al. doesn't change the
229  // direction of the lepton 3-momentum, it *does* change the direction of the
230  // 3-momentum transfer, and so the correction should be applied *before*
231  // rotating coordinates into a frame where q3VecTilde lies along the positive
232  // z axis.
233  TVector3 leptonMomCoulomb3 = (! fCoulomb ) ? leptonMom3
234  : plLocal * leptonMom3 * (1. / leptonMom3.Mag());
235  TVector3 q3VecTilde = neutrinoMom3 - leptonMomCoulomb3;
236 
237  // Find the rotation angle needed to put q3VecTilde along z
238  TVector3 zvec(0.0, 0.0, 1.0);
239  TVector3 rot = ( q3VecTilde.Cross(zvec) ).Unit(); // Vector to rotate about
240  // Angle between the z direction and q
241  double angle = zvec.Angle( q3VecTilde );
242 
243  // Handle the edge case where q3VecTilde is along -z, so the
244  // cross product above vanishes
245  if ( q3VecTilde.Perp() == 0. && q3VecTilde.Z() < 0. ) {
246  rot = TVector3(0., 1., 0.);
247  angle = kPi;
248  }
249 
250  // Rotate if the rotation vector is not 0
251  if ( rot.Mag() >= kASmallNum ) {
252 
253  neutrinoMom3.Rotate(angle,rot);
254  neutrinoMom.SetVect(neutrinoMom3);
255 
256  leptonMom3.Rotate(angle,rot);
257  leptonMom.SetVect(leptonMom3);
258 
259  inNucleonMom3.Rotate(angle,rot);
260  inNucleonMom.SetVect(inNucleonMom3);
261  inNucleonMomOnShell.SetVect(inNucleonMom3);
262 
263  outNucleonMom3.Rotate(angle,rot);
264  outNucleonMom.SetVect(outNucleonMom3);
265 
266  }
267 
268  // Calculate q and qTilde
269  TLorentzVector qP4 = neutrinoMom - leptonMom;
270  TLorentzVector qTildeP4(0., 0., q3VecTilde.Mag(), q0Tilde);
271 
272  double Q2 = -1. * qP4.Mag2();
273  double Q2tilde = -1. * qTildeP4.Mag2();
274 
275  // Store Q2tilde in the interaction so that we get the correct
276  // values of the form factors (according to the de Forest prescription)
277  interaction->KinePtr()->SetQ2(Q2tilde);
278 
279  double q2 = -Q2tilde;
280 
281  // Check that q2 < 0 (accounting for rounding errors)
282  if ( q2 >= kASmallNum ) {
283  LOG("Nieves", pWARN) << "q2 >= 0, returning xsec = 0.0";
284  return 0.0;
285  }
286 
287  // Calculate form factors
288  fFormFactors.Calculate( interaction );
289 
290  // Now that the form factors have been calculated, store Q2
291  // in the event instead of Q2tilde
292  interaction->KinePtr()->SetQ2( Q2 );
293 
294  // Do the contraction of the leptonic and hadronic tensors. See the
295  // RPA-corrected expressions for the hadronic tensor elements in appendix A
296  // of Phys. Rev. C 70, 055503 (2004). Note that the on-shell 4-momentum of
297  // the initial struck nucleon should be used in the calculation, as well as
298  // the effective 4-momentum transfer q tilde (corrected for the nucleon
299  // binding energy and Coulomb effects)
300  double LmunuAnumuResult = LmunuAnumu(neutrinoMom, inNucleonMomOnShell,
301  leptonMom, qTildeP4, mNucleon, is_neutrino, target,
302  interaction->TestBit( kIAssumeFreeNucleon ));
303 
304  // Calculate xsec
305  double xsec = Gfactor*coulombFactor*LmunuAnumuResult;
306 
307  // Apply the factor that arises from elimination of the energy-conserving
308  // delta function
309  xsec *= genie::utils::EnergyDeltaFunctionSolutionQEL( *interaction );
310 
311  // Apply given scaling factor
312  xsec *= fXSecScale;
313 
314  LOG("Nieves",pDEBUG) << "TESTING: RPA=" << fRPA
315  << ", Coulomb=" << fCoulomb
316  << ", q2 = " << q2 << ", xsec = " << xsec;
317 
318  //----- The algorithm computes dxsec/dQ2 or kPSQELEvGen
319  // Check whether variable tranformation is needed
320  if ( kps != kPSQELEvGen ) {
321 
322  // Compute the appropriate Jacobian for transformation to the requested
323  // phase space
324  double J = utils::kinematics::Jacobian(interaction, kPSQELEvGen, kps);
325 
326 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
327  LOG("Nieves", pDEBUG)
328  << "Jacobian for transformation to: "
329  << KinePhaseSpace::AsString(kps) << ", J = " << J;
330 #endif
331  xsec *= J;
332  }
333 
334  // Number of scattering centers in the target
335  int nucpdgc = target.HitNucPdg();
336  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
337 
338  xsec *= NNucl; // nuclear xsec
339 
340  return xsec;
341 }
342 //____________________________________________________________________________
344 {
345  // If we're using the new spline generation method (which integrates
346  // over the kPSQELEvGen phase space used by QELEventGenerator) then
347  // let the cross section integrator do all of the work. It's smart
348  // enough to handle free nucleon vs. nuclear targets, different
349  // nuclear models (including the local Fermi gas model), etc.
350  if ( fXSecIntegrator->Id().Name() == "genie::NewQELXSec" ) {
351  return fXSecIntegrator->Integrate(this, in);
352  }
353  else {
354  LOG("Nieves", pFATAL) << "Splines for the Nieves CCQE model must be"
355  << " generated using genie::NewQELXSec";
356  std::exit(1);
357  }
358 }
359 //____________________________________________________________________________
360 bool NievesQELCCPXSec::ValidProcess(const Interaction * interaction) const
361 {
362  if(interaction->TestBit(kISkipProcessChk)) return true;
363 
364  const InitialState & init_state = interaction->InitState();
365  const ProcessInfo & proc_info = interaction->ProcInfo();
366 
367  if(!proc_info.IsQuasiElastic()) return false;
368 
369  int nuc = init_state.Tgt().HitNucPdg();
370  int nu = init_state.ProbePdg();
371 
372  bool isP = pdg::IsProton(nuc);
373  bool isN = pdg::IsNeutron(nuc);
374  bool isnu = pdg::IsNeutrino(nu);
375  bool isnub = pdg::IsAntiNeutrino(nu);
376 
377  bool prcok = proc_info.IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
378  if(!prcok) return false;
379 
380  return true;
381 }
382 //____________________________________________________________________________
384 {
385  Algorithm::Configure(config);
386  this->LoadConfig();
387 }
388 //____________________________________________________________________________
390 {
391  Algorithm::Configure(config);
392  this->LoadConfig();
393 }
394 //____________________________________________________________________________
396 {
397  double thc;
398  GetParam( "CabibboAngle", thc ) ;
399  fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
400 
401  // Cross section scaling factor
402  GetParam( "QEL-CC-XSecScale", fXSecScale ) ;
403 
404  // hbarc for unit conversion, GeV*fm
406 
407  // load QEL form factors model
408  fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
409  this->SubAlg("FormFactorsAlg") );
411  fFormFactors.SetModel( fFormFactorsModel ); // <-- attach algorithm
412 
413  // load XSec Integrator
414  fXSecIntegrator = dynamic_cast<const XSecIntegratorI*>(
415  this->SubAlg("XSec-Integrator") );
417 
418  // Load settings for RPA and Coulomb effects
419 
420  // RPA corrections will not affect a free nucleon
421  GetParamDef("RPA", fRPA, true ) ;
422 
423  // Coulomb Correction- adds a correction factor, and alters outgoing lepton
424  // 3-momentum magnitude (but not direction)
425  // Correction only becomes sizeable near threshold and/or for heavy nuclei
426  GetParamDef( "Coulomb", fCoulomb, true ) ;
427 
428  LOG("Nieves", pNOTICE) << "RPA=" << fRPA << ", useCoulomb=" << fCoulomb;
429 
430  // Get nuclear model for use in Integral()
431  RgKey nuclkey = "IntegralNuclearModel";
432  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
434 
435  // Check if the model is a local Fermi gas
437 
438  if(!fLFG){
439  // get the Fermi momentum table for relativistic Fermi gas
440  GetParam( "FermiMomentumTable", fKFTableName ) ;
441 
442  fKFTable = 0;
444  fKFTable = kftp->GetTable( fKFTableName );
445  assert( fKFTable );
446  }
447 
448  // TESTING CODE
449  GetParamDef( "PrintDebugData", fCompareNievesTensors, false ) ;
450  // END TESTING CODE
451 
452  // Nuclear radius parameter (R = R0*A^(1/3)) to use when computing
453  // the maximum radius to use to integrate the Coulomb potential
454  GetParam("NUCL-R0", fR0) ; // fm
455 
456  std::string temp_mode;
457  GetParamDef( "RmaxMode", temp_mode, std::string("VertexGenerator") ) ;
458 
459  // Translate the string setting the Rmax mode to the appropriate
460  // enum value, or complain if one couldn't be found
461  if ( temp_mode == "VertexGenerator" ) {
463  }
464  else if ( temp_mode == "Nieves" ) {
466  }
467  else {
468  LOG("Nieves", pFATAL) << "Unrecognized setting \"" << temp_mode
469  << "\" requested for the RmaxMode parameter in the"
470  << " configuration for NievesQELCCPXSec";
471  gAbortingInErr = true;
472  std::exit(1);
473  }
474 
475  // Method to use to calculate the binding energy of the initial hit nucleon when
476  // generating splines
477  std::string temp_binding_mode;
478  GetParamDef( "IntegralNucleonBindingMode", temp_binding_mode, std::string("UseNuclearModel") );
480 
481  // Cutoff energy for integrating over nucleon momentum distribution (above this
482  // lab-frame probe energy, the effects of Fermi motion and binding energy
483  // are taken to be negligible for computing the total cross section)
484  GetParamDef("IntegralNuclearInfluenceCutoffEnergy", fEnergyCutOff, 2.5 ) ;
485 
486  // Get PauliBlocker for possible use in XSec()
487  fPauliBlocker = dynamic_cast<const PauliBlocker*>( this->SubAlg("PauliBlockerAlg") );
489 
490  // Decide whether or not it should be used in XSec()
491  GetParamDef( "DoPauliBlocking", fDoPauliBlocking, true );
492 }
493 //___________________________________________________________________________
494 void NievesQELCCPXSec::CNCTCLimUcalc(TLorentzVector qTildeP4,
495  double M, double r, bool is_neutrino, bool tgtIsNucleus, int tgt_pdgc,
496  int A, int Z, int N, bool hitNucIsProton, double & CN, double & CT, double & CL,
497  double & imaginaryU, double & t0, double & r00, bool assumeFreeNucleon) const
498 {
499  if ( tgtIsNucleus && !assumeFreeNucleon ) {
500  double dq = qTildeP4.Vect().Mag();
501  double dq2 = TMath::Power(dq,2);
502  double q2 = 1 * qTildeP4.Mag2();
503  //Terms for polarization coefficients CN,CT, and CL
504  double hbarc2 = TMath::Power(fhbarc,2);
505  double c0 = 0.380/fhbarc;//Constant for CN in natural units
506 
507  //Density gives the nuclear density, normalized to 1
508  //Input radius r must be in fm
509  double rhop = nuclear::Density(r,A)*Z;
510  double rhon = nuclear::Density(r,A)*N;
511  double rho = rhop + rhon;
512  double rho0 = A*nuclear::Density(0,A);
513 
514  double fPrime = (0.33*rho/rho0+0.45*(1-rho/rho0))*c0;
515 
516  // Get Fermi momenta
517  double kF1, kF2;
518  if(fLFG){
519  if(hitNucIsProton){
520  kF1 = TMath::Power(3*kPi2*rhop, 1.0/3.0) *fhbarc;
521  kF2 = TMath::Power(3*kPi2*rhon, 1.0/3.0) *fhbarc;
522  }else{
523  kF1 = TMath::Power(3*kPi2*rhon, 1.0/3.0) *fhbarc;
524  kF2 = TMath::Power(3*kPi2*rhop, 1.0/3.0) *fhbarc;
525  }
526  }else{
527  if(hitNucIsProton){
528  kF1 = fKFTable->FindClosestKF(tgt_pdgc, kPdgProton);
529  kF2 = fKFTable->FindClosestKF(tgt_pdgc, kPdgNeutron);
530  }else{
531  kF1 = fKFTable->FindClosestKF(tgt_pdgc, kPdgNeutron);
532  kF2 = fKFTable->FindClosestKF(tgt_pdgc, kPdgProton);
533  }
534  }
535 
536  double kF = TMath::Power(1.5*kPi2*rho, 1.0/3.0) *fhbarc;
537 
538  std::complex<double> imU(relLindhardIm(qTildeP4.E(),dq,kF1,kF2,
539  M,is_neutrino,t0,r00));
540 
541  imaginaryU = imag(imU);
542 
543  std::complex<double> relLin(0,0),udel(0,0);
544 
545  // By comparison with Nieves' fortran code
546  if(imaginaryU < 0.){
547  relLin = relLindhard(qTildeP4.E(),dq,kF,M,is_neutrino,imU);
548  udel = deltaLindhard(qTildeP4.E(),dq,rho,kF);
549  }
550  std::complex<double> relLinTot(relLin + udel);
551 
552  /* CRho = 2
553  DeltaRho = 2500 MeV, (2.5 GeV)^2 = 6.25 GeV^2
554  mRho = 770 MeV, (0.770 GeV)^2 = 0.5929 GeV^2
555  g' = 0.63 */
556  double Vt = 0.08*4*kPi/kPionMass2 *
557  (2* TMath::Power((6.25-0.5929)/(6.25-q2),2)*dq2/(q2-0.5929) + 0.63);
558  /* f^2/4/Pi = 0.08
559  DeltaSubPi = 1200 MeV, (1.2 GeV)^2 = 1.44 GeV^2
560  g' = 0.63 */
561  double Vl = 0.08*4*kPi/kPionMass2 *
562  (TMath::Power((1.44-kPionMass2)/(1.44-q2),2)*dq2/(q2-kPionMass2)+0.63);
563 
564  CN = 1.0/TMath::Power(abs(1.0-fPrime*relLin/hbarc2),2);
565 
566  CT = 1.0/TMath::Power(abs(1.0-relLinTot*Vt),2);
567  CL = 1.0/TMath::Power(abs(1.0-relLinTot*Vl),2);
568  }
569  else {
570  //Polarization Coefficients: all equal to 1.0 for free nucleon
571  CN = 1.0;
572  CT = 1.0;
573  CL = 1.0;
574  imaginaryU = 0.0;
575  }
576 }
577 //____________________________________________________________________________
578 // Gives the imaginary part of the relativistic lindhard function in GeV^2
579 // and sets the values of t0 and r00
580 std::complex<double> NievesQELCCPXSec::relLindhardIm(double q0, double dq,
581  double kFn, double kFp,
582  double M,
583  bool isNeutrino,
584  double & t0,
585  double & r00) const
586 {
587  double M2 = TMath::Power(M,2);
588  double EF1,EF2;
589  if(isNeutrino){
590  EF1 = TMath::Sqrt(M2+TMath::Power(kFn,2)); //EFn
591  EF2 = TMath::Sqrt(M2+TMath::Power(kFp,2)); //EFp
592  }else{
593  EF1 = TMath::Sqrt(M2+TMath::Power(kFp,2)); //EFp
594  EF2 = TMath::Sqrt(M2+TMath::Power(kFn,2)); //EFn
595  }
596 
597  double q2 = TMath::Power(q0,2) - TMath::Power(dq,2);
598  double a = (-q0+dq*TMath::Sqrt(1-4.0*M2/q2))/2.0;
599  double epsRP = TMath::Max(TMath::Max(M,EF2-q0),a);
600 
601  // Other theta functions for q are handled by nuclear suppression
602  // That is, q0>0 and -q2>0 are always handled, and q0>EF2-EF1 is
603  // handled if pauli blocking is on, because otherwise the final
604  // nucleon would be below the fermi sea
605  //if(fNievesSuppression && !interaction->TestBit(kIAssumeFreeNucleon )
606  //&& !EF1-epsRP<0){
607  //LOG("Nieves", pINFO) << "Average value of E(p) above Fermi sea";
608  //return 0;
609  //}else{
610  t0 = 0.5*(EF1+epsRP);
611  r00 = (TMath::Power(EF1,2)+TMath::Power(epsRP,2)+EF1*epsRP)/3.0;
612  std::complex<double> result(0.0,-M2/2.0/kPi/dq*(EF1-epsRP));
613  return result;
614  //}
615 }
616 //____________________________________________________________________________
617 //Following obtained from fortran code by J Nieves, which contained the following comment:
618 /*
619  NUCLEON relativistic Lindhard Function
620  Same normalization as ULIN
621  Real part
622  taken from Eur.Phys.J.A25:299-318,2005 (Barbaro et al)
623  Eq. 61
624 
625  Im. part: Juan.
626 
627  INPUT: Real*8
628  q0:Energy [fm]
629  qm: modulus 3mom [fm]
630  kf: Fermi mom [fm]
631 
632  OUTPUT: Complex*16 [fm]
633 
634  USES: ruLinRelX, relLindhardIm
635  */
636 //Takes inputs in GeV (with imU in GeV^2), and gives output in GeV^2
637 std::complex<double> NievesQELCCPXSec::relLindhard(double q0gev,
638  double dqgev, double kFgev, double M,
639  bool isNeutrino,
640  std::complex<double> /* relLindIm */) const
641 {
642  double q0 = q0gev/fhbarc;
643  double qm = dqgev/fhbarc;
644  double kf = kFgev/fhbarc;
645  double m = M/fhbarc;
646 
647  if(q0>qm){
648  LOG("Nieves", pWARN) << "relLindhard() failed";
649  return 0.0;
650  }
651 
652  std::complex<double> RealLinRel(ruLinRelX(q0,qm,kf,m)+ruLinRelX(-q0,qm,kf,m));
653  double t0,r00;
654  std::complex<double> ImLinRel(relLindhardIm(q0gev,dqgev,kFgev,kFgev,M,isNeutrino,t0,r00));
655  //Units of GeV^2
656  return(RealLinRel*TMath::Power(fhbarc,2) + 2.0*ImLinRel);
657 }
658 //____________________________________________________________________________
659 //Inputs assumed to be in natural units
660 std::complex<double> NievesQELCCPXSec::ruLinRelX(double q0, double qm,
661  double kf, double m) const
662 {
663  double q02 = TMath::Power(q0, 2);
664  double qm2 = TMath::Power(qm, 2);
665  double kf2 = TMath::Power(kf, 2);
666  double m2 = TMath::Power(m, 2);
667  double m4 = TMath::Power(m, 4);
668 
669  double ef = TMath::Sqrt(m2+kf2);
670  double q2 = q02-qm2;
671  double q4 = TMath::Power(q2,2);
672  double ds = TMath::Sqrt(1.0-4.0*m2/q2);
673  double L1 = log((kf+ef)/m);
674  std::complex<double> uL2(
675  TMath::Log(TMath::Abs(
676  (ef + q0 - TMath::Sqrt(m2+TMath::Power(kf-qm,2)))/
677  (ef + q0 - TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))) +
678  TMath::Log(TMath::Abs(
679  (ef + q0 + TMath::Sqrt(m2 + TMath::Power(kf - qm,2)))/
680  (ef + q0 + TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))));
681 
682  std::complex<double> uL3(
683  TMath::Log(TMath::Abs((TMath::Power(2*kf + q0*ds,2)-qm2)/
684  (TMath::Power(2*kf - q0*ds,2)-qm2))) +
685  TMath::Log(TMath::Abs((TMath::Power(kf-ef*ds,2) - (4*m4*qm2)/q4)/
686  (TMath::Power(kf+ef*ds,2) - (4*m4*qm2)/q4))));
687 
688  std::complex<double> RlinrelX(-L1/(16.0*kPi2)+
689  uL2*(2.0*ef+q0)/(32.0*kPi2*qm)-
690  uL3*ds/(64.0*kPi2));
691 
692  return RlinrelX*16.0*m2;
693 }
694 //____________________________________________________________________________
695 //Following obtained from fortran code by J Nieves, which contained the following comment:
696 /*
697  complex Lindhard function for symmetric nuclear matter:
698  from Appendix of
699  E.Oset et al Phys. Rept. 188:79, 1990
700  formula A.4
701 
702  input variables:
703  q_zero [fm^-1] : Energy
704  q_mod [fm^-1] : Momentum
705  rho [fm^3] : Nuclear density
706  k_fermi[fm^-1] : Fermi momentum
707 
708  All variables are real*8
709 
710  output variable:
711  delta_lind [fm^-2]
712 
713  ATTENTION!!!
714  Only works properly for real q_zero,
715  if q_zero has an imaginary part calculates the L. function
716  assuming Gamma= 0.
717  Therefore this subroutine provides two different functions
718  depending on whether q_zero is real or not!!!!!!!!!!!
719 */
720 std::complex<double> NievesQELCCPXSec::deltaLindhard(double q0,
721  double dq, double rho, double kF) const
722 {
723  double q_zero = q0/fhbarc;
724  double q_mod = dq/fhbarc;
725  double k_fermi = kF/fhbarc;
726  //Divide by hbarc in order to use natural units (rho is already in the correct units)
727 
728  //m = 939/197.3, md = 1232/197.3, mpi = 139/197.3
729  double m = 4.7592;
730  double md = 6.2433;
731  double mpi = 0.7045;
732 
733  double fdel_f = 2.13;
734  double wr = md-m;
735  double gamma = 0;
736  double gammap = 0;
737 
738  double q_zero2 = TMath::Power(q_zero, 2);
739  double q_mod2 = TMath::Power(q_mod, 2);
740  double k_fermi2 = TMath::Power(k_fermi, 2);
741 
742  double m2 = TMath::Power(m, 2);
743  double m4 = TMath::Power(m, 4);
744  double mpi2 = TMath::Power(mpi, 2);
745  double mpi4 = TMath::Power(mpi, 4);
746 
747  double fdel_f2 = TMath::Power(fdel_f, 2);
748 
749  //For the current code q_zero is always real
750  //If q_zero can have an imaginary part then only the real part is used
751  //until z and zp are calculated
752 
753  double s = m2+q_zero2-q_mod2+
754  2.0*q_zero *TMath::Sqrt(m2+3.0/5.0*k_fermi2);
755 
756  if(s>TMath::Power(m+mpi,2)){
757  double srot = TMath::Sqrt(s);
758  double qcm = TMath::Sqrt(TMath::Power(s,2)+mpi4+m4-2.0*(s*mpi2+s*m2+
759  mpi2*m2)) /(2.0*srot);
760  gamma = 1.0/3.0 * 1.0/(4.0*kPi) * fdel_f2*
761  TMath::Power(qcm,3)/srot*(m+TMath::Sqrt(m2+TMath::Power(qcm,2)))/mpi2;
762  }
763  double sp = m2+q_zero2-q_mod2-
764  2.0*q_zero *TMath::Sqrt(m2+3.0/5.0*k_fermi2);
765 
766 
767  if(sp > TMath::Power(m+mpi,2)){
768  double srotp = TMath::Sqrt(sp);
769  double qcmp = TMath::Sqrt(TMath::Power(sp,2)+mpi4+m4-2.0*(sp*mpi2+sp*m2+
770  mpi2*m2))/(2.0*srotp);
771  gammap = 1.0/3.0 * 1.0/(4.0*kPi) * fdel_f2*
772  TMath::Power(qcmp,3)/srotp*(m+TMath::Sqrt(m2+TMath::Power(qcmp,2)))/mpi2;
773  }
774  //}//End if statement
775  const std::complex<double> iNum(0,1.0);
776 
777  std::complex<double> z(md/(q_mod*k_fermi)*(q_zero-q_mod2/(2.0*md)
778  -wr +iNum*gamma/2.0));
779  std::complex<double> zp(md/(q_mod*k_fermi)*(-q_zero-q_mod2/(2.0*md)
780  -wr +iNum*gammap/2.0));
781 
782  std::complex<double> pzeta(0.0);
783  if(abs(z) > 50.0){
784  pzeta = 2.0/(3.0*z)+2.0/(15.0*z*z*z);
785  }else if(abs(z) < TMath::Power(10.0,-2)){
786  pzeta = 2.0*z-2.0/3.0*z*z*z-iNum*kPi/2.0*(1.0-z*z);
787  }else{
788  pzeta = z + (1.0-z*z) * log((z+1.0)/(z-1.0))/2.0;
789  }
790 
791  std::complex<double> pzetap(0);
792  if(abs(zp) > 50.0){
793  pzetap = 2.0/(3.0*zp)+2.0/(15.0*zp*zp*zp);
794  }else if(abs(zp) < TMath::Power(10.0,-2)){
795  pzetap = 2.0*zp-2.0/3.0*zp*zp*zp-iNum*kPi/2.0*(1.0-zp*zp);
796  }else{
797  pzetap = zp+ (1.0-zp*zp) * log((zp+1.0)/(zp-1.0))/2.0;
798  }
799 
800  //Multiply by hbarc^2 to give answer in units of GeV^2
801  return 2.0/3.0 * rho * md/(q_mod*k_fermi) * (pzeta +pzetap) * fdel_f2 *
802  TMath::Power(fhbarc,2);
803 }
804 
805 //____________________________________________________________________________
806 // Gives coulomb potential in units of GeV
807 double NievesQELCCPXSec::vcr(const Target * target, double Rcurr) const{
808  if(target->IsNucleus()){
809  int A = target->A();
810  int Z = target->Z();
811  double Rmax = 0.;
812 
813  if ( fCoulombRmaxMode == kMatchNieves ) {
814  // Rmax calculated using formula from Nieves' fortran code and default
815  // charge and neutron matter density parameters from NuclearUtils.cxx
816  if (A > 20) {
817  double c = TMath::Power(A,0.35), z = 0.54;
818  Rmax = c + 9.25*z;
819  }
820  else {
821  // c = 1.75 for A <= 20
822  Rmax = TMath::Sqrt(20.0)*1.75;
823  }
824  }
826  // TODO: This solution is fragile. If the formula used by VertexGenerator
827  // changes, then this one will need to change too. Switch to using
828  // a common function to get Rmax for both.
829  Rmax = 3. * fR0 * std::pow(A, 1./3.);
830  }
831  else {
832  LOG("Nieves", pFATAL) << "Unrecognized setting for fCoulombRmaxMode encountered"
833  << " in NievesQELCCPXSec::vcr()";
834  gAbortingInErr = true;
835  std::exit(1);
836  }
837 
838  if(Rcurr >= Rmax){
839  LOG("Nieves",pNOTICE) << "Radius greater than maximum radius for coulomb corrections."
840  << " Integrating to max radius.";
841  Rcurr = Rmax;
842  }
843 
844  ROOT::Math::IBaseFunctionOneDim * func = new
848 
849  double abstol = 1; // We mostly care about relative tolerance;
850  double reltol = 1E-4;
851  int nmaxeval = 100000;
852  ROOT::Math::Integrator ig(*func,ig_type,abstol,reltol,nmaxeval);
853  double result = ig.Integral(0,Rmax);
854  delete func;
855 
856  // Multiply by Z to normalize densities to number of protons
857  // Multiply by hbarc to put result in GeV instead of fm
858  return -kAem*4*kPi*result*fhbarc;
859  }else{
860  // If target is not a nucleus the potential will be 0
861  return 0.0;
862  }
863 }
864 //____________________________________________________________________________
865 int NievesQELCCPXSec::leviCivita(int input[]) const{
866  int copy[4] = {input[0],input[1],input[2],input[3]};
867  int permutations = 0;
868  int temp;
869 
870  for(int i=0;i<4;i++){
871  for(int j=i+1;j<4;j++){
872  //If any two elements are equal return 0
873  if(input[i] == input[j])
874  return 0;
875  //If a larger number is before a smaller one, use permutations
876  //(exchanges of two adjacent elements) to move the smaller element
877  //so it is before the larger element, eg 2341->2314->2134->1234
878  if(copy[i]>copy[j]){
879  temp = copy[j];
880  for(int k=j;k>i;k--){
881  copy[k] = copy[k-1];
882  permutations++;
883  }
884  copy[i] = temp;
885  }
886  }
887  }
888 
889  if(permutations % 2 == 0){
890  return 1;
891  }else{
892  return -1;
893  }
894 }
895 //____________________________________________________________________________
896 // Calculates the constraction of the leptonic and hadronic tensors. The
897 // expressions used here are valid in a frame in which the
898 // initial nucleus is at rest, and qTilde must be in the z direction.
899 double NievesQELCCPXSec::LmunuAnumu(const TLorentzVector neutrinoMom,
900 const TLorentzVector inNucleonMomOnShell, const TLorentzVector leptonMom,
901 const TLorentzVector qTildeP4, double M, bool is_neutrino,
902 const Target& target, bool assumeFreeNucleon) const
903 {
904  double r = target.HitNucPosition();
905  bool tgtIsNucleus = target.IsNucleus();
906  int tgt_pdgc = target.Pdg();
907  int A = target.A();
908  int Z = target.Z();
909  int N = target.N();
910  bool hitNucIsProton = pdg::IsProton( target.HitNucPdg() );
911 
912  const double k[4] = {neutrinoMom.E(),neutrinoMom.Px(),neutrinoMom.Py(),neutrinoMom.Pz()};
913  const double kPrime[4] = {leptonMom.E(),leptonMom.Px(),
914  leptonMom.Py(),leptonMom.Pz()};
915 
916  double q2 = qTildeP4.Mag2();
917 
918  const double q[4] = {qTildeP4.E(),qTildeP4.Px(),qTildeP4.Py(),qTildeP4.Pz()};
919  double q0 = q[0];
920  double dq = TMath::Sqrt(TMath::Power(q[1],2)+
921  TMath::Power(q[2],2)+TMath::Power(q[3],2));
922 
923  int sign = (is_neutrino) ? 1 : -1;
924 
925  // Get the QEL form factors (were calculated before this method was called)
926  double F1V = 0.5*fFormFactors.F1V();
927  double xiF2V = 0.5*fFormFactors.xiF2V();
928  double FA = -fFormFactors.FA();
929  // According to Nieves' paper, Fp = 2.0*M*FA/(kPionMass2-q2), but Llewelyn-
930  // Smith uses Fp = 2.0*M^2*FA/(kPionMass2-q2), so I divide by M
931  // This gives units of GeV^-1
932  double Fp = -1.0/M*fFormFactors.Fp();
933 
934 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
935  LOG("Nieves", pDEBUG) << "\n" << fFormFactors;
936 #endif
937 
938  // Calculate auxiliary parameters
939  double M2 = TMath::Power(M, 2);
940  double FA2 = TMath::Power(FA, 2);
941  double Fp2 = TMath::Power(Fp, 2);
942  double F1V2 = TMath::Power(F1V, 2);
943  double xiF2V2 = TMath::Power(xiF2V, 2);
944  double q02 = TMath::Power(q[0], 2);
945  double dq2 = TMath::Power(dq, 2);
946  double q4 = TMath::Power(q2, 2);
947 
948  double t0,r00;
949  double CN=1.,CT=1.,CL=1.,imU=0;
950  CNCTCLimUcalc(qTildeP4, M, r, is_neutrino, tgtIsNucleus,
951  tgt_pdgc, A, Z, N, hitNucIsProton, CN, CT, CL, imU,
952  t0, r00, assumeFreeNucleon);
953 
954  if ( imU > kASmallNum )
955  return 0.;
956 
957  if ( !fRPA || assumeFreeNucleon ) {
958  CN = 1.0;
959  CT = 1.0;
960  CL = 1.0;
961  }
962 
963  double tulin[4] = {0.,0.,0.,0.};
964  double rulin[4][4] = { {0.,0.,0.,0.},
965  {0.,0.,0.,0.},
966  {0.,0.,0.,0.},
967  {0.,0.,0.,0.} };
968 
969  // TESTING CODE:
971  // Use average values for initial momentum to calculate A, as given
972  // in Appendix B of Nieves' paper. T gives average values of components
973  // of p, and R gives the average value of two components multiplied
974  // together
975  double t3 = (0.5*q2 + q0*t0)/dq; // Average pz
976 
977  // Vector of p
978 
979  tulin[0] = t0;
980  tulin[3] = t3;
981 
982  // R is a 4x4 matrix, with R[mu][nu] is the average
983  // value of p[mu]*p[nu]
984  double aR = r00-M2;
985  double bR = (q4+4.0*r00*q02+4.0*q2*q0*t0)/(4.0*dq2);
986  double gamma = (aR-bR)/2.0;
987  double delta = (-aR+3.0*bR)/2.0/dq2;
988 
989  double r03 = (0.5*q2*t0 + q0*r00)/dq; // Average E(p)*pz
990 
991  rulin[0][0] = r00;
992  rulin[0][3] = r03;
993  rulin[1][1] = gamma;
994  rulin[2][2] = gamma;
995  rulin[3][0] = r03;
996  rulin[3][3] = gamma+delta*dq2; // END TESTING CODE
997  }
998  else {
999  // For normal code execulation, tulin is the initial nucleon momentum
1000  tulin[0] = inNucleonMomOnShell.E();
1001  tulin[1] = inNucleonMomOnShell.Px();
1002  tulin[2] = inNucleonMomOnShell.Py();
1003  tulin[3] = inNucleonMomOnShell.Pz();
1004 
1005  for(int i=0; i<4; i++)
1006  for(int j=0; j<4; j++)
1007  rulin[i][j] = tulin[i]*tulin[j];
1008  }
1009 
1010  //Additional constants and variables
1011  const int g[4][4] = {{1,0,0,0},{0,-1,0,0},{0,0,-1,0},{0,0,0,-1}};
1012  const std::complex<double> iNum(0,1);
1013  int leviCivitaIndexArray[4];
1014  double imaginaryPart = 0;
1015 
1016  std::complex<double> sum(0.0,0.0);
1017 
1018  double kPrimek = k[0]*kPrime[0]-k[1]*kPrime[1]-k[2]*kPrime[2]-k[3]*kPrime[3];
1019 
1020  std::complex<double> Lmunu(0.0,0.0),Lnumu(0.0,0.0);
1021  std::complex<double> Amunu(0.0,0.0),Anumu(0.0,0.0);
1022 
1023  // Calculate LmunuAnumu by iterating over mu and nu
1024  // In each iteration, add LmunuAnumu to sum if mu=nu, and add
1025  // LmunuAnumu + LnumuAmunu if mu != nu, since we start nu at mu
1026  double axx=0.,azz=0.,a0z=0.,a00=0.,axy=0.;
1027  for(int mu=0;mu<4;mu++){
1028  for(int nu=mu;nu<4;nu++){
1029  imaginaryPart = 0;
1030  if(mu == nu){
1031  //if mu==nu then levi-civita = 0, so imaginary part = 0
1032  Lmunu = g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu]-g[mu][nu]*kPrimek;
1033  }else{
1034  //if mu!=nu, then g[mu][nu] = 0
1035  //This same leviCivitaIndex array can be used in the else portion when
1036  //calculating Anumu
1037  leviCivitaIndexArray[0] = mu;
1038  leviCivitaIndexArray[1] = nu;
1039  for(int a=0;a<4;a++){
1040  for(int b=0;b<4;b++){
1041  leviCivitaIndexArray[2] = a;
1042  leviCivitaIndexArray[3] = b;
1043  imaginaryPart += - leviCivita(leviCivitaIndexArray)*kPrime[a]*k[b];
1044  }
1045  }
1046  //real(Lmunu) is symmetric, and imag(Lmunu) is antisymmetric
1047  //std::complex<double> num(g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu],imaginaryPart);
1048  Lmunu = g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu] + iNum*imaginaryPart;
1049  Lnumu = g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu]+g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu ]- iNum*imaginaryPart;
1050  } // End Lmunu calculation
1051 
1052  if(mu ==0 && nu == 0){
1053  Amunu = 16.0*F1V2*(2.0*rulin[0][0]*CN+2.0*q[0]*tulin[0]+q2/2.0)+
1054  2.0*q2*xiF2V2*
1055  (4.0-4.0*rulin[0][0]/M2-4.0*q[0]*tulin[0]/M2-q02*(4.0/q2+1.0/M2)) +
1056  4.0*FA2*(2.0*rulin[0][0]+2.0*q[0]*tulin[0]+(q2/2.0-2.0*M2))-
1057  (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*q02-16.0*F1V*xiF2V*(-q2+q02)*CN;
1058  a00 = real(Amunu); // TESTING CODE
1059  sum += Lmunu*Amunu;
1060  }else if(mu == 0 && nu == 3){
1061  Amunu = 16.0*F1V2*((2.0*rulin[0][3]+tulin[0]*dq)*CN+tulin[3]*q[0])+
1062  2.0*q2*xiF2V2*
1063  (-4.0*rulin[0][3]/M2-2.0*(dq*tulin[0]+q[0]*tulin[3])/M2-dq*q[0]*(4.0/q2+1.0/M2))+
1064  4.0*FA2*((2.0*rulin[0][3]+dq*tulin[0])*CL+q[0]*tulin[3])-
1065  (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq*q[0]-
1066  16.0*F1V*xiF2V*dq*q[0];
1067  a0z= real(Amunu); // TESTING CODE
1068  Anumu = Amunu;
1069  sum += Lmunu*Anumu + Lnumu*Amunu;
1070  }else if(mu == 3 && nu == 3){
1071  Amunu = 16.0*F1V2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-q2/2.0)+
1072  2.0*q2*xiF2V2*(-4.0-4.0*rulin[3][3]/M2-4.0*dq*tulin[3]/M2-dq2*(4.0/q2+1.0/M2))+
1073  4.0*FA2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-(q2/2.0-2.0*CL*M2))-
1074  (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq2-
1075  16.0*F1V*xiF2V*(q2+dq2);
1076  azz = real(Amunu); // TESTING CODE
1077  sum += Lmunu*Amunu;
1078  }else if(mu ==1 && nu == 1){
1079  Amunu = 16.0*F1V2*(2.0*rulin[1][1]-q2/2.0)+
1080  2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[1][1]/M2) +
1081  4.0*FA2*(2.0*rulin[1][1]-(q2/2.0-2.0*CT*M2))-
1082  16.0*F1V*xiF2V*CT*q2;
1083  axx = real(Amunu); // TESTING CODE
1084  sum += Lmunu*Amunu;
1085  }else if(mu == 2 && nu == 2){
1086  // Ayy not explicitly listed in paper. This is included so rotating the
1087  // coordinates of k and k' about the z-axis does not change the xsec.
1088  Amunu = 16.0*F1V2*(2.0*rulin[2][2]-q2/2.0)+
1089  2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[2][2]/M2) +
1090  4.0*FA2*(2.0*rulin[2][2]-(q2/2.0-2.0*CT*M2))-
1091  16.0*F1V*xiF2V*CT*q2;
1092  sum += Lmunu*Amunu;
1093  }else if(mu ==1 && nu == 2){
1094  Amunu = sign*16.0*iNum*FA*(xiF2V+F1V)*(-dq*tulin[0]*CT + q[0]*tulin[3]);
1095  Anumu = -Amunu; // Im(A) is antisymmetric
1096  axy = imag(Amunu); // TESTING CODE
1097  sum += Lmunu*Anumu+Lnumu*Amunu;
1098  }
1099  // All other terms will be 0 because the initial nucleus is at rest and
1100  // qTilde is in the z direction
1101 
1102  } // End loop over nu
1103  } // End loop over mu
1104 
1105  // TESTING CODE
1107  // get tmu
1108  double tmugev = leptonMom.E() - leptonMom.Mag();
1109  // Print Q2, form factors, and tensor elts
1110  std::ofstream ffstream;
1111  ffstream.open(fTensorsOutFile, std::ios_base::app);
1112  if(q0 > 0){
1113  ffstream << -q2 << "\t" << q[0] << "\t" << dq
1114  << "\t" << axx << "\t" << azz << "\t" << a0z
1115  << "\t" << a00 << "\t" << axy << "\t"
1116  << CT << "\t" << CL << "\t" << CN << "\t"
1117  << tmugev << "\t" << imU << "\t"
1118  << F1V << "\t" << xiF2V << "\t"
1119  << FA << "\t" << Fp << "\t"
1120  << tulin[0] << "\t"<< tulin[1] << "\t"
1121  << tulin[2] << "\t"<< tulin[3] << "\t"
1122  << rulin[0][0]<< "\t"<< rulin[0][1]<< "\t"
1123  << rulin[0][2]<< "\t"<< rulin[0][3]<< "\t"
1124  << rulin[1][0]<< "\t"<< rulin[1][1]<< "\t"
1125  << rulin[1][2]<< "\t"<< rulin[1][3]<< "\t"
1126  << rulin[2][0]<< "\t"<< rulin[2][1]<< "\t"
1127  << rulin[2][2]<< "\t"<< rulin[2][3]<< "\t"
1128  << rulin[3][0]<< "\t"<< rulin[3][1]<< "\t"
1129  << rulin[3][2]<< "\t"<< rulin[3][3]<< "\t"
1130  << fVc << "\t" << fCoulombFactor << "\t";
1131  ffstream << "\n";
1132  }
1133  ffstream.close();
1134  }
1135  // END TESTING CODE
1136 
1137  // Since the real parts of A and L are both symmetric and the imaginary
1138  // parts are antisymmetric, the contraction should be real
1139  if ( imag(sum) > kASmallNum )
1140  LOG("Nieves",pWARN) << "Imaginary part of tensor contraction is nonzero "
1141  << "in QEL XSec, real(sum) = " << real(sum)
1142  << "imag(sum) = " << imag(sum);
1143 
1144  return real(sum);
1145 }
1146 
1147 //___________________________________________________________________________
1148 // Auxiliary scalar function for internal integration
1149 //____________________________________________________________________________
1151  double Rcurr, int A, int Z):
1152 ROOT::Math::IBaseFunctionOneDim()
1153 {
1154  fRcurr = Rcurr;
1155  fA = A;
1156  fZ = Z;
1157 }
1158 //____________________________________________________________________________
1160 {
1161 
1162 }
1163 //____________________________________________________________________________
1165 {
1166  return 1;
1167 }
1168 //____________________________________________________________________________
1170 {
1171  double rhop = fZ*nuclear::Density(rin,fA);
1172  if(rin<fRcurr){
1173  return rhop*rin*rin/fRcurr;
1174  }else{
1175  return rhop*rin;
1176  }
1177 }
1178 //____________________________________________________________________________
1179 ROOT::Math::IBaseFunctionOneDim *
1181 {
1183 }
1184 //____________________________________________________________________________
1185 
1186 //____________________________________________________________________________
1187 //
1188 // NOTE: THE REMAINING IS TESTING CODE
1189 //
1190 // This method prints the tensor elements (as well as various inputs) for
1191 // different kinematics. The tensor elements can then be compared to the
1192 // output of Nieves' fortran code.
1193 //
1194 // The results of this code will only agree exactlly with Nieves' fortran
1195 // if Dipole form factors are set (in UserPhysicsOptions).
1196 //
1198  const {
1199  Interaction * interaction = new Interaction(*in); // copy in
1200 
1201  // Set input values here
1202  double ein = 0.2;
1203  double ctl = 0.5;
1204  double rmaxfrac = 0.25;
1205 
1206  bool carbon = false; // true -> C12, false -> Pb208
1207 
1208  if(fRPA)
1209  fTensorsOutFile = "gen.RPA";
1210  else
1211  fTensorsOutFile = "gen.noRPA";
1212 
1213  // Calculate radius
1214  bool klave;
1215  double rp,ap,rn,an;
1216  if(carbon){
1217  klave = true;
1218  rp = 1.692;
1219  ap = 1.082;
1220  rn = 1.692;
1221  an = 1.082;
1222  }else{
1223  // Pb208
1224  klave = false;
1225  rp = 6.624;
1226  ap = 0.549;
1227  rn = 6.890;
1228  an = 0.549;
1229  }
1230  double rmax;
1231  if(!klave)
1232  rmax = TMath::Max(rp,rn) + 9.25*TMath::Max(ap,an);
1233  else
1234  rmax = TMath::Sqrt(20.0)*TMath::Max(rp,rn);
1235  double r = rmax * rmaxfrac;
1236 
1237  // Relevant objects and parameters
1238  //const Kinematics & kinematics = interaction -> Kine();
1239  const InitialState & init_state = interaction -> InitState();
1240  const Target & target = init_state.Tgt();
1241 
1242  // Parameters required for LmunuAnumu
1243  double M = target.HitNucMass();
1244  double ml = interaction->FSPrimLepton()->Mass();
1245  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
1246 
1247  // Iterate over lepton energy (which then affects q, which is passed to
1248  // LmunuAnumu using in and out NucleonMom
1249  double delta = (ein-0.025)/100.0;
1250  for(int it=0;it<100;it++){
1251  double tmu = it*delta;
1252  double eout = ml + tmu;
1253  double pout = TMath::Sqrt(eout*eout-ml*ml);
1254 
1255  double pin = TMath::Sqrt(ein*ein); // Assume massless neutrinos
1256 
1257  double q0 = ein-eout;
1258  double dq = TMath::Sqrt(pin*pin+pout*pout-2.0*ctl*pin*pout);
1259  double q2 = q0*q0-dq*dq;
1260  interaction->KinePtr()->SetQ2(-q2);
1261 
1262  // When this method is called, inNucleonMomOnShell is unused.
1263  // I can thus provide the calculated values using a null vector for
1264  // inNucleonMomOnShell. I also need to put qTildeP4 in the z direction, as
1265  // Nieves does in his paper.
1266  TLorentzVector qTildeP4(0, 0, dq, q0);
1267  TLorentzVector inNucleonMomOnShell(0,0,0,0);
1268 
1269  // neutrinoMom and leptonMom only directly affect the leptonic tensor, which
1270  // we are not calculating now. Use them to transfer q.
1271  TLorentzVector neutrinoMom(0,0,pout+dq,eout+q0);
1272  TLorentzVector leptonMom(0,0,pout,eout);
1273 
1274  if(fCoulomb){ // Use same steps as in XSec()
1275  // Coulomb potential
1276  double Vc = vcr(& target, r);
1277  fVc = Vc;
1278 
1279  // Outgoing lepton energy and momentum including coulomb potential
1280  int sign = is_neutrino ? 1 : -1;
1281  double El = leptonMom.E();
1282  double ElLocal = El - sign*Vc;
1283  if(ElLocal - ml <= 0.0){
1284  LOG("Nieves",pINFO) << "Event should be rejected. Coulomb effects "
1285  << "push kinematics below threshold";
1286  return;
1287  }
1288  double plLocal = TMath::Sqrt(ElLocal*ElLocal-ml*ml);
1289 
1290  // Correction factor
1291  double coulombFactor= plLocal*ElLocal/leptonMom.Vect().Mag()/El;
1292  fCoulombFactor = coulombFactor; // Store and print
1293  }
1294 
1295  // TODO: apply Coulomb correction to 3-momentum transfer dq
1296 
1297  fFormFactors.Calculate(interaction);
1298  LmunuAnumu(neutrinoMom, inNucleonMomOnShell, leptonMom, qTildeP4,
1299  M, is_neutrino, target, false);
1300  }
1301  return;
1302 } // END TESTING CODE
1303 //____________________________________________________________________________
Cross Section Calculation Interface.
const double kPi
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Basic constants.
Double_t angle
Definition: plot.C:86
std::complex< double > ruLinRelX(double q0, double qm, double kf, double m) const
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:31
bool IsWeakCC(void) const
bool fRPA
use RPA corrections
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
const XML_Char * target
Definition: expat.h:268
void Configure(const Registry &config)
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int HitNucPdg(void) const
Definition: Target.cxx:321
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
double delta
Definition: runWimpSim.h:98
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
double Density(double r, int A, double ring=0.)
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
const FermiMomentumTable * fKFTable
int RecoilNucleonPdg(void) const
recoil nucleon pdg
static const double fermi
Definition: Units.h:63
int A(void) const
Definition: Target.h:71
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
double HitNucMass(void) const
Definition: Target.cxx:250
T sqrt(T number)
Definition: d0nt_math.hpp:156
#define pFATAL
Definition: Messenger.h:57
bool IsNucleus(void) const
Definition: Target.cxx:289
constexpr T pow(T x)
Definition: pow.h:75
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, bool assumeFreeNucleon) const
int Pdg(void) const
Definition: Target.h:72
static FermiMomentumTablePool * Instance(void)
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
Examines whether the generated event should be Pauli blocked. Is a concerete implementation of the Ev...
Definition: PauliBlocker.h:36
const XSecIntegratorI * fXSecIntegrator
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:67
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
void abs(TH1 *hist)
Definition: config.py:1
QELFormFactors fFormFactors
bool fCoulomb
use Coulomb corrections
enum genie::EKinePhaseSpace KinePhaseSpace_t
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
float abs(float number)
Definition: d0nt_math.hpp:39
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
std::complex< double > deltaLindhard(double q0gev, double dqgev, double rho, double kFgev) const
virtual NuclearModel_t ModelType(const Target &) const =0
const QELFormFactorsModelI * fFormFactorsModel
static const double kAem
Definition: Constants.h:57
double fhbarc
hbar*c in GeV*fm
Float_t Z
Definition: plot.C:38
Double_t q2[12][num]
Definition: f2_nu.C:137
double fCos8c2
cos^2(cabibbo angle)
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
const XML_Char * s
Definition: expat.h:262
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:92
double GetFermiMomentum(const Target &tgt, int pdg_Nf, double radius=0.0) const
Get the Fermi momentum needed to check Pauli blocking.
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
const NuclearModelI * fNuclModel
Nuclear Model for integration.
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:66
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
string Name(void) const
Definition: AlgId.h:45
Float_t E
Definition: plot.C:20
const FermiMomentumTable * GetTable(string name)
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
Definition: QELUtils.cxx:51
static string AsString(KinePhaseSpace_t kps)
bool fCompareNievesTensors
print tensors
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
const double a
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
double Integral(const Interaction *i) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
bool fDoPauliBlocking
Whether to apply Pauli blocking in XSec()
double func(double x, double y)
int Z(void) const
Definition: Target.h:69
const double j
Definition: BetheBloch.cxx:29
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:63
double fXSecScale
external xsec scaling factor
TVector3 Unit() const
Misc GENIE control constants.
z
Definition: test.py:28
double xiF2V(void) const
Get the computed form factor xi*F2V.
A very simple service to remember what detector we&#39;re working in.
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
static const double kLightSpeed
Definition: Constants.h:32
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
int kf2[14]
Definition: runWimpSim.h:94
double vcr(const Target *target, double r) const
static constexpr Double_t m2
Definition: Munits.h:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TString fTensorsOutFile
file to print tensors to
static const double A
Definition: Units.h:82
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
ifstream in
Definition: comparison.C:7
int N(void) const
Definition: Target.h:70
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
double HitNucPosition(void) const
Definition: Target.h:90
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:195
app
Definition: demo.py:189
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
exit(0)
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
int leviCivita(int input[]) const
void CompareNievesTensors(const Interaction *i) const
const int kPdgProton
Definition: PDGCodes.h:65
double F1V(void) const
Get the computed form factor F1V.
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
std::complex< double > relLindhard(double q0gev, double dqgev, double kFgev, double M, bool isNeutrino, std::complex< double > relLindIm) const
static const double kGF2
Definition: Constants.h:60
Double_t sum
Definition: plot.C:31
double Fp(void) const
Get the computed form factor Fp.
double LmunuAnumu(const TLorentzVector neutrinoMom, const TLorentzVector inNucleonMom, const TLorentzVector leptonMom, const TLorentzVector outNucleonMom, double M, bool is_neutrino, const Target &target, bool assumeFreeNucleon) const
void kinematics()
Definition: kinematics.C:10
bool gAbortingInErr
Definition: Messenger.cxx:56
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
static const double kPlankConstant
Definition: Constants.h:33
const int kPdgNeutron
Definition: PDGCodes.h:67
double FA(void) const
Get the computed form factor FA.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
double Density(double r)
Definition: PREM.cxx:26
static const double kPi2
Definition: Constants.h:39
Root of GENIE utility namespaces.
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
def sign(x)
Definition: canMan.py:197
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
std::complex< double > relLindhardIm(double q0gev, double dqgev, double kFngev, double kFpgev, double M, bool isNeutrino, double &t0, double &r00) const
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
static const double kPionMass2
Definition: Constants.h:87
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353