44 using std::ostringstream;
46 using namespace genie;
88 fEventGenList =
"Default";
115 fUnphysEventMask->SetBitNumber(
i,
true);
119 <<
" -> 0) : " << *fUnphysEventMask;
124 if (fUnphysEventMask)
delete fUnphysEventMask;
125 if (fInitState)
delete fInitState;
126 if (fEvGenList)
delete fEvGenList;
127 if (fIntSelector)
delete fIntSelector;
128 if (fIntGenMap)
delete fIntGenMap;
129 if (fXSecSumSpl)
delete fXSecSumSpl;
151 mesg <<
"Configuring event generation driver for initial state: `" 153 <<
"' using event generator list: `" 154 << fEventGenList <<
"'.";
159 this -> BuildInitialState (init_state);
160 this -> BuildGeneratorList ();
161 this -> BuildInteractionGeneratorMap ();
162 this -> BuildInteractionSelector ();
164 LOG(
"GEVGDriver",
pINFO) <<
"Done configuring. \n";
169 LOG(
"GEVGDriver",
pINFO) <<
"Setting the initial state";
171 if(fInitState)
delete fInitState;
174 this->AssertIsValidInitState();
183 <<
"Building the event generator list (specified list name: " 184 << fEventGenList <<
")";
196 <<
"Building the interaction -> generator associations...";
200 fIntGenMap->BuildMap(*fInitState);
202 string mesgh =
"Interaction -> Generator assignments for Initial State: ";
211 LOG(
"GEVGDriver",
pINFO) <<
"Building the interaction selector...";
215 if(fIntSelector)
delete fIntSelector;
217 algf->
AdoptAlgorithm(
"genie::PhysInteractionSelector",
"Default"));
222 *fUnphysEventMask =
mask;
226 <<
" -> 0) : " << *fUnphysEventMask;
232 LOG(
"GEVGDriver",
pINFO) <<
"Creating the initial state";
240 <<
"Selecting an Interaction & Bootstraping the EventRecord";
241 fCurrentRecord = fIntSelector->SelectInteraction(fIntGenMap, nu4p);
243 if(!fCurrentRecord) {
245 <<
"No interaction could be selected for: " 246 << init_state.
AsString() <<
" at E = " << nu4p.E() <<
" GeV";
251 LOG(
"GEVGDriver",
pDEBUG) <<
"Getting the selected interaction";
252 Interaction * interaction = fCurrentRecord->Summary();
263 LOG(
"GEVGDriver",
pINFO) <<
"Finding an appropriate EventGenerator";
280 string mesg =
"Requesting from event generation thread: " +
281 evgen->
Id().
Key() +
" to generate the selected interaction";
286 fCurrentRecord->SetUnphysEventMask(*fUnphysEventMask);
295 bool unphys = fCurrentRecord->IsUnphysical();
297 LOG(
"GEVGDriver",
pINFO) <<
"Returning the current event!";
299 return fCurrentRecord;
301 LOG(
"GEVGDriver",
pWARN) <<
"An unphysical event was generated...";
303 bool accept = fCurrentRecord->Accept();
306 <<
"The generated unphysical event is accepted by the user";
308 return fCurrentRecord;
312 <<
"The generated unphysical event is rejected";
313 delete fCurrentRecord;
319 <<
"Attempting to regenerate the event...";
323 <<
"Could not produce a physical event after " 325 delete fCurrentRecord;
341 <<
"Interaction->Generator Map has not being built yet!";
352 <<
"Setting event generator list: " << listname;
354 fEventGenList = listname;
361 LOG(
"GEVGDriver",
pWARN) <<
"Null interaction!!";
366 <<
"Interaction->Generator Map has not being built yet!";
379 LOG(
"GEVGDriver",
pDEBUG) <<
"Computing the cross section sum";
391 InteractionList::const_iterator intliter;
392 for(intliter = ilst.begin(); intliter != ilst.end(); ++intliter) {
400 <<
"Compute cross section for interaction: \n" << code;
404 fIntGenMap->FindGenerator(interaction)->CrossSectionAlg();
409 bool spline_exists = xssl->
SplineExists(xsec_alg, interaction);
410 if (spline_exists && fUseSplines) {
411 double E = nup4.Energy();
414 xsec = xsec_alg->
Integral(interaction);
416 xsec = TMath::Max(0., xsec);
421 <<
"\nInteraction = " << code
422 <<
"\nCross Section " 423 << (fUseSplines ?
"*interpolated*" :
"*computed*")
432 << pdglib->
Find(fInitState->ProbePdg())->
GetName() <<
"+" 433 << pdglib->
Find(fInitState->Tgt().Pdg())->
GetName() <<
"->X, " 434 <<
"E = " << nup4.Energy() <<
" GeV)" 435 << (fUseSplines ?
"*interpolated*" :
"*computed*")
442 int nk,
double Emin,
double Emax,
bool inlogE)
451 <<
"Creating spline (sum-xsec = f(" << ((inlogE) ?
"logE" :
"E")
452 <<
") in E = [" << Emin <<
", " << Emax <<
"] using " << nk <<
" knots";
455 LOG(
"GEVGDriver",
pFATAL) <<
"You haven't loaded any splines!! ";
458 assert(Emin<Emax && Emin>0 && nk>2);
460 double logEmin=0, logEmax=0,
dE=0;
462 double *
E =
new double[nk];
463 double *
xsec =
new double[nk];
466 logEmin = TMath::Log(Emin);
467 logEmax = TMath::Log(Emax);
468 dE = (logEmax-logEmin)/(nk-1);
470 dE = (Emax-Emin)/(nk-1);
473 TLorentzVector
p4(0,0,0,0);
475 for(
int i=0;
i<nk;
i++) {
476 double e = (inlogE) ? TMath::Exp(logEmin +
i*
dE) : Emin +
i*
dE;
477 p4.SetPxPyPzE(0.,0.,e,e);
478 double xs = this->XSecSum(p4);
483 if (fXSecSumSpl)
delete fXSecSumSpl;
484 fXSecSumSpl =
new Spline(nk, E, xsec);
494 if (!fUseSplines)
return 0;
502 fIntGenMap->FindGenerator(interaction)->CrossSectionAlg();
539 InteractionList::const_iterator intliter;
540 for(intliter = ilst.begin(); intliter != ilst.end(); ++intliter) {
547 fIntGenMap->FindGenerator(interaction)->CrossSectionAlg();
551 bool spl_exists = xsl->
SplineExists(xsec_alg, interaction);
554 fUseSplines = fUseSplines && spl_exists;
559 <<
"Null cross-section algorithm! Can not load cross-section spline.";
564 <<
"Null interaction! Can not load cross-section spline.";
568 <<
"*** At least a spline (algorithm: " 569 << xsec_alg->
Id().
Key() <<
", interaction: " 570 << interaction->
AsString() <<
" doesn't exist. " 571 <<
"Reverting back to not using splines";
585 <<
"Creating (missing) splines with [UseLogE: " 586 << ((useLogE) ?
"ON]" :
"OFF]");
591 EventGeneratorList::const_iterator evgliter;
592 InteractionList::iterator intliter;
595 for(evgliter = fEvGenList->begin();
596 evgliter != fEvGenList->end(); ++evgliter) {
600 <<
"Querying [" << evgen->
Id().
Key()
601 <<
"] for its InteractionList";
623 <<
"Refusing to exceed validity range: Emax = " <<
Emax;
625 emax = TMath::Min(emax,Emax);
635 nknots = (
int) (15 * TMath::Log10(emax-Emin));
637 nknots = TMath::Max(nknots,30);
641 for(intliter = ilst->begin(); intliter != ilst->end(); ++intliter) {
647 SLOG(
"GEVGDriver",
pINFO) <<
"Need xsec spline for " << code;
653 <<
"The spline wasn't loaded at initialization. " 654 <<
"I can build it now but it might take a while...";
655 xsl->
CreateSpline(alg, interaction, nknots, Emin, emax);
657 SLOG(
"GEVGDriver",
pDEBUG) <<
"Spline was found";
679 EventGeneratorList::const_iterator evgliter;
682 for(evgliter = fEvGenList->begin();
683 evgliter != fEvGenList->end(); ++evgliter) {
692 E.
min = TMath::Min(E.
min, Emin);
693 E.
max = TMath::Max(E.
max, Emax);
702 int ppdgc = fInitState->ProbePdg();
710 <<
"\n\n *********************** GEVGDriver ***************************";
712 int ppdg = fInitState->ProbePdg();
713 int tgtpdg = fInitState->Tgt().Pdg();
715 stream <<
"\n |---o Probe PDG-code ......: " << ppdg;
716 stream <<
"\n |---o Target PDG-code .....: " << tgtpdg;
718 stream <<
"\n |---o Using cross section splines is turned " 720 stream <<
"\n |---o Unphysical event filter mask (" 723 stream <<
"\n *********************************************************\n";
Cross Section Calculation Interface.
void Print(ostream &stream) const
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
void SetProbeP4(const TLorentzVector &P4)
THE MAIN GENIE PROJECT NAMESPACE
virtual void ProcessEventRecord(GHepRecord *event_rec) const =0
EventGeneratorList * AssembleGeneratorList()
void CreateSpline(const XSecAlgorithmI *alg, const Interaction *i, int nknots=-1, double e_min=-1, double e_max=-1)
Defines the InteractionListGeneratorI interface. Concrete implementations of this interface generate ...
virtual const InteractionListGeneratorI * IntListGenerator(void) const =0
A simple [min,max] interval for doubles.
const EventGeneratorI * FindGenerator(const Interaction *interaction) const
bool SplineExists(const XSecAlgorithmI *alg, const Interaction *i) const
A numeric analysis tool class for interpolating 1-D functions.
Defines the EventGeneratorI interface.
bool IsDarkMatter(int pdgc)
void BuildInteractionGeneratorMap(void)
static XSecSplineList * Instance()
double Evaluate(double x) const
string AsString(void) const
virtual const GVldContext & ValidityContext(void) const =0
Range1D_t ValidEnergyRange(void) const
string BoolAsIOString(bool b)
void GenerateEvent(string mesg)
Summary information for an interaction.
An Interaction -> EventGeneratorI associative container. The container is being built for the loaded ...
std::string GetName(int i)
static unsigned int NFlags(void)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void SetEventGeneratorList(string listname)
const InteractionList * Interactions(void) const
double XSecSum(const TLorentzVector &nup4)
virtual double Integral(const Interaction *i) const =0
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
string AsString(void) const
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Algorithm * AdoptAlgorithm(const AlgId &algid) const
virtual InteractionList * CreateInteractionList(const InitialState &init) const =0
Misc GENIE control constants.
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
void BuildGeneratorList(void)
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Defines the InteractionSelectorI interface to be implemented by algorithms selecting interactions to ...
static RunningThreadInfo * Instance(void)
void BuildInteractionSelector(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
static PDGLibrary * Instance(void)
static AlgFactory * Instance()
const Spline * XSecSpline(const Interaction *interaction) const
Singleton class to load & serve a TDatabasePDG.
ostream & operator<<(ostream &stream, const AlgConfigPool &config_pool)
void UpdateRunningThread(const EventGeneratorI *evg)
void Configure(string mesg)
TParticlePDG * Find(int pdgc)
void SetLogE(bool on)
set opt to build splines as f(E) or as f(logE)
void Configure(int nu_pdgc, int Z, int A)
assert(nhit_max >=nhit_nbins)
A vector of Interaction objects.
InitialState * InitStatePtr(void) const
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
void CreateSplines(int nknots=-1, double emax=-1, bool inLogE=true)
The GENIE Algorithm Factory.
void UseGeneratorList(const EventGeneratorList *list)
void SetUnphysEventMask(const TBits &mask)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
const Spline * GetSpline(const XSecAlgorithmI *alg, const Interaction *i) const
List of cross section vs energy splines.
Module to generate only pions from cosmic rays.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Keep info on the event generation thread currently on charge. This is used so that event generation m...
void AssertIsValidInitState(void) const
Assembles a list of all the EventGeneratorI subclasses that can be employed during a neutrino event g...
Initial State information.
EventRecord * GenerateEvent(const TLorentzVector &nu4p)
void BuildInitialState(const InitialState &init_state)
static const unsigned int kRecursiveModeMaxDepth