50 using namespace genie;
92 }
else if (
fXSecModel->
Id().
Name() ==
"genie::NievesSimoVacasMECPXSec2016") {
103 "ProcessEventRecord >> Cannot calculate kinematics for " <<
127 const TLorentzVector & p4cluster = *(cluster->
P4());
128 const TLorentzVector & p4tgt = *(target ->
P4());
130 const TLorentzVector
p4 = p4tgt - p4cluster;
131 const TLorentzVector v4(0.,0.,0., 0.);
133 int momidx =
event->TargetNucleusPosition();
169 GHepParticle * remnant_nucleus =
event->RemnantNucleus();
177 tgt.SetHitNucPdg(pdgv[0]);
180 tgt.SetHitNucPdg(pdgv[1]);
185 <<
"1st nucleon (code = " << pdgv[0] <<
") generated momentum: (" 186 << p3a.Px() <<
", " << p3a.Py() <<
", " << p3a.Pz() <<
"), " 187 <<
"|p| = " << p3a.Mag();
189 <<
"2nd nucleon (code = " << pdgv[1] <<
") generated momentum: (" 190 << p3b.Px() <<
", " << p3b.Py() <<
", " << p3b.Pz() <<
"), " 191 <<
"|p| = " << p3b.Mag();
195 TVector3
p3 = p3a + p3b;
198 <<
"di-nucleon cluster momentum: (" 199 << p3.Px() <<
", " << p3.Py() <<
", " << p3.Pz() <<
"), " 200 <<
"|p| = " << p3.Mag();
213 TLorentzVector p4nclust ( p3.Px(), p3.Py(), p3.Pz(), EN );
214 TLorentzVector p4remnant (-1*p3.Px(), -1*p3.Py(), -1*p3.Pz(), Mi-EN);
221 event->Summary()->InitStatePtr()->TgtPtr()->SetHitNucP4(p4nclust);
249 double dQ2 = (Q2max-Q2min) / (nq-1);
250 double dW = (Wmax-
Wmin ) / (nw-1);
252 for(
int iw=0; iw<nw; iw++) {
253 for(
int iq=0; iq<nq; iq++) {
254 double Q2 = Q2min + iq*dQ2;
255 double W = Wmin + iw*dW;
259 xsec_max = TMath::Max(xsec, xsec_max);
262 LOG(
"MEC",
pNOTICE) <<
"xsec_max (E = " << Ev <<
" GeV) = " << xsec_max;
266 unsigned int iter = 0;
272 <<
"Couldn't select a valid W, Q^2 pair after " 273 << iter <<
" iterations";
274 event->EventFlags()->SetBitNumber(
kKineGenErr,
true);
276 exception.
SetReason(
"Couldn't select kinematics");
282 double gQ2 = Q2min + (Q2max-Q2min) * rnd->
RndKine().Rndm();
283 double gW = Wmin + (Wmax -
Wmin ) * rnd->
RndKine().Rndm();
291 double t = xsec_max * rnd->
RndKine().Rndm();
293 accept = (t < J*
xsec);
297 LOG(
"MEC",
pINFO) <<
"Selected: Q^2 = " << gQ2 <<
", W = " << gW;
306 LOG(
"MEC",
pINFO) <<
"x = " << gx <<
", y = " << gy;
328 const TLorentzVector & pnuc4 = init_state.
Tgt().
HitNucP4();
329 TVector3
beta = pnuc4.BoostVector();
333 TLorentzVector * p4v =
event->Probe()->GetP4();
334 p4v->Boost(-1.*beta);
337 double Q2 = interaction->
Kine().
Q2(
true);
338 double y = interaction->
Kine().
y(
true);
342 LOG(
"MEC",
pNOTICE) <<
"neutrino energy = " << Ev;
344 double ml2 = TMath::Power(ml,2);
348 double El = (1-
y)*Ev;
349 double plp = El - 0.5*(Q2+ml2)/Ev;
350 double plt =
TMath::Sqrt(TMath::Max(0.,El*El-plp*plp-ml2));
353 <<
"fsl: E = " << El <<
", |p//| = " << plp <<
", |pT| = " << plt;
357 double phi = 2*
kPi * rnd->
RndLep().Rndm();
358 double pltx = plt * TMath::Cos(phi);
359 double plty = plt * TMath::Sin(phi);
363 TVector3 unit_nudir = p4v->Vect().Unit();
366 TVector3 p3l(pltx,plty,plp);
367 p3l.RotateUz(unit_nudir);
370 TLorentzVector p4l(p3l,El);
379 TLorentzVector v4(*event->
Probe()->
X4());
381 int momidx =
event->ProbePosition();
391 TLorentzVector *
tmp=nucleon_cluster->
GetP4();
392 TLorentzVector p4cluster(*tmp);
398 TLorentzVector p4v(*neutrino->
P4());
403 TLorentzVector p4l(*fsl->P4());
406 TLorentzVector q = p4v - p4l;
409 TLorentzVector p4cluster_recoil = p4cluster + q;
412 LOG(
"MEC",
pINFO) <<
"Interaction summary";
413 LOG(
"MEC",
pINFO) << *
event->Summary();
414 int recoil_nucleon_cluster_pdg =
event->Summary()->RecoilNucleonPdg();
417 TLorentzVector v4(*neutrino->
X4());
422 2, -1, -1, -1, p4cluster_recoil, v4);
430 LOG(
"MEC",
pINFO) <<
"Decaying nucleon cluster...";
433 int nucleon_cluster_id = 5;
434 GHepParticle * nucleon_cluster =
event->Particle(nucleon_cluster_id);
439 LOG(
"MEC",
pINFO) <<
"Decay product IDs: " << pdgv;
442 vector<int>::const_iterator pdg_iter;
444 double * mass =
new double[pdgv.size()];
446 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
447 int pdgc = *pdg_iter;
454 <<
"Performing a phase space decay to " 455 << pdgv.size() <<
" particles / total mass = " <<
sum;
457 TLorentzVector * p4d = nucleon_cluster->
GetP4();
458 TLorentzVector * v4d = nucleon_cluster->
GetX4();
467 <<
" *** Phase space decay is not permitted \n" 468 <<
" Total particle mass = " << sum <<
"\n" 477 exception.
SetReason(
"Decay not permitted kinematically");
485 for(
int idec=0; idec<200; idec++) {
487 wmax = TMath::Max(wmax,w);
493 <<
"Max phase space gen. weight = " << wmax;
496 bool accept_decay=
false;
505 <<
"Couldn't generate an unweighted phase space decay after " 506 << itry <<
" attempts";
514 exception.
SetReason(
"Couldn't select decay after N attempts");
522 <<
"Decay weight = " << w <<
" > max decay weight = " << wmax;
524 double gw = wmax * rnd->
RndDec().Rndm();
525 accept_decay = (gw<=
w);
528 <<
"Decay weight = " << w <<
" / R = " << gw
529 <<
" - accepted: " << accept_decay;
534 TLorentzVector v4(*v4d);
537 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
538 int pdgc = *pdg_iter;
540 event->AddParticle(pdgc, ist, nucleon_cluster_id,-1,-1,-1, *p4fin, v4);
552 bool allowdup =
true;
572 <<
"Unknown di-nucleon cluster PDG code (" << pdgc <<
")";
584 int FullDeltaNodelta = 1;
590 double XSecMaxPar1 = 2.2504;
591 double XSecMaxPar2 = 9.41158;
602 TLorentzVector v4(*event->
Probe()->
X4());
603 TLorentzVector tempp4(0.,0.,0.,0.);
607 double CosthMax = 1.0;
608 double CosthMin = -1.0;
625 TMax = Enu - LepMass;
635 TMin =
TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu -
fQ3Max), 2)) - LepMass;
641 double NuclearAfactorXSecMax = 1.0;
647 if (NuclearA < 12) NuclearAfactorXSecMax *= NuclearA / 12.0;
648 else NuclearAfactorXSecMax *= TMath::Power(NuclearA/12.0, 1.4);
651 LOG(
"MEC",
pERROR) <<
"Trying to scale XSecMax for larger nuclei, but " 652 << TgtPDG <<
" isn't a nucleus?";
661 unsigned int iter = 0;
669 <<
"Couldn't select a valid Tmu, CosTheta pair after " 670 << iter <<
" iterations";
671 event->EventFlags()->SetBitNumber(
kKineGenErr,
true);
673 exception.
SetReason(
"Couldn't select lepton kinematics");
679 T = TMin + (TMax-TMin)*rnd->
RndKine().Rndm();
680 Costh = CosthMin + (CosthMax-CosthMin)*rnd->
RndKine().Rndm();
684 Q3 =
TMath::Sqrt(Plep*Plep + Enu*Enu - 2.0 * Plep * Enu * Costh);
703 if (FullDeltaNodelta == 1){
709 double XSecMax = 1.35 * TMath::Power(10.0, XSecMaxPar1 * TMath::Log10(Enu) - XSecMaxPar2);
710 if (NuclearA > 12) XSecMax *= NuclearAfactorXSecMax;
712 LOG(
"MEC",
pDEBUG) <<
" T, Costh: " << T <<
", " << Costh ;
737 if (XSec > XSecMax) {
738 LOG(
"MEC",
pERROR) <<
"XSec is > XSecMax for nucleus " << TgtPDG <<
" " 739 << XSec <<
" > " << XSecMax
740 <<
" don't let this happen.";
743 accept = XSec > XSecMax*rnd->
RndKine().Rndm();
744 LOG(
"MEC",
pINFO) <<
"Xsec, Max, Accept: " << XSec <<
", " 745 << XSecMax <<
", " << accept;
758 double myrand = rnd->
RndKine().Rndm();
759 double pnFraction = XSecPN / XSec;
760 LOG(
"MEC",
pDEBUG) <<
"Test for pn: xsec_pn = " << XSecPN
761 <<
"; xsec = " << XSec
762 <<
"; pn_fraction = " << pnFraction
763 <<
"; random number val = " << myrand;
765 if (myrand <= pnFraction) {
768 1, -1, -1, -1, tempp4, v4);
772 if (rnd->
RndKine().Rndm() <= XSecDeltaPN / XSecPN) {
780 1, -1, -1, -1, tempp4, v4);
785 1, -1, -1, -1, tempp4, v4);
791 (XSecDelta - XSecDeltaPN) / (XSec - XSecPN)) {
824 double PlepZ = Plep * Costh;
825 double PlepXY = Plep *
TMath::Sqrt(1. - TMath::Power(Costh,2));
828 double phi= 2 *
kPi * rnd->
RndLep().Rndm();
830 double PlepX = PlepXY * TMath::Cos(phi);
831 double PlepY = PlepXY * TMath::Sin(phi);
835 TVector3 unit_nudir =
event->Probe()->P4()->Vect().Unit();
836 TVector3 p3l(PlepX, PlepY, PlepZ);
837 p3l.RotateUz(unit_nudir);
840 Elep =
TMath::Sqrt(LepMass*LepMass + PlepX*PlepX + PlepY*PlepY + PlepZ*PlepZ);
841 TLorentzVector p4l(p3l,Elep);
845 int momidx =
event->ProbePosition();
851 double gy = Q0 / Enu;
866 LOG(
"MEC",
pDEBUG) <<
"~~~ LEPTON DONE ~~~";
874 LOG(
"MEC",
pDEBUG) <<
"Generate Initial Hadrons - Start";
879 TLorentzVector p4nu(*neutrino->
P4());
884 TLorentzVector p4l(*fsl->P4());
887 TLorentzVector Q4 = p4nu - p4l;
893 GHepParticle * initial_nucleon_cluster =
event->HitNucleon();
894 assert(initial_nucleon_cluster);
895 GHepParticle * remnant_nucleus =
event->RemnantNucleus();
915 unsigned int iter = 0;
917 int initial_nucleon_cluster_pdg = initial_nucleon_cluster->
Pdg();
918 int final_nucleon_cluster_pdg = 0;
919 if (neutrino->
Pdg() > 0) {
927 LOG(
"MEC",
pERROR) <<
"Wrong pdg for a CC neutrino MEC interaction" 928 << initial_nucleon_cluster->
Pdg();
931 else if (neutrino->
Pdg() < 0) {
939 LOG(
"MEC",
pERROR) <<
"Wrong pdg for a CC anti-neutrino MEC interaction" 940 << initial_nucleon_cluster->
Pdg();
944 TLorentzVector p4initial_cluster;
945 TLorentzVector p4final_cluster;
946 TLorentzVector p4remnant_nucleus;
947 double removalenergy1;
948 double removalenergy2;
960 <<
"Couldn't select a valid W, Q^2 pair after " 961 << iter <<
" iterations";
962 event->EventFlags()->SetBitNumber(
kKineGenErr,
true);
964 exception.
SetReason(
"Couldn't select initial hadron kinematics");
975 tgt.SetHitNucPdg(pdgv[0]);
979 tgt.SetHitNucPdg(pdgv[1]);
987 TVector3 p3i = p31i + p32i;
991 p4initial_cluster.SetPxPyPzE(p3i.Px(),p3i.Py(),p3i.Pz(),
energy);
994 TLorentzVector tLVebind(0., 0., 0., -1.0 * (removalenergy1 + removalenergy2));
1009 p4final_cluster = p4initial_cluster + Q4 + tLVebind;
1013 if (p4final_cluster.M() <
1026 initial_nucleon_cluster->
SetMomentum(p4initial_cluster);
1030 remnant_nucleus->
SetMomentum(-1.0*p4initial_cluster.Px(),
1031 -1.0*p4initial_cluster.Py(),
1032 -1.0*p4initial_cluster.Pz(),
1033 Mi - p4initial_cluster.E() + removalenergy1 + removalenergy2);
1040 TLorentzVector v4(*neutrino->
X4());
1044 2, -1, -1, -1, p4final_cluster, v4);
1053 event->AddParticle(p1);
1073 RgKey nuclkey =
"NuclearModel";
void ProcessEventRecord(GHepRecord *event) const
TLorentzVector * GetX4(void) const
void SelectNSVLeptonKinematics(GHepRecord *event) const
double RemovalEnergy(void) const
const NuclearModelI * fNuclModel
double J(double q0, double q3, double Enu, double ml)
THE MAIN GENIE PROJECT NAMESPACE
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
static const double kNucleonMass
double Q2(const Interaction *const i)
static RandomGen * Instance()
Access instance.
int HitNucPdg(void) const
void SetQ2(double Q2, bool selected=false)
const TLorentzVector * P4(void) const
Kinematics * KinePtr(void) const
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the 'Visito...
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 Configure(const Registry &config)
Defines the EventGeneratorI interface.
int IonPdgCodeToA(int pdgc)
Generated/set kinematical variables for an event.
string P4AsString(const TLorentzVector *p)
static const unsigned int kMaxUnweightDecayIterations
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
::xsd::cxx::tree::exception< char > exception
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
A singleton holding random number generator classes. All random number generation in GENIE should tak...
void SetResonance(Resonance_t res)
double XYtoW(double Ev, double M, double x, double y)
const TVector3 & Momentum3(void) const
virtual GHepParticle * Probe(void) const
double y(bool selected=false) const
void GenerateNSVInitialHadrons(GHepRecord *event) const
Summary information for an interaction.
const TLorentzVector & HitNucP4(void) const
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...
void SetFSLeptonP4(const TLorentzVector &p4)
void RecoilNucleonCluster(GHepRecord *event) const
void DecayNucleonCluster(GHepRecord *event) const
TLorentzVector * GetP4(void) const
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
const Kinematics & Kine(void) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
virtual void Configure(const Registry &config)
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Misc GENIE control constants.
PDGCodeList NucleonClusterConstituents(int pdgc) const
void SelectEmpiricalKinematics(GHepRecord *event) const
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Setx(double x, bool selected=false)
void SetReturnStep(int s)
XclsTag * ExclTagPtr(void) const
void SetKV(KineVar_t kv, double value)
static RunningThreadInfo * Instance(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
void SetW(double W, bool selected=false)
void SwitchOnFastForward(void)
static PDGLibrary * Instance(void)
void SetReason(string reason)
double Q2YtoX(double Ev, double M, double Q2, double y)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
const TLorentzVector * X4(void) const
A registry. Provides the container for algorithm configuration parameters.
void SetHitNucPdg(int pdgc)
void Sety(double y, bool selected=false)
Target * TgtPtr(void) const
TParticlePDG * Find(int pdgc)
int IonPdgCode(int A, int Z)
const XSecAlgorithmI * fXSecModel
void SwitchOnStepBack(void)
assert(nhit_max >=nhit_nbins)
void SetHadSystP4(const TLorentzVector &p4)
InitialState * InitStatePtr(void) const
virtual bool GenerateNucleon(const Target &) const =0
const InitialState & InitState(void) const
void AddTargetRemnant(GHepRecord *event) const
const ProcessInfo & ProcInfo(void) const
void GenerateFermiMomentum(GHepRecord *event) const
double Q2(bool selected=false) const
void ClearRunningValues(void)
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
static const unsigned int kRjMaxIterations
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const EventGeneratorI * RunningThread(void)
double ProbeE(RefFrame_t rf) const
GENIE's GHEP MC event record.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TRandom3 & RndDec(void) const
rnd number generator used by decay models
STDHEP-like event record entry that can fit a particle or a nucleus.
Keep info on the event generation thread currently on charge. This is used so that event generation m...
void AddFinalStateLepton(GHepRecord *event) const
Root of GENIE utility namespaces.
void push_back(int pdg_code)
TGenPhaseSpace fPhaseSpaceGenerator
Initial State information.
const Algorithm * SubAlg(const RgKey ®istry_key) const