AivazisCharmPXSecLO.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::AivazisCharmPXSecLO")
43 {
44 
45 }
46 //____________________________________________________________________________
48 XSecAlgorithmI("genie::AivazisCharmPXSecLO", 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() || interaction->ProcInfo().IsDarkMatter()) return 0;
65 
66  //----- get init-state & kinematical parameters
67  const Kinematics & kinematics = interaction -> Kine();
68  const InitialState & init_state = interaction -> InitState();
69  const Target & target = init_state.Tgt();
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 kinematic & auxiliary parameters
87  double E = init_state.ProbeE(kRfHitNucRest);
88  double x = kinematics.x();
89  double y = kinematics.y();
90  double x2 = TMath::Power(x, 2);
91  double Mnuc = target.HitNucMass();
92  double Mnuc2 = TMath::Power(Mnuc, 2);
93  double Q2 = 2*Mnuc*E*x*y;
94  double inverse_eta = 0.5/x + TMath::Sqrt( 0.25/x2 + Mnuc2/Q2 );
95  double eta = 1 / inverse_eta;
96  double xi = eta * (1 + fMc2/Q2);
97  double coshpsi = (2-y)/y; // hyperbolic-cosine(psi)
98  double sinh2psi = TMath::Power(coshpsi, 2) - 1;
99 
100  //----- make sure that the mass-corrected x is in physical region
101  if(xi<=0 || xi>1) return 0;
102 
103  //----- Calculate the PDFs
104  PDF pdfs;
105  pdfs.SetModel(fPDFModel); // <-- attach algorithm
106  pdfs.Calculate(xi, Q2); // <-- calculate
107 
108  //----- proton pdfs
109  double us = pdfs.UpSea()/xi;
110  double uv = pdfs.UpValence()/xi;
111  double ds = pdfs.DownSea()/xi;
112  double dv = pdfs.DownValence()/xi;
113  double s = pdfs.Strange()/xi;
114 
115  //----- if the hit nucleon is a neutron, swap u<->d
116  double tmp;
117  if (isN) {
118  tmp = uv; uv = dv; dv = tmp;
119  tmp = us; us = ds; ds = tmp;
120  }
121 
122  //----- if a hit quark has been set then switch off other contributions
123  if(qset) {
124  if(isnub) { bool pass = (isdb||issb)&&sea; if(!pass) return 0; }
125  if(isnu) { bool pass = isd||(iss&&sea); if(!pass) return 0; }
126  dv = ( isd && !sea) ? dv : 0.;
127  ds = ( (isd||isdb) && sea) ? ds : 0.;
128  s = ( (iss||issb) && sea) ? s : 0.;
129  }
130  //----- in case of a \bar{v}+N calculation (quark not set) zero the d/val contribution
131  if(isnub) dv=0;
132 
133  //----- calculate the cross section
134  double Gw2 = TMath::Power((kGF/kSqrt2)*(1+Q2/kMw2), 2);
135  double f1 = TMath::Power( (1+coshpsi)/2, 2);
136  double f2 = 0.25 * (fMc2/Q2) * sinh2psi;
137  double xsec_0 = 2 * Gw2 * (y*Q2/kPi) * (f1+f2);
138  double xsec_d = xsec_0 * fVcd2 * (dv+ds);
139  double xsec_s = xsec_0 * fVcs2 * s;
140  double xsec = xsec_d + xsec_s;
141 
142 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
143  double W2 = Mnuc2 + 2*Mnuc*E*y*(1-x);
144  double W = TMath::Max(0., TMath::Sqrt(W2));
145  LOG("DISCharmXSec", pDEBUG)
146  << "\n dxsec[DISCharm,FreeN]/dxdy (E= " << E
147  << ", x= " << x << ", y= " << y
148  << ", W= " << W << ", Q2 = " << Q2 << ") = " << xsec;
149 #endif
150 
151  //----- The algorithm computes d^2xsec/dxdy
152  // Check whether variable tranformation is needed
153  if(kps!=kPSxyfE) {
154  double J = utils::kinematics::Jacobian(interaction,kPSxyfE,kps);
155  xsec *= J;
156  }
157 
158  //----- If requested return the free nucleon xsec even for input nuclear tgt
159  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
160 
161  //----- Nuclear cross section (simple scaling here)
162  int NNucl = (isP) ? target.Z() : target.N();
163  xsec *= NNucl;
164 
165  return xsec;
166 }
167 //____________________________________________________________________________
168 double AivazisCharmPXSecLO::Integral(const Interaction * interaction) const
169 {
170  double xsec = fXSecIntegrator->Integrate(this,interaction);
171  return xsec;
172 }
173 //____________________________________________________________________________
174 bool AivazisCharmPXSecLO::ValidProcess(const Interaction * interaction) const
175 {
176  if(interaction->TestBit(kISkipProcessChk)) return true;
177 
178  const XclsTag & xcls = interaction->ExclTag();
179  const InitialState & init_state = interaction->InitState();
180  const ProcessInfo & proc_info = interaction->ProcInfo();
181 
182  if(!proc_info.IsDeepInelastic()) return false;
183  if(!proc_info.IsWeak()) return false;
184 
185  bool is_inclusive_charm = (xcls.IsCharmEvent() && xcls.IsInclusiveCharm());
186  if(!is_inclusive_charm) return false;
187 
188  int nu = init_state.ProbePdg();
189  int nuc = init_state.Tgt().HitNucPdg();
190 
191  if (!pdg::IsProton(nuc) && !pdg::IsNeutron(nuc)) return false;
192  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
193 
194  return true;
195 }
196 //____________________________________________________________________________
198 {
199  Algorithm::Configure(config);
200  this->LoadConfig();
201 }
202 //____________________________________________________________________________
203 void AivazisCharmPXSecLO::Configure(string param_set)
204 {
205  Algorithm::Configure(param_set);
206  this->LoadConfig();
207 }
208 //____________________________________________________________________________
210 {
211 
212  // read mc, Vcd, Vcs from config or set defaults
213  GetParam( "Charm-Mass", fMc ) ;
214  GetParam( "CKM-Vcd", fVcd ) ;
215  GetParam( "CKM-Vcs", fVcs ) ;
216 
217  fMc2 = TMath::Power(fMc, 2);
218  fVcd2 = TMath::Power(fVcd, 2);
219  fVcs2 = TMath::Power(fVcs, 2);
220 
221  // load PDF set
222  fPDFModel = dynamic_cast<const PDFModelI *> (this->SubAlg("PDF-Set"));
223  assert(fPDFModel);
224 
225  // load XSec Integrator
227  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
229 }
230 //____________________________________________________________________________
Cross Section Calculation Interface.
const double kPi
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
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
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
Float_t f2
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?
const XSecIntegratorI * fXSecIntegrator
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
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
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
Float_t f1
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
void Configure(const Registry &config)
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
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
bool IsDarkMatter(void) const
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
double Integral(const Interaction *i) const
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:63
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
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 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