GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RSPPResonanceSelector.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 
27 
28 using std::vector;
29 using std::ostringstream;
30 
31 using namespace genie;
32 
33 //___________________________________________________________________________
35 HadronicSystemGenerator("genie::RSPPResonanceSelector")
36 {
37 
38 }
39 //___________________________________________________________________________
41 HadronicSystemGenerator("genie::RSPPResonanceSelector", config)
42 {
43 
44 }
45 //___________________________________________________________________________
47 {
48 
49 }
50 //___________________________________________________________________________
52 {
53  //-- select a baryon resonance
54  Resonance_t res = this->SelectResonance(evrec);
55  assert(res != kNoResonance);
56 
57  //-- add the resonance at the event summary
58  Interaction * interaction = evrec->Summary();
59  interaction->ExclTagPtr()->SetResonance(res);
60 
61  //-- add an entry at the GHep event record & the event summary
62  this->AddResonance(evrec);
63 }
64 //___________________________________________________________________________
66 {
67 // Select a baryon resonance from a list of considered resonances based on
68 // their differential d^2xsec/dWdQ^2 cross sections
69 
70  LOG("RESSelector", pNOTICE) << "Selecting a baryon resonance";
71 
72  Interaction * interaction = evrec->Summary();
73 
74  //-- Figure out what the resonance charge should be.
75  int q_res = this->ResonanceCharge(evrec);
76 
77  //-- Use selected kinematics
78  interaction->KinePtr()->UseSelectedKinematics();
79 
80  //-- Trust kinematics and process type already set.
81  interaction->SetBit(kISkipProcessChk);
82  interaction->SetBit(kISkipKinematicChk);
83 
84  //-- Access cross section algorithm for running thread
85  // The algorithm must be able to compute the RES contribution to
86  // the specified RES/SPP channel
88  const EventGeneratorI * evg = rtinfo->RunningThread();
89  const XSecAlgorithmI * xsecalg = evg->CrossSectionAlg();
90 
91  //-- Loop over all considered baryon resonances and compute the double
92  // differential cross section for the selected kinematical variables
93 
94  double xsec_sum = 0;
95  unsigned int nres = fResList.NResonances();
96  vector<double> xsec_vec(nres);
97 
98  for(unsigned int ires = 0; ires < nres; ires++) {
99 
100  //-- Current resonance
101  Resonance_t res = fResList.ResonanceId(ires);
102 
103  //-- Set the current resonance at the interaction summary
104  // compute the differential cross section d^2xsec/dWdQ^2
105  // (do it only for resonances that can conserve charge)
106  interaction->ExclTagPtr()->SetResonance(res);
107 
108  double xsec = 0;
109  bool skip = (q_res==2 && !utils::res::IsDelta(res));
110 
111  if(!skip) xsec = xsecalg->XSec(interaction,kPSWQ2fE);
112  else {
113  SLOG("RESSelector", pNOTICE)
114  << "RES: " << utils::res::AsString(res)
115  << " would not conserve charge -- skipping it";
116  }
117  //-- For the ith resonance store the sum of (xsec) * (breit-wigner)
118  // for the resonances in the range [0,i]
119  xsec_sum += xsec;
120  xsec_vec[ires] = xsec_sum;
121 
122  SLOG("RESSelector", pNOTICE)
123  << "Resonances (0->" << ires << "): "
124  << "Sum{ BW(W) * d^2xsec(E,W,Q^2)/dWd*Q^2 } = " << xsec_sum;
125  } // res
126 
127  //-- Reset 'trust' bits
128  interaction->ResetBit(kISkipProcessChk);
129  interaction->ResetBit(kISkipKinematicChk);
130 
131  //-- Reset running kinematics
132  interaction->KinePtr()->ClearRunningValues();
133 
134  //-- Use the computed differential cross sections to select a resonance
135  RandomGen * rnd = RandomGen::Instance();
136  double R = xsec_sum * rnd->RndGen().Rndm();
137 
138  SLOG("RESSelector", pDEBUG) << "R = " << R;
139 
140  for(unsigned int ires = 0; ires < nres; ires++) {
141  SLOG("RESSelector", pDEBUG)
142  << "SUM-XSEC(0->" << ires <<") = " << xsec_vec[ires];
143 
144  if(R < xsec_vec[ires]) {
145  Resonance_t sres = fResList.ResonanceId(ires); // selected RES.
146  LOG("RESSelector", pNOTICE)
147  << "Selected RES = " << utils::res::AsString(sres);
148  return sres;
149  }
150  }
151  LOG("RESSelector", pERROR) << "** Failed to select a resonance";
152  return kNoResonance;
153 }
154 //___________________________________________________________________________
156 {
157  // compute RES p4 = p4(neutrino) + p4(hit nucleon) - p4(primary lepton)
158  TLorentzVector p4 = this->Hadronic4pLAB(evrec);
159 
160  //-- Determine the RES pdg code (from the selected Resonance_t & charge)
161  Interaction * interaction = evrec->Summary();
162  Resonance_t res = interaction->ExclTag().Resonance();
163  int charge = this->ResonanceCharge(evrec);
164  int pdgc = utils::res::PdgCode(res,charge);
165 
166  LOG("RESSelector", pNOTICE)
167  << "Adding RES with PDGC = " << pdgc << ", Q = " << charge;
168 
169  //-- Add the resonance at the EventRecord
171  int mom = evrec->HitNucleonPosition();
172 
173  evrec->AddParticle(
174  pdgc, ist, mom,-1,-1,-1, p4.Px(),p4.Py(),p4.Pz(),p4.E(), 0,0,0,0);
175 }
176 //___________________________________________________________________________
178 {
179  Algorithm::Configure(config);
180  this->LoadConfig();
181 }
182 //____________________________________________________________________________
183 void RSPPResonanceSelector::Configure(string param_set)
184 {
185  Algorithm::Configure(param_set);
186  this->LoadConfig();
187 }
188 //____________________________________________________________________________
190 {
191  fResList.Clear();
192  string resonances = "";
193  this->GetParam("ResonanceNameList", resonances);
194  SLOG("RESSelector", pDEBUG) << "Resonance list: " << resonances;
195 
196  fResList.DecodeFromNameList(resonances);
197  LOG("RESSelector", pINFO) << fResList;
198 }
199 //____________________________________________________________________________
bool IsDelta(Resonance_t res)
is it a Delta resonance?
Cross Section Calculation Interface.
#define pERROR
Definition: Messenger.h:59
Resonance_t SelectResonance(GHepRecord *event_rec) const
TLorentzVector Hadronic4pLAB(GHepRecord *event_rec) const
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
void DecodeFromNameList(string list, string delimiter=",")
Defines the EventGeneratorI interface.
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:418
enum genie::EResonance Resonance_t
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
unsigned int NResonances(void) const
void SetResonance(Resonance_t res)
Definition: XclsTag.cxx:128
Summary information for an interaction.
Definition: Interaction.h:56
void UseSelectedKinematics(void)
Definition: Kinematics.cxx:359
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -&gt; PDG code
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void ProcessEventRecord(GHepRecord *event_rec) const
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
#define pINFO
Definition: Messenger.h:62
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
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:65
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:80
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
const char * AsString(Resonance_t res)
resonance id -&gt; string
Abstract class. Is used to pass some commonly recurring methods to all concrete implementations of th...
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
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:84
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
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
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
Resonance_t ResonanceId(unsigned int ires) const