GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PaisQELLambdaPXSec.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  Hugh Gallagher
7  Tufts University
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
14 #include "Framework/Conventions/GBuild.h"
31 
32 using namespace genie;
33 using namespace genie::constants;
34 using namespace genie::utils;
35 
36 //____________________________________________________________________________
38 XSecAlgorithmI("genie::PaisQELLambdaPXSec")
39 {
40 
41 }
42 //____________________________________________________________________________
44 XSecAlgorithmI("genie::PaisQELLambdaPXSec", config)
45 {
46 
47 }
48 //____________________________________________________________________________
50 {
51 
52 }
53 //____________________________________________________________________________
55  const Interaction * interaction, KinePhaseSpace_t kps) const
56 {
57  if(! this -> ValidProcess (interaction) ) return 0.;
58  if(! this -> ValidKinematics (interaction) ) return 0.;
59 
60  //----- get kinematics & init state - compute auxiliary vars
61  const Kinematics & kinematics = interaction->Kine();
62  const InitialState & init_state = interaction->InitState();
63  const Target & target = init_state.Tgt();
64 
65  //neutrino energy & momentum transfer
66  double E = init_state.ProbeE(kRfHitNucRest);
67  double E2 = E * E;
68  double q2 = kinematics.q2();
69 
70 
71  //resonance mass & nucleon mass
72  double Mnuc = target.HitNucMass();
73  double Mnuc2 = TMath::Power(Mnuc,2);
74 
75  //----- Calculate the differential cross section dxsec/dQ^2
76  double Gf = kGF2 / (2*kPi);
77  double ml = interaction->FSPrimLepton()->Mass();
78  double ml2 = TMath::Power(ml,2);
79  double M1 = Mnuc;
80  double M2 = (this)->MHyperon(interaction);
81  double v = (TMath::Power(M2,2) - Mnuc2 - q2) / (2*Mnuc);
82  double v2 = TMath::Power(v,2);
83  double s = Mnuc2 + 2*Mnuc*E;
84  double u = Mnuc2 + ml2 + 2*v*Mnuc - 2*Mnuc*E;
85 
86 // xsec term changes sign for antineutrinos
87  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
88  int sign = (is_neutrino) ? -1 : 1;
89 
90 // Calculate the QEL form factors
91  fFormFactors.Calculate(interaction);
92 
93  double F1V = fFormFactors.F1V();
94  double xiF2V = fFormFactors.xiF2V();
95  double FA = fFormFactors.FA();
96 // double Fp = fFormFactors.Fp();
97 
98 // calculate w coefficients
99  //start with Mass terms
100  double Mp = M2 + M1;
101  double Mm = M2 - M1;
102  double Mm2 = TMath::Power(Mm, 2);
103  double Mp2 = TMath::Power(Mp, 2);
104 
105  //Powers of Form Factors
106  double FA2 = TMath::Power(FA, 2);
107 // double FA3 = 0;
108 
109  //Calculate W terms
110 
111  double w1 = (Mm2 - q2)/(4*Mnuc2)*TMath::Power((F1V + xiF2V), 2) + (Mp2 - q2)/(4*Mnuc2) * FA2;
112  double w2 = FA2 + TMath::Power((F1V + xiF2V - Mp * xiF2V / (2 * Mnuc)), 2) - q2 / Mnuc2 * TMath::Power((xiF2V / 2), 2);
113  double w3 = 2 * FA * (F1V + xiF2V);
114 
115  double xsec = Gf*fSin8c2 / (16*Mnuc2*E2) * (-8*Mnuc2*q2*w1 - 4*(Mnuc2*v2 - q2)*w2 - sign*2*(s - u)*q2*w3 + (s-u)*(s-u)*w2);
116  xsec = TMath::Max(xsec,0.);
117 
118  //----- The algorithm computes dxsec/dQ2
119  // Check whether variable tranformation is needed
120  if(kps!=kPSQ2fE) {
121  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
122  xsec *= J;
123  }
124 
125  //----- If requested return the free nucleon xsec even for input nuclear tgt
126  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
127 
128  //----- Nuclear cross section (simple scaling here)
129  int nuc = target.HitNucPdg();
130  int NNucl = (pdg::IsProton(nuc)) ? target.Z() : target.N();
131  xsec *= NNucl;
132 
133  return xsec;
134 }
135 //____________________________________________________________________________
136 double PaisQELLambdaPXSec::MHyperon(const Interaction * interaction) const
137 {
138  const XclsTag & xcls = interaction->ExclTag();
139 
140  int pdgc = xcls.StrangeHadronPdg();
141  double MR = PDGLibrary::Instance()->Find(pdgc)->Mass();
142  return MR;
143 }
144 //____________________________________________________________________________
145 double PaisQELLambdaPXSec::Integral(const Interaction * interaction) const
146 {
147 
148  double xsec = fXSecIntegrator->Integrate(this,interaction);
149  return xsec;
150 }
151 //____________________________________________________________________________
153  const Interaction * interaction) const
154 {
155  // Make sure we are dealing with one of the following channels:
156  // v + n --> mu+ + Sigma^{-}
157  // v + p --> mu+ + Lambda^{0}
158  // v + p --> mu+ + Sigma^{0}
159 
160  if(interaction->TestBit(kISkipProcessChk)) return true;
161 
162  const XclsTag & xcls = interaction->ExclTag();
163  const InitialState & init_state = interaction->InitState();
164  const ProcessInfo & proc_info = interaction->ProcInfo();
165 
166  bool is_exclusive_strange = (xcls.IsStrangeEvent() && !xcls.IsInclusiveStrange());
167  if(!is_exclusive_strange) return false;
168 
169  if(!proc_info.IsQuasiElastic()) return false;
170  if(!proc_info.IsWeak()) return false;
171 
172  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
173  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );
174 
175  int pdgc = xcls.StrangeHadronPdg();
176 
177  bool can_handle = (
178  (pdgc == kPdgSigmaM && isN) || /* v + n -> l + #Sigma^{-} */
179  (pdgc == kPdgLambda && isP) || /* v + p -> l + #Lambda^{0} */
180  (pdgc == kPdgSigma0 && isP) /* v + p -> l + #Sigma^{0} */
181  );
182 
183  return can_handle;
184 }
185 //____________________________________________________________________________
187  const Interaction * interaction) const
188 {
189  if(interaction->TestBit(kISkipKinematicChk)) return true;
190 
191  const InitialState & init_state = interaction->InitState();
192  double E = init_state.ProbeE(kRfHitNucRest);
193 
194  //resonance, final state primary lepton & nucleon mass
195  double MR = this -> MHyperon (interaction);
196  double ml = interaction->FSPrimLepton()->Mass();
197  double Mnuc = init_state.Tgt().HitNucP4Ptr()->M();
198  double Mnuc2 = TMath::Power(Mnuc,2);
199 
200  //resonance threshold
201  double ER = ( TMath::Power(MR+ml,2) - Mnuc2 ) / (2*Mnuc);
202 
203  if(E <= ER) return false;
204 
205  return true;
206 }
207 //____________________________________________________________________________
209 {
210  Algorithm::Configure(config);
211  this->LoadConfig();
212 }
213 //____________________________________________________________________________
214 void PaisQELLambdaPXSec::Configure(string param_set)
215 {
216  Algorithm::Configure(param_set);
217  this->LoadConfig();
218 }
219 //____________________________________________________________________________
221 {
222 
223  double thc ;
224  GetParam( "CabibboAngle", thc ) ;
225  fSin8c2 = TMath::Power(TMath::Sin(thc), 2);
226 
227  // load QEL form factors model
228  fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
229  this->SubAlg("FormFactorsAlg"));
230  assert(fFormFactorsModel);
231  fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
232 
233  // load XSec Integrator
235  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
236  assert(fXSecIntegrator);
237 }
238 //____________________________________________________________________________
Cross Section Calculation Interface.
bool IsWeak(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
Cross Section Integrator Interface.
const int kPdgLambda
Definition: PDGCodes.h:85
void Configure(const Registry &config)
int HitNucPdg(void) const
Definition: Target.cxx:304
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
double HitNucMass(void) const
Definition: Target.cxx:233
bool IsStrangeEvent(void) const
Definition: XclsTag.h:53
static constexpr double s
Definition: Units.h:95
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
const int kPdgSigma0
Definition: PDGCodes.h:88
enum genie::EKinePhaseSpace KinePhaseSpace_t
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
Summary information for an interaction.
Definition: Interaction.h:56
double MHyperon(const Interaction *interaction) const
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double q2(bool selected=false) const
Definition: Kinematics.cxx:141
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
int Z(void) const
Definition: Target.h:68
double xiF2V(void) const
Get the computed form factor xi*F2V.
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
const int kPdgSigmaM
Definition: PDGCodes.h:89
double Integral(const Interaction *i) const
bool IsInclusiveStrange(void) const
Definition: XclsTag.cxx:71
int N(void) const
Definition: Target.h:69
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
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
double F1V(void) const
Get the computed form factor F1V.
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
double FA(void) const
Get the computed form factor FA.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
const XSecIntegratorI * fXSecIntegrator
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
const QELFormFactorsModelI * fFormFactorsModel
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345