G4MuonDecayChannelWithSpin.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuonDecayChannelWithSpin.cc

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
// ------------------------------------------------------------
//      GEANT 4 class header file
//
//      History:
//               17 August 2004 P.Gumplinger and T.MacPhail
//               samples Michel spectrum including 1st order
//               radiative corrections
//               Reference: Florian Scheck "Muon Physics", in Physics Reports
//                          (Review Section of Physics Letters) 44, No. 4 (1978)
//                          187-248. North-Holland Publishing Company, Amsterdam
//                          at page 210 cc.
//
//                          W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
//
// ------------------------------------------------------------
//
#include "G4MuonDecayChannelWithSpin.hh"

#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"

#include "G4DecayProducts.hh"
#include "G4LorentzVector.hh"

G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin()
  : G4MuonDecayChannel()
{
}

G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4String& theParentName, 
                   G4double        theBR)
  : G4MuonDecayChannel(theParentName,theBR)
{
}

G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin()
{
}

G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4MuonDecayChannelWithSpin &right):
  G4MuonDecayChannel(right)
{
}

G4MuonDecayChannelWithSpin & G4MuonDecayChannelWithSpin::operator=(const G4MuonDecayChannelWithSpin & 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 !=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;
}


G4DecayProducts *G4MuonDecayChannelWithSpin::DecayIt(G4double) 
{
  // 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-energy; 
  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;
    G4double TT = daughterparticle0->GetTotalEnergy()
      + daughterparticle1->GetTotalEnergy()
      + daughterparticle2->GetTotalEnergy(); 
    G4cout << "e  " << daughterparticle0->GetTotalEnergy()/MeV << G4endl; 
    G4cout << "nu1" << daughterparticle1->GetTotalEnergy()/MeV << G4endl; 
    G4cout << "nu2" << daughterparticle2->GetTotalEnergy()/MeV << G4endl; 
    G4cout << "total" << (TT-parentmass)/keV << G4endl;
    if (GetVerboseLevel()>2) { products->DumpInfo(); }
  }
#endif

  return products;
}

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

  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;
}