Geant4  10.01.p03
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 82270 2014-06-13 13:51:56Z 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 
61 
63 
65  :kinematics_name(""),
66  rbranch(0.0),
67  numberOfDaughters(0),
68  parent_name(0), daughters_name(0),
69  rangeMass(1.0),
70  particletable(0),
71  verboseLevel(1)
72 {
74  G4MT_parent = 0;
75  G4MT_daughters = 0;
76  G4MT_parent_mass = 0.0;
79 
80  // set pointer to G4ParticleTable (static and singleton object)
82 }
83 
85  :kinematics_name(aName),
86  rbranch(0.0),
87  numberOfDaughters(0),
88  parent_name(0), daughters_name(0),
89  rangeMass(1.0),
90  particletable(0),
91  verboseLevel(Verbose)
92 {
94  G4MT_parent = 0;
95  G4MT_daughters = 0;
96  G4MT_parent_mass = 0.0;
99 
100  // set pointer to G4ParticleTable (static and singleton object)
102 }
103 
105  const G4String& theParentName,
106  G4double theBR,
107  G4int theNumberOfDaughters,
108  const G4String& theDaughterName1,
109  const G4String& theDaughterName2,
110  const G4String& theDaughterName3,
111  const G4String& theDaughterName4 )
112  :kinematics_name(aName),
113  rbranch(theBR),
114  numberOfDaughters(theNumberOfDaughters),
115  parent_name(0), daughters_name(0),
116  rangeMass(1.0),
117  particletable(0),
118  verboseLevel(1)
119 {
121  G4MT_parent = 0;
122  G4MT_daughters = 0;
123  G4MT_parent_mass = 0.0;
126 
127  // set pointer to G4ParticleTable (static and singleton object)
129 
130  // parent name
131  parent_name = new G4String(theParentName);
132 
133  // cleate array
135  for (G4int index=0;index<numberOfDaughters;index++) daughters_name[index]=0;
136 
137  // daughters' name
138  if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
139  if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
140  if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
141  if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
142 }
143 
145 {
147 
149  verboseLevel = right.verboseLevel;
150  rbranch = right.rbranch;
151  rangeMass = right.rangeMass;
152 
153  // copy parent name
154  parent_name = new G4String(*right.parent_name);
155  G4MT_parent = 0;
156  G4MT_parent_mass = 0.0;
157 
158  //create array
160 
161  daughters_name =0;
162  if ( numberOfDaughters >0 ) {
164  //copy daughters name
165  for (G4int index=0; index < numberOfDaughters; index++)
166  {
167  daughters_name[index] = new G4String(*right.daughters_name[index]);
168  }
169  }
170 
171  //
173  G4MT_daughters = 0;
175 
176  // particle table
178 }
179 
181 {
182  if (this != &right) {
184  verboseLevel = right.verboseLevel;
185  rbranch = right.rbranch;
186  rangeMass = right.rangeMass;
187 
188  // copy parent name
189  parent_name = new G4String(*right.parent_name);
190 
191  // clear daughters_name array
193 
194  // recreate array
196  if ( numberOfDaughters >0 ) {
199  //copy daughters name
200  for (G4int index=0; index < numberOfDaughters; index++) {
201  daughters_name[index] = new G4String(*right.daughters_name[index]);
202  }
203  }
204  }
205 
206  //
207  G4MT_parent = 0;
208  G4MT_daughters = 0;
209  G4MT_parent_mass = 0.0;
212 
213  // particle table
215 
216  return *this;
217 }
218 
220 {
222  if (parent_name != 0) delete parent_name;
223  parent_name = 0;
224  if (G4MT_daughters_mass != 0) delete [] G4MT_daughters_mass;
226  if (G4MT_daughters_width != 0) delete [] G4MT_daughters_width;
228 }
229 
234 
236 {
237  if ( daughters_name != 0) {
238  if (numberOfDaughters>0) {
239 #ifdef G4VERBOSE
240  if (verboseLevel>1) {
241  G4cout << "G4VDecayChannel::ClearDaughtersName "
242  << " for " << *parent_name << G4endl;
243  }
244 #endif
245  for (G4int index=0; index < numberOfDaughters; index++) {
246  if (daughters_name[index] != 0) delete daughters_name[index];
247  }
248  }
249  delete [] daughters_name;
250  daughters_name = 0;
251  }
252  //
253  if (G4MT_daughters != 0) delete [] G4MT_daughters;
254  if (G4MT_daughters_mass != 0) delete [] G4MT_daughters_mass;
255  if (G4MT_daughters_width != 0) delete [] G4MT_daughters_width;
257  G4MT_daughters = 0;
260 
261  numberOfDaughters = 0;
262 }
263 
265 {
266  if (size >0) {
267  // remove old contents
269  // cleate array
270  daughters_name = new G4String*[size];
271  for (G4int index=0;index<size;index++) daughters_name[index]=0;
272  numberOfDaughters = size;
273  }
274 }
275 
277  const G4String &particle_name)
278 {
279  // check numberOfDaughters is positive
280  if (numberOfDaughters<=0) {
281 #ifdef G4VERBOSE
282  if (verboseLevel>0) {
283  G4cout << "G4VDecayChannel::SetDaughter: "
284  << "Number of daughters is not defined" << G4endl;
285  }
286 #endif
287  return;
288  }
289 
290  // check existence of daughters_name array
291  if (daughters_name == 0) {
292  // cleate array
294  for (G4int index=0;index<numberOfDaughters;index++) {
295  daughters_name[index]=0;
296  }
297  }
298 
299  // check an index
300  if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
301 #ifdef G4VERBOSE
302  if (verboseLevel>0) {
303  G4cout << "G4VDecayChannel::SetDaughter"
304  << "index out of range " << anIndex << G4endl;
305  }
306 #endif
307  } else {
308  // delete the old name if it exists
309  if (daughters_name[anIndex]!=0) delete daughters_name[anIndex];
310  // fill the name
311  daughters_name[anIndex] = new G4String(particle_name);
312  // refill the array of daughters[] if it exists
313  if (G4MT_daughters != 0) FillDaughters();
314 #ifdef G4VERBOSE
315  if (verboseLevel>1) {
316  G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
317  G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
318  }
319 #endif
320  }
321 }
322 
323 void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
324 {
325  if (parent_type != 0) SetDaughter(anIndex, parent_type->GetParticleName());
326 }
327 
329 {
330  G4int index;
331 
332 #ifdef G4VERBOSE
333  if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
334 #endif
335  if (G4MT_daughters != 0) {
336  delete [] G4MT_daughters;
337  G4MT_daughters = 0;
338  }
339 
340  // parent mass
341  if (G4MT_parent == 0) FillParent();
342  G4double parentmass = G4MT_parent->GetPDGMass();
343 
344  //
345  G4double sumofdaughtermass = 0.0;
346  G4double sumofdaughterwidthsq = 0.0;
347 
348  if ((numberOfDaughters <=0) || (daughters_name == 0) ){
349 #ifdef G4VERBOSE
350  if (verboseLevel>0) {
351  G4cout << "G4VDecayChannel::FillDaughters "
352  << "[ " << G4MT_parent->GetParticleName() << " ]"
353  << "numberOfDaughters is not defined yet";
354  }
355 #endif
356  G4MT_daughters = 0;
357  G4Exception("G4VDecayChannel::FillDaughters",
358  "PART011", FatalException,
359  "Can not fill daughters: numberOfDaughters is not defined yet");
360  }
361 
362  //create and set the array of pointers to daughter particles
364  if (G4MT_daughters_mass != 0) delete [] G4MT_daughters_mass;
365  if (G4MT_daughters_width != 0) delete [] G4MT_daughters_width;
368  // loop over all daughters
369  for (index=0; index < numberOfDaughters; index++) {
370  if (daughters_name[index] == 0) {
371  // daughter name is not defined
372 #ifdef G4VERBOSE
373  if (verboseLevel>0) {
374  G4cout << "G4VDecayChannel::FillDaughters "
375  << "[ " << G4MT_parent->GetParticleName() << " ]"
376  << index << "-th daughter is not defined yet" << G4endl;
377  }
378 #endif
379  G4MT_daughters[index] = 0;
380  G4Exception("G4VDecayChannel::FillDaughters",
381  "PART011", FatalException,
382  "Can not fill daughters: name of a daughter is not defined yet");
383  }
384  //search daughter particles in the particle table
386  if (G4MT_daughters[index] == 0) {
387  // can not find the daughter particle
388 #ifdef G4VERBOSE
389  if (verboseLevel>0) {
390  G4cout << "G4VDecayChannel::FillDaughters "
391  << "[ " << G4MT_parent->GetParticleName() << " ]"
392  << index << ":" << *daughters_name[index]
393  << " is not defined !!" << G4endl;
394  G4cout << " The BR of this decay mode is set to zero " << G4endl;
395  }
396 #endif
397  SetBR(0.0);
398  return;
399  }
400 #ifdef G4VERBOSE
401  if (verboseLevel>1) {
402  G4cout << index << ":" << *daughters_name[index];
403  G4cout << ":" << G4MT_daughters[index] << G4endl;
404  }
405 #endif
406  G4MT_daughters_mass[index] = G4MT_daughters[index]->GetPDGMass();
407  G4double d_width = G4MT_daughters[index]->GetPDGWidth();
408  G4MT_daughters_width[index] = d_width;
409  sumofdaughtermass += G4MT_daughters[index]->GetPDGMass();
410  sumofdaughterwidthsq += d_width*d_width;
411  } // end loop over all daughters
412 
413  // check sum of daghter mass
414  G4double widthMass = std::sqrt(G4MT_parent->GetPDGWidth()*G4MT_parent->GetPDGWidth()+sumofdaughterwidthsq);
415  if ( (G4MT_parent->GetParticleType() != "nucleus") &&
416  (sumofdaughtermass > parentmass + rangeMass*widthMass) ){
417  // !!! illegal mass !!!
418 #ifdef G4VERBOSE
419  if (GetVerboseLevel()>0) {
420  G4cout << "G4VDecayChannel::FillDaughters "
421  << "[ " << G4MT_parent->GetParticleName() << " ]"
422  << " Energy/Momentum conserevation breaks " <<G4endl;
423  if (GetVerboseLevel()>1) {
424  G4cout << " parent:" << *parent_name
425  << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
426  for (index=0; index < numberOfDaughters; index++){
427  G4cout << " daughter " << index << ":" << *daughters_name[index]
428  << " mass:" << G4MT_daughters[index]->GetPDGMass()/GeV
429  << "[GeV/c/c]" <<G4endl;
430  }
431  }
432  }
433 #endif
434  }
435 }
436 
437 
439 {
440  if (parent_name == 0) {
441  // parent name is not defined
442 #ifdef G4VERBOSE
443  if (verboseLevel>0) {
444  G4cout << "G4VDecayChannel::FillParent "
445  << ": parent name is not defined !!" << G4endl;
446  }
447 #endif
448  G4MT_parent = 0;
449  G4Exception("G4VDecayChannel::FillParent()",
450  "PART012", FatalException,
451  "Can not fill parent: parent name is not defined yet");
452  return;
453  }
454  // search parent particle in the particle table
456  if (G4MT_parent == 0) {
457  // parent particle does not exist
458 #ifdef G4VERBOSE
459  if (verboseLevel>0) {
460  G4cout << "G4VDecayChannel::FillParent "
461  << *parent_name << " does not exist !!" << G4endl;
462  }
463 #endif
464  G4Exception("G4VDecayChannel::FillParent()",
465  "PART012", FatalException,
466  "Can not fill parent: parent does not exist");
467  return;
468  }
470 }
471 
473 {
474  if (parent_type != 0) SetParent(parent_type->GetParticleName());
475 }
476 
478 {
479  // determine angular momentum
480 
481  // fill pointers to daughter particles if not yet set
482  if (G4MT_daughters == 0) FillDaughters();
483 
484  const G4int PiSpin = G4MT_parent->GetPDGiSpin();
485  const G4int PParity = G4MT_parent->GetPDGiParity();
486  if (2==numberOfDaughters) { // up to now we can only handle two particle decays
487  const G4int D1iSpin = G4MT_daughters[0]->GetPDGiSpin();
488  const G4int D1Parity = G4MT_daughters[0]->GetPDGiParity();
489  const G4int D2iSpin = G4MT_daughters[1]->GetPDGiSpin();
490  const G4int D2Parity = G4MT_daughters[1]->GetPDGiParity();
491  const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
492  const G4int MaxiSpin = D1iSpin + D2iSpin;
493  const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int
494  G4int lMin;
495 #ifdef G4VERBOSE
496  if (verboseLevel>1) {
497  G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
498  G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
499  }
500 #endif
501  for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings
502  lMin = std::abs(PiSpin-j)/2;
503 #ifdef G4VERBOSE
504  if (verboseLevel>1)
505  G4cout << "-> checking 2*j=" << j << G4endl;
506 #endif
507  for (G4int l=lMin; l<=lMax; l++) {
508 #ifdef G4VERBOSE
509  if (verboseLevel>1)
510  G4cout << " checking l=" << l << G4endl;
511 #endif
512  if (l%2==0) {
513  if (PParity == D1Parity*D2Parity) { // check parity for this l
514  return l;
515  }
516  } else {
517  if (PParity == -1*D1Parity*D2Parity) { // check parity for this l
518  return l;
519  }
520  }
521  }
522  }
523  } else {
524  G4Exception("G4VDecayChannel::GetAngularMomentum",
525  "PART111", JustWarning,
526  "Sorry, can't handle 3 particle decays (up to now)");
527  return 0;
528  }
529  G4Exception ("G4VDecayChannel::GetAngularMomentum",
530  "PART111", JustWarning,
531  "Can't find angular momentum for this decay");
532  return 0;
533 }
534 
536 {
537  G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
538  G4cout << " : " ;
539  for (G4int index=0; index < numberOfDaughters; index++)
540  {
541  if(daughters_name[index] != 0) {
542  G4cout << " " << *(daughters_name[index]);
543  } else {
544  G4cout << " not defined ";
545  }
546  }
547  G4cout << G4endl;
548 }
549 
551 {
552  return noName;
553 }
554 
555 #include "Randomize.hh"
557 {
558  if (maxDev >rangeMass) maxDev = rangeMass;
559  if (maxDev <-1.*rangeMass) return massPDG; // can not calculate
560 
561  G4double x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
562  G4double y = G4UniformRand();
563  while ( y * (width*width*x*x + massPDG*massPDG*width*width) > massPDG*massPDG*width*width ){
564  x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
565  y = G4UniformRand();
566  }
567  G4double mass = massPDG + x*width;
568  return mass;
569 }
const G4String & GetNoName() const
static const G4String noName
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
void SetBR(G4double value)
G4double * G4MT_daughters_mass
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
#define width
G4double G4MT_parent_mass
G4ParticleTable * particletable
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4int GetAngularMomentum()
G4double GetPDGWidth() const
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
G4String kinematics_name
void SetNumberOfDaughters(G4int value)
G4String * parent_name
static const double GeV
Definition: G4SIunits.hh:196
const G4String & GetParticleType() const
G4VDecayChannel & operator=(const G4VDecayChannel &)
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)
#define G4endl
Definition: G4ios.hh:61
G4double DynamicalMass(G4double massPDG, G4double width, G4double maxDev=+1.) const
double G4double
Definition: G4Types.hh:76
G4String ** daughters_name