ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VDecayChannel.cc
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 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, based on object model of
33 // 27 July 1996 H.Kurashige
34 // 30 May 1997 H.Kurashige
35 // 23 Mar. 2000 H.Weber : add GetAngularMomentum
36 // ------------------------------------------------------------
37 
38 #include "G4ParticleDefinition.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4ParticleTable.hh"
41 #include "G4DecayTable.hh"
42 #include "G4DecayProducts.hh"
43 #include "G4VDecayChannel.hh"
44 #include "G4AutoLock.hh"
45 
47 
49  :kinematics_name(""),
50  rbranch(0.0),
51  numberOfDaughters(0),
52  parent_name(nullptr),
53  daughters_name(nullptr),
54  rangeMass(2.5),
55  parent_polarization(),
56  particletable(nullptr),
57  verboseLevel(1)
58 {
59  G4MT_parent = nullptr;
60  G4MT_daughters = nullptr;
61  G4MT_parent_mass = 0.0;
62  G4MT_daughters_mass = nullptr;
63  G4MT_daughters_width = nullptr;
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(nullptr),
74  daughters_name(nullptr),
75  rangeMass(2.5),
76  parent_polarization(),
77  particletable(nullptr),
78  verboseLevel(Verbose)
79 {
80  G4MT_parent = nullptr;
81  G4MT_daughters = nullptr;
82  G4MT_parent_mass = 0.0;
83  G4MT_daughters_mass = nullptr;
84  G4MT_daughters_width = nullptr;
85 
86  // set pointer to G4ParticleTable (static and singleton object)
88 }
89 
91  const G4String& theParentName,
92  G4double theBR,
93  G4int theNumberOfDaughters,
94  const G4String& theDaughterName1,
95  const G4String& theDaughterName2,
96  const G4String& theDaughterName3,
97  const G4String& theDaughterName4 )
98  :kinematics_name(aName),
99  rbranch(theBR),
100  numberOfDaughters(theNumberOfDaughters),
101  parent_name(nullptr),
102  daughters_name(nullptr),
103  rangeMass(1.0),
104  parent_polarization(),
105  particletable(0),
106  verboseLevel(1)
107 {
108  G4MT_parent = nullptr;
109  G4MT_daughters = nullptr;
110  G4MT_parent_mass = 0.0;
111  G4MT_daughters_mass = nullptr;
112  G4MT_daughters_width = nullptr;
113 
114  // set pointer to G4ParticleTable (static and singleton object)
116 
117  // parent name
118  parent_name = new G4String(theParentName);
119 
120  // cleate array
122  for (G4int index=0;index<numberOfDaughters;index++) daughters_name[index]=0;
123 
124  // daughters' name
125  if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
126  if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
127  if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
128  if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
129 
130  if (rbranch <0. ) rbranch = 0.0;
131  else if (rbranch >1.0 ) rbranch = 1.0;
132 }
133 
135 {
137  verboseLevel = right.verboseLevel;
138  rbranch = right.rbranch;
139  rangeMass = right.rangeMass;
140 
141  // copy parent name
142  parent_name = new G4String(*right.parent_name);
143  G4MT_parent = nullptr;
144  G4MT_parent_mass = 0.0;
145 
146  //create array
148 
149  daughters_name =nullptr;
150  if ( numberOfDaughters >0 ) {
152  //copy daughters name
153  for (G4int index=0; index < numberOfDaughters; index++){
154  daughters_name[index] = new G4String(*right.daughters_name[index]);
155  }
156  }
157 
158  //
159  G4MT_daughters_mass = nullptr;
160  G4MT_daughters = nullptr;
161  G4MT_daughters_width = nullptr;
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 ) {
189  if (daughters_name != nullptr) ClearDaughtersName();
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 = nullptr;
200  G4MT_daughters = nullptr;
201  G4MT_parent_mass = 0.0;
202  G4MT_daughters_mass = nullptr;
203  G4MT_daughters_width = nullptr;
204 
205  // particle table
207 
210 
211  return *this;
212 }
213 
215 {
217  if (parent_name != nullptr) delete parent_name;
218  parent_name = nullptr;
219  if (G4MT_daughters_mass != nullptr) delete [] G4MT_daughters_mass;
220  G4MT_daughters_mass =nullptr;
221  if (G4MT_daughters_width != nullptr) delete [] G4MT_daughters_width;
222  G4MT_daughters_width = nullptr;
225 }
226 
228 {
230  if ( daughters_name != nullptr) {
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] != nullptr) delete daughters_name[index];
240  }
241  }
242  delete [] daughters_name;
243  daughters_name = nullptr;
244  }
245  //
246  if (G4MT_daughters != nullptr) delete [] G4MT_daughters;
247  if (G4MT_daughters_mass != nullptr) delete [] G4MT_daughters_mass;
248  if (G4MT_daughters_width != nullptr) delete [] G4MT_daughters_width;
249  G4MT_daughters_width = nullptr;
250  G4MT_daughters = nullptr;
251  G4MT_daughters_mass = nullptr;
252 
253  numberOfDaughters = 0;
254 }
255 
257 {
258  if (size >0) {
259  // remove old contents
261  // cleate array
262  daughters_name = new G4String*[size];
263  for (G4int index=0;index<size;index++) daughters_name[index]=nullptr;
264  numberOfDaughters = size;
265  }
266 }
267 
269  const G4String &particle_name)
270 {
271  // check numberOfDaughters is positive
272  if (numberOfDaughters<=0) {
273 #ifdef G4VERBOSE
274  if (verboseLevel>0) {
275  G4cout << "G4VDecayChannel::SetDaughter: "
276  << "Number of daughters is not defined" << G4endl;
277  }
278 #endif
279  return;
280  }
281 
282  //ANDREA:-> Feb 25 2016
283  // An analysis of this code, shows that this method is called
284  // only in the constructor of derived classes.
285  // The general idea of this method is probably to support
286  // the possibility to re-define daughters on the fly, however
287  // this design is extremely problematic for MT mode, we thus
288  // require (as practically happens) that the method is called only
289  // at construction, i.e. when G4MT_daugheters == 0
290  // moreover this method can be called only after SetNumberOfDaugthers
291  // has been called (see previous if), in such a case daughters_name != 0
292  if ( daughters_name == nullptr ) {
293  G4Exception("G4VDecayChannel::SetDaughter","PART112",FatalException,
294  "Trying to add a daughter without specifying number of secondaries, useSetNumberOfDaughters first");
295  return;
296  }
297  if ( G4MT_daughters != nullptr ) {
298  G4Exception("G4VDecayChannel::SetDaughter","PART111",FatalException,
299  "Trying to modify a daughter of a decay channel, but decay channel already has daughters.");
300  return;
301  }
302  //<-:ANDREA
303 
304  // check an index
305  if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
306 #ifdef G4VERBOSE
307  if (verboseLevel>0) {
308  G4cout << "G4VDecayChannel::SetDaughter"
309  << "index out of range " << anIndex << G4endl;
310  }
311 #endif
312  } else {
313  // fill the name
314  daughters_name[anIndex] = new G4String(particle_name);
315 #ifdef G4VERBOSE
316  if (verboseLevel>1) {
317  G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
318  G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
319  }
320 #endif
321  }
322 }
323 
324 void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
325 {
326  if (parent_type != nullptr) SetDaughter(anIndex, parent_type->GetParticleName());
327 }
328 
330 {
331  G4AutoLock lock(&daughtersMutex);
332  //Double check, check again if another thread has already filled this, in
333  //case do not need to do anything
334  if ( G4MT_daughters != nullptr ) return;
335 
336  G4int index;
337 
338 #ifdef G4VERBOSE
339  if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
340 #endif
341  if (G4MT_daughters != nullptr) {
342  delete [] G4MT_daughters;
343  G4MT_daughters = nullptr;
344  }
345 
346  // parent mass
348  G4double parentmass = G4MT_parent->GetPDGMass();
349 
350  //
351  G4double sumofdaughtermass = 0.0;
352  G4double sumofdaughterwidthsq = 0.0;
353 
354  if ((numberOfDaughters <=0) || (daughters_name == nullptr) ){
355 #ifdef G4VERBOSE
356  if (verboseLevel>0) {
357  G4cout << "G4VDecayChannel::FillDaughters "
358  << "[ " << G4MT_parent->GetParticleName() << " ]"
359  << "numberOfDaughters is not defined yet";
360  }
361 #endif
362  G4MT_daughters = nullptr;
363  G4Exception("G4VDecayChannel::FillDaughters",
364  "PART011", FatalException,
365  "Can not fill daughters: numberOfDaughters is not defined yet");
366  }
367 
368  //create and set the array of pointers to daughter particles
370  if (G4MT_daughters_mass != nullptr) delete [] G4MT_daughters_mass;
371  if (G4MT_daughters_width != nullptr) delete [] G4MT_daughters_width;
374  // loop over all daughters
375  for (index=0; index < numberOfDaughters; index++) {
376  if (daughters_name[index] == nullptr) {
377  // daughter name is not defined
378 #ifdef G4VERBOSE
379  if (verboseLevel>0) {
380  G4cout << "G4VDecayChannel::FillDaughters "
381  << "[ " << G4MT_parent->GetParticleName() << " ]"
382  << index << "-th daughter is not defined yet" << G4endl;
383  }
384 #endif
385  G4MT_daughters[index] = nullptr;
386  G4Exception("G4VDecayChannel::FillDaughters",
387  "PART011", FatalException,
388  "Can not fill daughters: name of a daughter is not defined yet");
389  }
390  //search daughter particles in the particle table
392  if (G4MT_daughters[index] == nullptr ) {
393  // can not find the daughter particle
394 #ifdef G4VERBOSE
395  if (verboseLevel>0) {
396  G4cout << "G4VDecayChannel::FillDaughters "
397  << "[ " << G4MT_parent->GetParticleName() << " ]"
398  << index << ":" << *daughters_name[index]
399  << " is not defined !!" << G4endl;
400  G4cout << " The BR of this decay mode is set to zero " << G4endl;
401  }
402 #endif
403  SetBR(0.0);
404  return;
405  }
406 #ifdef G4VERBOSE
407  if (verboseLevel>1) {
408  G4cout << index << ":" << *daughters_name[index];
409  G4cout << ":" << G4MT_daughters[index] << G4endl;
410  }
411 #endif
412  G4MT_daughters_mass[index] = G4MT_daughters[index]->GetPDGMass();
413  G4double d_width = G4MT_daughters[index]->GetPDGWidth();
414  G4MT_daughters_width[index] = d_width;
415  sumofdaughtermass += G4MT_daughters[index]->GetPDGMass();
416  sumofdaughterwidthsq += d_width*d_width;
417  } // end loop over all daughters
418 
419  // check sum of daghter mass
420  G4double widthMass = std::sqrt(G4MT_parent->GetPDGWidth()*G4MT_parent->GetPDGWidth()+sumofdaughterwidthsq);
421  if ( (G4MT_parent->GetParticleType() != "nucleus") &&
422  (numberOfDaughters !=1) &&
423  (sumofdaughtermass > parentmass + rangeMass*widthMass) ){
424  // !!! illegal mass !!!
425 #ifdef G4VERBOSE
426  if (GetVerboseLevel()>0) {
427  G4cout << "G4VDecayChannel::FillDaughters "
428  << "[ " << G4MT_parent->GetParticleName() << " ]"
429  << " Energy/Momentum conserevation breaks " <<G4endl;
430  if (GetVerboseLevel()>1) {
431  G4cout << " parent:" << *parent_name
432  << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
433  for (index=0; index < numberOfDaughters; index++){
434  G4cout << " daughter " << index << ":" << *daughters_name[index]
435  << " mass:" << G4MT_daughters[index]->GetPDGMass()/GeV
436  << "[GeV/c/c]" <<G4endl;
437  }
438  }
439  }
440 #endif
441  }
442 }
443 
444 
446 {
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 != nullptr ) return;
451 
452  if (parent_name == nullptr) {
453  // parent name is not defined
454 #ifdef G4VERBOSE
455  if (verboseLevel>0) {
456  G4cout << "G4VDecayChannel::FillParent "
457  << ": parent name is not defined !!" << G4endl;
458  }
459 #endif
460  G4MT_parent = nullptr;
461  G4Exception("G4VDecayChannel::FillParent()",
462  "PART012", FatalException,
463  "Can not fill parent: parent name is not defined yet");
464  return;
465  }
466  // search parent particle in the particle table
468  if (G4MT_parent == nullptr) {
469  // parent particle does not exist
470 #ifdef G4VERBOSE
471  if (verboseLevel>0) {
472  G4cout << "G4VDecayChannel::FillParent "
473  << *parent_name << " does not exist !!" << G4endl;
474  }
475 #endif
476  G4Exception("G4VDecayChannel::FillParent()",
477  "PART012", FatalException,
478  "Can not fill parent: parent does not exist");
479  return;
480  }
482 }
483 
485 {
486  if (parent_type != nullptr) SetParent(parent_type->GetParticleName());
487 }
488 
490 {
491  // determine angular momentum
492 
493  // fill pointers to daughter particles if not yet set
495 
496  const G4int PiSpin = G4MT_parent->GetPDGiSpin();
497  const G4int PParity = G4MT_parent->GetPDGiParity();
498  if (2==numberOfDaughters) { // up to now we can only handle two particle decays
499  const G4int D1iSpin = G4MT_daughters[0]->GetPDGiSpin();
500  const G4int D1Parity = G4MT_daughters[0]->GetPDGiParity();
501  const G4int D2iSpin = G4MT_daughters[1]->GetPDGiSpin();
502  const G4int D2Parity = G4MT_daughters[1]->GetPDGiParity();
503  const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
504  const G4int MaxiSpin = D1iSpin + D2iSpin;
505  const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int
506  G4int lMin;
507 #ifdef G4VERBOSE
508  if (verboseLevel>1) {
509  G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
510  G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
511  }
512 #endif
513  for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings
514  lMin = std::abs(PiSpin-j)/2;
515 #ifdef G4VERBOSE
516  if (verboseLevel>1)
517  G4cout << "-> checking 2*j=" << j << G4endl;
518 #endif
519  for (G4int l=lMin; l<=lMax; l++) {
520 #ifdef G4VERBOSE
521  if (verboseLevel>1)
522  G4cout << " checking l=" << l << G4endl;
523 #endif
524  if (l%2==0) {
525  if (PParity == D1Parity*D2Parity) { // check parity for this l
526  return l;
527  }
528  } else {
529  if (PParity == -1*D1Parity*D2Parity) { // check parity for this l
530  return l;
531  }
532  }
533  }
534  }
535  } else {
536  G4Exception("G4VDecayChannel::GetAngularMomentum",
537  "PART111", JustWarning,
538  "Sorry, can't handle 3 particle decays (up to now)");
539  return 0;
540  }
541  G4Exception ("G4VDecayChannel::GetAngularMomentum",
542  "PART111", JustWarning,
543  "Can't find angular momentum for this decay");
544  return 0;
545 }
546 
548 {
549  G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
550  G4cout << " : " ;
551  for (G4int index=0; index < numberOfDaughters; index++){
552  if(daughters_name[index] != nullptr) {
553  G4cout << " " << *(daughters_name[index]);
554  } else {
555  G4cout << " not defined ";
556  }
557  }
558  G4cout << G4endl;
559 }
560 
562 {
563  return noName;
564 }
565 
566 #include "Randomize.hh"
568 {
569  if (width<=0.0) return massPDG;
570  if (maxDev >rangeMass) maxDev = rangeMass;
571  if (maxDev <=-1.*rangeMass) return massPDG; // can not calculate
572 
573  G4double x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
575  const size_t MAX_LOOP=10000;
576  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
577  if ( y * (width*width*x*x + massPDG*massPDG*width*width) <= massPDG*massPDG*width*width ) break;
578  x = G4UniformRand()*(maxDev+rangeMass) - rangeMass;
579  y = G4UniformRand();
580  }
581  G4double mass = massPDG + x*width;
582  return mass;
583 }
584 
586 {
587  G4double sumOfDaughterMassMin=0.0;
590  // skip one body decay
591  if (numberOfDaughters==1) return true;
592 
593  for (G4int index=0; index < numberOfDaughters; index++) {
594  sumOfDaughterMassMin +=
596  }
597  return (parentMass >= sumOfDaughterMassMin);
598 }
599 
601 {
602  rbranch = value;
603  if (rbranch <0. ) rbranch = 0.0;
604  else if (rbranch >1.0 ) rbranch = 1.0;
605 }
606