source: trunk/source/particles/management/src/G4VDecayChannel.cc@ 1260

Last change on this file since 1260 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 14.7 KB
RevLine 
[824]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//
[1196]27// $Id: G4VDecayChannel.cc,v 1.19 2009/08/17 14:52:19 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
[824]29//
30//
31// ------------------------------------------------------------
32// GEANT 4 class header file
33//
34// History: first implementation, based on object model of
35// 27 July 1996 H.Kurashige
36// 30 May 1997 H.Kurashige
37// 23 Mar. 2000 H.Weber : add GetAngularMomentum
38// ------------------------------------------------------------
39
40#include "G4ParticleDefinition.hh"
41#include "G4ParticleTable.hh"
42#include "G4DecayTable.hh"
43#include "G4DecayProducts.hh"
44#include "G4VDecayChannel.hh"
45
46const G4String G4VDecayChannel::noName = " ";
47
48G4VDecayChannel::G4VDecayChannel(const G4String &aName, G4int Verbose)
49 :kinematics_name(aName),
50 rbranch(0.0),
51 numberOfDaughters(0),
52 parent_name(0), daughters_name(0),
53 particletable(0),
54 parent(0), daughters(0),
55 parent_mass(0.0), daughters_mass(0),
56 verboseLevel(Verbose)
57{
58 // set pointer to G4ParticleTable (static and singleton object)
59 particletable = G4ParticleTable::GetParticleTable();
60}
61
62G4VDecayChannel::G4VDecayChannel(const G4String &aName,
63 const G4String& theParentName,
64 G4double theBR,
65 G4int theNumberOfDaughters,
66 const G4String& theDaughterName1,
67 const G4String& theDaughterName2,
68 const G4String& theDaughterName3,
69 const G4String& theDaughterName4 )
70 :kinematics_name(aName),
71 rbranch(theBR),
72 numberOfDaughters(theNumberOfDaughters),
73 parent_name(0), daughters_name(0),
74 particletable(0),
75 parent(0), daughters(0),
76 parent_mass(0.0), daughters_mass(0),
77 verboseLevel(1)
78{
79 // set pointer to G4ParticleTable (static and singleton object)
80 particletable = G4ParticleTable::GetParticleTable();
81
82 // parent name
83 parent_name = new G4String(theParentName);
84
85 // cleate array
86 daughters_name = new G4String*[numberOfDaughters];
87 for (G4int index=0;index<numberOfDaughters;index++) daughters_name[index]=0;
88
89 // daughters' name
90 if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
91 if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
92 if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
93 if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
94}
95
96
97
98G4VDecayChannel::G4VDecayChannel(const G4VDecayChannel &right)
99{
100 kinematics_name = right.kinematics_name;
101 verboseLevel = right.verboseLevel;
102 rbranch = right.rbranch;
103
104 // copy parent name
105 parent_name = new G4String(*right.parent_name);
106 parent = 0;
107 parent_mass = 0.0;
108
109 //create array
110 numberOfDaughters = right.numberOfDaughters;
111
112 if ( numberOfDaughters >0 ) {
113 daughters_name = new G4String*[numberOfDaughters];
114 //copy daughters name
115 for (G4int index=0; index < numberOfDaughters; index++)
116 {
117 daughters_name[index] = new G4String(*right.daughters_name[index]);
118 }
119 }
120
121 //
122 daughters_mass = 0;
123 daughters = 0;
124
125 // particle table
126 particletable = G4ParticleTable::GetParticleTable();
127}
128
129G4VDecayChannel & G4VDecayChannel::operator=(const G4VDecayChannel &right)
130{
131 if (this != &right) {
132 kinematics_name = right.kinematics_name;
133 verboseLevel = right.verboseLevel;
134 rbranch = right.rbranch;
135
136 // copy parent name
137 parent_name = new G4String(*right.parent_name);
138
139 // clear daughters_name array
140 ClearDaughtersName();
141
142 // recreate array
143 numberOfDaughters = right.numberOfDaughters;
144 if ( numberOfDaughters >0 ) {
145 daughters_name = new G4String*[numberOfDaughters];
146 //copy daughters name
147 for (G4int index=0; index < numberOfDaughters; index++) {
148 daughters_name[index] = new G4String(*right.daughters_name[index]);
149 }
150 }
151 }
152
153 //
154 parent = 0;
155 daughters = 0;
156 parent_mass = 0.0;
157 daughters_mass = 0;
158
159 // particle table
160 particletable = G4ParticleTable::GetParticleTable();
161
162 return *this;
163}
164
165
166G4VDecayChannel::~G4VDecayChannel()
167{
[1196]168 ClearDaughtersName();
[824]169 if (parent_name != 0) delete parent_name;
[1196]170 parent_name = 0;
[824]171 if (daughters_mass != 0) delete [] daughters_mass;
[1196]172 daughters_mass =0;
[824]173}
174
175void G4VDecayChannel::ClearDaughtersName()
176{
177 if ( daughters_name != 0) {
178 if (numberOfDaughters>0) {
179#ifdef G4VERBOSE
180 if (verboseLevel>1) {
[1196]181 G4cerr << "G4VDecayChannel::ClearDaughtersName "
182 << " for " << *parent_name << G4endl;
[824]183 }
184#endif
185 for (G4int index=0; index < numberOfDaughters; index++) {
186 if (daughters_name[index] != 0) delete daughters_name[index];
187 }
188 }
189 delete [] daughters_name;
190 daughters_name = 0;
191 }
192 //
193 if (daughters != 0) delete [] daughters;
194 if (daughters_mass != 0) delete [] daughters_mass;
195 daughters = 0;
196 daughters_mass = 0;
197
198 numberOfDaughters = 0;
199}
200
201void G4VDecayChannel::SetNumberOfDaughters(G4int size)
202{
203 if (size >0) {
204 // remove old contents
205 ClearDaughtersName();
206 // cleate array
207 daughters_name = new G4String*[size];
208 for (G4int index=0;index<size;index++) daughters_name[index]=0;
209 numberOfDaughters = size;
210 }
211}
212
213void G4VDecayChannel::SetDaughter(G4int anIndex,
214 const G4String &particle_name)
215{
216 // check numberOfDaughters is positive
217 if (numberOfDaughters<=0) {
218#ifdef G4VERBOSE
219 if (verboseLevel>0) {
220 G4cout << "G4VDecayChannel::SetDaughter: ";
221 G4cout << "Number of daughters is not defined" << G4endl;
222 }
223#endif
224 return;
225 }
226
227 // check existence of daughters_name array
228 if (daughters_name == 0) {
229 // cleate array
230 daughters_name = new G4String*[numberOfDaughters];
231 for (G4int index=0;index<numberOfDaughters;index++) {
232 daughters_name[index]=0;
233 }
234 }
235
236 // check an index
237 if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
238#ifdef G4VERBOSE
239 if (verboseLevel>0) {
240 G4cout << "G4VDecayChannel::SetDaughter";
241 G4cout << "index out of range " << anIndex << G4endl;
242 }
243#endif
244 } else {
245 // delete the old name if it exists
246 if (daughters_name[anIndex]!=0) delete daughters_name[anIndex];
247 // fill the name
248 daughters_name[anIndex] = new G4String(particle_name);
249 // refill the array of daughters[] if it exists
250 if (daughters != 0) FillDaughters();
251#ifdef G4VERBOSE
252 if (verboseLevel>1) {
253 G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
254 G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
255 }
256#endif
257 }
258}
259
260void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
261{
262 if (parent_type != 0) SetDaughter(anIndex, parent_type->GetParticleName());
263}
264
265void G4VDecayChannel::FillDaughters()
266{
267 G4int index;
268
269#ifdef G4VERBOSE
270 if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
271#endif
272 if (daughters != 0) delete [] daughters;
273
274 // parent mass
275 if (parent == 0) FillParent();
276 G4double parentmass = parent->GetPDGMass();
277
278 //
279 G4double sumofdaughtermass = 0.0;
280 if ((numberOfDaughters <=0) || (daughters_name == 0) ){
281#ifdef G4VERBOSE
282 if (verboseLevel>0) {
283 G4cout << "G4VDecayChannel::FillDaughters ";
284 G4cout << "[ " << parent->GetParticleName() << " ]";
285 G4cout << "numberOfDaughters is not defined yet";
286 }
287#endif
288 daughters = 0;
289 G4Exception("G4VDecayChannel::FillDaughters",
290 "can not fill daughters", FatalException,
291 "numberOfDaughters is not defined yet");
292 }
293
294 //create and set the array of pointers to daughter particles
295 daughters = new G4ParticleDefinition*[numberOfDaughters];
296 if (daughters_mass != 0) delete [] daughters_mass;
297 daughters_mass = new G4double[numberOfDaughters];
298 // loop over all daughters
299 for (index=0; index < numberOfDaughters; index++) {
300 if (daughters_name[index] == 0) {
301 // daughter name is not defined
302#ifdef G4VERBOSE
303 if (verboseLevel>0) {
304 G4cout << "G4VDecayChannel::FillDaughters ";
305 G4cout << "[ " << parent->GetParticleName() << " ]";
306 G4cout << index << "-th daughter is not defined yet" << G4endl;
307 }
308#endif
309 daughters[index] = 0;
310 G4Exception("G4VDecayChannel::FillDaughters",
311 "can not fill daughters", FatalException,
312 "name of a daughter is not defined yet");
313 }
314 //search daughter particles in the particle table
315 daughters[index] = particletable->FindParticle(*daughters_name[index]);
316 if (daughters[index] == 0) {
317 // can not find the daughter particle
318#ifdef G4VERBOSE
319 if (verboseLevel>0) {
320 G4cout << "G4VDecayChannel::FillDaughters ";
321 G4cout << "[ " << parent->GetParticleName() << " ]";
322 G4cout << index << ":" << *daughters_name[index];
323 G4cout << " is not defined !!" << G4endl;
324 G4cout << " The BR of this decay mode is set to zero " << G4endl;
325 }
326#endif
327 SetBR(0.0);
328 return;
329 }
330#ifdef G4VERBOSE
331 if (verboseLevel>1) {
332 G4cout << index << ":" << *daughters_name[index];
333 G4cout << ":" << daughters[index] << G4endl;
334 }
335#endif
336 daughters_mass[index] = daughters[index]->GetPDGMass();
337 sumofdaughtermass += daughters[index]->GetPDGMass();
338 } // end loop over all daughters
339
340 // check sum of daghter mass
341 G4double widthMass = parent->GetPDGWidth();
342 if ( (parent->GetParticleType() != "nucleus") &&
343 (sumofdaughtermass > parentmass + 5*widthMass) ){
344 // !!! illegal mass !!!
345#ifdef G4VERBOSE
346 if (GetVerboseLevel()>0) {
347 G4cout << "G4VDecayChannel::FillDaughters ";
348 G4cout << "[ " << parent->GetParticleName() << " ]";
349 G4cout << " Energy/Momentum conserevation breaks " <<G4endl;
350 if (GetVerboseLevel()>1) {
351 G4cout << " parent:" << *parent_name;
352 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
353 for (index=0; index < numberOfDaughters; index++){
354 G4cout << " daughter " << index << ":" << *daughters_name[index];
355 G4cout << " mass:" << daughters[index]->GetPDGMass()/GeV;
356 G4cout << "[GeV/c/c]" <<G4endl;
357 }
358 }
359 }
360#endif
361 }
362}
363
364
365void G4VDecayChannel::FillParent()
366{
367 if (parent_name == 0) {
368 // parent name is not defined
369#ifdef G4VERBOSE
370 if (verboseLevel>0) {
371 G4cout << "G4VDecayChannel::FillParent ";
372 G4cout << ": parent name is not defined !!" << G4endl;
373 }
374#endif
375 parent = 0;
376 G4Exception("G4VDecayChannel::FillParent()",
377 "can not fill parent", FatalException,
378 "parent name is not defined yet");
379 }
380 // search parent particle in the particle table
381 parent = particletable->FindParticle(*parent_name);
382 if (parent == 0) {
383 // parent particle does not exist
384#ifdef G4VERBOSE
385 if (verboseLevel>0) {
386 G4cout << "G4VDecayChannel::FillParent ";
387 G4cout << *parent_name << " does not exist !!" << G4endl;
388 }
389#endif
390 G4Exception("G4VDecayChannel::FillParent()",
391 "can not fill parent", FatalException,
392 "parent does not exist");
393 }
394 parent_mass = parent->GetPDGMass();
395}
396
397void G4VDecayChannel::SetParent(const G4ParticleDefinition * parent_type)
398{
399 if (parent_type != 0) SetParent(parent_type->GetParticleName());
400}
401
402G4int G4VDecayChannel::GetAngularMomentum()
403{
404 // determine angular momentum
405
406 // fill pointers to daughter particles if not yet set
407 if (daughters == 0) FillDaughters();
408
409 const G4int PiSpin = parent->GetPDGiSpin();
410 const G4int PParity = parent->GetPDGiParity();
411 if (2==numberOfDaughters) { // up to now we can only handle two particle decays
412 const G4int D1iSpin = daughters[0]->GetPDGiSpin();
413 const G4int D1Parity = daughters[0]->GetPDGiParity();
414 const G4int D2iSpin = daughters[1]->GetPDGiSpin();
415 const G4int D2Parity = daughters[1]->GetPDGiParity();
416 const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
417 const G4int MaxiSpin = D1iSpin + D2iSpin;
418 const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int
419 G4int lMin;
420#ifdef G4VERBOSE
421 if (verboseLevel>1) {
422 G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
423 G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
424 }
425#endif
426 for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings
427 lMin = std::abs(PiSpin-j)/2;
428#ifdef G4VERBOSE
429 if (verboseLevel>1)
430 G4cout << "-> checking 2*j=" << j << G4endl;
431#endif
432 for (G4int l=lMin; l<=lMax; l++) {
433#ifdef G4VERBOSE
434 if (verboseLevel>1)
435 G4cout << " checking l=" << l << G4endl;
436#endif
437 if (l%2==0) {
438 if (PParity == D1Parity*D2Parity) { // check parity for this l
439 return l;
440 }
441 } else {
442 if (PParity == -1*D1Parity*D2Parity) { // check parity for this l
443 return l;
444 }
445 }
446 }
447 }
448 } else {
449 G4Exception("G4VDecayChannel::GetAngularMomentum",
450 "can not calculate", JustWarning,
451 "Sorry, can't handle 3 particle decays (up to now)");
452 return 0;
453 }
454 G4Exception ("G4VDecayChannel::GetAngularMomentum",
455 "can not calculate", JustWarning,
456 "Can't find angular momentum for this decay!");
457 return 0;
458}
459
460void G4VDecayChannel::DumpInfo()
461{
462 G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
463 G4cout << " : " ;
464 for (G4int index=0; index < numberOfDaughters; index++)
465 {
466 if(daughters_name[index] != 0) {
467 G4cout << " " << *(daughters_name[index]);
468 } else {
469 G4cout << " not defined ";
470 }
471 }
472 G4cout << G4endl;
473}
474
475const G4String& G4VDecayChannel::GetNoName() const
476{
477 return noName;
478}
Note: See TracBrowser for help on using the repository browser.