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