SlowRsclCharmDISPXSecLO.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 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
16 #include "Framework/Conventions/GBuild.h"
29 #include "Framework/Utils/Range1.h"
30 
31 using namespace genie;
32 using namespace genie::constants;
33 
34 //____________________________________________________________________________
36 XSecAlgorithmI("genie::SlowRsclCharmDISPXSecLO")
37 {
38 
39 }
40 //____________________________________________________________________________
42 XSecAlgorithmI("genie::SlowRsclCharmDISPXSecLO", config)
43 {
44 
45 }
46 //____________________________________________________________________________
48 {
49 
50 }
51 //____________________________________________________________________________
53  const Interaction * interaction, KinePhaseSpace_t kps) const
54 {
55  if(! this -> ValidProcess (interaction) ) return 0.;
56  if(! this -> ValidKinematics (interaction) ) return 0.;
57 
58  if(interaction->ProcInfo().IsWeakNC()) return 0;
59 
60  //----- get kinematics & init-state parameters
61  const Kinematics & kinematics = interaction->Kine();
62  const InitialState & init_state = interaction->InitState();
63  const Target & target = init_state.Tgt();
64 
65  double Mnuc = target.HitNucMass();
66  double E = init_state.ProbeE(kRfHitNucRest);
67  double x = kinematics.x();
68  double y = kinematics.y();
69  double Q2 = 2*Mnuc*E*x*y;
70 
71  //----- get target information (hit nucleon and quark)
72  int nu = init_state.ProbePdg();
73  int nuc = target.HitNucPdg();
74  bool isP = pdg::IsProton (nuc);
75  bool isN = pdg::IsNeutron(nuc);
76  bool qset = target.HitQrkIsSet();
77  int qpdg = (qset) ? target.HitQrkPdg() : 0;
78  bool sea = (qset) ? target.HitSeaQrk() : false;
79  bool isd = (qset) ? pdg::IsDQuark (qpdg) : false;
80  bool iss = (qset) ? pdg::IsSQuark (qpdg) : false;
81  bool isdb = (qset) ? pdg::IsAntiDQuark (qpdg) : false;
82  bool issb = (qset) ? pdg::IsAntiSQuark (qpdg) : false;
83  bool isnu = pdg::IsNeutrino(nu);
84  bool isnub = pdg::IsAntiNeutrino(nu);
85 
86  //----- compute slow rescaling variable & check its value
87  double xi = x * (1 + fMc2/Q2);
88  if(xi<=0 || xi>1) return 0.;
89 
90  //----- calculate the PDFs
91  PDF pdfs;
92  pdfs.SetModel(fPDFModel); // <-- attach algorithm
93  pdfs.Calculate(xi, Q2); // <-- calculate
94 
95  //----- proton pdfs
96  double us = pdfs.UpSea()/xi;
97  double uv = pdfs.UpValence()/xi;
98  double ds = pdfs.DownSea()/xi;
99  double dv = pdfs.DownValence()/xi;
100  double s = pdfs.Strange()/xi;
101 
102  //----- if the hit nucleon is a neutron, swap u<->d
103  double tmp;
104  if (isN) {
105  tmp = uv; uv = dv; dv = tmp;
106  tmp = us; us = ds; ds = tmp;
107  }
108 
109  //----- if a hit quark has been set then switch off other contributions
110  if(qset) {
111  if(isnub) { bool pass = (isdb||issb)&&sea; if(!pass) return 0; }
112  if(isnu) { bool pass = isd||(iss&&sea); if(!pass) return 0; }
113  dv = ( isd && !sea) ? dv : 0.;
114  ds = ( (isd||isdb) && sea) ? ds : 0.;
115  s = ( (iss||issb) && sea) ? s : 0.;
116  }
117  //----- in case of a \bar{v}+N calculation (quark not set) zero the d/val contribution
118  if(isnub) dv=0;
119 
120  //----- calculate cross section
121  double Gw2 = TMath::Power( (kGF/kSqrt2)*(1+Q2/kMw2), 2);
122  double xsec_0 = Gw2 * 2*Q2/(y*kPi) * (y + xi*(1-y)/x);
123  double xsec_d = xsec_0 * fVcd2 * (dv+ds);
124  double xsec_s = xsec_0 * fVcs2 * s;
125  double xsec = xsec_d + xsec_s;
126 
127 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
128  double Mnuc2 = TMath::Power(Mnuc, 2);
129  double W2 = Mnuc2 + 2*Mnuc*E*y*(1-x);
130  double W = TMath::Max(0., TMath::Sqrt(W2));
131  LOG("DISCharmXSec", pDEBUG)
132  << "\n dxsec[DISCharm,FreeN]/dxdy (E= " << E
133  << ", x= " << x << ", y= " << y
134  << ", W= " << W << ", Q2 = " << Q2 << ") = " << xsec;
135 #endif
136 
137  //----- The algorithm computes d^2xsec/dxdy
138  // Check whether variable tranformation is needed
139  if(kps!=kPSxyfE) {
140  double J = utils::kinematics::Jacobian(interaction,kPSxyfE,kps);
141  xsec *= J;
142  }
143 
144  //----- If requested return the free nucleon xsec even for input nuclear tgt
145  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
146 
147  //----- Nuclear cross section (simple scaling here)
148  int NNucl = (isP) ? target.Z() : target.N();
149  xsec *= NNucl;
150 
151  return xsec;
152 }
153 //____________________________________________________________________________
155 {
156  double xsec = fXSecIntegrator->Integrate(this,interaction);
157  return xsec;
158 }
159 //____________________________________________________________________________
161  const Interaction * interaction) const
162 {
163  if(interaction->TestBit(kISkipProcessChk)) return true;
164 
165  const XclsTag & xcls = interaction->ExclTag();
166  const InitialState & init_state = interaction->InitState();
167  const ProcessInfo & proc_info = interaction->ProcInfo();
168 
169  if(!proc_info.IsDeepInelastic()) return false;
170  if(!proc_info.IsWeak()) return false;
171 
172  bool is_inclusive_charm = (xcls.IsCharmEvent() && xcls.IsInclusiveCharm());
173  if(!is_inclusive_charm) return false;
174 
175  int nu = init_state.ProbePdg();
176  int nuc = init_state.Tgt().HitNucPdg();
177 
178  if (!pdg::IsProton(nuc) && !pdg::IsNeutron(nuc)) return false;
179  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
180 
181  return true;
182 }
183 //____________________________________________________________________________
185 {
186  Algorithm::Configure(config);
187  this->LoadConfig();
188 }
189 //____________________________________________________________________________
191 {
192  Algorithm::Configure(param_set);
193  this->LoadConfig();
194 }
195 //____________________________________________________________________________
197 {
198  // read mc, Vcd, Vcs from config or set defaults
199  GetParam( "Charm-Mass", fMc ) ;
200  GetParam( "CKM-Vcd", fVcd ) ;
201  GetParam( "CKM-Vcs", fVcs ) ;
202 
203  fMc2 = TMath::Power(fMc, 2);
204  fVcd2 = TMath::Power(fVcd, 2);
205  fVcs2 = TMath::Power(fVcs, 2);
206 
207  // load PDF set
208  fPDFModel = dynamic_cast<const PDFModelI *> (this->SubAlg("PDF-Set"));
209  assert(fPDFModel);
210 
211  // load XSec Integrator
213  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
214  assert(fXSecIntegrator);
215 }
216 //____________________________________________________________________________
Cross Section Calculation Interface.
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Basic constants.
bool IsWeak(void) const
bool HitSeaQrk(void) const
Definition: Target.cxx:299
static const double kSqrt2
Definition: Constants.h:115
static constexpr double us
Definition: Units.h:97
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
static const double kMw2
Definition: Constants.h:93
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
int HitNucPdg(void) const
Definition: Target.cxx:304
double DownValence(void) const
Definition: PDF.h:51
int HitQrkPdg(void) const
Definition: Target.cxx:242
double HitNucMass(void) const
Definition: Target.cxx:233
A class to store PDFs.
Definition: PDF.h:37
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double x(bool selected=false) const
Definition: Kinematics.cxx:99
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:273
double UpSea(void) const
Definition: PDF.h:52
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:303
enum genie::EKinePhaseSpace KinePhaseSpace_t
Pure abstract base class. Defines the PDFModelI interface to be implemented by wrapper classes to exi...
Definition: PDFModelI.h:28
bool IsAntiDQuark(int pdgc)
Definition: PDGUtils.cxx:298
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
double y(bool selected=false) const
Definition: Kinematics.cxx:112
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
double DownSea(void) const
Definition: PDF.h:53
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
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 IsProton(int pdgc)
Definition: PDGUtils.cxx:333
bool IsWeakNC(void) const
double Strange(void) const
Definition: PDF.h:54
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
static Config * config
Definition: config.cpp:1054
HLTPathStatus const pass
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:115
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
void SetModel(const PDFModelI *model)
Definition: PDF.cxx:42
int Z(void) const
Definition: Target.h:68
static const double kGF
Definition: Constants.h:58
string tmp
Definition: languages.py:63
void Calculate(double x, double q2)
Definition: PDF.cxx:49
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:84
int N(void) const
Definition: Target.h:69
double UpValence(void) const
Definition: PDF.h:50
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
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
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:54
void Configure(const Registry &config)
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:268
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
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
list x
Definition: train.py:276
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
double Integral(const Interaction *i) const
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:37
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
const XSecIntegratorI * fXSecIntegrator
static QCString * s
Definition: config.cpp:1042
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