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

Last change on this file since 1152 was 992, checked in by garnier, 17 years ago

fichiers oublies
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//
27// $Id: G4VDecayChannel.cc,v 1.18 2006/06/29 19:26:20 gunter Exp $
[992]28// GEANT4 tag $Name: geant4-09-02-ref-02 $
[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{
168 if (parent_name != 0) delete parent_name;
169 ClearDaughtersName();
170 if (daughters_mass != 0) delete [] daughters_mass;
171}
172
173void G4VDecayChannel::ClearDaughtersName()
174{
175 if ( daughters_name != 0) {
176 if (numberOfDaughters>0) {
177#ifdef G4VERBOSE
178 if (verboseLevel>1) {
179 G4cout << "G4VDecayChannel::ClearDaughtersName ";
180 G4cout << "clear all daughters " << G4endl;
181 }
182#endif
183 for (G4int index=0; index < numberOfDaughters; index++) {
184 if (daughters_name[index] != 0) delete daughters_name[index];
185 }
186 }
187 delete [] daughters_name;
188 daughters_name = 0;
189 }
190 //
191 if (daughters != 0) delete [] daughters;
192 if (daughters_mass != 0) delete [] daughters_mass;
193 daughters = 0;
194 daughters_mass = 0;
195
196 numberOfDaughters = 0;
197}
198
199void G4VDecayChannel::SetNumberOfDaughters(G4int size)
200{
201 if (size >0) {
202 // remove old contents
203 ClearDaughtersName();
204 // cleate array
205 daughters_name = new G4String*[size];
206 for (G4int index=0;index<size;index++) daughters_name[index]=0;
207 numberOfDaughters = size;
208 }
209}
210
211void G4VDecayChannel::SetDaughter(G4int anIndex,
212 const G4String &particle_name)
213{
214 // check numberOfDaughters is positive
215 if (numberOfDaughters<=0) {
216#ifdef G4VERBOSE
217 if (verboseLevel>0) {
218 G4cout << "G4VDecayChannel::SetDaughter: ";
219 G4cout << "Number of daughters is not defined" << G4endl;
220 }
221#endif
222 return;
223 }
224
225 // check existence of daughters_name array
226 if (daughters_name == 0) {
227 // cleate array
228 daughters_name = new G4String*[numberOfDaughters];
229 for (G4int index=0;index<numberOfDaughters;index++) {
230 daughters_name[index]=0;
231 }
232 }
233
234 // check an index
235 if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
236#ifdef G4VERBOSE
237 if (verboseLevel>0) {
238 G4cout << "G4VDecayChannel::SetDaughter";
239 G4cout << "index out of range " << anIndex << G4endl;
240 }
241#endif
242 } else {
243 // delete the old name if it exists
244 if (daughters_name[anIndex]!=0) delete daughters_name[anIndex];
245 // fill the name
246 daughters_name[anIndex] = new G4String(particle_name);
247 // refill the array of daughters[] if it exists
248 if (daughters != 0) FillDaughters();
249#ifdef G4VERBOSE
250 if (verboseLevel>1) {
251 G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
252 G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
253 }
254#endif
255 }
256}
257
258void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
259{
260 if (parent_type != 0) SetDaughter(anIndex, parent_type->GetParticleName());
261}
262
263void G4VDecayChannel::FillDaughters()
264{
265 G4int index;
266
267#ifdef G4VERBOSE
268 if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
269#endif
270 if (daughters != 0) delete [] daughters;
271
272 // parent mass
273 if (parent == 0) FillParent();
274 G4double parentmass = parent->GetPDGMass();
275
276 //
277 G4double sumofdaughtermass = 0.0;
278 if ((numberOfDaughters <=0) || (daughters_name == 0) ){
279#ifdef G4VERBOSE
280 if (verboseLevel>0) {
281 G4cout << "G4VDecayChannel::FillDaughters ";
282 G4cout << "[ " << parent->GetParticleName() << " ]";
283 G4cout << "numberOfDaughters is not defined yet";
284 }
285#endif
286 daughters = 0;
287 G4Exception("G4VDecayChannel::FillDaughters",
288 "can not fill daughters", FatalException,
289 "numberOfDaughters is not defined yet");
290 }
291
292 //create and set the array of pointers to daughter particles
293 daughters = new G4ParticleDefinition*[numberOfDaughters];
294 if (daughters_mass != 0) delete [] daughters_mass;
295 daughters_mass = new G4double[numberOfDaughters];
296 // loop over all daughters
297 for (index=0; index < numberOfDaughters; index++) {
298 if (daughters_name[index] == 0) {
299 // daughter name is not defined
300#ifdef G4VERBOSE
301 if (verboseLevel>0) {
302 G4cout << "G4VDecayChannel::FillDaughters ";
303 G4cout << "[ " << parent->GetParticleName() << " ]";
304 G4cout << index << "-th daughter is not defined yet" << G4endl;
305 }
306#endif
307 daughters[index] = 0;
308 G4Exception("G4VDecayChannel::FillDaughters",
309 "can not fill daughters", FatalException,
310 "name of a daughter is not defined yet");
311 }
312 //search daughter particles in the particle table
313 daughters[index] = particletable->FindParticle(*daughters_name[index]);
314 if (daughters[index] == 0) {
315 // can not find the daughter particle
316#ifdef G4VERBOSE
317 if (verboseLevel>0) {
318 G4cout << "G4VDecayChannel::FillDaughters ";
319 G4cout << "[ " << parent->GetParticleName() << " ]";
320 G4cout << index << ":" << *daughters_name[index];
321 G4cout << " is not defined !!" << G4endl;
322 G4cout << " The BR of this decay mode is set to zero " << G4endl;
323 }
324#endif
325 SetBR(0.0);
326 return;
327 }
328#ifdef G4VERBOSE
329 if (verboseLevel>1) {
330 G4cout << index << ":" << *daughters_name[index];
331 G4cout << ":" << daughters[index] << G4endl;
332 }
333#endif
334 daughters_mass[index] = daughters[index]->GetPDGMass();
335 sumofdaughtermass += daughters[index]->GetPDGMass();
336 } // end loop over all daughters
337
338 // check sum of daghter mass
339 G4double widthMass = parent->GetPDGWidth();
340 if ( (parent->GetParticleType() != "nucleus") &&
341 (sumofdaughtermass > parentmass + 5*widthMass) ){
342 // !!! illegal mass !!!
343#ifdef G4VERBOSE
344 if (GetVerboseLevel()>0) {
345 G4cout << "G4VDecayChannel::FillDaughters ";
346 G4cout << "[ " << parent->GetParticleName() << " ]";
347 G4cout << " Energy/Momentum conserevation breaks " <<G4endl;
348 if (GetVerboseLevel()>1) {
349 G4cout << " parent:" << *parent_name;
350 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
351 for (index=0; index < numberOfDaughters; index++){
352 G4cout << " daughter " << index << ":" << *daughters_name[index];
353 G4cout << " mass:" << daughters[index]->GetPDGMass()/GeV;
354 G4cout << "[GeV/c/c]" <<G4endl;
355 }
356 }
357 }
358#endif
359 }
360}
361
362
363void G4VDecayChannel::FillParent()
364{
365 if (parent_name == 0) {
366 // parent name is not defined
367#ifdef G4VERBOSE
368 if (verboseLevel>0) {
369 G4cout << "G4VDecayChannel::FillParent ";
370 G4cout << ": parent name is not defined !!" << G4endl;
371 }
372#endif
373 parent = 0;
374 G4Exception("G4VDecayChannel::FillParent()",
375 "can not fill parent", FatalException,
376 "parent name is not defined yet");
377 }
378 // search parent particle in the particle table
379 parent = particletable->FindParticle(*parent_name);
380 if (parent == 0) {
381 // parent particle does not exist
382#ifdef G4VERBOSE
383 if (verboseLevel>0) {
384 G4cout << "G4VDecayChannel::FillParent ";
385 G4cout << *parent_name << " does not exist !!" << G4endl;
386 }
387#endif
388 G4Exception("G4VDecayChannel::FillParent()",
389 "can not fill parent", FatalException,
390 "parent does not exist");
391 }
392 parent_mass = parent->GetPDGMass();
393}
394
395void G4VDecayChannel::SetParent(const G4ParticleDefinition * parent_type)
396{
397 if (parent_type != 0) SetParent(parent_type->GetParticleName());
398}
399
400G4int G4VDecayChannel::GetAngularMomentum()
401{
402 // determine angular momentum
403
404 // fill pointers to daughter particles if not yet set
405 if (daughters == 0) FillDaughters();
406
407 const G4int PiSpin = parent->GetPDGiSpin();
408 const G4int PParity = parent->GetPDGiParity();
409 if (2==numberOfDaughters) { // up to now we can only handle two particle decays
410 const G4int D1iSpin = daughters[0]->GetPDGiSpin();
411 const G4int D1Parity = daughters[0]->GetPDGiParity();
412 const G4int D2iSpin = daughters[1]->GetPDGiSpin();
413 const G4int D2Parity = daughters[1]->GetPDGiParity();
414 const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
415 const G4int MaxiSpin = D1iSpin + D2iSpin;
416 const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int
417 G4int lMin;
418#ifdef G4VERBOSE
419 if (verboseLevel>1) {
420 G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
421 G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
422 }
423#endif
424 for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings
425 lMin = std::abs(PiSpin-j)/2;
426#ifdef G4VERBOSE
427 if (verboseLevel>1)
428 G4cout << "-> checking 2*j=" << j << G4endl;
429#endif
430 for (G4int l=lMin; l<=lMax; l++) {
431#ifdef G4VERBOSE
432 if (verboseLevel>1)
433 G4cout << " checking l=" << l << G4endl;
434#endif
435 if (l%2==0) {
436 if (PParity == D1Parity*D2Parity) { // check parity for this l
437 return l;
438 }
439 } else {
440 if (PParity == -1*D1Parity*D2Parity) { // check parity for this l
441 return l;
442 }
443 }
444 }
445 }
446 } else {
447 G4Exception("G4VDecayChannel::GetAngularMomentum",
448 "can not calculate", JustWarning,
449 "Sorry, can't handle 3 particle decays (up to now)");
450 return 0;
451 }
452 G4Exception ("G4VDecayChannel::GetAngularMomentum",
453 "can not calculate", JustWarning,
454 "Can't find angular momentum for this decay!");
455 return 0;
456}
457
458void G4VDecayChannel::DumpInfo()
459{
460 G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
461 G4cout << " : " ;
462 for (G4int index=0; index < numberOfDaughters; index++)
463 {
464 if(daughters_name[index] != 0) {
465 G4cout << " " << *(daughters_name[index]);
466 } else {
467 G4cout << " not defined ";
468 }
469 }
470 G4cout << G4endl;
471}
472
473const G4String& G4VDecayChannel::GetNoName() const
474{
475 return noName;
476}
Note: See TracBrowser for help on using the repository browser.