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..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 @@ -328,12 +284,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 { @@ -341,139 +291,78 @@ 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); + 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, 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, 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 *= ( (EAsString(); - LOG("RESKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec; - LOG("RESKinematics", pDEBUG) << "Computed using " << fXSecModel->Id(); -#endif - return max_xsec; } //___________________________________________________________________________