| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4PhaseSpaceDecayChannel class implementati <<
27 // 23 //
28 // Author: H.Kurashige, 27 July 1996 << 24 // $Id: G4PhaseSpaceDecayChannel.cc,v 1.7 2001/10/15 09:58:36 kurasige Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-04-00 $
30 << 26 //
31 #include "G4PhaseSpaceDecayChannel.hh" << 27 //
>> 28 // ------------------------------------------------------------
>> 29 // GEANT 4 class header file
>> 30 //
>> 31 // History: first implementation, based on object model of
>> 32 // 27 July 1996 H.Kurashige
>> 33 // 20 Oct 1996 H.Kurashige
>> 34 // 30 May 1997 H.Kurashige
>> 35 // ------------------------------------------------------------
32 36
33 #include "G4DecayProducts.hh" <<
34 #include "G4LorentzRotation.hh" <<
35 #include "G4LorentzVector.hh" <<
36 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh"
37 #include "G4PhysicalConstants.hh" << 38 #include "G4DecayProducts.hh"
38 #include "G4SystemOfUnits.hh" <<
39 #include "G4VDecayChannel.hh" 39 #include "G4VDecayChannel.hh"
>> 40 #include "G4PhaseSpaceDecayChannel.hh"
40 #include "Randomize.hh" 41 #include "Randomize.hh"
>> 42 #include "G4LorentzVector.hh"
>> 43 #include "G4LorentzRotation.hh"
>> 44
41 45
42 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayCha 46 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
43 : G4VDecayChannel("Phase Space", Verbose) << 47 :G4VDecayChannel("Phase Space", Verbose)
44 {} << 48 {
45 49
46 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayCha << 50 }
47 <<
48 <<
49 <<
50 <<
51 <<
52 <<
53 : G4VDecayChannel("Phase Space", theParentNa <<
54 theDaughterName2, theDaugh <<
55 {} <<
56 51
57 G4DecayProducts* G4PhaseSpaceDecayChannel::Dec << 52 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(
>> 53 const G4String& theParentName,
>> 54 G4double theBR,
>> 55 G4int theNumberOfDaughters,
>> 56 const G4String& theDaughterName1,
>> 57 const G4String& theDaughterName2,
>> 58 const G4String& theDaughterName3,
>> 59 const G4String& theDaughterName4 )
>> 60 :G4VDecayChannel("Phase Space",
>> 61 theParentName,theBR,
>> 62 theNumberOfDaughters,
>> 63 theDaughterName1,
>> 64 theDaughterName2,
>> 65 theDaughterName3,
>> 66 theDaughterName4)
58 { 67 {
59 #ifdef G4VERBOSE <<
60 if (GetVerboseLevel() > 1) G4cout << "G4Phas <<
61 #endif <<
62 68
63 G4DecayProducts* products = nullptr; << 69 }
64 <<
65 CheckAndFillParent(); <<
66 CheckAndFillDaughters(); <<
67 70
68 if (parentMass > 0.0) << 71 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
69 current_parent_mass.Put(parentMass); << 72 {
70 else << 73 }
71 current_parent_mass.Put(G4MT_parent_mass); <<
72 74
73 switch (numberOfDaughters) { << 75 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
74 case 0: << 76 {
75 #ifdef G4VERBOSE 77 #ifdef G4VERBOSE
76 if (GetVerboseLevel() > 0) { << 78 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
77 G4cout << "G4PhaseSpaceDecayChannel::D <<
78 G4cout << " daughters not defined " << <<
79 } <<
80 #endif 79 #endif
81 break; << 80
82 case 1: << 81 G4DecayProducts * products = 0;
83 products = OneBodyDecayIt(); << 82
84 break; << 83 if (parent == 0) FillParent();
85 case 2: << 84 if (daughters == 0) FillDaughters();
86 products = TwoBodyDecayIt(); << 85
87 break; << 86 if (parentMass >0.0) parent_mass = parentMass;
88 case 3: << 87 switch (numberOfDaughters){
89 products = ThreeBodyDecayIt(); << 88 case 0:
90 break; << 89 #ifdef G4VERBOSE
91 default: << 90 if (GetVerboseLevel()>0) {
92 products = ManyBodyDecayIt(); << 91 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
93 break; << 92 G4cout << " daughters not defined " <<G4endl;
94 } << 93 }
95 #ifdef G4VERBOSE << 94 #endif
96 if ((products == nullptr) && (GetVerboseLeve << 95 break;
97 G4cout << "G4PhaseSpaceDecayChannel::Decay << 96 case 1:
98 G4cout << *parent_name << " cannot decay " << 97 products = OneBodyDecayIt();
>> 98 break;
>> 99 case 2:
>> 100 products = TwoBodyDecayIt();
>> 101 break;
>> 102 case 3:
>> 103 products = ThreeBodyDecayIt();
>> 104 break;
>> 105 default:
>> 106 products = ManyBodyDecayIt();
>> 107 break;
>> 108 }
>> 109 #ifdef G4VERBOSE
>> 110 if ((products == 0) && (GetVerboseLevel()>0)) {
>> 111 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
>> 112 G4cout << *parent_name << " can not decay " << G4endl;
99 DumpInfo(); 113 DumpInfo();
100 } 114 }
101 #endif 115 #endif
102 return products; 116 return products;
103 } 117 }
104 118
105 G4DecayProducts* G4PhaseSpaceDecayChannel::One << 119 G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
106 { 120 {
107 #ifdef G4VERBOSE 121 #ifdef G4VERBOSE
108 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 122 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
109 #endif 123 #endif
110 // parent mass <<
111 G4double parentmass = current_parent_mass.Ge <<
112 124
113 // create parent G4DynamicParticle at rest << 125 //create parent G4DynamicParticle at rest
114 G4ThreeVector dummy; 126 G4ThreeVector dummy;
115 auto parentparticle = new G4DynamicParticle( << 127 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
116 // create G4Decayproducts << 128 //create G4Decayproducts
117 auto products = new G4DecayProducts(*parentp << 129 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
118 delete parentparticle; 130 delete parentparticle;
119 131
120 // create daughter G4DynamicParticle at rest << 132 //create daughter G4DynamicParticle at rest
121 auto daughterparticle = new G4DynamicParticl << 133 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0);
122 if (useGivenDaughterMass) daughterparticle-> <<
123 products->PushProducts(daughterparticle); 134 products->PushProducts(daughterparticle);
124 135
125 #ifdef G4VERBOSE 136 #ifdef G4VERBOSE
126 if (GetVerboseLevel() > 1) { << 137 if (GetVerboseLevel()>1) {
127 G4cout << "G4PhaseSpaceDecayChannel::OneBo << 138 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
128 G4cout << " create decay products in rest << 139 G4cout << " create decay products in rest frame " <<G4endl;
129 products->DumpInfo(); << 140 products->DumpInfo();
130 } 141 }
131 #endif 142 #endif
132 return products; 143 return products;
133 } 144 }
134 145
135 G4DecayProducts* G4PhaseSpaceDecayChannel::Two << 146 G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
136 { 147 {
137 #ifdef G4VERBOSE 148 #ifdef G4VERBOSE
138 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 149 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
139 #endif 150 #endif
140 // parent mass 151 // parent mass
141 G4double parentmass = current_parent_mass.Ge << 152 G4double parentmass = parent_mass;
>> 153
>> 154 //daughters'mass
>> 155 G4double daughtermass[2];
>> 156 G4double daughtermomentum;
>> 157 daughtermass[0] = daughters_mass[0];
>> 158 daughtermass[1] = daughters_mass[1];
142 159
143 // daughters'mass, width << 160 //create parent G4DynamicParticle at rest
144 G4double daughtermass[2], daughterwidth[2]; <<
145 daughtermass[0] = G4MT_daughters_mass[0]; <<
146 daughtermass[1] = G4MT_daughters_mass[1]; <<
147 daughterwidth[0] = G4MT_daughters_width[0]; <<
148 daughterwidth[1] = G4MT_daughters_width[1]; <<
149 <<
150 // create parent G4DynamicParticle at rest <<
151 G4ThreeVector dummy; 161 G4ThreeVector dummy;
152 auto parentparticle = new G4DynamicParticle( << 162 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
153 // create G4Decayproducts << 163 //create G4Decayproducts
154 auto products = new G4DecayProducts(*parentp << 164 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
155 delete parentparticle; 165 delete parentparticle;
156 166
157 if (!useGivenDaughterMass) { << 167 //calculate daughter momentum
158 G4bool withWidth = (daughterwidth[0] > 1.0 << 168 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
159 || (daughterwidth[1] > << 169 G4double costheta = 2.*G4UniformRand()-1.0;
160 if (withWidth) { << 170 G4double sintheta = sqrt((1.0 - costheta)*(1.0 + costheta));
161 G4double sumofdaughterwidthsq = << 171 G4double phi = 2.0*M_PI*G4UniformRand()*rad;
162 daughterwidth[0] * daughterwidth[0] + << 172 G4ThreeVector direction(sintheta*cos(phi),sintheta*sin(phi),costheta);
163 // check parent mass and daughter mass << 173
164 G4double maxDev = << 174 //create daughter G4DynamicParticle
165 (parentmass - daughtermass[0] - daught << 175 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum);
166 if (maxDev <= -1.0 * rangeMass) { <<
167 #ifdef G4VERBOSE <<
168 if (GetVerboseLevel() > 0) { <<
169 G4cout << "G4PhaseSpaceDecayChannel: <<
170 << "Sum of daughter mass is l <<
171 G4cout << "Parent :" << G4MT_parent- <<
172 << current_parent_mass.Get() <<
173 G4cout << "Daughter 1 :" << G4MT_dau <<
174 << daughtermass[0] / GeV << G <<
175 G4cout << "Daughter 2:" << G4MT_daug <<
176 << daughtermass[1] / GeV << G <<
177 } <<
178 #endif <<
179 G4Exception("G4PhaseSpaceDecayChannel: <<
180 "Cannot create decay produ <<
181 larger than parent mass!" <<
182 return products; <<
183 } <<
184 G4double dm1 = daughtermass[0]; <<
185 if (daughterwidth[0] > 0.) dm1 = Dynamic <<
186 G4double dm2 = daughtermass[1]; <<
187 if (daughterwidth[1] > 0.) dm2 = Dynamic <<
188 while (dm1 + dm2 > parentmass) // Loop <<
189 { <<
190 dm1 = DynamicalMass(daughtermass[0], d <<
191 dm2 = DynamicalMass(daughtermass[1], d <<
192 } <<
193 daughtermass[0] = dm1; <<
194 daughtermass[1] = dm2; <<
195 } <<
196 } <<
197 else { <<
198 // use given daughter mass; <<
199 daughtermass[0] = givenDaughterMasses[0]; <<
200 daughtermass[1] = givenDaughterMasses[1]; <<
201 } <<
202 if (parentmass < daughtermass[0] + daughterm <<
203 #ifdef G4VERBOSE <<
204 if (GetVerboseLevel() > 0) { <<
205 G4cout << "G4PhaseSpaceDecayChannel::Two <<
206 << "Sum of daughter mass is large <<
207 G4cout << "Parent :" << G4MT_parent->Get <<
208 << current_parent_mass.Get() / Ge <<
209 G4cout << "Daughter 1 :" << G4MT_daughte <<
210 << daughtermass[0] / GeV << G4end <<
211 G4cout << "Daughter 2:" << G4MT_daughter <<
212 << daughtermass[1] / GeV << G4end <<
213 if (useGivenDaughterMass) { <<
214 G4cout << "Daughter Mass is given." << <<
215 } <<
216 } <<
217 #endif <<
218 G4Exception("G4PhaseSpaceDecayChannel::Two <<
219 "Cannot create decay products: <<
220 larger than parent mass!"); <<
221 return products; <<
222 } <<
223 <<
224 // calculate daughter momentum <<
225 G4double daughtermomentum = Pmx(parentmass, <<
226 <<
227 G4double costheta = 2. * G4UniformRand() - 1 <<
228 G4double sintheta = std::sqrt((1.0 - costhet <<
229 G4double phi = twopi * G4UniformRand() * rad <<
230 G4ThreeVector direction(sintheta * std::cos( <<
231 <<
232 // create daughter G4DynamicParticle <<
233 G4double Ekin = std::sqrt(daughtermomentum * <<
234 - daughtermass[0]; <<
235 auto daughterparticle = <<
236 new G4DynamicParticle(G4MT_daughters[0], d <<
237 products->PushProducts(daughterparticle); 176 products->PushProducts(daughterparticle);
238 Ekin = std::sqrt(daughtermomentum * daughter << 177 daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum));
239 - daughtermass[1]; <<
240 daughterparticle = <<
241 new G4DynamicParticle(G4MT_daughters[1], - <<
242 products->PushProducts(daughterparticle); 178 products->PushProducts(daughterparticle);
243 179
244 #ifdef G4VERBOSE 180 #ifdef G4VERBOSE
245 if (GetVerboseLevel() > 1) { << 181 if (GetVerboseLevel()>1) {
246 G4cout << "G4PhaseSpaceDecayChannel::TwoBo << 182 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
247 G4cout << " Create decay products in rest << 183 G4cout << " create decay products in rest frame " <<G4endl;
248 products->DumpInfo(); << 184 products->DumpInfo();
249 } 185 }
250 #endif 186 #endif
251 return products; 187 return products;
252 } 188 }
253 189
254 G4DecayProducts* G4PhaseSpaceDecayChannel::Thr << 190 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
>> 191 // algorism of this code is originally written in GDECA3 of GEANT3
255 { 192 {
256 // Algorithm of this code originally written <<
257 <<
258 #ifdef G4VERBOSE 193 #ifdef G4VERBOSE
259 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 194 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
260 #endif 195 #endif
261 // parent mass 196 // parent mass
262 G4double parentmass = current_parent_mass.Ge << 197 G4double parentmass = parent_mass;
263 // daughters'mass << 198 //daughters'mass
264 G4double daughtermass[3], daughterwidth[3]; << 199 G4double daughtermass[3];
265 G4double sumofdaughtermass = 0.0; 200 G4double sumofdaughtermass = 0.0;
266 G4double sumofdaughterwidthsq = 0.0; << 201 for (G4int index=0; index<3; index++){
267 G4bool withWidth = false; << 202 daughtermass[index] = daughters_mass[index];
268 for (G4int index = 0; index < 3; ++index) { <<
269 daughtermass[index] = G4MT_daughters_mass[ <<
270 sumofdaughtermass += daughtermass[index]; 203 sumofdaughtermass += daughtermass[index];
271 daughterwidth[index] = G4MT_daughters_widt <<
272 sumofdaughterwidthsq += daughterwidth[inde <<
273 withWidth = withWidth || (daughterwidth[in <<
274 } 204 }
275 << 205
276 // create parent G4DynamicParticle at rest << 206 //create parent G4DynamicParticle at rest
277 G4ThreeVector dummy; 207 G4ThreeVector dummy;
278 auto parentparticle = new G4DynamicParticle( << 208 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
279 << 209 //create G4Decayproducts
280 // create G4Decayproducts << 210 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
281 auto products = new G4DecayProducts(*parentp <<
282 delete parentparticle; 211 delete parentparticle;
283 212
284 if (!useGivenDaughterMass) { << 213 //calculate daughter momentum
285 if (withWidth) { << 214 // Generate two
286 G4double maxDev = (parentmass - sumofdau <<
287 if (maxDev <= -1.0 * rangeMass) { <<
288 #ifdef G4VERBOSE <<
289 if (GetVerboseLevel() > 0) { <<
290 G4cout << "G4PhaseSpaceDecayChannel: <<
291 << "Sum of daughter mass is l <<
292 G4cout << "Parent :" << G4MT_parent- <<
293 << current_parent_mass.Get() <<
294 G4cout << "Daughter 1 :" << G4MT_dau <<
295 << daughtermass[0] / GeV << G <<
296 G4cout << "Daughter 2:" << G4MT_daug <<
297 << daughtermass[1] / GeV << G <<
298 G4cout << "Daughter 3:" << G4MT_daug <<
299 << daughtermass[2] / GeV << G <<
300 } <<
301 #endif <<
302 G4Exception("G4PhaseSpaceDecayChannel: <<
303 "Cannot create decay produ <<
304 is larger than parent mas <<
305 return products; <<
306 } <<
307 G4double dm1 = daughtermass[0]; <<
308 if (daughterwidth[0] > 0.) dm1 = Dynamic <<
309 G4double dm2 = daughtermass[1]; <<
310 if (daughterwidth[1] > 0.) dm2 = Dynamic <<
311 G4double dm3 = daughtermass[2]; <<
312 if (daughterwidth[2] > 0.) dm3 = Dynamic <<
313 while (dm1 + dm2 + dm3 > parentmass) // <<
314 { <<
315 dm1 = DynamicalMass(daughtermass[0], d <<
316 dm2 = DynamicalMass(daughtermass[1], d <<
317 dm3 = DynamicalMass(daughtermass[2], d <<
318 } <<
319 daughtermass[0] = dm1; <<
320 daughtermass[1] = dm2; <<
321 daughtermass[2] = dm3; <<
322 sumofdaughtermass = dm1 + dm2 + dm3; <<
323 } <<
324 } <<
325 else { <<
326 // use given daughter mass; <<
327 daughtermass[0] = givenDaughterMasses[0]; <<
328 daughtermass[1] = givenDaughterMasses[1]; <<
329 daughtermass[2] = givenDaughterMasses[2]; <<
330 sumofdaughtermass = daughtermass[0] + daug <<
331 } <<
332 <<
333 if (sumofdaughtermass > parentmass) { <<
334 #ifdef G4VERBOSE <<
335 if (GetVerboseLevel() > 0) { <<
336 G4cout << "G4PhaseSpaceDecayChannel::Thr <<
337 << "Sum of daughter mass is large <<
338 G4cout << "Parent :" << G4MT_parent->Get <<
339 << current_parent_mass.Get() / Ge <<
340 G4cout << "Daughter 1 :" << G4MT_daughte <<
341 << daughtermass[0] / GeV << G4end <<
342 G4cout << "Daughter 2:" << G4MT_daughter <<
343 << daughtermass[1] / GeV << G4end <<
344 G4cout << "Daughter 3:" << G4MT_daughter <<
345 << daughtermass[2] / GeV << G4end <<
346 } <<
347 if (useGivenDaughterMass) { <<
348 G4cout << "Daughter Mass is given." << G <<
349 } <<
350 #endif <<
351 G4Exception("G4PhaseSpaceDecayChannel::Thr <<
352 "Can not create decay products <<
353 is larger than parent mass!") <<
354 return products; <<
355 } <<
356 <<
357 // calculate daughter momentum <<
358 // Generate two <<
359 G4double rd1, rd2, rd; 215 G4double rd1, rd2, rd;
360 G4double daughtermomentum[3]; 216 G4double daughtermomentum[3];
361 G4double momentummax = 0.0, momentumsum = 0. << 217 G4double momentummax=0.0, momentumsum = 0.0;
362 G4double energy; 218 G4double energy;
363 const std::size_t MAX_LOOP = 10000; <<
364 219
365 for (std::size_t loop_counter = 0; loop_coun << 220 do {
366 rd1 = G4UniformRand(); 221 rd1 = G4UniformRand();
367 rd2 = G4UniformRand(); 222 rd2 = G4UniformRand();
368 if (rd2 > rd1) { 223 if (rd2 > rd1) {
369 rd = rd1; << 224 rd = rd1;
370 rd1 = rd2; 225 rd1 = rd2;
371 rd2 = rd; 226 rd2 = rd;
372 } << 227 }
373 momentummax = 0.0; 228 momentummax = 0.0;
374 momentumsum = 0.0; 229 momentumsum = 0.0;
375 // daughter 0 230 // daughter 0
376 energy = rd2 * (parentmass - sumofdaughter << 231 energy = rd2*(parentmass - sumofdaughtermass);
377 daughtermomentum[0] = std::sqrt(energy * e << 232 daughtermomentum[0] = sqrt(energy*energy + 2.0*energy* daughtermass[0]);
378 if (daughtermomentum[0] > momentummax) mom << 233 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
379 momentumsum += daughtermomentum[0]; << 234 momentumsum += daughtermomentum[0];
380 // daughter 1 235 // daughter 1
381 energy = (1. - rd1) * (parentmass - sumofd << 236 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
382 daughtermomentum[1] = std::sqrt(energy * e << 237 daughtermomentum[1] = sqrt(energy*energy + 2.0*energy* daughtermass[1]);
383 if (daughtermomentum[1] > momentummax) mom << 238 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
384 momentumsum += daughtermomentum[1]; << 239 momentumsum += daughtermomentum[1];
385 // daughter 2 240 // daughter 2
386 energy = (rd1 - rd2) * (parentmass - sumof << 241 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
387 daughtermomentum[2] = std::sqrt(energy * e << 242 daughtermomentum[2] = sqrt(energy*energy + 2.0*energy* daughtermass[2]);
388 if (daughtermomentum[2] > momentummax) mom << 243 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
389 momentumsum += daughtermomentum[2]; << 244 momentumsum += daughtermomentum[2];
390 if (momentummax <= momentumsum - momentumm << 245 } while (momentummax > momentumsum - momentummax );
391 } <<
392 <<
393 // output message 246 // output message
394 #ifdef G4VERBOSE 247 #ifdef G4VERBOSE
395 if (GetVerboseLevel() > 1) { << 248 if (GetVerboseLevel()>1) {
396 G4cout << " daughter 0:" << daughtermo << 249 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
397 G4cout << " daughter 1:" << daughtermo << 250 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
398 G4cout << " daughter 2:" << daughtermo << 251 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
399 G4cout << " momentum sum:" << momentumsu << 252 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
400 } 253 }
401 #endif 254 #endif
402 255
403 // create daughter G4DynamicParticle << 256 //create daughter G4DynamicParticle
404 G4double costheta, sintheta, phi, sinphi, co << 257 G4double costheta, sintheta, phi, sinphi, cosphi;
405 G4double costhetan, sinthetan, phin, sinphin << 258 G4double costhetan, sinthetan, phin, sinphin, cosphin;
406 costheta = 2. * G4UniformRand() - 1.0; << 259 costheta = 2.*G4UniformRand()-1.0;
407 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 260 sintheta = sqrt((1.0-costheta)*(1.0+costheta));
408 phi = twopi * G4UniformRand() * rad; << 261 phi = 2.0*M_PI*G4UniformRand()*rad;
409 sinphi = std::sin(phi); << 262 sinphi = sin(phi);
410 cosphi = std::cos(phi); << 263 cosphi = cos(phi);
411 << 264 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
412 G4ThreeVector direction0(sintheta * cosphi, << 265 G4DynamicParticle * daughterparticle
413 G4double Ekin = << 266 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
414 std::sqrt(daughtermomentum[0] * daughtermo <<
415 - daughtermass[0]; <<
416 auto daughterparticle = <<
417 new G4DynamicParticle(G4MT_daughters[0], d <<
418 products->PushProducts(daughterparticle); 267 products->PushProducts(daughterparticle);
419 268
420 costhetan = (daughtermomentum[1] * daughterm << 269 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
421 - daughtermomentum[0] * daughte << 270 sinthetan = sqrt((1.0-costhetan)*(1.0+costhetan));
422 / (2.0 * daughtermomentum[2] * d << 271 phin = 2.0*M_PI*G4UniformRand()*rad;
423 sinthetan = std::sqrt((1.0 - costhetan) * (1 << 272 sinphin = sin(phin);
424 phin = twopi * G4UniformRand() * rad; << 273 cosphin = cos(phin);
425 sinphin = std::sin(phin); <<
426 cosphin = std::cos(phin); <<
427 G4ThreeVector direction2; 274 G4ThreeVector direction2;
428 direction2.setX(sinthetan * cosphin * costhe << 275 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
429 + costhetan * sintheta * cos << 276 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
430 direction2.setY(sinthetan * cosphin * costhe << 277 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
431 + costhetan * sintheta * sin << 278 daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
432 direction2.setZ(-sinthetan * cosphin * sinth <<
433 G4ThreeVector pmom = daughtermomentum[2] * d <<
434 Ekin = std::sqrt(pmom.mag2() + daughtermass[ <<
435 daughterparticle = <<
436 new G4DynamicParticle(G4MT_daughters[2], p <<
437 products->PushProducts(daughterparticle); 279 products->PushProducts(daughterparticle);
438 280
439 pmom = (direction0 * daughtermomentum[0] + d << 281 daughterparticle =
440 * (-1.0); << 282 new G4DynamicParticle(
441 Ekin = std::sqrt(pmom.mag2() + daughtermass[ << 283 daughters[1],
442 daughterparticle = << 284 (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)
443 new G4DynamicParticle(G4MT_daughters[1], p << 285 );
444 products->PushProducts(daughterparticle); 286 products->PushProducts(daughterparticle);
445 287
446 #ifdef G4VERBOSE 288 #ifdef G4VERBOSE
447 if (GetVerboseLevel() > 1) { << 289 if (GetVerboseLevel()>1) {
448 G4cout << "G4PhaseSpaceDecayChannel::Three << 290 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
449 G4cout << " create decay products in rest << 291 G4cout << " create decay products in rest frame " <<G4endl;
450 products->DumpInfo(); << 292 products->DumpInfo();
451 } 293 }
452 #endif 294 #endif
453 return products; 295 return products;
454 } 296 }
455 297
456 G4DecayProducts* G4PhaseSpaceDecayChannel::Man << 298 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
>> 299 // algorithm of this code is originally written in FORTRAN by M.Asai
>> 300 //******************************************************************
>> 301 // NBODY
>> 302 // N-body phase space Monte-Carlo generator
>> 303 // Makoto Asai
>> 304 // Hiroshima Institute of Technology
>> 305 // (asai@kekvax.kek.jp)
>> 306 // Revised release : 19/Apr/1995
>> 307 //
457 { 308 {
458 // Algorithm of this code originally written <<
459 // ***************************************** <<
460 // NBODY - N-body phase space Monte-Carlo ge <<
461 // Makoto Asai - Hiroshima Institute of Tech <<
462 // Revised release : 19/Apr/1995 <<
463 <<
464 G4int index, index2; 309 G4int index, index2;
465 310
>> 311 //return value
>> 312 G4DecayProducts *products;
>> 313
466 #ifdef G4VERBOSE 314 #ifdef G4VERBOSE
467 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 315 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
468 #endif 316 #endif
469 <<
470 // parent mass 317 // parent mass
471 G4double parentmass = current_parent_mass.Ge << 318 G4double parentmass = parent_mass;
472 << 319 //daughters'mass
473 // parent particle << 320 G4double *daughtermass = new G4double[numberOfDaughters];
474 G4ThreeVector dummy; <<
475 auto parentparticle = new G4DynamicParticle( <<
476 <<
477 // create G4Decayproducts <<
478 auto products = new G4DecayProducts(*parentp <<
479 delete parentparticle; <<
480 <<
481 // daughters'mass <<
482 auto daughtermass = new G4double[numberOfDau <<
483 <<
484 G4double sumofdaughtermass = 0.0; 321 G4double sumofdaughtermass = 0.0;
485 for (index = 0; index < numberOfDaughters; + << 322 for (index=0; index<numberOfDaughters; index++){
486 if (!useGivenDaughterMass) { << 323 daughtermass[index] = daughters_mass[index];
487 daughtermass[index] = G4MT_daughters_mas <<
488 } <<
489 else { <<
490 daughtermass[index] = givenDaughterMasse <<
491 } <<
492 sumofdaughtermass += daughtermass[index]; 324 sumofdaughtermass += daughtermass[index];
493 } 325 }
494 if (sumofdaughtermass > parentmass) { << 326
495 #ifdef G4VERBOSE << 327 //Calculate daughter momentum
496 if (GetVerboseLevel() > 0) { << 328 G4double *daughtermomentum = new G4double[numberOfDaughters];
497 G4cout << "G4PhaseSpaceDecayChannel::Man << 329 G4ThreeVector direction;
498 << "Sum of daughter mass is large << 330 G4DynamicParticle **daughterparticle;
499 G4cout << "Parent :" << G4MT_parent->Get << 331 G4double *sm = new G4double[numberOfDaughters];
500 << current_parent_mass.Get() / Ge <<
501 G4cout << "Daughter 1 :" << G4MT_daughte <<
502 << daughtermass[0] / GeV << G4end <<
503 G4cout << "Daughter 2:" << G4MT_daughter <<
504 << daughtermass[1] / GeV << G4end <<
505 G4cout << "Daughter 3:" << G4MT_daughter <<
506 << daughtermass[2] / GeV << G4end <<
507 G4cout << "Daughter 4:" << G4MT_daughter <<
508 << daughtermass[3] / GeV << G4end <<
509 G4cout << "Daughter 5:" << G4MT_daughter <<
510 << daughtermass[4] / GeV << G4end <<
511 } <<
512 #endif <<
513 G4Exception("G4PhaseSpaceDecayChannel::Man <<
514 "Can not create decay products <<
515 is larger than parent mass!") <<
516 delete[] daughtermass; <<
517 return products; <<
518 } <<
519 <<
520 // Calculate daughter momentum <<
521 auto daughtermomentum = new G4double[numberO <<
522 G4ThreeVector direction; <<
523 G4DynamicParticle** daughterparticle; <<
524 auto sm = new G4double[numberOfDaughters]; <<
525 G4double tmas; 332 G4double tmas;
526 G4double weight = 1.0; 333 G4double weight = 1.0;
527 G4int numberOfTry = 0; << 334 G4int numberOfTry = 0;
528 const std::size_t MAX_LOOP = 10000; << 335
529 << 336 do {
530 for (std::size_t loop_counter = 0; loop_coun << 337 //Generate rundom number in descending order
531 // Generate rundom number in descending or <<
532 G4double temp; 338 G4double temp;
533 auto rd = new G4double[numberOfDaughters]; << 339 G4double *rd = new G4double[numberOfDaughters];
534 rd[0] = 1.0; 340 rd[0] = 1.0;
535 for (index = 1; index < numberOfDaughters << 341 for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
536 rd[index] = G4UniformRand(); << 342 rd[ numberOfDaughters -1] = 0.0;
537 } << 343 for(index =1; index < numberOfDaughters -1; index++) {
538 rd[numberOfDaughters - 1] = 0.0; << 344 for(index2 = index+1; index2 < numberOfDaughters; index2++) {
539 for (index = 1; index < numberOfDaughters << 345 if (rd[index] < rd[index2]){
540 for (index2 = index + 1; index2 < number << 346 temp = rd[index];
541 if (rd[index] < rd[index2]) { << 347 rd[index] = rd[index2];
542 temp = rd[index]; <<
543 rd[index] = rd[index2]; <<
544 rd[index2] = temp; 348 rd[index2] = temp;
545 } 349 }
546 } 350 }
547 } 351 }
548 << 352
549 // calculate virtual mass << 353 //calcurate virtual mass
550 tmas = parentmass - sumofdaughtermass; << 354 tmas = parentmass - sumofdaughtermass;
551 temp = sumofdaughtermass; << 355 temp = sumofdaughtermass;
552 for (index = 0; index < numberOfDaughters; << 356 for(index =0; index < numberOfDaughters; index++) {
553 sm[index] = rd[index] * tmas + temp; << 357 sm[index] = rd[index]*tmas + temp;
554 temp -= daughtermass[index]; 358 temp -= daughtermass[index];
555 if (GetVerboseLevel() > 1) { << 359 if (GetVerboseLevel()>1) {
556 G4cout << index << " rundom number:" << 360 G4cout << index << " rundom number:" << rd[index];
557 G4cout << " virtual mass:" << sm[ind << 361 G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
558 } 362 }
559 } 363 }
560 delete[] rd; << 364 delete [] rd;
561 365
562 // Calculate daughter momentum << 366 //Calculate daughter momentum
563 weight = 1.0; 367 weight = 1.0;
564 G4bool smOK = true; << 368 index =numberOfDaughters-1;
565 for (index = 0; index < numberOfDaughters << 369 daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
566 smOK = (sm[index] - daughtermass[index] <<
567 } <<
568 if (!smOK) continue; <<
569 <<
570 index = numberOfDaughters - 1; <<
571 daughtermomentum[index] = Pmx(sm[index - 1 <<
572 #ifdef G4VERBOSE 370 #ifdef G4VERBOSE
573 if (GetVerboseLevel() > 1) { << 371 if (GetVerboseLevel()>1) {
574 G4cout << " daughter " << index << " 372 G4cout << " daughter " << index << ":" << *daughters_name[index];
575 G4cout << " momentum:" << daughtermoment << 373 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
576 } 374 }
577 #endif 375 #endif
578 for (index = numberOfDaughters - 2; index << 376 for(index =numberOfDaughters-2; index>=0; index--) {
579 // calculate << 377 // calculate
580 daughtermomentum[index] = Pmx(sm[index], << 378 daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
581 if (daughtermomentum[index] < 0.0) { << 379 if(daughtermomentum[index] < 0.0) {
582 // !!! illegal momentum !!! 380 // !!! illegal momentum !!!
583 #ifdef G4VERBOSE 381 #ifdef G4VERBOSE
584 if (GetVerboseLevel() > 0) { << 382 if (GetVerboseLevel()>0) {
585 G4cout << "G4PhaseSpaceDecayChannel: << 383 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
586 G4cout << " Cannot calculate dau << 384 G4cout << " can not calculate daughter momentum " <<G4endl;
587 G4cout << " parent:" << *parent_ 385 G4cout << " parent:" << *parent_name;
588 G4cout << " mass:" << parentmass / G << 386 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
589 G4cout << " daughter " << index 387 G4cout << " daughter " << index << ":" << *daughters_name[index];
590 G4cout << " mass:" << daughtermass[i << 388 G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
591 G4cout << " mass:" << daughtermoment << 389 G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
592 if (useGivenDaughterMass) { <<
593 G4cout << "Daughter Mass is given. <<
594 } <<
595 } 390 }
596 #endif 391 #endif
597 delete[] sm; << 392 delete [] sm;
598 delete[] daughtermass; << 393 delete [] daughtermass;
599 delete[] daughtermomentum; << 394 delete [] daughtermomentum;
600 delete products; << 395 return 0; // Error detection
601 G4Exception("G4PhaseSpaceDecayChannel: << 396
602 "Can not create decay prod << 397 } else {
603 is larger than parent mas << 398 // calculate weight of this events
604 return nullptr; // Error detection << 399 weight *= daughtermomentum[index]/sm[index];
605 } << 400 if (GetVerboseLevel()>1) {
606 << 401 G4cout << " daughter " << index << ":" << *daughters_name[index];
607 // calculate weight of this events << 402 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
608 weight *= daughtermomentum[index] / sm[i << 403 }
609 #ifdef G4VERBOSE <<
610 if (GetVerboseLevel() > 1) { <<
611 G4cout << " daughter " << index << <<
612 G4cout << " momentum:" << daughtermome <<
613 } 404 }
614 #endif <<
615 } 405 }
616 << 406 if (GetVerboseLevel()>1) {
617 #ifdef G4VERBOSE << 407 G4cout << " weight: " << weight <<G4endl;
618 if (GetVerboseLevel() > 1) { <<
619 G4cout << " weight: " << weight << G4 <<
620 } 408 }
621 #endif << 409
622 <<
623 // exit if number of Try exceeds 100 410 // exit if number of Try exceeds 100
624 if (++numberOfTry > 100) { << 411 if (numberOfTry++ >100) {
625 #ifdef G4VERBOSE << 412 if (GetVerboseLevel()>0) {
626 if (GetVerboseLevel() > 0) { << 413 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
627 G4cout << "G4PhaseSpaceDecayChannel::M << 414 G4cout << " can not determine Decay Kinematics " << G4endl;
628 G4cout << "Cannot determine Decay Kine <<
629 G4cout << "parent : " << *parent_name <<
630 G4cout << "daughters : "; <<
631 for (index = 0; index < numberOfDaught <<
632 G4cout << *daughters_name[index] << <<
633 } <<
634 G4cout << G4endl; <<
635 } 415 }
636 #endif << 416 delete [] sm;
637 G4Exception("G4PhaseSpaceDecayChannel::M << 417 delete [] daughtermass;
638 "Cannot decay: Decay Kinemat << 418 delete [] daughtermomentum;
639 << 419 return 0; // Error detection
640 delete[] sm; <<
641 delete[] daughtermass; <<
642 delete[] daughtermomentum; <<
643 delete products; <<
644 return nullptr; // Error detection <<
645 } 420 }
646 if (weight < G4UniformRand()) break; << 421 } while ( weight > G4UniformRand());
647 } <<
648 422
649 #ifdef G4VERBOSE 423 #ifdef G4VERBOSE
650 if (GetVerboseLevel() > 1) { << 424 if (GetVerboseLevel()>1) {
651 G4cout << "Start calculation of daughters << 425 G4cout << "Start calulation of daughters momentum vector "<<G4endl;
652 } 426 }
653 #endif << 427 #endif
654 <<
655 G4double costheta, sintheta, phi; 428 G4double costheta, sintheta, phi;
656 G4double beta; 429 G4double beta;
657 daughterparticle = new G4DynamicParticle*[nu 430 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
658 431
659 index = numberOfDaughters - 2; << 432 index = numberOfDaughters -2;
660 costheta = 2. * G4UniformRand() - 1.0; << 433 costheta = 2.*G4UniformRand()-1.0;
661 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 434 sintheta = sqrt((1.0-costheta)*(1.0+costheta));
662 phi = twopi * G4UniformRand() * rad; << 435 phi = 2.0*M_PI*G4UniformRand()*rad;
663 direction.setZ(costheta); 436 direction.setZ(costheta);
664 direction.setY(sintheta * std::sin(phi)); << 437 direction.setY(sintheta*sin(phi));
665 direction.setX(sintheta * std::cos(phi)); << 438 direction.setX(sintheta*cos(phi));
666 daughterparticle[index] = << 439 daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] );
667 new G4DynamicParticle(G4MT_daughters[index << 440 daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
668 daughterparticle[index + 1] = << 441
669 new G4DynamicParticle(G4MT_daughters[index << 442 for (index = numberOfDaughters -3; index >= 0; index--) {
670 << 443 //calculate momentum direction
671 for (index = numberOfDaughters - 3; index >= << 444 costheta = 2.*G4UniformRand()-1.0;
672 // calculate momentum direction << 445 sintheta = sqrt((1.0-costheta)*(1.0+costheta));
673 costheta = 2. * G4UniformRand() - 1.0; << 446 phi = 2.0*M_PI*G4UniformRand()*rad;
674 sintheta = std::sqrt((1.0 - costheta) * (1 <<
675 phi = twopi * G4UniformRand() * rad; <<
676 direction.setZ(costheta); 447 direction.setZ(costheta);
677 direction.setY(sintheta * std::sin(phi)); << 448 direction.setY(sintheta*sin(phi));
678 direction.setX(sintheta * std::cos(phi)); << 449 direction.setX(sintheta*cos(phi));
679 450
680 // boost already created particles << 451 // boost already created particles
681 beta = daughtermomentum[index]; 452 beta = daughtermomentum[index];
682 beta /= << 453 beta /= sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
683 std::sqrt(daughtermomentum[index] * daug << 454 for (G4int index2 = index+1; index2<numberOfDaughters; index2++) {
684 for (index2 = index + 1; index2 < numberOf <<
685 G4LorentzVector p4; 455 G4LorentzVector p4;
686 // make G4LorentzVector for secondaries 456 // make G4LorentzVector for secondaries
687 p4 = daughterparticle[index2]->Get4Momen 457 p4 = daughterparticle[index2]->Get4Momentum();
688 458
689 // boost secondaries to new frame << 459 // boost secondaries to new frame
690 p4.boost(direction.x() * beta, direction << 460 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
691 461
692 // change energy/momentum 462 // change energy/momentum
693 daughterparticle[index2]->Set4Momentum(p 463 daughterparticle[index2]->Set4Momentum(p4);
694 } 464 }
695 // create daughter G4DynamicParticle << 465 //create daughter G4DynamicParticle
696 daughterparticle[index] = << 466 daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index]));
697 new G4DynamicParticle(G4MT_daughters[ind <<
698 } 467 }
699 468
700 // add daughters to G4Decayproducts << 469 //create G4Decayproducts
701 for (index = 0; index < numberOfDaughters; + << 470 G4DynamicParticle *parentparticle;
>> 471 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
>> 472 parentparticle = new G4DynamicParticle( parent, direction, 0.0);
>> 473 products = new G4DecayProducts(*parentparticle);
>> 474 delete parentparticle;
>> 475 for (index = 0; index<numberOfDaughters; index++) {
702 products->PushProducts(daughterparticle[in 476 products->PushProducts(daughterparticle[index]);
703 } 477 }
704 <<
705 #ifdef G4VERBOSE 478 #ifdef G4VERBOSE
706 if (GetVerboseLevel() > 1) { << 479 if (GetVerboseLevel()>1) {
707 G4cout << "G4PhaseSpaceDecayChannel::ManyB << 480 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
708 G4cout << " create decay products in rest << 481 G4cout << " create decay products in rest frame " << G4endl;
709 products->DumpInfo(); 482 products->DumpInfo();
710 } 483 }
711 #endif 484 #endif
712 << 485 delete [] daughterparticle;
713 delete[] daughterparticle; << 486 delete [] daughtermomentum;
714 delete[] daughtermomentum; << 487 delete [] daughtermass;
715 delete[] daughtermass; << 488 delete [] sm;
716 delete[] sm; << 489
717 <<
718 return products; 490 return products;
719 } 491 }
720 492
721 G4bool G4PhaseSpaceDecayChannel::SetDaughterMa <<
722 { <<
723 for (G4int idx = 0; idx < numberOfDaughters; <<
724 givenDaughterMasses[idx] = masses[idx]; <<
725 } <<
726 useGivenDaughterMass = true; <<
727 return useGivenDaughterMass; <<
728 } <<
729 493
730 G4bool G4PhaseSpaceDecayChannel::SampleDaughte << 494 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
731 { 495 {
732 useGivenDaughterMass = false; << 496 // calcurate momentum of daughter particles in two-body decay
733 return useGivenDaughterMass; << 497 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
>> 498 if (ppp>0) return sqrt(ppp);
>> 499 else return -1.;
734 } 500 }
735 501
736 G4bool G4PhaseSpaceDecayChannel::IsOKWithParen <<
737 { <<
738 if (!useGivenDaughterMass) return G4VDecayCh <<
739 502
740 CheckAndFillParent(); <<
741 CheckAndFillDaughters(); <<
742 503
743 G4double sumOfDaughterMassMin = 0.0; <<
744 for (G4int index = 0; index < numberOfDaught <<
745 sumOfDaughterMassMin += givenDaughterMasse <<
746 } <<
747 return (parentMass >= sumOfDaughterMassMin); <<
748 } <<
749 <<
750 G4double G4PhaseSpaceDecayChannel::Pmx(G4doubl <<
751 { <<
752 // calcurate momentum of daughter particles <<
753 G4double ppp = (e + p1 + p2) * (e + p1 - p2) <<
754 if (ppp > 0) return std::sqrt(ppp); <<
755 return -1.; <<
756 } <<
757 504