PaisQELLambdaPXSec.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: Hugh Gallagher
8  Tufts University
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14 */
15 //____________________________________________________________________________
16 
17 #include <TMath.h>
18 
37 
38 using namespace genie;
39 using namespace genie::constants;
40 using namespace genie::utils;
41 
42 //____________________________________________________________________________
44 XSecAlgorithmI("genie::PaisQELLambdaPXSec")
45 {
46 
47 }
48 //____________________________________________________________________________
50 XSecAlgorithmI("genie::PaisQELLambdaPXSec", config)
51 {
52 
53 }
54 //____________________________________________________________________________
56 {
57 
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  //----- get kinematics & init state - compute auxiliary vars
67  const Kinematics & kinematics = interaction->Kine();
68  const InitialState & init_state = interaction->InitState();
69  const Target & target = init_state.Tgt();
70 
71  //neutrino energy & momentum transfer
72  double E = init_state.ProbeE(kRfHitNucRest);
73  double E2 = E * E;
74  double q2 = kinematics.q2();
75 
76 
77  //resonance mass & nucleon mass
78  double Mnuc = target.HitNucMass();
79  double Mnuc2 = TMath::Power(Mnuc,2);
80 
81  //----- Calculate the differential cross section dxsec/dQ^2
82  double Gf = kGF2 / (2*kPi);
83  double ml = interaction->FSPrimLepton()->Mass();
84  double ml2 = TMath::Power(ml,2);
85  double M1 = Mnuc;
86  double M2 = (this)->MHyperon(interaction);
87  double v = (TMath::Power(M2,2) - Mnuc2 - q2) / (2*Mnuc);
88  double v2 = TMath::Power(v,2);
89  double s = Mnuc2 + 2*Mnuc*E;
90  double u = Mnuc2 + ml2 + 2*v*Mnuc - 2*Mnuc*E;
91 
92 // xsec term changes sign for antineutrinos
93  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
94  int sign = (is_neutrino) ? -1 : 1;
95 
96 // Calculate the QEL form factors
97  fFormFactors.Calculate(interaction);
98 
99  double F1V = fFormFactors.F1V();
100  double xiF2V = fFormFactors.xiF2V();
101  double FA = fFormFactors.FA();
102 // double Fp = fFormFactors.Fp();
103 
104 // calculate w coefficients
105  //start with Mass terms
106  double Mp = M2 + M1;
107  double Mm = M2 - M1;
108  double Mm2 = TMath::Power(Mm, 2);
109  double Mp2 = TMath::Power(Mp, 2);
110 
111  //Powers of Form Factors
112  double FA2 = TMath::Power(FA, 2);
113 // double FA3 = 0;
114 
115  //Calculate W terms
116 
117  double w1 = (Mm2 - q2)/(4*Mnuc2)*TMath::Power((F1V + xiF2V), 2) + (Mp2 - q2)/(4*Mnuc2) * FA2;
118  double w2 = FA2 + TMath::Power((F1V + xiF2V - Mp * xiF2V / (2 * Mnuc)), 2) - q2 / Mnuc2 * TMath::Power((xiF2V / 2), 2);
119  double w3 = 2 * FA * (F1V + xiF2V);
120 
121  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);
122  xsec = TMath::Max(xsec,0.);
123 
124  //----- The algorithm computes dxsec/dQ2
125  // Check whether variable tranformation is needed
126  if(kps!=kPSQ2fE) {
127  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
128  xsec *= J;
129  }
130 
131  //----- If requested return the free nucleon xsec even for input nuclear tgt
132  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
133 
134  //----- Nuclear cross section (simple scaling here)
135  int nuc = target.HitNucPdg();
136  int NNucl = (pdg::IsProton(nuc)) ? target.Z() : target.N();
137  xsec *= NNucl;
138 
139  return xsec;
140 }
141 //____________________________________________________________________________
142 double PaisQELLambdaPXSec::MHyperon(const Interaction * interaction) const
143 {
144  const XclsTag & xcls = interaction->ExclTag();
145 
146  int pdgc = xcls.StrangeHadronPdg();
147  double MR = PDGLibrary::Instance()->Find(pdgc)->Mass();
148  return MR;
149 }
150 //____________________________________________________________________________
151 double PaisQELLambdaPXSec::Integral(const Interaction * interaction) const
152 {
153 
154  double xsec = fXSecIntegrator->Integrate(this,interaction);
155  return xsec;
156 }
157 //____________________________________________________________________________
159  const Interaction * interaction) const
160 {
161  // Make sure we are dealing with one of the following channels:
162  // v + n --> mu+ + Sigma^{-}
163  // v + p --> mu+ + Lambda^{0}
164  // v + p --> mu+ + Sigma^{0}
165 
166  if(interaction->TestBit(kISkipProcessChk)) return true;
167 
168  const XclsTag & xcls = interaction->ExclTag();
169  const InitialState & init_state = interaction->InitState();
170  const ProcessInfo & proc_info = interaction->ProcInfo();
171 
172  bool is_exclusive_strange = (xcls.IsStrangeEvent() && !xcls.IsInclusiveStrange());
173  if(!is_exclusive_strange) return false;
174 
175  if(!proc_info.IsQuasiElastic()) return false;
176  if(!proc_info.IsWeak()) return false;
177 
178  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
179  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );
180 
181  int pdgc = xcls.StrangeHadronPdg();
182 
183  bool can_handle = (
184  (pdgc == kPdgSigmaM && isN) || /* v + n -> l + #Sigma^{-} */
185  (pdgc == kPdgLambda && isP) || /* v + p -> l + #Lambda^{0} */
186  (pdgc == kPdgSigma0 && isP) /* v + p -> l + #Sigma^{0} */
187  );
188 
189  return can_handle;
190 }
191 //____________________________________________________________________________
193  const Interaction * interaction) const
194 {
195  if(interaction->TestBit(kISkipKinematicChk)) return true;
196 
197  const InitialState & init_state = interaction->InitState();
198  double E = init_state.ProbeE(kRfHitNucRest);
199 
200  //resonance, final state primary lepton & nucleon mass
201  double MR = this -> MHyperon (interaction);
202  double ml = interaction->FSPrimLepton()->Mass();
203  double Mnuc = init_state.Tgt().HitNucP4Ptr()->M();
204  double Mnuc2 = TMath::Power(Mnuc,2);
205 
206  //resonance threshold
207  double ER = ( TMath::Power(MR+ml,2) - Mnuc2 ) / (2*Mnuc);
208 
209  if(E <= ER) return false;
210 
211  return true;
212 }
213 //____________________________________________________________________________
215 {
216  Algorithm::Configure(config);
217  this->LoadConfig();
218 }
219 //____________________________________________________________________________
220 void PaisQELLambdaPXSec::Configure(string param_set)
221 {
222  Algorithm::Configure(param_set);
223  this->LoadConfig();
224 }
225 //____________________________________________________________________________
227 {
228 
229  double thc ;
230  GetParam( "CabibboAngle", thc ) ;
231  fSin8c2 = TMath::Power(TMath::Sin(thc), 2);
232 
233  // load QEL form factors model
234  fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
235  this->SubAlg("FormFactorsAlg"));
237  fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
238 
239  // load XSec Integrator
241  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
243 }
244 //____________________________________________________________________________
Cross Section Calculation Interface.
const double kPi
Basic constants.
bool IsWeak(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
const XML_Char * target
Definition: expat.h:268
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
const int kPdgLambda
Definition: PDGCodes.h:69
void Configure(const Registry &config)
int HitNucPdg(void) const
Definition: Target.cxx:321
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
double HitNucMass(void) const
Definition: Target.cxx:250
bool IsStrangeEvent(void) const
Definition: XclsTag.h:51
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
const int kPdgSigma0
Definition: PDGCodes.h:72
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:37
Double_t q2[12][num]
Definition: f2_nu.C:137
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
const XML_Char * s
Definition: expat.h:262
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:151
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
int StrangeHadronPdg(void) const
Definition: XclsTag.h:53
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
Float_t E
Definition: plot.C:20
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
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
double xiF2V(void) const
Get the computed form factor xi*F2V.
Double_t xsec[nknots]
Definition: testXsec.C:47
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:73
double Integral(const Interaction *i) const
bool IsInclusiveStrange(void) const
Definition: XclsTag.cxx:80
int N(void) const
Definition: Target.h:70
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
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:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
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.
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
static const double kGF2
Definition: Constants.h:60
void kinematics()
Definition: kinematics.C:10
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?
Root of GENIE utility namespaces.
def sign(x)
Definition: canMan.py:197
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:49
const QELFormFactorsModelI * fFormFactorsModel
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353