Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VDecayChannel.cc
Go to the documentation of this file.
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 //
27 // $Id: G4VDecayChannel.cc 98352 2016-07-08 08:21:00Z gcosmo $
28 //
29 //
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 // ------------------------------------------------------------
38 
39 #include "G4ParticleDefinition.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ParticleTable.hh"
42 #include "G4DecayTable.hh"
43 #include "G4DecayProducts.hh"
44 #include "G4VDecayChannel.hh"
45 #include "G4AutoLock.hh"
46 
48 
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;
64 
65  // set pointer to G4ParticleTable (static and singleton object)
67 }
68 
70  :kinematics_name(aName),
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;
86 
87  // set pointer to G4ParticleTable (static and singleton object)
89 }
90 
92  const G4String& theParentName,
93  G4double theBR,
94  G4int theNumberOfDaughters,
95  const G4String& theDaughterName1,
96  const G4String& theDaughterName2,
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;
115 
116  // set pointer to G4ParticleTable (static and singleton object)
118 
119  // parent name
120  parent_name = new G4String(theParentName);
121 
122  // cleate array
124  for (G4int index=0;index<numberOfDaughters;index++) daughters_name[index]=0;
125 
126  // daughters' name
127  if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
128  if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
129  if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
130  if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
131 }
132 
134 {
136  verboseLevel = right.verboseLevel;
137  rbranch = right.rbranch;
138  rangeMass = right.rangeMass;
139 
140  // copy parent name
141  parent_name = new G4String(*right.parent_name);
142  G4MT_parent = 0;
143  G4MT_parent_mass = 0.0;
144 
145  //create array
147 
148  daughters_name =0;
149  if ( numberOfDaughters >0 ) {
151  //copy daughters name
152  for (G4int index=0; index < numberOfDaughters; index++)
153  {
154  daughters_name[index] = new G4String(*right.daughters_name[index]);
155  }
156  }
157 
158  //
160  G4MT_daughters = 0;
162 
163  // particle table
165 
167 
170 }
171 
173 {
174  if (this != &right) {
176  verboseLevel = right.verboseLevel;
177  rbranch = right.rbranch;
178  rangeMass = right.rangeMass;
180  // copy parent name
181  parent_name = new G4String(*right.parent_name);
182 
183  // clear daughters_name array
185 
186  // recreate array
188  if ( numberOfDaughters >0 ) {
191  //copy daughters name
192  for (G4int index=0; index < numberOfDaughters; index++) {
193  daughters_name[index] = new G4String(*right.daughters_name[index]);
194  }
195  }
196  }
197 
198  //
199  G4MT_parent = 0;
200  G4MT_daughters = 0;
201  G4MT_parent_mass = 0.0;
204 
205  // particle table
207 
210 
211  return *this;
212 }
213 
215 {
217  if (parent_name != 0) delete parent_name;
218  parent_name = 0;
219  if (G4MT_daughters_mass != 0) delete [] G4MT_daughters_mass;
221  if (G4MT_daughters_width != 0) delete [] G4MT_daughters_width;
225 }
226 
228 {
230  if ( daughters_name != 0) {
231  if (numberOfDaughters>0) {
232 #ifdef G4VERBOSE
233  if (verboseLevel>1) {
234  G4cout << "G4VDecayChannel::ClearDaughtersName "
235  << " for " << *parent_name << G4endl;
236  }
237 #endif
238  for (G4int index=0; index < numberOfDaughters; index++) {
239  if (daughters_name[index] != 0) delete daughters_name[index];
240  }
241  }
242  delete [] daughters_name;
243  daughters_name = 0;
244  }
245  //
246  if (G4MT_daughters != 0) delete [] G4MT_daughters;
247  if (G4MT_daughters_mass != 0) delete [] G4MT_daughters_mass;
248  if (G4MT_daughters_width != 0) delete [] G4MT_daughters_width;
250  G4MT_daughters = 0;
253 
254  numberOfDaughters = 0;
255 }
256 
258 {
259  if (size >0) {
260  // remove old contents
262  // cleate array
263  daughters_name = new G4String*[size];
264  for (G4int index=0;index<size;index++) daughters_name[index]=0;
265  numberOfDaughters = size;
266  }
267 }
268 
270  const G4String &particle_name)
271 {
272  // check numberOfDaughters is positive
273  if (numberOfDaughters<=0) {
274 #ifdef G4VERBOSE
275  if (verboseLevel>0) {
276  G4cout << "G4VDecayChannel::SetDaughter: "
277  << "Number of daughters is not defined" << G4endl;
278  }
279 #endif
280  return;
281  }
282 
283  //ANDREA:-> Feb 25 2016
284  // An analysis of this code, shows that this method is called
285  // only in the constructor of derived classes.
286  // The general idea of this method is probably to support
287  // the possibility to re-define daughters on the fly, however
288  // this design is extremely problematic for MT mode, we thus
289  // require (as practically happens) that the method is called only
290  // at construction, i.e. when G4MT_daugheters == 0
291  // moreover this method can be called only after SetNumberOfDaugthers
292  // has been called (see previous if), in such a case daughters_name != 0
293  if ( daughters_name == 0 ) {
294  G4Exception("G4VDecayChannel::SetDaughter","PART112",FatalException,
295  "Trying to add a daughter without specifying number of secondaries, useSetNumberOfDaughters first");
296  return;
297  }
298  if ( G4MT_daughters != 0 ) {
299  G4Exception("G4VDecayChannel::SetDaughter","PART111",FatalException,
300  "Trying to modify a daughter of a decay channel, but decay channel already has daughters.");
301  return;
302  }
303  //<-:ANDREA
304 
305  // check an index
306  if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
307 #ifdef G4VERBOSE
308  if (verboseLevel>0) {
309  G4cout << "G4VDecayChannel::SetDaughter"
310  << "index out of range " << anIndex << G4endl;
311  }
312 #endif
313  } else {
314  // fill the name
315  daughters_name[anIndex] = new G4String(particle_name);
316 #ifdef G4VERBOSE
317  if (verboseLevel>1) {
318  G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
319  G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
320  }
321 #endif
322  }
323 }
324 
325 void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
326 {
327  if (parent_type != 0) SetDaughter(anIndex, parent_type->GetParticleName());
328 }
329 
330 void G4VDecayChannel::FillDaughters()
331 {
332  G4AutoLock lock(&daughtersMutex);
333  //Double check, check again if another thread has already filled this, in
334  //case do not need to do anything
335  if ( G4MT_daughters != 0 ) return;
336 
337  G4int index;
338 
339 #ifdef G4VERBOSE
340  if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
341 #endif
342  if (G4MT_daughters != 0) {
343  delete [] G4MT_daughters;
344  G4MT_daughters = 0;
345  }
346 
347  // parent mass
349  G4double parentmass = G4MT_parent->GetPDGMass();
350 
351  //
352  G4double sumofdaughtermass = 0.0;
353  G4double sumofdaughterwidthsq = 0.0;
354 
355  if ((numberOfDaughters <=0) || (daughters_name == 0) ){
356 #ifdef G4VERBOSE
357  if (verboseLevel>0) {
358  G4cout << "G4VDecayChannel::FillDaughters "
359  << "[ " << G4MT_parent->GetParticleName() << " ]"
360  << "numberOfDaughters is not defined yet";
361  }
362 #endif
363  G4MT_daughters = 0;
364  G4Exception("G4VDecayChannel::FillDaughters",
365  "PART011", FatalException,
366  "Can not fill daughters: numberOfDaughters is not defined yet");
367  }
368 
369  //create and set the array of pointers to daughter particles
371  if (G4MT_daughters_mass != 0) delete [] G4MT_daughters_mass;
372  if (G4MT_daughters_width != 0) delete [] G4MT_daughters_width;
375  // loop over all daughters
376  for (index=0; index < numberOfDaughters; index++) {
377  if (daughters_name[index] == 0) {
378  // daughter name is not defined
379 #ifdef G4VERBOSE
380  if (verboseLevel>0) {
381  G4cout << "G4VDecayChannel::FillDaughters "
382  << "[ " << G4MT_parent->GetParticleName() << " ]"
383  << index << "-th daughter is not defined yet" << G4endl;
384  }
385 #endif
386  G4MT_daughters[index] = 0;
387  G4Exception("G4VDecayChannel::FillDaughters",
388  "PART011", FatalException,
389  "Can not fill daughters: name of a daughter is not defined yet");
390  }
391  //search daughter particles in the particle table
393  if (G4MT_daughters[index] == 0) {
394  // can not find the daughter particle
395 #ifdef G4VERBOSE
396  if (verboseLevel>0) {
397  G4cout << "G4VDecayChannel::FillDaughters "
398  << "[ " << G4MT_parent->GetParticleName() << " ]"
399  << index << ":" << *daughters_name[index]
400  << " is not defined !!" << G4endl;
401  G4cout << " The BR of this decay mode is set to zero " << G4endl;
402  }
403 #endif
404  SetBR(0.0);
405  return;
406  }
407 #ifdef G4VERBOSE
408  if (verboseLevel>1) {
409  G4cout << index << ":" << *daughters_name[index];
410  G4cout << ":" << G4MT_daughters[index] << G4endl;
411  }
412 #endif
413  G4MT_daughters_mass[index] = G4MT_daughters[index]->GetPDGMass();
414  G4double d_width = G4MT_daughters[index]->GetPDGWidth();
415  G4MT_daughters_width[index] = d_width;
416  sumofdaughtermass += G4MT_daughters[index]->GetPDGMass();
417  sumofdaughterwidthsq += d_width*d_width;
418  } // end loop over all daughters
419 
420  // check sum of daghter mass
421  G4double widthMass = std::sqrt(G4MT_parent->GetPDGWidth()*G4MT_parent->GetPDGWidth()+sumofdaughterwidthsq);
422  if ( (G4MT_parent->GetParticleType() != "nucleus") &&
423  (numberOfDaughters !=1) &&
424  (sumofdaughtermass > parentmass + rangeMass*widthMass) ){
425  // !!! illegal mass !!!
426 #ifdef G4VERBOSE
427  if (GetVerboseLevel()>0) {
428  G4cout << "G4VDecayChannel::FillDaughters "
429  << "[ " << G4MT_parent->GetParticleName() << " ]"
430  << " Energy/Momentum conserevation breaks " <<G4endl;
431  if (GetVerboseLevel()>1) {
432  G4cout << " parent:" << *parent_name
433  << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
434  for (index=0; index < numberOfDaughters; index++){
435  G4cout << " daughter " << index << ":" << *daughters_name[index]
436  << " mass:" << G4MT_daughters[index]->GetPDGMass()/GeV
437  << "[GeV/c/c]" <<G4endl;
438  }
439  }
440  }
441 #endif
442  }
443 }
444 
445 
446 void G4VDecayChannel::FillParent()
447 {
448  G4AutoLock lock(&parentMutex);
449  //Double check, check again if another thread has already filled this, in
450  //case do not need to do anything
451  if ( G4MT_parent != 0 ) return;
452 
453  if (parent_name == 0) {
454  // parent name is not defined
455 #ifdef G4VERBOSE
456  if (verboseLevel>0) {
457  G4cout << "G4VDecayChannel::FillParent "
458  << ": parent name is not defined !!" << G4endl;
459  }
460 #endif
461  G4MT_parent = 0;
462  G4Exception("G4VDecayChannel::FillParent()",
463  "PART012", FatalException,
464  "Can not fill parent: parent name is not defined yet");
465  return;
466  }
467  // search parent particle in the particle table
469  if (G4MT_parent == 0) {
470  // parent particle does not exist
471 #ifdef G4VERBOSE
472  if (verboseLevel>0) {
473  G4cout << "G4VDecayChannel::FillParent "
474  << *parent_name << " does not exist !!" << G4endl;
475  }
476 #endif
477  G4Exception("G4VDecayChannel::FillParent()",
478  "PART012", FatalException,
479  "Can not fill parent: parent does not exist");
480  return;
481  }
483 }
484 
486 {
487  if (parent_type != 0) SetParent(parent_type->GetParticleName());
488 }
489 
491 {
492  // determine angular momentum
493 
494  // fill pointers to daughter particles if not yet set
496 
497  const G4int PiSpin = G4MT_parent->GetPDGiSpin();
498  const G4int PParity = G4MT_parent->GetPDGiParity();
499  if (2==numberOfDaughters) { // up to now we can only handle two particle decays
500  const G4int D1iSpin = G4MT_daughters[0]->GetPDGiSpin();
501  const G4int D1Parity = G4MT_daughters[0]->GetPDGiParity();
502  const G4int D2iSpin = G4MT_daughters[1]->GetPDGiSpin();
503  const G4int D2Parity = G4MT_daughters[1]->GetPDGiParity();
504  const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
505  const G4int MaxiSpin = D1iSpin + D2iSpin;
506  const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int
507  G4int lMin;
508 #ifdef G4VERBOSE
509  if (verboseLevel>1) {
510  G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
511  G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
512  }
513 #endif
514  for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings
515  lMin = std::abs(PiSpin-j)/2;
516 #ifdef G4VERBOSE
517  if (verboseLevel>1)
518  G4cout << "-> checking 2*j=" << j << G4endl;
519 #endif
520  for (G4int l=lMin; l<=lMax; l++) {
521 #ifdef G4VERBOSE
522  if (verboseLevel>1)
523  G4cout << " checking l=" << l << G4endl;
524 #endif
525  if (l%2==0) {
526  if (PParity == D1Parity*D2Parity) { // check parity for this l
527  return l;
528  }
529  } else {
530  if (PParity == -1*D1Parity*D2Parity) { // check parity for this l
531  return l;
532  }
533  }
534  }
535  }
536  } else {
537  G4Exception("G4VDecayChannel::GetAngularMomentum",
538  "PART111", JustWarning,
539  "Sorry, can't handle 3 particle decays (up to now)");
540  return 0;
541  }
542  G4Exception ("G4VDecayChannel::GetAngularMomentum",
543  "PART111", JustWarning,
544  "Can't find angular momentum for this decay");
545  return 0;
546 }
547 
549 {
550  G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
551  G4cout << " : " ;
552  for (G4int index=0; index < numberOfDaughters; index++){
553  if(daughters_name[index] != 0) {
554  G4cout << " " << *(daughters_name[index]);
555  } else {
556  G4cout << " not defined ";
557  }
558  }
559  G4cout << G4endl;
560 }
561 
562 const G4String& G4VDecayChannel::GetNoName() const
563 {
564  return noName;
565 }
566 
567 #include "Randomize.hh"
569 {
570  if (width<=0.0) return massPDG;
571  if (maxDev >rangeMass) maxDev = rangeMass;
572  if (maxDev <=-1.*rangeMass) return massPDG; // can not calculate
573 
574  G4double x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
575  G4double y = G4UniformRand();
576  const size_t MAX_LOOP=10000;
577  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
578  if ( y * (width*width*x*x + massPDG*massPDG*width*width) <= massPDG*massPDG*width*width ) break;
579  x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
580  y = G4UniformRand();
581  }
582  G4double mass = massPDG + x*width;
583  return mass;
584 }
585 
587 {
588  G4double sumOfDaughterMassMin=0.0;
591  // skip one body decay
592  if (numberOfDaughters==1) return true;
593 
594  for (G4int index=0; index < numberOfDaughters; index++) {
595  sumOfDaughterMassMin +=
597  }
598  return (parentMass >= sumOfDaughterMassMin);
599 }
void CheckAndFillDaughters()
static const G4String noName
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
void SetBR(G4double value)
G4double * G4MT_daughters_mass
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
#define G4MUTEXINIT(mutex)
Definition: G4Threading.hh:177
#define width
G4double G4MT_parent_mass
virtual G4bool IsOKWithParentMass(G4double parentMass)
tuple x
Definition: test.py:50
G4ParticleTable * particletable
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
const G4String & GetParticleName() const
G4int GetAngularMomentum()
G4double GetPDGWidth() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4String kinematics_name
void SetNumberOfDaughters(G4int value)
G4String * parent_name
const G4String & GetParticleType() const
G4VDecayChannel & operator=(const G4VDecayChannel &)
G4ThreeVector parent_polarization
virtual ~G4VDecayChannel()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4int GetVerboseLevel() const
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4double * G4MT_daughters_width
void SetParent(const G4ParticleDefinition *particle_type)
void SetDaughter(G4int anIndex, const G4ParticleDefinition *particle_type)
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
G4double DynamicalMass(G4double massPDG, G4double width, G4double maxDev=+1.) const
double G4double
Definition: G4Types.hh:76
#define G4MUTEXDESTROY(mutex)
Definition: G4Threading.hh:178
G4String ** daughters_name