GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AlvarezRusoCOHPiPXSec.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  Steve Dennis
7  University of Warwick, Rutherford Appleton Laboratory,
8 */
9 //____________________________________________________________________________
10 
11 #include <iostream>
12 
13 #include <TMath.h>
14 
17 #include "Framework/Conventions/GBuild.h"
27 
32 
33 
34 using namespace genie;
35 using namespace genie::constants;
36 using namespace genie::utils;
37 
38 using namespace alvarezruso;
39 
40 //____________________________________________________________________________
42 XSecAlgorithmI("genie::AlvarezRusoCOHPiPXSec")
43 {
44  fMultidiff = NULL;
45  fLastInteraction = NULL;
46 }
47 //____________________________________________________________________________
49 XSecAlgorithmI("genie::AlvarezRusoCOHPiPXSec", config)
50 {
51  fMultidiff = NULL;
52  fLastInteraction = NULL;
53 }
54 //____________________________________________________________________________
56 {
57  if (fMultidiff) delete fMultidiff;
58 }
59 //____________________________________________________________________________
61  const Interaction * interaction, KinePhaseSpace_t kps) const
62 {
63  if(! this -> ValidProcess (interaction) ) return 0.;
64  if(! this -> ValidKinematics (interaction) ) return 0.;
65 
66  const Kinematics & kinematics = interaction -> Kine();
67  const InitialState & init_state = interaction -> InitState();
68 
69  int A = init_state.Tgt().A(); // mass number
70  int Z = init_state.Tgt().Z(); // atomic number
71  double E_nu = init_state.ProbeE(kRfLab); // neutrino energy
72 
73  const TLorentzVector p4_lep = kinematics.FSLeptonP4();
74  const TLorentzVector p4_pi = kinematics.HadSystP4();
75  double E_lep = p4_lep.E();
76 
77  if (fLastInteraction!=interaction) {
78  if (fMultidiff != NULL) {
79  delete fMultidiff;
80  fMultidiff = NULL;
81  }
82 
83  current_t current;
84  if ( interaction->ProcInfo().IsWeakCC() ) {
85  current = kCC;
86  }
87  else if ( interaction->ProcInfo().IsWeakNC() ) {
88  current = kNC;
89  }
90  else {
91  LOG("AlvarezRusoCohPi",pDEBUG)<<"Unknown current for AlvarezRuso implementation";
92  return 0.;
93  }
94 
95  flavour_t flavour;
96  if ( init_state.ProbePdg() == 12 || init_state.ProbePdg() == -12) {
97  flavour=kE;
98  }
99  else if ( init_state.ProbePdg() == 14 || init_state.ProbePdg() == -14) {
100  flavour=kMu;
101  }
102  else if ( init_state.ProbePdg() == 16 || init_state.ProbePdg() == -16) {
103  flavour=kTau;
104  }
105  else {
106  LOG("AlvarezRusoCohPi",pDEBUG)<<"Unknown probe for AlvarezRuso implementation";
107  return 0.;
108  }
109 
110  nutype_t nutype;
111  if ( init_state.ProbePdg() > 0) {
112  nutype = kNu;
113  } else {
114  nutype = kAntiNu;
115  }
116 
117  fMultidiff = new AlvarezRusoCOHPiPDXSec(Z, A ,current, flavour, nutype);
118  fLastInteraction = interaction;
119  }
120 
121  double xsec = fMultidiff->DXSec(E_nu, E_lep, p4_lep.Theta(), p4_lep.Phi(), p4_pi.Theta(), p4_pi.Phi());
122  xsec = xsec * 1E-38 * units::cm2;
123 
124  if (kps != kPSElOlOpifE) {
125  xsec *= utils::kinematics::Jacobian(interaction, kPSElOlOpifE, kps );
126  }
127 
128  return (xsec);
129 }
130 //____________________________________________________________________________
131 double AlvarezRusoCOHPiPXSec::Integral(const Interaction * interaction) const
132 {
133  double xsec = fXSecIntegrator->Integrate(this,interaction);
134  return xsec;
135 }
136 //____________________________________________________________________________
137 bool AlvarezRusoCOHPiPXSec::ValidProcess(const Interaction * interaction) const
138 {
139  if(interaction->TestBit(kISkipProcessChk)) return true;
140 
141  const InitialState & init_state = interaction->InitState();
142  const ProcessInfo & proc_info = interaction->ProcInfo();
143  const Target & target = init_state.Tgt();
144 
145  int nu = init_state.ProbePdg();
146 
147  if (!proc_info.IsCoherentProduction()) return false;
148  if (!proc_info.IsWeak()) return false;
149  if (target.HitNucIsSet()) return false;
150  if (!(target.A()>1)) return false;
151  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
152 
153  return true;
154 }
155 //____________________________________________________________________________
157 {
158  Algorithm::Configure(config);
159  this->LoadConfig();
160 }
161 //____________________________________________________________________________
163 {
164  Algorithm::Configure(config);
165  this->LoadConfig();
166 }
167 //____________________________________________________________________________
169 {
170  //AlgConfigPool * confp = AlgConfigPool::Instance();
171  //const Registry * gc = confp->GlobalParameterList();
172 
173  /*fRo = fConfig->GetDoubleDef("COH-Ro", gc->GetDouble("COH-Ro"));
174  fMa = fConfig->GetDoubleDef("Ma", gc->GetDouble("COH-Ma"));
175  fa4 = fConfig->GetDoubleDef("a4", gc->GetDouble("COHAR-a4"));
176  fa5 = fConfig->GetDoubleDef("a5", gc->GetDouble("COHAR-a5"));
177  fb4 = fConfig->GetDoubleDef("b4", gc->GetDouble("COHAR-b4"));
178  fb5 = fConfig->GetDoubleDef("b5", gc->GetDouble("COHAR-b5"));
179  ffPi = fConfig->GetDoubleDef("fPi", gc->GetDouble("COHAR-fPi"));
180  ffStar = fConfig->GetDoubleDef("fStar", gc->GetDouble("COHAR-fStar"));*/
181 
182 
183  //-- load the differential cross section integrator
185  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
186  assert(fXSecIntegrator);
187 
188 }
189 //____________________________________________________________________________
Cross Section Calculation Interface.
bool IsWeak(void) const
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
Cross Section Integrator Interface.
int A(void) const
Definition: Target.h:70
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double Integral(const Interaction *i) const
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
bool IsCoherentProduction(void) const
enum genie::EKinePhaseSpace KinePhaseSpace_t
void Configure(const Registry &config)
Summary information for an interaction.
Definition: Interaction.h:56
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
bool IsWeakNC(void) const
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
double DXSec(const double E_nu_, const double E_l_, const double theta_l_, const double phi_l_, const double theta_pi_, const double phi_pi_)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
static constexpr double cm2
Definition: Units.h:69
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
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 HitNucIsSet(void) const
Definition: Target.cxx:283
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
const Target & Tgt(void) const
Definition: InitialState.h:66
const XSecIntegratorI * fXSecIntegrator
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
const Interaction * fLastInteraction
alvarezruso::AlvarezRusoCOHPiPDXSec * fMultidiff
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345