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