| 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 // G4VDecayChannel 26 // G4VDecayChannel
27 // 27 //
28 // Author: H.Kurashige, 27 July 1996 << 28 // Author: H.Kurashige, 27 July 1996
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4VDecayChannel.hh" <<
32 <<
33 #include "G4AutoLock.hh" <<
34 #include "G4DecayProducts.hh" <<
35 #include "G4DecayTable.hh" <<
36 #include "G4ParticleDefinition.hh" 31 #include "G4ParticleDefinition.hh"
37 #include "G4ParticleTable.hh" <<
38 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
>> 33 #include "G4ParticleTable.hh"
>> 34 #include "G4DecayTable.hh"
>> 35 #include "G4DecayProducts.hh"
>> 36 #include "G4VDecayChannel.hh"
>> 37 #include "G4AutoLock.hh"
39 #include "Randomize.hh" 38 #include "Randomize.hh"
40 39
41 const G4String G4VDecayChannel::noName = " "; 40 const G4String G4VDecayChannel::noName = " ";
42 41
>> 42 // --------------------------------------------------------------------
43 G4VDecayChannel::G4VDecayChannel() 43 G4VDecayChannel::G4VDecayChannel()
>> 44 : parent_polarization()
44 { 45 {
45 // set pointer to G4ParticleTable (static an 46 // set pointer to G4ParticleTable (static and singleton object)
46 particletable = G4ParticleTable::GetParticle 47 particletable = G4ParticleTable::GetParticleTable();
47 } 48 }
48 49
>> 50 // --------------------------------------------------------------------
49 G4VDecayChannel::G4VDecayChannel(const G4Strin 51 G4VDecayChannel::G4VDecayChannel(const G4String& aName, G4int verbose)
50 : kinematics_name(aName), verboseLevel(verbo << 52 : kinematics_name(aName), parent_polarization(), verboseLevel(verbose)
51 { 53 {
52 // set pointer to G4ParticleTable (static an 54 // set pointer to G4ParticleTable (static and singleton object)
53 particletable = G4ParticleTable::GetParticle 55 particletable = G4ParticleTable::GetParticleTable();
54 } 56 }
55 57
56 G4VDecayChannel::G4VDecayChannel(const G4Strin << 58 // --------------------------------------------------------------------
57 G4double theB << 59 G4VDecayChannel::G4VDecayChannel(const G4String& aName,
58 const G4Strin << 60 const G4String& theParentName,
59 const G4Strin << 61 G4double theBR,
60 const G4Strin << 62 G4int theNumberOfDaughters,
61 : kinematics_name(aName), rbranch(theBR), nu << 63 const G4String& theDaughterName1,
>> 64 const G4String& theDaughterName2,
>> 65 const G4String& theDaughterName3,
>> 66 const G4String& theDaughterName4,
>> 67 const G4String& theDaughterName5 )
>> 68 : kinematics_name(aName), rbranch(theBR), parent_polarization(),
>> 69 numberOfDaughters(theNumberOfDaughters)
62 { 70 {
63 // set pointer to G4ParticleTable (static an 71 // set pointer to G4ParticleTable (static and singleton object)
64 particletable = G4ParticleTable::GetParticle 72 particletable = G4ParticleTable::GetParticleTable();
65 73
66 // parent name 74 // parent name
67 parent_name = new G4String(theParentName); 75 parent_name = new G4String(theParentName);
68 76
69 // cleate array 77 // cleate array
70 daughters_name = new G4String*[numberOfDaugh 78 daughters_name = new G4String*[numberOfDaughters];
71 for (G4int index = 0; index < numberOfDaught << 79 for (G4int index=0; index<numberOfDaughters; ++index)
>> 80 {
72 daughters_name[index] = nullptr; 81 daughters_name[index] = nullptr;
73 } 82 }
74 83
75 // daughters' name 84 // daughters' name
76 if (numberOfDaughters > 0) daughters_name[0] << 85 if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
77 if (numberOfDaughters > 1) daughters_name[1] << 86 if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
78 if (numberOfDaughters > 2) daughters_name[2] << 87 if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
79 if (numberOfDaughters > 3) daughters_name[3] << 88 if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
80 if (numberOfDaughters > 4) daughters_name[4] << 89 if (numberOfDaughters>4) daughters_name[4] = new G4String(theDaughterName5);
81 << 90
82 if (rbranch < 0.) << 91 if (rbranch <0. ) rbranch = 0.0;
83 rbranch = 0.0; << 92 else if (rbranch >1.0 ) rbranch = 1.0;
84 else if (rbranch > 1.0) <<
85 rbranch = 1.0; <<
86 } 93 }
87 94
>> 95 // --------------------------------------------------------------------
88 G4VDecayChannel::G4VDecayChannel(const G4VDeca 96 G4VDecayChannel::G4VDecayChannel(const G4VDecayChannel& right)
89 { 97 {
90 kinematics_name = right.kinematics_name; 98 kinematics_name = right.kinematics_name;
91 verboseLevel = right.verboseLevel; 99 verboseLevel = right.verboseLevel;
92 rbranch = right.rbranch; 100 rbranch = right.rbranch;
93 rangeMass = right.rangeMass; 101 rangeMass = right.rangeMass;
94 102
95 // copy parent name 103 // copy parent name
96 parent_name = new G4String(*right.parent_nam 104 parent_name = new G4String(*right.parent_name);
97 105
98 // create array 106 // create array
99 numberOfDaughters = right.numberOfDaughters; 107 numberOfDaughters = right.numberOfDaughters;
100 108
101 daughters_name = nullptr; 109 daughters_name = nullptr;
102 if (numberOfDaughters > 0) { << 110 if ( numberOfDaughters >0 )
>> 111 {
103 daughters_name = new G4String*[numberOfDau 112 daughters_name = new G4String*[numberOfDaughters];
104 // copy daughters name 113 // copy daughters name
105 for (G4int index = 0; index < numberOfDaug << 114 for (G4int index=0; index<numberOfDaughters; ++index)
>> 115 {
106 daughters_name[index] = new G4String(*ri 116 daughters_name[index] = new G4String(*right.daughters_name[index]);
107 } 117 }
108 } 118 }
109 119
110 // particle table 120 // particle table
111 particletable = G4ParticleTable::GetParticle 121 particletable = G4ParticleTable::GetParticleTable();
112 122
113 parent_polarization = right.parent_polarizat 123 parent_polarization = right.parent_polarization;
114 124
115 G4MUTEXINIT(daughtersMutex); 125 G4MUTEXINIT(daughtersMutex);
116 G4MUTEXINIT(parentMutex); 126 G4MUTEXINIT(parentMutex);
117 } 127 }
118 128
>> 129 // --------------------------------------------------------------------
119 G4VDecayChannel& G4VDecayChannel::operator=(co 130 G4VDecayChannel& G4VDecayChannel::operator=(const G4VDecayChannel& right)
120 { 131 {
121 if (this != &right) { << 132 if (this != &right)
>> 133 {
122 kinematics_name = right.kinematics_name; 134 kinematics_name = right.kinematics_name;
123 verboseLevel = right.verboseLevel; 135 verboseLevel = right.verboseLevel;
124 rbranch = right.rbranch; 136 rbranch = right.rbranch;
125 rangeMass = right.rangeMass; << 137 rangeMass = right.rangeMass;
126 parent_polarization = right.parent_polariz 138 parent_polarization = right.parent_polarization;
127 // copy parent name 139 // copy parent name
128 delete parent_name; 140 delete parent_name;
129 parent_name = new G4String(*right.parent_n 141 parent_name = new G4String(*right.parent_name);
130 142
131 // clear daughters_name array 143 // clear daughters_name array
132 ClearDaughtersName(); 144 ClearDaughtersName();
133 145
134 // recreate array 146 // recreate array
135 numberOfDaughters = right.numberOfDaughter 147 numberOfDaughters = right.numberOfDaughters;
136 if (numberOfDaughters > 0) { << 148 if ( numberOfDaughters >0 )
>> 149 {
137 daughters_name = new G4String*[numberOfD 150 daughters_name = new G4String*[numberOfDaughters];
138 // copy daughters name 151 // copy daughters name
139 for (G4int index = 0; index < numberOfDa << 152 for (G4int index=0; index<numberOfDaughters; ++index)
>> 153 {
140 daughters_name[index] = new G4String(* 154 daughters_name[index] = new G4String(*right.daughters_name[index]);
141 } 155 }
142 } 156 }
143 } 157 }
144 158
145 // particle table 159 // particle table
146 particletable = G4ParticleTable::GetParticle 160 particletable = G4ParticleTable::GetParticleTable();
147 161
148 G4MUTEXINIT(daughtersMutex); 162 G4MUTEXINIT(daughtersMutex);
149 G4MUTEXINIT(parentMutex); 163 G4MUTEXINIT(parentMutex);
150 164
151 return *this; 165 return *this;
152 } 166 }
153 167
>> 168 // --------------------------------------------------------------------
154 G4VDecayChannel::~G4VDecayChannel() 169 G4VDecayChannel::~G4VDecayChannel()
155 { 170 {
156 ClearDaughtersName(); 171 ClearDaughtersName();
157 delete parent_name; 172 delete parent_name;
158 parent_name = nullptr; 173 parent_name = nullptr;
159 delete[] G4MT_daughters_mass; << 174 delete [] G4MT_daughters_mass;
160 G4MT_daughters_mass = nullptr; 175 G4MT_daughters_mass = nullptr;
161 delete[] G4MT_daughters_width; << 176 delete [] G4MT_daughters_width;
162 G4MT_daughters_width = nullptr; 177 G4MT_daughters_width = nullptr;
163 G4MUTEXDESTROY(daughtersMutex); 178 G4MUTEXDESTROY(daughtersMutex);
164 G4MUTEXDESTROY(parentMutex); 179 G4MUTEXDESTROY(parentMutex);
165 } << 180 }
166 181
>> 182 // --------------------------------------------------------------------
167 void G4VDecayChannel::ClearDaughtersName() 183 void G4VDecayChannel::ClearDaughtersName()
168 { 184 {
169 G4AutoLock l(&daughtersMutex); 185 G4AutoLock l(&daughtersMutex);
170 if (daughters_name != nullptr) { << 186 if ( daughters_name != nullptr)
171 if (numberOfDaughters > 0) { << 187 {
>> 188 if (numberOfDaughters>0)
>> 189 {
172 #ifdef G4VERBOSE 190 #ifdef G4VERBOSE
173 if (verboseLevel > 1) { << 191 if (verboseLevel>1)
>> 192 {
174 G4cout << "G4VDecayChannel::ClearDaugh 193 G4cout << "G4VDecayChannel::ClearDaughtersName() "
175 << " for " << *parent_name << G 194 << " for " << *parent_name << G4endl;
176 } 195 }
177 #endif 196 #endif
178 for (G4int index = 0; index < numberOfDa << 197 for (G4int index=0; index<numberOfDaughters; ++index)
>> 198 {
179 delete daughters_name[index]; 199 delete daughters_name[index];
180 } 200 }
181 } 201 }
182 delete[] daughters_name; << 202 delete [] daughters_name;
183 daughters_name = nullptr; 203 daughters_name = nullptr;
184 } 204 }
185 205
186 delete[] G4MT_daughters; << 206 delete [] G4MT_daughters;
187 delete[] G4MT_daughters_mass; << 207 delete [] G4MT_daughters_mass;
188 delete[] G4MT_daughters_width; << 208 delete [] G4MT_daughters_width;
189 G4MT_daughters_width = nullptr; 209 G4MT_daughters_width = nullptr;
190 G4MT_daughters = nullptr; 210 G4MT_daughters = nullptr;
191 G4MT_daughters_mass = nullptr; 211 G4MT_daughters_mass = nullptr;
192 212
193 numberOfDaughters = 0; 213 numberOfDaughters = 0;
194 } 214 }
195 215
>> 216 // --------------------------------------------------------------------
196 void G4VDecayChannel::SetNumberOfDaughters(G4i 217 void G4VDecayChannel::SetNumberOfDaughters(G4int size)
197 { 218 {
198 if (size > 0) { << 219 if (size >0)
>> 220 {
199 // remove old contents 221 // remove old contents
200 ClearDaughtersName(); 222 ClearDaughtersName();
201 // cleate array 223 // cleate array
202 daughters_name = new G4String*[size]; 224 daughters_name = new G4String*[size];
203 for (G4int index = 0; index < size; ++inde << 225 for (G4int index=0; index<size; ++index)
>> 226 {
204 daughters_name[index] = nullptr; 227 daughters_name[index] = nullptr;
205 } 228 }
206 numberOfDaughters = size; 229 numberOfDaughters = size;
207 } 230 }
208 } 231 }
209 232
>> 233 // --------------------------------------------------------------------
210 void G4VDecayChannel::SetDaughter(G4int anInde 234 void G4VDecayChannel::SetDaughter(G4int anIndex, const G4String& particle_name)
211 { 235 {
212 // check numberOfDaughters is positive 236 // check numberOfDaughters is positive
213 if (numberOfDaughters <= 0) { << 237 if (numberOfDaughters<=0)
>> 238 {
214 #ifdef G4VERBOSE 239 #ifdef G4VERBOSE
215 if (verboseLevel > 0) { << 240 if (verboseLevel>0)
>> 241 {
216 G4cout << "G4VDecayChannel::SetDaughter( 242 G4cout << "G4VDecayChannel::SetDaughter() - "
217 << "Number of daughters is not de 243 << "Number of daughters is not defined" << G4endl;
218 } 244 }
219 #endif 245 #endif
220 return; 246 return;
221 } 247 }
222 248
223 // An analysis of this code, shows that this 249 // An analysis of this code, shows that this method is called
224 // only in the constructor of derived classe 250 // only in the constructor of derived classes.
225 // The general idea of this method is probab 251 // The general idea of this method is probably to support
226 // the possibility to re-define daughters on 252 // the possibility to re-define daughters on the fly, however
227 // this design is extremely problematic for 253 // this design is extremely problematic for MT mode, we thus
228 // require (as practically happens) that the 254 // require (as practically happens) that the method is called only
229 // at construction, i.e. when G4MT_daughters 255 // at construction, i.e. when G4MT_daughters == 0
230 // moreover this method can be called only a 256 // moreover this method can be called only after SetNumberOfDaugthers()
231 // has been called (see previous if), in suc 257 // has been called (see previous if), in such a case daughters_name != nullptr
232 // 258 //
233 if (daughters_name == nullptr) { << 259 if ( daughters_name == nullptr )
>> 260 {
234 G4Exception("G4VDecayChannel::SetDaughter( 261 G4Exception("G4VDecayChannel::SetDaughter()", "PART112", FatalException,
235 "Trying to add a daughter with << 262 "Trying to add a daughter without specifying number of secondaries!");
236 return; 263 return;
237 } 264 }
238 if (G4MT_daughters != nullptr) { << 265 if ( G4MT_daughters != nullptr )
>> 266 {
239 G4Exception("G4VDecayChannel::SetDaughter( 267 G4Exception("G4VDecayChannel::SetDaughter()", "PART111", FatalException,
240 "Trying to modify a daughter o 268 "Trying to modify a daughter of a decay channel, \
241 but decay channel already has 269 but decay channel already has daughters.");
242 return; 270 return;
243 } 271 }
244 272
245 // check an index << 273 // check an index
246 if ((anIndex < 0) || (anIndex >= numberOfDau << 274 if ( (anIndex<0) || (anIndex>=numberOfDaughters) )
>> 275 {
247 #ifdef G4VERBOSE 276 #ifdef G4VERBOSE
248 if (verboseLevel > 0) { << 277 if (verboseLevel>0)
>> 278 {
249 G4cout << "G4VDecayChannel::SetDaughter( 279 G4cout << "G4VDecayChannel::SetDaughter() - "
250 << "index out of range " << anInd 280 << "index out of range " << anIndex << G4endl;
251 } 281 }
252 #endif 282 #endif
253 } 283 }
254 else { << 284 else
>> 285 {
255 // fill the name 286 // fill the name
256 daughters_name[anIndex] = new G4String(par 287 daughters_name[anIndex] = new G4String(particle_name);
257 #ifdef G4VERBOSE 288 #ifdef G4VERBOSE
258 if (verboseLevel > 1) { << 289 if (verboseLevel>1)
259 G4cout << "G4VDecayChannel::SetDaughter[ << 290 {
260 G4cout << daughters_name[anIndex] << ":" << 291 G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
>> 292 G4cout << daughters_name[anIndex] << ":"
>> 293 << *daughters_name[anIndex]<<G4endl;
261 } 294 }
262 #endif 295 #endif
263 } 296 }
264 } 297 }
265 298
266 void G4VDecayChannel::SetDaughter(G4int anInde << 299 // --------------------------------------------------------------------
>> 300 void G4VDecayChannel::
>> 301 SetDaughter(G4int anIndex, const G4ParticleDefinition* parent_type)
267 { 302 {
268 if (parent_type != nullptr) SetDaughter(anIn << 303 if (parent_type != nullptr)
>> 304 SetDaughter(anIndex, parent_type->GetParticleName());
269 } 305 }
270 306
>> 307 // --------------------------------------------------------------------
271 void G4VDecayChannel::FillDaughters() 308 void G4VDecayChannel::FillDaughters()
272 { 309 {
273 G4AutoLock lock(&daughtersMutex); 310 G4AutoLock lock(&daughtersMutex);
274 311
275 // Double check, check again if another thre 312 // Double check, check again if another thread has already filled this, in
276 // case do not need to do anything 313 // case do not need to do anything
277 if (G4MT_daughters != nullptr) return; << 314 if ( G4MT_daughters != nullptr ) return;
278 315
279 G4int index; 316 G4int index;
280 << 317
281 #ifdef G4VERBOSE 318 #ifdef G4VERBOSE
282 if (verboseLevel > 1) G4cout << "G4VDecayCha << 319 if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" << G4endl;
283 #endif 320 #endif
284 if (G4MT_daughters != nullptr) { << 321 if (G4MT_daughters != nullptr)
285 delete[] G4MT_daughters; << 322 {
>> 323 delete [] G4MT_daughters;
286 G4MT_daughters = nullptr; 324 G4MT_daughters = nullptr;
287 } 325 }
288 326
289 // parent mass 327 // parent mass
290 CheckAndFillParent(); 328 CheckAndFillParent();
291 G4double parentmass = G4MT_parent->GetPDGMas 329 G4double parentmass = G4MT_parent->GetPDGMass();
292 330
293 // 331 //
294 G4double sumofdaughtermass = 0.0; 332 G4double sumofdaughtermass = 0.0;
295 G4double sumofdaughterwidthsq = 0.0; 333 G4double sumofdaughterwidthsq = 0.0;
296 334
297 if ((numberOfDaughters <= 0) || (daughters_n << 335 if ((numberOfDaughters <=0) || (daughters_name == nullptr) )
>> 336 {
298 #ifdef G4VERBOSE 337 #ifdef G4VERBOSE
299 if (verboseLevel > 0) { << 338 if (verboseLevel>0)
>> 339 {
300 G4cout << "G4VDecayChannel::FillDaughter 340 G4cout << "G4VDecayChannel::FillDaughters() - "
301 << "[ " << G4MT_parent->GetPartic 341 << "[ " << G4MT_parent->GetParticleName() << " ]"
302 << "numberOfDaughters is not defi 342 << "numberOfDaughters is not defined yet";
303 } 343 }
304 #endif 344 #endif
305 G4MT_daughters = nullptr; 345 G4MT_daughters = nullptr;
306 G4Exception("G4VDecayChannel::FillDaughter << 346 G4Exception("G4VDecayChannel::FillDaughters()",
>> 347 "PART011", FatalException,
307 "Cannot fill daughters: number 348 "Cannot fill daughters: numberOfDaughters is not defined yet");
308 } << 349 }
309 350
310 // create and set the array of pointers to d 351 // create and set the array of pointers to daughter particles
311 G4MT_daughters = new G4ParticleDefinition*[n 352 G4MT_daughters = new G4ParticleDefinition*[numberOfDaughters];
312 delete[] G4MT_daughters_mass; << 353 if (G4MT_daughters_mass != nullptr) delete [] G4MT_daughters_mass;
313 delete[] G4MT_daughters_width; << 354 if (G4MT_daughters_width != nullptr) delete [] G4MT_daughters_width;
314 G4MT_daughters_mass = new G4double[numberOfD 355 G4MT_daughters_mass = new G4double[numberOfDaughters];
315 G4MT_daughters_width = new G4double[numberOf 356 G4MT_daughters_width = new G4double[numberOfDaughters];
316 // loop over all daughters 357 // loop over all daughters
317 for (index = 0; index < numberOfDaughters; + << 358 for (index=0; index<numberOfDaughters; ++index)
318 if (daughters_name[index] == nullptr) { << 359 {
>> 360 if (daughters_name[index] == nullptr)
>> 361 {
319 // daughter name is not defined 362 // daughter name is not defined
320 #ifdef G4VERBOSE 363 #ifdef G4VERBOSE
321 if (verboseLevel > 0) { << 364 if (verboseLevel>0)
>> 365 {
322 G4cout << "G4VDecayChannel::FillDaught 366 G4cout << "G4VDecayChannel::FillDaughters() - "
323 << "[ " << G4MT_parent->GetPart << 367 << "[ " << G4MT_parent->GetParticleName() << " ]"
324 << "-th daughter is not defined << 368 << index << "-th daughter is not defined yet" << G4endl;
325 } 369 }
326 #endif 370 #endif
327 G4MT_daughters[index] = nullptr; 371 G4MT_daughters[index] = nullptr;
328 G4Exception("G4VDecayChannel::FillDaught << 372 G4Exception("G4VDecayChannel::FillDaughters()",
>> 373 "PART011", FatalException,
329 "Cannot fill daughters: name 374 "Cannot fill daughters: name of daughter is not defined yet");
330 } << 375 }
331 // search daughter particles in the partic << 376 // search daughter particles in the particle table
332 G4MT_daughters[index] = particletable->Fin 377 G4MT_daughters[index] = particletable->FindParticle(*daughters_name[index]);
333 if (G4MT_daughters[index] == nullptr) { << 378 if (G4MT_daughters[index] == nullptr )
>> 379 {
334 // cannot find the daughter particle 380 // cannot find the daughter particle
335 #ifdef G4VERBOSE 381 #ifdef G4VERBOSE
336 if (verboseLevel > 0) { << 382 if (verboseLevel>0)
>> 383 {
337 G4cout << "G4VDecayChannel::FillDaught 384 G4cout << "G4VDecayChannel::FillDaughters() - "
338 << "[ " << G4MT_parent->GetPart << 385 << "[ " << G4MT_parent->GetParticleName() << " ]"
339 << *daughters_name[index] << " << 386 << index << ":" << *daughters_name[index]
>> 387 << " is not defined !!" << G4endl;
340 G4cout << " The BR of this decay mode 388 G4cout << " The BR of this decay mode is set to zero." << G4endl;
341 } 389 }
342 #endif 390 #endif
343 SetBR(0.0); 391 SetBR(0.0);
344 return; 392 return;
345 } 393 }
346 #ifdef G4VERBOSE 394 #ifdef G4VERBOSE
347 if (verboseLevel > 1) { << 395 if (verboseLevel>1)
>> 396 {
348 G4cout << index << ":" << *daughters_nam 397 G4cout << index << ":" << *daughters_name[index];
349 G4cout << ":" << G4MT_daughters[index] < 398 G4cout << ":" << G4MT_daughters[index] << G4endl;
350 } 399 }
351 #endif 400 #endif
352 G4MT_daughters_mass[index] = G4MT_daughter 401 G4MT_daughters_mass[index] = G4MT_daughters[index]->GetPDGMass();
353 G4double d_width = G4MT_daughters[index]-> 402 G4double d_width = G4MT_daughters[index]->GetPDGWidth();
354 G4MT_daughters_width[index] = d_width; 403 G4MT_daughters_width[index] = d_width;
355 sumofdaughtermass += G4MT_daughters[index] 404 sumofdaughtermass += G4MT_daughters[index]->GetPDGMass();
356 sumofdaughterwidthsq += d_width * d_width; << 405 sumofdaughterwidthsq += d_width*d_width;
357 } // end loop over all daughters 406 } // end loop over all daughters
358 407
359 // check sum of daghter mass 408 // check sum of daghter mass
360 G4double widthMass = << 409 G4double widthMass = std::sqrt(G4MT_parent->GetPDGWidth()
361 std::sqrt(G4MT_parent->GetPDGWidth() * G4M << 410 *G4MT_parent->GetPDGWidth()
362 if ((G4MT_parent->GetParticleType() != "nucl << 411 +sumofdaughterwidthsq);
363 && (sumofdaughtermass > parentmass + ran << 412 if ( (G4MT_parent->GetParticleType() != "nucleus")
>> 413 && (numberOfDaughters !=1)
>> 414 && (sumofdaughtermass > parentmass + rangeMass*widthMass) )
364 { 415 {
365 // !!! illegal mass !!! 416 // !!! illegal mass !!!
366 #ifdef G4VERBOSE 417 #ifdef G4VERBOSE
367 if (GetVerboseLevel() > 0) { << 418 if (GetVerboseLevel()>0)
>> 419 {
368 G4cout << "G4VDecayChannel::FillDaughter 420 G4cout << "G4VDecayChannel::FillDaughters() - "
369 << "[ " << G4MT_parent->GetPartic 421 << "[ " << G4MT_parent->GetParticleName() << " ]"
370 << " Energy/Momentum consereva 422 << " Energy/Momentum conserevation breaks " << G4endl;
371 if (GetVerboseLevel() > 1) { << 423 if (GetVerboseLevel()>1)
372 G4cout << " parent:" << *parent_nam << 424 {
373 << G4endl; << 425 G4cout << " parent:" << *parent_name
374 for (index = 0; index < numberOfDaught << 426 << " mass:" << parentmass/GeV << "[GeV/c/c]" << G4endl;
375 G4cout << " daughter " << index << 427 for (index=0; index<numberOfDaughters; ++index)
376 << " mass:" << G4MT_daughters << 428 {
377 } << 429 G4cout << " daughter " << index << ":" << *daughters_name[index]
378 } << 430 << " mass:" << G4MT_daughters[index]->GetPDGMass()/GeV
379 } << 431 << "[GeV/c/c]" << G4endl;
>> 432 }
>> 433 }
>> 434 }
380 #endif 435 #endif
381 } << 436 }
382 } 437 }
383 438
>> 439 // --------------------------------------------------------------------
384 void G4VDecayChannel::FillParent() 440 void G4VDecayChannel::FillParent()
385 { 441 {
386 G4AutoLock lock(&parentMutex); 442 G4AutoLock lock(&parentMutex);
387 443
388 // Double check, check again if another thre 444 // Double check, check again if another thread has already filled this, in
389 // case do not need to do anything 445 // case do not need to do anything
390 if (G4MT_parent != nullptr) return; << 446 if ( G4MT_parent != nullptr ) return;
391 447
392 if (parent_name == nullptr) { << 448 if (parent_name == nullptr)
>> 449 {
393 // parent name is not defined 450 // parent name is not defined
394 #ifdef G4VERBOSE 451 #ifdef G4VERBOSE
395 if (verboseLevel > 0) { << 452 if (verboseLevel>0)
>> 453 {
396 G4cout << "G4VDecayChannel::FillParent() 454 G4cout << "G4VDecayChannel::FillParent() - "
397 << "parent name is not defined !! 455 << "parent name is not defined !!" << G4endl;
398 } 456 }
399 #endif 457 #endif
400 G4MT_parent = nullptr; 458 G4MT_parent = nullptr;
401 G4Exception("G4VDecayChannel::FillParent() << 459 G4Exception("G4VDecayChannel::FillParent()",
402 "Cannot fill parent: parent na << 460 "PART012", FatalException,
>> 461 "Cannot fill parent: parent name is not defined yet");
403 return; 462 return;
404 } 463 }
405 464
406 // search parent particle in the particle ta 465 // search parent particle in the particle table
407 G4MT_parent = particletable->FindParticle(*p 466 G4MT_parent = particletable->FindParticle(*parent_name);
408 if (G4MT_parent == nullptr) { << 467 if (G4MT_parent == nullptr)
>> 468 {
409 // parent particle does not exist 469 // parent particle does not exist
410 #ifdef G4VERBOSE 470 #ifdef G4VERBOSE
411 if (verboseLevel > 0) { << 471 if (verboseLevel>0)
412 G4cout << "G4VDecayChannel::FillParent() << 472 {
413 << G4endl; << 473 G4cout << "G4VDecayChannel::FillParent() - "
>> 474 << *parent_name << " does not exist !!" << G4endl;
414 } 475 }
415 #endif 476 #endif
416 G4Exception("G4VDecayChannel::FillParent() << 477 G4Exception("G4VDecayChannel::FillParent()",
>> 478 "PART012", FatalException,
417 "Cannot fill parent: parent do 479 "Cannot fill parent: parent does not exist");
418 return; 480 return;
419 } 481 }
420 G4MT_parent_mass = G4MT_parent->GetPDGMass() 482 G4MT_parent_mass = G4MT_parent->GetPDGMass();
421 } 483 }
422 484
>> 485 // --------------------------------------------------------------------
423 void G4VDecayChannel::SetParent(const G4Partic 486 void G4VDecayChannel::SetParent(const G4ParticleDefinition* parent_type)
424 { 487 {
425 if (parent_type != nullptr) SetParent(parent 488 if (parent_type != nullptr) SetParent(parent_type->GetParticleName());
426 } 489 }
427 490
>> 491 // --------------------------------------------------------------------
428 G4int G4VDecayChannel::GetAngularMomentum() 492 G4int G4VDecayChannel::GetAngularMomentum()
429 { 493 {
430 // determine angular momentum 494 // determine angular momentum
431 495
432 // fill pointers to daughter particles if no << 496 // fill pointers to daughter particles if not yet set
433 CheckAndFillDaughters(); 497 CheckAndFillDaughters();
434 498
435 const G4int PiSpin = G4MT_parent->GetPDGiSpi 499 const G4int PiSpin = G4MT_parent->GetPDGiSpin();
436 const G4int PParity = G4MT_parent->GetPDGiPa 500 const G4int PParity = G4MT_parent->GetPDGiParity();
437 if (2 == numberOfDaughters) // up to now we << 501 if (2==numberOfDaughters) // up to now we can only handle two particle decays
438 { 502 {
439 const G4int D1iSpin = G4MT_daughters[0]->G << 503 const G4int D1iSpin = G4MT_daughters[0]->GetPDGiSpin();
440 const G4int D1Parity = G4MT_daughters[0]-> 504 const G4int D1Parity = G4MT_daughters[0]->GetPDGiParity();
441 const G4int D2iSpin = G4MT_daughters[1]->G << 505 const G4int D2iSpin = G4MT_daughters[1]->GetPDGiSpin();
442 const G4int D2Parity = G4MT_daughters[1]-> 506 const G4int D2Parity = G4MT_daughters[1]->GetPDGiParity();
443 const G4int MiniSpin = std::abs(D1iSpin - << 507 const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
444 const G4int MaxiSpin = D1iSpin + D2iSpin; 508 const G4int MaxiSpin = D1iSpin + D2iSpin;
445 const G4int lMax = (PiSpin + D1iSpin + D2i << 509 const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is always int
446 G4int lMin; 510 G4int lMin;
447 #ifdef G4VERBOSE 511 #ifdef G4VERBOSE
448 if (verboseLevel > 1) { << 512 if (verboseLevel>1)
449 G4cout << "iSpin: " << PiSpin << " -> " << 513 {
450 G4cout << "2*jmin, 2*jmax, lmax " << Min << 514 G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin
451 } << 515 << " + " << D2iSpin << G4endl;
452 #endif << 516 G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin
453 for (G4int j = MiniSpin; j <= MaxiSpin; j << 517 << " " << lMax << G4endl;
454 { // spin couplings << 518 }
455 lMin = std::abs(PiSpin - j) / 2; << 519 #endif
456 #ifdef G4VERBOSE << 520 for (G4int j=MiniSpin; j<=MaxiSpin; j+=2) // loop over all possible
457 if (verboseLevel > 1) G4cout << "-> chec << 521 { // spin couplings
458 #endif << 522 lMin = std::abs(PiSpin-j)/2;
459 for (G4int l = lMin; l <= lMax; ++l) { << 523 #ifdef G4VERBOSE
460 #ifdef G4VERBOSE << 524 if (verboseLevel>1)
461 if (verboseLevel > 1) G4cout << " chec << 525 G4cout << "-> checking 2*j=" << j << G4endl;
462 #endif << 526 #endif
463 if (l % 2 == 0) { << 527 for (G4int l=lMin; l<=lMax; ++l)
464 if (PParity == D1Parity * D2Parity) << 528 {
>> 529 #ifdef G4VERBOSE
>> 530 if (verboseLevel>1)
>> 531 G4cout << " checking l=" << l << G4endl;
>> 532 #endif
>> 533 if (l%2==0)
>> 534 {
>> 535 if (PParity == D1Parity*D2Parity) // check parity for this l
465 return l; 536 return l;
466 } 537 }
467 else { << 538 else
468 if (PParity == -1 * D1Parity * D2Par << 539 {
>> 540 if (PParity == -1*D1Parity*D2Parity) // check parity for this l
469 return l; 541 return l;
470 } 542 }
471 } 543 }
472 } 544 }
473 } 545 }
474 else { << 546 else
475 G4Exception("G4VDecayChannel::GetAngularMo << 547 {
>> 548 G4Exception("G4VDecayChannel::GetAngularMomentum()",
>> 549 "PART111", JustWarning,
476 "Sorry, can't handle 3 particl 550 "Sorry, can't handle 3 particle decays (up to now)");
477 return 0; 551 return 0;
478 } 552 }
479 G4Exception("G4VDecayChannel::GetAngularMome << 553 G4Exception ("G4VDecayChannel::GetAngularMomentum()",
480 "Can't find angular momentum for << 554 "PART111", JustWarning,
>> 555 "Can't find angular momentum for this decay");
481 return 0; 556 return 0;
482 } 557 }
483 558
>> 559 // --------------------------------------------------------------------
484 void G4VDecayChannel::DumpInfo() 560 void G4VDecayChannel::DumpInfo()
485 { 561 {
486 G4cout << " BR: " << rbranch << " [" << ki 562 G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
487 G4cout << " : "; << 563 G4cout << " : " ;
488 for (G4int index = 0; index < numberOfDaught << 564 for (G4int index=0; index<numberOfDaughters; ++index)
489 if (daughters_name[index] != nullptr) { << 565 {
>> 566 if(daughters_name[index] != nullptr)
>> 567 {
490 G4cout << " " << *(daughters_name[index] 568 G4cout << " " << *(daughters_name[index]);
491 } 569 }
492 else { << 570 else
>> 571 {
493 G4cout << " not defined "; 572 G4cout << " not defined ";
494 } 573 }
495 } 574 }
496 G4cout << G4endl; 575 G4cout << G4endl;
497 } 576 }
498 577
>> 578 // --------------------------------------------------------------------
499 const G4String& G4VDecayChannel::GetNoName() c 579 const G4String& G4VDecayChannel::GetNoName() const
500 { 580 {
501 return noName; 581 return noName;
502 } 582 }
503 583
504 G4double G4VDecayChannel::DynamicalMass(G4doub << 584 // --------------------------------------------------------------------
505 { << 585 G4double G4VDecayChannel::
506 if (width <= 0.0) return massPDG; << 586 DynamicalMass(G4double massPDG, G4double width, G4double maxDev ) const
507 if (maxDev > rangeMass) maxDev = rangeMass; << 587 {
508 if (maxDev <= -1. * rangeMass) return massPD << 588 if (width<=0.0) return massPDG;
509 << 589 if (maxDev >rangeMass) maxDev = rangeMass;
510 G4double x = G4UniformRand() * (maxDev + ran << 590 if (maxDev <=-1.*rangeMass) return massPDG; // cannot calculate
>> 591
>> 592 G4double x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
511 G4double y = G4UniformRand(); 593 G4double y = G4UniformRand();
512 const std::size_t MAX_LOOP = 10000; << 594 const std::size_t MAX_LOOP=10000;
513 for (std::size_t loop_counter = 0; loop_coun << 595 for (std::size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter)
514 if (y * (width * width * x * x + massPDG * << 596 {
515 <= massPDG * massPDG * width * width) << 597 if ( y * (width*width*x*x + massPDG*massPDG*width*width)
516 break; << 598 <= massPDG*massPDG*width*width ) break;
517 x = G4UniformRand() * (maxDev + rangeMass) << 599 x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
518 y = G4UniformRand(); 600 y = G4UniformRand();
519 } 601 }
520 G4double mass = massPDG + x * width; << 602 G4double mass = massPDG + x*width;
521 return mass; 603 return mass;
522 } 604 }
523 << 605
>> 606 // --------------------------------------------------------------------
524 G4bool G4VDecayChannel::IsOKWithParentMass(G4d 607 G4bool G4VDecayChannel::IsOKWithParentMass(G4double parentMass)
525 { 608 {
526 G4double sumOfDaughterMassMin = 0.0; 609 G4double sumOfDaughterMassMin = 0.0;
527 CheckAndFillParent(); 610 CheckAndFillParent();
528 CheckAndFillDaughters(); 611 CheckAndFillDaughters();
529 // skip one body decay 612 // skip one body decay
530 if (numberOfDaughters == 1) return true; << 613 if (numberOfDaughters==1) return true;
531 << 614
532 for (G4int index = 0; index < numberOfDaught << 615 for (G4int index=0; index<numberOfDaughters; ++index)
533 sumOfDaughterMassMin += G4MT_daughters_mas << 616 {
>> 617 sumOfDaughterMassMin +=
>> 618 G4MT_daughters_mass[index] - rangeMass*G4MT_daughters_width[index];
534 } 619 }
535 return (parentMass >= sumOfDaughterMassMin); << 620 return (parentMass >= sumOfDaughterMassMin);
536 } 621 }
537 622
538 void G4VDecayChannel::SetBR(G4double value) << 623 // --------------------------------------------------------------------
539 { << 624 void G4VDecayChannel::SetBR(G4double value)
540 rbranch = value; << 625 {
541 if (rbranch < 0.) << 626 rbranch = value;
542 rbranch = 0.0; << 627 if (rbranch <0. ) rbranch = 0.0;
543 else if (rbranch > 1.0) << 628 else if (rbranch >1.0 ) rbranch = 1.0;
544 rbranch = 1.0; <<
545 } 629 }
546 630