From 311955c281f1f06d44119e684c7416adf6c79dd8 Mon Sep 17 00:00:00 2001 From: Igor Kakorin Date: Sat, 10 Jun 2023 19:29:29 +0300 Subject: [PATCH 1/7] fix RESKinematicsGenerator --- src/Framework/Utils/KineUtils.cxx | 6 +++--- .../EventGen/RESKinematicsGenerator.cxx | 19 +++++-------------- 2 files changed, 8 insertions(+), 17 deletions(-) diff --git a/src/Framework/Utils/KineUtils.cxx b/src/Framework/Utils/KineUtils.cxx index 23c8710c1..5b5917134 100644 --- a/src/Framework/Utils/KineUtils.cxx +++ b/src/Framework/Utils/KineUtils.cxx @@ -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 @@ -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; diff --git a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx index f2355f3db..996ec11b4 100644 --- a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx +++ b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx @@ -124,7 +124,6 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const double gW = 0; // current hadronic invariant mass double gQ2 = 0; // current momentum transfer - double gQD2 = 0; // tranformed Q2 to take out dipole form if(fGenerateUniformly) { //-- Generate a W uniformly in the kinematically allowed range. @@ -158,14 +157,9 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const 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 << "]"; + << "Q^2: [" << Q2min << ", " << Q2max << "]"; #endif double mR, gR; if(!interaction->ExclTag().KnownResonance()) { @@ -181,7 +175,7 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const << "(m,g) = (" << mR << ", " << gR << "), max(xsec,W) = (" << xsec_max << ", " << W.max << ")"; #endif - fEnvelope->SetRange(QD2min,W.min,QD2max,W.max); // range + fEnvelope->SetRange(Q2min,W.min,Q2max,W.max); // range fEnvelope->SetParameter(0, mR); // resonance mass fEnvelope->SetParameter(1, gR); // resonance width fEnvelope->SetParameter(2, xsec_max); // max differential xsec @@ -189,10 +183,8 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const }// first pass // Generate W,QD2 using the 2-D envelope as PDF - fEnvelope->GetRandom2(gQD2,gW); + fEnvelope->GetRandom2(gQ2,gW); - // QD2 -> Q2 - gQ2 = utils::kinematics::QD2toQ2(gQD2); } // uniformly over phase space? LOG("RESKinematics", pINFO) << "Trying: W = " << gW << ", Q2 = " << gQ2; @@ -208,9 +200,8 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const if(!fGenerateUniformly) { // unified neutrino / electron scattering - double max = fEnvelope->Eval(gQD2, gW); + double max = fEnvelope->Eval(gQ2, gW); double t = max * rnd->RndKine().Rndm(); - double J = kinematics::Jacobian(interaction,kPSWQ2fE,kPSWQD2fE); this->AssertXSecLimits(interaction, xsec, max); @@ -218,7 +209,7 @@ void RESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const LOG("RESKinematics", pDEBUG) << "xsec= " << xsec << ", J= " << J << ", Rnd= " << t; #endif - accept = (t < J*xsec); + accept = (t < xsec); } // charged lepton or neutrino scattering? else { accept = (xsec>0); From c3d8d2e67e0657e8cf31b7e34a57e57d72a334da Mon Sep 17 00:00:00 2001 From: Igor Kakorin Date: Thu, 15 Jun 2023 15:53:51 +0300 Subject: [PATCH 2/7] fix max searching --- .../EventGen/RESKinematicsGenerator.cxx | 172 ++++++------------ 1 file changed, 52 insertions(+), 120 deletions(-) diff --git a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx index f2355f3db..3d1670673 100644 --- a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx +++ b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx @@ -328,152 +328,84 @@ 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; iq2KinePtr()->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; iq2bKinePtr()->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; iwPhaseSpace(); + 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; iwKinePtr()->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; iq2KinePtr()->SetQ2(Q2); - double xsec = fXSecModel->XSec(interaction, kPSWQ2fE); - 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; iq2bKinePtr()->SetQ2(Q2); - xsec = fXSecModel->XSec(interaction, kPSWQ2fE); - LOG("RESKinematics", pDEBUG) + + for(int iq2=0; iq2KinePtr()->SetQ2(Q2); + double xsec = fXSecModel->XSec(interaction, kPSWQD2fE); + LOG("RESKinematics", pDEBUG) << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec; - max_xsec = TMath::Max(xsec, max_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; iq2bKinePtr()->SetQ2(Q2); + xsec = fXSecModel->XSec(interaction, kPSWQD2fE); + LOG("RESKinematics", pDEBUG) + << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec; + max_xsec = TMath::Max(xsec, max_xsec); + } + break; } - 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 *= ( (EAsString(); - LOG("RESKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec; - LOG("RESKinematics", pDEBUG) << "Computed using " << fXSecModel->Id(); -#endif - return max_xsec; } //___________________________________________________________________________ From 2b12e18be0836126d5b6c2d2e6725907f332432e Mon Sep 17 00:00:00 2001 From: Igor Kakorin Date: Thu, 15 Jun 2023 15:54:24 +0300 Subject: [PATCH 3/7] fix Messenger.xml --- config/Messenger.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/Messenger.xml b/config/Messenger.xml index 52b6ff705..fbe6559b2 100644 --- a/config/Messenger.xml +++ b/config/Messenger.xml @@ -224,7 +224,7 @@ ERROR FATAL - WARN + NOTICE NOTICE NOTICE NOTICE From 454cf401e2e98da4330d9edda5b75d2cba64fe2a Mon Sep 17 00:00:00 2001 From: Igor Kakorin Date: Thu, 15 Jun 2023 15:56:18 +0300 Subject: [PATCH 4/7] fix Messenger.xml --- config/Messenger.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/Messenger.xml b/config/Messenger.xml index 52b6ff705..fbe6559b2 100644 --- a/config/Messenger.xml +++ b/config/Messenger.xml @@ -224,7 +224,7 @@ ERROR FATAL - WARN + NOTICE NOTICE NOTICE NOTICE From 0751e8490c3e89cba7e3a2645c4e798c75e3c986 Mon Sep 17 00:00:00 2001 From: Igor Kakorin Date: Fri, 16 Jun 2023 14:44:42 +0300 Subject: [PATCH 5/7] delete unused code --- .../EventGen/RESKinematicsGenerator.cxx | 151 +++++------------- 1 file changed, 42 insertions(+), 109 deletions(-) diff --git a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx index 996ec11b4..bad89eb9e 100644 --- a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx +++ b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx @@ -319,12 +319,6 @@ 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 { @@ -332,99 +326,43 @@ double RESKinematicsGenerator::ComputeMaxXSec( 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; iq2KinePtr()->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; iq2bKinePtr()->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; iwPhaseSpace(); + 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; iwKinePtr()->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; iq2KinePtr()->SetQ2(Q2); double xsec = fXSecModel->XSec(interaction, kPSWQ2fE); @@ -433,38 +371,33 @@ double RESKinematicsGenerator::ComputeMaxXSec( 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; iq2bKinePtr()->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; iq2bKinePtr()->SetQ2(Q2); + xsec = fXSecModel->XSec(interaction, kPSWQ2fE); + 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 *= ( (EAsString(); - LOG("RESKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec; - LOG("RESKinematics", pDEBUG) << "Computed using " << fXSecModel->Id(); -#endif - return max_xsec; } //___________________________________________________________________________ From 9b86fb2ccfb69d6604010b324fb4a42768f710c4 Mon Sep 17 00:00:00 2001 From: Igor Kakorin Date: Fri, 16 Jun 2023 15:16:38 +0300 Subject: [PATCH 6/7] fix layout --- .../EventGen/RESKinematicsGenerator.cxx | 59 ++++++++++--------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx index 3d1670673..2e2d2c85f 100644 --- a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx +++ b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx @@ -334,7 +334,8 @@ double RESKinematicsGenerator::ComputeMaxXSec( Resonance_t res = interaction->ExclTag().Resonance(); md=res::Mass(res); } - + + // ** 2-D Scan const KPhaseSpace & kps = interaction->PhaseSpace(); Range1D_t rW = kps.WLim(); @@ -371,36 +372,36 @@ double RESKinematicsGenerator::ComputeMaxXSec( for(int iq2=0; iq2KinePtr()->SetQ2(Q2); - 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; iq2bKinePtr()->SetQ2(Q2); - xsec = fXSecModel->XSec(interaction, kPSWQD2fE); - LOG("RESKinematics", pDEBUG) - << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec; - max_xsec = TMath::Max(xsec, max_xsec); - } - break; - } + double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2); + interaction->KinePtr()->SetQ2(Q2); + 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; iq2bKinePtr()->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 - + // Apply safety factor, since value retrieved from the cache might // correspond to a slightly different energy // Apply larger safety factor for smaller energies. From a9a99a9880160413c39507d1d5458fb309331403 Mon Sep 17 00:00:00 2001 From: Igor Kakorin Date: Fri, 16 Jun 2023 23:42:33 +0300 Subject: [PATCH 7/7] delete envelope, max search in QD2 --- .../EventGen/RESKinematicsGenerator.cxx | 118 ++++++------------ 1 file changed, 37 insertions(+), 81 deletions(-) diff --git a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx index 2e2d2c85f..d46a83d81 100644 --- a/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx +++ b/src/Physics/Resonance/EventGen/RESKinematicsGenerator.cxx @@ -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) { LOG("RESKinematics", pWARN) << "No available phase space"; @@ -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 @@ -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; @@ -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 @@ -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