G4MuonRadiativeDecayChannelWithSpin.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// G4MuonRadiativeDecayChannelWithSpin class implementation
//
// References:
// - TRIUMF/TWIST Technote TN-55:
//   "Radiative muon decay" by P. Depommier and A. Vacheret
// - Yoshitaka Kuno and Yasuhiro Okada
//   "Muon Decays and Physics Beyond the Standard Model"
//   Rev. Mod. Phys. 73, 151 (2001)
//
// Author: P.Gumplinger - Triumf, 25 July 2007
// Revision: D.Mingming - Center for HEP, Tsinghua Univ., 10 August 2011
// --------------------------------------------------------------------
 
#include "G4MuonRadiativeDecayChannelWithSpin.hh"
 
#include "G4DecayProducts.hh"
#include "G4LorentzVector.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
 
G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin(
  const G4String& theParentName, G4double theBR)
  : 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
  }
}
 
G4MuonRadiativeDecayChannelWithSpin&
G4MuonRadiativeDecayChannelWithSpin::operator=(const G4MuonRadiativeDecayChannelWithSpin& right)
{
  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;
}
 
G4DecayProducts* G4MuonRadiativeDecayChannelWithSpin::DecayIt(G4double)
{
#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;
}
 
G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu, G4double x, G4double y,
                                                   G4double cthetaE, G4double cthetaG,
                                                   G4double cthetaGE)
{
  G4double mu = 105.65;
  G4double me = 0.511;
  G4double rho = 0.75;
  G4double del = 0.75;
  G4double eps = 0.0;
  G4double kap = 0.0;
  G4double ksi = 1.0;
 
  G4double delta = 1 - cthetaGE;
 
  // Calculation of the functions f(x,y)
 
  G4double f_1s = 12.0
                  * ((y * y) * (1.0 - y) + x * y * (2.0 - 3.0 * y) + 2.0 * (x * x) * (1.0 - 2.0 * y)
                     - 2.0 * (x * x * x));
  G4double f0s = 6.0
                 * (-x * y * (2.0 - 3.0 * (y * y)) - 2.0 * (x * x) * (1.0 - y - 3.0 * (y * y))
                    + 2.0 * (x * x * x) * (1.0 + 2.0 * y));
  G4double f1s =
    3.0 * ((x * x) * y * (2.0 - 3.0 * y - 3.0 * (y * y)) - (x * x * x) * y * (4.0 + 3.0 * y));
  G4double f2s = 1.5 * ((x * x * x) * (y * y) * (2.0 + y));
 
  G4double f_1se = 12.0 * (x * y * (1.0 - y) + (x * x) * (2.0 - 3.0 * y) - 2.0 * (x * x * x));
  G4double f0se = 6.0 * (-(x * x) * (2.0 - y - 2.0 * (y * y)) + (x * x * x) * (2.0 + 3.0 * y));
  G4double f1se = -3.0 * (x * x * x) * y * (2.0 + y);
  G4double f2se = 0.0;
 
  G4double f_1sg = 12.0 * ((y * y) * (1.0 - y) + x * y * (1.0 - 2.0 * y) - (x * x) * y);
  G4double f0sg =
    6.0 * (-x * (y * y) * (2.0 - 3.0 * y) - (x * x) * y * (1.0 - 4.0 * y) + (x * x * x) * y);
  G4double f1sg = 3.0 * ((x * x) * (y * y) * (1.0 - 3.0 * y) - 2.0 * (x * x * x) * (y * y));
  G4double f2sg = 1.5 * (x * x * x) * (y * y * y);
 
  G4double f_1v = 8.0
                  * ((y * y) * (3.0 - 2.0 * y) + 6.0 * x * y * (1.0 - y)
                     + 2.0 * (x * x) * (3.0 - 4.0 * y) - 4.0 * (x * x * x));
  G4double f0v = 8.0
                 * (-x * y * (3.0 - y - (y * y)) - (x * x) * (3.0 - y - 4.0 * (y * y))
                    + 2.0 * (x * x * x) * (1.0 + 2.0 * y));
  G4double f1v =
    2.0 * ((x * x) * y * (6.0 - 5.0 * y - 2.0 * (y * y)) - 2.0 * (x * x * x) * y * (4.0 + 3.0 * y));
  G4double f2v = 2.0 * (x * x * x) * (y * y) * (2.0 + y);
 
  G4double f_1ve =
    8.0 * (x * y * (1.0 - 2.0 * y) + 2.0 * (x * x) * (1.0 - 3.0 * y) - 4.0 * (x * x * x));
  G4double f0ve =
    4.0 * (-(x * x) * (2.0 - 3.0 * y - 4.0 * (y * y)) + 2.0 * (x * x * x) * (2.0 + 3.0 * y));
  G4double f1ve = -4.0 * (x * x * x) * y * (2.0 + y);
  G4double f2ve = 0.0;
 
  G4double f_1vg = 8.0 * ((y * y) * (1.0 - 2.0 * y) + x * y * (1.0 - 4.0 * y) - 2.0 * (x * x) * y);
  G4double f0vg =
    4.0 * (2.0 * x * (y * y) * (1.0 + y) - (x * x) * y * (1.0 - 4.0 * y) + 2.0 * (x * x * x) * y);
  G4double f1vg = 2.0 * ((x * x) * (y * y) * (1.0 - 2.0 * y) - 4.0 * (x * x * x) * (y * y));
  G4double f2vg = 2.0 * (x * x * x) * (y * y * y);
 
  G4double f_1t = 8.0
                  * ((y * y) * (3.0 - y) + 3.0 * x * y * (2.0 - y) + 2.0 * (x * x) * (3.0 - 2.0 * y)
                     - 2.0 * (x * x * x));
  G4double f0t = 4.0
                 * (-x * y * (6.0 + (y * y)) - 2.0 * (x * x) * (3.0 + y - 3.0 * (y * y))
                    + 2.0 * (x * x * x) * (1.0 + 2.0 * y));
  G4double f1t =
    2.0 * ((x * x) * y * (6.0 - 5.0 * y + (y * y)) - (x * x * x) * y * (4.0 + 3.0 * y));
  G4double f2t = (x * x * x) * (y * y) * (2.0 + y);
 
  G4double f_1te = -8.0 * (x * y * (1.0 + 3.0 * y) + (x * x) * (2.0 + 3.0 * y) + 2.0 * (x * x * x));
  G4double f0te = 4.0 * ((x * x) * (2.0 + 3.0 * y + 4.0 * (y * y)) + (x * x * x) * (2.0 + 3.0 * y));
  G4double f1te = -2.0 * (x * x * x) * y * (2.0 + y);
  G4double f2te = 0.0;
 
  G4double f_1tg = -8.0 * ((y * y) * (1.0 + y) + x * y + (x * x) * y);
  G4double f0tg = 4.0 * (x * (y * y) * (2.0 - y) + (x * x) * y * (1.0 + 2.0 * y) + (x * x * x) * y);
  G4double f1tg = -2.0 * ((x * x) * (y * y) * (1.0 - y) + 2.0 * (x * x * x) * y);
  G4double f2tg = (x * x * x) * (y * y * y);
 
  G4double term = delta + 2.0 * (me * me) / ((mu * mu) * (x * x));
  term = 1.0 / term;
 
  G4double nss = term * f_1s + f0s + delta * f1s + (delta * delta) * f2s;
  G4double nv = term * f_1v + f0v + delta * f1v + (delta * delta) * f2v;
  G4double nt = term * f_1t + f0t + delta * f1t + (delta * delta) * f2t;
 
  G4double nse = term * f_1se + f0se + delta * f1se + (delta * delta) * f2se;
  G4double nve = term * f_1ve + f0ve + delta * f1ve + (delta * delta) * f2ve;
  G4double nte = term * f_1te + f0te + delta * f1te + (delta * delta) * f2te;
 
  G4double nsg = term * f_1sg + f0sg + delta * f1sg + (delta * delta) * f2sg;
  G4double nvg = term * f_1vg + f0vg + delta * f1vg + (delta * delta) * f2vg;
  G4double ntg = term * f_1tg + f0tg + delta * f1tg + (delta * delta) * f2tg;
 
  G4double term1 = nv;
  G4double term2 = 2.0 * nss + nv - nt;
  G4double term3 = 2.0 * nss - 2.0 * nv + nt;
 
  G4double term1e = 1.0 / 3.0 * (1.0 - 4.0 / 3.0 * del);
  G4double term2e = 2.0 * nse + 5.0 * nve - nte;
  G4double term3e = 2.0 * nse - 2.0 * nve + nte;
 
  G4double term1g = 1.0 / 3.0 * (1.0 - 4.0 / 3.0 * del);
  G4double term2g = 2.0 * nsg + 5.0 * nvg - ntg;
  G4double term3g = 2.0 * nsg - 2.0 * nvg + ntg;
 
  G4double som00 = term1 + (1.0 - 4.0 / 3.0 * rho) * term2 + eps * term3;
  G4double som01 = Pmu * ksi
                   * (cthetaE * (nve - term1e * term2e + kap * term3e)
                      + cthetaG * (nvg - term1g * term2g + kap * term3g));
 
  G4double som0 = (som00 + som01) / y;
  som0 = fine_structure_const / 8. / (twopi * twopi * twopi) * som0;
 
  return som0;
}