| 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 = nullptr;
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 G4DynamicParticle* aSecondary = (*products)[index];
97 const G4ParticleDefinition* aSecondaryDe 97 const G4ParticleDefinition* aSecondaryDef = aSecondary->GetDefinition();
98 98
99 if (aSecondaryDef == aMuonPlus || 99 if (aSecondaryDef == aMuonPlus ||
100 aSecondaryDef == aMuonMinus ) { 100 aSecondaryDef == aMuonMinus ) {
101 // Muon+ or Muon- 101 // Muon+ or Muon-
102 aMuon = aSecondary; 102 aMuon = aSecondary;
103 emu = aSecondary->GetTotalEnergy(); 103 emu = aSecondary->GetTotalEnergy();
104 p_muon = aSecondary->GetMomentum(); 104 p_muon = aSecondary->GetMomentum();
105 } else if (aSecondaryDef == aNeutrinoMu 105 } else if (aSecondaryDef == aNeutrinoMu ||
106 aSecondaryDef == aAntiNeutrinoMu ) { 106 aSecondaryDef == aAntiNeutrinoMu ) {
107 // Muon-Neutrino / Muon-Anti-Neutri 107 // Muon-Neutrino / Muon-Anti-Neutrino
108 eneutrino = aSecondary->GetTotalEnergy(); 108 eneutrino = aSecondary->GetTotalEnergy();
109 p_neutrino = aSecondary->GetMomentum(); 109 p_neutrino = aSecondary->GetMomentum();
110 } 110 }
111 } 111 }
112 } 112 }
113 113
114 // This routine deals only with decays with 114 // This routine deals only with decays with a
115 // muon and mu-(anti)neutrinos in the final 115 // muon and mu-(anti)neutrinos in the final state
116 if (aMuon == nullptr) return; 116 if (aMuon == nullptr) return;
117 if (eneutrino==0||emu==0) return; 117 if (eneutrino==0||emu==0) return;
118 118
119 G4ThreeVector spin(0,0,0); 119 G4ThreeVector spin(0,0,0);
120 120
121 const G4DynamicParticle* theParentParticle = 121 const G4DynamicParticle* theParentParticle = products->GetParentParticle();
122 122
123 G4double amass = theParentParticle->GetMass( 123 G4double amass = theParentParticle->GetMass();
124 G4double emmu = aMuonPlus->GetPDGMass(); 124 G4double emmu = aMuonPlus->GetPDGMass();
125 125
126 if (numberOfSecondaries == 2 ) { 126 if (numberOfSecondaries == 2 ) {
127 G4double scale = - (eneutrino - ( p_muon 127 G4double scale = - (eneutrino - ( p_muon * p_neutrino )/(emu+emmu));
128 128
129 p_muon = scale * p_muon; 129 p_muon = scale * p_muon;
130 p_neutrino = emmu * p_neutrino; 130 p_neutrino = emmu * p_neutrino;
131 spin = p_muon + p_neutrino; 131 spin = p_muon + p_neutrino;
132 132
133 scale = 2./(amass*amass-emmu*emmu); 133 scale = 2./(amass*amass-emmu*emmu);
134 spin = scale * spin; 134 spin = scale * spin;
135 135
136 if (aParticle->GetCharge() < 0.0) spin = 136 if (aParticle->GetCharge() < 0.0) spin = -spin;
137 137
138 } else { 138 } else {
139 spin = G4RandomDirection(); 139 spin = G4RandomDirection();
140 140
141 } 141 }
142 142
143 spin = spin.unit(); 143 spin = spin.unit();
144 144
145 aMuon->SetPolarization(spin.x(),spin.y(),spi 145 aMuon->SetPolarization(spin.x(),spin.y(),spin.z());
146 146
147 return; 147 return;
148 } 148 }
149 149
150 void G4PionDecayMakeSpin::ProcessDescription(s 150 void G4PionDecayMakeSpin::ProcessDescription(std::ostream& outFile) const
151 { 151 {
152 outFile << GetProcessName() 152 outFile << GetProcessName()
153 << ": Decay of mesons that can decay into 153 << ": Decay of mesons that can decay into a muon \n"
154 << " i.e. pi+, pi-, K+, K- and K0_long \n" 154 << " i.e. pi+, pi-, K+, K- and K0_long \n"
155 << " kinematics of daughters are dertermin 155 << " kinematics of daughters are dertermined by DecayChannels \n"
156 << " polarization of daughter particles ar 156 << " polarization of daughter particles are take into account. \n";
157 } 157 }
158 158
159 159