#include <G4MuonRadiativeDecayChannelWithSpin.hh>

Inheritance diagram for G4MuonRadiativeDecayChannelWithSpin:

Public Member Functions
	G4MuonRadiativeDecayChannelWithSpin (const G4String &theParentName, G4double theBR)

virtual	~G4MuonRadiativeDecayChannelWithSpin ()

virtual G4DecayProducts *	DecayIt (G4double)

void	SetPolarization (G4ThreeVector)

const G4ThreeVector &	GetPolarization () const

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="")

virtual	~G4VDecayChannel ()

G4int	operator== (const G4VDecayChannel &right) const

G4int	operator!= (const G4VDecayChannel &right) const

G4int	operator< (const G4VDecayChannel &right) const

virtual G4DecayProducts *	DecayIt (G4double parentMass=-1.0)=0

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 ()

Protected Member Functions
	G4MuonRadiativeDecayChannelWithSpin (const G4MuonRadiativeDecayChannelWithSpin &)

G4MuonRadiativeDecayChannelWithSpin &	operator= (const G4MuonRadiativeDecayChannelWithSpin &)

Protected Member Functions inherited from G4VDecayChannel
	G4VDecayChannel ()

	G4VDecayChannel (const G4VDecayChannel &)

G4VDecayChannel &	operator= (const G4VDecayChannel &)

void	ClearDaughtersName ()

void	FillDaughters ()

void	FillParent ()

Additional Inherited Members
Protected Attributes inherited from G4VDecayChannel
G4String	kinematics_name

G4double	rbranch

G4int	numberOfDaughters

G4String *	parent_name

G4String **	daughters_name

G4ParticleTable *	particletable

G4ParticleDefinition *	parent

G4ParticleDefinition **	daughters

G4double	parent_mass

G4double *	daughters_mass

G4int	verboseLevel

Static Protected Attributes inherited from G4VDecayChannel
static const G4String	noName = " "

Detailed Description

Definition at line 61 of file G4MuonRadiativeDecayChannelWithSpin.hh.

Constructor & Destructor Documentation

◆ G4MuonRadiativeDecayChannelWithSpin() [1/2]

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

Definition at line 59 of file G4MuonRadiativeDecayChannelWithSpin.cc.

         : G4VDecayChannel("Radiative Muon Decay",1),
           parent_polarization()
{
  // 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:: constructor :";
      G4cout << " parent particle is not muon but ";
      G4cout << theParentName << G4endl;
    }
#endif
  }
}

Referenced by G4MuonRadiativeDecayChannelWithSpin().

◆ ~G4MuonRadiativeDecayChannelWithSpin()

G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin ( )

virtual

Definition at line 93 of file G4MuonRadiativeDecayChannelWithSpin.cc.

94{

95}

◆ G4MuonRadiativeDecayChannelWithSpin() [2/2]

G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin ( const G4MuonRadiativeDecayChannelWithSpin & right )

protected

Definition at line 97 of file G4MuonRadiativeDecayChannelWithSpin.cc.

                                                                                                                        :
  G4VDecayChannel(right)
{
  parent_polarization = right.parent_polarization;
}

Member Function Documentation

◆ DecayIt()

G4DecayProducts * G4MuonRadiativeDecayChannelWithSpin::DecayIt ( G4double )

virtual

if(i<10000000)goto leap1:

Implements G4VDecayChannel.

Definition at line 132 of file G4MuonRadiativeDecayChannelWithSpin.cc.

{
 
#ifdef G4VERBOSE
  if (GetVerboseLevel()>1) 
                 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
#endif
 
  if (parent == 0) FillParent();  
  if (daughters == 0) FillDaughters();
 
  // parent mass
  G4double parentmass = parent->GetPDGMass();
 
  G4double EMMU = parentmass;
 
  //daughters'mass
  G4double daughtermass[4]; 
  G4double sumofdaughtermass = 0.0;
  for (G4int index=0; index<4; index++){
    daughtermass[index] = daughters[index]->GetPDGMass();
    sumofdaughtermass += daughtermass[index];
  }
 
  G4double EMASS = daughtermass[0];
 
  //create parent G4DynamicParticle at rest
  G4ThreeVector dummy;
  G4DynamicParticle * parentparticle = 
                               new G4DynamicParticle( parent, dummy, 0.0);
  //create G4Decayproducts
  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
  delete parentparticle;
 
  G4int i = 0;
 
  G4double eps = EMASS/EMMU;
 
  G4double som0, Qsqr, x, y, xx, yy, zz;
  G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
 
  do {
 
//     leap1:
 
     i++;
 
//     leap2:
 
     do {
//
//--------------------------------------------------------------------------
//      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;
//
//-----------------------------------------------------------------------
//
//      Maybe certain restrictions on the kinematical variables:
//
////      if (cthetaE    > 0.01)goto leap2;
////      if (cthetaG    > 0.01)goto leap2;
////      if (std::fabs(x-0.5) > 0.5 )goto leap2;
////      if (std::fabs(y-0.5) > 0.5 )goto leap2;
//
//-----------------------------------------------------------------------
//
//      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;
 
        Qsqr  = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps));
//
//-----------------------------------------------------------------------
//
//      Check the kinematics.
//
     } while ( Qsqr<0.0 || Qsqr>1.0 );
//
////   G4cout << x << " " << y << " " <<  beta << " " << Qsqr << G4endl;
//
//   Do the calculation for -1 muon polarization (i.e. mu+)
//
     G4double Pmu = -1.0;
     if (GetParentName() == "mu-")Pmu = +1.0;
//
//   and for Fronsdal
//
//-----------------------------------------------------------------------
//
     som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE);
//
////     if(som0<0.0){
////       G4cout << " som0 < 0 in Fronsdal " << som0 
////              << " at event " << i << G4endl;
////       G4cout << Pmu << " " << x << " " << y << " " 
////              << cthetaE << " " << cthetaG << " "
////              << cthetaGE << " " << som0 << G4endl;
////     }
//
//-----------------------------------------------------------------------
//
////     G4cout << x << " " << y << " " << som0 << G4endl;
//
//----------------------------------------------------------------------
//
//   Sample the decay rate
//
 
  } while (G4UniformRand()*177.0 > som0);
 
///   if(i<10000000)goto leap1:
//
//-----------------------------------------------------------------------
//
  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);
 
  G4DynamicParticle * daughterparticle0 
    = new G4DynamicParticle( 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);
 
  G4DynamicParticle * daughterparticle1
    = new G4DynamicParticle( 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*(1.0 - (x+y)/2.0);
 
  G4double vmass   = std::sqrt((energy2-
                                (daughtermomentum[0]+daughtermomentum[1]))*
                               (energy2+
                                (daughtermomentum[0]+daughtermomentum[1])));
  G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2;
  beta = -1.0 * std::min(beta,0.99);
 
  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);
 
  G4DynamicParticle * daughterparticle2
    = new G4DynamicParticle( daughters[2], direction2*(vmass/2.));
  G4DynamicParticle * daughterparticle3
    = new G4DynamicParticle( daughters[3], direction2*(-1.0*vmass/2.));
 
  // boost to the muon rest frame
 
  G4ThreeVector direction34(direction0.x()+direction1.x(),
                            direction0.y()+direction1.y(),
                            direction0.z()+direction1.z());
  direction34 = direction34.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;
    products->DumpInfo();
  }
#endif
  return products;
}

◆ GetPolarization()

const G4ThreeVector & G4MuonRadiativeDecayChannelWithSpin::GetPolarization ( ) const

inline

Definition at line 113 of file G4MuonRadiativeDecayChannelWithSpin.hh.

{
  return parent_polarization;
}

◆ operator=()

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

protected

Definition at line 103 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 !=0) 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;
}

◆ SetPolarization()

void G4MuonRadiativeDecayChannelWithSpin::SetPolarization ( G4ThreeVector polar )

inline

Definition at line 107 of file G4MuonRadiativeDecayChannelWithSpin.hh.

{
  parent_polarization = polar;
}

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

Public Member Functions

Protected Member Functions

Additional Inherited Members