| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // ----------------------------------------------------------------
28 // GEANT 4 class header file 28 // GEANT 4 class header file
29 // 29 //
30 // History: first implementation, A. Feli 30 // History: first implementation, A. Feliciello, 21st May 1998
31 // 31 //
32 // Note: this class is a generalization o 32 // Note: this class is a generalization of the
33 // G4PhaseSpaceDecayChannel one 33 // G4PhaseSpaceDecayChannel one
34 // ------------------------------------------- 34 // ----------------------------------------------------------------
35 35
36 #include "G4ParticleDefinition.hh" 36 #include "G4ParticleDefinition.hh"
37 #include "G4DecayProducts.hh" 37 #include "G4DecayProducts.hh"
38 #include "G4VDecayChannel.hh" 38 #include "G4VDecayChannel.hh"
39 #include "G4GeneralPhaseSpaceDecay.hh" 39 #include "G4GeneralPhaseSpaceDecay.hh"
40 #include "G4PhysicalConstants.hh" 40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh"
42 #include "Randomize.hh" 42 #include "Randomize.hh"
43 #include "G4LorentzVector.hh" 43 #include "G4LorentzVector.hh"
44 #include "G4LorentzRotation.hh" 44 #include "G4LorentzRotation.hh"
45 #include "G4ios.hh" 45 #include "G4ios.hh"
46 46
47 47
48 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD 48 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(G4int Verbose) :
49 G4VDecayChannel("Pha 49 G4VDecayChannel("Phase Space", Verbose),
50 parentmass(0.), theD 50 parentmass(0.), theDaughterMasses(0)
51 { 51 {
52 if (GetVerboseLevel()>1) G4cout << "G4Genera 52 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
53 } 53 }
54 54
55 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD 55 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
56 G4double 56 G4double theBR,
57 G4int 57 G4int theNumberOfDaughters,
58 const G4Strin 58 const G4String& theDaughterName1,
59 const G4Strin 59 const G4String& theDaughterName2,
60 const G4Strin 60 const G4String& theDaughterName3) :
61 G4VDecayCha 61 G4VDecayChannel("Phase Space",
62 theParentName,theBR, 62 theParentName,theBR,
63 theNumberOfDaughters, 63 theNumberOfDaughters,
64 theDaughterName1, 64 theDaughterName1,
65 theDaughterName2, 65 theDaughterName2,
66 theDaughterName3), 66 theDaughterName3),
67 theDaughterMasses(0) 67 theDaughterMasses(0)
68 { 68 {
69 if (GetVerboseLevel()>1) G4cout << "G4Genera 69 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
70 70
71 // Set the parent particle (resonance) mas 71 // Set the parent particle (resonance) mass to the (default) PDG vale
72 if (G4MT_parent != NULL) 72 if (G4MT_parent != NULL)
73 { 73 {
74 parentmass = G4MT_parent->GetPDGMass(); 74 parentmass = G4MT_parent->GetPDGMass();
75 } else { 75 } else {
76 parentmass=0.; 76 parentmass=0.;
77 } 77 }
78 78
79 } 79 }
80 80
81 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD 81 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
82 82 G4double theParentMass,
83 G4double 83 G4double theBR,
84 G4int 84 G4int theNumberOfDaughters,
85 const G4Strin 85 const G4String& theDaughterName1,
86 const G4Strin 86 const G4String& theDaughterName2,
87 const G4Strin 87 const G4String& theDaughterName3) :
88 G4VDecayCha 88 G4VDecayChannel("Phase Space",
89 theParentName,theBR, 89 theParentName,theBR,
90 theNumberOfDaughters, 90 theNumberOfDaughters,
91 theDaughterName1, 91 theDaughterName1,
92 theDaughterName2, 92 theDaughterName2,
93 theDaughterName3), 93 theDaughterName3),
94 parentmass(theParentMass), 94 parentmass(theParentMass),
95 theDaughterMasses(0) 95 theDaughterMasses(0)
96 { 96 {
97 if (GetVerboseLevel()>1) G4cout << "G4Genera 97 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
98 } 98 }
99 99
100 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD 100 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
101 101 G4double theParentMass,
102 G4double 102 G4double theBR,
103 G4int 103 G4int theNumberOfDaughters,
104 const G4Strin 104 const G4String& theDaughterName1,
105 const G4Strin 105 const G4String& theDaughterName2,
106 const G4Strin 106 const G4String& theDaughterName3,
107 const G4double *masses) : 107 const G4double *masses) :
108 G4VDecayCha 108 G4VDecayChannel("Phase Space",
109 theParentName,theBR, 109 theParentName,theBR,
110 theNumberOfDaughters, 110 theNumberOfDaughters,
111 theDaughterName1, 111 theDaughterName1,
112 theDaughterName2, 112 theDaughterName2,
113 theDaughterName3), 113 theDaughterName3),
114 parentmass(theParentMass), 114 parentmass(theParentMass),
115 theDaughterMasses(masses) 115 theDaughterMasses(masses)
116 { 116 {
117 if (GetVerboseLevel()>1) G4cout << "G4Genera 117 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
118 } 118 }
119 119
120 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceD 120 G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay(const G4String& theParentName,
121 121 G4double theParentMass,
122 G4double 122 G4double theBR,
123 G4int 123 G4int theNumberOfDaughters,
124 const G4Strin 124 const G4String& theDaughterName1,
125 const G4Strin 125 const G4String& theDaughterName2,
126 const G4Strin 126 const G4String& theDaughterName3,
127 const G4Strin 127 const G4String& theDaughterName4,
128 const G4double *masses) : 128 const G4double *masses) :
129 G4VDecayCha 129 G4VDecayChannel("Phase Space",
130 theParentName,theBR, 130 theParentName,theBR,
131 theNumberOfDaughters, 131 theNumberOfDaughters,
132 theDaughterName1, 132 theDaughterName1,
133 theDaughterName2, 133 theDaughterName2,
134 theDaughterName3, 134 theDaughterName3,
135 theDaughterName4), 135 theDaughterName4),
136 parentmass(theParentMass), 136 parentmass(theParentMass),
137 theDaughterMasses(masses) 137 theDaughterMasses(masses)
138 { 138 {
139 if (GetVerboseLevel()>1) G4cout << "G4Genera 139 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
140 } 140 }
141 141
142 G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpace 142 G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay()
143 { 143 {
144 } 144 }
145 145
146 G4DecayProducts *G4GeneralPhaseSpaceDecay::Dec 146 G4DecayProducts *G4GeneralPhaseSpaceDecay::DecayIt(G4double)
147 { 147 {
148 if (GetVerboseLevel()>1) G4cout << "G4Genera 148 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
149 G4DecayProducts * products = NULL; 149 G4DecayProducts * products = NULL;
150 150
151 CheckAndFillParent(); 151 CheckAndFillParent();
152 CheckAndFillDaughters(); 152 CheckAndFillDaughters();
153 153
154 switch (numberOfDaughters){ 154 switch (numberOfDaughters){
155 case 0: 155 case 0:
156 if (GetVerboseLevel()>0) { 156 if (GetVerboseLevel()>0) {
157 G4cout << "G4GeneralPhaseSpaceDecay::Dec 157 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
158 G4cout << " daughters not defined " <<G4 158 G4cout << " daughters not defined " <<G4endl;
159 } 159 }
160 break; 160 break;
161 case 1: 161 case 1:
162 products = OneBodyDecayIt(); 162 products = OneBodyDecayIt();
163 break; 163 break;
164 case 2: 164 case 2:
165 products = TwoBodyDecayIt(); 165 products = TwoBodyDecayIt();
166 break; 166 break;
167 case 3: 167 case 3:
168 products = ThreeBodyDecayIt(); 168 products = ThreeBodyDecayIt();
169 break; 169 break;
170 default: 170 default:
171 products = ManyBodyDecayIt(); 171 products = ManyBodyDecayIt();
172 break; 172 break;
173 } 173 }
174 if ((products == NULL) && (GetVerboseLevel() 174 if ((products == NULL) && (GetVerboseLevel()>0)) {
175 G4cout << "G4GeneralPhaseSpaceDecay::Decay 175 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
176 G4cout << *parent_name << " can not decay 176 G4cout << *parent_name << " can not decay " << G4endl;
177 DumpInfo(); 177 DumpInfo();
178 } 178 }
179 return products; 179 return products;
180 } 180 }
181 181
182 G4DecayProducts *G4GeneralPhaseSpaceDecay::One 182 G4DecayProducts *G4GeneralPhaseSpaceDecay::OneBodyDecayIt()
183 { 183 {
184 if (GetVerboseLevel()>1) G4cout << "G4Genera 184 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt()"<<G4endl;
185 185
186 // G4double daughtermass = daughters[0]->GetP 186 // G4double daughtermass = daughters[0]->GetPDGMass();
187 187
188 //create parent G4DynamicParticle at rest 188 //create parent G4DynamicParticle at rest
189 G4ParticleMomentum dummy; 189 G4ParticleMomentum dummy;
190 G4DynamicParticle * parentparticle = new G4D 190 G4DynamicParticle * parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
191 191
192 //create G4Decayproducts 192 //create G4Decayproducts
193 G4DecayProducts *products = new G4DecayProdu 193 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
194 delete parentparticle; 194 delete parentparticle;
195 195
196 //create daughter G4DynamicParticle at rest 196 //create daughter G4DynamicParticle at rest
197 G4DynamicParticle * daughterparticle = new G 197 G4DynamicParticle * daughterparticle = new G4DynamicParticle(G4MT_daughters[0], dummy, 0.0);
198 products->PushProducts(daughterparticle); 198 products->PushProducts(daughterparticle);
199 199
200 if (GetVerboseLevel()>1) 200 if (GetVerboseLevel()>1)
201 { 201 {
202 G4cout << "G4GeneralPhaseSpaceDecay::OneB 202 G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt ";
203 G4cout << " create decay products in res 203 G4cout << " create decay products in rest frame " <<G4endl;
204 products->DumpInfo(); 204 products->DumpInfo();
205 } 205 }
206 return products; 206 return products;
207 } 207 }
208 208
209 G4DecayProducts *G4GeneralPhaseSpaceDecay::Two 209 G4DecayProducts *G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()
210 { 210 {
211 if (GetVerboseLevel()>1) G4cout << "G4Genera 211 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()"<<G4endl;
212 212
213 //daughters'mass 213 //daughters'mass
214 G4double daughtermass[2]; 214 G4double daughtermass[2];
215 G4double daughtermomentum; 215 G4double daughtermomentum;
216 if ( theDaughterMasses ) 216 if ( theDaughterMasses )
217 { 217 {
218 daughtermass[0]= *(theDaughterMasses); 218 daughtermass[0]= *(theDaughterMasses);
219 daughtermass[1] = *(theDaughterMasses+1); 219 daughtermass[1] = *(theDaughterMasses+1);
220 } else { 220 } else {
221 daughtermass[0] = G4MT_daughters[0]->GetP 221 daughtermass[0] = G4MT_daughters[0]->GetPDGMass();
222 daughtermass[1] = G4MT_daughters[1]->GetP 222 daughtermass[1] = G4MT_daughters[1]->GetPDGMass();
223 } 223 }
224 224
225 // G4double sumofdaughtermass = daughtermass 225 // G4double sumofdaughtermass = daughtermass[0] + daughtermass[1];
226 226
227 //create parent G4DynamicParticle at rest 227 //create parent G4DynamicParticle at rest
228 G4ParticleMomentum dummy; 228 G4ParticleMomentum dummy;
229 G4DynamicParticle * parentparticle = new G4D 229 G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
230 230
231 //create G4Decayproducts @@GF why dummy par 231 //create G4Decayproducts @@GF why dummy parentparticle?
232 G4DecayProducts *products = new G4DecayProdu 232 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
233 delete parentparticle; 233 delete parentparticle;
234 234
235 //calculate daughter momentum 235 //calculate daughter momentum
236 daughtermomentum = Pmx(parentmass,daughterma 236 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
237 G4double costheta = 2.*G4UniformRand()-1.0; 237 G4double costheta = 2.*G4UniformRand()-1.0;
238 G4double sintheta = std::sqrt((1.0 - costhet 238 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
239 G4double phi = twopi*G4UniformRand()*rad; 239 G4double phi = twopi*G4UniformRand()*rad;
240 G4ParticleMomentum direction(sintheta*std::c 240 G4ParticleMomentum direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
241 241
242 //create daughter G4DynamicParticle 242 //create daughter G4DynamicParticle
243 G4double Etotal= std::sqrt(daughtermass[0]*d 243 G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum);
244 G4DynamicParticle * daughterparticle = new G 244 G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0],Etotal, direction*daughtermomentum);
245 products->PushProducts(daughterparticle); 245 products->PushProducts(daughterparticle);
246 Etotal= std::sqrt(daughtermass[1]*daughterma 246 Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum);
247 daughterparticle = new G4DynamicParticle( G4 247 daughterparticle = new G4DynamicParticle( G4MT_daughters[1],Etotal, direction*(-1.0*daughtermomentum));
248 products->PushProducts(daughterparticle); 248 products->PushProducts(daughterparticle);
249 249
250 if (GetVerboseLevel()>1) 250 if (GetVerboseLevel()>1)
251 { 251 {
252 G4cout << "G4GeneralPhaseSpaceDecay::TwoB 252 G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ";
253 G4cout << " create decay products in res 253 G4cout << " create decay products in rest frame " <<G4endl;
254 products->DumpInfo(); 254 products->DumpInfo();
255 } 255 }
256 return products; 256 return products;
257 } 257 }
258 258
259 G4DecayProducts *G4GeneralPhaseSpaceDecay::Thr 259 G4DecayProducts *G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()
260 // algorism of this code is originally written 260 // algorism of this code is originally written in GDECA3 of GEANT3
261 { 261 {
262 if (GetVerboseLevel()>1) G4cout << "G4Genera 262 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()"<<G4endl;
263 263
264 //daughters'mass 264 //daughters'mass
265 G4double daughtermass[3]; 265 G4double daughtermass[3];
266 G4double sumofdaughtermass = 0.0; 266 G4double sumofdaughtermass = 0.0;
267 for (G4int index=0; index<3; index++) 267 for (G4int index=0; index<3; index++)
268 { 268 {
269 if ( theDaughterMasses ) 269 if ( theDaughterMasses )
270 { 270 {
271 daughtermass[index]= *(theDaughterMas 271 daughtermass[index]= *(theDaughterMasses+index);
272 } else { 272 } else {
273 daughtermass[index] = G4MT_daughters[ 273 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
274 } 274 }
275 sumofdaughtermass += daughtermass[index]; 275 sumofdaughtermass += daughtermass[index];
276 } 276 }
277 277
278 //create parent G4DynamicParticle at rest 278 //create parent G4DynamicParticle at rest
279 G4ParticleMomentum dummy; 279 G4ParticleMomentum dummy;
280 G4DynamicParticle * parentparticle = new G4D 280 G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
281 281
282 //create G4Decayproducts 282 //create G4Decayproducts
283 G4DecayProducts *products = new G4DecayProdu 283 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
284 delete parentparticle; 284 delete parentparticle;
285 285
286 //calculate daughter momentum 286 //calculate daughter momentum
287 // Generate two 287 // Generate two
288 G4double rd1, rd2, rd; 288 G4double rd1, rd2, rd;
289 G4double daughtermomentum[3]; 289 G4double daughtermomentum[3];
290 G4double momentummax=0.0, momentumsum = 0.0; 290 G4double momentummax=0.0, momentumsum = 0.0;
291 G4double energy; 291 G4double energy;
292 const G4int maxNumberOfLoops = 10000; 292 const G4int maxNumberOfLoops = 10000;
293 G4int loopCounter = 0; 293 G4int loopCounter = 0;
294 294
295 do 295 do
296 { 296 {
297 rd1 = G4UniformRand(); 297 rd1 = G4UniformRand();
298 rd2 = G4UniformRand(); 298 rd2 = G4UniformRand();
299 if (rd2 > rd1) 299 if (rd2 > rd1)
300 { 300 {
301 rd = rd1; 301 rd = rd1;
302 rd1 = rd2; 302 rd1 = rd2;
303 rd2 = rd; 303 rd2 = rd;
304 } 304 }
305 momentummax = 0.0; 305 momentummax = 0.0;
306 momentumsum = 0.0; 306 momentumsum = 0.0;
307 // daughter 0 307 // daughter 0
308 308
309 energy = rd2*(parentmass - sumofdaughterm 309 energy = rd2*(parentmass - sumofdaughtermass);
310 daughtermomentum[0] = std::sqrt(energy*en 310 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
311 if ( daughtermomentum[0] >momentummax )mo 311 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
312 momentumsum += daughtermomentum[0]; 312 momentumsum += daughtermomentum[0];
313 313
314 // daughter 1 314 // daughter 1
315 energy = (1.-rd1)*(parentmass - sumofdaug 315 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
316 daughtermomentum[1] = std::sqrt(energy*en 316 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
317 if ( daughtermomentum[1] >momentummax )mo 317 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
318 momentumsum += daughtermomentum[1]; 318 momentumsum += daughtermomentum[1];
319 319
320 // daughter 2 320 // daughter 2
321 energy = (rd1-rd2)*(parentmass - sumofdau 321 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
322 daughtermomentum[2] = std::sqrt(energy*en 322 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
323 if ( daughtermomentum[2] >momentummax )mo 323 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
324 momentumsum += daughtermomentum[2]; 324 momentumsum += daughtermomentum[2];
325 } while ( ( momentummax > momentumsum - m 325 } while ( ( momentummax > momentumsum - momentummax ) && /* Loop checking, 02.11.2015, A.Ribon */
326 ++loopCounter < maxNumberOfLoops 326 ++loopCounter < maxNumberOfLoops );
327 if ( loopCounter >= maxNumberOfLoops ) { 327 if ( loopCounter >= maxNumberOfLoops ) {
328 G4ExceptionDescription ed; 328 G4ExceptionDescription ed;
329 ed << " Failed sampling after maxNumberO 329 ed << " Failed sampling after maxNumberOfLoops attempts : forced exit" << G4endl;
330 G4Exception( " G4GeneralPhaseSpaceDecay: 330 G4Exception( " G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ", "HAD_PHASESPACE_001", FatalException, ed );
331 } 331 }
332 332
333 // output message 333 // output message
334 if (GetVerboseLevel()>1) { 334 if (GetVerboseLevel()>1) {
335 G4cout << " daughter 0:" << daughtermo 335 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
336 G4cout << " daughter 1:" << daughtermo 336 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
337 G4cout << " daughter 2:" << daughtermo 337 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
338 G4cout << " momentum sum:" << momentumsu 338 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
339 } 339 }
340 340
341 //create daughter G4DynamicParticle 341 //create daughter G4DynamicParticle
342 G4double costheta, sintheta, phi, sinphi, co 342 G4double costheta, sintheta, phi, sinphi, cosphi;
343 G4double costhetan, sinthetan, phin, sinphin 343 G4double costhetan, sinthetan, phin, sinphin, cosphin;
344 costheta = 2.*G4UniformRand()-1.0; 344 costheta = 2.*G4UniformRand()-1.0;
345 sintheta = std::sqrt((1.0-costheta)*(1.0+cos 345 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
346 phi = twopi*G4UniformRand()*rad; 346 phi = twopi*G4UniformRand()*rad;
347 sinphi = std::sin(phi); 347 sinphi = std::sin(phi);
348 cosphi = std::cos(phi); 348 cosphi = std::cos(phi);
349 G4ParticleMomentum direction0(sintheta*cosph 349 G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
350 G4double Etotal=std::sqrt( daughtermass[0]*d 350 G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]);
351 G4DynamicParticle * daughterparticle 351 G4DynamicParticle * daughterparticle
352 = new G4DynamicParticle( G4MT_daughte 352 = new G4DynamicParticle( G4MT_daughters[0], Etotal, direction0*daughtermomentum[0]);
353 products->PushProducts(daughterparticle); 353 products->PushProducts(daughterparticle);
354 354
355 costhetan = (daughtermomentum[1]*daughtermom 355 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
356 sinthetan = std::sqrt((1.0-costhetan)*(1.0+c 356 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
357 phin = twopi*G4UniformRand()*rad; 357 phin = twopi*G4UniformRand()*rad;
358 sinphin = std::sin(phin); 358 sinphin = std::sin(phin);
359 cosphin = std::cos(phin); 359 cosphin = std::cos(phin);
360 G4ParticleMomentum direction2; 360 G4ParticleMomentum direction2;
361 direction2.setX( sinthetan*cosphin*costheta* 361 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
362 direction2.setY( sinthetan*cosphin*costheta* 362 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
363 direction2.setZ( -sinthetan*cosphin*sintheta 363 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
364 Etotal=std::sqrt( daughtermass[2]*daughterma 364 Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.mag2());
365 daughterparticle = new G4DynamicParticle( G4 365 daughterparticle = new G4DynamicParticle( G4MT_daughters[2],Etotal, direction2*(daughtermomentum[2]/direction2.mag()));
366 products->PushProducts(daughterparticle); 366 products->PushProducts(daughterparticle);
367 G4ThreeVector mom=(direction0*daughtermoment 367 G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
368 Etotal= std::sqrt( daughtermass[1]*daughterm 368 Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.mag2() );
369 daughterparticle = 369 daughterparticle =
370 new G4DynamicParticle(G4MT_daughters[1] 370 new G4DynamicParticle(G4MT_daughters[1], Etotal, mom);
371 products->PushProducts(daughterparticle); 371 products->PushProducts(daughterparticle);
372 372
373 if (GetVerboseLevel()>1) { 373 if (GetVerboseLevel()>1) {
374 G4cout << "G4GeneralPhaseSpaceDecay::Thre 374 G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ";
375 G4cout << " create decay products in res 375 G4cout << " create decay products in rest frame " <<G4endl;
376 products->DumpInfo(); 376 products->DumpInfo();
377 } 377 }
378 return products; 378 return products;
379 } 379 }
380 380
381 G4DecayProducts *G4GeneralPhaseSpaceDecay::Man 381 G4DecayProducts *G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()
382 // algorism of this code is originally written 382 // algorism of this code is originally written in FORTRAN by M.Asai
383 //******************************************** 383 //*****************************************************************
384 // NBODY 384 // NBODY
385 // N-body phase space Monte-Carlo generator 385 // N-body phase space Monte-Carlo generator
386 // Makoto Asai 386 // Makoto Asai
387 // Hiroshima Institute of Technology 387 // Hiroshima Institute of Technology
388 // (asai@kekvax.kek.jp) 388 // (asai@kekvax.kek.jp)
389 // Revised release : 19/Apr/1995 389 // Revised release : 19/Apr/1995
390 // 390 //
391 { 391 {
392 //return value 392 //return value
393 G4DecayProducts *products; 393 G4DecayProducts *products;
394 394
395 if (GetVerboseLevel()>1) G4cout << "G4Genera 395 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()"<<G4endl;
396 396
397 //daughters'mass 397 //daughters'mass
398 G4double *daughtermass = new G4double[number 398 G4double *daughtermass = new G4double[numberOfDaughters];
399 G4double sumofdaughtermass = 0.0; 399 G4double sumofdaughtermass = 0.0;
400 for (G4int index=0; index<numberOfDaughters; 400 for (G4int index=0; index<numberOfDaughters; index++){
401 daughtermass[index] = G4MT_daughters[index 401 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
402 sumofdaughtermass += daughtermass[index]; 402 sumofdaughtermass += daughtermass[index];
403 } 403 }
404 404
405 //Calculate daughter momentum 405 //Calculate daughter momentum
406 G4double *daughtermomentum = new G4double[nu 406 G4double *daughtermomentum = new G4double[numberOfDaughters];
407 G4ParticleMomentum direction; 407 G4ParticleMomentum direction;
408 G4DynamicParticle **daughterparticle; 408 G4DynamicParticle **daughterparticle;
409 G4double *sm = new G4double[numberOfDaughter 409 G4double *sm = new G4double[numberOfDaughters];
410 G4double tmas; 410 G4double tmas;
411 G4double weight = 1.0; 411 G4double weight = 1.0;
412 G4int numberOfTry = 0; 412 G4int numberOfTry = 0;
413 G4int index1; 413 G4int index1;
414 414
415 do { 415 do {
416 //Generate rundom number in descending ord 416 //Generate rundom number in descending order
417 G4double temp; 417 G4double temp;
418 G4double *rd = new G4double[numberOfDaught 418 G4double *rd = new G4double[numberOfDaughters];
419 rd[0] = 1.0; 419 rd[0] = 1.0;
420 for(index1 =1; index1 < numberOfDaughters 420 for(index1 =1; index1 < numberOfDaughters -1; index1++)
421 rd[index1] = G4UniformRand(); 421 rd[index1] = G4UniformRand();
422 rd[ numberOfDaughters -1] = 0.0; 422 rd[ numberOfDaughters -1] = 0.0;
423 for(index1 =1; index1 < numberOfDaughters 423 for(index1 =1; index1 < numberOfDaughters -1; index1++) {
424 for(G4int index2 = index1+1; index2 < nu 424 for(G4int index2 = index1+1; index2 < numberOfDaughters; index2++) {
425 if (rd[index1] < rd[index2]){ 425 if (rd[index1] < rd[index2]){
426 temp = rd[index1]; 426 temp = rd[index1];
427 rd[index1] = rd[index2]; 427 rd[index1] = rd[index2];
428 rd[index2] = temp; 428 rd[index2] = temp;
429 } 429 }
430 } 430 }
431 } 431 }
432 432
433 //calcurate virtual mass 433 //calcurate virtual mass
434 tmas = parentmass - sumofdaughtermass; 434 tmas = parentmass - sumofdaughtermass;
435 temp = sumofdaughtermass; 435 temp = sumofdaughtermass;
436 for(index1 =0; index1 < numberOfDaughters; 436 for(index1 =0; index1 < numberOfDaughters; index1++) {
437 sm[index1] = rd[index1]*tmas + temp; 437 sm[index1] = rd[index1]*tmas + temp;
438 temp -= daughtermass[index1]; 438 temp -= daughtermass[index1];
439 if (GetVerboseLevel()>1) { 439 if (GetVerboseLevel()>1) {
440 G4cout << index1 << " rundom number:" 440 G4cout << index1 << " rundom number:" << rd[index1];
441 G4cout << " virtual mass:" << sm[ind 441 G4cout << " virtual mass:" << sm[index1]/GeV << "[GeV/c/c]" <<G4endl;
442 } 442 }
443 } 443 }
444 delete [] rd; 444 delete [] rd;
445 445
446 //Calculate daughter momentum 446 //Calculate daughter momentum
447 weight = 1.0; 447 weight = 1.0;
448 index1 =numberOfDaughters-1; 448 index1 =numberOfDaughters-1;
449 daughtermomentum[index1]= Pmx( sm[index1-1 449 daughtermomentum[index1]= Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]);
450 if (GetVerboseLevel()>1) { 450 if (GetVerboseLevel()>1) {
451 G4cout << " daughter " << index1 << 451 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
452 G4cout << " momentum:" << daughtermoment 452 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
453 } 453 }
454 for(index1 =numberOfDaughters-2; index1>=0 454 for(index1 =numberOfDaughters-2; index1>=0; index1--) {
455 // calculate 455 // calculate
456 daughtermomentum[index1]= Pmx( sm[index1 456 daughtermomentum[index1]= Pmx( sm[index1],daughtermass[index1], sm[index1 +1]);
457 if(daughtermomentum[index1] < 0.0) { 457 if(daughtermomentum[index1] < 0.0) {
458 // !!! illegal momentum !!! 458 // !!! illegal momentum !!!
459 if (GetVerboseLevel()>0) { 459 if (GetVerboseLevel()>0) {
460 G4cout << "G4GeneralPhaseSpaceDecay: 460 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
461 G4cout << " can not calculate da 461 G4cout << " can not calculate daughter momentum " <<G4endl;
462 G4cout << " parent:" << *parent_ 462 G4cout << " parent:" << *parent_name;
463 G4cout << " mass:" << parentmass/GeV 463 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
464 G4cout << " daughter " << index1 464 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
465 G4cout << " mass:" << daughtermass[i 465 G4cout << " mass:" << daughtermass[index1]/GeV << "[GeV/c/c]" ;
466 G4cout << " mass:" << daughtermoment 466 G4cout << " mass:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
467 } 467 }
468 delete [] sm; 468 delete [] sm;
469 delete [] daughtermass; 469 delete [] daughtermass;
470 delete [] daughtermomentum; 470 delete [] daughtermomentum;
471 return NULL; // Error detection 471 return NULL; // Error detection
472 472
473 } else { 473 } else {
474 // calculate weight of this events 474 // calculate weight of this events
475 weight *= daughtermomentum[index1]/sm 475 weight *= daughtermomentum[index1]/sm[index1];
476 if (GetVerboseLevel()>1) { 476 if (GetVerboseLevel()>1) {
477 G4cout << " daughter " << index1 477 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
478 G4cout << " momentum:" << daughtermo 478 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
479 } 479 }
480 } 480 }
481 } 481 }
482 if (GetVerboseLevel()>1) { 482 if (GetVerboseLevel()>1) {
483 G4cout << " weight: " << weight <<G4e 483 G4cout << " weight: " << weight <<G4endl;
484 } 484 }
485 485
486 // exit if number of Try exceeds 100 486 // exit if number of Try exceeds 100
487 if (numberOfTry++ >100) { 487 if (numberOfTry++ >100) {
488 if (GetVerboseLevel()>0) { 488 if (GetVerboseLevel()>0) {
489 G4cout << "G4GeneralPhaseSpaceDecay::M 489 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: ";
490 G4cout << " can not determine Decay Kinemati 490 G4cout << " can not determine Decay Kinematics " << G4endl;
491 } 491 }
492 delete [] sm; 492 delete [] sm;
493 delete [] daughtermass; 493 delete [] daughtermass;
494 delete [] daughtermomentum; 494 delete [] daughtermomentum;
495 return NULL; // Error detection 495 return NULL; // Error detection
496 } 496 }
497 } while ( weight > G4UniformRand()); /* Loo 497 } while ( weight > G4UniformRand()); /* Loop checking, 02.11.2015, A.Ribon */
498 if (GetVerboseLevel()>1) { 498 if (GetVerboseLevel()>1) {
499 G4cout << "Start calculation of daughter 499 G4cout << "Start calculation of daughters momentum vector "<<G4endl;
500 } 500 }
501 501
502 G4double costheta, sintheta, phi; 502 G4double costheta, sintheta, phi;
503 G4double beta; 503 G4double beta;
504 daughterparticle = new G4DynamicParticle*[nu 504 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
505 505
506 index1 = numberOfDaughters -2; 506 index1 = numberOfDaughters -2;
507 costheta = 2.*G4UniformRand()-1.0; 507 costheta = 2.*G4UniformRand()-1.0;
508 sintheta = std::sqrt((1.0-costheta)*(1.0+cos 508 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
509 phi = twopi*G4UniformRand()*rad; 509 phi = twopi*G4UniformRand()*rad;
510 direction.setZ(costheta); 510 direction.setZ(costheta);
511 direction.setY(sintheta*std::sin(phi)); 511 direction.setY(sintheta*std::sin(phi));
512 direction.setX(sintheta*std::cos(phi)); 512 direction.setX(sintheta*std::cos(phi));
513 daughterparticle[index1] = new G4DynamicPart 513 daughterparticle[index1] = new G4DynamicParticle( G4MT_daughters[index1], direction*daughtermomentum[index1] );
514 daughterparticle[index1+1] = new G4DynamicPa 514 daughterparticle[index1+1] = new G4DynamicParticle( G4MT_daughters[index1+1], direction*(-1.0*daughtermomentum[index1]) );
515 515
516 for (index1 = numberOfDaughters -3; index1 516 for (index1 = numberOfDaughters -3; index1 >= 0; index1--) {
517 //calculate momentum direction 517 //calculate momentum direction
518 costheta = 2.*G4UniformRand()-1.0; 518 costheta = 2.*G4UniformRand()-1.0;
519 sintheta = std::sqrt((1.0-costheta)*(1.0+c 519 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
520 phi = twopi*G4UniformRand()*rad; 520 phi = twopi*G4UniformRand()*rad;
521 direction.setZ(costheta); 521 direction.setZ(costheta);
522 direction.setY(sintheta*std::sin(phi)); 522 direction.setY(sintheta*std::sin(phi));
523 direction.setX(sintheta*std::cos(phi)); 523 direction.setX(sintheta*std::cos(phi));
524 524
525 // boost already created particles 525 // boost already created particles
526 beta = daughtermomentum[index1]; 526 beta = daughtermomentum[index1];
527 beta /= std::sqrt( daughtermomentum[index1 527 beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] );
528 for (G4int index2 = index1+1; index2<numbe 528 for (G4int index2 = index1+1; index2<numberOfDaughters; index2++) {
529 G4LorentzVector p4; 529 G4LorentzVector p4;
530 // make G4LorentzVector for secondaries 530 // make G4LorentzVector for secondaries
531 p4 = daughterparticle[index2]->Get4Momen 531 p4 = daughterparticle[index2]->Get4Momentum();
532 532
533 // boost secondaries to new frame 533 // boost secondaries to new frame
534 p4.boost( direction.x()*beta, direction. 534 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
535 535
536 // change energy/momentum 536 // change energy/momentum
537 daughterparticle[index2]->Set4Momentum(p 537 daughterparticle[index2]->Set4Momentum(p4);
538 } 538 }
539 //create daughter G4DynamicParticle 539 //create daughter G4DynamicParticle
540 daughterparticle[index1]= new G4DynamicPar 540 daughterparticle[index1]= new G4DynamicParticle( G4MT_daughters[index1], direction*(-1.0*daughtermomentum[index1]));
541 } 541 }
542 542
543 //create G4Decayproducts 543 //create G4Decayproducts
544 G4DynamicParticle *parentparticle; 544 G4DynamicParticle *parentparticle;
545 direction.setX(1.0); direction.setY(0.0); d 545 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
546 parentparticle = new G4DynamicParticle( G4MT 546 parentparticle = new G4DynamicParticle( G4MT_parent, direction, 0.0);
547 products = new G4DecayProducts(*parentpartic 547 products = new G4DecayProducts(*parentparticle);
548 delete parentparticle; 548 delete parentparticle;
549 for (index1 = 0; index1<numberOfDaughters; i 549 for (index1 = 0; index1<numberOfDaughters; index1++) {
550 products->PushProducts(daughterparticle[in 550 products->PushProducts(daughterparticle[index1]);
551 } 551 }
552 if (GetVerboseLevel()>1) { 552 if (GetVerboseLevel()>1) {
553 G4cout << "G4GeneralPhaseSpaceDecay::ManyB 553 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
554 G4cout << " create decay products in rest 554 G4cout << " create decay products in rest frame " << G4endl;
555 products->DumpInfo(); 555 products->DumpInfo();
556 } 556 }
557 557
558 delete [] daughterparticle; 558 delete [] daughterparticle;
559 delete [] daughtermomentum; 559 delete [] daughtermomentum;
560 delete [] daughtermass; 560 delete [] daughtermass;
561 delete [] sm; 561 delete [] sm;
562 562
563 return products; 563 return products;
564 } 564 }
565 565
566 566
567 567
568 568
569 569
570 570