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