GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BertuzzoDNuCOHPXSec.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  Author: Iker de Icaza <i.de-icaza-astiz \at sussex.ac.uk>
7  University of Sussex
8 
9  Costas Andreopoulos <c.andreopoulos \at cern.ch>
10  University of Liverpool
11 */
12 //____________________________________________________________________________
13 
14 
15 #include <TMath.h>
16 #include <Math/Integrator.h>
17 
28 
29 using namespace genie;
30 using namespace genie::utils;
31 using namespace genie::constants;
32 
33 //____________________________________________________________________________
35 XSecAlgorithmI("genie::BertuzzoDNuCOHPXSec")
36 {
37 
38 }
39 //____________________________________________________________________________
41 XSecAlgorithmI("genie::BertuzzoDNuCOHPXSec", config)
42 {
43 
44 }
45 //____________________________________________________________________________
47 {
48 
49 }
50 //____________________________________________________________________________
52  const Interaction * interaction, KinePhaseSpace_t kps) const
53 {
54  if(! this -> ValidProcess (interaction) ) return 0.;
55  if(! this -> ValidKinematics (interaction) ) return 0.;
56 
57  const InitialState & init_state = interaction -> InitState();
58  const Kinematics & kinematics = interaction -> Kine();
59  const Target & target = init_state.Tgt();
60 
61  // User inputs to the calculation
62  const int nu_pdg = init_state.ProbePdg();
63  const double E = init_state.ProbeE(kRfLab); // neutrino energy, units: GeV
64  const double Q2 = kinematics.Q2(); // momentum transfer, units: GeV^2
65  const double TE = kinematics.HadSystP4().E(); // energy of the target
66  const unsigned int Z = target.Z(); // number of protons
67 
68  // select the mixing depending on the incoming neutrino
69  unsigned short nu_i = 3;
70  if( pdg::IsNuE( TMath::Abs( nu_pdg ) ) ) nu_i = 0;
71  else if ( pdg::IsNuMu( TMath::Abs( nu_pdg ) ) ) nu_i = 1;
72  else if ( pdg::IsNuTau( TMath::Abs( nu_pdg ) ) ) nu_i = 2;
73  const double DTheta2 = fMixing2s[nu_i] * fMixing2s[3] ;
74 
75  // Target atomic mass number and mass calculated from inputs
76  const double M = target.Mass(); // units: GeV
77 
78  const double FF = fFF->FormFactor(Q2, target);
79  const double TT = TE - M;
80 
81  // auxiliary variables
82  const double E2 = E * E;
83  const double Z2 = Z * Z;
84  const double FF2 = FF * FF;
85  const double TTDiff = TT - M;
86 
87  const double const_factor = 2* constants::kPi * constants::kAem ;
88  const double model_params = fEps2 * DTheta2 * fAlpha_D ;
89 
90  const double num_fact1 = FF2 * Z2;
91  const double num_fact21 = fDNuMass2 * (TTDiff - 2.*E);
92  const double num_fact22 = 2. * M * (2.*E2 - 2.*TT*E + TT*TTDiff);
93  const double den_fact1 = 1. / (E2);
94  const double den_fact2 = TMath::Power((fDMediatorMass2 + 2.*TT*M), -2.);
95 
96  if(kps == kPSEDNufE) {
97  const double xsec = const_factor * model_params * num_fact1 *
98  (num_fact21 + num_fact22) * den_fact1 * den_fact2;
99 
100  return xsec;
101  }
102  return 0.;
103 }
104 //____________________________________________________________________________
105 double BertuzzoDNuCOHPXSec::Integral(const Interaction * interaction) const
106 {
107  double xsec = fXSecIntegrator->Integrate(this,interaction);
108  return xsec;
109 }
110 //____________________________________________________________________________
111 bool BertuzzoDNuCOHPXSec::ValidProcess(const Interaction * interaction) const
112 {
113  if(interaction->TestBit(kISkipProcessChk)) return true;
114 
115  const ProcessInfo & proc_info = interaction->ProcInfo();
116  if ( ! proc_info.IsCoherentElastic() ) return false;
117  if ( ! proc_info.IsDarkNeutralCurrent() ) return false ;
118 
119  const InitialState & init_state = interaction->InitState();
120  if ( ! pdg::IsNeutrino( TMath::Abs( init_state.ProbePdg() ) ) ) return false ;
121 
122  const Target & target = init_state.Tgt();
123  if( ! target.IsNucleus() ) return false ;
124 
125  return true;
126 }
127 //____________________________________________________________________________
128 bool BertuzzoDNuCOHPXSec::ValidKinematics(const Interaction* interaction) const
129 {
130  if(interaction->TestBit(kISkipKinematicChk)) return true;
131 
132  if(!interaction->PhaseSpace().IsAboveThreshold()) return false;
133 
134  const double E = interaction->InitState().ProbeE(kRfLab);
135  const double M = interaction->InitState().Tgt().Mass();
136  const TLorentzVector& DNu = interaction->Kine().FSLeptonP4();
137 
138  const double tl = DNu.E()*(M+E) - E*M - 0.5*fDNuMass2;
139  const double tr = E * DNu.P();
140 
141  if(tl < -1.*tr) return false;
142  if(tl > tr) return false;
143 
144  return true;
145 
146 }
147 //____________________________________________________________________________
149 {
150  Algorithm::Configure(config);
151  this->LoadConfig();
152 }
153 //____________________________________________________________________________
155 {
156  Algorithm::Configure(config);
157  this->LoadConfig();
158 }
159 //____________________________________________________________________________
161 {
162 
163  bool good_configuration = true ;
164 
165  double DKineticMixing = 0.;
166  this->GetParam("Dark-KineticMixing", DKineticMixing);
167  fEps2 = DKineticMixing * DKineticMixing;
168 
169  bool force_unitarity = false ;
170  GetParam( "Dark-Mixing-ForceUnitarity", force_unitarity ) ;
171 
172  unsigned int n_min_mixing = force_unitarity ? 3 : 4 ;
173 
174  std::vector<double> DMixing2s; // |U_{\alpha 4}|^2
175  this->GetParamVect("Dark-Mixings2", DMixing2s);
176 
177  // check whether we have enough mixing elements
178  if ( DMixing2s.size () < n_min_mixing ) {
179  good_configuration = false ;
180  LOG("BertuzzoDNuCOH", pERROR )
181  << "Not enough mixing elements specified, only specified "
182  << DMixing2s.size() << " / " << n_min_mixing ;
183  }
184 
185  double tot_mix = 0.;
186  for( unsigned int i = 0; i < n_min_mixing ; ++i ) {
187  if ( DMixing2s[i] < 0. ) {
188  good_configuration = false ;
189  LOG("BertuzzoDNuCOH", pERROR )
190  << "Mixing " << i << " non positive: " << DMixing2s[i] ;
191  continue ;
192  }
193  tot_mix += fMixing2s[i] = DMixing2s[i] ;
194  }
195 
196  if ( force_unitarity ) {
197  fMixing2s[3] = 1. - tot_mix ;
198  }
199  if ( DMixing2s[3] < 0. ) {
200  good_configuration = false ;
201  LOG("BertuzzoDNuCOH", pERROR )
202  << "Mixing D4 non positive: " << DMixing2s[3] ;
203  }
204 
205  this->GetParam("Dark-Alpha", fAlpha_D);
206 
207  fDNuMass = 0.;
208  this->GetParam("Dark-NeutrinoMass", fDNuMass);
210 
211  fDMediatorMass = 0.;
212  this->GetParam("Dark-MediatorMass", fDMediatorMass);
214 
216  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
217  assert(fXSecIntegrator);
218  fFF = dynamic_cast<const EngelFormFactor *> (this->SubAlg("FormFactor"));
219  assert(fFF);
220 
221  if ( ! good_configuration ) {
222  LOG("BertuzzoDNuCOH", pFATAL ) << "Wrong configuration. Exiting" ;
223  exit ( 78 ) ;
224  }
225 
226 }
227 //____________________________________________________________________________
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
bool IsNuTau(int pdgc)
Definition: PDGUtils.cxx:168
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
Cross Section Integrator Interface.
Form Factor for BertuzzoDNuCOHXSec...
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
bool IsDarkNeutralCurrent(void) const
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
bool IsNuE(int pdgc)
Definition: PDGUtils.cxx:158
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Mass(void) const
Definition: Target.cxx:224
void Configure(const Registry &config) override
std::array< double, 4 > fMixing2s
bool IsNuMu(int pdgc)
Definition: PDGUtils.cxx:163
Summary information for an interaction.
Definition: Interaction.h:56
double FormFactor(const double Q, const Target &target) const
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double Integral(const Interaction *i) const override
bool IsCoherentElastic(void) const
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
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
bool ValidProcess(const Interaction *i) const override
Can this cross section algorithm handle the input process?
int Z(void) const
Definition: Target.h:68
const XSecIntegratorI * fXSecIntegrator
cross section integrator
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
bool ValidKinematics(const Interaction *i) const override
Is the input kinematical point a physically allowed one?
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
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 EngelFormFactor * fFF
Engel Form Factor algorithm.
double XSec(const Interaction *i, KinePhaseSpace_t k) const override
Compute the cross section for the input interaction.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345