GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
QPMDISStrucFuncBase.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 
9  This GENIE code was adapted from the neugen3 code co-authored by Donna Naples
10  (Pittsburgh U.), Hugh Gallagher (Tufts U), and Costas Andreopoulos (RAL)
11 
12  A fix was installed (Aug 12, 2014) by Brian Tice (Rochester) so that
13  the nuclear modification to the pdf should be calculated in terms
14  of the experimental x, not the rescaled x. The same goes for R(x,Q2).
15 
16  A fix of the scaling variable used for the relations between structure
17  functions was installed by C. Bronner and J. Morrison Jun 06, 2016
18  after it was confirmed by A. Bodek that x and not the modified
19  scaling variable should be used there.
20 
21  Changes required to implement the GENIE Boosted Dark Matter module
22  were installed by Josh Berger (Univ. of Wisconsin)
23 */
24 //____________________________________________________________________________
25 
26 #include <TMath.h>
27 
29 #include "Framework/Conventions/GBuild.h"
39 
40 using namespace genie;
41 using namespace genie::constants;
42 
43 //____________________________________________________________________________
46 {
47  this->InitPDF();
48 }
49 //____________________________________________________________________________
52 {
53  this->InitPDF();
54 }
55 //____________________________________________________________________________
56 QPMDISStrucFuncBase::QPMDISStrucFuncBase(string name, string config):
57 DISStructureFuncModelI(name, config)
58 {
59  this->InitPDF();
60 }
61 //____________________________________________________________________________
63 {
64  delete fPDF;
65  delete fPDFc;
66 }
67 //____________________________________________________________________________
69 {
70  Algorithm::Configure(config);
71  this->LoadConfig();
72 }
73 //____________________________________________________________________________
74 void QPMDISStrucFuncBase::Configure(string param_set)
75 {
76  Algorithm::Configure(param_set);
77  this->LoadConfig();
78 }
79 //____________________________________________________________________________
81 {
82  LOG("DISSF", pDEBUG) << "Loading configuration...";
83 
84  //-- pdf
85  const PDFModelI * pdf_model =
86  dynamic_cast<const PDFModelI *> (this->SubAlg("PDF-Set"));
87  fPDF -> SetModel(pdf_model);
88  fPDFc -> SetModel(pdf_model);
89 
90  //-- get CKM elements
91  GetParam( "CKM-Vcd", fVcd ) ;
92  GetParam( "CKM-Vcs", fVcs ) ;
93  GetParam( "CKM-Vud", fVud ) ;
94  GetParam( "CKM-Vus", fVus ) ;
95 
96  fVcd2 = TMath::Power( fVcd, 2 );
97  fVcs2 = TMath::Power( fVcs, 2 );
98  fVud2 = TMath::Power( fVud, 2 );
99  fVus2 = TMath::Power( fVus, 2 );
100 
101  //-- charm mass
102  GetParam( "Charm-Mass", fMc ) ;
103 
104  //-- min Q2 for PDF evaluation
105  GetParam( "PDF-Q2min", fQ2min ) ;
106 
107  //-- include R (~FL)?
108  GetParam( "IncludeR", fIncludeR ) ;
109 
110  //-- include nuclear factor (shadowing / anti-shadowing / ...)
111  GetParam( "IncludeNuclMod", fIncludeNuclMod ) ;
112 
113  //-- Use 2016 SF relation corrections
114  GetParam( "Use2016Corrections", fUse2016Corrections ) ;
115 
116  //-- Set min for relation between 2xF1 and F2
117  GetParam( "LowQ2CutoffF1F2", fLowQ2CutoffF1F2 ) ;
118 
119  //-- turn charm production off?
120  GetParamDef( "Charm-Prod-Off", fCharmOff, false ) ;
121 
122  //-- weinberg angle
123  double thw ;
124  GetParam( "WeinbergAngle", thw ) ;
125  fSin2thw = TMath::Power(TMath::Sin(thw), 2);
126 
127  LOG("DISSF", pDEBUG) << "Done loading configuration";
128 }
129 //____________________________________________________________________________
131 {
132  // evaluated at:
133  fPDF = new PDF(); // x = computed (+/-corrections) scaling var, Q2
134  fPDFc = new PDF(); // x = computed charm slow re-scaling var, Q2
135 }
136 //____________________________________________________________________________
137 void QPMDISStrucFuncBase::Calculate(const Interaction * interaction) const
138 {
139  // Reset mutable members
140  fF1 = 0;
141  fF2 = 0;
142  fF3 = 0;
143  fF4 = 0;
144  fF5 = 0;
145  fF6 = 0;
146 
147  // Get process info & perform various checks
148  const ProcessInfo & proc_info = interaction->ProcInfo();
149  const InitialState & init_state = interaction->InitState();
150  const Target & tgt = init_state.Tgt();
151 
152  int nuc_pdgc = tgt.HitNucPdg();
153  int probe_pdgc = init_state.ProbePdg();
154  bool is_p = pdg::IsProton ( nuc_pdgc );
155  bool is_n = pdg::IsNeutron ( nuc_pdgc );
156  bool is_nu = pdg::IsNeutrino ( probe_pdgc );
157  bool is_nubar = pdg::IsAntiNeutrino ( probe_pdgc );
158  bool is_lepton = pdg::IsLepton ( probe_pdgc );
159  bool is_dm = pdg::IsDarkMatter ( probe_pdgc );
160  bool is_CC = proc_info.IsWeakCC();
161  bool is_NC = proc_info.IsWeakNC();
162  bool is_EM = proc_info.IsEM();
163  bool is_dmi = proc_info.IsDarkMatter();
164 
165  if ( !is_lepton && !is_dm ) return;
166  if ( !is_p && !is_n ) return;
167  if ( tgt.N() == 0 && is_n ) return;
168  if ( tgt.Z() == 0 && is_p ) return;
169 
170  // Flags switching on/off quark contributions so that this algorithm can be
171  // used for both l + N -> l' + X, and l + q -> l' + q' level calculations
172 
173  double switch_uv = 1.;
174  double switch_us = 1.;
175  double switch_ubar = 1.;
176  double switch_dv = 1.;
177  double switch_ds = 1.;
178  double switch_dbar = 1.;
179  double switch_s = 1.;
180  double switch_sbar = 1.;
181  double switch_c = 1.;
182  double switch_cbar = 1.;
183 
184  if(tgt.HitQrkIsSet()) {
185 
186  switch_uv = 0.;
187  switch_us = 0.;
188  switch_ubar = 0.;
189  switch_dv = 0.;
190  switch_ds = 0.;
191  switch_dbar = 0.;
192  switch_s = 0.;
193  switch_sbar = 0.;
194  switch_c = 0.;
195  switch_cbar = 0.;
196 
197  int qpdg = tgt.HitQrkPdg();
198  bool sea = tgt.HitSeaQrk();
199 
200  bool is_u = pdg::IsUQuark (qpdg);
201  bool is_ubar = pdg::IsAntiUQuark (qpdg);
202  bool is_d = pdg::IsDQuark (qpdg);
203  bool is_dbar = pdg::IsAntiDQuark (qpdg);
204  bool is_s = pdg::IsSQuark (qpdg);
205  bool is_sbar = pdg::IsAntiSQuark (qpdg);
206  bool is_c = pdg::IsCQuark (qpdg);
207  bool is_cbar = pdg::IsAntiCQuark (qpdg);
208 
209  if (!sea && is_u ) { switch_uv = 1; }
210  else if ( sea && is_u ) { switch_us = 1; }
211  else if ( sea && is_ubar) { switch_ubar = 1; }
212  else if (!sea && is_d ) { switch_dv = 1; }
213  else if ( sea && is_d ) { switch_ds = 1; }
214  else if ( sea && is_dbar) { switch_dbar = 1; }
215  else if ( sea && is_s ) { switch_s = 1; }
216  else if ( sea && is_sbar) { switch_sbar = 1; }
217  else if ( sea && is_c ) { switch_c = 1; }
218  else if ( sea && is_cbar) { switch_cbar = 1; }
219  else return;
220 
221  // make sure user inputs make sense
222  if(is_nu && is_CC && is_u ) return;
223  if(is_nu && is_CC && is_c ) return;
224  if(is_nu && is_CC && is_dbar) return;
225  if(is_nu && is_CC && is_sbar) return;
226  if(is_nubar && is_CC && is_ubar) return;
227  if(is_nubar && is_CC && is_cbar) return;
228  if(is_nubar && is_CC && is_d ) return;
229  if(is_nubar && is_CC && is_s ) return;
230  }
231 
232  // Compute PDFs [both at (scaling-var,Q2) and (slow-rescaling-var,Q2)
233  // Applying all PDF K-factors abd scaling variable corrections
234 
235  this -> CalcPDFs (interaction);
236 
237  //
238  // Compute structure functions for the EM, NC and CC cases
239  //
240 
241  double F2val=0, xF3val=0;
242 
243  // *** NEUTRAL CURRENT
244 
245  // Include DM in NC
246  if(is_NC || is_dmi) {
247 
248  if(!is_nu && !is_nubar && !is_dm) return;
249 
250  double GL = (is_nu) ? ( 0.5 - (2./3.)*fSin2thw) : ( - (2./3.)*fSin2thw); // clu
251  double GR = (is_nu) ? ( - (2./3.)*fSin2thw) : ( 0.5 - (2./3.)*fSin2thw); // cru
252  double GLp = (is_nu) ? (-0.5 + (1./3.)*fSin2thw) : ( (1./3.)*fSin2thw); // cld
253  double GRp = (is_nu) ? ( (1./3.)*fSin2thw) : (-0.5 + (1./3.)*fSin2thw); // crd
254  // Set the couplings to up and down quarks to be axial for DM
255  if (is_dm) {
256  GL = -1.;
257  GR = 1.;
258  GLp = -1.;
259  GRp = 1.;
260  }
261  double gvu = GL + GR;
262  double gau = GL - GR;
263  double gvd = GLp + GRp;
264  double gad = GLp - GRp;
265  double gvu2 = TMath::Power(gvu, 2.);
266  double gau2 = TMath::Power(gau, 2.);
267  double gvd2 = TMath::Power(gvd, 2.);
268  double gad2 = TMath::Power(gad, 2.);
269 
270  double q2 = (switch_uv * fuv + switch_us * fus + switch_c * fc) * (gvu2+gau2) +
271  (switch_dv * fdv + switch_ds * fds + switch_s * fs) * (gvd2+gad2);
272  double q3 = (switch_uv * fuv + switch_us * fus + switch_c * fc) * (2*gvu*gau) +
273  (switch_dv * fdv + switch_ds * fds + switch_s * fs) * (2*gvd*gad);
274 
275  double qb2 = (switch_ubar * fus + switch_cbar * fc) * (gvu2+gau2) +
276  (switch_dbar * fds + switch_sbar * fs) * (gvd2+gad2);
277  double qb3 = (switch_ubar * fus + switch_cbar * fc) * (2*gvu*gau) +
278  (switch_dbar * fds + switch_sbar * fs) * (2*gvd*gad);
279 
280 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
281  LOG("DISSF", pINFO) << "f2 : q = " << q2 << ", bar{q} = " << qb2;
282  LOG("DISSF", pINFO) << "xf3: q = " << q3 << ", bar{q} = " << qb3;
283 #endif
284 
285  F2val = q2+qb2;
286  xF3val = q3-qb3;
287  }
288 
289  // *** CHARGED CURRENT
290 
291  if(is_CC) {
292  double q=0, qbar=0;
293 
294  if (is_nu) {
295  q = ( switch_dv * fdv + switch_ds * fds ) * fVud2 +
296  ( switch_s * fs ) * fVus2 +
297  ( switch_dv * fdv_c + switch_ds * fds_c ) * fVcd2 +
298  ( switch_s * fs_c ) * fVcs2;
299 
300  qbar = ( switch_ubar * fus ) * fVud2 +
301  ( switch_ubar * fus ) * fVus2 +
302  ( switch_cbar * fc_c ) * fVcd2 +
303  ( switch_cbar * fc_c ) * fVcs2;
304  }
305  else
306  if (is_nubar) {
307  q = ( switch_uv * fuv + switch_us * fus ) * fVud2 +
308  ( switch_uv * fuv + switch_us * fus ) * fVus2 +
309  ( switch_c * fc_c ) * fVcd2 +
310  ( switch_c * fc_c ) * fVcs2;
311 
312  qbar = ( switch_dbar * fds_c ) * fVcd2 +
313  ( switch_dbar * fds ) * fVud2 +
314  ( switch_sbar * fs ) * fVus2 +
315  ( switch_sbar * fs_c ) * fVcs2;
316  }
317  else {
318  return;
319  }
320 
321  F2val = 2*(q+qbar);
322  xF3val = 2*(q-qbar);
323  }
324 
325  // *** ELECTROMAGNETIC
326 
327  if(is_EM) {
328 
329  if(!pdg::IsChargedLepton(probe_pdgc)) return;
330 
331  double sq23 = TMath::Power(2./3., 2.);
332  double sq13 = TMath::Power(1./3., 2.);
333 
334  double qu = sq23 * ( switch_uv * fuv + switch_us * fus );
335  double qd = sq13 * ( switch_dv * fdv + switch_ds * fds );
336  double qs = sq13 * ( switch_s * fs );
337  double qbu = sq23 * ( switch_ubar * fus );
338  double qbd = sq13 * ( switch_dbar * fds );
339  double qbs = sq13 * ( switch_sbar * fs );
340 
341  double q = qu + qd + qs;
342  double qbar = qbu + qbd + qbs;
343 
344  F2val = q + qbar;;
345  xF3val = 0.;
346 
347  }
348 
349  double Q2val = this->Q2 (interaction);
350  double x = this->ScalingVar(interaction);
351  double f = this->NuclMod (interaction); // nuclear modification
352  double r = this->R (interaction); // R ~ FL
353 
354 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
355  LOG("DISSF", pDEBUG) << "Nucl. mod = " << f;
356  LOG("DISSF", pDEBUG) << "R(=FL/2xF1) = " << r;
357 #endif
358 
359  if(fUse2016Corrections) {
360  //It was confirmed by A.Bodek that the modified scaling variable
361  //should just be used to compute the strucure functions F2 and xF3,
362  //but that the usual Bjorken x should be used for the relations
363  //between the structure functions.
364  //For the same reason remove the freezing of Q2 at 0.8 for those relations,
365  //although it has not been explicitly asked to A.Bodek if it should be done.
366 
367  const Kinematics & kinematics = interaction->Kine();
368  double bjx = kinematics.x();
369 
370  double a = TMath::Power(bjx,2.) / TMath::Max(Q2val, fLowQ2CutoffF1F2);
371  double c = (1. + 4. * kNucleonMass2 * a) / (1.+r);
372 
373  fF3 = f * xF3val/bjx;
374  fF2 = f * F2val;
375  fF1 = fF2 * 0.5*c/bjx;
376  fF5 = fF2/bjx; // Albright-Jarlskog relation
377  fF4 = 0.; // Nucl.Phys.B 84, 467 (1975)
378  }
379  else {
380  double a = TMath::Power(x,2.) / TMath::Max(Q2val, fLowQ2CutoffF1F2);
381  double c = (1. + 4. * kNucleonMass2 * a) / (1.+r);
382  //double a = TMath::Power(x,2.) / Q2val;
383  //double c = (1. + 4. * kNucleonMass * a) / (1.+r);
384 
385  fF3 = f * xF3val / x;
386  fF2 = f * F2val;
387  fF1 = fF2 * 0.5 * c / x;
388  fF5 = fF2 / x; // Albright-Jarlskog relation
389  fF4 = 0.; // Nucl.Phys.B 84, 467 (1975)
390  }
391 
392 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
393  LOG("DISSF", pDEBUG)
394  << "F1-F5 = "
395  << fF1 << ", " << fF2 << ", " << fF3 << ", " << fF4 << ", " << fF5;
396 #endif
397 }
398 //____________________________________________________________________________
399 double QPMDISStrucFuncBase::Q2(const Interaction * interaction) const
400 {
401 // Return Q2 from the kinematics or, if not set, compute it from x,y
402 // The x might be corrected
403 
404  const Kinematics & kinematics = interaction->Kine();
405 
406  // if Q2 (or q2) is set then prefer this value
407  if (kinematics.KVSet(kKVQ2) || kinematics.KVSet(kKVq2)) {
408  double Q2val = kinematics.Q2();
409  return Q2val;
410  }
411  // if Q2 was not set, then compute it from x,y,Ev,Mnucleon
412  if (kinematics.KVSet(kKVy)) {
413  const InitialState & init_state = interaction->InitState();
414  double Mn = init_state.Tgt().HitNucP4Ptr()->M(); // could be off-shell
415  //double x = this->ScalingVar(interaction); // could be redefined
416  double x = kinematics.x();
417  double y = kinematics.y();
418  double Ev = init_state.ProbeE(kRfHitNucRest);
419  double Q2val = 2*Mn*Ev*x*y;
420  return Q2val;
421  }
422  LOG("DISSF", pERROR) << "Could not compute Q2!";
423  return 0;
424 }
425 //____________________________________________________________________________
426 double QPMDISStrucFuncBase::ScalingVar(const Interaction* interaction) const
427 {
428 // The scaling variable is set to the normal Bjorken x.
429 // Override DISStructureFuncModel::ScalingVar() to compute corrections
430 
431  return interaction->Kine().x();
432 }
433 //____________________________________________________________________________
435  double & kuv, double & kdv, double & kus, double & kds) const
436 {
437 // This is an abstract class: no model-specific correction
438 // The PDF scaling variables are set to 1
439 // Override this method to compute model-dependent corrections
440 
441  kuv = 1.;
442  kdv = 1.;
443  kus = 1.;
444  kds = 1.;
445 }
446 //____________________________________________________________________________
447 double QPMDISStrucFuncBase::NuclMod(const Interaction * interaction) const
448 {
449 // Nuclear modification to Fi
450 // The scaling variable can be overwritten to include corrections
451 
452  if( interaction->TestBit(kIAssumeFreeNucleon) ) return 1.0;
453  if( interaction->TestBit(kINoNuclearCorrection) ) return 1.0;
454 
455  double f = 1.;
456  if(fIncludeNuclMod) {
457  const Target & tgt = interaction->InitState().Tgt();
458 
459 // The x used for computing the DIS Nuclear correction factor should be the
460 // experimental x, not the rescaled x or off-shell-rest-frame version of x
461 // (i.e. selected x). Since we do not have access to experimental x at this
462 // point in the calculation, just use selected x.
463  const Kinematics & kine = interaction->Kine();
464  double x = kine.x();
465  int A = tgt.A();
467 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
468  LOG("DISSF", pDEBUG) << "Nuclear factor for x of " << x << " = " << f;
469 #endif
470  }
471 
472  return f;
473 }
474 //____________________________________________________________________________
475 double QPMDISStrucFuncBase::R(const Interaction * interaction) const
476 {
477 // Computes R ( ~ longitudinal structure function FL = R * 2xF1)
478 // The scaling variable can be overwritten to include corrections
479 
480 // The x used for computing the DIS Nuclear correction factor should be the
481 // experimental x, not the rescaled x or off-shell-rest-frame version of x
482 // (i.e. selected x). Since we do not have access to experimental x at this
483 // point in the calculation, just use selected x.
484  if(fIncludeR) {
485  const Kinematics & kine = interaction->Kine();
486  double x = kine.x();
487 // double x = this->ScalingVar(interaction);
488  double Q2val = this->Q2(interaction);
489  double Rval = utils::phys::RWhitlow(x, Q2val);
490  return Rval;
491  }
492  return 0;
493 }
494 //____________________________________________________________________________
495 void QPMDISStrucFuncBase::CalcPDFs(const Interaction * interaction) const
496 {
497  // Clean-up previous calculation
498  fPDF -> Reset();
499  fPDFc -> Reset();
500 
501  // Get the kinematical variables x,Q2 (could include corrections)
502  double x = this->ScalingVar(interaction);
503  double Q2val = this->Q2(interaction);
504 
505  // Get the hit nucleon mass (could be off-shell)
506  const Target & tgt = interaction->InitState().Tgt();
507  double M = tgt.HitNucP4().M();
508 
509  // Get the Q2 for which PDFs will be evaluated
510  double Q2pdf = TMath::Max(Q2val, fQ2min);
511 
512  // Compute PDFs at (x,Q2)
513 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
514  LOG("DISSF", pDEBUG) << "Calculating PDFs @ x = " << x << ", Q2 = " << Q2pdf;
515 #endif
516  fPDF->Calculate(x, Q2pdf);
517 
518  // Check whether it is above charm threshold
519  bool above_charm =
521  if(above_charm) {
522 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
523  LOG("DISSF", pDEBUG)
524  << "The event is above the charm threshold (mcharm = " << fMc << ")";
525 #endif
526  if(fCharmOff) {
527  LOG("DISSF", pINFO) << "Charm production is turned off";
528  } else {
529  // compute the slow rescaling var
530  double xc = utils::kinematics::SlowRescalingVar(x, Q2val, M, fMc);
531  if(xc<0 || xc>1) {
532  LOG("DISSF", pINFO) << "Unphys. slow rescaling var: xc = " << xc;
533  } else {
534  // compute PDFs at (xc,Q2)
535 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
536  LOG("DISSF", pDEBUG)
537  << "Calculating PDFs @ xc (slow rescaling) = " << x << ", Q2 = " << Q2val;
538 #endif
539  fPDFc->Calculate(xc, Q2pdf);
540  }
541  }// charm off?
542  }//above charm thr?
543  else {
544  LOG("DISSF", pDEBUG)
545  << "The event is below the charm threshold (mcharm = " << fMc << ")";
546  }
547 
548  // Compute the K factors
549  double kval_u = 1.;
550  double kval_d = 1.;
551  double ksea_u = 1.;
552  double ksea_d = 1.;
553 
554  this->KFactors(interaction, kval_u, kval_d, ksea_u, ksea_d);
555 
556 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
557  LOG("DISSF", pDEBUG) << "K-Factors:";
558  LOG("DISSF", pDEBUG) << "U: Kval = " << kval_u << ", Ksea = " << ksea_u;
559  LOG("DISSF", pDEBUG) << "D: Kval = " << kval_d << ", Ksea = " << ksea_d;
560 #endif
561 
562  // Apply the K factors
563  //
564  // Always scale d pdfs with d kfactors and u pdfs with u kfactors.
565  // Don't swap the applied kfactors for neutrons.
566  // Debdatta & Donna noted (Sep.2006) that a similar swap in the neugen
567  // implementation was the cause of the difference in nu and nubar F2
568  //
569  fPDF->ScaleUpValence (kval_u);
570  fPDF->ScaleDownValence (kval_d);
571  fPDF->ScaleUpSea (ksea_u);
572  fPDF->ScaleDownSea (ksea_d);
573  fPDF->ScaleStrange (ksea_d);
574  fPDF->ScaleCharm (ksea_u);
575  if(above_charm) {
576  fPDFc->ScaleUpValence (kval_u);
577  fPDFc->ScaleDownValence (kval_d);
578  fPDFc->ScaleUpSea (ksea_u);
579  fPDFc->ScaleDownSea (ksea_d);
580  fPDFc->ScaleStrange (ksea_d);
581  fPDFc->ScaleCharm (ksea_u);
582  }
583 
584  // Rules of thumb
585  // ---------------------------------------
586  // - For W+ exchange use: -1/3|e| quarks and -2/3|e| antiquarks
587  // - For W- exchange use: 2/3|e| quarks and 1/3|e| antiquarks
588  // - For each qi -> qj transition multiply with the (ij CKM element)^2
589  // - Use isospin symmetry to get neutron's u,d from proton's u,d
590  // -- neutron d = proton u
591  // -- neutron u = proton d
592  // - Use u = usea + uvalence. Same for d
593  // - For s,c use q=qbar
594  // - For t,b use q=qbar=0
595 
596  fuv = fPDF -> UpValence();
597  fus = fPDF -> UpSea();
598  fdv = fPDF -> DownValence();
599  fds = fPDF -> DownSea();
600  fs = fPDF -> Strange();
601  fc = 0.;
602  fuv_c = fPDFc -> UpValence(); // will be 0 if < charm threshold
603  fus_c = fPDFc -> UpSea(); // ...
604  fdv_c = fPDFc -> DownValence(); // ...
605  fds_c = fPDFc -> DownSea(); // ...
606  fs_c = fPDFc -> Strange(); // ...
607  fc_c = fPDFc -> Charm(); // ...
608 
609  // The above are the proton parton density function. Get the PDFs for the
610  // hit nucleon (p or n) by swapping u<->d if necessary
611 
612  int nuc_pdgc = tgt.HitNucPdg();
613  bool isP = pdg::IsProton (nuc_pdgc);
614  bool isN = pdg::IsNeutron (nuc_pdgc);
615  assert(isP || isN);
616 
617  double tmp = 0;
618  if (isN) { // swap u <-> d
619  tmp = fuv; fuv = fdv; fdv = tmp;
620  tmp = fus; fus = fds; fds = tmp;
621  tmp = fuv_c; fuv_c = fdv_c; fdv_c = tmp;
622  tmp = fus_c; fus_c = fds_c; fds_c = tmp;
623  }
624 
625 }
626 //____________________________________________________________________________
Pure Abstract Base Class. Defines the DISStructureFuncModelI interface to be implemented by any algor...
bool HitSeaQrk(void) const
Definition: Target.cxx:299
bool IsWeakCC(void) const
double fVcd
CKM element Vcd used.
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
double fQ2min
min Q^2 allowed for PDFs: PDF(Q2&lt;Q2min):=PDF(Q2min)
bool IsUQuark(int pdgc)
Definition: PDGUtils.cxx:266
int HitNucPdg(void) const
Definition: Target.cxx:304
int A(void) const
Definition: Target.h:70
PDF * fPDFc
computed PDFs @ (slow-rescaling-var,Q2)
int HitQrkPdg(void) const
Definition: Target.cxx:242
A class to store PDFs.
Definition: PDF.h:37
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double fVud
CKM element Vud used.
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
double x(bool selected=false) const
Definition: Kinematics.cxx:99
bool fIncludeNuclMod
include nuclear factor (shadowing, anti-shadowing,...)?
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:276
double fMc
charm mass used
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:306
Pure abstract base class. Defines the PDFModelI interface to be implemented by wrapper classes to exi...
Definition: PDFModelI.h:28
bool IsAntiDQuark(int pdgc)
Definition: PDGUtils.cxx:301
PDF * fPDF
computed PDFs @ (x,Q2)
bool fIncludeR
include R (~FL) in DIS SF calculation?
void ScaleDownValence(double kscale)
Definition: PDF.cxx:86
double y(bool selected=false) const
Definition: Kinematics.cxx:112
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
virtual double NuclMod(const Interaction *i) const
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
double SlowRescalingVar(double x, double Q2, double M, double mc)
Definition: KineUtils.cxx:1257
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
bool IsWeakNC(void) const
const UInt_t kINoNuclearCorrection
if set, inhibit nuclear corrections
Definition: Interaction.h:51
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool KVSet(KineVar_t kv) const
Definition: Kinematics.cxx:317
virtual double Q2(const Interaction *i) const
static constexpr double A
Definition: Units.h:74
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const double a
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
int Z(void) const
Definition: Target.h:68
bool IsAboveCharmThreshold(double x, double Q2, double M, double mc)
Definition: KineUtils.cxx:1238
bool fCharmOff
turn charm production off?
#define pINFO
Definition: Messenger.h:62
static const double kNucleonMass2
virtual double ScalingVar(const Interaction *i) const
void Configure(const Registry &config)
void Calculate(double x, double q2)
Definition: PDF.cxx:49
double fLowQ2CutoffF1F2
Set min for relation between 2xF1 and F2.
bool IsEM(void) const
virtual void Calculate(const Interaction *interaction) const
Calculate the structure functions F1-F6 for the input interaction.
double fVus
CKM element Vcs used.
double fVcs
CKM element Vcs used.
virtual void KFactors(const Interaction *i, double &kuv, double &kdv, double &kus, double &kds) const
void ScaleUpValence(double kscale)
Definition: PDF.cxx:81
int N(void) const
Definition: Target.h:69
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
bool IsCQuark(int pdgc)
Definition: PDGUtils.cxx:281
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
bool IsDarkMatter(void) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
bool IsAntiCQuark(int pdgc)
Definition: PDGUtils.cxx:311
virtual double R(const Interaction *i) const
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:271
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
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:66
virtual void CalcPDFs(const Interaction *i) const
bool fUse2016Corrections
Use 2016 SF relation corrections.
void ScaleUpSea(double kscale)
Definition: PDF.cxx:91
void ScaleCharm(double kscale)
Definition: PDF.cxx:106
double DISNuclFactor(double x, int A)
double RWhitlow(double x, double Q2)
Definition: PhysUtils.cxx:75
void ScaleDownSea(double kscale)
Definition: PDF.cxx:96
double ProbeE(RefFrame_t rf) const
bool IsLepton(int pdgc)
Definition: PDGUtils.cxx:86
bool IsAntiUQuark(int pdgc)
Definition: PDGUtils.cxx:296
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
void ScaleStrange(double kscale)
Definition: PDF.cxx:101
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345