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