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