36 #include "Math/Minimizer.h" 37 #include "Math/Factory.h" 56 using namespace genie;
83 <<
"Generating kinematics uniformly over the allowed phase space";
94 }
else if (
fXSecModel->
Id().
Name() ==
"genie::BergerSehgalFMCOHPiPXSec2015") {
101 "ProcessEventRecord >> Cannot calculate kinematics for " <<
122 double xsec_max = this->
MaxXSec(evrec);
132 const double dy = ymax -
ymin;
135 const double dQ2 = Q2max - Q2min;
137 unsigned int iter = 0;
149 gy = ymin + dy * rnd->
RndKine().Rndm();
153 "Trying: Q^2 = " <<
gQ2 <<
", y = " << gy;
163 accept = (xsec_max * rnd->
RndKine().Rndm() <
xsec);
168 <<
"Selected: Q^2 = " <<
gQ2 <<
", y = " << gy;
179 double Epi2 = TMath::Power(Epi,2);
181 double gx = interaction->
KinePtr()->
x();
183 double tA = 1. + xME - 0.5*pme2;
186 double tmin = 2*Epi2 * (tA-tB);
187 double tmax = 2*Epi2 * (tA+tB);
188 double A = (double) init_state.
Tgt().
A();
189 double A13 = TMath::Power(A,1./3.);
191 double R2 = TMath::Power(R,2.);
192 double b = 0.33333 * R2;
193 double tsum = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
194 double rt = tsum * rnd->
RndKine().Rndm();
195 double gt = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/
b;
236 double xsec_max = this->
MaxXSec(evrec);
246 const double dy = ymax -
ymin;
249 const double dQ2 = Q2max - Q2min;
252 const double dt = tmax - tmin;
256 unsigned int iter = 0;
268 gy = ymin + dy * rnd->
RndKine().Rndm();
273 "Trying: Q^2 = " <<
gQ2 <<
", y = " << gy <<
", t = " <<
gt;
283 accept = (xsec_max * rnd->
RndKine().Rndm() <
xsec);
288 <<
"Selected: Q^2 = " <<
gQ2 <<
", y = " << gy <<
", t = " <<
gt;
339 const double dx = xmax -
xmin;
340 const double dy = ymax -
ymin;
344 unsigned int iter = 0;
346 double xsec=-1, gx=-1, gy=-1;
354 gx = xmin + dx * rnd->
RndKine().Rndm();
355 gy = ymin + dy * rnd->
RndKine().Rndm();
361 LOG(
"COHKinematics",
pNOTICE) <<
"Initializing the sampling envelope";
363 fEnvelope->SetRange(xmin,ymin,xmax,ymax);
372 LOG(
"COHKinematics",
pINFO) <<
"Trying: x = " << gx <<
", y = " << gy;
383 double t = max * rnd->
RndKine().Rndm();
386 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 388 <<
"xsec= " << xsec <<
", J= 1, Rnd= " <<
t;
398 LOG(
"COHKinematics",
pNOTICE) <<
"Selected: x = "<< gx <<
", y = "<< gy;
409 double Epi2 = TMath::Power(Epi,2);
412 double tA = 1. + xME - 0.5*pme2;
414 double tmin = 2*Epi2 * (tA-tB);
415 double tmax = 2*Epi2 * (tA+tB);
416 double A = (double) init_state.
Tgt().
A();
417 double A13 = TMath::Power(A,1./3.);
419 double R2 = TMath::Power(R,2.);
420 double b = 0.33333 * R2;
421 double tsum = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
422 double rt = tsum * rnd->
RndKine().Rndm();
423 double gt = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/
b;
426 <<
"Selected: t = "<< gt <<
", from ["<< tmin <<
", "<< tmax <<
"]";
432 double totxsec = evrec->
XSec();
433 double wght = (vol/totxsec)*xsec;
434 LOG(
"COHKinematics",
pNOTICE) <<
"Kinematics wght = "<< wght;
438 LOG(
"COHKinematics",
pNOTICE) <<
"Current event wght = " << wght;
465 LOG(
"COHKinematics",
pNOTICE) <<
"Using AlvarezRuso Model";
484 const double E_l_min = interaction->
FSPrimLepton()->Mass();
487 const double ctheta_l_min = 0.4;
490 const double ctheta_pi_min = 0.4;
491 const double ctheta_pi_max = 1.0 -
kASmallNum;
496 const double d_E_l = E_l_max - E_l_min;
497 const double d_ctheta_l = ctheta_l_max - ctheta_l_min;
498 const double d_ctheta_pi = ctheta_pi_max - ctheta_pi_min;
499 const double d_phi = phi_max - phi_min;
502 unsigned int iter = 0;
504 double xsec=-1, g_E_l=-1, g_theta_l=-1, g_phi_l=-1, g_theta_pi=-1, g_phi_pi=-1;
505 double g_ctheta_l, g_ctheta_pi;
512 g_E_l = E_l_min + d_E_l * rnd->
RndKine().Rndm();
513 g_ctheta_l = ctheta_l_min + d_ctheta_l * rnd->
RndKine().Rndm();
514 g_ctheta_pi = ctheta_pi_min + d_ctheta_pi * rnd->
RndKine().Rndm();
515 g_phi_l = phi_min + d_phi * rnd->
RndKine().Rndm();
517 g_phi_pi = g_phi_l + (phi_min + d_phi * rnd->
RndKine().Rndm());
518 g_theta_l = TMath::ACos(g_ctheta_l);
519 g_theta_pi = TMath::ACos(g_ctheta_pi);
521 LOG(
"COHKinematics",
pINFO) <<
"Trying: Lep(" <<g_E_l <<
", " <<
522 g_theta_l <<
", " << g_phi_l <<
") Pi(" <<
523 g_theta_pi <<
", " << g_phi_pi <<
")";
525 this->
SetKinematics(g_E_l, g_theta_l, g_phi_l, g_theta_pi, g_phi_pi,
526 interaction, interaction->
KinePtr());
533 double t = xsec_max * rnd->
RndKine().Rndm();
535 LOG(
"COHKinematics",
pINFO) <<
"Got: xsec = " << xsec <<
", t = " <<
536 t <<
" (max_xsec = " << xsec_max <<
")";
539 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 541 <<
"xsec= " << xsec <<
", J= 1, Rnd= " <<
t;
551 LOG(
"COHKinematics",
pNOTICE) <<
"Selected: Lepton(" <<
552 g_E_l <<
", " << g_theta_l <<
", " <<
553 g_phi_l <<
") Pion(" << g_theta_pi <<
", " << g_phi_pi <<
")";
556 double theta_l = g_theta_l;
557 double theta_pi = g_theta_pi;
558 double phi_l = g_phi_l;
559 double phi_pi = g_phi_pi;
561 double E_nu = P4_nu.E();
562 double E_pi= E_nu-E_l;
564 double m_pi = this->
pionMass(interaction);
567 TVector3 lepton_3vector = TVector3(0,0,0);
568 lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
569 TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
572 TVector3 pion_3vector = TVector3(0,0,0);
573 pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
574 TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
576 TLorentzVector q = P4_nu - P4_lep;
577 double Q2 = -q.Mag2();
579 double y = E_pi/E_nu;
581 double t = TMath::Abs( (q - P4_pion).
Mag2() );
587 double vol = d_E_l*d_ctheta_l*d_phi*d_ctheta_pi*d_phi;
588 double totxsec = evrec->
XSec();
589 double wght = (vol/totxsec)*xsec;
590 LOG(
"COHKinematics",
pNOTICE) <<
"Kinematics wght = "<< wght;
594 LOG(
"COHKinematics",
pNOTICE) <<
"Current event wght = " << wght;
616 const double theta_l,
618 const double theta_pi,
624 double E_nu = P4_nu.E();
625 double E_pi= E_nu-E_l;
637 TVector3 lepton_3vector = TVector3(0,0,0);
638 lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
639 TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
645 TVector3 pion_3vector = TVector3(0,0,0);
646 pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
647 TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
649 double Q2 = -(P4_nu-P4_lep).
Mag2();
651 double y = E_pi/E_nu;
669 double E_nu = P4_nu.E();
670 double E_pi= E_nu-E_l;
695 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 697 <<
"Scanning the allowed phase space {K} for the max(dxsec/d{K})";
699 double max_xsec = 0.;
702 }
else if ((
fXSecModel->
Id().
Name() ==
"genie::BergerSehgalCOHPiPXSec2015")) {
704 }
else if ((
fXSecModel->
Id().
Name() ==
"genie::BergerSehgalFMCOHPiPXSec2015")) {
711 "ComputeMaxXSec >> Cannot calculate max cross-section for " <<
719 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 721 SLOG(
"COHKinematics",
pDEBUG) <<
"Max xsec in phase space = " << max_xsec;
739 const double logQ2max = TMath::Log10(Q2r.
max);
740 const double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
742 for(
int i=0;
i<NQ2;
i++) {
743 double Q2 = TMath::Power(10, logQ2min +
i * dlogQ2);
748 (yr.
max > 1) || (yr.
min < 0)) {
751 const double logymin = TMath::Log10(yr.
min);
752 const double logymax = TMath::Log10(yr.
max);
753 const double dlogy = (logymax - logymin) /(Ny-1);
755 for(
int j=0;
j<Ny;
j++) {
756 double gy = TMath::Power(10, logymin +
j * dlogy);
764 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 766 <<
"xsec(Q2= " << Q2 <<
", y= " << gy <<
", t = " <<
gt <<
") = " <<
xsec;
768 max_xsec = TMath::Max(max_xsec, xsec);
788 const double logQ2max = TMath::Log10(Q2r.
max);
789 const double logtmin = TMath::Log10(
kASmallNum);
791 const double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
792 const double dlogt = (logtmax - logtmin) /(Nt-1);
794 for(
int i=0;
i<NQ2;
i++) {
795 double Q2 = TMath::Power(10, logQ2min +
i * dlogQ2);
800 (yr.
max > 1) || (yr.
min < 0)) {
803 const double logymin = TMath::Log10(yr.
min);
804 const double logymax = TMath::Log10(yr.
max);
805 const double dlogy = (logymax - logymin) /(Ny-1);
807 for(
int j=0;
j<Ny;
j++) {
808 double gy = TMath::Power(10, logymin +
j * dlogy);
810 for(
int k=0; k<Nt; k++) {
811 double gt = TMath::Power(10, logtmin + k * dlogt);
817 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 819 <<
"xsec(Q2= " << Q2 <<
", y= " << gy <<
", t = " << gt <<
") = " <<
xsec;
821 max_xsec = TMath::Max(max_xsec, xsec);
840 const double logxmin = TMath::Log10(1
E-5);
841 const double logxmax = TMath::Log10(1.0);
842 const double logymin = TMath::Log10(y.
min);
843 const double logymax = TMath::Log10(y.
max);
845 const double dlogx = (logxmax - logxmin) /(Nx-1);
846 const double dlogy = (logymax - logymin) /(Ny-1);
848 for(
int i=0;
i<Nx;
i++) {
849 double gx = TMath::Power(10, logxmin +
i * dlogx);
850 for(
int j=0;
j<Ny;
j++) {
851 double gy = TMath::Power(10, logymin +
j * dlogy);
854 if(Ev>1.0 && Q2>0.01)
continue;
860 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 862 <<
"xsec(x= " << gx <<
", y= " << gy <<
") = " <<
xsec;
864 max_xsec = TMath::Max(max_xsec, xsec);
878 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 880 <<
"Scanning the allowed phase space {K} for the max(dxsec/d{K})";
882 double max_xsec = 0.;
888 ROOT::Math::Minimizer *
min = ROOT::Math::Factory::CreateMinimizer(
"Minuit2");
892 min->SetFunction( f );
893 min->SetMaxFunctionCalls(10000);
894 min->SetTolerance(0.05);
898 const unsigned int n_el = 100;
899 const double d_el = (max_el - min_el) /
double(n_el - 1);
902 const double max_thetal =
kPi / 4.0;
903 const unsigned int n_thetal = 10;
904 const double d_thetal = (max_thetal - min_thetal) /
double(n_thetal - 1);
907 const double max_thetapi =
kPi / 2.0;
908 const unsigned int n_thetapi = 10;
909 const double d_thetapi = (max_thetapi - min_thetapi) /
double(n_thetapi - 1);
915 const unsigned int n_phipi = 10;
916 const double d_phipi = (max_phipi - min_phipi) /
double(n_phipi - 1);
918 min->SetLimitedVariable ( 0 ,
"E_lep" , max_el -kASmallNum , d_el , min_el , max_el );
919 min->SetLimitedVariable ( 1 ,
"theta_l" , min_thetal +kASmallNum , d_thetal , min_thetal , max_thetal );
920 min->SetLimitedVariable ( 2 ,
"theta_pi" , min_thetapi+kASmallNum , d_thetapi , min_thetapi, max_thetapi );
921 min->SetLimitedVariable ( 3 ,
"phi_pi" , min_phipi +kASmallNum , d_phipi , min_phipi , max_phipi );
924 max_xsec = -min->MinValue();
930 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 932 SLOG(
"COHKinematics",
pDEBUG) <<
"Max xsec in phase space = " << max_xsec;
966 <<
"*** Could not select valid kinematics after " 967 << iters <<
" iterations";
970 exception.
SetReason(
"Couldn't select kinematics");
1016 gROOT->GetListOfFunctions()->Remove(
fEnvelope);
double COHImportanceSamplingEnvelope(double *x, double *par)
virtual double MaxXSec(GHepRecord *evrec) const
const KPhaseSpace & PhaseSpace(void) const
virtual void SetWeight(double wght)
bool IsWeakCC(void) const
std::map< std::string, double > xmax
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double MaxXSec_ReinSehgal(const Interaction *in) const
THE MAIN GENIE PROJECT NAMESPACE
double fRo
nuclear scale parameter
static const double kNucleonMass
double Q2(const Interaction *const i)
virtual Interaction * Summary(void) const
static RandomGen * Instance()
Access instance.
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
void CalculateKin_BergerSehgalFM(GHepRecord *event_rec) const
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
static const double fermi
A simple [min,max] interval for doubles.
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
static const double kPi0Mass
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
Defines the EventGeneratorI interface.
void CalculateKin_AlvarezRuso(GHepRecord *event_rec) const
Generated/set kinematical variables for an event.
void Configure(const Registry &config)
virtual double Weight(void) const
double x(bool selected=false) const
::xsd::cxx::tree::exception< char > exception
Range1D_t YLim(void) const
y limits
void SetKinematics(const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction, Kinematics *kinematics) const
Range1D_t Q2Lim(void) const
Q2 limits.
Abstract class. Provides a data caching mechanism for for concrete implementations of the EventRecord...
A singleton holding random number generator classes. All random number generation in GENIE should tak...
double MaxXSec_AlvarezRuso(const Interaction *in) const
string AsString(void) const
const XSecAlgorithmI * fXSecModel
void CalculateKin_ReinSehgal(GHepRecord *event_rec) const
void ProcessEventRecord(GHepRecord *event_rec) const
Summary information for an interaction.
double fTMax
upper bound for t = (q - p_pi)^2
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)
double fQ2Min
lower bound of integration for Q^2 in Berger-Sehgal Model
virtual void Configure(const Registry &config)
void Sett(double t, bool selected=false)
static const double kASmallNum
void CalculateKin_BergerSehgal(GHepRecord *event_rec) const
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Misc GENIE control constants.
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static float min(const float a, const float b, const float c)
void SetFactor(double factor)
void Setx(double x, bool selected=false)
void UpdateWQ2FromXY(const Interaction *in)
static RunningThreadInfo * Instance(void)
double ComputeMaxXSec(const Interaction *in) const
virtual const AlgId & Id(void) const
Get algorithm ID.
void SetW(double W, bool selected=false)
void SwitchOnFastForward(void)
double Energy(const Interaction *in) const
void SetReason(string reason)
virtual TBits * EventFlags(void) const
static const double kPionMass
Var Sqrt(const Var &v)
Use to take sqrt of a var.
A registry. Provides the container for algorithm configuration parameters.
void Sety(double y, bool selected=false)
~COHKinematicsGenerator()
double MaxXSec_BergerSehgal(const Interaction *in) const
void UpdateXFromQ2Y(const Interaction *in)
bool CheckKinematics(const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction) const
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
double pionMass(const Interaction *in) const
assert(nhit_max >=nhit_nbins)
void SetHadSystP4(const TLorentzVector &p4)
virtual double XSec(void) const
InitialState * InitStatePtr(void) const
TF2 * fEnvelope
2-D envelope used for importance sampling
const InitialState & InitState(void) const
const ProcessInfo & ProcInfo(void) const
void ClearRunningValues(void)
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
static const unsigned int kRjMaxIterations
bool gt(unsigned short int a, unsigned short int b)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const EventGeneratorI * RunningThread(void)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
double ProbeE(RefFrame_t rf) const
GENIE's GHEP MC event record.
double MaxXSec_BergerSehgalFM(const Interaction *in) const
Keep info on the event generation thread currently on charge. This is used so that event generation m...
Root of GENIE utility namespaces.
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
void throwOnTooManyIterations(unsigned int iters, GHepRecord *evrec) const
Initial State information.
static const double kPionMass2