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