All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
GSLXSecFunc.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2020, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
7  University of Liverpool & STFC Rutherford Appleton Laboratory
8 
9  Changes required to implement the GENIE Boosted Dark Matter module
10  were installed by Josh Berger (Univ. of Wisconsin)
11 
12  Changes required to implement the GENIE Dark Neutrino module
13  were installed by Iker de Icaza (Univ. of Sussex)
14 */
15 //____________________________________________________________________________
16 
17 #include <cassert>
18 
19 #include <TMath.h>
20 
23 #include "Framework/Conventions/GBuild.h"
33 
34 using namespace genie;
35 
36 //____________________________________________________________________________
38  const XSecAlgorithmI * m, const Interaction * i, double scale) :
39 ROOT::Math::IBaseFunctionOneDim(),
40 fModel(m),
41 fInteraction(i),
42 fScale(scale)
43 {
44 
45 }
47 {
48 
49 }
50 unsigned int genie::utils::gsl::dXSec_dQ2_E::NDim(void) const
51 {
52  return 1;
53 }
55 {
56 // inputs:
57 // Q2 [GeV^2]
58 // outputs:
59 // differential cross section [10^-38 cm^2 / GeV^2]
60 //
61  double Q2 = xin;
62  fInteraction->KinePtr()->SetQ2(Q2);
63  double xsec = fModel->XSec(fInteraction, kPSQ2fE);
64 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
65  LOG("GSLXSecFunc", pDEBUG) << "xsec(Q2 = " << Q2 << ") = " << xsec;
66 #endif
67  return fScale*xsec/(1E-38 * units::cm2);
68 }
69 ROOT::Math::IBaseFunctionOneDim *
71 {
72  return
74 }
75 //____________________________________________________________________________
77  const XSecAlgorithmI * m, const Interaction * i) :
78 ROOT::Math::IBaseFunctionOneDim(),
79 fModel(m),
80 fInteraction(i)
81 {
82 
83 }
85 {
86 
87 }
88 unsigned int genie::utils::gsl::dXSec_dy_E::NDim(void) const
89 {
90  return 1;
91 }
92 double genie::utils::gsl::dXSec_dy_E::DoEval(double xin) const
93 {
94 // inputs:
95 // y [-]
96 // outputs:
97 // differential cross section [10^-38 cm^2]
98 //
99  double y = xin;
100  fInteraction->KinePtr()->Sety(y);
101  double xsec = fModel->XSec(fInteraction, kPSyfE);
102 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
103  LOG("GXSecFunc", pDEBUG) << "xsec(y = " << y << ") = " << xsec;
104 #endif
105  return xsec/(1E-38 * units::cm2);
106 }
107 ROOT::Math::IBaseFunctionOneDim *
109 {
110  return
112 }
113 //____________________________________________________________________________
115  const XSecAlgorithmI * m, const Interaction * i,
116  double DNuMass, double scale) :
117  ROOT::Math::IBaseFunctionOneDim(),
118  fModel(m),
119  fInteraction(i),
120  fDNuMass(DNuMass),
121  fScale(scale)
122 {
123 
124 }
126 {
127 
128 }
130 {
131  return 1;
132 }
134 {
135 // inputs:
136 // DNuEnergy [GeV]
137 // outputs:
138 // differential cross section [10^-38 cm^2 / GeV]
139 //
140 
141  double DNuEnergy = xin;
142  double fDNuMass2 = fDNuMass*fDNuMass;
143 
144  Kinematics * kinematics = fInteraction->KinePtr();
145  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
146  double E_nu = P4_nu->E();
147  double M_target = fInteraction->InitState().Tgt().Mass();
148 
149  double ETimesM = E_nu * M_target;
150  double EPlusM = E_nu + M_target;
151 
152  double p_DNu = TMath::Sqrt(DNuEnergy*DNuEnergy - fDNuMass2);
153  double cos_theta_DNu = (DNuEnergy*(EPlusM) - ETimesM - 0.5*fDNuMass2) / (E_nu * p_DNu);
154  double theta_DNu = TMath::ACos(cos_theta_DNu);
155  TVector3 DNu_3vector = TVector3(0,0,0);
156  DNu_3vector.SetMagThetaPhi(p_DNu, theta_DNu, 0.);
157  TLorentzVector P4_DNu = TLorentzVector(DNu_3vector, DNuEnergy);
158  kinematics->SetFSLeptonP4(P4_DNu);
159 
160  TVector3 target_3vector = P4_nu->Vect() - DNu_3vector;
161  double E_target = E_nu + M_target - DNuEnergy;
162  TLorentzVector P4_target = TLorentzVector(target_3vector , E_target);
163  kinematics->SetHadSystP4(P4_target);
164  kinematics->SetQ2(2.*M_target*(E_target-M_target));
165 
166  delete P4_nu;
167  double xsec = fModel->XSec(fInteraction, kPSEDNufE);
168 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
169  LOG("GSLXSecFunc", pDEBUG) << "xsec(DNuEnergy = " << DNuEnergy << ") = " << xsec;
170 #endif
171 
172  return fScale*xsec/(1E-38 * units::cm2);
173 }
174 ROOT::Math::IBaseFunctionOneDim *
176 {
177  return
179 }
181 {
182  // look at valid angles section at tech note
183  const double E = fInteraction->InitState().ProbeE(kRfLab);
184  const double M = fInteraction->InitState().Tgt().Mass();
185  const double M2 = M * M;
186  double fDNuMass2 = fDNuMass*fDNuMass;
187 
188  const double A = M2 + 2.*M*E;
189  const double B = (M+E) * (E*M + 0.5*fDNuMass2);
190  const double C = E*E *(M2 + fDNuMass2) + E*M*fDNuMass2 + 0.25*fDNuMass2*fDNuMass2;
191  const double D = sqrt(B*B - A*C);
192 
193  Range1D_t DNuEnergy((B - D)/A, (B + D)/A);
194  return DNuEnergy;
195 
196 }
197 //____________________________________________________________________________
199  const XSecAlgorithmI * m, const Interaction * i) :
200 ROOT::Math::IBaseFunctionMultiDim(),
201 fModel(m),
202 fInteraction(i)
203 {
204 
205 }
207 {
208 
209 }
211 {
212  return 2;
213 }
214 double genie::utils::gsl::d2XSec_dxdy_E::DoEval(const double * xin) const
215 {
216 // inputs:
217 // x [-]
218 // y [-]
219 // outputs:
220 // differential cross section [10^-38 cm^2]
221 //
222  double x = xin[0];
223  double y = xin[1];
224  fInteraction->KinePtr()->Setx(x);
225  fInteraction->KinePtr()->Sety(y);
227  double xsec = fModel->XSec(fInteraction, kPSxyfE);
228  return xsec/(1E-38 * units::cm2);
229 }
230 ROOT::Math::IBaseFunctionMultiDim *
232 {
233  return
235 }
236 //____________________________________________________________________________
238  const XSecAlgorithmI * m, const Interaction * i) :
239 ROOT::Math::IBaseFunctionMultiDim(),
240 fModel(m),
241 fInteraction(i)
242 {
243 
244 }
246 {
247 
248 }
250 {
251  return 2;
252 }
253 double genie::utils::gsl::d2XSec_dQ2dy_E::DoEval(const double * xin) const
254 {
255 // inputs:
256 // Q2 [-]
257 // y [-]
258 // outputs:
259 // differential cross section [10^-38 cm^2]
260 //
261  double Q2 = xin[0];
262  double y = xin[1];
263  fInteraction->KinePtr()->SetQ2(Q2);
264  fInteraction->KinePtr()->Sety(y);
266  double xsec = fModel->XSec(fInteraction, kPSQ2yfE);
267  return xsec/(1E-38 * units::cm2);
268 }
269 ROOT::Math::IBaseFunctionMultiDim *
271 {
272  return
274 }
275 //____________________________________________________________________________
277  const XSecAlgorithmI * m, const Interaction * i, double scale) :
278 ROOT::Math::IBaseFunctionMultiDim(),
279 fModel(m),
280 fInteraction(i),
281 fScale(scale)
282 {
283 
284 }
286 {
287 
288 }
290 {
291  return 2;
292 }
294 {
295 // inputs:
296 // log10(x) [-]
297 // log10(Q2) [-]
298 // outputs:
299 // differential cross section [10^-38 cm^2]
300 //
301  fInteraction->KinePtr()->Setx(TMath::Power(10,xin[0]));
302  fInteraction->KinePtr()->SetQ2(TMath::Power(10,xin[1]));
304  double xsec = fModel->XSec(fInteraction, kPSlog10xlog10Q2fE);
305  return fScale*xsec/(1E-38 * units::cm2);
306 }
307 ROOT::Math::IBaseFunctionMultiDim *
309 {
310  return
312 }
313 //____________________________________________________________________________
315  const XSecAlgorithmI * m, const Interaction * i) :
316 ROOT::Math::IBaseFunctionMultiDim(),
317 fModel(m),
318 fInteraction(i)
319 {
320 
321 }
323 {
324 
325 }
327 {
328  return 3;
329 }
330 double genie::utils::gsl::d2XSec_dQ2dydt_E::DoEval(const double * xin) const
331 {
332 // inputs:
333 // Q2 [-]
334 // y [-]
335 // t [-]
336 // outputs:
337 // differential cross section [10^-38 cm^2]
338 //
339  //double E = fInteraction->InitState().ProbeE(kRfLab);
340  double Q2 = xin[0];
341  double y = xin[1];
342  double t = xin[2];
343  fInteraction->KinePtr()->SetQ2(Q2);
344  fInteraction->KinePtr()->Sety(y);
345  fInteraction->KinePtr()->Sett(t);
347  double xsec = fModel->XSec(fInteraction, kPSQ2yfE);
348  return xsec/(1E-38 * units::cm2);
349 }
350 ROOT::Math::IBaseFunctionMultiDim *
352 {
353  return
355 }
356 //____________________________________________________________________________
358  const XSecAlgorithmI * m, const Interaction * i) :
359 ROOT::Math::IBaseFunctionMultiDim(),
360 fModel(m),
361 fInteraction(i)
362 {
363 
364 }
366 {
367 
368 }
370 {
371  return 3;
372 }
373 double genie::utils::gsl::d3XSec_dxdydt_E::DoEval(const double * xin) const
374 {
375 // inputs:
376 // x [-]
377 // y [-]
378 // t [-]
379 // outputs:
380 // differential cross section [10^-38 cm^2]
381 //
382  //double E = fInteraction->InitState().ProbeE(kRfLab);
383  double x = xin[0];
384  double y = xin[1];
385  double t = xin[2];
386  fInteraction->KinePtr()->Setx(x);
387  fInteraction->KinePtr()->Sety(y);
388  fInteraction->KinePtr()->Sett(t);
389  double xsec = fModel->XSec(fInteraction, kPSxytfE);
390  return xsec/(1E-38 * units::cm2);
391 }
392 ROOT::Math::IBaseFunctionMultiDim *
394 {
395  return
397 }
398 //____________________________________________________________________________
400  const XSecAlgorithmI * m, const Interaction * i) :
401 ROOT::Math::IBaseFunctionMultiDim(),
402 fModel(m),
403 fInteraction(i)
404 {
405 
406 }
408 {
409 
410 }
412 {
413  return 2;
414 }
415 double genie::utils::gsl::d2XSec_dWdQ2_E::DoEval(const double * xin) const
416 {
417 // inputs:
418 // W [GeV]
419 // Q2 [GeV^2]
420 // outputs:
421 // differential cross section [10^-38 cm^2/GeV^3]
422 //
423  double W = xin[0];
424  double Q2 = xin[1];
425  fInteraction->KinePtr()->SetW(W);
426  fInteraction->KinePtr()->SetQ2(Q2);
429  double x=0,y=0;
431  double M = fInteraction->InitState().Tgt().HitNucP4Ptr()->M();
432 
433  kinematics::WQ2toXY(E,M,W,Q2,x,y);
434  fInteraction->KinePtr()->Setx(x);
435  fInteraction->KinePtr()->Sety(y);
436  }
437  double xsec = fModel->XSec(fInteraction, kPSWQ2fE);
438  return xsec/(1E-38 * units::cm2);
439 }
440 ROOT::Math::IBaseFunctionMultiDim *
442 {
443  return
445 }
446 //____________________________________________________________________________
448  const XSecAlgorithmI * m, const Interaction * i, double x) :
449 ROOT::Math::IBaseFunctionOneDim(),
450 fModel(m),
451 fInteraction(i),
452 fx(x)
453 {
454 
455 }
457 {
458 
459 }
461 {
462  return 1;
463 }
465 {
466 // inputs:
467 // y [-]
468 // outputs:
469 // differential cross section [10^-38 cm^2]
470 //
471  double y = xin;
473  fInteraction->KinePtr()->Sety(y);
474  double xsec = fModel->XSec(fInteraction, kPSxyfE);
475  return xsec/(1E-38 * units::cm2);
476 }
477 ROOT::Math::IBaseFunctionOneDim *
479 {
480  return
482 }
483 //____________________________________________________________________________
485  const XSecAlgorithmI * m, const Interaction * i, double y) :
486 ROOT::Math::IBaseFunctionOneDim(),
487 fModel(m),
488 fInteraction(i),
489 fy(y)
490 {
491 
492 }
494 {
495 
496 }
498 {
499  return 1;
500 }
502 {
503 // inputs:
504 // x [-]
505 // outputs:
506 // differential cross section [10^-38 cm^2]
507 //
508  double x = xin;
509  fInteraction->KinePtr()->Setx(x);
511  double xsec = fModel->XSec(fInteraction, kPSxyfE);
512  return xsec/(1E-38 * units::cm2);
513 }
514 ROOT::Math::IBaseFunctionOneDim *
516 {
517  return
519 }
520 //____________________________________________________________________________
522  const XSecAlgorithmI * m, const Interaction * i, double W) :
523 ROOT::Math::IBaseFunctionOneDim(),
524 fModel(m),
525 fInteraction(i),
526 fW(W)
527 {
528 
529 }
531 {
532 
533 }
535 {
536  return 1;
537 }
539 {
540 // inputs:
541 // Q2 [GeV^2]
542 // outputs:
543 // differential cross section [10^-38 cm^2/GeV^3]
544 //
545  double Q2 = xin;
547  fInteraction->KinePtr()->SetQ2(Q2);
548  double xsec = fModel->XSec(fInteraction, kPSWQ2fE);
549  return xsec/(1E-38 * units::cm2);
550 }
551 ROOT::Math::IBaseFunctionOneDim *
553 {
554  return
556 }
557 //____________________________________________________________________________
559  const XSecAlgorithmI * m, const Interaction * i, double Q2) :
560 ROOT::Math::IBaseFunctionOneDim(),
561 fModel(m),
562 fInteraction(i),
563 fQ2(Q2)
564 {
565 
566 }
568 {
569 
570 }
572 {
573  return 1;
574 }
576 {
577 // inputs:
578 // W [GeV]
579 // outputs:
580 // differential cross section [10^-38 cm^2/GeV^3]
581 //
582  double W = xin;
583  fInteraction->KinePtr()->SetW(W);
585  double xsec = fModel->XSec(fInteraction,kPSWQ2fE);
586  return xsec/(1E-38 * units::cm2);
587 }
588 ROOT::Math::IBaseFunctionOneDim *
590 {
591  return
593 }
594 //____________________________________________________________________________
595 //
596 // This just returns the 5-D differential cross section
597 //
599  const XSecAlgorithmI * m, const Interaction * i) :
600 ROOT::Math::IBaseFunctionMultiDim(),
601 fModel(m),
602 fInteraction(i),
603 flip(false)
604 {
605 
606 }
608 {
609 }
610 
611 unsigned int genie::utils::gsl::d5XSecAR::NDim(void) const
612 {
613  return 5;
614 }
615 double genie::utils::gsl::d5XSecAR::DoEval(const double * xin) const
616 {
617 // inputs:
618 // x [-]
619 // y [-]
620 // outputs:
621 // differential cross section [10^-38 cm^2]
622 //
623 
624  Kinematics * kinematics = fInteraction->KinePtr();
625  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
626  double E_nu = P4_nu->E();
627 
628  double E_l = xin[0];
629  double theta_l = xin[1];
630  double phi_l = xin[2];
631  double theta_pi = xin[3];
632  double phi_pi = xin[4];
633 
634  double E_pi= E_nu-E_l;
635 
636  double y = E_pi/E_nu;
637 
638  double m_l = fInteraction->FSPrimLepton()->Mass();
639  if (E_l < m_l) {
640  return 0.;
641  }
642 
643  double m_pi;
644  if ( fInteraction->ProcInfo().IsWeakCC() ) {
645  m_pi = constants::kPionMass;
646  }
647  else {
648  m_pi = constants::kPi0Mass;
649  }
650 
651 
652  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
653  TVector3 lepton_3vector = TVector3(0,0,0);
654  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
655  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
656 
657  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
658  TVector3 pion_3vector = TVector3(0,0,0);
659  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
660  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
661 
662  double Q2 = -(*P4_nu-P4_lep).Mag2();
663 
664  double x = Q2/(2*E_pi*constants::kNucleonMass);
665 
667 
668  if ( x < xlim.min || x > xlim.max ) {
669  return 0.;
670  }
671 
672  kinematics->Setx(x);
673  kinematics->Sety(y);
675 
676  kinematics->SetFSLeptonP4(P4_lep );
677  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
678 
679  double xsec = fModel->XSec(fInteraction);
680  if (xsec>0 && flip) {
681  xsec = xsec*-1.0;
682  }
683  delete P4_nu;
684  //return xsec/(1E-38 * units::cm2);
685  return xsec;
686 }
687 
688 ROOT::Math::IBaseFunctionMultiDim *
690 {
691  return
693 }
694 
695 //____________________________________________________________________________
696 //
697 // This is the original 5-D cross-section that Steve D coded
698 //
700  const XSecAlgorithmI * m, const Interaction * i) :
701 ROOT::Math::IBaseFunctionMultiDim(),
702 fModel(m),
703 fInteraction(i)
704 {
705 
706 }
708 {
709 
710 }
712 {
713  return 5;
714 }
716 {
717 // inputs:
718 // x [-]
719 // y [-]
720 // outputs:
721 // differential cross section [10^-38 cm^2]
722 //
723  Kinematics * kinematics = fInteraction->KinePtr();
724  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
725  double E_nu = P4_nu->E();
726 
727  double E_l = xin[0];
728  double theta_l = xin[1];
729  double phi_l = xin[2];
730  double theta_pi = xin[3];
731  double phi_pi = xin[4];
732 
733  double E_pi= E_nu-E_l;
734 
735  double y = E_pi/E_nu;
736 
737  double m_l = fInteraction->FSPrimLepton()->Mass();
738  if (E_l < m_l) {
739  return 0.;
740  }
741 
742  double m_pi;
743  if ( fInteraction->ProcInfo().IsWeakCC() ) {
744  m_pi = constants::kPionMass;
745  }
746  else {
747  m_pi = constants::kPi0Mass;
748  }
749 
750  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
751  TVector3 lepton_3vector = TVector3(0,0,0);
752  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
753  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
754 
755  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
756  TVector3 pion_3vector = TVector3(0,0,0);
757  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
758  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
759 
760  double Q2 = -(*P4_nu-P4_lep).Mag2();
761 
762  double x = Q2/(2*E_pi*constants::kNucleonMass);
763 
765 
766  if ( x < xlim.min || x > xlim.max ) {
767  return 0.;
768  }
769 
770  kinematics->Setx(x);
771  kinematics->Sety(y);
773 
774  kinematics->SetFSLeptonP4(P4_lep );
775  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
776 
777  delete P4_nu;
778 
779  double xsec = fModel->XSec(fInteraction)*TMath::Sin(theta_l)*TMath::Sin(theta_pi);
780  return xsec/(1E-38 * units::cm2);
781 }
782 
783 ROOT::Math::IBaseFunctionMultiDim *
785 {
786  return
788 }
789 
790 //____________________________________________________________________________
791 //
792 // This is the same as the 5d space that Steve D coded,
793 // But we remove the integration of phi_l
794 //
795 
797  const XSecAlgorithmI * m, const Interaction * i) :
798 ROOT::Math::IBaseFunctionMultiDim(),
799 fModel(m),
800 fInteraction(i)
801 {
802 
803 }
805 {
806 
807 }
809 {
810  return 4;
811 }
813 {
814 // inputs:
815 // El [GeV]
816 // theta l [rad]
817 // theta pi [rad]
818 // phi pi [rad]
819 // outputs:
820 // differential cross section [10^-38 cm^2]
821 //
822  Kinematics * kinematics = fInteraction->KinePtr();
823  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
824  double E_nu = P4_nu->E();
825 
826  double E_l = xin[0];
827  double theta_l = xin[1];
828  double phi_l = 0.0;
829  double theta_pi = xin[2];
830  double phi_pi = xin[3];
831 
832  double sin_theta_l = TMath::Sin(theta_l);
833  double sin_theta_pi = TMath::Sin(theta_pi);
834 
835  double E_pi= E_nu-E_l;
836 
837  double y = E_pi/E_nu;
838 
839  double m_l = fInteraction->FSPrimLepton()->Mass();
840  if (E_l < m_l) {
841  return 0.;
842  }
843 
844  double m_pi;
845  if ( fInteraction->ProcInfo().IsWeakCC() ) {
846  m_pi = constants::kPionMass;
847  }
848  else {
849  m_pi = constants::kPi0Mass;
850  }
851 
852  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
853  TVector3 lepton_3vector = TVector3(0,0,0);
854  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
855  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
856 
857  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
858  TVector3 pion_3vector = TVector3(0,0,0);
859  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
860  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
861 
862  double Q2 = -(*P4_nu-P4_lep).Mag2();
863 
864  double x = Q2/(2*E_pi*constants::kNucleonMass);
865 
867 
868  if ( x < xlim.min || x > xlim.max ) {
869  return 0.;
870  }
871 
872  kinematics->Setx(x);
873  kinematics->Sety(y);
875 
876  kinematics->SetFSLeptonP4(P4_lep );
877  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
878 
879  delete P4_nu;
880 
881  double xsec = sin_theta_l * sin_theta_pi * fModel->XSec(fInteraction,kPSElOlTpifE);
882  return fFactor * xsec/(1E-38 * units::cm2);
883 }
884 ROOT::Math::IBaseFunctionMultiDim *
886 {
887  return
889 }
891 {
892  fFactor = factor;
893 }
895 {
896  return fFactor;
897 }
898 //____________________________________________________________________________
900  const XSecAlgorithmI * m, const Interaction * i) :
901 ROOT::Math::IBaseFunctionMultiDim(),
902 fModel(m),
903 fInteraction(i),
904 fElep(-1)
905 {
906 
907 }
909 {
910 
911 }
913 {
914  return 3;
915 }
917 {
918 // inputs:
919 // theta l [rad]
920 // phi_l [rad]
921 // phi pi [rad]
922 // outputs:
923 // differential cross section [10^-38 cm^2]
924 //
925  Kinematics * kinematics = fInteraction->KinePtr();
926  const TLorentzVector * P4_nu = fInteraction->InitStatePtr()->GetProbeP4(kRfLab);
927  double E_nu = P4_nu->E();
928 
929  double E_l = fElep;
930 
931  double theta_l = xin[0];
932  double phi_l = xin[1];
933  double theta_pi = xin[2];
934  double phi_pi = 0;
935 
936  double sin_theta_l = TMath::Sin(theta_l);
937  double sin_theta_pi = TMath::Sin(theta_pi);
938 
939  double E_pi= E_nu-E_l;
940 
941  double y = E_pi/E_nu;
942 
943  double m_l = fInteraction->FSPrimLepton()->Mass();
944  if (E_l < m_l) {
945  return 0.;
946  }
947 
948  double m_pi;
949  if ( fInteraction->ProcInfo().IsWeakCC() ) {
950  m_pi = constants::kPionMass;
951  }
952  else {
953  m_pi = constants::kPi0Mass;
954  }
955 
956  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
957  TVector3 lepton_3vector = TVector3(0,0,0);
958  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
959  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
960 
961  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
962  TVector3 pion_3vector = TVector3(0,0,0);
963  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
964  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
965 
966  double Q2 = -(*P4_nu-P4_lep).Mag2();
967 
968  double x = Q2/(2*E_pi*constants::kNucleonMass);
969 
971 
972  if ( x < xlim.min || x > xlim.max ) {
973  return 0.;
974  }
975 
976  kinematics->Setx(x);
977  kinematics->Sety(y);
979 
980  kinematics->SetFSLeptonP4(P4_lep );
981  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
982 
983  delete P4_nu;
984 
985  double xsec = (sin_theta_l * sin_theta_pi) * fModel->XSec(fInteraction,kPSElOlTpifE);
986  return xsec/(1E-38 * units::cm2);
987 }
990 {
992  out->SetE_lep(fElep);
993  return out;
994 }
995 //____________________________________________________________________________
997 {
998  fElep = E_lepton;
999 }
1000 //____________________________________________________________________________
1002  const XSecAlgorithmI * m, const Interaction * i,
1003  string gsl_nd_integrator_type, double gsl_relative_tolerance,
1004  unsigned int max_n_calls) :
1005 ROOT::Math::IBaseFunctionOneDim(),
1006 fModel(m),
1007 fInteraction(i),
1008 integrator(utils::gsl::IntegrationNDimTypeFromString(gsl_nd_integrator_type),1, gsl_relative_tolerance, max_n_calls),
1009 fGSLIntegratorType(gsl_nd_integrator_type),
1010 fGSLRelTol(gsl_relative_tolerance),
1011 fGSLMaxCalls(max_n_calls)
1012 {
1014 
1015  integrator.SetRelTolerance(gsl_relative_tolerance);
1016  integrator.SetFunction(*func);
1017 
1019 
1022 }
1024 {
1025  delete func;
1026 }
1028 {
1029  double Elep = xin;
1030  func->SetE_lep(Elep);
1031  double xsec = integrator.Integral(&kine_min[0], &kine_max[0]) ;
1032  LOG("GSLXSecFunc",pINFO) << "dXSec_dElep_AR >> "<<func->NDim()<<"d integral done. (Elep = " <<Elep<< " , dxsec/dElep = "<<xsec << ")";
1033  return xsec;
1034 }
1037 {
1038  return
1040 }
1041 //____________________________________________________________________________
1043  const ROOT::Math::IBaseFunctionMultiDim * fn,
1044  bool * ifLog, double * mins, double * maxes) :
1045  fFn(fn),
1046  fIfLog(ifLog),
1047  fMins(mins),
1048  fMaxes(maxes)
1049 {
1050 }
1052 {
1053 }
1054 //____________________________________________________________________________
1055 
1056 // ROOT::Math::IBaseFunctionMultiDim interface
1058 {
1059  return fFn->NDim();
1060 }
1061 //____________________________________________________________________________
1062 double genie::utils::gsl::dXSec_Log_Wrapper::DoEval (const double * xin) const
1063 {
1064  double * toEval = new double[this->NDim()];
1065  double a,b,x;
1066  for (unsigned int i = 0 ; i < this->NDim() ; i++ )
1067  {
1068  if (fIfLog[i]) {
1069  a = fMins[i];
1070  b = fMaxes[i];
1071  x = xin[i];
1072  toEval[i] = a + (b-a)/(constants::kNapierConst-1.) * (exp(x/(b-a)) - 1.);
1073  }
1074  else {
1075  toEval[i] = xin[i];
1076  }
1077  }
1078  double val = (*fFn)(toEval);
1079  delete[] toEval;
1080  return val;
1081 }
1082 ROOT::Math::IBaseFunctionMultiDim * genie::utils::gsl::dXSec_Log_Wrapper::Clone (void) const
1083 {
1084  return new dXSec_Log_Wrapper(fFn,fIfLog,fMins,fMaxes);
1085 }
1086 
1087 //____________________________________________________________________________
Cross Section Calculation Interface.
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:75
double m_pi
double DoEval(const double *xin) const
unsigned int NDim(void) const
Definition: GSLXSecFunc.cxx:88
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
static const double kNapierConst
Definition: Constants.h:43
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:121
double DoEval(const double *xin) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:99
double DoEval(double xin) const
Definition: GSLXSecFunc.cxx:92
bool IsWeakCC(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:207
Range1D_t IntegrationRange(void) const
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
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:77
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:59
unsigned int NDim(void) const
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
d3Xsec_dOmegaldThetapi * Clone(void) const
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
dXSec_dQ2_E(const XSecAlgorithmI *m, const Interaction *i, double scale=1.)
Definition: GSLXSecFunc.cxx:37
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:42
static const double kPi0Mass
Definition: Constants.h:74
unsigned int NDim(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:271
double DoEval(double xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double DoEval(const double *xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:53
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:292
d5XSecAR(const XSecAlgorithmI *m, const Interaction *i)
double DoEval(const double *xin) const
unsigned int NDim(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:76
double Mass(void) const
Definition: Target.cxx:224
const Interaction * fInteraction
Definition: GSLXSecFunc.h:315
double DoEval(double xin) const
Definition: GSLXSecFunc.cxx:54
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
void UpdateWYFromXQ2(const Interaction *in)
Definition: KineUtils.cxx:1313
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:164
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:50
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:248
dXSec_dEDNu_E(const XSecAlgorithmI *m, const Interaction *i, double DNuMass, double scale=1.)
Summary information for an interaction.
Definition: Interaction.h:56
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
unsigned int NDim(void) const
double DoEval(double xin) const
QAsciiDict< Entry > fn
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
Definition: GSLXSecFunc.cxx:70
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
static constexpr double cm2
Definition: Units.h:69
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:314
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:227
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1119
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:291
double DoEval(const double *xin) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:206
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:185
const Interaction * fInteraction
Definition: GSLXSecFunc.h:165
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
dXSec_dy_E(const XSecAlgorithmI *m, const Interaction *i)
Definition: GSLXSecFunc.cxx:76
static const double kASmallNum
Definition: Controls.h:40
void SetE_lep(double E_lepton) const
#define pINFO
Definition: Messenger.h:62
unsigned int NDim(void) const
d2XSec_dQ2dy_E(const XSecAlgorithmI *m, const Interaction *i)
const Interaction * fInteraction
Definition: GSLXSecFunc.h:122
double DoEval(const double *xin) const
unsigned int NDim(void) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const Interaction * fInteraction
Definition: GSLXSecFunc.h:338
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
unsigned int NDim(void) const
const Interaction * fInteraction
Definition: GSLXSecFunc.h:437
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
d2XSec_dxdy_Ey(const XSecAlgorithmI *m, const Interaction *i, double x)
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:98
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1277
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
const Interaction * fInteraction
Definition: GSLXSecFunc.h:293
const Interaction * fInteraction
Definition: GSLXSecFunc.h:249
d2XSec_dxdy_E(const XSecAlgorithmI *m, const Interaction *i)
static const double kPionMass
Definition: Constants.h:73
const genie::utils::gsl::d3Xsec_dOmegaldThetapi * func
Definition: GSLXSecFunc.h:439
double DoEval(const double *xin) const
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:337
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:270
d5Xsec_dEldOmegaldOmegapi(const XSecAlgorithmI *m, const Interaction *i)
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
void UpdateXFromQ2Y(const Interaction *in)
Definition: KineUtils.cxx:1331
unsigned int NDim(void) const
E
Definition: 018_def.c:13
Definition: utils.py:1
d4Xsec_dEldThetaldOmegapi(const XSecAlgorithmI *m, const Interaction *i)
const Interaction * fInteraction
Definition: GSLXSecFunc.h:186
unsigned int NDim(void) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
d2XSec_dlog10xdlog10Q2_E(const XSecAlgorithmI *m, const Interaction *i, double scale=1.)
#define A
Definition: memgrp.cpp:38
static bool * b
Definition: config.cpp:1043
d2XSec_dWdQ2_E(const XSecAlgorithmI *m, const Interaction *i)
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
const Interaction * fInteraction
Definition: GSLXSecFunc.h:228
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
list x
Definition: train.py:276
unsigned int NDim(void) const
double DoEval(const double *xin) const
const ROOT::Math::IBaseFunctionMultiDim * fFn
Definition: GSLXSecFunc.h:469
const Target & Tgt(void) const
Definition: InitialState.h:66
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
d2XSec_dWdQ2_EW(const XSecAlgorithmI *m, const Interaction *i, double W)
const Interaction * fInteraction
Definition: GSLXSecFunc.h:54
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
ROOT::Math::IntegratorMultiDim integrator
Definition: GSLXSecFunc.h:441
bool IsDarkMatterDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
dXSec_dElep_AR * Clone(void) const
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:37
const XSecAlgorithmI * fModel
Definition: GSLXSecFunc.h:436
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
#define pDEBUG
Definition: Messenger.h:63
unsigned int NDim(void) const