00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4ParticleTable.hh"
00042 #include "G4DecayTable.hh"
00043 #include "G4DecayProducts.hh"
00044 #include "G4VDecayChannel.hh"
00045
00046 const G4String G4VDecayChannel::noName = " ";
00047
00048 G4VDecayChannel::G4VDecayChannel()
00049 :kinematics_name(""),
00050 rbranch(0.0),
00051 numberOfDaughters(0),
00052 parent_name(0), daughters_name(0),
00053 particletable(0),
00054 parent(0), daughters(0),
00055 parent_mass(0.0), daughters_mass(0),
00056 verboseLevel(1)
00057 {
00058
00059 particletable = G4ParticleTable::GetParticleTable();
00060 }
00061
00062
00063 G4VDecayChannel::G4VDecayChannel(const G4String &aName, G4int Verbose)
00064 :kinematics_name(aName),
00065 rbranch(0.0),
00066 numberOfDaughters(0),
00067 parent_name(0), daughters_name(0),
00068 particletable(0),
00069 parent(0), daughters(0),
00070 parent_mass(0.0), daughters_mass(0),
00071 verboseLevel(Verbose)
00072 {
00073
00074 particletable = G4ParticleTable::GetParticleTable();
00075 }
00076
00077 G4VDecayChannel::G4VDecayChannel(const G4String &aName,
00078 const G4String& theParentName,
00079 G4double theBR,
00080 G4int theNumberOfDaughters,
00081 const G4String& theDaughterName1,
00082 const G4String& theDaughterName2,
00083 const G4String& theDaughterName3,
00084 const G4String& theDaughterName4 )
00085 :kinematics_name(aName),
00086 rbranch(theBR),
00087 numberOfDaughters(theNumberOfDaughters),
00088 parent_name(0), daughters_name(0),
00089 particletable(0),
00090 parent(0), daughters(0),
00091 parent_mass(0.0), daughters_mass(0),
00092 verboseLevel(1)
00093 {
00094
00095 particletable = G4ParticleTable::GetParticleTable();
00096
00097
00098 parent_name = new G4String(theParentName);
00099
00100
00101 daughters_name = new G4String*[numberOfDaughters];
00102 for (G4int index=0;index<numberOfDaughters;index++) daughters_name[index]=0;
00103
00104
00105 if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
00106 if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
00107 if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
00108 if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
00109 }
00110
00111
00112
00113 G4VDecayChannel::G4VDecayChannel(const G4VDecayChannel &right)
00114 {
00115 kinematics_name = right.kinematics_name;
00116 verboseLevel = right.verboseLevel;
00117 rbranch = right.rbranch;
00118
00119
00120 parent_name = new G4String(*right.parent_name);
00121 parent = 0;
00122 parent_mass = 0.0;
00123
00124
00125 numberOfDaughters = right.numberOfDaughters;
00126
00127 daughters_name =0;
00128 if ( numberOfDaughters >0 ) {
00129 daughters_name = new G4String*[numberOfDaughters];
00130
00131 for (G4int index=0; index < numberOfDaughters; index++)
00132 {
00133 daughters_name[index] = new G4String(*right.daughters_name[index]);
00134 }
00135 }
00136
00137
00138 daughters_mass = 0;
00139 daughters = 0;
00140
00141
00142 particletable = G4ParticleTable::GetParticleTable();
00143 }
00144
00145 G4VDecayChannel & G4VDecayChannel::operator=(const G4VDecayChannel &right)
00146 {
00147 if (this != &right) {
00148 kinematics_name = right.kinematics_name;
00149 verboseLevel = right.verboseLevel;
00150 rbranch = right.rbranch;
00151
00152
00153 parent_name = new G4String(*right.parent_name);
00154
00155
00156 ClearDaughtersName();
00157
00158
00159 numberOfDaughters = right.numberOfDaughters;
00160 if ( numberOfDaughters >0 ) {
00161 if (daughters_name !=0) ClearDaughtersName();
00162 daughters_name = new G4String*[numberOfDaughters];
00163
00164 for (G4int index=0; index < numberOfDaughters; index++) {
00165 daughters_name[index] = new G4String(*right.daughters_name[index]);
00166 }
00167 }
00168 }
00169
00170
00171 parent = 0;
00172 daughters = 0;
00173 parent_mass = 0.0;
00174 daughters_mass = 0;
00175
00176
00177 particletable = G4ParticleTable::GetParticleTable();
00178
00179 return *this;
00180 }
00181
00182
00183 G4VDecayChannel::~G4VDecayChannel()
00184 {
00185 ClearDaughtersName();
00186 if (parent_name != 0) delete parent_name;
00187 parent_name = 0;
00188 if (daughters_mass != 0) delete [] daughters_mass;
00189 daughters_mass =0;
00190 }
00191
00192 void G4VDecayChannel::ClearDaughtersName()
00193 {
00194 if ( daughters_name != 0) {
00195 if (numberOfDaughters>0) {
00196 #ifdef G4VERBOSE
00197 if (verboseLevel>1) {
00198 G4cerr << "G4VDecayChannel::ClearDaughtersName "
00199 << " for " << *parent_name << G4endl;
00200 }
00201 #endif
00202 for (G4int index=0; index < numberOfDaughters; index++) {
00203 if (daughters_name[index] != 0) delete daughters_name[index];
00204 }
00205 }
00206 delete [] daughters_name;
00207 daughters_name = 0;
00208 }
00209
00210 if (daughters != 0) delete [] daughters;
00211 if (daughters_mass != 0) delete [] daughters_mass;
00212 daughters = 0;
00213 daughters_mass = 0;
00214
00215 numberOfDaughters = 0;
00216 }
00217
00218 void G4VDecayChannel::SetNumberOfDaughters(G4int size)
00219 {
00220 if (size >0) {
00221
00222 ClearDaughtersName();
00223
00224 daughters_name = new G4String*[size];
00225 for (G4int index=0;index<size;index++) daughters_name[index]=0;
00226 numberOfDaughters = size;
00227 }
00228 }
00229
00230 void G4VDecayChannel::SetDaughter(G4int anIndex,
00231 const G4String &particle_name)
00232 {
00233
00234 if (numberOfDaughters<=0) {
00235 #ifdef G4VERBOSE
00236 if (verboseLevel>0) {
00237 G4cerr << "G4VDecayChannel::SetDaughter: "
00238 << "Number of daughters is not defined" << G4endl;
00239 }
00240 #endif
00241 return;
00242 }
00243
00244
00245 if (daughters_name == 0) {
00246
00247 daughters_name = new G4String*[numberOfDaughters];
00248 for (G4int index=0;index<numberOfDaughters;index++) {
00249 daughters_name[index]=0;
00250 }
00251 }
00252
00253
00254 if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
00255 #ifdef G4VERBOSE
00256 if (verboseLevel>0) {
00257 G4cerr << "G4VDecayChannel::SetDaughter"
00258 << "index out of range " << anIndex << G4endl;
00259 }
00260 #endif
00261 } else {
00262
00263 if (daughters_name[anIndex]!=0) delete daughters_name[anIndex];
00264
00265 daughters_name[anIndex] = new G4String(particle_name);
00266
00267 if (daughters != 0) FillDaughters();
00268 #ifdef G4VERBOSE
00269 if (verboseLevel>1) {
00270 G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
00271 G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
00272 }
00273 #endif
00274 }
00275 }
00276
00277 void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
00278 {
00279 if (parent_type != 0) SetDaughter(anIndex, parent_type->GetParticleName());
00280 }
00281
00282 void G4VDecayChannel::FillDaughters()
00283 {
00284 G4int index;
00285
00286 #ifdef G4VERBOSE
00287 if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
00288 #endif
00289 if (daughters != 0) {
00290 delete [] daughters;
00291 daughters = 0;
00292 }
00293
00294
00295 if (parent == 0) FillParent();
00296 G4double parentmass = parent->GetPDGMass();
00297
00298
00299 G4double sumofdaughtermass = 0.0;
00300 if ((numberOfDaughters <=0) || (daughters_name == 0) ){
00301 #ifdef G4VERBOSE
00302 if (verboseLevel>0) {
00303 G4cerr << "G4VDecayChannel::FillDaughters "
00304 << "[ " << parent->GetParticleName() << " ]"
00305 << "numberOfDaughters is not defined yet";
00306 }
00307 #endif
00308 daughters = 0;
00309 G4Exception("G4VDecayChannel::FillDaughters",
00310 "PART011", FatalException,
00311 "Can not fill daughters: numberOfDaughters is not defined yet");
00312 }
00313
00314
00315 daughters = new G4ParticleDefinition*[numberOfDaughters];
00316 if (daughters_mass != 0) delete [] daughters_mass;
00317 daughters_mass = new G4double[numberOfDaughters];
00318
00319 for (index=0; index < numberOfDaughters; index++) {
00320 if (daughters_name[index] == 0) {
00321
00322 #ifdef G4VERBOSE
00323 if (verboseLevel>0) {
00324 G4cerr << "G4VDecayChannel::FillDaughters "
00325 << "[ " << parent->GetParticleName() << " ]"
00326 << index << "-th daughter is not defined yet" << G4endl;
00327 }
00328 #endif
00329 daughters[index] = 0;
00330 G4Exception("G4VDecayChannel::FillDaughters",
00331 "PART011", FatalException,
00332 "Can not fill daughters: name of a daughter is not defined yet");
00333 }
00334
00335 daughters[index] = particletable->FindParticle(*daughters_name[index]);
00336 if (daughters[index] == 0) {
00337
00338 #ifdef G4VERBOSE
00339 if (verboseLevel>0) {
00340 G4cerr << "G4VDecayChannel::FillDaughters "
00341 << "[ " << parent->GetParticleName() << " ]"
00342 << index << ":" << *daughters_name[index]
00343 << " is not defined !!" << G4endl;
00344 G4cerr << " The BR of this decay mode is set to zero " << G4endl;
00345 }
00346 #endif
00347 SetBR(0.0);
00348 return;
00349 }
00350 #ifdef G4VERBOSE
00351 if (verboseLevel>1) {
00352 G4cout << index << ":" << *daughters_name[index];
00353 G4cout << ":" << daughters[index] << G4endl;
00354 }
00355 #endif
00356 daughters_mass[index] = daughters[index]->GetPDGMass();
00357 sumofdaughtermass += daughters[index]->GetPDGMass();
00358 }
00359
00360
00361 G4double widthMass = parent->GetPDGWidth();
00362 if ( (parent->GetParticleType() != "nucleus") &&
00363 (sumofdaughtermass > parentmass + 5*widthMass) ){
00364
00365 #ifdef G4VERBOSE
00366 if (GetVerboseLevel()>0) {
00367 G4cerr << "G4VDecayChannel::FillDaughters "
00368 << "[ " << parent->GetParticleName() << " ]"
00369 << " Energy/Momentum conserevation breaks " <<G4endl;
00370 if (GetVerboseLevel()>1) {
00371 G4cerr << " parent:" << *parent_name
00372 << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
00373 for (index=0; index < numberOfDaughters; index++){
00374 G4cerr << " daughter " << index << ":" << *daughters_name[index]
00375 << " mass:" << daughters[index]->GetPDGMass()/GeV
00376 << "[GeV/c/c]" <<G4endl;
00377 }
00378 }
00379 }
00380 #endif
00381 }
00382 }
00383
00384
00385 void G4VDecayChannel::FillParent()
00386 {
00387 if (parent_name == 0) {
00388
00389 #ifdef G4VERBOSE
00390 if (verboseLevel>0) {
00391 G4cerr << "G4VDecayChannel::FillParent "
00392 << ": parent name is not defined !!" << G4endl;
00393 }
00394 #endif
00395 parent = 0;
00396 G4Exception("G4VDecayChannel::FillParent()",
00397 "PART012", FatalException,
00398 "Can not fill parent: parent name is not defined yet");
00399 }
00400
00401 parent = particletable->FindParticle(*parent_name);
00402 if (parent == 0) {
00403
00404 #ifdef G4VERBOSE
00405 if (verboseLevel>0) {
00406 G4cerr << "G4VDecayChannel::FillParent "
00407 << *parent_name << " does not exist !!" << G4endl;
00408 }
00409 #endif
00410 G4Exception("G4VDecayChannel::FillParent()",
00411 "PART012", FatalException,
00412 "Can not fill parent: parent does not exist");
00413 }
00414 parent_mass = parent->GetPDGMass();
00415 }
00416
00417 void G4VDecayChannel::SetParent(const G4ParticleDefinition * parent_type)
00418 {
00419 if (parent_type != 0) SetParent(parent_type->GetParticleName());
00420 }
00421
00422 G4int G4VDecayChannel::GetAngularMomentum()
00423 {
00424
00425
00426
00427 if (daughters == 0) FillDaughters();
00428
00429 const G4int PiSpin = parent->GetPDGiSpin();
00430 const G4int PParity = parent->GetPDGiParity();
00431 if (2==numberOfDaughters) {
00432 const G4int D1iSpin = daughters[0]->GetPDGiSpin();
00433 const G4int D1Parity = daughters[0]->GetPDGiParity();
00434 const G4int D2iSpin = daughters[1]->GetPDGiSpin();
00435 const G4int D2Parity = daughters[1]->GetPDGiParity();
00436 const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
00437 const G4int MaxiSpin = D1iSpin + D2iSpin;
00438 const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2;
00439 G4int lMin;
00440 #ifdef G4VERBOSE
00441 if (verboseLevel>1) {
00442 G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
00443 G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
00444 }
00445 #endif
00446 for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){
00447 lMin = std::abs(PiSpin-j)/2;
00448 #ifdef G4VERBOSE
00449 if (verboseLevel>1)
00450 G4cout << "-> checking 2*j=" << j << G4endl;
00451 #endif
00452 for (G4int l=lMin; l<=lMax; l++) {
00453 #ifdef G4VERBOSE
00454 if (verboseLevel>1)
00455 G4cout << " checking l=" << l << G4endl;
00456 #endif
00457 if (l%2==0) {
00458 if (PParity == D1Parity*D2Parity) {
00459 return l;
00460 }
00461 } else {
00462 if (PParity == -1*D1Parity*D2Parity) {
00463 return l;
00464 }
00465 }
00466 }
00467 }
00468 } else {
00469 G4Exception("G4VDecayChannel::GetAngularMomentum",
00470 "PART111", JustWarning,
00471 "Sorry, can't handle 3 particle decays (up to now)");
00472 return 0;
00473 }
00474 G4Exception ("G4VDecayChannel::GetAngularMomentum",
00475 "PART111", JustWarning,
00476 "Can't find angular momentum for this decay");
00477 return 0;
00478 }
00479
00480 void G4VDecayChannel::DumpInfo()
00481 {
00482 G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
00483 G4cout << " : " ;
00484 for (G4int index=0; index < numberOfDaughters; index++)
00485 {
00486 if(daughters_name[index] != 0) {
00487 G4cout << " " << *(daughters_name[index]);
00488 } else {
00489 G4cout << " not defined ";
00490 }
00491 }
00492 G4cout << G4endl;
00493 }
00494
00495 const G4String& G4VDecayChannel::GetNoName() const
00496 {
00497 return noName;
00498 }