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