Skip to content

RMC phase-space approximation #1879

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

Open
michaelmackenzie wants to merge 1 commit into
base: main
Choose a base branch
from
Open

RMC phase-space approximation #1879

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
5 changes: 3 additions & 2 deletions GlobalConstantsService/data/globalConstants_01.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,13 +139,14 @@ double physicsParams.C.capture.RMCInternalRate = 0.0069; //measurement of intern

//Ce+ endpoint energy
double physicsParams.Al.ePlusEndpointEnergy = 92.32;
double physicsParams.Ti.ePlusEndpointEnergy = 98.89;

// parameters for RMC spectrum shape
double physicsParams.Al.RMCbindingEnergyFit = 0.464;
double physicsParams.Al.RMCrecoilEnergyFit = 0.220;
double physicsParams.Al.RMCdeltaMassFit = 3.121;


double physicsParams.Al.RMCKMaxKnockout.0 = 101.87; // 0 nucleon knockout
double physicsParams.Al.RMCKMaxKnockout.1 = 95.45; // 1 neutron knockout

// Shanker interpolation coefficients - Phys. Rev. D 25, 1847 (1982)
// - The Shanker formula includes one D, E and F coefficient for a
Expand Down
9 changes: 9 additions & 0 deletions GlobalConstantsService/inc/PhysicsParams.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,13 @@ namespace mu2e
double getRMCdeltaMassFit(targetMat material = "") const {
return doubleOrThrow(_RMCdeltaMassFit,material,"RMCdeltaMassFit");
}
double getRMCKMaxKnockout(targetMat material = "", int nKnockout = 0) const {
if(nKnockout == 0)
return doubleOrThrow(_RMCKMaxKnockout0n,material,std::string("RMCKMaxKnockout.0"));
if(nKnockout == 1)
return doubleOrThrow(_RMCKMaxKnockout1n,material,std::string("RMCKMaxKnockout.1"));
throw cet::exception("PhysicsParams") << "Greater than 1 knockout RMC not yet implemented!";
}

PhysicsParams( SimpleConfig const& config );

Expand Down Expand Up @@ -220,6 +227,8 @@ namespace mu2e
std::map<targetMat, double> _RMCbindingEnergyFit;
std::map<targetMat, double> _RMCrecoilEnergyFit;
std::map<targetMat, double> _RMCdeltaMassFit;
std::map<targetMat, double> _RMCKMaxKnockout0n;
std::map<targetMat, double> _RMCKMaxKnockout1n;


inline targetMat checkMaterial( const targetMat& material ) const {
Expand Down
4 changes: 4 additions & 0 deletions GlobalConstantsService/src/PhysicsParams.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,10 @@ namespace mu2e {
if(config.hasName(name)) _RMCrecoilEnergyFit[material] = config.getDouble(name);
name = "physicsParams."+material+".RMCdeltaMassFit";
if(config.hasName(name)) _RMCdeltaMassFit[material] = config.getDouble(name);
name = "physicsParams."+material+".RMCKMaxKnockout.0";
if(config.hasName(name)) _RMCKMaxKnockout0n[material] = config.getDouble(name);
name = "physicsParams."+material+".RMCKMaxKnockout.1";
if(config.hasName(name)) _RMCKMaxKnockout1n[material] = config.getDouble(name);

}

Expand Down
6 changes: 5 additions & 1 deletion Mu2eUtilities/inc/MuonCaptureSpectrum.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ namespace mu2e {

public:

enum enum_type { Flat , RMC };
enum enum_type { Flat , RMC , PhaseSpace };
enum enum_type_2D { Flat2D, KrollWadaJoseph };

MuonCaptureSpectrum(){}
Expand All @@ -32,6 +32,8 @@ namespace mu2e {

MuonCaptureSpectrum(bool kMaxUserSet, double kMaxUser, double kMaxMax,
CLHEP::RandFlat* randFlat = 0, RandomUnitSphere* randomUnitSphere = 0);
MuonCaptureSpectrum(double kMax, int nKnockout,
CLHEP::RandFlat* randFlat = 0, RandomUnitSphere* randomUnitSphere = 0);

~MuonCaptureSpectrum(){}

Expand All @@ -45,6 +47,7 @@ namespace mu2e {

double getFlat (double e, double x = 0., double y = 0.) const ;
double getRMCSpectrum(double e, bool kMaxUserSet, double kMaxUser, double kMaxMax) const;
double getRMCPhaseSpaceSpectrum(double e, double kMax, int nKnockout) const;

void getElecPosiVectors(double energy, CLHEP::HepLorentzVector& mome, CLHEP::HepLorentzVector& momp) const;

Expand All @@ -60,6 +63,7 @@ namespace mu2e {
bool _kMaxUserSet;
double _kMaxUser;
double _kMaxMax;
int _nKnockout; // for phase space approximation

CLHEP::RandFlat* _rnFlat;
RandomUnitSphere* _rnUnitSphere;
Expand Down
27 changes: 17 additions & 10 deletions Mu2eUtilities/src/BinnedSpectrum.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,23 +62,30 @@ namespace mu2e {
this->initialize<ConversionSpectrum>(elow,ehi,bin,ehi,bin);
}else if (spectrumShape == "ejectedProtons") {
// should be kinetic energy
double elow = 0.;
double elow = psphys.get<double>("elow", 0.);
// cut off at muon mass
double ehi = GlobalConstantsHandle<ParticleDataList>()->particle(PDGCode::mu_minus).mass();
double ehi = psphys.get<double>("ehi", GlobalConstantsHandle<ParticleDataList>()->particle(PDGCode::mu_minus).mass());
double bin = (ehi - elow)/psphys.get<unsigned>("nbins");
this->initialize<EjectedProtonSpectrum>(elow, ehi, bin);
}else if (spectrumShape == "RMC") {
double elow = psphys.get<double>("elow");
double ehi = psphys.get<double>("ehi");
double res = psphys.get<double>("spectrumResolution");
bool kMaxUserSet = psphys.get<bool> ("kMaxUserSet",false);
double kMaxUser = psphys.get<double>("kMaxUser",0);
const double bindingEnergyFit = GlobalConstantsHandle<PhysicsParams>()->getRMCbindingEnergyFit("Al");
const double recoilEnergyFit = GlobalConstantsHandle<PhysicsParams>()->getRMCrecoilEnergyFit("Al");
const double deltaMassFit = GlobalConstantsHandle<PhysicsParams>()->getRMCdeltaMassFit("Al");
const double mmu = GlobalConstantsHandle<ParticleDataList>()->particle(PDGCode::mu_minus).mass();
const double kMaxMax =mmu - bindingEnergyFit - recoilEnergyFit - deltaMassFit;
this->initialize<MuonCaptureSpectrum>(elow, ehi, res, kMaxUserSet, kMaxUser, kMaxMax);
const bool phase_space = psphys.get<bool>("phaseSpaceShape", false);
if(phase_space) {
const int nKnockout = psphys.get<int>("nKnockout");
const double kMax = GlobalConstantsHandle<PhysicsParams>()->getRMCKMaxKnockout("Al", nKnockout);
this->initialize<MuonCaptureSpectrum>(elow, ehi, res, kMax, nKnockout);
} else { // Closure approximation
const double bindingEnergyFit = GlobalConstantsHandle<PhysicsParams>()->getRMCbindingEnergyFit("Al");
const double recoilEnergyFit = GlobalConstantsHandle<PhysicsParams>()->getRMCrecoilEnergyFit("Al");
const double deltaMassFit = GlobalConstantsHandle<PhysicsParams>()->getRMCdeltaMassFit("Al");
const double mmu = GlobalConstantsHandle<ParticleDataList>()->particle(PDGCode::mu_minus).mass();
const double kMaxMax =mmu - bindingEnergyFit - recoilEnergyFit - deltaMassFit;
bool kMaxUserSet = psphys.get<bool> ("kMaxUserSet",false);
double kMaxUser = psphys.get<double>("kMaxUser",0);
this->initialize<MuonCaptureSpectrum>(elow, ehi, res, kMaxUserSet, kMaxUser, kMaxMax);
}
}else if (spectrumShape == "Bistirlich") {
double elow = psphys.get<double>("elow");
double ehi = psphys.get<double>("ehi");
Expand Down
42 changes: 41 additions & 1 deletion Mu2eUtilities/src/MuonCaptureSpectrum.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ namespace mu2e {
_kMaxUserSet (false ),
_kMaxUser (0. ),
_kMaxMax (0. ),
_nKnockout (0 ),
_rnFlat (rnFlat ),
_rnUnitSphere(rnUnitSphere )
{
Expand All @@ -55,6 +56,7 @@ namespace mu2e {
_kMaxUserSet (kMaxUserSet ),
_kMaxUser (kMaxUser ),
_kMaxMax (kMaxMax ),
_nKnockout (0 ),
_rnFlat (rnFlat ),
_rnUnitSphere(rnUnitSphere )
{
Expand All @@ -65,6 +67,25 @@ namespace mu2e {
_MN = GlobalConstantsHandle<PhysicsParams>()->getAtomicMass("Al");
}

MuonCaptureSpectrum::MuonCaptureSpectrum(double kMax, int nKnockout,
CLHEP::RandFlat* rnFlat, RandomUnitSphere* rnUnitSphere):
_spectrum (PhaseSpace ),
_spectrum2D (KrollWadaJoseph ),
_kMaxUserSet (true ),
_kMaxUser (kMax ),
_kMaxMax (kMax ),
_nKnockout (nKnockout ),
_rnFlat (rnFlat ),
_rnUnitSphere(rnUnitSphere )
{
GlobalConstantsHandle<ParticleDataList> pdt;

_me = pdt->particle(PDGCode::e_minus ).mass();
_mmu = pdt->particle(PDGCode::mu_minus).mass();
_MN = GlobalConstantsHandle<PhysicsParams>()->getAtomicMass("Al");

}

double MuonCaptureSpectrum::getWeight(double E) const {

double weight(0.);
Expand All @@ -73,7 +94,8 @@ namespace mu2e {
else if ( _spectrum == RMC ) weight = getRMCSpectrum( E , _kMaxUserSet, _kMaxUser, _kMaxMax);
*/

weight = getRMCSpectrum(E, _kMaxUserSet,_kMaxUser,_kMaxMax);
if (_spectrum == RMC ) weight = getRMCSpectrum(E, _kMaxUserSet,_kMaxUser,_kMaxMax);
else if(_spectrum == PhaseSpace) weight = getRMCPhaseSpaceSpectrum(E, _kMaxUser, _nKnockout);

// std::cout << "spectrum is " << _spectrum << std::endl;

Expand Down Expand Up @@ -130,6 +152,24 @@ namespace mu2e {
}


//=======================================================
// Phase-space approximation to the photon energy spectrum
// -
// - see doc-db 57178
//=======================================================

double MuonCaptureSpectrum::getRMCPhaseSpaceSpectrum(double e, double kMax, int nKnockout) const {
if ( e >= kMax ) return 0.;
if ( e <= 0. ) return 0.;

const double x = e/kMax;
const double power = 2. + 1.5*nKnockout;
const double norm = (power + 1.) * (power + 2.) / kMax;
const double weight = norm * x * std::pow(1. - x, power);
return weight;
}


//=======================================================
// Analytic expression for electron/positron energies via
// - Kroll and Wada, Phys. Rev. 98, 1355 (1955)
Expand Down