Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::EmpiricalMECPXSec2015 Class Reference

Computes the MEC differential cross section. Is a concrete implementation of the XSecAlgorithmI interface.
. More...

#include <EmpiricalMECPXSec2015.h>

Inheritance diagram for genie::EmpiricalMECPXSec2015:
genie::XSecAlgorithmI genie::Algorithm

Public Member Functions

 EmpiricalMECPXSec2015 ()
 
 EmpiricalMECPXSec2015 (string config)
 
virtual ~EmpiricalMECPXSec2015 ()
 
double XSec (const Interaction *i, KinePhaseSpace_t k) const
 Compute the cross section for the input interaction. More...
 
double Integral (const Interaction *i) const
 
bool ValidProcess (const Interaction *i) const
 Can this cross section algorithm handle the input process? More...
 
void Configure (const Registry &config)
 
void Configure (string param_set)
 
- Public Member Functions inherited from genie::XSecAlgorithmI
virtual ~XSecAlgorithmI ()
 
virtual bool ValidKinematics (const Interaction *i) const
 Is the input kinematical point a physically allowed one? More...
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

void LoadConfig (void)
 

Private Attributes

double fMq2d
 toy model param: `mass' in dipole (Q2 - dependence) form factor (GeV) More...
 
double fMass
 toy model param: peak of W distribution (GeV) More...
 
double fWidth
 toy model param: width of W distribution (GeV) More...
 
double fMECAPower
 power of A relative to carbon More...
 
double fFracPN_NC
 toy model param: fraction of nucleon pairs that are pn, not nn or pp More...
 
double fFracPN_CC
 toy model param: fraction of nucleon pairs that are pn, not nn or pp More...
 
double fFracPN_EM
 toy model param: fraction of nucleon pairs that are pn, not nn or pp More...
 
double fFracCCQE
 empirical model param: MEC cross section is taken to be this fraction of CCQE cross section More...
 
double fFracNCQE
 empirical model param: MEC cross section is taken to be this fraction of NCQE cross section More...
 
double fFracEMQE
 empirical model param: MEC cross section is taken to be this fraction of Rosenbluth xs More...
 
const XSecAlgorithmIfXSecAlgCCQE
 cross section algorithm for CCQE More...
 
const XSecAlgorithmIfXSecAlgNCQE
 cross section algorithm for NCQE More...
 
const XSecAlgorithmIfXSecAlgEMQE
 cross section algorithm for EMQE More...
 
const XSecIntegratorIfXSecIntegrator
 Integrator used for reweighting. More...
 
bool fIntegrateForReweighting
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::XSecAlgorithmI
 XSecAlgorithmI ()
 
 XSecAlgorithmI (string name)
 
 XSecAlgorithmI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Computes the MEC differential cross section. Is a concrete implementation of the XSecAlgorithmI interface.
.

Author
Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

Steve Dytman <dytman+ pitt.edu> Pittsburgh University

May 05, 2009

Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 30 of file EmpiricalMECPXSec2015.h.

Constructor & Destructor Documentation

EmpiricalMECPXSec2015::EmpiricalMECPXSec2015 ( )

Definition at line 34 of file EmpiricalMECPXSec2015.cxx.

34  :
35 XSecAlgorithmI("genie::EmpiricalMECPXSec2015")
36 {
37 
38 }
EmpiricalMECPXSec2015::EmpiricalMECPXSec2015 ( string  config)

Definition at line 40 of file EmpiricalMECPXSec2015.cxx.

40  :
41 XSecAlgorithmI("genie::EmpiricalMECPXSec2015", config)
42 {
43 
44 }
static Config * config
Definition: config.cpp:1054
EmpiricalMECPXSec2015::~EmpiricalMECPXSec2015 ( )
virtual

Definition at line 46 of file EmpiricalMECPXSec2015.cxx.

47 {
48 
49 }

Member Function Documentation

void EmpiricalMECPXSec2015::Configure ( const Registry config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 353 of file EmpiricalMECPXSec2015.cxx.

354 {
355  Algorithm::Configure(config);
356  this->LoadConfig();
357 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void EmpiricalMECPXSec2015::Configure ( string  config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 359 of file EmpiricalMECPXSec2015.cxx.

360 {
362 
363  Registry* algos = AlgConfigPool::Instance() -> GlobalParameterList() ;
364  string global_key_head = "XSecModel@genie::EventGenerator/" ;
365  string local_key_head = "XSecModel-" ;
366 
367  Registry r( "EmpiricalMECPXSec2015_specific", false ) ;
368  r.Set( local_key_head + "QEL-NC", algos -> GetAlg( global_key_head + "QEL-NC") ) ;
369  r.Set( local_key_head + "QEL-CC", algos -> GetAlg( global_key_head + "QEL-CC") ) ;
370  r.Set( local_key_head + "QEL-EM", algos -> GetAlg( global_key_head + "QEL-EM") ) ;
371 
373 
374  this->LoadConfig();
375 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
static AlgConfigPool * Instance()
double EmpiricalMECPXSec2015::Integral ( const Interaction i) const
virtual

Integrate the model over the kinematic phase space available to the input interaction (kinematical cuts can be included)

Implements genie::XSecAlgorithmI.

Definition at line 203 of file EmpiricalMECPXSec2015.cxx.

204 {
205  // Normally the empirical MEC splines are computed as a fraction of the CCQE
206  // ones. In cases where we want to integrate the differential cross section
207  // directly (e.g., for reweighting via the XSecShape_CCMEC dial), do that
208  // instead.
209  if ( fIntegrateForReweighting ) {
210  return fXSecIntegrator->Integrate( this, interaction );
211  }
212 
213 // Calculate the CCMEC cross section as a fraction of the CCQE cross section
214 // for the given nuclear target at the given energy.
215 // Alternative strategy is to calculate the MEC cross section as the difference
216 // of CCQE cross section for two different M_A values (eg ~1.3 GeV and ~1.0 GeV)
217 // Include hit-object combinatorial factor? Would yield different A-dependence
218 // for MEC and QE.
219 //
220 
221  bool iscc = interaction->ProcInfo().IsWeakCC();
222  bool isnc = interaction->ProcInfo().IsWeakNC();
223  bool isem = interaction->ProcInfo().IsEM();
224 
225  int nupdg = interaction->InitState().ProbePdg();
226  int tgtpdg = interaction->InitState().Tgt().Pdg();
227  double E = interaction->InitState().ProbeE(kRfLab);
228  int nucleon_cluster_pdg = interaction->InitState().Tgt().HitNucPdg();
229  double Z=interaction->InitState().Tgt().Z();
230  double A=interaction->InitState().Tgt().A();
231  double N=A-Z;
232 
233  if(iscc) {
234 
235  int nucpdg = 0;
236  // neutrino CC: calculate the CCQE cross section resetting the
237  // hit nucleon cluster to neutron
238  if(pdg::IsNeutrino(nupdg)) {
239  nucpdg = kPdgNeutron;
240  }
241  // anti-neutrino CC: calculate the CCQE cross section resetting the
242  // hit nucleon cluster to proton
243  else
244  if(pdg::IsAntiNeutrino(nupdg)) {
245  nucpdg = kPdgProton;
246  }
247  else {
248  exit(1);
249  }
250 
251  // Create a tmp QE process
252  Interaction * in = Interaction::QELCC(tgtpdg,nucpdg,nupdg,E);
253 
254  // Calculate cross section for the QE process
255  double xsec = fXSecAlgCCQE->Integral(in);
256 
257  // Add A dependence which is not known from theory
258  double fFracADep = 1.;
259  if(A>=12) fFracADep = TMath::Power((N/6.),fMECAPower-1.);
260 
261  // Use tunable fraction
262  // FFracCCQE is fraction of QE going to MEC
263  // fFracCCQE_cluster is fraction of MEC going to each NN pair
264  // double fFracCCQE = fFracCCQElo;
265 
266  double fFracCCQE_cluster=0.;
267  if(pdg::IsNeutrino(nupdg) && nucleon_cluster_pdg==2000000201) fFracCCQE_cluster= fFracPN_CC; //n+p
268  if(pdg::IsNeutrino(nupdg) && nucleon_cluster_pdg==2000000200) fFracCCQE_cluster= 1.0-fFracPN_CC; //n+n
269  if(pdg::IsAntiNeutrino(nupdg) && nucleon_cluster_pdg==2000000201) fFracCCQE_cluster= fFracPN_CC; //n+p
270  if(pdg::IsAntiNeutrino(nupdg) && nucleon_cluster_pdg==2000000202) fFracCCQE_cluster= 1.0-fFracPN_CC; //p+p
271 
272 
273  xsec *= fFracCCQE*fFracCCQE_cluster*fFracADep;
274 
275  // Use gross combinatorial factor (number of 2-nucleon targets over number
276  // of 1-nucleon targets) : (A-1)/2
277  // double combfact = (in->InitState().Tgt().A()-1)/2.;
278  // xsec *= combfact;
279 
280  delete in;
281  return xsec;
282  }
283 
284  else if(isnc) {
285  int nucpdg = kPdgProton;
286  // Create a tmp QE process
287  Interaction * inp = Interaction::QELNC(tgtpdg,nucpdg,nupdg,E);
288  nucpdg = kPdgNeutron;
289  // Create a tmp QE process
290  Interaction * inn = Interaction::QELNC(tgtpdg,nucpdg,nupdg,E);
291 
292  // Calculate cross section for the QE process - avg of p and n - best for isoscalar nuclei
293  double xsec = (Z*fXSecAlgNCQE->Integral(inp) + N*fXSecAlgNCQE->Integral(inn))/A;
294 
295  // Add A dependence which is not known from theory
296  double fFracADep = 1.;
297  if(A>=12) fFracADep = TMath::Power((A/12.),fMECAPower-1.);
298 
299  // Use tunable fraction
300  // FFracNCQE is fraction of QE going to MEC
301  // fFracNCQE_cluster is fraction of MEC going to each NN pair
302  double fFracNCQE_cluster=0.;
303  if(nucleon_cluster_pdg==2000000200) fFracNCQE_cluster= 0.5*(1-fFracPN_NC); //n+n
304  if(nucleon_cluster_pdg==2000000201) fFracNCQE_cluster= fFracPN_NC; //n+p
305  if(nucleon_cluster_pdg==2000000202) fFracNCQE_cluster= 0.5*(1-fFracPN_NC); //p+p
306  xsec *= fFracNCQE*fFracNCQE_cluster*fFracADep;
307  delete inn;
308  delete inp;
309  return xsec;
310  }
311 
312  else if(isem) {
313  int nucpdg = kPdgProton;
314  // Create a tmp QE process
315  Interaction * inp = Interaction::QELEM(tgtpdg,nucpdg,nupdg,E);
316  nucpdg = kPdgNeutron;
317  // Create a tmp QE process
318  Interaction * inn = Interaction::QELEM(tgtpdg,nucpdg,nupdg,E);
319 
320  // Calculate cross section for the QE process - avg of p and n - best for isoscalar nuclei
321  double xsec = (Z*fXSecAlgEMQE->Integral(inp) + N*fXSecAlgEMQE->Integral(inn))/A;
322 
323  // Add A dependence which is not known from theory, data wants high A suppression
324  double fFracADep = 1.;
325  if(A>=12) fFracADep = TMath::Power((A/12.),fMECAPower-1.);
326 
327  // Use tunable fraction
328  // FFracEMQE is fraction of QE going to MEC
329  // fFracEMQE_cluster is fraction of MEC going to each NN pair
330  double fFracEMQE_cluster=0.;
331  if(nucleon_cluster_pdg==2000000200) fFracEMQE_cluster= 0.5*(1-fFracPN_EM); //n+n
332  if(nucleon_cluster_pdg==2000000201) fFracEMQE_cluster= fFracPN_EM; //n+p
333  if(nucleon_cluster_pdg==2000000202) fFracEMQE_cluster= 0.5*(1-fFracPN_EM); //p+p
334  xsec *= fFracEMQE*fFracEMQE_cluster*fFracADep;
335  delete inn;
336  delete inp;
337  return xsec;
338  }
339 
340  return 0;
341 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
double fFracPN_NC
toy model param: fraction of nucleon pairs that are pn, not nn or pp
const XSecAlgorithmI * fXSecAlgCCQE
cross section algorithm for CCQE
const XSecAlgorithmI * fXSecAlgNCQE
cross section algorithm for NCQE
static Interaction * QELCC(int tgt, int nuc, int probe, double E=0)
const XSecIntegratorI * fXSecIntegrator
Integrator used for reweighting.
Summary information for an interaction.
Definition: Interaction.h:56
double fFracPN_CC
toy model param: fraction of nucleon pairs that are pn, not nn or pp
static Interaction * QELNC(int tgt, int nuc, int probe, double E=0)
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
virtual double Integral(const Interaction *i) const =0
double fFracPN_EM
toy model param: fraction of nucleon pairs that are pn, not nn or pp
const XSecAlgorithmI * fXSecAlgEMQE
cross section algorithm for EMQE
static Interaction * QELEM(int tgt, int nuc, int probe, double E=0)
double fFracEMQE
empirical model param: MEC cross section is taken to be this fraction of Rosenbluth xs ...
double fFracNCQE
empirical model param: MEC cross section is taken to be this fraction of NCQE cross section ...
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
#define A
Definition: memgrp.cpp:38
double fFracCCQE
empirical model param: MEC cross section is taken to be this fraction of CCQE cross section ...
const int kPdgProton
Definition: PDGCodes.h:81
double fMECAPower
power of A relative to carbon
const int kPdgNeutron
Definition: PDGCodes.h:83
void EmpiricalMECPXSec2015::LoadConfig ( void  )
private

Definition at line 377 of file EmpiricalMECPXSec2015.cxx.

378 {
379  fXSecAlgCCQE = 0;
380  fXSecAlgNCQE = 0;
381  fXSecAlgEMQE = 0;
382 
383  GetParam( "EmpiricalMEC-Mq2d", fMq2d ) ;
384  GetParam( "EmpiricalMEC-Mass", fMass ) ;
385  GetParam( "EmpiricalMEC-Width", fWidth ) ;
386  GetParam( "EmpiricalMEC-APower", fMECAPower ) ;
387 
388  GetParam( "EmpiricalMEC-FracPN_NC", fFracPN_NC ) ;
389  GetParam( "EmpiricalMEC-FracPN_CC", fFracPN_CC ) ;
390  GetParam( "EmpiricalMEC-FracPN_EM", fFracPN_EM ) ;
391 
392  GetParam( "EmpiricalMEC-FracCCQE", fFracCCQE ) ;
393  GetParam( "EmpiricalMEC-FracNCQE", fFracNCQE ) ;
394  GetParam( "EmpiricalMEC-FracEMQE", fFracEMQE ) ;
395 
396  string key_head = "XSecModel-" ;
397 
398  fXSecAlgNCQE =
399  dynamic_cast<const XSecAlgorithmI *> ( this -> SubAlg( key_head + "QEL-NC" ) ) ;
400  assert(fXSecAlgNCQE);
401 
402  fXSecAlgCCQE =
403  dynamic_cast<const XSecAlgorithmI *> ( this -> SubAlg( key_head + "QEL-CC" ) ) ;
404  assert(fXSecAlgCCQE);
405 
406  fXSecAlgEMQE =
407  dynamic_cast<const XSecAlgorithmI *> ( this -> SubAlg( key_head + "QEL-EM" ) ) ;
408  assert(fXSecAlgEMQE);
409 
410  // Get the "fast" configuration of MECXSec. This will be used when integrating
411  // the total cross section for reweighting purposes.
413  fXSecIntegrator = dynamic_cast<const XSecIntegratorI*>( algf->AdoptAlgorithm(
414  "genie::MECXSec", "Fast") );
415  assert( fXSecIntegrator );
416 
417  fIntegrateForReweighting = false;
418  GetParamDef( "IntegrateForReweighting", fIntegrateForReweighting, false );
419 }
Cross Section Calculation Interface.
Cross Section Integrator Interface.
double fFracPN_NC
toy model param: fraction of nucleon pairs that are pn, not nn or pp
const XSecAlgorithmI * fXSecAlgCCQE
cross section algorithm for CCQE
const XSecAlgorithmI * fXSecAlgNCQE
cross section algorithm for NCQE
const XSecIntegratorI * fXSecIntegrator
Integrator used for reweighting.
double fMass
toy model param: peak of W distribution (GeV)
double fFracPN_CC
toy model param: fraction of nucleon pairs that are pn, not nn or pp
double fWidth
toy model param: width of W distribution (GeV)
Algorithm * AdoptAlgorithm(const AlgId &algid) const
Definition: AlgFactory.cxx:116
double fMq2d
toy model param: `mass&#39; in dipole (Q2 - dependence) form factor (GeV)
double fFracPN_EM
toy model param: fraction of nucleon pairs that are pn, not nn or pp
const XSecAlgorithmI * fXSecAlgEMQE
cross section algorithm for EMQE
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
double fFracEMQE
empirical model param: MEC cross section is taken to be this fraction of Rosenbluth xs ...
double fFracNCQE
empirical model param: MEC cross section is taken to be this fraction of NCQE cross section ...
double fFracCCQE
empirical model param: MEC cross section is taken to be this fraction of CCQE cross section ...
double fMECAPower
power of A relative to carbon
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
bool EmpiricalMECPXSec2015::ValidProcess ( const Interaction i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 343 of file EmpiricalMECPXSec2015.cxx.

344 {
345  if(interaction->TestBit(kISkipProcessChk)) return true;
346 
347  const ProcessInfo & proc_info = interaction->ProcInfo();
348  if(!proc_info.IsMEC()) return false;
349 
350  return true;
351 }
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsMEC(void) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double EmpiricalMECPXSec2015::XSec ( const Interaction i,
KinePhaseSpace_t  k 
) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 51 of file EmpiricalMECPXSec2015.cxx.

53 {
54 
55  // If we've been asked for the kPSTlctl phase space, get W and Q^2 from those
56  // variables. You actually need the lepton phi set in order to do the
57  // conversion (see the "important note" in src/Framework/Utils/KineUtils.cxx
58  // about the kPSTlctl --> kPSWQ2fE Jacobian) when Fermi motion is properly
59  // taken into account. For now, I neglect Fermi motion as a stopgap solution.
60  // - S. Gardiner, 29 July 2020
61  if ( kps == kPSTlctl ) {
62 
63  // {Tl, ctl} --> {W, Q2}
64 
65  // Probe properties (mass, energy, momentum)
66  const InitialState& init_state = interaction->InitState();
67  double mv = init_state.Probe()->Mass();
68  double Ev = init_state.ProbeE( kRfLab );
69  double pv = std::sqrt( std::max(0., Ev*Ev - mv*mv) );
70 
71  // Invariant mass of the initial hit nucleon
72  const TLorentzVector& hit_nuc_P4 = init_state.Tgt().HitNucP4();
73  double M = hit_nuc_P4.M();
74 
75  // Outgoing lepton mass
76  double ml = interaction->FSPrimLepton()->Mass();
77 
78  // Lab-frame lepton kinetic energy and scattering cosine
79  double Tl = interaction->Kine().GetKV( kKVTl );
80  double ctl = interaction->Kine().GetKV( kKVctl );
81 
82  // Q^2 from Tl and ctl
83  double El = Tl + ml;
84  double pl = std::sqrt( Tl*Tl + 2.*ml*Tl );
85  double Q2 = -mv*mv - ml*ml + 2.*Ev*El - 2.*pv*pl*ctl;
86 
87  // Energy transfer
88  double omega = Ev - El;
89 
90  double W = std::sqrt( std::max(0., M*M - Q2 + 2.*omega*M) );
91 
92  interaction->KinePtr()->SetW( W );
93  interaction->KinePtr()->SetQ2( Q2 );
94  }
95 
96 
97 // meson exchange current contribution depends a lot on QE model.
98 // This is an empirical model in development, not used in default event generation.
99 
100 // if(! this -> ValidProcess (interaction) ) return 0.;
101 // if(! this -> ValidKinematics (interaction) ) return 0.;
102  bool iscc = interaction->ProcInfo().IsWeakCC();
103  bool isnc = interaction->ProcInfo().IsWeakNC();
104  bool isem = interaction->ProcInfo().IsEM();
105 
106  const Kinematics & kinematics = interaction -> Kine();
107  double W = kinematics.W();
108  double Q2 = kinematics.Q2();
109  //LOG("MEC", pINFO) << "W, Q2 trial= " << W << " " << Q2 ;
110 
111  //
112  // Do a check whether W,Q2 is allowed. Return 0 otherwise.
113  //
114  double Ev = interaction->InitState().ProbeE(kRfHitNucRest); // kRfLab
115  int nucleon_cluster_pdg = interaction->InitState().Tgt().HitNucPdg();
116  double M2n = PDGLibrary::Instance()->Find(nucleon_cluster_pdg)-> Mass(); // nucleon cluster mass
117  double M2n2 = M2n*M2n;
118  double ml = interaction->FSPrimLepton()->Mass();
119  Range1D_t Wlim = isem ? genie::utils::kinematics::electromagnetic::InelWLim(Ev, ml, M2n) : genie::utils::kinematics::InelWLim(Ev, M2n, ml);
120 
121  //LOG("MEC", pINFO) << "Ev, ml, M2n = " << Ev << " " << ml << " " << M2n;
122  //LOG("MEC", pINFO) << "Wlim= " << Wlim.min << " " <<Wlim.max ;
123  if(W < Wlim.min || W > Wlim.max)
124  {double xsec = 0.;
125  return xsec;
126  }
127  //use proper Q2 limit from Controls.h
128  Range1D_t Q2lim = isem ? genie::utils::kinematics::electromagnetic::InelQ2Lim_W(Ev, ml, M2n, W) : genie::utils::kinematics::InelQ2Lim_W (Ev, M2n, ml, W, kMinQ2Limit);
129 
130  //LOG("MEC", pINFO) << "Q2lim= " << Q2lim.min << " " <<Q2lim.max ;
131  if(Q2 < Q2lim.min || Q2 > Q2lim.max)
132  {double xsec = 0.;
133  return xsec;
134  }
135 
136  //get x and y
137  double x = 0.;
138  double y = 0.;
139  genie::utils::kinematics::WQ2toXY(Ev,M2n,W,Q2,x,y);
140  // LOG("MEC", pINFO) << "x = " << x << ", y = " << y;
141  // double Tmu = (1.-y)*Ev; // UNUSED - comment to quiet compiler warnings
142 
143  // Calculate d^2xsec/dWdQ2 - first form factor which is common to both
144  double Wdep = TMath::Gaus(W, fMass, fWidth);
145  double Q2dep = Q2*TMath::Power((1+Q2/fMq2d),-8.);
146  // double nudep = TMath::Power(Tmu,2.5);
147  // LOG("MEC", pINFO) << "Tmu = " << Tmu << ", nudep = " << nudep;
148  double FF2 = Wdep * Q2dep;// * nudep;
149  //LOG("MEC", pINFO) << "form factor = " << FF2 << ", Q2 = " << Q2 << ", W = " << W;
150 
151 // using formulas in Bodek and Budd for (e,e') inclusive cross section
152  double xsec = 1.;
153  if(isem) {
154  // Calculate scattering angle
155  //
156  // Q^2 = 4 * E^2 * sin^2 (theta/2) / ( 1 + 2 * (E/M) * sin^2(theta/2) ) =>
157  // sin^2 (theta/2) = MQ^2 / (4ME^2 - 2EQ^2)
158 
159  double E = Ev;
160  double E2 = E*E;
161  double sin2_halftheta = M2n*Q2 / (4*M2n*E2 - 2*E*Q2);
162  // double sin4_halftheta = TMath::Power(sin2_halftheta, 2.);
163  double cos2_halftheta = 1.-sin2_halftheta;
164  // double cos_halftheta = TMath::Sqrt(cos2_halftheta);
165  double tan2_halftheta = sin2_halftheta/cos2_halftheta;
166  double Q4 = Q2*Q2;
167 
168  // Calculate tau and the virtual photon polarization (epsilon)
169  double tau = Q2/(4*M2n2);
170  // double epsilon = 1. / (1. + 2.*(tau/x))*tan2_halftheta); //different than RosenbluthPXSec.cxx
171 
172  // Calculate the scattered lepton energy
173  double Ep = E / (1. + 2.*(E/M2n)*sin2_halftheta);
174  double Ep2 = Ep*Ep;
175 
176  //calculate cross section - d2sig/dOmega dE for purely transverse process
177  xsec = 4*kAem2*Ep2*cos2_halftheta/Q4 * FF2 * (tau/(1+tau) +2*tau*tan2_halftheta);
178  }
179  // use BB formula which seems to be same as Llewlyn-Smith
180  // note B term is only difference between nu and antinu, so both same here
181  else if(isnc||iscc){
182  double tau = Q2/(4*M2n2);
183  double tau2 = tau*tau;
184  double smufac = 4*M2n*Ev - Q2 - ml*ml;
185  double A = (ml*ml+Q2)/M2n2 * (tau*(1+tau) - tau2*(1-tau)+4*tau2)/TMath::Power(1+tau,2.) * FF2;
186  double C = tau/4/(1+tau) * FF2;
187  xsec = A + smufac*smufac*C; // CC or NC case - Llewelyn-Smith for transverse vector process.
188  }
189  // Check whether variable tranformation is needed
190  if ( kps!=kPSWQ2fE && xsec != 0. ) {
192 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
193  LOG("MEC", pDEBUG)
194  << "Jacobian for transformation to: "
195  << KinePhaseSpace::AsString(kps) << ", J = " << J;
196 #endif
197  xsec *= J;
198  }
199 
200  return xsec;
201 }
double W(bool selected=false) const
Definition: Kinematics.cxx:157
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
Range1D_t InelWLim(double El, double ml, double M)
Definition: KineUtils.cxx:534
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
Range1D_t InelWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:345
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
TParticlePDG * Probe(void) const
double Mass(Resonance_t res)
resonance mass (GeV)
double fMass
toy model param: peak of W distribution (GeV)
static const double kMinQ2Limit
Definition: Controls.h:41
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:366
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static const double kAem2
Definition: Constants.h:57
static string AsString(KinePhaseSpace_t kps)
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1119
double fWidth
toy model param: width of W distribution (GeV)
static int max(int a, int b)
double fMq2d
toy model param: `mass&#39; in dipole (Q2 - dependence) form factor (GeV)
double max
Definition: Range1.h:53
Range1D_t InelQ2Lim_W(double El, double ml, double M, double W)
Definition: KineUtils.cxx:556
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
E
Definition: 018_def.c:13
Definition: utils.py:1
Definition: 018_def.c:13
#define A
Definition: memgrp.cpp:38
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
list x
Definition: train.py:276
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
const Target & Tgt(void) const
Definition: InitialState.h:66
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

double genie::EmpiricalMECPXSec2015::fFracCCQE
private

empirical model param: MEC cross section is taken to be this fraction of CCQE cross section

Definition at line 60 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fFracEMQE
private

empirical model param: MEC cross section is taken to be this fraction of Rosenbluth xs

Definition at line 62 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fFracNCQE
private

empirical model param: MEC cross section is taken to be this fraction of NCQE cross section

Definition at line 61 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fFracPN_CC
private

toy model param: fraction of nucleon pairs that are pn, not nn or pp

Definition at line 57 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fFracPN_EM
private

toy model param: fraction of nucleon pairs that are pn, not nn or pp

Definition at line 58 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fFracPN_NC
private

toy model param: fraction of nucleon pairs that are pn, not nn or pp

Definition at line 56 of file EmpiricalMECPXSec2015.h.

bool genie::EmpiricalMECPXSec2015::fIntegrateForReweighting
private

Whether to integrate in the usual way (false) or in "reweighting mode" (true)

Definition at line 73 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fMass
private

toy model param: peak of W distribution (GeV)

Definition at line 52 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fMECAPower
private

power of A relative to carbon

Definition at line 54 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fMq2d
private

toy model param: `mass' in dipole (Q2 - dependence) form factor (GeV)

Definition at line 51 of file EmpiricalMECPXSec2015.h.

double genie::EmpiricalMECPXSec2015::fWidth
private

toy model param: width of W distribution (GeV)

Definition at line 53 of file EmpiricalMECPXSec2015.h.

const XSecAlgorithmI* genie::EmpiricalMECPXSec2015::fXSecAlgCCQE
private

cross section algorithm for CCQE

Definition at line 64 of file EmpiricalMECPXSec2015.h.

const XSecAlgorithmI* genie::EmpiricalMECPXSec2015::fXSecAlgEMQE
private

cross section algorithm for EMQE

Definition at line 66 of file EmpiricalMECPXSec2015.h.

const XSecAlgorithmI* genie::EmpiricalMECPXSec2015::fXSecAlgNCQE
private

cross section algorithm for NCQE

Definition at line 65 of file EmpiricalMECPXSec2015.h.

const XSecIntegratorI* genie::EmpiricalMECPXSec2015::fXSecIntegrator
private

Integrator used for reweighting.

Definition at line 69 of file EmpiricalMECPXSec2015.h.


The documentation for this class was generated from the following files: