GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NievesSimoVacasMECPXSec2016.cxx
Go to the documentation of this file.
1 //_________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Original code contributed by J. Schwehr, D. Cherdack, R. Gran
7  Substantial code refactorizations by the core GENIE group.
8 */
9 //_________________________________________________________________________
10 
24 
25 using namespace genie;
26 using namespace genie::constants;
27 
28 //_________________________________________________________________________
30 XSecAlgorithmI("genie::NievesSimoVacasMECPXSec2016")
31 {
32 
33 }
34 //_________________________________________________________________________
36 XSecAlgorithmI("genie::NievesSimoVacasMECPXSec2016", config)
37 {
38 
39 }
40 //_________________________________________________________________________
42 {
43 
44 }
45 //_________________________________________________________________________
47  const Interaction * interaction, KinePhaseSpace_t kps) const
48 {
49  // If {W,Q2} have been supplied instead, compute {Tl, ctl}
50  // NOTE: The expressions used here neglect Fermi motion and
51  // should eventually be revisited. See the "important note"
52  // in src/Framework/Utils/KineUtils.cxx about the
53  // Jacobian for transforming {W,Q2} --> {Tl, ctl}.
54  // - S. Gardiner, 29 July 2020
55  if ( kps == kPSWQ2fE ) {
56 
57  double Q2 = interaction->Kine().GetKV( kKVQ2 );
58  double W = interaction->Kine().GetKV( kKVW );
59 
60  // Probe properties (mass, energy, momentum)
61  const InitialState& init_state = interaction->InitState();
62  double mv = init_state.Probe()->Mass();
63  double Ev = init_state.ProbeE( kRfLab );
64  double pv = std::sqrt( std::max(0., Ev*Ev - mv*mv) );
65 
66  // Invariant mass of the initial hit nucleon
67  const TLorentzVector& hit_nuc_P4 = init_state.Tgt().HitNucP4();
68  double M = hit_nuc_P4.M();
69 
70  // Get the outgoing lepton kinetic energy
71  double ml = interaction->FSPrimLepton()->Mass();
72  double Tl = Ev - ml - ( (W*W + Q2 - M*M) / (2.*M) );
73 
74  // Get the outgoing lepton scattering cosine
75  double El = Tl + ml;
76  double pl = std::sqrt( std::max(0., El*El - ml*ml) );
77  double ctl = ( 2.*Ev*El - Q2 - mv*mv - ml*ml ) / ( 2. * pv * pl );
78 
79  // Set Tl, ctl in the interaction
80  interaction->KinePtr()->SetKV( kKVTl, Tl );
81  interaction->KinePtr()->SetKV( kKVctl, ctl );
82  }
83 
84  // This function returns d2sigma/(dTmu dcos_mu) in GeV^(-3)
85  int target_pdg = interaction->InitState().Tgt().Pdg();
86 
87  int A_request = pdg::IonPdgCodeToA(target_pdg);
88  int Z_request = pdg::IonPdgCodeToZ(target_pdg);
89 
90  // To generate cross-sections for nuclei other than those with hadron
91  // tensors we need to pull both the full cross-section and
92  // the pn initial state fraction.
93  // Non-isoscalar nuclei are beyond the original published Valencia model
94  // and scale with A according to the number of pp, pn, or nn pairs
95  // the probe is expected to find.
96  // There is some by-hand optimization here, skipping the delta part when
97  // only the total cross-section is requested.
98  // Possible future models without a Delta had tensor would also use that
99  // flag to call this without computing the Delta part.
100 
101  // Try to look up a hadron tensor in the pool that is an exact match for
102  // the target nucleus. If an exact match cannot be found, decide upon a
103  // suitable substitute based on the mass number A and proton number Z.
104 
105  int tensor_pdg = target_pdg;
106 
107  /// \todo Replace these hard-coded replacements with an equivalent XML
108  /// configuration
110  {
111 
112  if ( A_request == 4 && Z_request == 2 ) {
113  tensor_pdg = kPdgTgtC12;
114  // This is for helium 4, but use carbon tensor
115  // the use of nuclear density parameterization is suspicious
116  // but some users (MINERvA) need something not nothing.
117  // The pn will be exactly 1/3, but pp and nn will be ~1/4
118  // Because the combinatorics are different.
119  // Could do lithium beryllium boron which you don't need
120  }
121  else if (A_request < 9) {
122  // refuse to do D, T, He3, Li, and some Be, B
123  // actually it would work technically, maybe except D, T
124  MAXLOG("NievesSimoVacasMEC", pWARN, 10)
125  << "Asked to scale to deuterium through boron "
126  << target_pdg << " nope, lets not do that.";
127  return 0;
128  }
129  else if (A_request >= 9 && A_request < 15) {
130  tensor_pdg = kPdgTgtC12;
131  //}
132  // could explicitly put in nitrogen for air
133  //else if ( A_request >= 14 && A < 15) { // AND CHANGE <=14 to <14.
134  // tensor_pdg = kPdgTgtN14;
135  }
136  else if (A_request >= 15 && A_request < 22) {
137  tensor_pdg = kPdgTgtO16;
138  }
139  else if (A_request >= 22 && A_request < 33) {
140  // of special interest, this gets Al27 and Si28
141  tensor_pdg = 1000140280;
142  }
143  else if (A_request >= 33 && A_request < 50) {
144  // of special interest, this gets Ar40 and Ti48
145  tensor_pdg = kPdgTgtCa40;
146  }
147  else if (A_request >= 50 && A_request < 90) {
148  // pseudoFe56, also covers many other ferrometals and Ge
149  tensor_pdg = 1000280560;
150  }
151  else if (A_request >= 90 && A_request < 160) {
152  // use Ba112 = PseudoCd. Row5 of Periodic table useless. Ag, Xe?
153  tensor_pdg = 1000561120;
154  }
155  else if (A_request >= 160) {
156  // use Rf208 = pseudoPb
157  tensor_pdg = 1001042080;
158  }
159  else {
160  MAXLOG("NievesSimoVacasMEC", pWARN, 10)
161  << "Asked to scale to a nucleus "
162  << target_pdg << " which we don't know yet.";
163  return 0;
164  }
165  }
166 
167  // Check that the input kinematical point is within the range
168  // in which hadron tensors are known (for chosen target)
169  double Ev = interaction->InitState().ProbeE(kRfLab);
170  double Tl = interaction->Kine().GetKV(kKVTl);
171  double costl = interaction->Kine().GetKV(kKVctl);
172  double ml = interaction->FSPrimLepton()->Mass();
173  double Q0 = interaction->Kine().GetKV(kKVQ0);
174  double Q3 = interaction->Kine().GetKV(kKVQ3);
175 
176  const LabFrameHadronTensorI* tensor
177  = dynamic_cast<const LabFrameHadronTensorI*>(
179 
180  // If retrieving the tensor failed, complain and return zero
181  if ( !tensor ) {
182  LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a"
183  " hadronic tensor for the nuclide " << tensor_pdg;
184  return 0.;
185  }
186 
187  // Assume for now that the range of validity for the "FullAll" hadron
188  // tensor is the same as for the partial hadron tensors
189  /// \todo Revisit this assumption, and perhaps implement something
190  /// more robust
191  double Q0min = tensor->q0Min();
192  double Q0max = tensor->q0Max();
193  double Q3min = tensor->qMagMin();
194  double Q3max = tensor->qMagMax();
195  if (Q0 < Q0min || Q0 > Q0max || Q3 < Q3min || Q3 > Q3max) {
196  return 0.0;
197  }
198 
199  // Get the Q-value needed to calculate the cross sections using the
200  // hadron tensor.
201  /// \todo Shouldn't we get this from the nuclear model?
202  int nu_pdg = interaction->InitState().ProbePdg();
203  double Q_value = genie::utils::mec::Qvalue(target_pdg, nu_pdg);
204 
205  // Apply Qvalue relative shift if needed:
206  if( fQvalueShifter ) Q_value += Q_value * fQvalueShifter -> Shift( interaction->InitState().Tgt() ) ;
207 
208  // By default, we will compute the full cross-section. If a resonance is
209  // set, we will calculate the part of the cross-section with an internal
210  // Delta line without a final state pion (usually called PPD for pioness
211  // Delta decay). If a {p,n} hit dinucleon was set we will calculate the
212  // cross-section for that component only (either full or PDD cross-section)
213  bool delta = interaction->ExclTag().KnownResonance();
214  bool pn = (interaction->InitState().Tgt().HitNucPdg() == kPdgClusterNP);
215 
216  double xsec_all = 0.;
217  double xsec_pn = 0.;
218  if ( delta ) {
219 
220  const LabFrameHadronTensorI* tensor_delta_all
221  = dynamic_cast<const LabFrameHadronTensorI*>(
223 
224  if ( !tensor_delta_all ) {
225  LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"DeltaAll\""
226  << " hadronic tensor for nuclide " << tensor_pdg;
227  return 0.;
228  }
229 
230  const LabFrameHadronTensorI* tensor_delta_pn
231  = dynamic_cast<const LabFrameHadronTensorI*>(
233 
234  if ( !tensor_delta_pn ) {
235  LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"Deltapn\""
236  << " hadronic tensor for nuclide " << tensor_pdg;
237  return 0.;
238  }
239 
240  xsec_all = tensor_delta_all->dSigma_dT_dCosTheta(interaction, Q_value);
241  xsec_pn = tensor_delta_pn->dSigma_dT_dCosTheta(interaction, Q_value);
242 
243  }
244  else {
245 
246  const LabFrameHadronTensorI* tensor_full_all
247  = dynamic_cast<const LabFrameHadronTensorI*>(
249 
250  if ( !tensor_full_all ) {
251  LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"FullAll\""
252  << " hadronic tensor for nuclide " << tensor_pdg;
253  return 0.;
254  }
255 
256  const LabFrameHadronTensorI* tensor_full_pn
257  = dynamic_cast<const LabFrameHadronTensorI*>(
259 
260  if ( !tensor_full_pn ) {
261  LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"Fullpn\""
262  << " hadronic tensor for nuclide " << tensor_pdg;
263  return 0.;
264  }
265 
266  xsec_all = tensor_full_all->dSigma_dT_dCosTheta(interaction, Q_value);
267  xsec_pn = tensor_full_pn->dSigma_dT_dCosTheta(interaction, Q_value);
268  }
269 
270  // We need to scale the cross section appropriately if
271  // we are using a hadronic tensor for a nuclide that is different
272  // from the actual target
273  bool need_to_scale = (target_pdg != tensor_pdg);
274 
275  // would need to trap and treat He3, T, D special here.
276  if ( need_to_scale ) {
277 
278  double PP = Z_request;
279  double NN = A_request - PP;
280  double P = pdg::IonPdgCodeToZ(tensor_pdg);
281  double N = pdg::IonPdgCodeToA(tensor_pdg) - P;
282 
283  double scale_pn = TMath::Sqrt( (PP*NN)/(P*N) );
284  double scale_pp = TMath::Sqrt( (PP * (PP - 1.)) / (P * (P - 1.)) );
285  double scale_nn = TMath::Sqrt( (NN * (NN - 1.)) / (N * (N - 1.)) );
286 
287  LOG("NievesSimoVacasMEC", pDEBUG)
288  << "Scale pn pp nn for (" << target_pdg << ", " << tensor_pdg << ")"
289  << " : " << scale_pn << " " << scale_pp << " " << scale_nn;
290 
291  // This is an approximation in at least three senses:
292  // 1. We are scaling from an isoscalar nucleus using p and n counting
293  // 2. We are not using the right qvalue in the had tensor
294  // 3. We are not scaling the Delta faster than the non-Delta.
295  // The guess is that these are good approximations.
296  // A test we could document is to scale from O16 to N14 or C12 using this
297  // algorithm and see how many percent deviation we see from the full
298  // calculation.
299  double temp_all = xsec_all;
300  double temp_pn = xsec_pn * scale_pn;
301  if (nu_pdg > 0) {
302  // matter neutrinos
303  temp_all = xsec_pn * scale_pn + (xsec_all - xsec_pn) * scale_nn;
304  }
305  else {
306  // antineutrinos
307  temp_all = xsec_pn * scale_pn + (xsec_all - xsec_pn) * scale_pp;
308  }
309 
310  xsec_all = temp_all;
311  xsec_pn = temp_pn;
312 
313  }
314 
315  // Choose the right kind of cross section ("all" or "pn") to return
316  // based on whether a {p, n} dinucleon was hit
317  double xsec = (pn) ? xsec_pn : xsec_all;
318 
319  // Apply given scaling factor
320  const ProcessInfo& proc_info = interaction->ProcInfo();
321  if( proc_info.IsWeakCC() ) xsec *= fXSecCCScale;
322  else if( proc_info.IsWeakNC() ) xsec *= fXSecNCScale;
323 
324  if( fMECScaleAlg ) xsec *= fMECScaleAlg->GetScaling( * interaction ) ;
325 
326  if ( kps != kPSTlctl && kps != kPSWQ2fE ) {
327  LOG("NievesSimoVacasMEC", pWARN)
328  << "Doesn't support transformation from "
329  << KinePhaseSpace::AsString(kPSTlctl) << " to "
330  << KinePhaseSpace::AsString(kps);
331  xsec = 0;
332  }
333  else if ( kps == kPSWQ2fE && xsec != 0. ) {
334  double J = utils::kinematics::Jacobian( interaction, kPSTlctl, kps );
335  xsec *= J;
336  }
337 
338  return xsec;
339 }
340 //_________________________________________________________________________
342  const Interaction * interaction) const
343 {
344  double xsec = fXSecIntegrator->Integrate(this,interaction);
345  return xsec;
346 }
347 //_________________________________________________________________________
349  const Interaction * interaction) const
350 {
351  if (interaction->TestBit(kISkipProcessChk)) return true;
352 
353  const ProcessInfo & proc_info = interaction->ProcInfo();
354  if (!proc_info.IsMEC()) {
355  return false;
356  }
357  return true;
358 }
359 //_________________________________________________________________________
361 {
362  Algorithm::Configure(config);
363  this->LoadConfig();
364 }
365 //____________________________________________________________________________
367 {
368  Algorithm::Configure(config);
369  this->LoadConfig();
370 }
371 //_________________________________________________________________________
373 {
374  bool good_config = true;
375 
376  // Cross section scaling factor
377  GetParam( "MEC-CC-XSecScale", fXSecCCScale ) ;
378  GetParam( "MEC-NC-XSecScale", fXSecNCScale ) ;
379 
380  fHadronTensorModel = dynamic_cast<const HadronTensorModelI *> ( this->SubAlg("HadronTensorAlg") );
381  if( !fHadronTensorModel ) {
382  good_config = false ;
383  LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required HadronTensorAlg does not exist. AlgID is : " << SubAlg("HadronTensorAlg")->Id() ;
384  }
385 
386  fXSecIntegrator = dynamic_cast<const XSecIntegratorI *> (this->SubAlg("NumericalIntegrationAlg"));
387  if( !fXSecIntegrator ) {
388  good_config = false ;
389  LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required NumericalIntegrationAlg does not exist. AlgID is : " << SubAlg("NumericalIntegrationAlg")->Id();
390  }
391 
392  // Read optional QvalueShifter:
393  fQvalueShifter = nullptr;
394  if( GetConfig().Exists("QvalueShifterAlg") ) {
395  fQvalueShifter = dynamic_cast<const QvalueShifter *> ( this->SubAlg("QvalueShifterAlg") );
396  if( !fQvalueShifter ) {
397  good_config = false ;
398  LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required QvalueShifterAlg does not exist. AlgID is : " << SubAlg("QvalueShifterAlg")->Id() ;
399  }
400  }
401 
402  // Read optional MECScaleVsW:
403  fMECScaleAlg = nullptr;
404  if( GetConfig().Exists("MECScaleAlg") ) {
405  fMECScaleAlg = dynamic_cast<const XSecScaleI *> ( this->SubAlg("MECScaleAlg") );
406  if( !fMECScaleAlg ) {
407  good_config = false ;
408  LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required MECScaleAlg cannot be casted. AlgID is : " << SubAlg("MECScaleAlg")->Id() ;
409  }
410  }
411 
412  if( ! good_config ) {
413  LOG("NievesSimoVacasMECPXSec2016", pERROR) << "Configuration has failed.";
414  exit(78) ;
415  }
416 
417 }
Cross Section Calculation Interface.
virtual double dSigma_dT_dCosTheta(const Interaction *interaction, double Q_value) const =0
virtual const HadronTensorI * GetTensor(int tensor_pdg, HadronTensorType_t type) const =0
bool IsWeakCC(void) const
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
#define pERROR
Definition: Messenger.h:59
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int HitNucPdg(void) const
Definition: Target.cxx:304
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
const int kPdgClusterNP
Definition: PDGCodes.h:215
bool KnownResonance(void) const
Definition: XclsTag.h:68
This class is responsible to compute a scaling factor for the XSec.
Definition: XSecScaleI.h:25
int Pdg(void) const
Definition: Target.h:71
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
TParticlePDG * Probe(void) const
Abstract interface for an object that computes the elements ( , , etc.) and structure functions ( ...
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double fXSecCCScale
external xsec scaling factor
enum genie::EKinePhaseSpace KinePhaseSpace_t
virtual double q0Max() const =0
#define MAXLOG(s, l, c)
Similar to LOG(stream,priority) but quits after &quot;maxcount&quot; messages.
Definition: Messenger.h:241
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
double fXSecNCScale
external xsec scaling factor
const int kPdgTgtCa40
Definition: PDGCodes.h:204
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
const int kPdgTgtO16
Definition: PDGCodes.h:203
double Qvalue(int targetpdg, int nupdg)
Definition: MECUtils.cxx:164
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
bool IsWeakNC(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static string AsString(KinePhaseSpace_t kps)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
#define pWARN
Definition: Messenger.h:60
virtual double qMagMax() const =0
virtual double q0Min() const =0
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsMEC(void) const
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
double Integral(const Interaction *i) const
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
virtual double GetScaling(const Interaction &) const =0
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
const int kPdgTgtC12
Definition: PDGCodes.h:202
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
Creates hadron tensor objects for use in cross section calculations.
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:55
virtual double qMagMin() const =0
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345