38 using namespace genie;
95 TLorentzVector v4(0,0,0,0);
105 TLorentzVector p4i(0,0,0,Mi);
106 event->AddParticle(ipdg,stis,-1,-1,-1,-1, p4i, v4);
111 TLorentzVector p4neut(0,0,0,mneut);
112 event->AddParticle(neutpdg,stdc,0,-1,-1,-1, p4neut, v4);
117 TLorentzVector p4n(0,0,0,mn);
118 event->AddParticle(dpdg,stdc, 0,-1,-1,-1, p4n, v4);
127 TLorentzVector p4f(0,0,0,Mf);
128 event->AddParticle(rpdg,stfs,0,-1,-1,-1, p4f, v4);
135 int A = initial_nucleus->
A();
143 double dA = (double)A;
144 double R = R0 * TMath::Power(dA, 1./3.);
147 <<
"Generating vertex according to a realistic nuclear density profile";
153 for(
double r = 0;
r < rmax;
r+=dr) {
159 TLorentzVector vtx(0,0,0,0);
160 unsigned int iter = 0;
167 <<
"Couldn't generate a vertex position after " << iter <<
" iterations";
169 exception.
SetReason(
"Couldn't generate vertex");
174 double r = rmax * rnd->
RndFsi().Rndm();
175 double t = ymax * rnd->
RndFsi().Rndm();
179 <<
"y = " << y <<
" > ymax = " << ymax <<
" for r = " << r <<
", A = " <<
A;
181 bool accept = (t <
y);
184 double cosphi = TMath::Cos(phi);
185 double sinphi = TMath::Sin(phi);
186 double costheta = -1 + 2 * rnd->
RndFsi().Rndm();
187 double sintheta =
TMath::Sqrt(1-costheta*costheta);
188 vtx.SetX(r*sintheta*cosphi);
189 vtx.SetY(r*sintheta*sinphi);
190 vtx.SetZ(r*costheta);
197 GHepParticle * oscillating_neutron =
event->Particle(1);
198 assert(oscillating_neutron);
202 GHepParticle * annihilation_nucleon =
event->Particle(2);
203 assert(annihilation_nucleon);
211 int A = initial_nucleus->
A();
216 GHepParticle * oscillating_neutron =
event->Particle(1);
217 GHepParticle * annihilation_nucleon =
event->Particle(2);
219 assert(oscillating_neutron);
220 assert(annihilation_nucleon);
234 double mass = oscillating_neutron->
Mass();
237 TLorentzVector
p4(p3, energy);
241 <<
"Generated neutron momentum: (" 242 << p3.Px() <<
", " << p3.Py() <<
", " << p3.Pz() <<
"), " 243 <<
"|p| = " << p3.Mag();
246 tgt.SetHitNucPdg(annihilation_nucleon->
Pdg());
252 mass = annihilation_nucleon->
Mass();
253 energy =
sqrt(
pow(mass,2) + p3.Mag2()) - w;
255 p4 = TLorentzVector(p3, energy);
259 <<
"Generated nucleon momentum: (" 260 << p3.Px() <<
", " << p3.Py() <<
", " << p3.Pz() <<
"), " 261 <<
"|p| = " << p3.Mag();
264 p3 = -1 * (oscillating_neutron->
GetP4()->Vect() + annihilation_nucleon->
GetP4()->Vect());
266 mass = remnant_nucleus->
Mass();
267 energy =
sqrt(
pow(mass,2) + p3.Mag2());
269 p4 = TLorentzVector(p3, energy);
276 LOG(
"NNBarOsc",
pINFO) <<
"Generating decay...";
279 LOG(
"NNBarOsc",
pINFO) <<
"Decay product IDs: " << pdgv;
280 assert ( pdgv.size() > 1);
282 LOG(
"NNBarOsc",
pINFO) <<
"Performing a phase space decay...";
286 vector<int>::const_iterator pdg_iter;
288 double * mass =
new double[pdgv.size()];
290 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
291 int pdgc = *pdg_iter;
298 <<
"Decaying N = " << pdgv.size() <<
" particles / total mass = " <<
sum;
299 int initial_nucleus_id = 0;
300 int oscillating_neutron_id = 1;
301 int annihilation_nucleon_id = 2;
304 GHepParticle * initial_nucleus =
event->Particle(initial_nucleus_id);
306 GHepParticle * oscillating_neutron =
event->Particle(oscillating_neutron_id);
307 assert(oscillating_neutron);
308 GHepParticle * annihilation_nucleon =
event->Particle(annihilation_nucleon_id);
309 assert(annihilation_nucleon);
315 TLorentzVector * p4_1 = oscillating_neutron->
GetP4();
316 TLorentzVector * p4_2 = annihilation_nucleon->
GetP4();
317 TLorentzVector * p4d =
new TLorentzVector(*p4_1 + *p4_2);
318 TVector3
boost = p4d->BoostVector();
322 TLorentzVector * v4d = annihilation_nucleon->
GetX4();
337 <<
"Not enough energy to generate decay products! Selecting a new final state.";
341 int initial_nucleus_pdg = initial_nucleus->
Pdg();
344 if (nucleus_pdg.size() != 10) {
346 <<
"Expecting the nuclear PDG code to be a 10-digit integer, but it is " << nucleus_pdg <<
", which is " 347 << nucleus_pdg.size() <<
" digits long. Drop me an email at jezhewes@gmail.com ; exiting...";
351 int n_nucleons = std::stoi(nucleus_pdg.substr(6,3)) - 1;
352 int n_protons = std::stoi(nucleus_pdg.substr(3,3));
353 double proton_frac = ((double)n_protons) / ((double) n_nucleons);
354 double neutron_frac = 1 - proton_frac;
357 const int n_modes = 16;
358 double br [n_modes] = { 0.010, 0.080, 0.100, 0.220,
359 0.360, 0.160, 0.070, 0.020,
360 0.015, 0.065, 0.110, 0.280,
361 0.070, 0.240, 0.100, 0.100 };
363 for (
int i = 0;
i < n_modes;
i++) {
366 br[
i] *= proton_frac;
369 br[
i] *= neutron_frac;
375 double p = rnd->
RndNum().Rndm();
379 for (
int j = 0;
j < n_modes;
j++) {
393 event->AttachSummary(interaction);
398 LOG(
"NNBarOsc",
pINFO) <<
"Decay product IDs: " << pdgv;
399 assert ( pdgv.size() > 1);
402 LOG(
"NNBarOsc",
pINFO) <<
"Performing a phase space decay...";
405 mass =
new double[pdgv.size()];
407 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
408 int pdgc = *pdg_iter;
415 <<
"Decaying N = " << pdgv.size() <<
" particles / total mass = " <<
sum;
424 <<
" *** Phase space decay is not permitted \n" 425 <<
" Total particle mass = " << sum <<
"\n" 433 exception.
SetReason(
"Decay not permitted kinematically");
441 for(
int i=0;
i<200;
i++) {
443 wmax = TMath::Max(wmax,w);
449 <<
"Max phase space gen. weight @ current hadronic system: " << wmax;
454 bool accept_decay=
false;
463 <<
"Couldn't generate an unweighted phase space decay after " 464 << itry <<
" attempts";
471 exception.
SetReason(
"Couldn't select decay after N attempts");
478 <<
"Decay weight = " << w <<
" > max decay weight = " << wmax;
480 double gw = wmax * rnd->
RndHadro().Rndm();
481 accept_decay = (gw<=
w);
484 <<
"Decay weight = " << w <<
" / R = " << gw
485 <<
" - accepted: " << accept_decay;
490 TLorentzVector v4(*v4d);
492 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
493 int pdgc = *pdg_iter;
498 event->AddParticle(pdgc, ist, oscillating_neutron_id,-1,-1,-1, *p4fin, v4);
527 RgKey nuclkey =
"NuclearModel";
TLorentzVector * GetX4(void) const
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
double RemovalEnergy(void) const
int AnnihilatingNucleonPdgCode(NNBarOscMode_t ndm)
THE MAIN GENIE PROJECT NAMESPACE
NNBarOscMode_t fCurrDecayMode
static RandomGen * Instance()
Access instance.
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the 'Visito...
double Density(double r, int A, double ring=0.)
GHepStatus_t DecayProductStatus(bool in_nucleus, int pdgc)
NNBarOscPrimaryVtxGenerator()
enum genie::EGHepStatus GHepStatus_t
int IonPdgCodeToA(int pdgc)
string P4AsString(const TLorentzVector *p)
static const unsigned int kMaxUnweightDecayIterations
~NNBarOscPrimaryVtxGenerator()
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
::xsd::cxx::tree::exception< char > exception
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
const TVector3 & Momentum3(void) const
Summary information for an interaction.
const NuclearModelI * fNuclModel
void SetPosition(const TLorentzVector &v4)
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
enum genie::ENNBarOscMode NNBarOscMode_t
TLorentzVector * GetP4(void) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
virtual void Configure(const Registry &config)
void AddInitialState(GHepRecord *event) const
int DecayMode(void) const
string AsString(NNBarOscMode_t ndm)
void Configure(const Registry &config)
TRandom3 & RndNum(void) const
rnd number generator used by MC integrators & other numerical methods
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
PDGCodeList DecayProductList(NNBarOscMode_t ndm)
void SwitchOnFastForward(void)
static PDGLibrary * Instance(void)
void SetReason(string reason)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
A registry. Provides the container for algorithm configuration parameters.
void SetHitNucPdg(int pdgc)
const XclsTag & ExclTag(void) const
TParticlePDG * Find(int pdgc)
void ProcessEventRecord(GHepRecord *event) const
int IonPdgCode(int A, int Z)
assert(nhit_max >=nhit_nbins)
virtual bool GenerateNucleon(const Target &) const =0
const InitialState & InitState(void) const
int IonPdgCodeToZ(int pdgc)
void GenerateOscillatingNeutronPosition(GHepRecord *event) const
const Target & Tgt(void) const
std::string to_string(ModuleType mt)
static const unsigned int kRjMaxIterations
GENIE's GHEP MC event record.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
TGenPhaseSpace fPhaseSpaceGenerator
void GenerateFermiMomentum(GHepRecord *event) const
static Interaction * NOsc(int tgt, int annihilation_mode=-1)
void GenerateDecayProducts(GHepRecord *event) const
void SetSeed(long int seed)
const Algorithm * SubAlg(const RgKey ®istry_key) const