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