GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PhysInteractionSelector.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 <vector>
12 #include <sstream>
13 #include <cstdlib>
14 #include <iomanip>
15 
16 #include <TMath.h>
17 #include <TLorentzVector.h>
18 
31 
32 using std::vector;
33 using std::endl;
34 using std::setw;
35 using std::setprecision;
36 using std::setfill;
37 using std::ostringstream;
38 using namespace genie;
39 //using namespace genie::units;
40 
41 //___________________________________________________________________________
43 InteractionSelectorI("genie::PhysInteractionSelector")
44 {
45 
46 }
47 //___________________________________________________________________________
49 InteractionSelectorI("genie::PhysInteractionSelector", config)
50 {
51 
52 }
53 //___________________________________________________________________________
55 {
56 
57 }
58 //___________________________________________________________________________
60  (const InteractionGeneratorMap * igmap, const TLorentzVector & p4) const
61 {
62  if(!igmap) {
63  LOG("IntSel", pERROR)
64  << "\n*** NULL InteractionGeneratorMap! Can't select interaction";
65  return 0;
66  }
67  if(igmap->size() <= 0) {
68  LOG("IntSel", pERROR)
69  << "\n*** Empty InteractionGeneratorMap! Can't select interaction";
70  return 0;
71  }
72 
73  // Get the list of spline objects
74  // Should have been constructed at the job initialization
75  XSecSplineList * xssl = 0;
76  if (fUseSplines) xssl = XSecSplineList::Instance();
77 
78  const InteractionList & ilst = igmap->GetInteractionList();
79  vector<double> xseclist(ilst.size());
80 
81  string istate = ilst[0]->InitState().AsString();
82  ostringstream msg;
83  msg << "Selecting an interaction for the given initial state = "
84  << istate << " at E = " << p4.E() << " GeV";
85 
86  LOG("IntSel", pNOTICE)
87  << utils::print::PrintFramedMesg(msg.str(), 0, '=');
88  LOG("IntSel", pNOTICE)
89  << "Computing xsecs for all relevant modeled interactions:";
90 
91  unsigned int i=0;
92  InteractionList::const_iterator intliter = ilst.begin();
93 
94  ostringstream xsec_table_printout;
95 
96  xsec_table_printout
97  << " |" << setfill('-') << setw(112) << "|" << endl
98  << " | " << setfill(' ') << setw(80) << "interaction"
99  << " | cross-section (1E-38*cm^2) |" << endl
100  << " |" << setfill('-') << setw(112) << "|" << endl;
101 
102  for( ; intliter != ilst.end(); ++intliter) {
103 
104  Interaction * interaction = new Interaction(**intliter);
105  interaction->InitStatePtr()->SetProbeP4(p4);
106 
107  SLOG("IntSel", pDEBUG)
108  << "Computing xsec for: \n " << interaction->AsString();
109 
110  // get the cross section for this interaction
111  const XSecAlgorithmI * xsec_alg =
112  igmap->FindGenerator(interaction)->CrossSectionAlg();
113  assert(xsec_alg);
114 
115  double xsec = 0; // cross section for this interaction
116 
117  bool spline_computed = xssl->SplineExists(xsec_alg, interaction);
118  bool eval = fUseSplines && spline_computed;
119  if (eval) {
120  const InitialState & init = interaction->InitState();
121  const ProcessInfo & proc = interaction->ProcInfo();
122  double E = init.ProbeE(kRfLab);
123  if(TMath::IsNaN(E)) {
124  BLOG("IntSel", pFATAL) << *interaction;
125  BLOG("IntSel", pFATAL) << "E = " << E;
126  abort();
127  }
128  const Spline * spl = xssl->GetSpline(xsec_alg,interaction);
129  if(spl->ClosestKnotValueIsZero(E,"-")) xsec = 0;
130  else xsec = spl->Evaluate(E);
131  } else {
132  xsec = xsec_alg->Integral(interaction);
133  }
134  TMath::Max(0., xsec);
135 /*
136  LOG("IntSel", pNOTICE)
137  << interaction->AsString()
138  << " --> xsec " << (eval ? "[**interp**]" : "[**calc**]")
139  << " = " << xsec/genie::units::cm2 << " cm^2";
140 */
141  xsec_table_printout
142  << " | " << setfill(' ') << setw(80) << interaction->AsString()
143  << " | " << setfill(' ') << setw(26) << xsec/(1E-38*genie::units::cm2)
144  << " | " << endl;
145 
146  xseclist[i++] = xsec;
147  delete interaction;
148 
149  } // loop over interaction that can be generated
150 
151  xsec_table_printout
152  << " |" << setfill('-') << setw(112) << "|" << endl;
153 
154  LOG("IntSel", pNOTICE)
155  << "\n" << xsec_table_printout.str();
156 
157  // select an interaction
158 
159  LOG("IntSel", pINFO)
160  << "Selecting an entry from the Interaction List";
161  double xsec_sum = 0;
162  for(unsigned int iint = 0; iint < xseclist.size(); iint++) {
163  xsec_sum += xseclist[iint];
164  xseclist[iint] = xsec_sum;
165 
166  SLOG("IntSel", pINFO)
167  << "Sum{xsec}(0->" << iint << ") = " << xsec_sum;
168  }
169  RandomGen * rnd = RandomGen::Instance();
170  double R = xsec_sum * rnd->RndISel().Rndm();
171 
172  LOG("IntSel", pINFO)
173  << "Generating Rndm (0. -> max = " << xsec_sum << ") = " << R;
174 
175  for(unsigned int iint = 0; iint < xseclist.size(); iint++) {
176 
177  SLOG("IntSel", pDEBUG)
178  << "Sum{xsec}(0->" << iint <<") = " << xseclist[iint];
179 
180  if( R < xseclist[iint] ) {
181  Interaction * selected_interaction = new Interaction (*ilst[iint]);
182  selected_interaction->InitStatePtr()->SetProbeP4(p4);
183 
184  // set the cross section for the selected interaction (just extract it
185  // from the array of summed xsecs rather than recomputing it)
186  double xsec_pedestal = (iint > 0) ? xseclist[iint-1] : 0.;
187  double xsec = xseclist[iint] - xsec_pedestal;
188  assert(xsec>0);
189 
190  LOG("IntSel", pNOTICE)
191  << "Selected interaction: " << selected_interaction->AsString();
192 
193  // bootstrap the event record
194  EventRecord * evrec = new EventRecord;
195  evrec->AttachSummary(selected_interaction);
196  evrec->SetXSec(xsec);
197 
198  return evrec;
199  }
200  }
201  LOG("IntSel", pERROR) << "Could not select interaction";
202  return 0;
203 }
204 //___________________________________________________________________________
206 {
207  Algorithm::Configure(config);
208  this->LoadConfigData();
209 }
210 //____________________________________________________________________________
212 {
213  Algorithm::Configure(param_set);
214  this->LoadConfigData();
215 }
216 //____________________________________________________________________________
218 {
219  //check whether the user prefers the cross sections to be calculated or
220  //evaluated from a spline object constructed at the job initialization
221  fUseSplines = false ;
222  GetParam( "UseStoredXSecs", fUseSplines ) ;
223 
224 }
225 //___________________________________________________________________________
Cross Section Calculation Interface.
virtual void SetXSec(double xsec)
Definition: GHepRecord.h:132
EventRecord * SelectInteraction(const InteractionGeneratorMap *igmp, const TLorentzVector &p4) const
implement the InteractionSelectorI interface
void SetProbeP4(const TLorentzVector &P4)
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
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:58
#define pFATAL
Definition: Messenger.h:56
static XSecSplineList * Instance()
double Evaluate(double x) const
Definition: Spline.cxx:363
#define BLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending no additional information.
Definition: Messenger.h:212
bool ClosestKnotValueIsZero(double x, Option_t *opt="-+") const
Definition: Spline.cxx:561
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
string AsString(void) const
virtual void AttachSummary(Interaction *interaction)
Definition: GHepRecord.cxx:99
Summary information for an interaction.
Definition: Interaction.h:56
An Interaction -&gt; 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:96
static constexpr double cm2
Definition: Units.h:69
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
#define pINFO
Definition: Messenger.h:62
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
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:65
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
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:164
TRandom3 & RndISel(void) const
rnd number generator used by interaction selectors
Definition: RandomGen.h:65
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
#define pNOTICE
Definition: Messenger.h:61
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:84
const Spline * GetSpline(const XSecAlgorithmI *alg, const Interaction *i) const
List of cross section vs energy splines.
double ProbeE(RefFrame_t rf) const
virtual double Integral(const Interaction *i) const =0
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63