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