PhysInteractionSelector.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 - January 25, 2005
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Dec 03, 2007 - CA
14  Setting the reference frame explicitly before computing the energy passed
15  to the spline objects rather than depending on the interaction doing the
16  right thing. Protect cross section eval. against NaN input as TSpline3
17  itself doesn't and its hard to diagnose problems from its actuall err mesg.
18  @ Jun 23, 2008 - CA
19  Protect against round off err / negative xsec
20 */
21 //____________________________________________________________________________
22 
23 #include <vector>
24 #include <sstream>
25 #include <cstdlib>
26 #include <iomanip>
27 
28 #include <TMath.h>
29 #include <TLorentzVector.h>
30 
43 
44 using std::vector;
45 using std::endl;
46 using std::setw;
47 using std::setprecision;
48 using std::setfill;
49 using std::ostringstream;
50 using namespace genie;
51 //using namespace genie::units;
52 
53 //___________________________________________________________________________
55 InteractionSelectorI("genie::PhysInteractionSelector")
56 {
57 
58 }
59 //___________________________________________________________________________
61 InteractionSelectorI("genie::PhysInteractionSelector", config)
62 {
63 
64 }
65 //___________________________________________________________________________
67 {
68 
69 }
70 //___________________________________________________________________________
72  (const InteractionGeneratorMap * igmap, const TLorentzVector & p4) const
73 {
74  if(!igmap) {
75  LOG("IntSel", pERROR)
76  << "\n*** NULL InteractionGeneratorMap! Can't select interaction";
77  return 0;
78  }
79  if(igmap->size() <= 0) {
80  LOG("IntSel", pERROR)
81  << "\n*** Empty InteractionGeneratorMap! Can't select interaction";
82  return 0;
83  }
84 
85  // Get the list of spline objects
86  // Should have been constructed at the job initialization
87  XSecSplineList * xssl = 0;
89 
90  const InteractionList & ilst = igmap->GetInteractionList();
91  vector<double> xseclist(ilst.size());
92 
93  string istate = ilst[0]->InitState().AsString();
94  ostringstream msg;
95  msg << "Selecting an interaction for the given initial state = "
96  << istate << " at E = " << p4.E() << " GeV";
97 
98  LOG("IntSel", pNOTICE)
99  << utils::print::PrintFramedMesg(msg.str(), 0, '=');
100  LOG("IntSel", pNOTICE)
101  << "Computing xsecs for all relevant modeled interactions:";
102 
103  unsigned int i=0;
104  InteractionList::const_iterator intliter = ilst.begin();
105 
106  ostringstream xsec_table_printout;
107 
108  xsec_table_printout
109  << " |" << setfill('-') << setw(112) << "|" << endl
110  << " | " << setfill(' ') << setw(80) << "interaction"
111  << " | cross-section (1E-38*cm^2) |" << endl
112  << " |" << setfill('-') << setw(112) << "|" << endl;
113 
114  for( ; intliter != ilst.end(); ++intliter) {
115 
116  Interaction * interaction = new Interaction(**intliter);
117  interaction->InitStatePtr()->SetProbeP4(p4);
118 
119  SLOG("IntSel", pDEBUG)
120  << "Computing xsec for: \n " << interaction->AsString();
121 
122  // get the cross section for this interaction
123  const XSecAlgorithmI * xsec_alg =
124  igmap->FindGenerator(interaction)->CrossSectionAlg();
125  assert(xsec_alg);
126 
127  double xsec = 0; // cross section for this interaction
128 
129  bool spline_computed = xssl->SplineExists(xsec_alg, interaction);
130  bool eval = fUseSplines && spline_computed;
131  if (eval) {
132  const InitialState & init = interaction->InitState();
133  const ProcessInfo & proc = interaction->ProcInfo();
134  // choose ref frame ('Lab' or 'Hit nucleon rest frame')
135  RefFrame_t frame =
136  (proc.IsCoherent() || proc.IsElectronScattering()) ?
138  double E = init.ProbeE(frame);
139  if(TMath::IsNaN(E)) {
140  BLOG("IntSel", pFATAL) << *interaction;
141  BLOG("IntSel", pFATAL) << "E = " << E;
142  abort();
143  }
144  const Spline * spl = xssl->GetSpline(xsec_alg,interaction);
145  if(spl->ClosestKnotValueIsZero(E,"-")) xsec = 0;
146  else xsec = spl->Evaluate(E);
147  } else {
148  xsec = xsec_alg->Integral(interaction);
149  }
150  TMath::Max(0., xsec);
151 /*
152  LOG("IntSel", pNOTICE)
153  << interaction->AsString()
154  << " --> xsec " << (eval ? "[**interp**]" : "[**calc**]")
155  << " = " << xsec/genie::units::cm2 << " cm^2";
156 */
157  xsec_table_printout
158  << " | " << setfill(' ') << setw(80) << interaction->AsString()
159  << " | " << setfill(' ') << setw(26) << xsec/(1E-38*genie::units::cm2)
160  << " | " << endl;
161 
162  xseclist[i++] = xsec;
163  delete interaction;
164 
165  } // loop over interaction that can be generated
166 
167  xsec_table_printout
168  << " |" << setfill('-') << setw(112) << "|" << endl;
169 
170  LOG("IntSel", pNOTICE)
171  << "\n" << xsec_table_printout.str();
172 
173  // select an interaction
174 
175  LOG("IntSel", pINFO)
176  << "Selecting an entry from the Interaction List";
177  double xsec_sum = 0;
178  for(unsigned int iint = 0; iint < xseclist.size(); iint++) {
179  xsec_sum += xseclist[iint];
180  xseclist[iint] = xsec_sum;
181 
182  SLOG("IntSel", pINFO)
183  << "Sum{xsec}(0->" << iint << ") = " << xsec_sum;
184  }
186  double R = xsec_sum * rnd->RndISel().Rndm();
187 
188  LOG("IntSel", pINFO)
189  << "Generating Rndm (0. -> max = " << xsec_sum << ") = " << R;
190 
191  for(unsigned int iint = 0; iint < xseclist.size(); iint++) {
192 
193  SLOG("IntSel", pDEBUG)
194  << "Sum{xsec}(0->" << iint <<") = " << xseclist[iint];
195 
196  if( R < xseclist[iint] ) {
197  Interaction * selected_interaction = new Interaction (*ilst[iint]);
198  selected_interaction->InitStatePtr()->SetProbeP4(p4);
199 
200  // set the cross section for the selected interaction (just extract it
201  // from the array of summed xsecs rather than recomputing it)
202  double xsec_pedestal = (iint > 0) ? xseclist[iint-1] : 0.;
203  double xsec = xseclist[iint] - xsec_pedestal;
204  assert(xsec>0);
205 
206  LOG("IntSel", pNOTICE)
207  << "Selected interaction: " << selected_interaction->AsString();
208 
209  // bootstrap the event record
210  EventRecord * evrec = new EventRecord;
211  evrec->AttachSummary(selected_interaction);
212  evrec->SetXSec(xsec);
213 
214  return evrec;
215  }
216  }
217  LOG("IntSel", pERROR) << "Could not select interaction";
218  return 0;
219 }
220 //___________________________________________________________________________
222 {
223  Algorithm::Configure(config);
224  this->LoadConfigData();
225 }
226 //____________________________________________________________________________
228 {
229  Algorithm::Configure(param_set);
230  this->LoadConfigData();
231 }
232 //____________________________________________________________________________
234 {
235  //check whether the user prefers the cross sections to be calculated or
236  //evaluated from a spline object constructed at the job initialization
237  fUseSplines = false ;
238  GetParam( "UseStoredXSecs", fUseSplines ) ;
239 
240 }
241 //___________________________________________________________________________
Cross Section Calculation Interface.
virtual void SetXSec(double xsec)
Definition: GHepRecord.h:133
EventRecord * SelectInteraction(const InteractionGeneratorMap *igmp, const TLorentzVector &p4) const
implement the InteractionSelectorI interface
void SetProbeP4(const TLorentzVector &P4)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const EventGeneratorI * FindGenerator(const Interaction *in) const
bool SplineExists(const XSecAlgorithmI *alg, const Interaction *i) const
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
#define pFATAL
Definition: Messenger.h:57
bool IsElectronScattering(void) const
Definition: config.py:1
enum genie::ERefFrame RefFrame_t
static XSecSplineList * Instance()
double Evaluate(double x) const
Definition: Spline.cxx:362
#define BLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending no additional information.
Definition: Messenger.h:215
static const double cm2
Definition: Units.h:77
bool ClosestKnotValueIsZero(double x, Option_t *opt="-+") const
Definition: Spline.cxx:555
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
string AsString(void) const
virtual void AttachSummary(Interaction *interaction)
Definition: GHepRecord.cxx:143
Summary information for an interaction.
Definition: Interaction.h:56
An Interaction -> EventGeneratorI associative container. The container is being built for the loaded ...
void Configure(const Registry &config)
#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
virtual double Integral(const Interaction *i) const =0
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
#define R(x)
#define pINFO
Definition: Messenger.h:63
Double_t xsec[nknots]
Definition: testXsec.C:47
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
Defines the InteractionSelectorI interface to be implemented by algorithms selecting interactions to ...
const InteractionList & GetInteractionList(void) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
assert(nhit_max >=nhit_nbins)
A vector of Interaction objects.
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
TRandom3 & RndISel(void) const
rnd number generator used by interaction selectors
Definition: RandomGen.h:66
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
#define pNOTICE
Definition: Messenger.h:62
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
#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 Spline * GetSpline(const XSecAlgorithmI *alg, const Interaction *i) const
List of cross section vs energy splines.
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97