Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // 27 // 28 28 29 #include "G4PionDecayMakeSpin.hh" 29 #include "G4PionDecayMakeSpin.hh" 30 30 31 #include "G4Decay.hh" 31 #include "G4Decay.hh" 32 #include "G4DecayProducts.hh" 32 #include "G4DecayProducts.hh" 33 33 34 #include "G4RandomDirection.hh" 34 #include "G4RandomDirection.hh" 35 35 36 // constructor 36 // constructor 37 37 38 G4PionDecayMakeSpin::G4PionDecayMakeSpin(const 38 G4PionDecayMakeSpin::G4PionDecayMakeSpin(const G4String& processName) 39 : G4Decay(proce 39 : G4Decay(processName) 40 { 40 { 41 // set Process Sub Type 41 // set Process Sub Type 42 SetProcessSubType(static_cast<int>(DECAY_Pio 42 SetProcessSubType(static_cast<int>(DECAY_PionMakeSpin)); 43 43 44 } 44 } 45 45 46 G4PionDecayMakeSpin::~G4PionDecayMakeSpin() { 46 G4PionDecayMakeSpin::~G4PionDecayMakeSpin() { } 47 47 48 void G4PionDecayMakeSpin::DaughterPolarization 48 void G4PionDecayMakeSpin::DaughterPolarization(const G4Track& aTrack, 49 G4DecayProd 49 G4DecayProducts* products) 50 { 50 { 51 // This routine deals only with particles t 51 // This routine deals only with particles that can decay into a muon 52 // pi+, pi-, K+, K- and K0_l 52 // pi+, pi-, K+, K- and K0_long 53 53 54 // get particle 54 // get particle 55 55 56 const G4DynamicParticle* aParticle = aTrack. 56 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 57 const G4ParticleDefinition* aParticleDef = a 57 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition(); 58 58 59 G4ParticleDefinition* aMuonPlus = 59 G4ParticleDefinition* aMuonPlus = 60 G4ParticleTable::GetParticleTable( 60 G4ParticleTable::GetParticleTable()->FindParticle("mu+"); 61 G4ParticleDefinition* aMuonMinus = 61 G4ParticleDefinition* aMuonMinus = 62 G4ParticleTable::GetParticleTable( 62 G4ParticleTable::GetParticleTable()->FindParticle("mu-"); 63 G4ParticleDefinition* aPionPlus = 63 G4ParticleDefinition* aPionPlus = 64 G4ParticleTable::GetParticleTable( 64 G4ParticleTable::GetParticleTable()->FindParticle("pi+"); 65 G4ParticleDefinition* aPionMinus = 65 G4ParticleDefinition* aPionMinus = 66 G4ParticleTable::GetParticleTable( 66 G4ParticleTable::GetParticleTable()->FindParticle("pi-"); 67 G4ParticleDefinition* aKaonPlus = 67 G4ParticleDefinition* aKaonPlus = 68 G4ParticleTable::GetParticleTable( 68 G4ParticleTable::GetParticleTable()->FindParticle("kaon+"); 69 G4ParticleDefinition* aKaonMinus = 69 G4ParticleDefinition* aKaonMinus = 70 G4ParticleTable::GetParticleTable( 70 G4ParticleTable::GetParticleTable()->FindParticle("kaon-"); 71 G4ParticleDefinition* aKaon0Long = 71 G4ParticleDefinition* aKaon0Long = 72 G4ParticleTable::GetParticleTable( 72 G4ParticleTable::GetParticleTable()->FindParticle("kaon0L"); 73 G4ParticleDefinition* aNeutrinoMu = 73 G4ParticleDefinition* aNeutrinoMu = 74 G4ParticleTable::GetParticleTable( 74 G4ParticleTable::GetParticleTable()->FindParticle("nu_mu"); 75 G4ParticleDefinition* aAntiNeutrinoMu = 75 G4ParticleDefinition* aAntiNeutrinoMu = 76 G4ParticleTable::GetParticleTable( 76 G4ParticleTable::GetParticleTable()->FindParticle("anti_nu_mu"); 77 77 78 if( aParticleDef == aPionPlus || 78 if( aParticleDef == aPionPlus || 79 aParticleDef == aPionMinus || 79 aParticleDef == aPionMinus || 80 aParticleDef == aKaonPlus || 80 aParticleDef == aKaonPlus || 81 aParticleDef == aKaonMinus || 81 aParticleDef == aKaonMinus || 82 aParticleDef == aKaon0Long ) { 82 aParticleDef == aKaon0Long ) { 83 } else { 83 } else { 84 return; 84 return; 85 } 85 } 86 86 87 G4DynamicParticle* aMuon = nullptr; << 87 G4DynamicParticle* aMuon = NULL; 88 88 89 G4double emu(0), eneutrino(0); 89 G4double emu(0), eneutrino(0); 90 G4ThreeVector p_muon, p_neutrino; 90 G4ThreeVector p_muon, p_neutrino; 91 91 92 G4int numberOfSecondaries = products->entrie 92 G4int numberOfSecondaries = products->entries(); 93 93 94 if (numberOfSecondaries > 0) { 94 if (numberOfSecondaries > 0) { 95 for (G4int index=0; index < numberOfSecond << 95 for (G4int index=0; index < numberOfSecondaries; index++) 96 G4DynamicParticle* aSecondary = (*produc << 96 { 97 const G4ParticleDefinition* aSecondaryDe << 97 G4DynamicParticle* aSecondary = (*products)[index]; 98 << 98 const G4ParticleDefinition* aSecondaryDef = aSecondary->GetDefinition(); 99 if (aSecondaryDef == aMuonPlus || << 99 100 aSecondaryDef == aMuonMinus ) { << 100 if (aSecondaryDef == aMuonPlus || 101 // Muon+ or Muon- << 101 aSecondaryDef == aMuonMinus ) { 102 aMuon = aSecondary; << 102 // Muon+ or Muon- 103 emu = aSecondary->GetTotalEnergy(); << 103 aMuon = aSecondary; 104 p_muon = aSecondary->GetMomentum(); << 104 emu = aSecondary->GetTotalEnergy(); 105 } else if (aSecondaryDef == aNeutrinoMu << 105 p_muon = aSecondary->GetMomentum(); 106 aSecondaryDef == aAntiNeutrinoMu ) { << 106 } else if (aSecondaryDef == aNeutrinoMu || 107 // Muon-Neutrino / Muon-Anti-Neutri << 107 aSecondaryDef == aAntiNeutrinoMu ) { 108 eneutrino = aSecondary->GetTotalEnergy(); << 108 // Muon-Neutrino / Muon-Anti-Neutrino 109 p_neutrino = aSecondary->GetMomentum(); << 109 eneutrino = aSecondary->GetTotalEnergy(); 110 } << 110 p_neutrino = aSecondary->GetMomentum(); >> 111 } 111 } 112 } 112 } 113 } 113 114 114 // This routine deals only with decays with 115 // This routine deals only with decays with a 115 // muon and mu-(anti)neutrinos in the final 116 // muon and mu-(anti)neutrinos in the final state 116 if (aMuon == nullptr) return; << 117 >> 118 if (!aMuon) return; 117 if (eneutrino==0||emu==0) return; 119 if (eneutrino==0||emu==0) return; 118 120 119 G4ThreeVector spin(0,0,0); 121 G4ThreeVector spin(0,0,0); 120 122 121 const G4DynamicParticle* theParentParticle = 123 const G4DynamicParticle* theParentParticle = products->GetParentParticle(); 122 124 123 G4double amass = theParentParticle->GetMass( 125 G4double amass = theParentParticle->GetMass(); 124 G4double emmu = aMuonPlus->GetPDGMass(); 126 G4double emmu = aMuonPlus->GetPDGMass(); 125 127 126 if (numberOfSecondaries == 2 ) { 128 if (numberOfSecondaries == 2 ) { >> 129 127 G4double scale = - (eneutrino - ( p_muon 130 G4double scale = - (eneutrino - ( p_muon * p_neutrino )/(emu+emmu)); 128 131 129 p_muon = scale * p_muon; 132 p_muon = scale * p_muon; 130 p_neutrino = emmu * p_neutrino; 133 p_neutrino = emmu * p_neutrino; 131 spin = p_muon + p_neutrino; 134 spin = p_muon + p_neutrino; 132 135 133 scale = 2./(amass*amass-emmu*emmu); 136 scale = 2./(amass*amass-emmu*emmu); 134 spin = scale * spin; 137 spin = scale * spin; 135 138 136 if (aParticle->GetCharge() < 0.0) spin = 139 if (aParticle->GetCharge() < 0.0) spin = -spin; 137 140 138 } else { 141 } else { >> 142 139 spin = G4RandomDirection(); 143 spin = G4RandomDirection(); 140 144 141 } 145 } 142 146 143 spin = spin.unit(); 147 spin = spin.unit(); 144 148 145 aMuon->SetPolarization(spin.x(),spin.y(),spi 149 aMuon->SetPolarization(spin.x(),spin.y(),spin.z()); 146 150 147 return; 151 return; 148 } 152 } 149 << 150 void G4PionDecayMakeSpin::ProcessDescription(s << 151 { << 152 outFile << GetProcessName() << 153 << ": Decay of mesons that can decay into << 154 << " i.e. pi+, pi-, K+, K- and K0_long \n" << 155 << " kinematics of daughters are dertermin << 156 << " polarization of daughter particles ar << 157 } << 158 << 159 153