RSPPResonanceSelector.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
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Mar 03, 2009 - CA
14  Moved into the new RES package from its previous location (EVGModules).
15  @ Jul 23, 2010 - CA
16  Use ResonanceCharge() from base class. Function removed from utils::res.
17 
18 */
19 //____________________________________________________________________________
20 
21 #include <vector>
22 #include <sstream>
23 
37 
38 using std::vector;
39 using std::ostringstream;
40 
41 using namespace genie;
42 
43 //___________________________________________________________________________
45 HadronicSystemGenerator("genie::RSPPResonanceSelector")
46 {
47 
48 }
49 //___________________________________________________________________________
51 HadronicSystemGenerator("genie::RSPPResonanceSelector", config)
52 {
53 
54 }
55 //___________________________________________________________________________
57 {
58 
59 }
60 //___________________________________________________________________________
62 {
63  //-- select a baryon resonance
64  Resonance_t res = this->SelectResonance(evrec);
65  assert(res != kNoResonance);
66 
67  //-- add the resonance at the event summary
68  Interaction * interaction = evrec->Summary();
69  interaction->ExclTagPtr()->SetResonance(res);
70 
71  //-- add an entry at the GHep event record & the event summary
72  this->AddResonance(evrec);
73 }
74 //___________________________________________________________________________
76 {
77 // Select a baryon resonance from a list of considered resonances based on
78 // their differential d^2xsec/dWdQ^2 cross sections
79 
80  LOG("RESSelector", pNOTICE) << "Selecting a baryon resonance";
81 
82  Interaction * interaction = evrec->Summary();
83 
84  //-- Figure out what the resonance charge should be.
85  int q_res = this->ResonanceCharge(evrec);
86 
87  //-- Use selected kinematics
88  interaction->KinePtr()->UseSelectedKinematics();
89 
90  //-- Trust kinematics and process type already set.
91  interaction->SetBit(kISkipProcessChk);
92  interaction->SetBit(kISkipKinematicChk);
93 
94  //-- Access cross section algorithm for running thread
95  // The algorithm must be able to compute the RES contribution to
96  // the specified RES/SPP channel
98  const EventGeneratorI * evg = rtinfo->RunningThread();
99  const XSecAlgorithmI * xsecalg = evg->CrossSectionAlg();
100 
101  //-- Loop over all considered baryon resonances and compute the double
102  // differential cross section for the selected kinematical variables
103 
104  double xsec_sum = 0;
105  unsigned int nres = fResList.NResonances();
106  vector<double> xsec_vec(nres);
107 
108  for(unsigned int ires = 0; ires < nres; ires++) {
109 
110  //-- Current resonance
111  Resonance_t res = fResList.ResonanceId(ires);
112 
113  //-- Set the current resonance at the interaction summary
114  // compute the differential cross section d^2xsec/dWdQ^2
115  // (do it only for resonances that can conserve charge)
116  interaction->ExclTagPtr()->SetResonance(res);
117 
118  double xsec = 0;
119  bool skip = (q_res==2 && !utils::res::IsDelta(res));
120 
121  if(!skip) xsec = xsecalg->XSec(interaction,kPSWQ2fE);
122  else {
123  SLOG("RESSelector", pNOTICE)
124  << "RES: " << utils::res::AsString(res)
125  << " would not conserve charge -- skipping it";
126  }
127  //-- For the ith resonance store the sum of (xsec) * (breit-wigner)
128  // for the resonances in the range [0,i]
129  xsec_sum += xsec;
130  xsec_vec[ires] = xsec_sum;
131 
132  SLOG("RESSelector", pNOTICE)
133  << "Resonances (0->" << ires << "): "
134  << "Sum{ BW(W) * d^2xsec(E,W,Q^2)/dWd*Q^2 } = " << xsec_sum;
135  } // res
136 
137  //-- Reset 'trust' bits
138  interaction->ResetBit(kISkipProcessChk);
139  interaction->ResetBit(kISkipKinematicChk);
140 
141  //-- Reset running kinematics
142  interaction->KinePtr()->ClearRunningValues();
143 
144  //-- Use the computed differential cross sections to select a resonance
146  double R = xsec_sum * rnd->RndGen().Rndm();
147 
148  SLOG("RESSelector", pDEBUG) << "R = " << R;
149 
150  for(unsigned int ires = 0; ires < nres; ires++) {
151  SLOG("RESSelector", pDEBUG)
152  << "SUM-XSEC(0->" << ires <<") = " << xsec_vec[ires];
153 
154  if(R < xsec_vec[ires]) {
155  Resonance_t sres = fResList.ResonanceId(ires); // selected RES.
156  LOG("RESSelector", pNOTICE)
157  << "Selected RES = " << utils::res::AsString(sres);
158  return sres;
159  }
160  }
161  LOG("RESSelector", pERROR) << "** Failed to select a resonance";
162  return kNoResonance;
163 }
164 //___________________________________________________________________________
166 {
167  // compute RES p4 = p4(neutrino) + p4(hit nucleon) - p4(primary lepton)
168  TLorentzVector p4 = this->Hadronic4pLAB(evrec);
169 
170  //-- Determine the RES pdg code (from the selected Resonance_t & charge)
171  Interaction * interaction = evrec->Summary();
172  Resonance_t res = interaction->ExclTag().Resonance();
173  int charge = this->ResonanceCharge(evrec);
174  int pdgc = utils::res::PdgCode(res,charge);
175 
176  LOG("RESSelector", pNOTICE)
177  << "Adding RES with PDGC = " << pdgc << ", Q = " << charge;
178 
179  //-- Add the resonance at the EventRecord
181  int mom = evrec->HitNucleonPosition();
182 
183  evrec->AddParticle(
184  pdgc, ist, mom,-1,-1,-1, p4.Px(),p4.Py(),p4.Pz(),p4.E(), 0,0,0,0);
185 }
186 //___________________________________________________________________________
188 {
189  Algorithm::Configure(config);
190  this->LoadConfig();
191 }
192 //____________________________________________________________________________
193 void RSPPResonanceSelector::Configure(string param_set)
194 {
195  Algorithm::Configure(param_set);
196  this->LoadConfig();
197 }
198 //____________________________________________________________________________
200 {
201  fResList.Clear();
202  string resonances = "";
203  this->GetParam("ResonanceNameList", resonances);
204  SLOG("RESSelector", pDEBUG) << "Resonance list: " << resonances;
205 
206  fResList.DecodeFromNameList(resonances);
207  LOG("RESSelector", pINFO) << fResList;
208 }
209 //____________________________________________________________________________
bool IsDelta(Resonance_t res)
is it a Delta resonance?
Cross Section Calculation Interface.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
Resonance_t SelectResonance(GHepRecord *event_rec) const
TLorentzVector Hadronic4pLAB(GHepRecord *event_rec) const
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
enum genie::EGHepStatus GHepStatus_t
void DecodeFromNameList(string list, string delimiter=",")
Defines the EventGeneratorI interface.
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:454
Definition: config.py:1
enum genie::EResonance Resonance_t
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
unsigned int NResonances(void) const
void SetResonance(Resonance_t res)
Definition: XclsTag.cxx:123
Summary information for an interaction.
Definition: Interaction.h:56
void UseSelectedKinematics(void)
Definition: Kinematics.cxx:369
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void ProcessEventRecord(GHepRecord *event_rec) const
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
#define R(x)
#define pINFO
Definition: Messenger.h:63
Resonance_t Resonance(void) const
Definition: XclsTag.h:62
Double_t xsec[nknots]
Definition: testXsec.C:47
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
XclsTag * ExclTagPtr(void) const
Definition: Interaction.h:77
static RunningThreadInfo * Instance(void)
void AddResonance(GHepRecord *event_rec) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
int ResonanceCharge(GHepRecord *event_rec) const
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:81
assert(nhit_max >=nhit_nbins)
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
const char * AsString(Resonance_t res)
resonance id -> string
Abstract class. Is used to pass some commonly recurring methods to all concrete implementations of th...
#define pNOTICE
Definition: Messenger.h:62
void ClearRunningValues(void)
Definition: Kinematics.cxx:357
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const EventGeneratorI * RunningThread(void)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
Keep info on the event generation thread currently on charge. This is used so that event generation m...
BaryonResList fResList
baryon resonances taken into account
void Configure(const Registry &config)
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
#define pDEBUG
Definition: Messenger.h:64
Resonance_t ResonanceId(unsigned int ires) const