Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix RESKinematicsGenerator #291

Merged
merged 9 commits into from
Jul 25, 2023
Merged

fix RESKinematicsGenerator #291

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
311955c
fix RESKinematicsGenerator
Jun 10, 2023
c3d8d2e
fix max searching
Jun 15, 2023
2b12e18
fix Messenger.xml
Jun 15, 2023
454cf40
fix Messenger.xml
Jun 15, 2023
0751e84
delete unused code
Jun 16, 2023
9b86fb2
fix layout
Jun 16, 2023
a9a99a9
delete envelope, max search in QD2
Jun 16, 2023
9043ca6
get rid of envelope. max search in QD2
Jun 20, 2023
ad8c81e
Merge branch 'GENIE-MC:master' into fix_RESKinematicsGenerator
idkakorin Jun 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions src/Framework/Utils/KineUtils.cxx
View file
Open in desktop
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
Open in desktop
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;
}
//___________________________________________________________________________