AlvarezRusoCOHPiPXSec.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: Steve Dennis
8  University of Warwick, Rutherford Appleton Laboratory,
9  October 5, 2012
10 
11  For the class documentation see the corresponding header file.
12 
13 */
14 //____________________________________________________________________________
15 
16 #include <iostream>
17 
18 #include <TMath.h>
19 
32 
37 
38 
39 using namespace genie;
40 using namespace genie::constants;
41 using namespace genie::utils;
42 
43 using namespace alvarezruso;
44 
45 //____________________________________________________________________________
47 XSecAlgorithmI("genie::AlvarezRusoCOHPiPXSec")
48 {
49  fMultidiff = NULL;
50  fLastInteraction = NULL;
51 }
52 //____________________________________________________________________________
54 XSecAlgorithmI("genie::AlvarezRusoCOHPiPXSec", config)
55 {
56  fMultidiff = NULL;
57  fLastInteraction = NULL;
58 }
59 //____________________________________________________________________________
61 {
62  if (fMultidiff) delete fMultidiff;
63 }
64 //____________________________________________________________________________
66  const Interaction * interaction, KinePhaseSpace_t kps) const
67 {
68  if(! this -> ValidProcess (interaction) ) return 0.;
69  if(! this -> ValidKinematics (interaction) ) return 0.;
70 
71  const Kinematics & kinematics = interaction -> Kine();
72  const InitialState & init_state = interaction -> InitState();
73 
74  int A = init_state.Tgt().A(); // mass number
75  int Z = init_state.Tgt().Z(); // atomic number
76  double E_nu = init_state.ProbeE(kRfLab); // neutrino energy
77 
78  const TLorentzVector p4_lep = kinematics.FSLeptonP4();
79  const TLorentzVector p4_pi = kinematics.HadSystP4();
80  double E_lep = p4_lep.E();
81 
82  if (fLastInteraction!=interaction) {
83  if (fMultidiff != NULL) {
84  delete fMultidiff;
85  fMultidiff = NULL;
86  }
87 
88  current_t current;
89  if ( interaction->ProcInfo().IsWeakCC() ) {
90  current = kCC;
91  }
92  else if ( interaction->ProcInfo().IsWeakNC() ) {
93  current = kNC;
94  }
95  else {
96  LOG("AlvarezRusoCohPi",pDEBUG)<<"Unknown current for AlvarezRuso implementation";
97  return 0.;
98  }
99 
100  flavour_t flavour;
101  if ( init_state.ProbePdg() == 12 || init_state.ProbePdg() == -12) {
102  flavour=kE;
103  }
104  else if ( init_state.ProbePdg() == 14 || init_state.ProbePdg() == -14) {
105  flavour=kMu;
106  }
107  else if ( init_state.ProbePdg() == 16 || init_state.ProbePdg() == -16) {
108  flavour=kTau;
109  }
110  else {
111  LOG("AlvarezRusoCohPi",pDEBUG)<<"Unknown probe for AlvarezRuso implementation";
112  return 0.;
113  }
114 
116  if ( init_state.ProbePdg() > 0) {
117  nutype = kNu;
118  } else {
119  nutype = kAntiNu;
120  }
121 
122  fMultidiff = new AlvarezRusoCOHPiPDXSec(Z, A ,current, flavour, nutype);
123  fLastInteraction = interaction;
124  }
125 
126  double xsec = fMultidiff->DXSec(E_nu, E_lep, p4_lep.Theta(), p4_lep.Phi(), p4_pi.Theta(), p4_pi.Phi());
127  xsec = xsec * 1E-38 * units::cm2;
128 
129  if (kps != kPSElOlOpifE) {
130  xsec *= utils::kinematics::Jacobian(interaction, kPSElOlOpifE, kps );
131  }
132 
133  return (xsec);
134 }
135 //____________________________________________________________________________
136 double AlvarezRusoCOHPiPXSec::Integral(const Interaction * interaction) const
137 {
138  double xsec = fXSecIntegrator->Integrate(this,interaction);
139  return xsec;
140 }
141 //____________________________________________________________________________
142 bool AlvarezRusoCOHPiPXSec::ValidProcess(const Interaction * interaction) const
143 {
144  if(interaction->TestBit(kISkipProcessChk)) return true;
145 
146  const InitialState & init_state = interaction->InitState();
147  const ProcessInfo & proc_info = interaction->ProcInfo();
148  const Target & target = init_state.Tgt();
149 
150  int nu = init_state.ProbePdg();
151 
152  if (!proc_info.IsCoherent()) return false;
153  if (!proc_info.IsWeak()) return false;
154  if (target.HitNucIsSet()) return false;
155  if (!(target.A()>1)) return false;
156  if (!pdg::IsNeutrino(nu) && !pdg::IsAntiNeutrino(nu)) return false;
157 
158  return true;
159 }
160 //____________________________________________________________________________
162 {
163  Algorithm::Configure(config);
164  this->LoadConfig();
165 }
166 //____________________________________________________________________________
168 {
169  Algorithm::Configure(config);
170  this->LoadConfig();
171 }
172 //____________________________________________________________________________
174 {
175  //AlgConfigPool * confp = AlgConfigPool::Instance();
176  //const Registry * gc = confp->GlobalParameterList();
177 
178  /*fRo = fConfig->GetDoubleDef("COH-Ro", gc->GetDouble("COH-Ro"));
179  fMa = fConfig->GetDoubleDef("Ma", gc->GetDouble("COH-Ma"));
180  fa4 = fConfig->GetDoubleDef("a4", gc->GetDouble("COHAR-a4"));
181  fa5 = fConfig->GetDoubleDef("a5", gc->GetDouble("COHAR-a5"));
182  fb4 = fConfig->GetDoubleDef("b4", gc->GetDouble("COHAR-b4"));
183  fb5 = fConfig->GetDoubleDef("b5", gc->GetDouble("COHAR-b5"));
184  ffPi = fConfig->GetDoubleDef("fPi", gc->GetDouble("COHAR-fPi"));
185  ffStar = fConfig->GetDoubleDef("fStar", gc->GetDouble("COHAR-fStar"));*/
186 
187 
188  //-- load the differential cross section integrator
190  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
192 
193 }
194 //____________________________________________________________________________
195 
Cross Section Calculation Interface.
Basic constants.
bool IsWeak(void) const
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
const XML_Char * target
Definition: expat.h:268
Antineutrinos-only.
Definition: IPrediction.h:50
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
int A(void) const
Definition: Target.h:71
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double Integral(const Interaction *i) const
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:67
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
void Configure(const Registry &config)
static const double cm2
Definition: Units.h:77
Float_t Z
Definition: plot.C:38
unsigned int nutype(unsigned int u)
Definition: runWimpSim.h:122
Summary information for an interaction.
Definition: Interaction.h:56
Charged-current interactions.
Definition: IPrediction.h:39
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
bool IsWeakNC(void) const
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:66
double DXSec(const double E_nu_, const double E_l_, const double theta_l_, const double phi_l_, const double theta_pi_, const double phi_pi_)
#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
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
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
int Z(void) const
Definition: Target.h:69
Double_t xsec[nknots]
Definition: testXsec.C:47
static const double A
Definition: Units.h:82
bool HitNucIsSet(void) const
Definition: Target.cxx:300
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
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
const Target & Tgt(void) const
Definition: InitialState.h:67
const XSecIntegratorI * fXSecIntegrator
void kinematics()
Definition: kinematics.C:10
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
const Interaction * fLastInteraction
Root of GENIE utility namespaces.
alvarezruso::AlvarezRusoCOHPiPDXSec * fMultidiff
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
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353