GSLXSecFunc.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  Changes required to implement the GENIE Boosted Dark Matter module
11  were installed by Josh Berger (Univ. of Wisconsin)
12 */
13 //____________________________________________________________________________
14 
15 #include <cassert>
16 
17 #include <TMath.h>
18 
31 
32 using namespace genie;
33 
34 //____________________________________________________________________________
36  const XSecAlgorithmI * m, const Interaction * i) :
37 ROOT::Math::IBaseFunctionOneDim(),
38 fModel(m),
39 fInteraction(i)
40 {
41 
42 }
44 {
45 
46 }
47 unsigned int genie::utils::gsl::dXSec_dQ2_E::NDim(void) const
48 {
49  return 1;
50 }
52 {
53 // inputs:
54 // Q2 [GeV^2]
55 // outputs:
56 // differential cross section [10^-38 cm^2 / GeV^2]
57 //
58  double Q2 = xin;
59  fInteraction->KinePtr()->SetQ2(Q2);
60  double xsec = fModel->XSec(fInteraction, kPSQ2fE);
61 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
62  LOG("GSLXSecFunc", pDEBUG) << "xsec(Q2 = " << Q2 << ") = " << xsec;
63 #endif
64  return xsec/(1E-38 * units::cm2);
65 }
66 ROOT::Math::IBaseFunctionOneDim *
68 {
69  return
71 }
72 //____________________________________________________________________________
74  const XSecAlgorithmI * m, const Interaction * i) :
75 ROOT::Math::IBaseFunctionOneDim(),
76 fModel(m),
77 fInteraction(i)
78 {
79 
80 }
82 {
83 
84 }
85 unsigned int genie::utils::gsl::dXSec_dy_E::NDim(void) const
86 {
87  return 1;
88 }
89 double genie::utils::gsl::dXSec_dy_E::DoEval(double xin) const
90 {
91 // inputs:
92 // y [-]
93 // outputs:
94 // differential cross section [10^-38 cm^2]
95 //
96  double y = xin;
97  fInteraction->KinePtr()->Sety(y);
98  double xsec = fModel->XSec(fInteraction, kPSyfE);
99 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
100  LOG("GXSecFunc", pDEBUG) << "xsec(y = " << y << ") = " << xsec;
101 #endif
102  return xsec/(1E-38 * units::cm2);
103 }
104 ROOT::Math::IBaseFunctionOneDim *
106 {
107  return
109 }
110 //____________________________________________________________________________
112  const XSecAlgorithmI * m, const Interaction * i) :
113 ROOT::Math::IBaseFunctionMultiDim(),
114 fModel(m),
115 fInteraction(i)
116 {
117 
118 }
120 {
121 
122 }
124 {
125  return 2;
126 }
127 double genie::utils::gsl::d2XSec_dxdy_E::DoEval(const double * xin) const
128 {
129 // inputs:
130 // x [-]
131 // y [-]
132 // outputs:
133 // differential cross section [10^-38 cm^2]
134 //
135  double x = xin[0];
136  double y = xin[1];
137  fInteraction->KinePtr()->Setx(x);
138  fInteraction->KinePtr()->Sety(y);
140  double xsec = fModel->XSec(fInteraction, kPSxyfE);
141  return xsec/(1E-38 * units::cm2);
142 }
143 ROOT::Math::IBaseFunctionMultiDim *
145 {
146  return
148 }
149 //____________________________________________________________________________
151  const XSecAlgorithmI * m, const Interaction * i) :
152 ROOT::Math::IBaseFunctionMultiDim(),
153 fModel(m),
154 fInteraction(i)
155 {
156 
157 }
159 {
160 
161 }
163 {
164  return 2;
165 }
166 double genie::utils::gsl::d2XSec_dQ2dy_E::DoEval(const double * xin) const
167 {
168 // inputs:
169 // Q2 [-]
170 // y [-]
171 // outputs:
172 // differential cross section [10^-38 cm^2]
173 //
174  double Q2 = xin[0];
175  double y = xin[1];
176  fInteraction->KinePtr()->SetQ2(Q2);
177  fInteraction->KinePtr()->Sety(y);
179  double xsec = fModel->XSec(fInteraction, kPSQ2yfE);
180  return xsec/(1E-38 * units::cm2);
181 }
182 ROOT::Math::IBaseFunctionMultiDim *
184 {
185  return
187 }
188 //____________________________________________________________________________
190  const XSecAlgorithmI * m, const Interaction * i) :
191 ROOT::Math::IBaseFunctionMultiDim(),
192 fModel(m),
193 fInteraction(i)
194 {
195 
196 }
198 {
199 
200 }
202 {
203  return 3;
204 }
205 double genie::utils::gsl::d2XSec_dQ2dydt_E::DoEval(const double * xin) const
206 {
207 // inputs:
208 // Q2 [-]
209 // y [-]
210 // t [-]
211 // outputs:
212 // differential cross section [10^-38 cm^2]
213 //
214  //double E = fInteraction->InitState().ProbeE(kRfLab);
215  double Q2 = xin[0];
216  double y = xin[1];
217  double t = xin[2];
218  fInteraction->KinePtr()->SetQ2(Q2);
219  fInteraction->KinePtr()->Sety(y);
220  fInteraction->KinePtr()->Sett(t);
222  double xsec = fModel->XSec(fInteraction, kPSQ2yfE);
223  return xsec/(1E-38 * units::cm2);
224 }
225 ROOT::Math::IBaseFunctionMultiDim *
227 {
228  return
230 }
231 //____________________________________________________________________________
233  const XSecAlgorithmI * m, const Interaction * i) :
234 ROOT::Math::IBaseFunctionMultiDim(),
235 fModel(m),
236 fInteraction(i)
237 {
238 
239 }
241 {
242 
243 }
245 {
246  return 3;
247 }
248 double genie::utils::gsl::d3XSec_dxdydt_E::DoEval(const double * xin) const
249 {
250 // inputs:
251 // x [-]
252 // y [-]
253 // t [-]
254 // outputs:
255 // differential cross section [10^-38 cm^2]
256 //
257  //double E = fInteraction->InitState().ProbeE(kRfLab);
258  double x = xin[0];
259  double y = xin[1];
260  double t = xin[2];
261  fInteraction->KinePtr()->Setx(x);
262  fInteraction->KinePtr()->Sety(y);
263  fInteraction->KinePtr()->Sett(t);
264  double xsec = fModel->XSec(fInteraction, kPSxytfE);
265  return xsec/(1E-38 * units::cm2);
266 }
267 ROOT::Math::IBaseFunctionMultiDim *
269 {
270  return
272 }
273 //____________________________________________________________________________
275  const XSecAlgorithmI * m, const Interaction * i) :
276 ROOT::Math::IBaseFunctionMultiDim(),
277 fModel(m),
278 fInteraction(i)
279 {
280 
281 }
283 {
284 
285 }
287 {
288  return 2;
289 }
290 double genie::utils::gsl::d2XSec_dWdQ2_E::DoEval(const double * xin) const
291 {
292 // inputs:
293 // W [GeV]
294 // Q2 [GeV^2]
295 // outputs:
296 // differential cross section [10^-38 cm^2/GeV^3]
297 //
298  double W = xin[0];
299  double Q2 = xin[1];
300  fInteraction->KinePtr()->SetW(W);
301  fInteraction->KinePtr()->SetQ2(Q2);
304  double x=0,y=0;
306  double M = fInteraction->InitState().Tgt().HitNucP4Ptr()->M();
307 
308  kinematics::WQ2toXY(E,M,W,Q2,x,y);
309  fInteraction->KinePtr()->Setx(x);
310  fInteraction->KinePtr()->Sety(y);
311  }
312  double xsec = fModel->XSec(fInteraction, kPSWQ2fE);
313  return xsec/(1E-38 * units::cm2);
314 }
315 ROOT::Math::IBaseFunctionMultiDim *
317 {
318  return
320 }
321 //____________________________________________________________________________
323  const XSecAlgorithmI * m, const Interaction * i, double x) :
324 ROOT::Math::IBaseFunctionOneDim(),
325 fModel(m),
326 fInteraction(i),
327 fx(x)
328 {
329 
330 }
332 {
333 
334 }
336 {
337  return 1;
338 }
340 {
341 // inputs:
342 // y [-]
343 // outputs:
344 // differential cross section [10^-38 cm^2]
345 //
346  double y = xin;
348  fInteraction->KinePtr()->Sety(y);
349  double xsec = fModel->XSec(fInteraction, kPSxyfE);
350  return xsec/(1E-38 * units::cm2);
351 }
352 ROOT::Math::IBaseFunctionOneDim *
354 {
355  return
357 }
358 //____________________________________________________________________________
360  const XSecAlgorithmI * m, const Interaction * i, double y) :
361 ROOT::Math::IBaseFunctionOneDim(),
362 fModel(m),
363 fInteraction(i),
364 fy(y)
365 {
366 
367 }
369 {
370 
371 }
373 {
374  return 1;
375 }
377 {
378 // inputs:
379 // x [-]
380 // outputs:
381 // differential cross section [10^-38 cm^2]
382 //
383  double x = xin;
384  fInteraction->KinePtr()->Setx(x);
386  double xsec = fModel->XSec(fInteraction, kPSxyfE);
387  return xsec/(1E-38 * units::cm2);
388 }
389 ROOT::Math::IBaseFunctionOneDim *
391 {
392  return
394 }
395 //____________________________________________________________________________
397  const XSecAlgorithmI * m, const Interaction * i, double W) :
398 ROOT::Math::IBaseFunctionOneDim(),
399 fModel(m),
400 fInteraction(i),
401 fW(W)
402 {
403 
404 }
406 {
407 
408 }
410 {
411  return 1;
412 }
414 {
415 // inputs:
416 // Q2 [GeV^2]
417 // outputs:
418 // differential cross section [10^-38 cm^2/GeV^3]
419 //
420  double Q2 = xin;
422  fInteraction->KinePtr()->SetQ2(Q2);
423  double xsec = fModel->XSec(fInteraction, kPSWQ2fE);
424  return xsec/(1E-38 * units::cm2);
425 }
426 ROOT::Math::IBaseFunctionOneDim *
428 {
429  return
431 }
432 //____________________________________________________________________________
434  const XSecAlgorithmI * m, const Interaction * i, double Q2) :
435 ROOT::Math::IBaseFunctionOneDim(),
436 fModel(m),
437 fInteraction(i),
438 fQ2(Q2)
439 {
440 
441 }
443 {
444 
445 }
447 {
448  return 1;
449 }
451 {
452 // inputs:
453 // W [GeV]
454 // outputs:
455 // differential cross section [10^-38 cm^2/GeV^3]
456 //
457  double W = xin;
458  fInteraction->KinePtr()->SetW(W);
460  double xsec = fModel->XSec(fInteraction,kPSWQ2fE);
461  return xsec/(1E-38 * units::cm2);
462 }
463 ROOT::Math::IBaseFunctionOneDim *
465 {
466  return
468 }
469 //____________________________________________________________________________
470 //
471 // This just returns the 5-D differential cross section
472 //
474  const XSecAlgorithmI * m, const Interaction * i) :
475 ROOT::Math::IBaseFunctionMultiDim(),
476 fModel(m),
477 fInteraction(i),
478 flip(false)
479 {
480 
481 }
483 {
484 }
485 
486 unsigned int genie::utils::gsl::d5XSecAR::NDim(void) const
487 {
488  return 5;
489 }
490 double genie::utils::gsl::d5XSecAR::DoEval(const double * xin) const
491 {
492 // inputs:
493 // x [-]
494 // y [-]
495 // outputs:
496 // differential cross section [10^-38 cm^2]
497 //
498 
500  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
501  double E_nu = P4_nu->E();
502 
503  double E_l = xin[0];
504  double theta_l = xin[1];
505  double phi_l = xin[2];
506  double theta_pi = xin[3];
507  double phi_pi = xin[4];
508 
509  double E_pi= E_nu-E_l;
510 
511  double y = E_pi/E_nu;
512 
513  double m_l = fInteraction->FSPrimLepton()->Mass();
514  if (E_l < m_l) {
515  return 0.;
516  }
517 
518  double m_pi;
519  if ( fInteraction->ProcInfo().IsWeakCC() ) {
520  m_pi = constants::kPionMass;
521  }
522  else {
523  m_pi = constants::kPi0Mass;
524  }
525 
526 
527  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
528  TVector3 lepton_3vector = TVector3(0,0,0);
529  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
530  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
531 
532  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
533  TVector3 pion_3vector = TVector3(0,0,0);
534  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
535  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
536 
537  double Q2 = -(*P4_nu-P4_lep).Mag2();
538 
539  double x = Q2/(2*E_pi*constants::kNucleonMass);
540 
542 
543  if ( x < xlim.min || x > xlim.max ) {
544  return 0.;
545  }
546 
547  kinematics->Setx(x);
548  kinematics->Sety(y);
550 
551  kinematics->SetFSLeptonP4(P4_lep );
552  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
553 
554  double xsec = fModel->XSec(fInteraction);
555  if (xsec>0 && flip) {
556  xsec = xsec*-1.0;
557  }
558  delete P4_nu;
559  //return xsec/(1E-38 * units::cm2);
560  return xsec;
561 }
562 
563 ROOT::Math::IBaseFunctionMultiDim *
565 {
566  return
568 }
569 
570 //____________________________________________________________________________
571 //
572 // This is the original 5-D cross-section that Steve D coded
573 //
575  const XSecAlgorithmI * m, const Interaction * i) :
576 ROOT::Math::IBaseFunctionMultiDim(),
577 fModel(m),
578 fInteraction(i)
579 {
580 
581 }
583 {
584 
585 }
587 {
588  return 5;
589 }
591 {
592 // inputs:
593 // x [-]
594 // y [-]
595 // outputs:
596 // differential cross section [10^-38 cm^2]
597 //
599  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
600  double E_nu = P4_nu->E();
601 
602  double E_l = xin[0];
603  double theta_l = xin[1];
604  double phi_l = xin[2];
605  double theta_pi = xin[3];
606  double phi_pi = xin[4];
607 
608  double E_pi= E_nu-E_l;
609 
610  double y = E_pi/E_nu;
611 
612  double m_l = fInteraction->FSPrimLepton()->Mass();
613  if (E_l < m_l) {
614  return 0.;
615  }
616 
617  double m_pi;
618  if ( fInteraction->ProcInfo().IsWeakCC() ) {
619  m_pi = constants::kPionMass;
620  }
621  else {
622  m_pi = constants::kPi0Mass;
623  }
624 
625  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
626  TVector3 lepton_3vector = TVector3(0,0,0);
627  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
628  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
629 
630  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
631  TVector3 pion_3vector = TVector3(0,0,0);
632  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
633  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
634 
635  double Q2 = -(*P4_nu-P4_lep).Mag2();
636 
637  double x = Q2/(2*E_pi*constants::kNucleonMass);
638 
640 
641  if ( x < xlim.min || x > xlim.max ) {
642  return 0.;
643  }
644 
645  kinematics->Setx(x);
646  kinematics->Sety(y);
648 
649  kinematics->SetFSLeptonP4(P4_lep );
650  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
651 
652  delete P4_nu;
653 
654  double xsec = fModel->XSec(fInteraction)*TMath::Sin(theta_l)*TMath::Sin(theta_pi);
655  return xsec/(1E-38 * units::cm2);
656 }
657 
658 ROOT::Math::IBaseFunctionMultiDim *
660 {
661  return
663 }
664 
665 //____________________________________________________________________________
666 //
667 // This is the same as the 5d space that Steve D coded,
668 // But we remove the integration of phi_l
669 //
670 
672  const XSecAlgorithmI * m, const Interaction * i) :
673 ROOT::Math::IBaseFunctionMultiDim(),
674 fModel(m),
675 fInteraction(i)
676 {
677 
678 }
680 {
681 
682 }
684 {
685  return 4;
686 }
688 {
689 // inputs:
690 // El [GeV]
691 // theta l [rad]
692 // theta pi [rad]
693 // phi pi [rad]
694 // outputs:
695 // differential cross section [10^-38 cm^2]
696 //
698  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
699  double E_nu = P4_nu->E();
700 
701  double E_l = xin[0];
702  double theta_l = xin[1];
703  double phi_l = 0.0;
704  double theta_pi = xin[2];
705  double phi_pi = xin[3];
706 
707  double sin_theta_l = TMath::Sin(theta_l);
708  double sin_theta_pi = TMath::Sin(theta_pi);
709 
710  double E_pi= E_nu-E_l;
711 
712  double y = E_pi/E_nu;
713 
714  double m_l = fInteraction->FSPrimLepton()->Mass();
715  if (E_l < m_l) {
716  return 0.;
717  }
718 
719  double m_pi;
720  if ( fInteraction->ProcInfo().IsWeakCC() ) {
721  m_pi = constants::kPionMass;
722  }
723  else {
724  m_pi = constants::kPi0Mass;
725  }
726 
727  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
728  TVector3 lepton_3vector = TVector3(0,0,0);
729  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
730  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
731 
732  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
733  TVector3 pion_3vector = TVector3(0,0,0);
734  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
735  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
736 
737  double Q2 = -(*P4_nu-P4_lep).Mag2();
738 
739  double x = Q2/(2*E_pi*constants::kNucleonMass);
740 
742 
743  if ( x < xlim.min || x > xlim.max ) {
744  return 0.;
745  }
746 
747  kinematics->Setx(x);
748  kinematics->Sety(y);
750 
751  kinematics->SetFSLeptonP4(P4_lep );
752  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
753 
754  delete P4_nu;
755 
756  double xsec = sin_theta_l * sin_theta_pi * fModel->XSec(fInteraction,kPSElOlTpifE);
757  return fFactor * xsec/(1E-38 * units::cm2);
758 }
759 ROOT::Math::IBaseFunctionMultiDim *
761 {
762  return
764 }
766 {
767  fFactor = factor;
768 }
770 {
771  return fFactor;
772 }
773 //____________________________________________________________________________
775  const XSecAlgorithmI * m, const Interaction * i) :
776 ROOT::Math::IBaseFunctionMultiDim(),
777 fModel(m),
778 fInteraction(i),
779 fElep(-1)
780 {
781 
782 }
784 {
785 
786 }
788 {
789  return 3;
790 }
792 {
793 // inputs:
794 // theta l [rad]
795 // phi_l [rad]
796 // phi pi [rad]
797 // outputs:
798 // differential cross section [10^-38 cm^2]
799 //
801  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
802  double E_nu = P4_nu->E();
803 
804  double E_l = fElep;
805 
806  double theta_l = xin[0];
807  double phi_l = xin[1];
808  double theta_pi = xin[2];
809  double phi_pi = 0;
810 
811  double sin_theta_l = TMath::Sin(theta_l);
812  double sin_theta_pi = TMath::Sin(theta_pi);
813 
814  double E_pi= E_nu-E_l;
815 
816  double y = E_pi/E_nu;
817 
818  double m_l = fInteraction->FSPrimLepton()->Mass();
819  if (E_l < m_l) {
820  return 0.;
821  }
822 
823  double m_pi;
824  if ( fInteraction->ProcInfo().IsWeakCC() ) {
825  m_pi = constants::kPionMass;
826  }
827  else {
828  m_pi = constants::kPi0Mass;
829  }
830 
831  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
832  TVector3 lepton_3vector = TVector3(0,0,0);
833  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
834  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
835 
836  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
837  TVector3 pion_3vector = TVector3(0,0,0);
838  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
839  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
840 
841  double Q2 = -(*P4_nu-P4_lep).Mag2();
842 
843  double x = Q2/(2*E_pi*constants::kNucleonMass);
844 
846 
847  if ( x < xlim.min || x > xlim.max ) {
848  return 0.;
849  }
850 
851  kinematics->Setx(x);
852  kinematics->Sety(y);
854 
855  kinematics->SetFSLeptonP4(P4_lep );
856  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
857 
858  delete P4_nu;
859 
860  double xsec = (sin_theta_l * sin_theta_pi) * fModel->XSec(fInteraction,kPSElOlTpifE);
861  return xsec/(1E-38 * units::cm2);
862 }
865 {
867  out->SetE_lep(fElep);
868  return out;
869 }
870 //____________________________________________________________________________
872 {
873  fElep = E_lepton;
874 }
875 //____________________________________________________________________________
877  const XSecAlgorithmI * m, const Interaction * i,
878  string gsl_nd_integrator_type, double gsl_relative_tolerance,
879  unsigned int max_n_calls) :
880 ROOT::Math::IBaseFunctionOneDim(),
881 fModel(m),
882 fInteraction(i),
883 integrator(utils::gsl::IntegrationNDimTypeFromString(gsl_nd_integrator_type),1, gsl_relative_tolerance, max_n_calls),
884 fGSLIntegratorType(gsl_nd_integrator_type),
885 fGSLRelTol(gsl_relative_tolerance),
886 fGSLMaxCalls(max_n_calls)
887 {
889 
890  integrator.SetRelTolerance(gsl_relative_tolerance);
891  integrator.SetFunction(*func);
892 
894 
897 }
899 {
900  delete func;
901 }
903 {
904  double Elep = xin;
905  func->SetE_lep(Elep);
906  double xsec = integrator.Integral(&kine_min[0], &kine_max[0]) ;
907  LOG("GSLXSecFunc",pINFO) << "dXSec_dElep_AR >> "<<func->NDim()<<"d integral done. (Elep = " <<Elep<< " , dxsec/dElep = "<<xsec << ")";
908  return xsec;
909 }
912 {
913  return
915 }
916 //____________________________________________________________________________
918  const ROOT::Math::IBaseFunctionMultiDim * fn,
919  bool * ifLog, double * mins, double * maxes) :
920  fFn(fn),
921  fIfLog(ifLog),
922  fMins(mins),
923  fMaxes(maxes)
924 {
925 }
927 {
928 }
929 
930 // ROOT::Math::IBaseFunctionMultiDim interface
932 {
933  return fFn->NDim();
934 }
935 double genie::utils::gsl::dXSec_Log_Wrapper::DoEval (const double * xin) const
936 {
937  double * toEval = new double[this->NDim()];
938  double a,b,x;
939  for (unsigned int i = 0 ; i < this->NDim() ; i++ )
940  {
941  if (fIfLog[i]) {
942  a = fMins[i];
943  b = fMaxes[i];
944  x = xin[i];
945  toEval[i] = a + (b-a)/(constants::kNapierConst-1.) * (exp(x/(b-a)) - 1.);
946  }
947  else {
948  toEval[i] = xin[i];
949  }
950  }
951  double val = (*fFn)(toEval);
952  delete[] toEval;
953  return val;
954 }
955 ROOT::Math::IBaseFunctionMultiDim * genie::utils::gsl::dXSec_Log_Wrapper::Clone (void) const
956 {
958 }
959 
960 //____________________________________________________________________________
961 
Cross Section Calculation Interface.
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:71
double DoEval(const double *xin) const
unsigned int NDim(void) const
Definition: GSLXSecFunc.cxx:85
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
static const double kNapierConst
Definition: Constants.h:44
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:92
double DoEval(const double *xin) const
double DoEval(double xin) const
Definition: GSLXSecFunc.cxx:89
bool IsWeakCC(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:156
dXSec_dQ2_E(const XSecAlgorithmI *m, const Interaction *i)
Definition: GSLXSecFunc.cxx:35
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
d2XSec_dQ2dydt_E(const XSecAlgorithmI *m, const Interaction *i)
double DoEval(double xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
double DoEval(const double *xin) const
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
unsigned int NDim(void) const
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
d3Xsec_dOmegaldThetapi * Clone(void) const
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
double DoEval(double xin) const
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
A simple [min,max] interval for doubles.
Definition: Range1.h:43
static const double kPi0Mass
Definition: Constants.h:75
unsigned int NDim(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:220
double DoEval(double xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double DoEval(const double *xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:50
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:241
d5XSecAR(const XSecAlgorithmI *m, const Interaction *i)
double DoEval(const double *xin) const
unsigned int NDim(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:72
const Interaction * fInteraction
Definition: GSLXSecFunc.h:264
static const double cm2
Definition: Units.h:77
double DoEval(double xin) const
Definition: GSLXSecFunc.cxx:51
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:113
double DoEval(double xin) const
double DoEval(const double *xin) const
d3Xsec_dOmegaldThetapi(const XSecAlgorithmI *m, const Interaction *i)
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
double DoEval(double xin) const
d2XSec_dxdy_Ex(const XSecAlgorithmI *m, const Interaction *i, double x)
unsigned int NDim(void) const
Definition: GSLXSecFunc.cxx:47
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:197
Summary information for an interaction.
Definition: Interaction.h:56
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
unsigned int NDim(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
Definition: GSLXSecFunc.cxx:67
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
Float_t E
Definition: plot.C:20
static keras::KerasModel * fModel
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:263
d3XSec_dxdydt_E(const XSecAlgorithmI *m, const Interaction *i)
const double a
unsigned int NDim(void) const
dXSec_Log_Wrapper(const ROOT::Math::IBaseFunctionMultiDim *fn, bool *ifLog, double *min, double *maxes)
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:176
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
d2XSec_dWdQ2_EQ2(const XSecAlgorithmI *m, const Interaction *i, double Q2)
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:301
double DoEval(const double *xin) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:155
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:134
const Interaction * fInteraction
Definition: GSLXSecFunc.h:114
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
dXSec_dy_E(const XSecAlgorithmI *m, const Interaction *i)
Definition: GSLXSecFunc.cxx:73
static const double kASmallNum
Definition: Controls.h:40
void SetE_lep(double E_lepton) const
#define pINFO
Definition: Messenger.h:63
unsigned int NDim(void) const
d2XSec_dQ2dy_E(const XSecAlgorithmI *m, const Interaction *i)
const Interaction * fInteraction
Definition: GSLXSecFunc.h:93
Double_t xsec[nknots]
Definition: testXsec.C:47
unsigned int NDim(void) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const Interaction * fInteraction
Definition: GSLXSecFunc.h:287
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:386
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:82
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
d2XSec_dxdy_Ey(const XSecAlgorithmI *m, const Interaction *i, double x)
float Mag2() const
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1185
double max
Definition: Range1.h:54
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
const Interaction * fInteraction
Definition: GSLXSecFunc.h:242
const Interaction * fInteraction
Definition: GSLXSecFunc.h:198
d2XSec_dxdy_E(const XSecAlgorithmI *m, const Interaction *i)
static const double kPionMass
Definition: Constants.h:74
const genie::utils::gsl::d3Xsec_dOmegaldThetapi * func
Definition: GSLXSecFunc.h:388
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double DoEval(const double *xin) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:286
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:219
d5Xsec_dEldOmegaldOmegapi(const XSecAlgorithmI *m, const Interaction *i)
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
void UpdateXFromQ2Y(const Interaction *in)
Definition: KineUtils.cxx:1221
unsigned int NDim(void) const
d4Xsec_dEldThetaldOmegapi(const XSecAlgorithmI *m, const Interaction *i)
const Interaction * fInteraction
Definition: GSLXSecFunc.h:135
const hit & b
Definition: hits.cxx:21
unsigned int NDim(void) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
d2XSec_dWdQ2_E(const XSecAlgorithmI *m, const Interaction *i)
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:317
const Interaction * fInteraction
Definition: GSLXSecFunc.h:177
double DoEval(const double *xin) const
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
Range1D_t XLim(void) const
x limits
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
unsigned int NDim(void) const
double DoEval(const double *xin) const
const ROOT::Math::IBaseFunctionMultiDim * fFn
Definition: GSLXSecFunc.h:418
const Target & Tgt(void) const
Definition: InitialState.h:67
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
d2XSec_dWdQ2_EW(const XSecAlgorithmI *m, const Interaction *i, double W)
const Interaction * fInteraction
Definition: GSLXSecFunc.h:51
ROOT::Math::IntegratorMultiDim integrator
Definition: GSLXSecFunc.h:390
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:87
void kinematics()
Definition: kinematics.C:10
#define W(x)
dXSec_dElep_AR * Clone(void) const
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:38
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:385
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
#define pDEBUG
Definition: Messenger.h:64
unsigned int NDim(void) const