GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AhrensDMELPXSec.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 
7  Author: Joshua Berger <jberger \at physics.wisc.edu>
8  University of Wisconsin-Madison
9 
10  Costas Andreopoulos <c.andreopoulos \at cern.ch>
11  University of Liverpool
12 */
13 //____________________________________________________________________________
14 
15 #include <TMath.h>
16 
30 #include "Framework/Utils/RunOpt.h"
31 
32 using namespace genie;
33 using namespace genie::utils;
34 using namespace genie::constants;
35 
36 //____________________________________________________________________________
38 XSecAlgorithmI("genie::AhrensDMELPXSec")
39 {
40 
41 }
42 //____________________________________________________________________________
44 XSecAlgorithmI("genie::AhrensDMELPXSec", config)
45 {
46 
47 }
48 //____________________________________________________________________________
50 {
51 
52 }
53 //____________________________________________________________________________
55  const Interaction * interaction, KinePhaseSpace_t kps) const
56 {
57  if(! this -> ValidProcess (interaction) ) return 0.;
58  if(! this -> ValidKinematics (interaction) ) return 0.;
59 
60  const InitialState & init_state = interaction -> InitState();
61  const Kinematics & kinematics = interaction -> Kine();
62  const Target & target = init_state.Tgt();
63 
64  LOG("AhrensDMEL", pDEBUG) << "Using v^" << fVelMode << " dependence";
65 
66  double E = init_state.ProbeE(kRfHitNucRest);
67  double ml = init_state.GetProbeP4(kRfHitNucRest)->M();
68  double Q2 = kinematics.Q2();
69  double M = target.HitNucMass();
70  double M2 = TMath::Power(M, 2.);
71  double E2 = TMath::Power(E, 2.);
72  double ml2 = TMath::Power(ml,2.);
73  LOG("AhrensDMEL", pNOTICE) << "Form factor masses are " << fMv2 << ", " << fMa2;
74  double qmv2 = TMath::Power(1 + Q2/fMv2, 2);
75  double qma2 = TMath::Power(1 + Q2/fMa2, 2);
76 
77  //-- handle terms changing sign for antineutrinos and isospin rotations
78  int nusign = 1;
79  int nucsign = 1;
80  int nupdgc = init_state.ProbePdg();
81  int nucpdgc = target.HitNucPdg();
82  if( pdg::IsAntiDarkMatter(nupdgc) ) nusign = -1;
83  if( pdg::IsNeutron(nucpdgc) ) nucsign = -1;
84 
85  LOG("AhrensDMEL", pNOTICE) << "Calculating for nuclear sign " << nucsign;
86 
87  //-- compute up quark form factor terms
88  double Geu = fQuV * (1.5 + nucsign*0.5) / qmv2;
89  double Gmu = fQuV * ((1.5 + nucsign*0.5) * fMuP + (1.5 - nucsign*0.5) * fMuN) / qmv2;
90  double FAu = fQuA * (nucsign > 0 ? fDelu : fDeld) / qma2;
91 
92  //-- compute down quark form factor terms
93  double Ged = fQdV * (1.5 - nucsign*0.5) / qmv2;
94  double Gmd = fQdV * ((1.5 - nucsign*0.5) * fMuP + (1.5 + nucsign*0.5) * fMuN) / qmv2;
95  double FAd = fQdA * (nucsign > 0 ? fDeld : fDelu) / qma2;
96 
97  //-- compute the induced pseudoscalar form factors
98  double pole3 = 4.0*M2 / (Q2 + fMpi2);
99  double pole0 = 4.0*M2 / (Q2 + fMeta2);
100  double FPu = 0.5 * (pole3 * (FAu - FAd) + pole0 * (FAu + FAd));
101  double FPd = 0.5 * (- pole3 * (FAu - FAd) + pole0 * (FAu + FAd));
102 
103  //-- compute strange quark form factor terms
104  double Ges = 0.0;
105  double Gms = 0.0;
106  double FAs = fQsA * fDels / qma2;
107  double FPs = 0.0; // fQsA * 2.0 * M2 * fDels / qmp2 / fMpi2;
108 
109  //-- compute form factors
110  double Ge = Geu + Ged + Ges;
111  double Gm = Gmu + Gmd + Gms;
112  double FA = FAu + FAd + FAs;
113  double FP = FPu + FPd + FPs;
114  double tau = 0.25 * Q2/M2;
115  double F1 = (Ge + tau * Gm) / (1.0 + tau);
116  double F2 = (Gm - Ge) / (1.0 + tau);
117  double F12 = TMath::Power(F1,2);
118  double F22 = TMath::Power(F2,2);
119  double FA2 = TMath::Power(FA,2);
120 
121  //-- compute the free nucleon cross section
122  double xsec = 0.;
123  double del = ml2 / M2;
124  double AT_F1F1 = 0.;
125  double AT_F2F2 = 0.;
126  double AT_FAFA = 0.;
127  double AT_F1F2 = 0.;
128  double AL = 0.;
129  double B = 0.;
130  double C = 0.;
131  if (fVelMode == 0) {
132  double QchiV2 = TMath::Power(fQchiV,2);
133  double QchiA2 = TMath::Power(fQchiA,2);
134  C = (QchiA2 + QchiV2) * (F12 + F22 * tau + FA2);
135  B = 8. * fQchiA * fQchiV * tau * FA * (F1 + F2);
136  AL = 16. * QchiA2 * del * TMath::Power(tau*(FA - 2.*FP*tau),2);
137  AT_F1F1 = QchiA2*(tau-1.)*(del+tau) + QchiV2*tau*(-del+tau-1);
138  AT_F2F2 = -tau*(QchiA2*(tau-1.)*(del+tau) + QchiV2*(del + (tau-1.)*tau));
139  AT_FAFA = (1.+tau)*(QchiA2*(del+tau) + QchiV2*(tau-del));
140  AT_F1F2 = 2.*tau*(2.*QchiA2*(del+tau) - QchiV2*(del-2.*tau));
141  }
142  else if (fVelMode == 2) {
143  double QchiS2 = TMath::Power(fQchiS,2);
144  C = QchiS2 * (F12 + F22 * tau + FA2);
145  AT_F1F1 = -QchiS2 * tau * (del + tau);
146  AT_F2F2 = AT_F1F1;
147  AT_FAFA = -QchiS2 * (tau + 1.) * (del + tau);
148  AT_F1F2 = 2.*AT_F1F1;
149  }
150  double smu = E/M - tau;
151  double MZ2 = TMath::Power(fMedMass,2);
152  double lon = TMath::Power(M2 / MZ2 + 0.25/tau,2);
153  LOG("AhrensDMEL", pDEBUG)
154  << "Using a mediator mass of " << fMedMass;
155  // double fd = 8*ml2*M2*tau*(2*M2*tau+MZ2) / (MZ2*MZ2*E2);
156  double gZp = fgZp;
157  double gZp4 = TMath::Power(gZp,4);
158  double prop = 1. / (Q2 + MZ2);
159  double prop2 = TMath::Power(prop,2);
160  double xsec0 = gZp4 * M2 * prop2 / (4. * kPi * (E2 - ml2));
161  xsec = xsec0 * (AL * lon + AT_F1F1 * F12 + AT_F2F2 * F22 + AT_FAFA * FA2 + AT_F1F2 * F1 * F2 + nusign * B * smu + C * smu * smu);
162 
163  LOG("AhrensDMEL", pDEBUG)
164  << "dXSec[vN,El]/dQ2 [FreeN](Ev = "<< E<< ", Q2 = "<< Q2 << ") = "<< xsec;
165 
166  //-- The algorithm computes dxsec/dQ2
167  // Check whether variable tranformation is needed
168  if(kps!=kPSQ2fE) {
169  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
170  xsec *= J;
171  }
172 
173  //-- if requested return the free nucleon xsec even for input nuclear tgt
174  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
175 
176  //-- compute nuclear suppression factor
177  // (R(Q2) is adapted from NeuGEN - see comments therein)
178  double R = nuclear::NuclQELXSecSuppression("Default", 0.5, interaction);
179 
180  //-- number of scattering centers in the target
181  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
182 
183  LOG("AhrensDMEL", pDEBUG)
184  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
185 
186  //-- compute nuclear cross section
187  xsec *= (R*NNucl);
188 
189  return xsec;
190 }
191 //____________________________________________________________________________
192 double AhrensDMELPXSec::Integral(const Interaction * interaction) const
193 {
194  double xsec = fXSecIntegrator->Integrate(this,interaction);
195  return xsec;
196 }
197 //____________________________________________________________________________
198 bool AhrensDMELPXSec::ValidProcess(const Interaction * interaction) const
199 {
200  if(interaction->TestBit(kISkipProcessChk)) return true;
201  return true;
202 }
203 //____________________________________________________________________________
205 {
206  Algorithm::Configure(config);
207  this->LoadConfig();
208 }
209 //____________________________________________________________________________
210 void AhrensDMELPXSec::Configure(string config)
211 {
212  Algorithm::Configure(config);
213  this->LoadConfig();
214 }
215 //____________________________________________________________________________
217 {
218  // dark matter couplings to mediator
219  double QchiL, QchiR;
220  this->GetParam( "DarkLeftCharge", QchiL ) ;
221  this->GetParam( "DarkRightCharge", QchiR ) ;
222  this->GetParam( "DarkScalarCharge", fQchiS ) ;
223  fQchiV = 0.5*(QchiL + QchiR);
224  fQchiA = 0.5*(- QchiL + QchiR);
225 
226  // quark couplings to mediator
227  double QuL, QuR, QdL, QdR, QsL, QsR;
228  this->GetParam( "UpLeftCharge", QuL ) ;
229  this->GetParam( "UpRightCharge", QuR ) ;
230  this->GetParam( "DownLeftCharge", QdL ) ;
231  this->GetParam( "DownRightCharge", QdR ) ;
232  this->GetParam( "StrangeLeftCharge", QsL ) ;
233  this->GetParam( "StrangeRightCharge", QsR ) ;
234  fQuV = 0.5*(QuL + QuR);
235  fQuA = 0.5*(- QuL + QuR);
236  fQdV = 0.5*(QdL + QdR);
237  fQdA = 0.5*(- QdL + QdR);
238  fQsV = 0.5*(QsL + QsR);
239  fQsA = 0.5*(- QsL + QsR);
240 
241  // axial and vector masses
242  double ma, mv, mp, mpi, meta ;
243  this->GetParam( "QEL-Ma", ma ) ;
244  this->GetParam( "QEL-Mv", mv ) ;
245  this->GetParam( "DMEL-Mp", mp ) ;
246  this->GetParam( "DMEL-Mpi", mpi ) ;
247  this->GetParam( "DMEL-Meta", meta ) ;
248  fMa2 = TMath::Power(ma,2);
249  fMv2 = TMath::Power(mv,2);
250  fMp2 = TMath::Power(mp,2);
251  fMpi2 = TMath::Power(mpi,2);
252  fMeta2 = TMath::Power(meta,2);
253 
254  // anomalous magnetic moments
255  this->GetParam( "AnomMagnMoment-P", fMuP ) ;
256  this->GetParam( "AnomMagnMoment-N", fMuN ) ;
257 
258  // Axial-vector spin charge
259  // Since we have a more general axial dependence,
260  // we need a more complex treatment than the usual model
261  this->GetParam( "AxialVectorSpin-u", fDelu );
262  this->GetParam( "AxialVectorSpin-d", fDeld );
263  this->GetParam( "AxialVectorSpin-s", fDels );
264 
265  // velocity dependence of interaction
266  this->GetParamDef("velocity-mode", fVelMode, 0 );
267 
268  // mediator coupling
269  this->GetParam("ZpCoupling", fgZp ) ;
270 
271  // mediator mass
273 
274  // load XSec Integrator
276  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
277  assert(fXSecIntegrator);
278 }
279 //____________________________________________________________________________
Cross Section Calculation Interface.
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
const int kPdgMediator
Definition: PDGCodes.h:220
double HitNucMass(void) const
Definition: Target.cxx:233
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
enum genie::EKinePhaseSpace KinePhaseSpace_t
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?
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsAntiDarkMatter(int pdgc)
Definition: PDGUtils.cxx:133
double Integral(const Interaction *i) const
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
int Z(void) const
Definition: Target.h:68
const XSecIntegratorI * fXSecIntegrator
int N(void) const
Definition: Target.h:69
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
#define pNOTICE
Definition: Messenger.h:61
bool GetParamDef(const RgKey &name, T &p, const T &def) const
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 XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void Configure(const Registry &config)
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
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