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