GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AivazisCharmPXSecLO.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
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::AivazisCharmPXSecLO")
37 {
38 
39 }
40 //____________________________________________________________________________
42 XSecAlgorithmI("genie::AivazisCharmPXSecLO", 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() || interaction->ProcInfo().IsDarkMatter()) return 0;
59 
60  //----- get init-state & kinematical parameters
61  const Kinematics & kinematics = interaction -> Kine();
62  const InitialState & init_state = interaction -> InitState();
63  const Target & target = init_state.Tgt();
64 
65  //----- get target information (hit nucleon and quark)
66  int nu = init_state.ProbePdg();
67  int nuc = target.HitNucPdg();
68  bool isP = pdg::IsProton (nuc);
69  bool isN = pdg::IsNeutron(nuc);
70  bool qset = target.HitQrkIsSet();
71  int qpdg = (qset) ? target.HitQrkPdg() : 0;
72  bool sea = (qset) ? target.HitSeaQrk() : false;
73  bool isd = (qset) ? pdg::IsDQuark (qpdg) : false;
74  bool iss = (qset) ? pdg::IsSQuark (qpdg) : false;
75  bool isdb = (qset) ? pdg::IsAntiDQuark (qpdg) : false;
76  bool issb = (qset) ? pdg::IsAntiSQuark (qpdg) : false;
77  bool isnu = pdg::IsNeutrino(nu);
78  bool isnub = pdg::IsAntiNeutrino(nu);
79 
80  //----- compute kinematic & auxiliary parameters
81  double E = init_state.ProbeE(kRfHitNucRest);
82  double x = kinematics.x();
83  double y = kinematics.y();
84  double x2 = TMath::Power(x, 2);
85  double Mnuc = target.HitNucMass();
86  double Mnuc2 = TMath::Power(Mnuc, 2);
87  double Q2 = 2*Mnuc*E*x*y;
88  double inverse_eta = 0.5/x + TMath::Sqrt( 0.25/x2 + Mnuc2/Q2 );
89  double eta = 1 / inverse_eta;
90  double xi = eta * (1 + fMc2/Q2);
91  double coshpsi = (2-y)/y; // hyperbolic-cosine(psi)
92  double sinh2psi = TMath::Power(coshpsi, 2) - 1;
93 
94  //----- make sure that the mass-corrected x is in physical region
95  if(xi<=0 || xi>1) return 0;
96 
97  //----- Calculate the PDFs
98  PDF pdfs;
99  pdfs.SetModel(fPDFModel); // <-- attach algorithm
100  pdfs.Calculate(xi, Q2); // <-- calculate
101 
102  //----- proton pdfs
103  double us = pdfs.UpSea()/xi;
104  double uv = pdfs.UpValence()/xi;
105  double ds = pdfs.DownSea()/xi;
106  double dv = pdfs.DownValence()/xi;
107  double s = pdfs.Strange()/xi;
108 
109  //----- if the hit nucleon is a neutron, swap u<->d
110  double tmp;
111  if (isN) {
112  tmp = uv; uv = dv; dv = tmp;
113  tmp = us; us = ds; ds = tmp;
114  }
115 
116  //----- if a hit quark has been set then switch off other contributions
117  if(qset) {
118  if(isnub) { bool pass = (isdb||issb)&&sea; if(!pass) return 0; }
119  if(isnu) { bool pass = isd||(iss&&sea); if(!pass) return 0; }
120  dv = ( isd && !sea) ? dv : 0.;
121  ds = ( (isd||isdb) && sea) ? ds : 0.;
122  s = ( (iss||issb) && sea) ? s : 0.;
123  }
124  //----- in case of a \bar{v}+N calculation (quark not set) zero the d/val contribution
125  if(isnub) dv=0;
126 
127  //----- calculate the cross section
128  double Gw2 = TMath::Power((kGF/kSqrt2)*(1+Q2/kMw2), 2);
129  double f1 = TMath::Power( (1+coshpsi)/2, 2);
130  double f2 = 0.25 * (fMc2/Q2) * sinh2psi;
131  double xsec_0 = 2 * Gw2 * (y*Q2/kPi) * (f1+f2);
132  double xsec_d = xsec_0 * fVcd2 * (dv+ds);
133  double xsec_s = xsec_0 * fVcs2 * s;
134  double xsec = xsec_d + xsec_s;
135 
136 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
137  double W2 = Mnuc2 + 2*Mnuc*E*y*(1-x);
138  double W = TMath::Max(0., TMath::Sqrt(W2));
139  LOG("DISCharmXSec", pDEBUG)
140  << "\n dxsec[DISCharm,FreeN]/dxdy (E= " << E
141  << ", x= " << x << ", y= " << y
142  << ", W= " << W << ", Q2 = " << Q2 << ") = " << xsec;
143 #endif
144 
145  //----- The algorithm computes d^2xsec/dxdy
146  // Check whether variable tranformation is needed
147  if(kps!=kPSxyfE) {
148  double J = utils::kinematics::Jacobian(interaction,kPSxyfE,kps);
149  xsec *= J;
150  }
151 
152  //----- If requested return the free nucleon xsec even for input nuclear tgt
153  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
154 
155  //----- Nuclear cross section (simple scaling here)
156  int NNucl = (isP) ? target.Z() : target.N();
157  xsec *= NNucl;
158 
159  return xsec;
160 }
161 //____________________________________________________________________________
162 double AivazisCharmPXSecLO::Integral(const Interaction * interaction) const
163 {
164  double xsec = fXSecIntegrator->Integrate(this,interaction);
165  return xsec;
166 }
167 //____________________________________________________________________________
168 bool AivazisCharmPXSecLO::ValidProcess(const Interaction * interaction) const
169 {
170  if(interaction->TestBit(kISkipProcessChk)) return true;
171 
172  const XclsTag & xcls = interaction->ExclTag();
173  const InitialState & init_state = interaction->InitState();
174  const ProcessInfo & proc_info = interaction->ProcInfo();
175 
176  if(!proc_info.IsDeepInelastic()) return false;
177  if(!proc_info.IsWeak()) return false;
178 
179  bool is_inclusive_charm = (xcls.IsCharmEvent() && xcls.IsInclusiveCharm());
180  if(!is_inclusive_charm) return false;
181 
182  int nu = init_state.ProbePdg();
183  int nuc = init_state.Tgt().HitNucPdg();
184 
185  if (!pdg::IsProton(nuc) && !pdg::IsNeutron(nuc)) return false;
186  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
187 
188  return true;
189 }
190 //____________________________________________________________________________
192 {
193  Algorithm::Configure(config);
194  this->LoadConfig();
195 }
196 //____________________________________________________________________________
197 void AivazisCharmPXSecLO::Configure(string param_set)
198 {
199  Algorithm::Configure(param_set);
200  this->LoadConfig();
201 }
202 //____________________________________________________________________________
204 {
205 
206  // read mc, Vcd, Vcs from config or set defaults
207  GetParam( "Charm-Mass", fMc ) ;
208  GetParam( "CKM-Vcd", fVcd ) ;
209  GetParam( "CKM-Vcs", fVcs ) ;
210 
211  fMc2 = TMath::Power(fMc, 2);
212  fVcd2 = TMath::Power(fVcd, 2);
213  fVcs2 = TMath::Power(fVcs, 2);
214 
215  // load PDF set
216  fPDFModel = dynamic_cast<const PDFModelI *> (this->SubAlg("PDF-Set"));
217  assert(fPDFModel);
218 
219  // load XSec Integrator
221  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
222  assert(fXSecIntegrator);
223 }
224 //____________________________________________________________________________
Cross Section Calculation Interface.
bool IsWeak(void) const
bool HitSeaQrk(void) const
Definition: Target.cxx:299
static constexpr double us
Definition: Units.h:97
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.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
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
static constexpr double s
Definition: Units.h:95
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double x(bool selected=false) const
Definition: Kinematics.cxx:99
bool IsSQuark(int pdgc)
Definition: PDGUtils.cxx:276
double UpSea(void) const
Definition: PDF.h:52
bool IsAntiSQuark(int pdgc)
Definition: PDGUtils.cxx:306
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:301
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 W(const Interaction *const i)
Definition: KineUtils.cxx:1101
double DownSea(void) const
Definition: PDF.h:53
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
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:336
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
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:118
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
void Calculate(double x, double q2)
Definition: PDF.cxx:49
bool IsDeepInelastic(void) const
Definition: ProcessInfo.cxx:89
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:69
double UpValence(void) const
Definition: PDF.h:50
bool HitQrkIsSet(void) const
Definition: Target.cxx:292
bool IsDarkMatter(void) const
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 Integral(const Interaction *i) const
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:54
bool IsDQuark(int pdgc)
Definition: PDGUtils.cxx:271
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
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
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:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345