G4MuonRadiativeDecayChannelWithSpin Class Reference

#include <G4MuonRadiativeDecayChannelWithSpin.hh>

Inheritance diagram for G4MuonRadiativeDecayChannelWithSpin:

Public Member Functions
	G4MuonRadiativeDecayChannelWithSpin (const G4String &theParentName, G4double theBR)
	~G4MuonRadiativeDecayChannelWithSpin () override=default
G4DecayProducts *	DecayIt (G4double) override
Public Member Functions inherited from G4VDecayChannel
	G4VDecayChannel (const G4String &aName, G4int Verbose=1)
	G4VDecayChannel (const G4String &aName, const G4String &theParentName, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="", const G4String &theDaughterName4="", const G4String &theDaughterName5="")
virtual	~G4VDecayChannel ()
G4bool	operator== (const G4VDecayChannel &r) const
G4bool	operator!= (const G4VDecayChannel &r) const
G4bool	operator< (const G4VDecayChannel &right) const
const G4String &	GetKinematicsName () const
G4double	GetBR () const
G4int	GetNumberOfDaughters () const
G4ParticleDefinition *	GetParent ()
G4ParticleDefinition *	GetDaughter (G4int anIndex)
G4int	GetAngularMomentum ()
const G4String &	GetParentName () const
const G4String &	GetDaughterName (G4int anIndex) const
G4double	GetParentMass () const
G4double	GetDaughterMass (G4int anIndex) const
void	SetParent (const G4ParticleDefinition *particle_type)
void	SetParent (const G4String &particle_name)
void	SetBR (G4double value)
void	SetNumberOfDaughters (G4int value)
void	SetDaughter (G4int anIndex, const G4ParticleDefinition *particle_type)
void	SetDaughter (G4int anIndex, const G4String &particle_name)
void	SetVerboseLevel (G4int value)
G4int	GetVerboseLevel () const
void	DumpInfo ()
G4double	GetRangeMass () const
void	SetRangeMass (G4double val)
virtual G4bool	IsOKWithParentMass (G4double parentMass)
void	SetPolarization (const G4ThreeVector &)
const G4ThreeVector &	GetPolarization () const

Protected Member Functions
	G4MuonRadiativeDecayChannelWithSpin (const G4MuonRadiativeDecayChannelWithSpin &)=default
G4MuonRadiativeDecayChannelWithSpin &	operator= (const G4MuonRadiativeDecayChannelWithSpin &)
Protected Member Functions inherited from G4VDecayChannel
	G4VDecayChannel ()
	G4VDecayChannel (const G4VDecayChannel &)
G4VDecayChannel &	operator= (const G4VDecayChannel &)
void	ClearDaughtersName ()
void	CheckAndFillDaughters ()
void	CheckAndFillParent ()
G4double	DynamicalMass (G4double massPDG, G4double width, G4double maxDev=1.0) const

Additional Inherited Members
Protected Attributes inherited from G4VDecayChannel
G4String	kinematics_name = ""
G4double	rbranch = 0.0
G4String *	parent_name = nullptr
G4String **	daughters_name = nullptr
G4double	rangeMass = 2.5
G4ThreeVector	parent_polarization
G4ParticleTable *	particletable = nullptr
G4ParticleDefinition *	G4MT_parent = nullptr
G4ParticleDefinition **	G4MT_daughters = nullptr
G4double	G4MT_parent_mass = 0.0
G4double *	G4MT_daughters_mass = nullptr
G4double *	G4MT_daughters_width = nullptr
G4Mutex	daughtersMutex
G4Mutex	parentMutex
G4int	numberOfDaughters = 0
G4int	verboseLevel = 1
Static Protected Attributes inherited from G4VDecayChannel
static const G4String	noName = " "

Detailed Description

Definition at line 53 of file G4MuonRadiativeDecayChannelWithSpin.hh.

Constructor & Destructor Documentation

G4MuonRadiativeDecayChannelWithSpin() [1/2]

G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin	(	const G4String &	theParentName,
		G4double	theBR )

Definition at line 47 of file G4MuonRadiativeDecayChannelWithSpin.cc.

  : G4VDecayChannel("Radiative Muon Decay", 1)
{
  // set names for daughter particles
  if (theParentName == "mu+") {
    SetBR(theBR);
    SetParent("mu+");
    SetNumberOfDaughters(4);
    SetDaughter(0, "e+");
    SetDaughter(1, "gamma");
    SetDaughter(2, "nu_e");
    SetDaughter(3, "anti_nu_mu");
  }
  else if (theParentName == "mu-") {
    SetBR(theBR);
    SetParent("mu-");
    SetNumberOfDaughters(4);
    SetDaughter(0, "e-");
    SetDaughter(1, "gamma");
    SetDaughter(2, "anti_nu_e");
    SetDaughter(3, "nu_mu");
  }
  else {
#ifdef G4VERBOSE
    if (GetVerboseLevel() > 0) {
      G4cout << "G4RadiativeMuonDecayChannel::G4RadiativeMuonDecayChannel():";
      G4cout << " parent particle is not muon but ";
      G4cout << theParentName << G4endl;
    }
#endif
  }
}

Referenced by G4MuonRadiativeDecayChannelWithSpin().

~G4MuonRadiativeDecayChannelWithSpin()

G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin ( )

overridedefault

G4MuonRadiativeDecayChannelWithSpin() [2/2]

G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin ( const G4MuonRadiativeDecayChannelWithSpin & )

protecteddefault

Member Function Documentation

DecayIt()

G4DecayProducts * G4MuonRadiativeDecayChannelWithSpin::DecayIt ( G4double )

overridevirtual

Implements G4VDecayChannel.

Definition at line 110 of file G4MuonRadiativeDecayChannelWithSpin.cc.

{
#ifdef G4VERBOSE
  if (GetVerboseLevel() > 1) G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt()";
#endif
 
  CheckAndFillParent();
  CheckAndFillDaughters();
 
  // parent mass
  G4double parentmass = G4MT_parent->GetPDGMass();
 
  G4double EMMU = parentmass;
 
  // daughters'mass
  G4double daughtermass[4];
  // G4double sumofdaughtermass = 0.0;
  for (G4int index = 0; index < 4; ++index) {
    daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
    // sumofdaughtermass += daughtermass[index];
  }
 
  G4double EMASS = daughtermass[0];
 
  // create parent G4DynamicParticle at rest
  G4ThreeVector dummy;
  auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
 
  // create G4Decayproducts
  auto products = new G4DecayProducts(*parentparticle);
  delete parentparticle;
 
  G4double eps = EMASS / EMMU;
 
  G4double som0, x, y, xx, yy, zz;
  G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
  const std::size_t MAX_LOOP = 10000;
 
  for (std::size_t loop_counter1 = 0; loop_counter1 < MAX_LOOP; ++loop_counter1) {
    for (std::size_t loop_counter2 = 0; loop_counter2 < MAX_LOOP; ++loop_counter2) {
      // -------------------------------------------------------------------
      // Build two vectors of random length and random direction, for the
      // positron and the photon.
      // x/y is the length of the vector, xx, yy and zz the components,
      // phi is the azimutal angle, theta the polar angle.
      // -------------------------------------------------------------------
 
      // For the positron
      //
      x = G4UniformRand();
 
      rn3dim(xx, yy, zz, x);
 
      if (std::fabs((xx * xx) + (yy * yy) + (zz * zz) - (x * x)) > 0.001) {
        G4cout << "Norm of x not correct" << G4endl;
      }
 
      phiE = atan4(xx, yy);
      cthetaE = zz / x;
      G4double sthetaE = std::sqrt((xx * xx) + (yy * yy)) / x;
 
      // What you get:
      //
      // x       = positron energy
      // phiE    = azimutal angle of positron momentum
      // cthetaE = cosine of polar angle of positron momentum
      // sthetaE = sine of polar angle of positron momentum
      //
      //// G4cout << " x, xx, yy, zz " << x  << " " << xx << " "
      ////                             << yy << " " << zz << G4endl;
      //// G4cout << " phiE, cthetaE, sthetaE " << phiE    << " "
      ////                                      << cthetaE << " "
      ////                                      << sthetaE << " " << G4endl;
 
      // For the photon
      //
      y = G4UniformRand();
 
      rn3dim(xx, yy, zz, y);
 
      if (std::fabs((xx * xx) + (yy * yy) + (zz * zz) - (y * y)) > 0.001) {
        G4cout << " Norm of y not correct " << G4endl;
      }
 
      phiG = atan4(xx, yy);
      cthetaG = zz / y;
      G4double sthetaG = std::sqrt((xx * xx) + (yy * yy)) / y;
 
      // What you get:
      //
      // y       = photon energy
      // phiG    = azimutal angle of photon momentum
      // cthetaG = cosine of polar angle of photon momentum
      // sthetaG = sine of polar angle of photon momentum
      //
      //// G4cout << " y, xx, yy, zz " << y  << " " << xx << " "
      ////                             << yy << " " << zz << G4endl;
      //// G4cout << " phiG, cthetaG, sthetaG " << phiG    << " "
      ////                                      << cthetaG << " "
      ////                                      << sthetaG << " " << G4endl;
 
      //      Calculate the angle between positron and photon (cosine)
      //
      cthetaGE = cthetaE * cthetaG + sthetaE * sthetaG * std::cos(phiE - phiG);
 
      //// G4cout << x << " " << cthetaE << " " << sthetaE << " "
      ////        << y << " " << cthetaG << " " << sthetaG << " "
      ////        << cthetaGE
 
      G4double term0 = eps * eps;
      G4double term1 = x * ((1.0 - eps) * (1.0 - eps)) + 2.0 * eps;
      G4double beta =
        std::sqrt(x * ((1.0 - eps) * (1.0 - eps)) * (x * ((1.0 - eps) * (1.0 - eps)) + 4.0 * eps))
        / term1;
      G4double delta = 1.0 - beta * cthetaGE;
 
      G4double term3 = y * (1.0 - (eps * eps));
      G4double term6 = term1 * delta * term3;
 
      G4double Qsqr = (1.0 - term1 - term3 + term0 + 0.5 * term6) / ((1.0 - eps) * (1.0 - eps));
 
      // Check the kinematics.
      //
      if (Qsqr >= 0.0 && Qsqr <= 1.0) break;
 
    }  // end loop count
 
    // Do the calculation for -1 muon polarization (i.e. mu+)
    //
    G4double Pmu = -1.0;
    if (GetParentName() == "mu-") {
      Pmu = +1.0;
    }
 
    som0 = fron(Pmu, x, y, cthetaE, cthetaG, cthetaGE);
 
    // Sample the decay rate
    //
    if (G4UniformRand() * 177.0 <= som0) break;
  }
 
  G4double E = EMMU / 2. * (x * ((1. - eps) * (1. - eps)) + 2. * eps);
  G4double G = EMMU / 2. * y * (1. - eps * eps);
 
  if (E < EMASS) E = EMASS;
 
  // calculate daughter momentum
  G4double daughtermomentum[4];
 
  daughtermomentum[0] = std::sqrt(E * E - EMASS * EMASS);
 
  G4double sthetaE = std::sqrt(1. - cthetaE * cthetaE);
  G4double cphiE = std::cos(phiE);
  G4double sphiE = std::sin(phiE);
 
  // Coordinates of the decay positron with respect to the muon spin
 
  G4double px = sthetaE * cphiE;
  G4double py = sthetaE * sphiE;
  G4double pz = cthetaE;
 
  G4ThreeVector direction0(px, py, pz);
 
  direction0.rotateUz(parent_polarization);
 
  auto daughterparticle0 =
    new G4DynamicParticle(G4MT_daughters[0], daughtermomentum[0] * direction0);
 
  products->PushProducts(daughterparticle0);
 
  daughtermomentum[1] = G;
 
  G4double sthetaG = std::sqrt(1. - cthetaG * cthetaG);
  G4double cphiG = std::cos(phiG);
  G4double sphiG = std::sin(phiG);
 
  // Coordinates of the decay gamma with respect to the muon spin
 
  px = sthetaG * cphiG;
  py = sthetaG * sphiG;
  pz = cthetaG;
 
  G4ThreeVector direction1(px, py, pz);
 
  direction1.rotateUz(parent_polarization);
 
  auto daughterparticle1 =
    new G4DynamicParticle(G4MT_daughters[1], daughtermomentum[1] * direction1);
 
  products->PushProducts(daughterparticle1);
 
  // daughter 3 ,4 (neutrinos)
  // create neutrinos in the C.M frame of two neutrinos
 
  G4double energy2 = parentmass - E - G;
 
  G4ThreeVector P34 = -1. * (daughtermomentum[0] * direction0 + daughtermomentum[1] * direction1);
  G4double vmass2 = energy2 * energy2 - P34.mag2();
  G4double vmass = std::sqrt(vmass2);
 
  G4double costhetan = 2. * G4UniformRand() - 1.0;
  G4double sinthetan = std::sqrt((1.0 - costhetan) * (1.0 + costhetan));
  G4double phin = twopi * G4UniformRand() * rad;
  G4double sinphin = std::sin(phin);
  G4double cosphin = std::cos(phin);
 
  G4ThreeVector direction2(sinthetan * cosphin, sinthetan * sinphin, costhetan);
 
  auto daughterparticle2 = new G4DynamicParticle(G4MT_daughters[2], direction2 * (vmass / 2.));
  auto daughterparticle3 =
    new G4DynamicParticle(G4MT_daughters[3], direction2 * (-1.0 * vmass / 2.));
 
  // boost to the muon rest frame
  G4double beta = P34.mag() / energy2;
  G4ThreeVector direction34 = P34.unit();
 
  G4LorentzVector p4 = daughterparticle2->Get4Momentum();
  p4.boost(direction34.x() * beta, direction34.y() * beta, direction34.z() * beta);
  daughterparticle2->Set4Momentum(p4);
 
  p4 = daughterparticle3->Get4Momentum();
  p4.boost(direction34.x() * beta, direction34.y() * beta, direction34.z() * beta);
  daughterparticle3->Set4Momentum(p4);
 
  products->PushProducts(daughterparticle2);
  products->PushProducts(daughterparticle3);
 
  daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
  daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
 
  // output message
#ifdef G4VERBOSE
  if (GetVerboseLevel() > 1) {
    G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt() -";
    G4cout << " create decay products in rest frame " << G4endl;
    G4double TT = daughterparticle0->GetTotalEnergy() + daughterparticle1->GetTotalEnergy()
                  + daughterparticle2->GetTotalEnergy() + daughterparticle3->GetTotalEnergy();
    G4cout << "e    :" << daughterparticle0->GetTotalEnergy() / MeV << G4endl;
    G4cout << "gamma:" << daughterparticle1->GetTotalEnergy() / MeV << G4endl;
    G4cout << "nu2  :" << daughterparticle2->GetTotalEnergy() / MeV << G4endl;
    G4cout << "nu2  :" << daughterparticle3->GetTotalEnergy() / MeV << G4endl;
    G4cout << "total:" << (TT - parentmass) / keV << G4endl;
    if (GetVerboseLevel() > 1) {
      products->DumpInfo();
    }
  }
#endif
 
  return products;
}

operator=()

G4MuonRadiativeDecayChannelWithSpin & G4MuonRadiativeDecayChannelWithSpin::operator= ( const G4MuonRadiativeDecayChannelWithSpin & right )

protected

Definition at line 82 of file G4MuonRadiativeDecayChannelWithSpin.cc.

{
  if (this != &right) {
    kinematics_name = right.kinematics_name;
    verboseLevel = right.verboseLevel;
    rbranch = right.rbranch;
 
    // copy parent name
    parent_name = new G4String(*right.parent_name);
 
    // clear daughters_name array
    ClearDaughtersName();
 
    // recreate array
    numberOfDaughters = right.numberOfDaughters;
    if (numberOfDaughters > 0) {
      if (daughters_name != nullptr) ClearDaughtersName();
      daughters_name = new G4String*[numberOfDaughters];
      // copy daughters name
      for (G4int index = 0; index < numberOfDaughters; ++index) {
        daughters_name[index] = new G4String(*right.daughters_name[index]);
      }
    }
    parent_polarization = right.parent_polarization;
  }
  return *this;
}

---

The documentation for this class was generated from the following files:

• G4MuonRadiativeDecayChannelWithSpin.hh
• G4MuonRadiativeDecayChannelWithSpin.cc