Skip to content

Commit

Permalink
Merge pull request #291 from igor144/fix_RESKinematicsGenerator
Browse files Browse the repository at this point in the history 
fix RESKinematicsGenerator
  • Loading branch information
@mroda88
mroda88 authored Jul 25, 2023
2 parents efac7f5 + ad8c81e commit 971a8b7
Show file tree
Hide file tree
Showing 2 changed files with 83 additions and 194 deletions.
6 changes: 3 additions & 3 deletions src/Framework/Utils/KineUtils.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1373,7 +1373,7 @@ double genie::utils::kinematics::RESImportanceSamplingEnvelope(
double * x, double * par)
{
//-- inputs
double QD2 = x[0]; // QD2 (Q2 transformed to take out the dipole form)
double Q2 = x[0]; // Q2
double W = x[1]; // invariant mass

//-- parameters
Expand Down Expand Up @@ -1424,8 +1424,8 @@ double genie::utils::kinematics::RESImportanceSamplingEnvelope(
}
}

if(QD2<0.5) {
func *= (1 - (0.5-QD2));
if(Q2>controls::kMQD2) {
func *= (0.5 + 1./(1 + Q2/controls::kMQD2));
}

return func;
Expand Down
271 changes: 80 additions & 191 deletions src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,6 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const
// (the physically allowed W's, unless an external cut is imposed)
const KPhaseSpace & kps = interaction->PhaseSpace();
Range1D_t W = kps.Limits(kKVW);
//costas says this is bad if(W.max>1.7) W.max=1.7;
//assert(W.min>=0. && W.min<W.max);

if(W.max <=0 || W.min>=W.max) {
LOG("RESKinematics", pWARN) << "No available phase space";
Expand All @@ -92,8 +90,6 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const

const InitialState & init_state = interaction -> InitState();
double E = init_state.ProbeE(kRfHitNucRest);
// double M = init_state.Tgt().HitNucP4().M();
// double ml = interaction->FSPrimLepton()->Mass();

//-- For the subsequent kinematic selection with the rejection method:
// Calculate the max differential cross section or retrieve it from the
Expand Down Expand Up @@ -126,73 +122,41 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const
double gQ2 = 0; // current momentum transfer
double gQD2 = 0; // tranformed Q2 to take out dipole form

if(fGenerateUniformly) {
if(fGenerateUniformly)
{
//-- Generate a W uniformly in the kinematically allowed range.
// For the generated W, compute the Q2 range and generate a value
// uniformly over that range
gW = W.min + dW * rnd->RndKine().Rndm();
interaction->KinePtr()->SetW(gW);
Range1D_t Q2 = kps.Q2Lim_W();
if(Q2.max<=0. || Q2.min>=Q2.max) continue;
gQ2 = Q2.min + (Q2.max-Q2.min) * rnd->RndKine().Rndm();

interaction->SetBit(kISkipKinematicChk);

} else {


// > neutrino scattering
// Selecting unweighted event kinematics using an importance sampling
// method. Q2 with be transformed to QD2 to take out the dipole form.
// An importance sampling envelope will be constructed for W.
// first pass, configure the sampling envelope
if(iter==1) {
LOG("RESKinematics", pINFO) << "Initializing the sampling envelope";
if(!fEnvelope) {
LOG("RESKinematics", pFATAL) << "Null sampling envelope!";
exit(1);
}
interaction->KinePtr()->SetW(W.min);
Range1D_t Q2 = kps.Q2Lim_W();
double Q2min = -99.;
if (is_em) { Q2min = Q2.min + kASmallNum; }
else { Q2min = 0 + kASmallNum; }
double Q2max = Q2.max - kASmallNum;

// In unweighted mode - use transform that takes out the dipole form
double QD2min = utils::kinematics::Q2toQD2(Q2max);
double QD2max = utils::kinematics::Q2toQD2(Q2min);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG)
<< "Q^2: [" << Q2min << ", " << Q2max << "] => "
<< "QD^2: [" << QD2min << ", " << QD2max << "]";
#endif
double mR, gR;
if(!interaction->ExclTag().KnownResonance()) {
mR = 1.2;
gR = 0.6;
} else {
Resonance_t res = interaction->ExclTag().Resonance();
mR = res::Mass(res);
gR = (E>mR) ? 0.220 : 0.400;
}
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG)
<< "(m,g) = (" << mR << ", " << gR
<< "), max(xsec,W) = (" << xsec_max << ", " << W.max << ")";
#endif
fEnvelope->SetRange(QD2min,W.min,QD2max,W.max); // range
fEnvelope->SetParameter(0, mR); // resonance mass
fEnvelope->SetParameter(1, gR); // resonance width
fEnvelope->SetParameter(2, xsec_max); // max differential xsec
fEnvelope->SetParameter(3, W.max); // kinematically allowed Wmax
}// first pass

// Generate W,QD2 using the 2-D envelope as PDF
fEnvelope->GetRandom2(gQD2,gW);

// QD2 -> Q2
gQ2 = utils::kinematics::QD2toQ2(gQD2);
}
else
{
// neutrino scattering
// Selecting unweighted event kinematics using an importance sampling
// method. Q2 with be transformed to QD2 to take out the dipole form.
interaction->KinePtr()->SetW(W.min);
Range1D_t Q2 = kps.Q2Lim_W();
double Q2min = -99.;
if (is_em)
Q2min = Q2.min + kASmallNum;
else
Q2min = 0 + kASmallNum;
double Q2max = Q2.max - kASmallNum;

// In unweighted mode - use transform that takes out the dipole form
double QD2min = utils::kinematics::Q2toQD2(Q2max);
double QD2max = utils::kinematics::Q2toQD2(Q2min);

gW = W.min + dW * rnd->RndKine().Rndm();
gQD2 = QD2min + (QD2max - QD2min) * rnd->RndKine().Rndm();

// QD2 -> Q2
gQ2 = utils::kinematics::QD2toQ2(gQD2);
} // uniformly over phase space?

LOG("RESKinematics", pINFO) << "Trying: W = " << gW << ", Q2 = " << gQ2;
Expand All @@ -202,25 +166,18 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const
interaction->KinePtr()->SetQ2(gQ2);

//-- Computing cross section for the current kinematics
xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
xsec = fXSecModel->XSec(interaction, kPSWQD2fE);

//-- Decide whether to accept the current kinematics
if(!fGenerateUniformly) {

// unified neutrino / electron scattering
double max = fEnvelope->Eval(gQD2, gW);
double t = max * rnd->RndKine().Rndm();
double J = kinematics::Jacobian(interaction,kPSWQ2fE,kPSWQD2fE);

this->AssertXSecLimits(interaction, xsec, max);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG)
<< "xsec= " << xsec << ", J= " << J << ", Rnd= " << t;
#endif
accept = (t < J*xsec);
} // charged lepton or neutrino scattering?
else {
if(!fGenerateUniformly)
{
// unified neutrino / electron scattering
double t = xsec_max * rnd->RndKine().Rndm();
this->AssertXSecLimits(interaction, xsec, xsec_max);
accept = (t < xsec);
} // charged lepton or neutrino scattering?
else
{
accept = (xsec>0);
} // uniformly over phase space

Expand All @@ -236,7 +193,6 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const
// note: hit nucleon can be off the mass-shell
double gx=-1, gy=-1;
double M = init_state.Tgt().HitNucP4().M();
//double M = init_state.Tgt().HitNucMass();
kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy);

// set the cross section for the selected kinematics
Expand Down Expand Up @@ -328,152 +284,85 @@ double RESKinematicsGenerator::ComputeMaxXSec(
bool is_em = interaction->ProcInfo().IsEM();
double Q2Thres = is_em ? utils::kinematics::electromagnetic::kMinQ2Limit : controls::kMinQ2Limit;

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG) << "Scanning phase space for E= " << E;
#endif
//bool scan1d = (E>1.0);
bool scan1d = false;

double md;
if(!interaction->ExclTag().KnownResonance()) md=1.23;
else {
Resonance_t res = interaction->ExclTag().Resonance();
md=res::Mass(res);
}

if(scan1d) {
// ** 1-D Scan
//
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG)
<< "Will search for max{xsec} along W(=MRes) = " << md;
#endif
// Set W around the value where d^2xsec/dWdQ^2 peaks
interaction->KinePtr()->SetW(md);

const KPhaseSpace & kps = interaction->PhaseSpace();
Range1D_t rQ2 = kps.Q2Lim_W();
if( rQ2.max < Q2Thres || rQ2.min <=0 ) return 0.;

int NQ2 = 25;
int NQ2b = 5;
double logQ2min = TMath::Log(rQ2.min+kASmallNum);
double logQ2max = TMath::Log(rQ2.max-kASmallNum);
double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);

double xseclast = -1;
bool increasing = true;

for(int iq2=0; iq2<NQ2; iq2++) {
double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2);
interaction->KinePtr()->SetQ2(Q2);
double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG)
<< "xsec(W= " << md << ", Q2= " << Q2 << ") = " << xsec;
#endif
max_xsec = TMath::Max(xsec, max_xsec);
increasing = xsec-xseclast>=0;
xseclast=xsec;

// once the cross section stops increasing, I reduce the step size and
// step backwards a little bit to handle cases that the max cross section
// is grossly underestimated (very peaky distribution & large step)
if(!increasing) {
dlogQ2/=NQ2b;
for(int iq2b=0; iq2b<NQ2b; iq2b++) {
Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
if(Q2 < rQ2.min) continue;
interaction->KinePtr()->SetQ2(Q2);
xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG)
<< "xsec(W= " << md << ", Q2= " << Q2 << ") = " << xsec;
#endif
max_xsec = TMath::Max(xsec, max_xsec);
}
break;
}
} // Q2

} else {

// ** 2-D Scan
//
const KPhaseSpace & kps = interaction->PhaseSpace();
Range1D_t rW = kps.WLim();

int NW = 20;
double Wmin = rW.min + kASmallNum;
double Wmax = rW.max - kASmallNum;

Wmax = TMath::Min(Wmax,fWcut);

Wmin = TMath::Max(Wmin, md-.3);
Wmax = TMath::Min(Wmax, md+.3);

if(Wmax-Wmin<0.05) { NW=1; Wmin=Wmax; }

double dW = (NW>1) ? (Wmax-Wmin)/(NW-1) : 0.;

for(int iw=0; iw<NW; iw++) {
// ** 2-D Scan
const KPhaseSpace & kps = interaction->PhaseSpace();
Range1D_t rW = kps.WLim();

int NW = 20;
double Wmin = rW.min + kASmallNum;
double Wmax = rW.max - kASmallNum;

Wmax = TMath::Min(Wmax,fWcut);

Wmin = TMath::Max(Wmin, md-.3);
Wmax = TMath::Min(Wmax, md+.3);

if(Wmax-Wmin<0.05) { NW=1; Wmin=Wmax; }

double dW = (NW>1) ? (Wmax-Wmin)/(NW-1) : 0.;

for(int iw=0; iw<NW; iw++)
{
double W = Wmin + iw*dW;
interaction->KinePtr()->SetW(W);

int NQ2 = 25;
int NQ2b = 4;

Range1D_t rQ2 = kps.Q2Lim_W();
if( rQ2.max < Q2Thres || rQ2.min <=0 ) continue;
if( rQ2.max-rQ2.min<0.02 ) {NQ2=5; NQ2b=3;}

double logQ2min = TMath::Log(rQ2.min+kASmallNum);
double logQ2max = TMath::Log(rQ2.max-kASmallNum);
double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
double xseclast = -1;
bool increasing = true;

for(int iq2=0; iq2<NQ2; iq2++) {

for(int iq2=0; iq2<NQ2; iq2++)
{
double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2);
interaction->KinePtr()->SetQ2(Q2);
double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
double xsec = fXSecModel->XSec(interaction, kPSWQD2fE);
LOG("RESKinematics", pDEBUG)
<< "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec;
max_xsec = TMath::Max(xsec, max_xsec);
increasing = xsec-xseclast>=0;
xseclast=xsec;

// once the cross section stops increasing, I reduce the step size and
// step backwards a little bit to handle cases that the max cross section
// is grossly underestimated (very peaky distribution & large step)
if(!increasing) {
dlogQ2/=NQ2b;
for(int iq2b=0; iq2b<NQ2b; iq2b++) {
Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
if(Q2 < rQ2.min) continue;
interaction->KinePtr()->SetQ2(Q2);
xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
LOG("RESKinematics", pDEBUG)
<< "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec;
max_xsec = TMath::Max(xsec, max_xsec);
}
break;
}
if(!increasing)
{
dlogQ2/=NQ2b;
for(int iq2b=0; iq2b<NQ2b; iq2b++)
{
Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
if(Q2 < rQ2.min) continue;
interaction->KinePtr()->SetQ2(Q2);
xsec = fXSecModel->XSec(interaction, kPSWQD2fE);
LOG("RESKinematics", pDEBUG)
<< "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec;
max_xsec = TMath::Max(xsec, max_xsec);
}
break;
}
} // Q2
}//W
}//2d scan
}//W

// Apply safety factor, since value retrieved from the cache might
// correspond to a slightly different energy
// Apply larger safety factor for smaller energies.
max_xsec *= ( (E<md) ? 2. : fSafetyFactor);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("RESKinematics", pDEBUG) << interaction->AsString();
LOG("RESKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
LOG("RESKinematics", pDEBUG) << "Computed using " << fXSecModel->Id();
#endif

return max_xsec;
}
//___________________________________________________________________________

0 comments on commit 971a8b7

Please sign in to comment.