NievesSimoVacasMECPXSec2016.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  For the class documentation see the corresponding header file.
8 */
9 //_________________________________________________________________________
10 
23 
24 using namespace genie;
25 using namespace genie::constants;
26 
27 //_________________________________________________________________________
29 XSecAlgorithmI("genie::NievesSimoVacasMECPXSec2016")
30 {
31 
32 }
33 //_________________________________________________________________________
35 XSecAlgorithmI("genie::NievesSimoVacasMECPXSec2016", config)
36 {
37 
38 }
39 //_________________________________________________________________________
41 {
42 
43 }
44 //_________________________________________________________________________
46  const Interaction * interaction, KinePhaseSpace_t kps) const
47 {
48 // The basic quantity calculated is d2sigma/(dTmu dcos_mu)
49 
50  // Get hadron tensor
52 
53  int targetpdg = interaction->InitState().Tgt().Pdg();
54  int Arequest = pdg::IonPdgCodeToA(targetpdg);
55  int Zrequest = pdg::IonPdgCodeToZ(targetpdg);
56 
57  // To generate cross-sections for nuclei other than those with hadron
58  // tensors we need to pull both the full cross-section and
59  // the pn initial state fraction.
60  // Non-isoscalar nuclei are beyond the original published Valencia model
61  // and scale with A according to the number of pp, pn, or nn pairs
62  // the probe is expected to find.
63  // There is some by-hand optimization here, skipping the delta part when
64  // only the total cross-section is requested.
65  // Possible future models without a Delta had tensor would also use that
66  // flag to call this without computing the Delta part.
67 
68  int tensorpdg = targetpdg;
69 
70  if( ! hadtensor->KnownTensor(tensorpdg) ) {
71 
72  // Decide which hadron tensor to use for different ranges of A.
73 
74  if( Arequest == 4 && Zrequest == 2 ){
75  tensorpdg = kPdgTgtC12;
76  // This is for helium 4, but use carbon tensor
77  // the use of nuclear density parameterization is suspicious
78  // but some users (MINERvA) need something not nothing.
79  // The pn will be exactly 1/3, but pp and nn will be ~1/4
80  // Because the combinatorics are different.
81  // Could do lithium beryllium boron which you don't need
82  }
83  else if (Arequest < 9){
84  // refuse to do D, T, He3, Li, and some Be, B
85  // actually it would work technically, maybe except D, T
86  MAXLOG("NievesSimoVacasMEC", pWARN, 10)
87  << "Asked to scale to deuterium through boron "
88  << targetpdg << " nope, lets not do that.";
89  return 0;
90  }
91  else if( Arequest >= 9 && Arequest < 15){
92  tensorpdg = kPdgTgtC12;
93  //}
94  // could explicitly put in nitrogen for air
95  //else if ( Arequest >= 14 && A < 15) { // AND CHANGE <=14 to <14.
96  // tensorpdg = kPdgTgtN14;
97  }
98  else if( Arequest >= 15 && Arequest < 22){
99  tensorpdg = kPdgTgtO16;
100  }
101  else if( Arequest >= 22 && Arequest < 33){
102  // of special interest, this gets Al27 and Si28
103  tensorpdg = 1000140280;
104  }
105  else if(Arequest >= 33 && Arequest < 50){
106  // of special interest, this gets Ar40 and Ti48
107  tensorpdg = kPdgTgtCa40;
108  }
109  else if( Arequest >= 50 && Arequest < 90){
110  // pseudoFe56, also covers many other ferrometals and Ge
111  tensorpdg = 1000280560;
112  }
113  else if( Arequest >= 90 && Arequest < 160){
114  // use Ba112 = PseudoCd. Row5 of Periodic table useless. Ag, Xe?
115  tensorpdg = 1000561120;
116  }
117  else if( Arequest >= 160 ){
118  // use Rf208 = pseudoPb
119  tensorpdg = 1001042080;
120  }
121  else {
122  MAXLOG("NievesSimoVacasMEC", pWARN, 10)
123  << "Asked to scale to a nucleus "
124  << targetpdg << " which we don't know yet.";
125  return 0;
126  }
127  }
128 
129  // Check that the input kinematical point is within the range
130  // in which hadron tensors are known (for chosen target)
131  double Ev = interaction->InitState().ProbeE(kRfLab);
132  double Tl = interaction->Kine().GetKV(kKVTl);
133  double costl = interaction->Kine().GetKV(kKVctl);
134  double ml = interaction->FSPrimLepton()->Mass();
135  double Q0 = 0;
136  double Q3 = 0;
137  genie::utils::mec::Getq0q3FromTlCostl(Tl, costl, Ev, ml, Q0, Q3);
138  const vector <genie::BLI2DNonUnifGrid *> &
139  tensor_table = hadtensor->TensorTable(
141  double Q0min = tensor_table[0]->XMin();
142  double Q0max = tensor_table[0]->XMax();
143  double Q3min = tensor_table[0]->YMin();
144  double Q3max = tensor_table[0]->YMax();
145  if(Q0 < Q0min || Q0 > Q0max || Q3 < Q3min || Q3 > Q3max) {
146  return 0.0;
147  }
148 
149  // By default, compute will compute the full cross-section.
150  // If a resonance is set, will will calculate the part of the cross-section
151  // with an internal Delta line without a final state pion (usually called PPD
152  // for pioness Delta decay).
153  // If a {p,n} hit dinucleon was set will calculate the cross-section
154  // for that component only (either full or PDD cross-section)
155  bool delta = interaction->ExclTag().KnownResonance();
156  bool pn = (interaction->InitState().Tgt().HitNucPdg() == kPdgClusterNP);
157  //LOG("NievesSimoVacasMEC", pDEBUG) << "delta: " << delta << ", pn: " << pn;
158 
159  double xsec_all = 0;
160  double xsec_pn = 0;
161  if(delta) {
162  xsec_all = genie::utils::mec::TensorContraction(interaction,
164  xsec_pn = genie::utils::mec::TensorContraction(interaction,
166  } else {
167  xsec_all = genie::utils::mec::TensorContraction(interaction,
169  xsec_pn = genie::utils::mec::TensorContraction(interaction,
171  }
172 
173  // now again, run this only if targetpdg != tensorpdg.
174  // would need to trap and treat He3, T, D special here.
175  if( ! hadtensor->KnownTensor(targetpdg) ) {
176  // if we need to scale, figure it out here.
177 
178  double PP = Zrequest;
179  double NN = Arequest - PP;
180  double P = pdg::IonPdgCodeToZ(tensorpdg);
181  double N = pdg::IonPdgCodeToA(tensorpdg) - P;
182 
183  double scale_pn = TMath::Sqrt( (PP*NN)/(P*N) );
184  double scale_pp = TMath::Sqrt( (PP * (PP - 1.)) / (P * (P - 1.)) );
185  double scale_nn = TMath::Sqrt( (NN * (NN - 1.)) / (N * (N - 1.)) );
186 
187  LOG("NievesSimoVacasMEC", pDEBUG)
188  << "Scale pn pp nn for (" << targetpdg << ", " << tensorpdg << ")"
189  << " : " << scale_pn << " " << scale_pp << " " << scale_nn;
190 
191  // this is an approximation in at least three senses.
192  // we are scaling from an isoscalar nucleus using p and n counting
193  // we are not using the right qvalue in the had tensor
194  // we are not scaling the Delta faster than the non-Delta.
195  // the guess is that these are good approximations.
196  // a test we could document is to scale from O16 to N14 or C12 using this algorithm
197  // and see how many percent deviation we see from the full calculation.
198 
199  double temp_all = xsec_all;
200  double temp_pn = xsec_pn * scale_pn;
201  int nupdg = interaction->InitState().ProbePdg();
202  if(nupdg > 0){
203  temp_all = xsec_pn * scale_pn + (xsec_all - xsec_pn) * scale_nn;
204  } else {
205  temp_all = xsec_pn * scale_pn + (xsec_all - xsec_pn) * scale_pp;
206  }
207  xsec_all = temp_all;
208  xsec_pn = temp_pn;
209  }
210 
211  double xsec = (pn) ? xsec_pn : xsec_all;
212 
213  // Apply given scaling factor
214  xsec *= fXSecScale;
215 
216  if(kps!=kPSTlctl) {
217  LOG("NievesSimoVacasMEC", pWARN)
218  << "Doesn't support transformation from "
219  << KinePhaseSpace::AsString(kPSTlctl) << " to "
220  << KinePhaseSpace::AsString(kps);
221  xsec = 0;
222  }
223 
224  return xsec;
225 }
226 //_________________________________________________________________________
228  const Interaction * interaction) const
229 {
230  double xsec = fXSecIntegrator->Integrate(this,interaction);
231  return xsec;
232 }
233 //_________________________________________________________________________
235  const Interaction * interaction) const
236 {
237  if (interaction->TestBit(kISkipProcessChk)) return true;
238 
239  const ProcessInfo & proc_info = interaction->ProcInfo();
240  if (!proc_info.IsMEC()) {
241  return false;
242  }
243  return true;
244 }
245 //_________________________________________________________________________
247 {
248  Algorithm::Configure(config);
249  this->LoadConfig();
250 }
251 //____________________________________________________________________________
253 {
254  Algorithm::Configure(config);
255  this->LoadConfig();
256 }
257 //_________________________________________________________________________
259 {
260  // Cross section scaling factor
261  GetParam( "MEC-CC-XSecScale", fXSecScale ) ;
262 
264  dynamic_cast<const XSecIntegratorI *> (
265  this->SubAlg("NumericalIntegrationAlg"));
267 
268 }
269 //_________________________________________________________________________
Cross Section Calculation Interface.
Basic constants.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
int HitNucPdg(void) const
Definition: Target.cxx:321
double delta
Definition: runWimpSim.h:98
const int kPdgClusterNP
Definition: PDGCodes.h:192
bool KnownResonance(void) const
Definition: XclsTag.h:61
static MECHadronTensor * Instance()
int Pdg(void) const
Definition: Target.h:72
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
#define MAXLOG(s, l, c)
Similar to LOG(stream,priority) but quits after "maxcount" messages.
Definition: Messenger.h:244
double fXSecScale
external xsec scaling factor
const int kPdgTgtCa40
Definition: PDGCodes.h:181
const int kPdgTgtO16
Definition: PDGCodes.h:180
#define P(a, b, c, d, e, x)
Summary information for an interaction.
Definition: Interaction.h:56
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...
Definition: Messenger.h:97
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:44
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
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)
Definition: MECUtils.cxx:115
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:333
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsMEC(void) const
double Integral(const Interaction *i) const
bool KnownTensor(int targetpdg)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
assert(nhit_max >=nhit_nbins)
const int kPdgTgtC12
Definition: PDGCodes.h:179
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:53
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
double TensorContraction(const Interaction *interaction, int tensor_pdg, MECHadronTensor::MECHadronTensorType_t tensor_type)
Definition: MECUtils.cxx:202
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define NN
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
#define pDEBUG
Definition: Messenger.h:64
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353