Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // G4VDecayChannel 27 // 28 // Author: H.Kurashige, 27 July 1996 29 // ------------------------------------------- 30 31 #include "G4VDecayChannel.hh" 32 33 #include "G4AutoLock.hh" 34 #include "G4DecayProducts.hh" 35 #include "G4DecayTable.hh" 36 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleTable.hh" 38 #include "G4SystemOfUnits.hh" 39 #include "Randomize.hh" 40 41 const G4String G4VDecayChannel::noName = " "; 42 43 G4VDecayChannel::G4VDecayChannel() 44 { 45 // set pointer to G4ParticleTable (static an 46 particletable = G4ParticleTable::GetParticle 47 } 48 49 G4VDecayChannel::G4VDecayChannel(const G4Strin 50 : kinematics_name(aName), verboseLevel(verbo 51 { 52 // set pointer to G4ParticleTable (static an 53 particletable = G4ParticleTable::GetParticle 54 } 55 56 G4VDecayChannel::G4VDecayChannel(const G4Strin 57 G4double theB 58 const G4Strin 59 const G4Strin 60 const G4Strin 61 : kinematics_name(aName), rbranch(theBR), nu 62 { 63 // set pointer to G4ParticleTable (static an 64 particletable = G4ParticleTable::GetParticle 65 66 // parent name 67 parent_name = new G4String(theParentName); 68 69 // cleate array 70 daughters_name = new G4String*[numberOfDaugh 71 for (G4int index = 0; index < numberOfDaught 72 daughters_name[index] = nullptr; 73 } 74 75 // daughters' name 76 if (numberOfDaughters > 0) daughters_name[0] 77 if (numberOfDaughters > 1) daughters_name[1] 78 if (numberOfDaughters > 2) daughters_name[2] 79 if (numberOfDaughters > 3) daughters_name[3] 80 if (numberOfDaughters > 4) daughters_name[4] 81 82 if (rbranch < 0.) 83 rbranch = 0.0; 84 else if (rbranch > 1.0) 85 rbranch = 1.0; 86 } 87 88 G4VDecayChannel::G4VDecayChannel(const G4VDeca 89 { 90 kinematics_name = right.kinematics_name; 91 verboseLevel = right.verboseLevel; 92 rbranch = right.rbranch; 93 rangeMass = right.rangeMass; 94 95 // copy parent name 96 parent_name = new G4String(*right.parent_nam 97 98 // create array 99 numberOfDaughters = right.numberOfDaughters; 100 101 daughters_name = nullptr; 102 if (numberOfDaughters > 0) { 103 daughters_name = new G4String*[numberOfDau 104 // copy daughters name 105 for (G4int index = 0; index < numberOfDaug 106 daughters_name[index] = new G4String(*ri 107 } 108 } 109 110 // particle table 111 particletable = G4ParticleTable::GetParticle 112 113 parent_polarization = right.parent_polarizat 114 115 G4MUTEXINIT(daughtersMutex); 116 G4MUTEXINIT(parentMutex); 117 } 118 119 G4VDecayChannel& G4VDecayChannel::operator=(co 120 { 121 if (this != &right) { 122 kinematics_name = right.kinematics_name; 123 verboseLevel = right.verboseLevel; 124 rbranch = right.rbranch; 125 rangeMass = right.rangeMass; 126 parent_polarization = right.parent_polariz 127 // copy parent name 128 delete parent_name; 129 parent_name = new G4String(*right.parent_n 130 131 // clear daughters_name array 132 ClearDaughtersName(); 133 134 // recreate array 135 numberOfDaughters = right.numberOfDaughter 136 if (numberOfDaughters > 0) { 137 daughters_name = new G4String*[numberOfD 138 // copy daughters name 139 for (G4int index = 0; index < numberOfDa 140 daughters_name[index] = new G4String(* 141 } 142 } 143 } 144 145 // particle table 146 particletable = G4ParticleTable::GetParticle 147 148 G4MUTEXINIT(daughtersMutex); 149 G4MUTEXINIT(parentMutex); 150 151 return *this; 152 } 153 154 G4VDecayChannel::~G4VDecayChannel() 155 { 156 ClearDaughtersName(); 157 delete parent_name; 158 parent_name = nullptr; 159 delete[] G4MT_daughters_mass; 160 G4MT_daughters_mass = nullptr; 161 delete[] G4MT_daughters_width; 162 G4MT_daughters_width = nullptr; 163 G4MUTEXDESTROY(daughtersMutex); 164 G4MUTEXDESTROY(parentMutex); 165 } 166 167 void G4VDecayChannel::ClearDaughtersName() 168 { 169 G4AutoLock l(&daughtersMutex); 170 if (daughters_name != nullptr) { 171 if (numberOfDaughters > 0) { 172 #ifdef G4VERBOSE 173 if (verboseLevel > 1) { 174 G4cout << "G4VDecayChannel::ClearDaugh 175 << " for " << *parent_name << G 176 } 177 #endif 178 for (G4int index = 0; index < numberOfDa 179 delete daughters_name[index]; 180 } 181 } 182 delete[] daughters_name; 183 daughters_name = nullptr; 184 } 185 186 delete[] G4MT_daughters; 187 delete[] G4MT_daughters_mass; 188 delete[] G4MT_daughters_width; 189 G4MT_daughters_width = nullptr; 190 G4MT_daughters = nullptr; 191 G4MT_daughters_mass = nullptr; 192 193 numberOfDaughters = 0; 194 } 195 196 void G4VDecayChannel::SetNumberOfDaughters(G4i 197 { 198 if (size > 0) { 199 // remove old contents 200 ClearDaughtersName(); 201 // cleate array 202 daughters_name = new G4String*[size]; 203 for (G4int index = 0; index < size; ++inde 204 daughters_name[index] = nullptr; 205 } 206 numberOfDaughters = size; 207 } 208 } 209 210 void G4VDecayChannel::SetDaughter(G4int anInde 211 { 212 // check numberOfDaughters is positive 213 if (numberOfDaughters <= 0) { 214 #ifdef G4VERBOSE 215 if (verboseLevel > 0) { 216 G4cout << "G4VDecayChannel::SetDaughter( 217 << "Number of daughters is not de 218 } 219 #endif 220 return; 221 } 222 223 // An analysis of this code, shows that this 224 // only in the constructor of derived classe 225 // The general idea of this method is probab 226 // the possibility to re-define daughters on 227 // this design is extremely problematic for 228 // require (as practically happens) that the 229 // at construction, i.e. when G4MT_daughters 230 // moreover this method can be called only a 231 // has been called (see previous if), in suc 232 // 233 if (daughters_name == nullptr) { 234 G4Exception("G4VDecayChannel::SetDaughter( 235 "Trying to add a daughter with 236 return; 237 } 238 if (G4MT_daughters != nullptr) { 239 G4Exception("G4VDecayChannel::SetDaughter( 240 "Trying to modify a daughter o 241 but decay channel already has 242 return; 243 } 244 245 // check an index 246 if ((anIndex < 0) || (anIndex >= numberOfDau 247 #ifdef G4VERBOSE 248 if (verboseLevel > 0) { 249 G4cout << "G4VDecayChannel::SetDaughter( 250 << "index out of range " << anInd 251 } 252 #endif 253 } 254 else { 255 // fill the name 256 daughters_name[anIndex] = new G4String(par 257 #ifdef G4VERBOSE 258 if (verboseLevel > 1) { 259 G4cout << "G4VDecayChannel::SetDaughter[ 260 G4cout << daughters_name[anIndex] << ":" 261 } 262 #endif 263 } 264 } 265 266 void G4VDecayChannel::SetDaughter(G4int anInde 267 { 268 if (parent_type != nullptr) SetDaughter(anIn 269 } 270 271 void G4VDecayChannel::FillDaughters() 272 { 273 G4AutoLock lock(&daughtersMutex); 274 275 // Double check, check again if another thre 276 // case do not need to do anything 277 if (G4MT_daughters != nullptr) return; 278 279 G4int index; 280 281 #ifdef G4VERBOSE 282 if (verboseLevel > 1) G4cout << "G4VDecayCha 283 #endif 284 if (G4MT_daughters != nullptr) { 285 delete[] G4MT_daughters; 286 G4MT_daughters = nullptr; 287 } 288 289 // parent mass 290 CheckAndFillParent(); 291 G4double parentmass = G4MT_parent->GetPDGMas 292 293 // 294 G4double sumofdaughtermass = 0.0; 295 G4double sumofdaughterwidthsq = 0.0; 296 297 if ((numberOfDaughters <= 0) || (daughters_n 298 #ifdef G4VERBOSE 299 if (verboseLevel > 0) { 300 G4cout << "G4VDecayChannel::FillDaughter 301 << "[ " << G4MT_parent->GetPartic 302 << "numberOfDaughters is not defi 303 } 304 #endif 305 G4MT_daughters = nullptr; 306 G4Exception("G4VDecayChannel::FillDaughter 307 "Cannot fill daughters: number 308 } 309 310 // create and set the array of pointers to d 311 G4MT_daughters = new G4ParticleDefinition*[n 312 delete[] G4MT_daughters_mass; 313 delete[] G4MT_daughters_width; 314 G4MT_daughters_mass = new G4double[numberOfD 315 G4MT_daughters_width = new G4double[numberOf 316 // loop over all daughters 317 for (index = 0; index < numberOfDaughters; + 318 if (daughters_name[index] == nullptr) { 319 // daughter name is not defined 320 #ifdef G4VERBOSE 321 if (verboseLevel > 0) { 322 G4cout << "G4VDecayChannel::FillDaught 323 << "[ " << G4MT_parent->GetPart 324 << "-th daughter is not defined 325 } 326 #endif 327 G4MT_daughters[index] = nullptr; 328 G4Exception("G4VDecayChannel::FillDaught 329 "Cannot fill daughters: name 330 } 331 // search daughter particles in the partic 332 G4MT_daughters[index] = particletable->Fin 333 if (G4MT_daughters[index] == nullptr) { 334 // cannot find the daughter particle 335 #ifdef G4VERBOSE 336 if (verboseLevel > 0) { 337 G4cout << "G4VDecayChannel::FillDaught 338 << "[ " << G4MT_parent->GetPart 339 << *daughters_name[index] << " 340 G4cout << " The BR of this decay mode 341 } 342 #endif 343 SetBR(0.0); 344 return; 345 } 346 #ifdef G4VERBOSE 347 if (verboseLevel > 1) { 348 G4cout << index << ":" << *daughters_nam 349 G4cout << ":" << G4MT_daughters[index] < 350 } 351 #endif 352 G4MT_daughters_mass[index] = G4MT_daughter 353 G4double d_width = G4MT_daughters[index]-> 354 G4MT_daughters_width[index] = d_width; 355 sumofdaughtermass += G4MT_daughters[index] 356 sumofdaughterwidthsq += d_width * d_width; 357 } // end loop over all daughters 358 359 // check sum of daghter mass 360 G4double widthMass = 361 std::sqrt(G4MT_parent->GetPDGWidth() * G4M 362 if ((G4MT_parent->GetParticleType() != "nucl 363 && (sumofdaughtermass > parentmass + ran 364 { 365 // !!! illegal mass !!! 366 #ifdef G4VERBOSE 367 if (GetVerboseLevel() > 0) { 368 G4cout << "G4VDecayChannel::FillDaughter 369 << "[ " << G4MT_parent->GetPartic 370 << " Energy/Momentum consereva 371 if (GetVerboseLevel() > 1) { 372 G4cout << " parent:" << *parent_nam 373 << G4endl; 374 for (index = 0; index < numberOfDaught 375 G4cout << " daughter " << index 376 << " mass:" << G4MT_daughters 377 } 378 } 379 } 380 #endif 381 } 382 } 383 384 void G4VDecayChannel::FillParent() 385 { 386 G4AutoLock lock(&parentMutex); 387 388 // Double check, check again if another thre 389 // case do not need to do anything 390 if (G4MT_parent != nullptr) return; 391 392 if (parent_name == nullptr) { 393 // parent name is not defined 394 #ifdef G4VERBOSE 395 if (verboseLevel > 0) { 396 G4cout << "G4VDecayChannel::FillParent() 397 << "parent name is not defined !! 398 } 399 #endif 400 G4MT_parent = nullptr; 401 G4Exception("G4VDecayChannel::FillParent() 402 "Cannot fill parent: parent na 403 return; 404 } 405 406 // search parent particle in the particle ta 407 G4MT_parent = particletable->FindParticle(*p 408 if (G4MT_parent == nullptr) { 409 // parent particle does not exist 410 #ifdef G4VERBOSE 411 if (verboseLevel > 0) { 412 G4cout << "G4VDecayChannel::FillParent() 413 << G4endl; 414 } 415 #endif 416 G4Exception("G4VDecayChannel::FillParent() 417 "Cannot fill parent: parent do 418 return; 419 } 420 G4MT_parent_mass = G4MT_parent->GetPDGMass() 421 } 422 423 void G4VDecayChannel::SetParent(const G4Partic 424 { 425 if (parent_type != nullptr) SetParent(parent 426 } 427 428 G4int G4VDecayChannel::GetAngularMomentum() 429 { 430 // determine angular momentum 431 432 // fill pointers to daughter particles if no 433 CheckAndFillDaughters(); 434 435 const G4int PiSpin = G4MT_parent->GetPDGiSpi 436 const G4int PParity = G4MT_parent->GetPDGiPa 437 if (2 == numberOfDaughters) // up to now we 438 { 439 const G4int D1iSpin = G4MT_daughters[0]->G 440 const G4int D1Parity = G4MT_daughters[0]-> 441 const G4int D2iSpin = G4MT_daughters[1]->G 442 const G4int D2Parity = G4MT_daughters[1]-> 443 const G4int MiniSpin = std::abs(D1iSpin - 444 const G4int MaxiSpin = D1iSpin + D2iSpin; 445 const G4int lMax = (PiSpin + D1iSpin + D2i 446 G4int lMin; 447 #ifdef G4VERBOSE 448 if (verboseLevel > 1) { 449 G4cout << "iSpin: " << PiSpin << " -> " 450 G4cout << "2*jmin, 2*jmax, lmax " << Min 451 } 452 #endif 453 for (G4int j = MiniSpin; j <= MaxiSpin; j 454 { // spin couplings 455 lMin = std::abs(PiSpin - j) / 2; 456 #ifdef G4VERBOSE 457 if (verboseLevel > 1) G4cout << "-> chec 458 #endif 459 for (G4int l = lMin; l <= lMax; ++l) { 460 #ifdef G4VERBOSE 461 if (verboseLevel > 1) G4cout << " chec 462 #endif 463 if (l % 2 == 0) { 464 if (PParity == D1Parity * D2Parity) 465 return l; 466 } 467 else { 468 if (PParity == -1 * D1Parity * D2Par 469 return l; 470 } 471 } 472 } 473 } 474 else { 475 G4Exception("G4VDecayChannel::GetAngularMo 476 "Sorry, can't handle 3 particl 477 return 0; 478 } 479 G4Exception("G4VDecayChannel::GetAngularMome 480 "Can't find angular momentum for 481 return 0; 482 } 483 484 void G4VDecayChannel::DumpInfo() 485 { 486 G4cout << " BR: " << rbranch << " [" << ki 487 G4cout << " : "; 488 for (G4int index = 0; index < numberOfDaught 489 if (daughters_name[index] != nullptr) { 490 G4cout << " " << *(daughters_name[index] 491 } 492 else { 493 G4cout << " not defined "; 494 } 495 } 496 G4cout << G4endl; 497 } 498 499 const G4String& G4VDecayChannel::GetNoName() c 500 { 501 return noName; 502 } 503 504 G4double G4VDecayChannel::DynamicalMass(G4doub 505 { 506 if (width <= 0.0) return massPDG; 507 if (maxDev > rangeMass) maxDev = rangeMass; 508 if (maxDev <= -1. * rangeMass) return massPD 509 510 G4double x = G4UniformRand() * (maxDev + ran 511 G4double y = G4UniformRand(); 512 const std::size_t MAX_LOOP = 10000; 513 for (std::size_t loop_counter = 0; loop_coun 514 if (y * (width * width * x * x + massPDG * 515 <= massPDG * massPDG * width * width) 516 break; 517 x = G4UniformRand() * (maxDev + rangeMass) 518 y = G4UniformRand(); 519 } 520 G4double mass = massPDG + x * width; 521 return mass; 522 } 523 524 G4bool G4VDecayChannel::IsOKWithParentMass(G4d 525 { 526 G4double sumOfDaughterMassMin = 0.0; 527 CheckAndFillParent(); 528 CheckAndFillDaughters(); 529 // skip one body decay 530 if (numberOfDaughters == 1) return true; 531 532 for (G4int index = 0; index < numberOfDaught 533 sumOfDaughterMassMin += G4MT_daughters_mas 534 } 535 return (parentMass >= sumOfDaughterMassMin); 536 } 537 538 void G4VDecayChannel::SetBR(G4double value) 539 { 540 rbranch = value; 541 if (rbranch < 0.) 542 rbranch = 0.0; 543 else if (rbranch > 1.0) 544 rbranch = 1.0; 545 } 546