ReinSehgalSPPXSec.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: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab - March 09, 2006
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14 */
15 //____________________________________________________________________________
16 
17 #include <TMath.h>
18 
27 #include "Framework/Utils/Cache.h"
30 
31 using namespace genie;
32 using namespace genie::constants;
33 
34 //____________________________________________________________________________
36 ReinSehgalRESXSecWithCache("genie::ReinSehgalSPPXSec")
37 {
38 
39 }
40 //____________________________________________________________________________
42 ReinSehgalRESXSecWithCache("genie::ReinSehgalSPPXSec", config)
43 {
44 
45 }
46 //____________________________________________________________________________
48 {
49 
50 }
51 //____________________________________________________________________________
53  const XSecAlgorithmI * model, const Interaction * interaction) const
54 {
55  if(! model->ValidProcess(interaction) ) return 0.;
56 
57  const KPhaseSpace & kps = interaction->PhaseSpace();
58  if(!kps.IsAboveThreshold()) {
59  LOG("COHXSec", pDEBUG) << "*** Below energy threshold";
60  return 0;
61  }
62 
64 
65  //-- Get 1pi exclusive channel
66  SppChannel_t spp_channel = SppChannel::FromInteraction(interaction);
67 
68  //-- Get cache
69  Cache * cache = Cache::Instance();
70 
71  const InitialState & init_state = interaction->InitState();
72  const ProcessInfo & proc_info = interaction->ProcInfo();
73  const Target & target = init_state.Tgt();
74 
76  int nucleon_pdgc = target.HitNucPdg();
77  int nu_pdgc = init_state.ProbePdg();
78 
79  // Get neutrino energy in the struck nucleon rest frame
80  double Ev = init_state.ProbeE(kRfHitNucRest);
81 
82  double xsec = 0;
83 
84  unsigned int nres = fResList.NResonances();
85  for(unsigned int ires = 0; ires < nres; ires++) {
86 
87  //-- Get next resonance from the resonance list
88  Resonance_t res = fResList.ResonanceId(ires);
89 
90  //-- Build a unique name for the cache branch
91  string key = this->CacheBranchName(res, it, nu_pdgc, nucleon_pdgc);
92  LOG("ReinSehgalSpp", pINFO)
93  << "Finding cache branch with key: " << key;
94  CacheBranchFx * cache_branch =
95  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
96 
97  if(!cache_branch) {
98  LOG("ReinSehgalSpp", pWARN)
99  << "No cached RES v-production data for input neutrino"
100  << " (pdgc: " << nu_pdgc << ")";
101  LOG("ReinSehgalSpp", pWARN)
102  << "Wait while computing/caching RES production xsec first...";
103 
104  this->CacheResExcitationXSec(interaction);
105 
106  LOG("ReinSehgalSpp", pINFO) << "Done caching resonance xsec data";
107  LOG("ReinSehgalSpp", pINFO)
108  << "Finding newly created cache branch with key: " << key;
109  cache_branch =
110  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
111  assert(cache_branch);
112  }
113  const CacheBranchFx & cbranch = (*cache_branch);
114 
115  //-- Get cached resonance neutrinoproduction xsec
116  // (If E>Emax, assume xsec = xsec(Emax) - but do not evaluate the
117  // cross section spline at the end of its energy range-)
118  double rxsec = (Ev<fEMax-1) ? cbranch(Ev) : cbranch(fEMax-1);
119 
120  //-- Get the BR for the (resonance) -> (exclusive final state)
121  double br = SppChannel::BranchingRatio(spp_channel, res);
122 
123  //-- Get the Isospin Clebsch-Gordon coefficient for the given resonance
124  // and exclusive final state
125  double igg = SppChannel::IsospinWeight(spp_channel, res);
126 
127  //-- Compute the weighted xsec
128  // (total weight = Breit-Wigner * BR * isospin Clebsch-Gordon)
129  double res_xsec_contrib = rxsec*br*igg;
130 
131  SLOG("ReinSehgalSpp", pINFO)
132  << "Contrib. from [" << utils::res::AsString(res) << "] = "
133  << "<Clebsch-Gordon = " << igg
134  << "> * <BR(->1pi) = " << br
135  << "> * <Breit-Wigner * d^2xsec/dQ^2dW = " << rxsec
136  << "> = " << res_xsec_contrib;
137 
138  //-- Add contribution of this resonance to the cross section
139  xsec += res_xsec_contrib;
140 
141  }//res
142 
143  SLOG("ReinSehgalSpp", pNOTICE)
144  << "XSec[SPP/" << SppChannel::AsString(spp_channel)
145  << "/free] (Ev = " << Ev << " GeV) = " << xsec;
146 
147  //-- If requested return the free nucleon xsec even for input nuclear tgt
148  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
149 
150  //-- number of scattering centers in the target
151  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
152 
153  xsec*=NNucl; // nuclear xsec
154 
155  return xsec;
156 }
157 //____________________________________________________________________________
159 {
160  Algorithm::Configure(config);
161  this->LoadConfig();
162 }
163 //____________________________________________________________________________
165 {
166  Algorithm::Configure(config);
167  this->LoadConfig();
168 }
169 //____________________________________________________________________________
171 {
172  // Get GSL integration type & relative tolerance
173  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
174  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
175  GetParamDef( "gsl-max-eval", fGSLMaxEval, 100000 ) ;
176 
177  // get upper E limit on res xsec spline (=f(E)) before assuming xsec=const
178  GetParamDef( "ESplineMax", fEMax, 100. ) ;
179  fEMax = TMath::Max(fEMax, 20.); // don't accept user Emax if less than 20 GeV
180 
181  // create the baryon resonance list specified in the config.
182  fResList.Clear();
183  string resonances ;
184  GetParam( "ResonanceNameList", resonances ) ;
185  fResList.DecodeFromNameList(resonances);
186 
187 }
188 //____________________________________________________________________________
Cross Section Calculation Interface.
static SppChannel_t FromInteraction(const Interaction *interaction)
Definition: SppChannel.h:277
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
InteractionType_t InteractionTypeId(void) const
string fGSLIntgType
name of GSL numerical integrator
const XML_Char * target
Definition: expat.h:268
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
int HitNucPdg(void) const
Definition: Target.cxx:321
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
void DecodeFromNameList(string list, string delimiter=",")
void Configure(const Registry &config)
Definition: config.py:1
static double IsospinWeight(SppChannel_t channel, Resonance_t res)
Definition: SppChannel.h:201
enum genie::EResonance Resonance_t
unsigned int NResonances(void) const
enum genie::ESppChannel SppChannel_t
static string AsString(SppChannel_t channel)
Definition: SppChannel.h:66
Summary information for an interaction.
Definition: Interaction.h:56
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
string CacheBranchName(Resonance_t r, InteractionType_t it, int nu, int nuc) const
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
Kinematical phase space.
Definition: KPhaseSpace.h:34
int Z(void) const
Definition: Target.h:69
#define pINFO
Definition: Messenger.h:63
void CacheResExcitationXSec(const Interaction *interaction) const
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:89
static double BranchingRatio(SppChannel_t, Resonance_t res)
Definition: SppChannel.h:245
GENIE Cache Memory.
Definition: Cache.h:39
int fGSLMaxEval
GSL max evaluations.
int N(void) const
Definition: Target.h:70
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
assert(nhit_max >=nhit_nbins)
An ABC that caches resonance neutrinoproduction cross sections on free nucleons according to the Rein...
const InitialState & InitState(void) const
Definition: Interaction.h:69
const char * AsString(Resonance_t res)
resonance id -> string
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
#define pNOTICE
Definition: Messenger.h:62
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:67
static Cache * Instance(void)
Definition: Cache.cxx:76
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:38
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
const XML_Char XML_Content * model
Definition: expat.h:151
double ProbeE(RefFrame_t rf) const
enum genie::EInteractionType InteractionType_t
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double fGSLRelTol
required relative tolerance (error)
Resonance_t ResonanceId(unsigned int ires) const