00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // 00028 00029 #include "G4PionDecayMakeSpin.hh" 00030 00031 #include "G4Decay.hh" 00032 #include "G4DecayProducts.hh" 00033 00034 #include "G4RandomDirection.hh" 00035 00036 // constructor 00037 00038 G4PionDecayMakeSpin::G4PionDecayMakeSpin(const G4String& processName) 00039 : G4Decay(processName) 00040 { 00041 // set Process Sub Type 00042 SetProcessSubType(static_cast<int>(DECAY_PionMakeSpin)); 00043 00044 } 00045 00046 G4PionDecayMakeSpin::~G4PionDecayMakeSpin() { } 00047 00048 void G4PionDecayMakeSpin::DaughterPolarization(const G4Track& aTrack, 00049 G4DecayProducts* products) 00050 { 00051 // This routine deals only with particles that can decay into a muon 00052 // pi+, pi-, K+, K- and K0_long 00053 00054 // get particle 00055 00056 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 00057 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition(); 00058 00059 G4ParticleDefinition* aMuonPlus = 00060 G4ParticleTable::GetParticleTable()->FindParticle("mu+"); 00061 G4ParticleDefinition* aMuonMinus = 00062 G4ParticleTable::GetParticleTable()->FindParticle("mu-"); 00063 G4ParticleDefinition* aPionPlus = 00064 G4ParticleTable::GetParticleTable()->FindParticle("pi+"); 00065 G4ParticleDefinition* aPionMinus = 00066 G4ParticleTable::GetParticleTable()->FindParticle("pi-"); 00067 G4ParticleDefinition* aKaonPlus = 00068 G4ParticleTable::GetParticleTable()->FindParticle("kaon+"); 00069 G4ParticleDefinition* aKaonMinus = 00070 G4ParticleTable::GetParticleTable()->FindParticle("kaon-"); 00071 G4ParticleDefinition* aKaon0Long = 00072 G4ParticleTable::GetParticleTable()->FindParticle("kaon0L"); 00073 G4ParticleDefinition* aNeutrinoMu = 00074 G4ParticleTable::GetParticleTable()->FindParticle("nu_mu"); 00075 G4ParticleDefinition* aAntiNeutrinoMu = 00076 G4ParticleTable::GetParticleTable()->FindParticle("anti_nu_mu"); 00077 00078 if( aParticleDef == aPionPlus || 00079 aParticleDef == aPionMinus || 00080 aParticleDef == aKaonPlus || 00081 aParticleDef == aKaonMinus || 00082 aParticleDef == aKaon0Long ) { 00083 } else { 00084 return; 00085 } 00086 00087 G4DynamicParticle* aMuon = NULL; 00088 00089 G4double emu(0), eneutrino(0); 00090 G4ThreeVector p_muon, p_neutrino; 00091 00092 G4int numberOfSecondaries = products->entries(); 00093 00094 if (numberOfSecondaries > 0) { 00095 for (G4int index=0; index < numberOfSecondaries; index++) 00096 { 00097 G4DynamicParticle* aSecondary = (*products)[index]; 00098 const G4ParticleDefinition* aSecondaryDef = aSecondary->GetDefinition(); 00099 00100 if (aSecondaryDef == aMuonPlus || 00101 aSecondaryDef == aMuonMinus ) { 00102 // Muon+ or Muon- 00103 aMuon = aSecondary; 00104 emu = aSecondary->GetTotalEnergy(); 00105 p_muon = aSecondary->GetMomentum(); 00106 } else if (aSecondaryDef == aNeutrinoMu || 00107 aSecondaryDef == aAntiNeutrinoMu ) { 00108 // Muon-Neutrino / Muon-Anti-Neutrino 00109 eneutrino = aSecondary->GetTotalEnergy(); 00110 p_neutrino = aSecondary->GetMomentum(); 00111 } 00112 } 00113 } 00114 00115 // This routine deals only with decays with a 00116 // muon and mu-(anti)neutrinos in the final state 00117 00118 if (!aMuon) return; 00119 if (eneutrino==0||emu==0) return; 00120 00121 G4ThreeVector spin(0,0,0); 00122 00123 const G4DynamicParticle* theParentParticle = products->GetParentParticle(); 00124 00125 G4double amass = theParentParticle->GetMass(); 00126 G4double emmu = aMuonPlus->GetPDGMass(); 00127 00128 if (numberOfSecondaries == 2 ) { 00129 00130 G4double scale = - (eneutrino - ( p_muon * p_neutrino )/(emu+emmu)); 00131 00132 p_muon = scale * p_muon; 00133 p_neutrino = emmu * p_neutrino; 00134 spin = p_muon + p_neutrino; 00135 00136 scale = 2./(amass*amass-emmu*emmu); 00137 spin = scale * spin; 00138 00139 if (aParticle->GetCharge() < 0.0) spin = -spin; 00140 00141 } else { 00142 00143 spin = G4RandomDirection(); 00144 00145 } 00146 00147 spin = spin.unit(); 00148 00149 aMuon->SetPolarization(spin.x(),spin.y(),spin.z()); 00150 00151 return; 00152 }