#include <G4MuonDecayChannelWithSpin.hh>

Inheritance diagram for G4MuonDecayChannelWithSpin:

Collaboration diagram for G4MuonDecayChannelWithSpin:

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

virtual	~G4MuonDecayChannelWithSpin ()

virtual G4DecayProducts *	DecayIt (G4double)

Public Member Functions inherited from G4MuonDecayChannel
	G4MuonDecayChannel (const G4String &theParentName, G4double theBR)

virtual	~G4MuonDecayChannel ()

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

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
	G4MuonDecayChannelWithSpin (const G4MuonDecayChannelWithSpin &)

G4MuonDecayChannelWithSpin &	operator= (const G4MuonDecayChannelWithSpin &)

Protected Member Functions inherited from G4MuonDecayChannel
	G4MuonDecayChannel (const G4MuonDecayChannel &)

G4MuonDecayChannel &	operator= (const G4MuonDecayChannel &)

	G4MuonDecayChannel ()

Protected Member Functions inherited from G4VDecayChannel
void	ClearDaughtersName ()

void	CheckAndFillDaughters ()

void	CheckAndFillParent ()

	G4VDecayChannel ()

	G4VDecayChannel (const G4VDecayChannel &)

G4VDecayChannel &	operator= (const G4VDecayChannel &)

G4double	DynamicalMass (G4double massPDG, G4double width, G4double maxDev=+1.) const

Private Member Functions
	G4MuonDecayChannelWithSpin ()

G4double	F_c (G4double x, G4double x0, G4double omega)

G4double	F_theta (G4double x, G4double x0, G4double omega)

G4double	R_c (G4double x, G4double omega)

Additional Inherited Members
Protected Attributes inherited from G4VDecayChannel
G4String	kinematics_name

G4double	rbranch

G4int	numberOfDaughters

G4String *	parent_name

G4String **	daughters_name

G4double	rangeMass

G4ThreeVector	parent_polarization

G4ParticleTable *	particletable

G4int	verboseLevel

G4ParticleDefinition *	G4MT_parent

G4ParticleDefinition **	G4MT_daughters

G4double	G4MT_parent_mass

G4double *	G4MT_daughters_mass

G4double *	G4MT_daughters_width

G4Mutex	daughtersMutex

G4Mutex	parentMutex

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

Detailed Description

Definition at line 50 of file G4MuonDecayChannelWithSpin.hh.

Constructor & Destructor Documentation

◆ G4MuonDecayChannelWithSpin() [1/3]

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

Definition at line 56 of file G4MuonDecayChannelWithSpin.cc.

   : G4MuonDecayChannel(theParentName,theBR)
 {
 }

◆ ~G4MuonDecayChannelWithSpin()

G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin ( )

virtual

Definition at line 62 of file G4MuonDecayChannelWithSpin.cc.

63 {

64 }

◆ G4MuonDecayChannelWithSpin() [2/3]

G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin ( const G4MuonDecayChannelWithSpin & right )

protected

Definition at line 66 of file G4MuonDecayChannelWithSpin.cc.

                                                                                              :
   G4MuonDecayChannel(right)
 {
 }

◆ G4MuonDecayChannelWithSpin() [3/3]

G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin ( )

private

Definition at line 51 of file G4MuonDecayChannelWithSpin.cc.

   : G4MuonDecayChannel()
 {
 }

Member Function Documentation

◆ DecayIt()

G4DecayProducts * G4MuonDecayChannelWithSpin::DecayIt ( G4double )

virtual

Reimplemented from G4MuonDecayChannel.

Definition at line 99 of file G4MuonDecayChannelWithSpin.cc.

 {
   // This version assumes V-A coupling with 1st order radiative correctons,
   //              the standard model Michel parameter values, but
   //              gives incorrect energy spectrum for neutrinos
 
 #ifdef G4VERBOSE
   if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
 #endif
 
   CheckAndFillParent();
   CheckAndFillDaughters();
 
   // parent mass
   G4double parentmass = G4MT_parent->GetPDGMass();
 
   G4double EMMU = parentmass;
 
   //daughters'mass
   G4double daughtermass[3]; 
   G4double sumofdaughtermass = 0.0;
   for (G4int index=0; index<3; index++){
     daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
     sumofdaughtermass += daughtermass[index];
   }
 
   G4double EMASS = daughtermass[0];
 
   //create parent G4DynamicParticle at rest
   G4ThreeVector dummy;
   G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
   //create G4Decayproducts
   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
   delete parentparticle;
 
   // calcurate electron energy
 
   G4double michel_rho   = 0.75; //Standard Model Michel rho
   G4double michel_delta = 0.75; //Standard Model Michel delta
   G4double michel_xsi   = 1.00; //Standard Model Michel xsi
   G4double michel_eta   = 0.00; //Standard Model eta
 
   G4double rndm, x, ctheta;
 
   G4double FG; 
   G4double FG_max = 2.00;
 
   G4double W_mue  = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
   G4double x0     =           EMASS/W_mue;
 
   G4double x0_squared = x0*x0;
   
   // ***************************************************
   //     x0 <= x <= 1.   and   -1 <= y <= 1
   //
   //     F(x,y) = f(x)*g(x,y);   g(x,y) = 1.+g(x)*y
   // ***************************************************
 
   // ***** sampling F(x,y) directly (brute force) *****
 
   const size_t MAX_LOOP=10000;
   for (size_t loop_count =0; loop_count <MAX_LOOP; ++loop_count){
 
     // Sample the positron energy by sampling from F
 
     rndm = G4UniformRand();
 
     x = x0 + rndm*(1.-x0);
 
     G4double x_squared = x*x;
 
     G4double F_IS, F_AS, G_IS, G_AS;
 
     F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
     F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
 
     G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
     G_IS = G_IS + michel_eta*(1.-x)*x0;
 
     G_AS = 3.*(michel_xsi-1.)*(1.-x);
     G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
     G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
 
     F_IS = F_IS + G_IS;
     F_AS = F_AS + G_AS;
 
     // *** Radiative Corrections ***
     const G4double omega =  std::log(EMMU/EMASS);
     G4double R_IS = F_c(x,x0,omega);
 
     G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
 
     // *** Radiative Corrections ***
 
     G4double R_AS = F_theta(x,x0,omega);
 
     rndm = G4UniformRand();
 
     ctheta = 2.*rndm-1.;
 
     G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
 
     FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
 
     if(FG>FG_max){
       G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl;
       FG_max = FG;
     }
 
     rndm = G4UniformRand();
 
     if (FG >= rndm*FG_max) break;
   }
 
   G4double energy = x * W_mue;
 
   rndm = G4UniformRand();
 
   G4double phi = twopi * rndm;
 
   if(energy < EMASS) energy = EMASS;
 
   // calculate daughter momentum
   G4double daughtermomentum[3];
 
   daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
 
   G4double stheta = std::sqrt(1.-ctheta*ctheta);
   G4double cphi = std::cos(phi);
   G4double sphi = std::sin(phi);
 
   //Coordinates of the decay positron with respect to the muon spin
 
   G4double px = stheta*cphi;
   G4double py = stheta*sphi;
   G4double pz = ctheta;
 
   G4ThreeVector direction0(px,py,pz);
 
   direction0.rotateUz(parent_polarization);
 
   G4DynamicParticle * daughterparticle0 
     = new G4DynamicParticle( G4MT_daughters[0], daughtermomentum[0]*direction0);
 
   products->PushProducts(daughterparticle0);
 
 
   // daughter 1 ,2 (neutrinos)
   // create neutrinos in the C.M frame of two neutrinos
   G4double energy2 = parentmass*(1.0 - x/2.0); 
   G4double vmass   = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
   G4double beta = -1.0*daughtermomentum[0]/energy2;
   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 direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
   G4DynamicParticle * daughterparticle1 
     = new G4DynamicParticle( G4MT_daughters[1], direction1*(vmass/2.));
   G4DynamicParticle * daughterparticle2
     = new G4DynamicParticle( G4MT_daughters[2], direction1*(-1.0*vmass/2.));
 
   // boost to the muon rest frame
   G4LorentzVector p4;
   p4 = daughterparticle1->Get4Momentum();
   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
   daughterparticle1->Set4Momentum(p4);
   p4 = daughterparticle2->Get4Momentum();
   p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
   daughterparticle2->Set4Momentum(p4);
   products->PushProducts(daughterparticle1);
   products->PushProducts(daughterparticle2);
   daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
   daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
 
   // output message
 #ifdef G4VERBOSE
   if (GetVerboseLevel()>1) {
     G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
     G4cout << "  create decay products in rest frame " <<G4endl;
     products->DumpInfo();
   }
 #endif
   return products;
 }

Here is the call graph for this function:

◆ F_c()

G4double G4MuonDecayChannelWithSpin::F_c	(	G4double	x,
		G4double	x0,
		G4double	omega
	)

inlineprivate

Definition at line 87 of file G4MuonDecayChannelWithSpin.hh.

 {
 
   G4double f_c;
 
   f_c = (5.+17.*x-34.*x*x)*(omega+std::log(x))-22.*x+34.*x*x;
   f_c = (1.-x)/(3.*x*x)*f_c;
   f_c = (6.-4.*x)*R_c(x,omega)+(6.-6.*x)*std::log(x) + f_c;
   f_c = (CLHEP::fine_structure_const/CLHEP::twopi) * (x*x-x0*x0) * f_c;
 
   return f_c;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ F_theta()

G4double G4MuonDecayChannelWithSpin::F_theta	(	G4double	x,
		G4double	x0,
		G4double	omega
	)

inlineprivate

Definition at line 100 of file G4MuonDecayChannelWithSpin.hh.

 {
   G4double f_theta;
 
   f_theta = (1.+x+34*x*x)*(omega+std::log(x))+3.-7.*x-32.*x*x;
   f_theta = f_theta + ((4.*(1.-x)*(1.-x))/x)*std::log(1.-x);
   f_theta = (1.-x)/(3.*x*x) * f_theta;
   f_theta = (2.-4.*x)*R_c(x,omega)+(2.-6.*x)*std::log(x)-f_theta;
   f_theta = (CLHEP::fine_structure_const/CLHEP::twopi) * (x*x-x0*x0) * f_theta;
 
   return f_theta;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ operator=()

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

protected

Definition at line 71 of file G4MuonDecayChannelWithSpin.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]);
       }
     }
   }
   return *this;
 }

Here is the call graph for this function:

◆ R_c()

G4double G4MuonDecayChannelWithSpin::R_c	(	G4double	x,
		G4double	omega
	)

private

Definition at line 287 of file G4MuonDecayChannelWithSpin.cc.

                                                                  {
 
   G4int n_max = (int)(100.*x);
 
   if(n_max<10)n_max=10;
 
   G4double L2 = 0.0;
 
   for(G4int n=1; n<=n_max; n++){
     L2 += std::pow(x,n)/(n*n);
   }
 
   G4double r_c;
 
   r_c = 2.*L2-(pi*pi/3.)-2.;
   r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x));
   r_c = r_c - std::log(x)*(2.*std::log(x)-1.);
   r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x);
 
   return r_c;
 }

Here is the caller graph for this function:

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

Public Member Functions

Protected Member Functions

Private Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ G4MuonDecayChannelWithSpin() [1/3]

◆ ~G4MuonDecayChannelWithSpin()

◆ G4MuonDecayChannelWithSpin() [2/3]

◆ G4MuonDecayChannelWithSpin() [3/3]

Member Function Documentation

◆ DecayIt()

◆ F_c()

◆ F_theta()

◆ operator=()

◆ R_c()