diff --git a/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc index bdb9345..41bc1e3 100644 --- a/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc +++ b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc @@ -1,190 +1,209 @@ -// -// ******************************************************************** -// * DISCLAIMER * -// * * -// * The following disclaimer summarizes all the specific disclaimers * -// * of contributors to this software. The specific disclaimers,which * -// * govern, are listed with their locations in: * -// * http://cern.ch/geant4/license * -// * * -// * 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. * -// * * -// * This code implementation is the intellectual property of the * -// * GEANT4 collaboration. * -// * By copying, distributing or modifying the Program (or any work * -// * based on the Program) you indicate your acceptance of this * -// * statement, and all its terms. * -// ******************************************************************** -// -// -// $Id: G4MuonDecayChannel.cc,v 1.11 2004/12/10 18:02:04 gcosmo Exp $ -// GEANT4 tag $Name: geant4-07-00-cand-05 $ -// -// -// ------------------------------------------------------------ -// GEANT 4 class header file -// -// History: first implementation, based on object model of -// 30 May 1997 H.Kurashige -// -// Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige -// ------------------------------------------------------------ - -#include "G4ParticleDefinition.hh" -#include "G4DecayProducts.hh" -#include "G4VDecayChannel.hh" -#include "G4MuonDecayChannel.hh" -#include "Randomize.hh" -#include "G4LorentzVector.hh" -#include "G4LorentzRotation.hh" - - -G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName, - G4double theBR) - :G4VDecayChannel("Muon Decay",1) -{ - // set names for daughter particles - if (theParentName == "mu+") { - SetBR(theBR); - SetParent("mu+"); - SetNumberOfDaughters(3); - SetDaughter(0, "e+"); - SetDaughter(1, "nu_e"); - SetDaughter(2, "anti_nu_mu"); - } else if (theParentName == "Mu") { - SetBR(theBR); - SetParent("Mu"); - SetNumberOfDaughters(3); - SetDaughter(0, "e+"); - SetDaughter(1, "nu_e"); - SetDaughter(2, "anti_nu_mu"); - } else if (theParentName == "mu-") { - SetBR(theBR); - SetParent("mu-"); - SetNumberOfDaughters(3); - SetDaughter(0, "e-"); - SetDaughter(1, "anti_nu_e"); - SetDaughter(2, "nu_mu"); - } else { -#ifdef G4VERBOSE - if (GetVerboseLevel()>0) { - G4cout << "G4MuonDecayChannel:: constructor :"; - G4cout << " parent particle is not muon but "; - G4cout << theParentName << G4endl; - } -#endif - } -} - -G4MuonDecayChannel::~G4MuonDecayChannel() -{ -} - -G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double) -{ - // this version neglects muon polarization - // assumes the pure V-A coupling - // gives incorrect energy spectrum for neutrinos -#ifdef G4VERBOSE - if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt "; -#endif - - if (parent == 0) FillParent(); - if (daughters == 0) FillDaughters(); - - // parent mass - G4double parentmass = parent->GetPDGMass(); - - //daughters'mass - G4double daughtermass[3]; - G4double sumofdaughtermass = 0.0; - for (G4int index=0; index<3; index++){ - daughtermass[index] = daughters[index]->GetPDGMass(); - sumofdaughtermass += daughtermass[index]; - } - - //create parent G4DynamicParticle at rest - G4ThreeVector dummy; - G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); - //create G4Decayproducts - G4DecayProducts *products = new G4DecayProducts(*parentparticle); - delete parentparticle; - - // calculate daughter momentum - G4double daughtermomentum[3]; - G4double energy; - // calcurate electron energy - G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass); - G4double x; - G4double r; - do { - do { - r = G4UniformRand(); - x = xmax*G4UniformRand(); - } while (r > (3.0 - 2.0*x)*x*x); - energy = x*parentmass/2.0 - daughtermass[0]; - } while (energy <0.0); - //create daughter G4DynamicParticle - // daughter 0 (electron) - daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]); - G4double costheta, sintheta, phi, sinphi, cosphi; - costheta = 2.*G4UniformRand()-1.0; - sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); - phi = twopi*G4UniformRand()*rad; - sinphi = std::sin(phi); - cosphi = std::cos(phi); - G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta); - G4DynamicParticle * daughterparticle - = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]); - products->PushProducts(daughterparticle); - - // daughter 1 ,2 (nutrinos) - // 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( daughters[1], direction1*(vmass/2.)); - G4DynamicParticle * daughterparticle2 - = new G4DynamicParticle( 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 << "G4MuonDecayChannel::DecayIt "; - G4cout << " create decay products in rest frame " <DumpInfo(); - } -#endif - return products; -} - - - - - - +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * 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. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4MuonDecayChannel.cc,v 1.13 2005/06/23 11:02:26 gcosmo Exp $ +// GEANT4 tag $Name: geant4-07-01 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 30 May 1997 H.Kurashige +// +// Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige +//2005 +// M. Melissas ( melissas AT cppm.in2p3.fr) +// J. Brunner ( brunner AT cppm.in2p3.fr) +// Adding V-A fluxes for neutrinos using a new algortithm : +// ------------------------------------------------------------ + +#include "G4ParticleDefinition.hh" +#include "G4DecayProducts.hh" +#include "G4VDecayChannel.hh" +#include "G4MuonDecayChannel.hh" +#include "Randomize.hh" +#include "G4LorentzVector.hh" +#include "G4LorentzRotation.hh" +#include "G4RotationMatrix.hh" + + +G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName, + G4double theBR) + :G4VDecayChannel("Muon Decay",1) +{ + // set names for daughter particles + if (theParentName == "mu+") { + SetBR(theBR); + SetParent("mu+"); + SetNumberOfDaughters(3); + SetDaughter(0, "e+"); + SetDaughter(1, "nu_e"); + SetDaughter(2, "anti_nu_mu"); + } else if (theParentName == "Mu") { + SetBR(theBR); + SetParent("Mu"); + SetNumberOfDaughters(3); + SetDaughter(0, "e+"); + SetDaughter(1, "nu_e"); + SetDaughter(2, "anti_nu_mu"); + } else if (theParentName == "mu-") { + SetBR(theBR); + SetParent("mu-"); + SetNumberOfDaughters(3); + SetDaughter(0, "e-"); + SetDaughter(1, "anti_nu_e"); + SetDaughter(2, "nu_mu"); + } else { +#ifdef G4VERBOSE + if (GetVerboseLevel()>0) { + G4cout << "G4MuonDecayChannel:: constructor :"; + G4cout << " parent particle is not muon but "; + G4cout << theParentName << G4endl; + } +#endif + } +} + +G4MuonDecayChannel::~G4MuonDecayChannel() +{ +} + +G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double) +{ + // this version neglects muon polarization,and electron mass + // assumes the pure V-A coupling + // the Neutrinos are correctly V-A. +#ifdef G4VERBOSE + if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt "; +#endif + + if (parent == 0) FillParent(); + if (daughters == 0) FillDaughters(); + + // parent mass + G4double parentmass = parent->GetPDGMass(); + + //daughters'mass + G4double daughtermass[3]; + G4double sumofdaughtermass = 0.0; + for (G4int index=0; index<3; index++){ + daughtermass[index] = daughters[index]->GetPDGMass(); + sumofdaughtermass += daughtermass[index]; + } + + //create parent G4DynamicParticle at rest + G4ThreeVector dummy; + G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); + //create G4Decayproducts + G4DecayProducts *products = new G4DecayProducts(*parentparticle); + delete parentparticle; + + // calculate daughter momentum + G4double daughtermomentum[3]; + // calcurate electron energy + G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass); + G4double x; + + G4double Ee,Ene; + + G4double gam; + G4double EMax=parentmass/2-daughtermass[0]; + + + //Generating Random Energy +do { + Ee=G4UniformRand(); + do{ + x=xmax*G4UniformRand(); + gam=G4UniformRand(); + }while (gam >x*(1.-x)); + Ene=x; + } while ( Ene < (1.-Ee)); + G4double Enm=(2.-Ee-Ene); + + + //initialisation of rotation parameters + + G4double costheta,sintheta,rphi,rtheta,rpsi; + costheta= 1.-2./Ee-2./Ene+2./Ene/Ee; + sintheta=sqrt(1.-costheta*costheta); + + + rphi=twopi*G4UniformRand()*rad; + rtheta=(acos(2.*G4UniformRand()-1.)); + rpsi=twopi*G4UniformRand()*rad; + + G4RotationMatrix *rot= new G4RotationMatrix(); + rot->set(rphi,rtheta,rpsi); + + //electron 0 + daughtermomentum[0]=sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]); + G4ThreeVector *direction0 =new G4ThreeVector(0.0,0.0,1.0); + + *direction0 *= *rot; + + G4DynamicParticle * daughterparticle = new G4DynamicParticle ( daughters[0], *direction0 * daughtermomentum[0]); + + products->PushProducts(daughterparticle); + + //electronic neutrino 1 + + daughtermomentum[1]=sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]); + G4ThreeVector *direction1 =new G4ThreeVector(sintheta,0.0,costheta); + + *direction1 *= *rot; + + G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( daughters[1], *direction1 * daughtermomentum[1]); + products->PushProducts(daughterparticle1); + + //muonnic neutrino 2 + + daughtermomentum[2]=sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]); + G4ThreeVector *direction2 =new G4ThreeVector(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta); + + *direction2 *= *rot; + + G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( daughters[2], + *direction2 * daughtermomentum[2]); + products->PushProducts(daughterparticle2); + + + + + // output message +#ifdef G4VERBOSE + if (GetVerboseLevel()>1) { + G4cout << "G4MuonDecayChannel::DecayIt "; + G4cout << " create decay products in rest frame " <DumpInfo(); + } +#endif + return products; +} + + + + + + diff --git a/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh index 0a3f6d7..a58dfa0 100644 --- a/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh +++ b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh @@ -22,7 +22,7 @@ // // // $Id: G4MuonDecayChannel.hh,v 1.5 2001/07/11 10:01:56 gunter Exp $ -// GEANT4 tag $Name: geant4-07-00-cand-01 $ +// GEANT4 tag $Name: geant4-07-01 $ // // // ------------------------------------------------------------