Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4PhaseSpaceDecayChannel class implementati << 27 // 23 // 28 // Author: H.Kurashige, 27 July 1996 << 24 // $Id: G4PhaseSpaceDecayChannel.cc,v 1.7 2001/10/15 09:58:36 kurasige Exp $ 29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-05-02-patch-01 $ 30 << 26 // 31 #include "G4PhaseSpaceDecayChannel.hh" << 27 // >> 28 // ------------------------------------------------------------ >> 29 // GEANT 4 class header file >> 30 // >> 31 // History: first implementation, based on object model of >> 32 // 27 July 1996 H.Kurashige >> 33 // 20 Oct 1996 H.Kurashige >> 34 // 30 May 1997 H.Kurashige >> 35 // ------------------------------------------------------------ 32 36 33 #include "G4DecayProducts.hh" << 34 #include "G4LorentzRotation.hh" << 35 #include "G4LorentzVector.hh" << 36 #include "G4ParticleDefinition.hh" 37 #include "G4ParticleDefinition.hh" 37 #include "G4PhysicalConstants.hh" << 38 #include "G4DecayProducts.hh" 38 #include "G4SystemOfUnits.hh" << 39 #include "G4VDecayChannel.hh" 39 #include "G4VDecayChannel.hh" >> 40 #include "G4PhaseSpaceDecayChannel.hh" 40 #include "Randomize.hh" 41 #include "Randomize.hh" >> 42 #include "G4LorentzVector.hh" >> 43 #include "G4LorentzRotation.hh" >> 44 41 45 42 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayCha 46 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose) 43 : G4VDecayChannel("Phase Space", Verbose) << 47 :G4VDecayChannel("Phase Space", Verbose) 44 {} << 48 { 45 49 46 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayCha << 50 } 47 << 48 << 49 << 50 << 51 << 52 << 53 : G4VDecayChannel("Phase Space", theParentNa << 54 theDaughterName2, theDaugh << 55 {} << 56 51 57 G4DecayProducts* G4PhaseSpaceDecayChannel::Dec << 52 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel( >> 53 const G4String& theParentName, >> 54 G4double theBR, >> 55 G4int theNumberOfDaughters, >> 56 const G4String& theDaughterName1, >> 57 const G4String& theDaughterName2, >> 58 const G4String& theDaughterName3, >> 59 const G4String& theDaughterName4 ) >> 60 :G4VDecayChannel("Phase Space", >> 61 theParentName,theBR, >> 62 theNumberOfDaughters, >> 63 theDaughterName1, >> 64 theDaughterName2, >> 65 theDaughterName3, >> 66 theDaughterName4) 58 { 67 { 59 #ifdef G4VERBOSE << 60 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 61 #endif << 62 68 63 G4DecayProducts* products = nullptr; << 69 } 64 << 65 CheckAndFillParent(); << 66 CheckAndFillDaughters(); << 67 70 68 if (parentMass > 0.0) << 71 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel() 69 current_parent_mass.Put(parentMass); << 72 { 70 else << 73 } 71 current_parent_mass.Put(G4MT_parent_mass); << 72 74 73 switch (numberOfDaughters) { << 75 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass) 74 case 0: << 76 { 75 #ifdef G4VERBOSE 77 #ifdef G4VERBOSE 76 if (GetVerboseLevel() > 0) { << 78 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt "; 77 G4cout << "G4PhaseSpaceDecayChannel::D << 78 G4cout << " daughters not defined " << << 79 } << 80 #endif 79 #endif 81 break; << 80 82 case 1: << 81 G4DecayProducts * products = 0; 83 products = OneBodyDecayIt(); << 82 84 break; << 83 if (parent == 0) FillParent(); 85 case 2: << 84 if (daughters == 0) FillDaughters(); 86 products = TwoBodyDecayIt(); << 85 87 break; << 86 if (parentMass >0.0) parent_mass = parentMass; 88 case 3: << 87 switch (numberOfDaughters){ 89 products = ThreeBodyDecayIt(); << 88 case 0: 90 break; << 89 #ifdef G4VERBOSE 91 default: << 90 if (GetVerboseLevel()>0) { 92 products = ManyBodyDecayIt(); << 91 G4cout << "G4PhaseSpaceDecayChannel::DecayIt "; 93 break; << 92 G4cout << " daughters not defined " <<G4endl; 94 } << 93 } 95 #ifdef G4VERBOSE << 94 #endif 96 if ((products == nullptr) && (GetVerboseLeve << 95 break; 97 G4cout << "G4PhaseSpaceDecayChannel::Decay << 96 case 1: 98 G4cout << *parent_name << " cannot decay " << 97 products = OneBodyDecayIt(); >> 98 break; >> 99 case 2: >> 100 products = TwoBodyDecayIt(); >> 101 break; >> 102 case 3: >> 103 products = ThreeBodyDecayIt(); >> 104 break; >> 105 default: >> 106 products = ManyBodyDecayIt(); >> 107 break; >> 108 } >> 109 #ifdef G4VERBOSE >> 110 if ((products == 0) && (GetVerboseLevel()>0)) { >> 111 G4cout << "G4PhaseSpaceDecayChannel::DecayIt "; >> 112 G4cout << *parent_name << " can not decay " << G4endl; 99 DumpInfo(); 113 DumpInfo(); 100 } 114 } 101 #endif 115 #endif 102 return products; 116 return products; 103 } 117 } 104 118 105 G4DecayProducts* G4PhaseSpaceDecayChannel::One << 119 G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt() 106 { 120 { 107 #ifdef G4VERBOSE 121 #ifdef G4VERBOSE 108 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 122 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl; 109 #endif 123 #endif 110 // parent mass << 111 G4double parentmass = current_parent_mass.Ge << 112 124 113 // create parent G4DynamicParticle at rest << 125 //create parent G4DynamicParticle at rest 114 G4ThreeVector dummy; 126 G4ThreeVector dummy; 115 auto parentparticle = new G4DynamicParticle( << 127 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 116 // create G4Decayproducts << 128 //create G4Decayproducts 117 auto products = new G4DecayProducts(*parentp << 129 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 118 delete parentparticle; 130 delete parentparticle; 119 131 120 // create daughter G4DynamicParticle at rest << 132 //create daughter G4DynamicParticle at rest 121 auto daughterparticle = new G4DynamicParticl << 133 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0); 122 if (useGivenDaughterMass) daughterparticle-> << 123 products->PushProducts(daughterparticle); 134 products->PushProducts(daughterparticle); 124 135 125 #ifdef G4VERBOSE 136 #ifdef G4VERBOSE 126 if (GetVerboseLevel() > 1) { << 137 if (GetVerboseLevel()>1) { 127 G4cout << "G4PhaseSpaceDecayChannel::OneBo << 138 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt "; 128 G4cout << " create decay products in rest << 139 G4cout << " create decay products in rest frame " <<G4endl; 129 products->DumpInfo(); << 140 products->DumpInfo(); 130 } 141 } 131 #endif 142 #endif 132 return products; 143 return products; 133 } 144 } 134 145 135 G4DecayProducts* G4PhaseSpaceDecayChannel::Two << 146 G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt() 136 { 147 { 137 #ifdef G4VERBOSE 148 #ifdef G4VERBOSE 138 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 149 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl; 139 #endif 150 #endif 140 // parent mass 151 // parent mass 141 G4double parentmass = current_parent_mass.Ge << 152 G4double parentmass = parent_mass; >> 153 >> 154 //daughters'mass >> 155 G4double daughtermass[2]; >> 156 G4double daughtermomentum; >> 157 daughtermass[0] = daughters_mass[0]; >> 158 daughtermass[1] = daughters_mass[1]; 142 159 143 // daughters'mass, width << 160 //create parent G4DynamicParticle at rest 144 G4double daughtermass[2], daughterwidth[2]; << 145 daughtermass[0] = G4MT_daughters_mass[0]; << 146 daughtermass[1] = G4MT_daughters_mass[1]; << 147 daughterwidth[0] = G4MT_daughters_width[0]; << 148 daughterwidth[1] = G4MT_daughters_width[1]; << 149 << 150 // create parent G4DynamicParticle at rest << 151 G4ThreeVector dummy; 161 G4ThreeVector dummy; 152 auto parentparticle = new G4DynamicParticle( << 162 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 153 // create G4Decayproducts << 163 //create G4Decayproducts 154 auto products = new G4DecayProducts(*parentp << 164 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 155 delete parentparticle; 165 delete parentparticle; 156 166 157 if (!useGivenDaughterMass) { << 167 //calculate daughter momentum 158 G4bool withWidth = (daughterwidth[0] > 1.0 << 168 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]); 159 || (daughterwidth[1] > << 169 G4double costheta = 2.*G4UniformRand()-1.0; 160 if (withWidth) { << 170 G4double sintheta = sqrt((1.0 - costheta)*(1.0 + costheta)); 161 G4double sumofdaughterwidthsq = << 171 G4double phi = 2.0*M_PI*G4UniformRand()*rad; 162 daughterwidth[0] * daughterwidth[0] + << 172 G4ThreeVector direction(sintheta*cos(phi),sintheta*sin(phi),costheta); 163 // check parent mass and daughter mass << 173 164 G4double maxDev = << 174 //create daughter G4DynamicParticle 165 (parentmass - daughtermass[0] - daught << 175 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum); 166 if (maxDev <= -1.0 * rangeMass) { << 167 #ifdef G4VERBOSE << 168 if (GetVerboseLevel() > 0) { << 169 G4cout << "G4PhaseSpaceDecayChannel: << 170 << "Sum of daughter mass is l << 171 G4cout << "Parent :" << G4MT_parent- << 172 << current_parent_mass.Get() << 173 G4cout << "Daughter 1 :" << G4MT_dau << 174 << daughtermass[0] / GeV << G << 175 G4cout << "Daughter 2:" << G4MT_daug << 176 << daughtermass[1] / GeV << G << 177 } << 178 #endif << 179 G4Exception("G4PhaseSpaceDecayChannel: << 180 "Cannot create decay produ << 181 larger than parent mass!" << 182 return products; << 183 } << 184 G4double dm1 = daughtermass[0]; << 185 if (daughterwidth[0] > 0.) dm1 = Dynamic << 186 G4double dm2 = daughtermass[1]; << 187 if (daughterwidth[1] > 0.) dm2 = Dynamic << 188 while (dm1 + dm2 > parentmass) // Loop << 189 { << 190 dm1 = DynamicalMass(daughtermass[0], d << 191 dm2 = DynamicalMass(daughtermass[1], d << 192 } << 193 daughtermass[0] = dm1; << 194 daughtermass[1] = dm2; << 195 } << 196 } << 197 else { << 198 // use given daughter mass; << 199 daughtermass[0] = givenDaughterMasses[0]; << 200 daughtermass[1] = givenDaughterMasses[1]; << 201 } << 202 if (parentmass < daughtermass[0] + daughterm << 203 #ifdef G4VERBOSE << 204 if (GetVerboseLevel() > 0) { << 205 G4cout << "G4PhaseSpaceDecayChannel::Two << 206 << "Sum of daughter mass is large << 207 G4cout << "Parent :" << G4MT_parent->Get << 208 << current_parent_mass.Get() / Ge << 209 G4cout << "Daughter 1 :" << G4MT_daughte << 210 << daughtermass[0] / GeV << G4end << 211 G4cout << "Daughter 2:" << G4MT_daughter << 212 << daughtermass[1] / GeV << G4end << 213 if (useGivenDaughterMass) { << 214 G4cout << "Daughter Mass is given." << << 215 } << 216 } << 217 #endif << 218 G4Exception("G4PhaseSpaceDecayChannel::Two << 219 "Cannot create decay products: << 220 larger than parent mass!"); << 221 return products; << 222 } << 223 << 224 // calculate daughter momentum << 225 G4double daughtermomentum = Pmx(parentmass, << 226 << 227 G4double costheta = 2. * G4UniformRand() - 1 << 228 G4double sintheta = std::sqrt((1.0 - costhet << 229 G4double phi = twopi * G4UniformRand() * rad << 230 G4ThreeVector direction(sintheta * std::cos( << 231 << 232 // create daughter G4DynamicParticle << 233 G4double Ekin = std::sqrt(daughtermomentum * << 234 - daughtermass[0]; << 235 auto daughterparticle = << 236 new G4DynamicParticle(G4MT_daughters[0], d << 237 products->PushProducts(daughterparticle); 176 products->PushProducts(daughterparticle); 238 Ekin = std::sqrt(daughtermomentum * daughter << 177 daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum)); 239 - daughtermass[1]; << 240 daughterparticle = << 241 new G4DynamicParticle(G4MT_daughters[1], - << 242 products->PushProducts(daughterparticle); 178 products->PushProducts(daughterparticle); 243 179 244 #ifdef G4VERBOSE 180 #ifdef G4VERBOSE 245 if (GetVerboseLevel() > 1) { << 181 if (GetVerboseLevel()>1) { 246 G4cout << "G4PhaseSpaceDecayChannel::TwoBo << 182 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "; 247 G4cout << " Create decay products in rest << 183 G4cout << " create decay products in rest frame " <<G4endl; 248 products->DumpInfo(); << 184 products->DumpInfo(); 249 } 185 } 250 #endif 186 #endif 251 return products; 187 return products; 252 } 188 } 253 189 254 G4DecayProducts* G4PhaseSpaceDecayChannel::Thr << 190 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt() >> 191 // algorism of this code is originally written in GDECA3 of GEANT3 255 { 192 { 256 // Algorithm of this code originally written << 257 << 258 #ifdef G4VERBOSE 193 #ifdef G4VERBOSE 259 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 194 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl; 260 #endif 195 #endif 261 // parent mass 196 // parent mass 262 G4double parentmass = current_parent_mass.Ge << 197 G4double parentmass = parent_mass; 263 // daughters'mass << 198 //daughters'mass 264 G4double daughtermass[3], daughterwidth[3]; << 199 G4double daughtermass[3]; 265 G4double sumofdaughtermass = 0.0; 200 G4double sumofdaughtermass = 0.0; 266 G4double sumofdaughterwidthsq = 0.0; << 201 for (G4int index=0; index<3; index++){ 267 G4bool withWidth = false; << 202 daughtermass[index] = daughters_mass[index]; 268 for (G4int index = 0; index < 3; ++index) { << 269 daughtermass[index] = G4MT_daughters_mass[ << 270 sumofdaughtermass += daughtermass[index]; 203 sumofdaughtermass += daughtermass[index]; 271 daughterwidth[index] = G4MT_daughters_widt << 272 sumofdaughterwidthsq += daughterwidth[inde << 273 withWidth = withWidth || (daughterwidth[in << 274 } 204 } 275 << 205 276 // create parent G4DynamicParticle at rest << 206 //create parent G4DynamicParticle at rest 277 G4ThreeVector dummy; 207 G4ThreeVector dummy; 278 auto parentparticle = new G4DynamicParticle( << 208 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 279 << 209 //create G4Decayproducts 280 // create G4Decayproducts << 210 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 281 auto products = new G4DecayProducts(*parentp << 282 delete parentparticle; 211 delete parentparticle; 283 212 284 if (!useGivenDaughterMass) { << 213 //calculate daughter momentum 285 if (withWidth) { << 214 // Generate two 286 G4double maxDev = (parentmass - sumofdau << 287 if (maxDev <= -1.0 * rangeMass) { << 288 #ifdef G4VERBOSE << 289 if (GetVerboseLevel() > 0) { << 290 G4cout << "G4PhaseSpaceDecayChannel: << 291 << "Sum of daughter mass is l << 292 G4cout << "Parent :" << G4MT_parent- << 293 << current_parent_mass.Get() << 294 G4cout << "Daughter 1 :" << G4MT_dau << 295 << daughtermass[0] / GeV << G << 296 G4cout << "Daughter 2:" << G4MT_daug << 297 << daughtermass[1] / GeV << G << 298 G4cout << "Daughter 3:" << G4MT_daug << 299 << daughtermass[2] / GeV << G << 300 } << 301 #endif << 302 G4Exception("G4PhaseSpaceDecayChannel: << 303 "Cannot create decay produ << 304 is larger than parent mas << 305 return products; << 306 } << 307 G4double dm1 = daughtermass[0]; << 308 if (daughterwidth[0] > 0.) dm1 = Dynamic << 309 G4double dm2 = daughtermass[1]; << 310 if (daughterwidth[1] > 0.) dm2 = Dynamic << 311 G4double dm3 = daughtermass[2]; << 312 if (daughterwidth[2] > 0.) dm3 = Dynamic << 313 while (dm1 + dm2 + dm3 > parentmass) // << 314 { << 315 dm1 = DynamicalMass(daughtermass[0], d << 316 dm2 = DynamicalMass(daughtermass[1], d << 317 dm3 = DynamicalMass(daughtermass[2], d << 318 } << 319 daughtermass[0] = dm1; << 320 daughtermass[1] = dm2; << 321 daughtermass[2] = dm3; << 322 sumofdaughtermass = dm1 + dm2 + dm3; << 323 } << 324 } << 325 else { << 326 // use given daughter mass; << 327 daughtermass[0] = givenDaughterMasses[0]; << 328 daughtermass[1] = givenDaughterMasses[1]; << 329 daughtermass[2] = givenDaughterMasses[2]; << 330 sumofdaughtermass = daughtermass[0] + daug << 331 } << 332 << 333 if (sumofdaughtermass > parentmass) { << 334 #ifdef G4VERBOSE << 335 if (GetVerboseLevel() > 0) { << 336 G4cout << "G4PhaseSpaceDecayChannel::Thr << 337 << "Sum of daughter mass is large << 338 G4cout << "Parent :" << G4MT_parent->Get << 339 << current_parent_mass.Get() / Ge << 340 G4cout << "Daughter 1 :" << G4MT_daughte << 341 << daughtermass[0] / GeV << G4end << 342 G4cout << "Daughter 2:" << G4MT_daughter << 343 << daughtermass[1] / GeV << G4end << 344 G4cout << "Daughter 3:" << G4MT_daughter << 345 << daughtermass[2] / GeV << G4end << 346 } << 347 if (useGivenDaughterMass) { << 348 G4cout << "Daughter Mass is given." << G << 349 } << 350 #endif << 351 G4Exception("G4PhaseSpaceDecayChannel::Thr << 352 "Can not create decay products << 353 is larger than parent mass!") << 354 return products; << 355 } << 356 << 357 // calculate daughter momentum << 358 // Generate two << 359 G4double rd1, rd2, rd; 215 G4double rd1, rd2, rd; 360 G4double daughtermomentum[3]; 216 G4double daughtermomentum[3]; 361 G4double momentummax = 0.0, momentumsum = 0. << 217 G4double momentummax=0.0, momentumsum = 0.0; 362 G4double energy; 218 G4double energy; 363 const std::size_t MAX_LOOP = 10000; << 364 219 365 for (std::size_t loop_counter = 0; loop_coun << 220 do { 366 rd1 = G4UniformRand(); 221 rd1 = G4UniformRand(); 367 rd2 = G4UniformRand(); 222 rd2 = G4UniformRand(); 368 if (rd2 > rd1) { 223 if (rd2 > rd1) { 369 rd = rd1; << 224 rd = rd1; 370 rd1 = rd2; 225 rd1 = rd2; 371 rd2 = rd; 226 rd2 = rd; 372 } << 227 } 373 momentummax = 0.0; 228 momentummax = 0.0; 374 momentumsum = 0.0; 229 momentumsum = 0.0; 375 // daughter 0 230 // daughter 0 376 energy = rd2 * (parentmass - sumofdaughter << 231 energy = rd2*(parentmass - sumofdaughtermass); 377 daughtermomentum[0] = std::sqrt(energy * e << 232 daughtermomentum[0] = sqrt(energy*energy + 2.0*energy* daughtermass[0]); 378 if (daughtermomentum[0] > momentummax) mom << 233 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0]; 379 momentumsum += daughtermomentum[0]; << 234 momentumsum += daughtermomentum[0]; 380 // daughter 1 235 // daughter 1 381 energy = (1. - rd1) * (parentmass - sumofd << 236 energy = (1.-rd1)*(parentmass - sumofdaughtermass); 382 daughtermomentum[1] = std::sqrt(energy * e << 237 daughtermomentum[1] = sqrt(energy*energy + 2.0*energy* daughtermass[1]); 383 if (daughtermomentum[1] > momentummax) mom << 238 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1]; 384 momentumsum += daughtermomentum[1]; << 239 momentumsum += daughtermomentum[1]; 385 // daughter 2 240 // daughter 2 386 energy = (rd1 - rd2) * (parentmass - sumof << 241 energy = (rd1-rd2)*(parentmass - sumofdaughtermass); 387 daughtermomentum[2] = std::sqrt(energy * e << 242 daughtermomentum[2] = sqrt(energy*energy + 2.0*energy* daughtermass[2]); 388 if (daughtermomentum[2] > momentummax) mom << 243 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2]; 389 momentumsum += daughtermomentum[2]; << 244 momentumsum += daughtermomentum[2]; 390 if (momentummax <= momentumsum - momentumm << 245 } while (momentummax > momentumsum - momentummax ); 391 } << 392 << 393 // output message 246 // output message 394 #ifdef G4VERBOSE 247 #ifdef G4VERBOSE 395 if (GetVerboseLevel() > 1) { << 248 if (GetVerboseLevel()>1) { 396 G4cout << " daughter 0:" << daughtermo << 249 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl; 397 G4cout << " daughter 1:" << daughtermo << 250 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl; 398 G4cout << " daughter 2:" << daughtermo << 251 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl; 399 G4cout << " momentum sum:" << momentumsu << 252 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl; 400 } 253 } 401 #endif 254 #endif 402 255 403 // create daughter G4DynamicParticle << 256 //create daughter G4DynamicParticle 404 G4double costheta, sintheta, phi, sinphi, co << 257 G4double costheta, sintheta, phi, sinphi, cosphi; 405 G4double costhetan, sinthetan, phin, sinphin << 258 G4double costhetan, sinthetan, phin, sinphin, cosphin; 406 costheta = 2. * G4UniformRand() - 1.0; << 259 costheta = 2.*G4UniformRand()-1.0; 407 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 260 sintheta = sqrt((1.0-costheta)*(1.0+costheta)); 408 phi = twopi * G4UniformRand() * rad; << 261 phi = 2.0*M_PI*G4UniformRand()*rad; 409 sinphi = std::sin(phi); << 262 sinphi = sin(phi); 410 cosphi = std::cos(phi); << 263 cosphi = cos(phi); 411 << 264 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta); 412 G4ThreeVector direction0(sintheta * cosphi, << 265 G4DynamicParticle * daughterparticle 413 G4double Ekin = << 266 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]); 414 std::sqrt(daughtermomentum[0] * daughtermo << 415 - daughtermass[0]; << 416 auto daughterparticle = << 417 new G4DynamicParticle(G4MT_daughters[0], d << 418 products->PushProducts(daughterparticle); 267 products->PushProducts(daughterparticle); 419 268 420 costhetan = (daughtermomentum[1] * daughterm << 269 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]); 421 - daughtermomentum[0] * daughte << 270 sinthetan = sqrt((1.0-costhetan)*(1.0+costhetan)); 422 / (2.0 * daughtermomentum[2] * d << 271 phin = 2.0*M_PI*G4UniformRand()*rad; 423 sinthetan = std::sqrt((1.0 - costhetan) * (1 << 272 sinphin = sin(phin); 424 phin = twopi * G4UniformRand() * rad; << 273 cosphin = cos(phin); 425 sinphin = std::sin(phin); << 426 cosphin = std::cos(phin); << 427 G4ThreeVector direction2; 274 G4ThreeVector direction2; 428 direction2.setX(sinthetan * cosphin * costhe << 275 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 429 + costhetan * sintheta * cos << 276 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 430 direction2.setY(sinthetan * cosphin * costhe << 277 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta); 431 + costhetan * sintheta * sin << 278 daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag())); 432 direction2.setZ(-sinthetan * cosphin * sinth << 433 G4ThreeVector pmom = daughtermomentum[2] * d << 434 Ekin = std::sqrt(pmom.mag2() + daughtermass[ << 435 daughterparticle = << 436 new G4DynamicParticle(G4MT_daughters[2], p << 437 products->PushProducts(daughterparticle); 279 products->PushProducts(daughterparticle); 438 280 439 pmom = (direction0 * daughtermomentum[0] + d << 281 daughterparticle = 440 * (-1.0); << 282 new G4DynamicParticle( 441 Ekin = std::sqrt(pmom.mag2() + daughtermass[ << 283 daughters[1], 442 daughterparticle = << 284 (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0) 443 new G4DynamicParticle(G4MT_daughters[1], p << 285 ); 444 products->PushProducts(daughterparticle); 286 products->PushProducts(daughterparticle); 445 287 446 #ifdef G4VERBOSE 288 #ifdef G4VERBOSE 447 if (GetVerboseLevel() > 1) { << 289 if (GetVerboseLevel()>1) { 448 G4cout << "G4PhaseSpaceDecayChannel::Three << 290 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "; 449 G4cout << " create decay products in rest << 291 G4cout << " create decay products in rest frame " <<G4endl; 450 products->DumpInfo(); << 292 products->DumpInfo(); 451 } 293 } 452 #endif 294 #endif 453 return products; 295 return products; 454 } 296 } 455 297 456 G4DecayProducts* G4PhaseSpaceDecayChannel::Man << 298 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt() >> 299 // algorithm of this code is originally written in FORTRAN by M.Asai >> 300 //****************************************************************** >> 301 // NBODY >> 302 // N-body phase space Monte-Carlo generator >> 303 // Makoto Asai >> 304 // Hiroshima Institute of Technology >> 305 // (asai@kekvax.kek.jp) >> 306 // Revised release : 19/Apr/1995 >> 307 // 457 { 308 { 458 // Algorithm of this code originally written << 459 // ***************************************** << 460 // NBODY - N-body phase space Monte-Carlo ge << 461 // Makoto Asai - Hiroshima Institute of Tech << 462 // Revised release : 19/Apr/1995 << 463 << 464 G4int index, index2; 309 G4int index, index2; 465 310 >> 311 //return value >> 312 G4DecayProducts *products; >> 313 466 #ifdef G4VERBOSE 314 #ifdef G4VERBOSE 467 if (GetVerboseLevel() > 1) G4cout << "G4Phas << 315 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl; 468 #endif 316 #endif 469 << 470 // parent mass 317 // parent mass 471 G4double parentmass = current_parent_mass.Ge << 318 G4double parentmass = parent_mass; 472 << 319 //daughters'mass 473 // parent particle << 320 G4double *daughtermass = new G4double[numberOfDaughters]; 474 G4ThreeVector dummy; << 475 auto parentparticle = new G4DynamicParticle( << 476 << 477 // create G4Decayproducts << 478 auto products = new G4DecayProducts(*parentp << 479 delete parentparticle; << 480 << 481 // daughters'mass << 482 auto daughtermass = new G4double[numberOfDau << 483 << 484 G4double sumofdaughtermass = 0.0; 321 G4double sumofdaughtermass = 0.0; 485 for (index = 0; index < numberOfDaughters; + << 322 for (index=0; index<numberOfDaughters; index++){ 486 if (!useGivenDaughterMass) { << 323 daughtermass[index] = daughters_mass[index]; 487 daughtermass[index] = G4MT_daughters_mas << 488 } << 489 else { << 490 daughtermass[index] = givenDaughterMasse << 491 } << 492 sumofdaughtermass += daughtermass[index]; 324 sumofdaughtermass += daughtermass[index]; 493 } 325 } 494 if (sumofdaughtermass > parentmass) { << 326 495 #ifdef G4VERBOSE << 327 //Calculate daughter momentum 496 if (GetVerboseLevel() > 0) { << 328 G4double *daughtermomentum = new G4double[numberOfDaughters]; 497 G4cout << "G4PhaseSpaceDecayChannel::Man << 329 G4ThreeVector direction; 498 << "Sum of daughter mass is large << 330 G4DynamicParticle **daughterparticle; 499 G4cout << "Parent :" << G4MT_parent->Get << 331 G4double *sm = new G4double[numberOfDaughters]; 500 << current_parent_mass.Get() / Ge << 501 G4cout << "Daughter 1 :" << G4MT_daughte << 502 << daughtermass[0] / GeV << G4end << 503 G4cout << "Daughter 2:" << G4MT_daughter << 504 << daughtermass[1] / GeV << G4end << 505 G4cout << "Daughter 3:" << G4MT_daughter << 506 << daughtermass[2] / GeV << G4end << 507 G4cout << "Daughter 4:" << G4MT_daughter << 508 << daughtermass[3] / GeV << G4end << 509 G4cout << "Daughter 5:" << G4MT_daughter << 510 << daughtermass[4] / GeV << G4end << 511 } << 512 #endif << 513 G4Exception("G4PhaseSpaceDecayChannel::Man << 514 "Can not create decay products << 515 is larger than parent mass!") << 516 delete[] daughtermass; << 517 return products; << 518 } << 519 << 520 // Calculate daughter momentum << 521 auto daughtermomentum = new G4double[numberO << 522 G4ThreeVector direction; << 523 G4DynamicParticle** daughterparticle; << 524 auto sm = new G4double[numberOfDaughters]; << 525 G4double tmas; 332 G4double tmas; 526 G4double weight = 1.0; 333 G4double weight = 1.0; 527 G4int numberOfTry = 0; << 334 G4int numberOfTry = 0; 528 const std::size_t MAX_LOOP = 10000; << 335 529 << 336 do { 530 for (std::size_t loop_counter = 0; loop_coun << 337 //Generate rundom number in descending order 531 // Generate rundom number in descending or << 532 G4double temp; 338 G4double temp; 533 auto rd = new G4double[numberOfDaughters]; << 339 G4double *rd = new G4double[numberOfDaughters]; 534 rd[0] = 1.0; 340 rd[0] = 1.0; 535 for (index = 1; index < numberOfDaughters << 341 for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand(); 536 rd[index] = G4UniformRand(); << 342 rd[ numberOfDaughters -1] = 0.0; 537 } << 343 for(index =1; index < numberOfDaughters -1; index++) { 538 rd[numberOfDaughters - 1] = 0.0; << 344 for(index2 = index+1; index2 < numberOfDaughters; index2++) { 539 for (index = 1; index < numberOfDaughters << 345 if (rd[index] < rd[index2]){ 540 for (index2 = index + 1; index2 < number << 346 temp = rd[index]; 541 if (rd[index] < rd[index2]) { << 347 rd[index] = rd[index2]; 542 temp = rd[index]; << 543 rd[index] = rd[index2]; << 544 rd[index2] = temp; 348 rd[index2] = temp; 545 } 349 } 546 } 350 } 547 } 351 } 548 << 352 549 // calculate virtual mass << 353 //calcurate virtual mass 550 tmas = parentmass - sumofdaughtermass; << 354 tmas = parentmass - sumofdaughtermass; 551 temp = sumofdaughtermass; << 355 temp = sumofdaughtermass; 552 for (index = 0; index < numberOfDaughters; << 356 for(index =0; index < numberOfDaughters; index++) { 553 sm[index] = rd[index] * tmas + temp; << 357 sm[index] = rd[index]*tmas + temp; 554 temp -= daughtermass[index]; 358 temp -= daughtermass[index]; 555 if (GetVerboseLevel() > 1) { << 359 if (GetVerboseLevel()>1) { 556 G4cout << index << " rundom number:" << 360 G4cout << index << " rundom number:" << rd[index]; 557 G4cout << " virtual mass:" << sm[ind << 361 G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl; 558 } 362 } 559 } 363 } 560 delete[] rd; << 364 delete [] rd; 561 365 562 // Calculate daughter momentum << 366 //Calculate daughter momentum 563 weight = 1.0; 367 weight = 1.0; 564 G4bool smOK = true; << 368 index =numberOfDaughters-1; 565 for (index = 0; index < numberOfDaughters << 369 daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]); 566 smOK = (sm[index] - daughtermass[index] << 567 } << 568 if (!smOK) continue; << 569 << 570 index = numberOfDaughters - 1; << 571 daughtermomentum[index] = Pmx(sm[index - 1 << 572 #ifdef G4VERBOSE 370 #ifdef G4VERBOSE 573 if (GetVerboseLevel() > 1) { << 371 if (GetVerboseLevel()>1) { 574 G4cout << " daughter " << index << " 372 G4cout << " daughter " << index << ":" << *daughters_name[index]; 575 G4cout << " momentum:" << daughtermoment << 373 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl; 576 } 374 } 577 #endif 375 #endif 578 for (index = numberOfDaughters - 2; index << 376 for(index =numberOfDaughters-2; index>=0; index--) { 579 // calculate << 377 // calculate 580 daughtermomentum[index] = Pmx(sm[index], << 378 daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]); 581 if (daughtermomentum[index] < 0.0) { << 379 if(daughtermomentum[index] < 0.0) { 582 // !!! illegal momentum !!! 380 // !!! illegal momentum !!! 583 #ifdef G4VERBOSE 381 #ifdef G4VERBOSE 584 if (GetVerboseLevel() > 0) { << 382 if (GetVerboseLevel()>0) { 585 G4cout << "G4PhaseSpaceDecayChannel: << 383 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt "; 586 G4cout << " Cannot calculate dau << 384 G4cout << " can not calculate daughter momentum " <<G4endl; 587 G4cout << " parent:" << *parent_ 385 G4cout << " parent:" << *parent_name; 588 G4cout << " mass:" << parentmass / G << 386 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl; 589 G4cout << " daughter " << index 387 G4cout << " daughter " << index << ":" << *daughters_name[index]; 590 G4cout << " mass:" << daughtermass[i << 388 G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ; 591 G4cout << " mass:" << daughtermoment << 389 G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl; 592 if (useGivenDaughterMass) { << 593 G4cout << "Daughter Mass is given. << 594 } << 595 } 390 } 596 #endif 391 #endif 597 delete[] sm; << 392 delete [] sm; 598 delete[] daughtermass; << 393 delete [] daughtermass; 599 delete[] daughtermomentum; << 394 delete [] daughtermomentum; 600 delete products; << 395 return 0; // Error detection 601 G4Exception("G4PhaseSpaceDecayChannel: << 396 602 "Can not create decay prod << 397 } else { 603 is larger than parent mas << 398 // calculate weight of this events 604 return nullptr; // Error detection << 399 weight *= daughtermomentum[index]/sm[index]; 605 } << 400 if (GetVerboseLevel()>1) { 606 << 401 G4cout << " daughter " << index << ":" << *daughters_name[index]; 607 // calculate weight of this events << 402 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl; 608 weight *= daughtermomentum[index] / sm[i << 403 } 609 #ifdef G4VERBOSE << 610 if (GetVerboseLevel() > 1) { << 611 G4cout << " daughter " << index << << 612 G4cout << " momentum:" << daughtermome << 613 } 404 } 614 #endif << 615 } 405 } 616 << 406 if (GetVerboseLevel()>1) { 617 #ifdef G4VERBOSE << 407 G4cout << " weight: " << weight <<G4endl; 618 if (GetVerboseLevel() > 1) { << 619 G4cout << " weight: " << weight << G4 << 620 } 408 } 621 #endif << 409 622 << 623 // exit if number of Try exceeds 100 410 // exit if number of Try exceeds 100 624 if (++numberOfTry > 100) { << 411 if (numberOfTry++ >100) { 625 #ifdef G4VERBOSE << 412 if (GetVerboseLevel()>0) { 626 if (GetVerboseLevel() > 0) { << 413 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: "; 627 G4cout << "G4PhaseSpaceDecayChannel::M << 414 G4cout << " can not determine Decay Kinematics " << G4endl; 628 G4cout << "Cannot determine Decay Kine << 629 G4cout << "parent : " << *parent_name << 630 G4cout << "daughters : "; << 631 for (index = 0; index < numberOfDaught << 632 G4cout << *daughters_name[index] << << 633 } << 634 G4cout << G4endl; << 635 } 415 } 636 #endif << 416 delete [] sm; 637 G4Exception("G4PhaseSpaceDecayChannel::M << 417 delete [] daughtermass; 638 "Cannot decay: Decay Kinemat << 418 delete [] daughtermomentum; 639 << 419 return 0; // Error detection 640 delete[] sm; << 641 delete[] daughtermass; << 642 delete[] daughtermomentum; << 643 delete products; << 644 return nullptr; // Error detection << 645 } 420 } 646 if (weight < G4UniformRand()) break; << 421 } while ( weight > G4UniformRand()); 647 } << 648 422 649 #ifdef G4VERBOSE 423 #ifdef G4VERBOSE 650 if (GetVerboseLevel() > 1) { << 424 if (GetVerboseLevel()>1) { 651 G4cout << "Start calculation of daughters << 425 G4cout << "Start calulation of daughters momentum vector "<<G4endl; 652 } 426 } 653 #endif << 427 #endif 654 << 655 G4double costheta, sintheta, phi; 428 G4double costheta, sintheta, phi; 656 G4double beta; 429 G4double beta; 657 daughterparticle = new G4DynamicParticle*[nu 430 daughterparticle = new G4DynamicParticle*[numberOfDaughters]; 658 431 659 index = numberOfDaughters - 2; << 432 index = numberOfDaughters -2; 660 costheta = 2. * G4UniformRand() - 1.0; << 433 costheta = 2.*G4UniformRand()-1.0; 661 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 434 sintheta = sqrt((1.0-costheta)*(1.0+costheta)); 662 phi = twopi * G4UniformRand() * rad; << 435 phi = 2.0*M_PI*G4UniformRand()*rad; 663 direction.setZ(costheta); 436 direction.setZ(costheta); 664 direction.setY(sintheta * std::sin(phi)); << 437 direction.setY(sintheta*sin(phi)); 665 direction.setX(sintheta * std::cos(phi)); << 438 direction.setX(sintheta*cos(phi)); 666 daughterparticle[index] = << 439 daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] ); 667 new G4DynamicParticle(G4MT_daughters[index << 440 daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) ); 668 daughterparticle[index + 1] = << 441 669 new G4DynamicParticle(G4MT_daughters[index << 442 for (index = numberOfDaughters -3; index >= 0; index--) { 670 << 443 //calculate momentum direction 671 for (index = numberOfDaughters - 3; index >= << 444 costheta = 2.*G4UniformRand()-1.0; 672 // calculate momentum direction << 445 sintheta = sqrt((1.0-costheta)*(1.0+costheta)); 673 costheta = 2. * G4UniformRand() - 1.0; << 446 phi = 2.0*M_PI*G4UniformRand()*rad; 674 sintheta = std::sqrt((1.0 - costheta) * (1 << 675 phi = twopi * G4UniformRand() * rad; << 676 direction.setZ(costheta); 447 direction.setZ(costheta); 677 direction.setY(sintheta * std::sin(phi)); << 448 direction.setY(sintheta*sin(phi)); 678 direction.setX(sintheta * std::cos(phi)); << 449 direction.setX(sintheta*cos(phi)); 679 450 680 // boost already created particles << 451 // boost already created particles 681 beta = daughtermomentum[index]; 452 beta = daughtermomentum[index]; 682 beta /= << 453 beta /= sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] ); 683 std::sqrt(daughtermomentum[index] * daug << 454 for (G4int index2 = index+1; index2<numberOfDaughters; index2++) { 684 for (index2 = index + 1; index2 < numberOf << 685 G4LorentzVector p4; 455 G4LorentzVector p4; 686 // make G4LorentzVector for secondaries 456 // make G4LorentzVector for secondaries 687 p4 = daughterparticle[index2]->Get4Momen 457 p4 = daughterparticle[index2]->Get4Momentum(); 688 458 689 // boost secondaries to new frame << 459 // boost secondaries to new frame 690 p4.boost(direction.x() * beta, direction << 460 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta); 691 461 692 // change energy/momentum 462 // change energy/momentum 693 daughterparticle[index2]->Set4Momentum(p 463 daughterparticle[index2]->Set4Momentum(p4); 694 } 464 } 695 // create daughter G4DynamicParticle << 465 //create daughter G4DynamicParticle 696 daughterparticle[index] = << 466 daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index])); 697 new G4DynamicParticle(G4MT_daughters[ind << 698 } 467 } 699 468 700 // add daughters to G4Decayproducts << 469 //create G4Decayproducts 701 for (index = 0; index < numberOfDaughters; + << 470 G4DynamicParticle *parentparticle; >> 471 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0); >> 472 parentparticle = new G4DynamicParticle( parent, direction, 0.0); >> 473 products = new G4DecayProducts(*parentparticle); >> 474 delete parentparticle; >> 475 for (index = 0; index<numberOfDaughters; index++) { 702 products->PushProducts(daughterparticle[in 476 products->PushProducts(daughterparticle[index]); 703 } 477 } 704 << 705 #ifdef G4VERBOSE 478 #ifdef G4VERBOSE 706 if (GetVerboseLevel() > 1) { << 479 if (GetVerboseLevel()>1) { 707 G4cout << "G4PhaseSpaceDecayChannel::ManyB << 480 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt "; 708 G4cout << " create decay products in rest << 481 G4cout << " create decay products in rest frame " << G4endl; 709 products->DumpInfo(); 482 products->DumpInfo(); 710 } 483 } 711 #endif 484 #endif 712 << 485 delete [] daughterparticle; 713 delete[] daughterparticle; << 486 delete [] daughtermomentum; 714 delete[] daughtermomentum; << 487 delete [] daughtermass; 715 delete[] daughtermass; << 488 delete [] sm; 716 delete[] sm; << 489 717 << 718 return products; 490 return products; 719 } 491 } 720 492 721 G4bool G4PhaseSpaceDecayChannel::SetDaughterMa << 722 { << 723 for (G4int idx = 0; idx < numberOfDaughters; << 724 givenDaughterMasses[idx] = masses[idx]; << 725 } << 726 useGivenDaughterMass = true; << 727 return useGivenDaughterMass; << 728 } << 729 493 730 G4bool G4PhaseSpaceDecayChannel::SampleDaughte << 494 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2) 731 { 495 { 732 useGivenDaughterMass = false; << 496 // calcurate momentum of daughter particles in two-body decay 733 return useGivenDaughterMass; << 497 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e); >> 498 if (ppp>0) return sqrt(ppp); >> 499 else return -1.; 734 } 500 } 735 501 736 G4bool G4PhaseSpaceDecayChannel::IsOKWithParen << 737 { << 738 if (!useGivenDaughterMass) return G4VDecayCh << 739 502 740 CheckAndFillParent(); << 741 CheckAndFillDaughters(); << 742 503 743 G4double sumOfDaughterMassMin = 0.0; << 744 for (G4int index = 0; index < numberOfDaught << 745 sumOfDaughterMassMin += givenDaughterMasse << 746 } << 747 return (parentMass >= sumOfDaughterMassMin); << 748 } << 749 << 750 G4double G4PhaseSpaceDecayChannel::Pmx(G4doubl << 751 { << 752 // calcurate momentum of daughter particles << 753 G4double ppp = (e + p1 + p2) * (e + p1 - p2) << 754 if (ppp > 0) return std::sqrt(ppp); << 755 return -1.; << 756 } << 757 504