| 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 // G4MuonDecayChannel class implementation <<
27 // 26 //
28 // Author: H.Kurashige, 30 May 1997 << 27 // $Id: G4MuonDecayChannel.cc 95906 2016-03-02 10:56:50Z gcosmo $
29 // Contributions: << 28 //
30 // - 2005 - M.Melissas, J.Brunner CPPM/IN2P3 << 29 //
31 // Added V-A fluxes for neutrinos using a ne << 30 // ------------------------------------------------------------
32 // ------------------------------------------- << 31 // GEANT 4 class header file
33 << 32 //
34 #include "G4MuonDecayChannel.hh" << 33 // History: first implementation, based on object model of
>> 34 // 30 May 1997 H.Kurashige
>> 35 //
>> 36 // - Fix bug in calculation of electron energy in DecayIt()
>> 37 // 28 Feb 01 H.Kurashige
>> 38 // - Adding V-A fluxes for neutrinos using a new algorithm, 2005
>> 39 // M. Melissas ( melissas AT cppm.in2p3.fr)
>> 40 // J. Brunner ( brunner AT cppm.in2p3.fr)
>> 41 // ------------------------------------------------------------
35 42
36 #include "G4DecayProducts.hh" <<
37 #include "G4LorentzRotation.hh" <<
38 #include "G4LorentzVector.hh" <<
39 #include "G4ParticleDefinition.hh" 43 #include "G4ParticleDefinition.hh"
40 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh"
41 #include "G4RotationMatrix.hh" <<
42 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
>> 46 #include "G4DecayProducts.hh"
43 #include "G4VDecayChannel.hh" 47 #include "G4VDecayChannel.hh"
>> 48 #include "G4MuonDecayChannel.hh"
44 #include "Randomize.hh" 49 #include "Randomize.hh"
>> 50 #include "G4LorentzVector.hh"
>> 51 #include "G4LorentzRotation.hh"
>> 52 #include "G4RotationMatrix.hh"
>> 53
>> 54 G4MuonDecayChannel::G4MuonDecayChannel()
>> 55 :G4VDecayChannel()
>> 56 {
>> 57 }
45 58
46 G4MuonDecayChannel::G4MuonDecayChannel(const G << 59 G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName,
47 : G4VDecayChannel("Muon Decay", 1) << 60 G4double theBR)
>> 61 :G4VDecayChannel("Muon Decay",1)
48 { 62 {
49 // set names for daughter particles 63 // set names for daughter particles
50 if (theParentName == "mu+") { 64 if (theParentName == "mu+") {
51 SetBR(theBR); 65 SetBR(theBR);
52 SetParent("mu+"); 66 SetParent("mu+");
53 SetNumberOfDaughters(3); 67 SetNumberOfDaughters(3);
54 SetDaughter(0, "e+"); 68 SetDaughter(0, "e+");
55 SetDaughter(1, "nu_e"); 69 SetDaughter(1, "nu_e");
56 SetDaughter(2, "anti_nu_mu"); 70 SetDaughter(2, "anti_nu_mu");
57 } << 71 } else if (theParentName == "mu-") {
58 else if (theParentName == "mu-") { <<
59 SetBR(theBR); 72 SetBR(theBR);
60 SetParent("mu-"); 73 SetParent("mu-");
61 SetNumberOfDaughters(3); 74 SetNumberOfDaughters(3);
62 SetDaughter(0, "e-"); 75 SetDaughter(0, "e-");
63 SetDaughter(1, "anti_nu_e"); 76 SetDaughter(1, "anti_nu_e");
64 SetDaughter(2, "nu_mu"); 77 SetDaughter(2, "nu_mu");
65 } << 78 } else {
66 else { <<
67 #ifdef G4VERBOSE 79 #ifdef G4VERBOSE
68 if (GetVerboseLevel() > 0) { << 80 if (GetVerboseLevel()>0) {
69 G4cout << "G4MuonDecayChannel:: construc 81 G4cout << "G4MuonDecayChannel:: constructor :";
70 G4cout << " parent particle is not muon 82 G4cout << " parent particle is not muon but ";
71 G4cout << theParentName << G4endl; 83 G4cout << theParentName << G4endl;
72 } 84 }
73 #endif 85 #endif
74 } 86 }
75 } 87 }
76 88
77 G4MuonDecayChannel& G4MuonDecayChannel::operat << 89 G4MuonDecayChannel::G4MuonDecayChannel(const G4MuonDecayChannel &right):
>> 90 G4VDecayChannel(right)
>> 91 {
>> 92 }
>> 93
>> 94 G4MuonDecayChannel::~G4MuonDecayChannel()
>> 95 {
>> 96 }
>> 97
>> 98 G4MuonDecayChannel & G4MuonDecayChannel::operator=(const G4MuonDecayChannel & right)
78 { 99 {
79 if (this != &right) { << 100 if (this != &right) {
80 kinematics_name = right.kinematics_name; 101 kinematics_name = right.kinematics_name;
81 verboseLevel = right.verboseLevel; 102 verboseLevel = right.verboseLevel;
82 rbranch = right.rbranch; 103 rbranch = right.rbranch;
83 104
84 // copy parent name 105 // copy parent name
85 parent_name = new G4String(*right.parent_n 106 parent_name = new G4String(*right.parent_name);
86 107
87 // clear daughters_name array 108 // clear daughters_name array
88 ClearDaughtersName(); 109 ClearDaughtersName();
89 110
90 // recreate array 111 // recreate array
91 numberOfDaughters = right.numberOfDaughter 112 numberOfDaughters = right.numberOfDaughters;
92 if (numberOfDaughters > 0) { << 113 if ( numberOfDaughters >0 ) {
93 if (daughters_name != nullptr) ClearDaug << 114 if (daughters_name !=0) ClearDaughtersName();
94 daughters_name = new G4String*[numberOfD 115 daughters_name = new G4String*[numberOfDaughters];
95 // copy daughters name << 116 //copy daughters name
96 for (G4int index = 0; index < numberOfDa << 117 for (G4int index=0; index < numberOfDaughters; index++) {
97 daughters_name[index] = new G4String(* << 118 daughters_name[index] = new G4String(*right.daughters_name[index]);
98 } 119 }
99 } 120 }
100 } 121 }
101 return *this; 122 return *this;
102 } 123 }
103 124
104 G4DecayProducts* G4MuonDecayChannel::DecayIt(G << 125 G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double)
105 { 126 {
106 // this version neglects muon polarization,a << 127 // this version neglects muon polarization,and electron mass
107 // assumes the pure V-A couplin 128 // assumes the pure V-A coupling
108 // the Neutrinos are correctly << 129 // the Neutrinos are correctly V-A.
109 <<
110 #ifdef G4VERBOSE 130 #ifdef G4VERBOSE
111 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 131 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt ";
112 #endif 132 #endif
113 133
114 CheckAndFillParent(); 134 CheckAndFillParent();
115 CheckAndFillDaughters(); 135 CheckAndFillDaughters();
116 << 136
117 // parent mass 137 // parent mass
118 G4double parentmass = G4MT_parent->GetPDGMas 138 G4double parentmass = G4MT_parent->GetPDGMass();
119 const G4int N_DAUGHTER = 3; << 139 const int N_DAUGHTER=3;
120 140
121 // daughters'mass << 141 //daughters'mass
122 G4double daughtermass[N_DAUGHTER]; << 142 G4double daughtermass[N_DAUGHTER];
123 // G4double sumofdaughtermass = 0.0; << 143 G4double sumofdaughtermass = 0.0;
124 for (G4int index = 0; index < N_DAUGHTER; ++ << 144 for (G4int index=0; index<N_DAUGHTER; index++){
125 daughtermass[index] = G4MT_daughters[index 145 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
126 // sumofdaughtermass += daughtermass[index << 146 sumofdaughtermass += daughtermass[index];
127 } 147 }
128 148
129 // create parent G4DynamicParticle at rest << 149 //create parent G4DynamicParticle at rest
130 G4ThreeVector dummy; 150 G4ThreeVector dummy;
131 auto parentparticle = new G4DynamicParticle( << 151 G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
132 // create G4Decayproducts << 152 //create G4Decayproducts
133 auto products = new G4DecayProducts(*parentp << 153 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
134 delete parentparticle; 154 delete parentparticle;
135 155
136 // calculate daughter momentum 156 // calculate daughter momentum
137 G4double daughtermomentum[N_DAUGHTER]; 157 G4double daughtermomentum[N_DAUGHTER];
138 // calculate electron energy << 158 // calcurate electron energy
139 G4double xmax = (1.0 + daughtermass[0] * dau << 159 G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
140 G4double x; 160 G4double x;
141 << 161
142 G4double Ee, Ene; << 162 G4double Ee,Ene;
143 << 163
144 G4double gam; 164 G4double gam;
145 G4double EMax = parentmass / 2 - daughtermas << 165 G4double EMax=parentmass/2-daughtermass[0];
146 << 166
147 const std::size_t MAX_LOOP = 1000; << 167 const size_t MAX_LOOP=1000;
148 // Generating Random Energy << 168 //Generating Random Energy
149 for (std::size_t loop1 = 0; loop1 < MAX_LOOP << 169 for (size_t loop1=0; loop1 <MAX_LOOP; ++loop1){
150 Ee = G4UniformRand(); << 170 Ee=G4UniformRand();
151 for (std::size_t loop2 = 0; loop2 < MAX_LO << 171 for (size_t loop2 =0; loop2<MAX_LOOP; ++loop2){
152 x = xmax * G4UniformRand(); << 172 x=xmax*G4UniformRand();
153 gam = G4UniformRand(); << 173 gam=G4UniformRand();
154 if (gam <= x * (1. - x)) break; << 174 if (gam <= x*(1.-x)) break;
155 x = xmax; 175 x = xmax;
156 } 176 }
157 Ene = x; << 177 Ene=x;
158 if (Ene >= (1. - Ee)) break; << 178 if ( Ene >= (1.-Ee)) break;
159 Ene = 1. - Ee; << 179 Ene = 1.-Ee;
160 } 180 }
161 G4double Enm = (2. - Ee - Ene); << 181 G4double Enm=(2.-Ee-Ene);
>> 182
162 183
163 // initialisation of rotation parameters << 184 //initialisation of rotation parameters
164 185
165 G4double costheta, sintheta, rphi, rtheta, r << 186 G4double costheta,sintheta,rphi,rtheta,rpsi;
166 costheta = 1. - 2. / Ee - 2. / Ene + 2. / En << 187 costheta= 1.-2./Ee-2./Ene+2./Ene/Ee;
167 sintheta = std::sqrt(1. - costheta * costhet << 188 sintheta=std::sqrt(1.-costheta*costheta);
>> 189
168 190
169 rphi = twopi * G4UniformRand() * rad; << 191 rphi=twopi*G4UniformRand()*rad;
170 rtheta = (std::acos(2. * G4UniformRand() - 1 << 192 rtheta=(std::acos(2.*G4UniformRand()-1.));
171 rpsi = twopi * G4UniformRand() * rad; << 193 rpsi=twopi*G4UniformRand()*rad;
172 194
173 G4RotationMatrix rot; 195 G4RotationMatrix rot;
174 rot.set(rphi, rtheta, rpsi); << 196 rot.set(rphi,rtheta,rpsi);
175 197
176 // electron 0 << 198 //electron 0
177 daughtermomentum[0] = std::sqrt(Ee * Ee * EM << 199 daughtermomentum[0]=std::sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]);
178 G4ThreeVector direction0(0.0, 0.0, 1.0); << 200 G4ThreeVector direction0(0.0,0.0,1.0);
179 201
180 direction0 *= rot; 202 direction0 *= rot;
181 203
182 auto daughterparticle = << 204 G4DynamicParticle * daughterparticle = new G4DynamicParticle ( G4MT_daughters[0], direction0 * daughtermomentum[0]);
183 new G4DynamicParticle(G4MT_daughters[0], d <<
184 205
185 products->PushProducts(daughterparticle); 206 products->PushProducts(daughterparticle);
>> 207
>> 208 //electronic neutrino 1
186 209
187 // electronic neutrino 1 << 210 daughtermomentum[1]=std::sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]);
188 << 211 G4ThreeVector direction1(sintheta,0.0,costheta);
189 daughtermomentum[1] = std::sqrt(Ene * Ene * <<
190 G4ThreeVector direction1(sintheta, 0.0, cost <<
191 212
192 direction1 *= rot; 213 direction1 *= rot;
193 214
194 auto daughterparticle1 = << 215 G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( G4MT_daughters[1], direction1 * daughtermomentum[1]);
195 new G4DynamicParticle(G4MT_daughters[1], d <<
196 products->PushProducts(daughterparticle1); 216 products->PushProducts(daughterparticle1);
197 217
198 // muonnic neutrino 2 << 218 //muonnic neutrino 2
199 << 219
200 daughtermomentum[2] = std::sqrt(Enm * Enm * << 220 daughtermomentum[2]=std::sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]);
201 G4ThreeVector direction2(-Ene / Enm * sinthe << 221 G4ThreeVector direction2(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta);
202 222
203 direction2 *= rot; 223 direction2 *= rot;
204 224
205 auto daughterparticle2 = << 225 G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( G4MT_daughters[2],
206 new G4DynamicParticle(G4MT_daughters[2], d << 226 direction2 * daughtermomentum[2]);
207 products->PushProducts(daughterparticle2); 227 products->PushProducts(daughterparticle2);
208 228
209 // output message << 229 // output message
210 #ifdef G4VERBOSE 230 #ifdef G4VERBOSE
211 if (GetVerboseLevel() > 1) { << 231 if (GetVerboseLevel()>1) {
212 G4cout << "G4MuonDecayChannel::DecayIt()"; << 232 G4cout << "G4MuonDecayChannel::DecayIt ";
213 G4cout << " create decay products in rest << 233 G4cout << " create decay products in rest frame " <<G4endl;
214 products->DumpInfo(); 234 products->DumpInfo();
215 } 235 }
216 #endif 236 #endif
217 return products; 237 return products;
218 } 238 }
219 239