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 // >> 27 // $Id: G4ChargeExchange.cc 91897 2015-08-10 09:55:12Z gcosmo $ >> 28 // >> 29 // 26 // G4 Model: Charge and strangness exchange ba 30 // G4 Model: Charge and strangness exchange based on G4LightMedia model 27 // 28 May 2006 V.Ivanchenko 31 // 28 May 2006 V.Ivanchenko 28 // 32 // 29 // Modified: 33 // Modified: 30 // 07-Jun-06 V.Ivanchenko fix problem of rotat 34 // 07-Jun-06 V.Ivanchenko fix problem of rotation of final state 31 // 25-Jul-06 V.Ivanchenko add 19 MeV low energ 35 // 25-Jul-06 V.Ivanchenko add 19 MeV low energy, below which S-wave is sampled 32 // 12-Jun-12 A.Ribon fix warnings of shadowed 36 // 12-Jun-12 A.Ribon fix warnings of shadowed variables 33 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log 37 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log and G4Pow 34 // 38 // 35 39 36 #include "G4ChargeExchange.hh" 40 #include "G4ChargeExchange.hh" 37 #include "G4ChargeExchangeXS.hh" << 38 #include "G4PhysicalConstants.hh" 41 #include "G4PhysicalConstants.hh" 39 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" 40 #include "G4ParticleTable.hh" 43 #include "G4ParticleTable.hh" 41 #include "G4ParticleDefinition.hh" 44 #include "G4ParticleDefinition.hh" 42 #include "G4IonTable.hh" 45 #include "G4IonTable.hh" 43 #include "Randomize.hh" 46 #include "Randomize.hh" 44 #include "G4NucleiProperties.hh" 47 #include "G4NucleiProperties.hh" 45 #include "G4DecayTable.hh" << 46 #include "G4VDecayChannel.hh" << 47 #include "G4DecayProducts.hh" << 48 #include "G4NistManager.hh" << 49 #include "G4Fragment.hh" << 50 #include "G4ExcitationHandler.hh" << 51 #include "G4ReactionProductVector.hh" << 52 48 53 #include "G4Exp.hh" 49 #include "G4Exp.hh" 54 #include "G4Log.hh" 50 #include "G4Log.hh" 55 #include "G4Pow.hh" 51 #include "G4Pow.hh" 56 52 57 #include "G4HadronicParameters.hh" << 58 #include "G4PhysicsModelCatalog.hh" << 59 53 60 namespace << 54 G4ChargeExchange::G4ChargeExchange() : G4HadronicInteraction("Charge Exchange") 61 { 55 { 62 constexpr G4int maxN = 1000; << 56 SetMinEnergy( 0.0*GeV ); 63 constexpr G4double emin = 2*136.9*CLHEP::MeV << 57 SetMaxEnergy( 100.*TeV ); 64 } << 65 58 66 G4ChargeExchange::G4ChargeExchange(G4ChargeExc << 59 lowEnergyRecoilLimit = 100.*keV; 67 : G4HadronicInteraction("ChargeExchange"), << 60 lowestEnergyLimit = 1.*MeV; 68 fXSection(ptr), fXSWeightFactor(1.0) << 61 69 { << 62 theProton = G4Proton::Proton(); 70 lowEnergyLimit = 1.*CLHEP::MeV; << 63 theNeutron = G4Neutron::Neutron(); 71 secID = G4PhysicsModelCatalog::GetModelID( " << 64 theAProton = G4AntiProton::AntiProton(); 72 nist = G4NistManager::Instance(); << 65 theANeutron = G4AntiNeutron::AntiNeutron(); 73 fHandler = new G4ExcitationHandler(); << 66 thePiPlus = G4PionPlus::PionPlus(); 74 if (nullptr != fXSection) { << 67 thePiMinus = G4PionMinus::PionMinus(); 75 fXSWeightFactor = 1.0/fXSection->GetCrossS << 68 thePiZero = G4PionZero::PionZero(); 76 } << 69 theKPlus = G4KaonPlus::KaonPlus(); >> 70 theKMinus = G4KaonMinus::KaonMinus(); >> 71 theK0S = G4KaonZeroShort::KaonZeroShort(); >> 72 theK0L = G4KaonZeroLong::KaonZeroLong(); >> 73 theL = G4Lambda::Lambda(); >> 74 theAntiL = G4AntiLambda::AntiLambda(); >> 75 theSPlus = G4SigmaPlus::SigmaPlus(); >> 76 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus(); >> 77 theSMinus = G4SigmaMinus::SigmaMinus(); >> 78 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus(); >> 79 theS0 = G4SigmaZero::SigmaZero(); >> 80 theAS0 = G4AntiSigmaZero::AntiSigmaZero(); >> 81 theXiMinus = G4XiMinus::XiMinus(); >> 82 theXi0 = G4XiZero::XiZero(); >> 83 theAXiMinus = G4AntiXiMinus::AntiXiMinus(); >> 84 theAXi0 = G4AntiXiZero::AntiXiZero(); >> 85 theOmega = G4OmegaMinus::OmegaMinus(); >> 86 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus(); >> 87 theD = G4Deuteron::Deuteron(); >> 88 theT = G4Triton::Triton(); >> 89 theA = G4Alpha::Alpha(); >> 90 theHe3 = G4He3::He3(); 77 } 91 } 78 92 79 G4ChargeExchange::~G4ChargeExchange() 93 G4ChargeExchange::~G4ChargeExchange() 80 { << 94 {} 81 delete fHandler; << 82 } << 83 95 84 G4HadFinalState* G4ChargeExchange::ApplyYourse 96 G4HadFinalState* G4ChargeExchange::ApplyYourself( 85 const G4HadProjectile& aTrack, G4Nucleus& 97 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus) 86 { 98 { 87 theParticleChange.Clear(); 99 theParticleChange.Clear(); 88 auto part = aTrack.GetDefinition(); << 100 const G4HadProjectile* aParticle = &aTrack; 89 G4double ekin = aTrack.GetKineticEnergy(); << 101 G4double ekin = aParticle->GetKineticEnergy(); 90 102 91 G4int A = targetNucleus.GetA_asInt(); 103 G4int A = targetNucleus.GetA_asInt(); 92 G4int Z = targetNucleus.GetZ_asInt(); 104 G4int Z = targetNucleus.GetZ_asInt(); 93 105 94 if (ekin <= lowEnergyLimit) { << 106 if(ekin <= lowestEnergyLimit || A < 3) { >> 107 theParticleChange.SetEnergyChange(ekin); >> 108 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 95 return &theParticleChange; 109 return &theParticleChange; 96 } 110 } 97 theParticleChange.SetWeightChange(fXSWeightF << 98 111 99 G4int projPDG = part->GetPDGEncoding(); << 112 G4double plab = aParticle->GetTotalMomentum(); 100 113 101 // for hydrogen targets and positive project << 102 // is not possible on proton, only on deuter << 103 if (1 == Z && (211 == projPDG || 321 == proj << 104 << 105 if (verboseLevel > 1) 114 if (verboseLevel > 1) 106 G4cout << "G4ChargeExchange for " << part- << 115 G4cout << "G4ChargeExchange::DoIt: Incident particle plab=" >> 116 << plab/GeV << " GeV/c " >> 117 << " ekin(MeV) = " << ekin/MeV << " " >> 118 << aParticle->GetDefinition()->GetParticleName() << G4endl; >> 119 >> 120 // Scattered particle referred to axis of incident particle >> 121 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); >> 122 >> 123 G4int N = A - Z; >> 124 G4int projPDG = theParticle->GetPDGEncoding(); >> 125 if (verboseLevel > 1) >> 126 G4cout << "G4ChargeExchange for " << theParticle->GetParticleName() 107 << " PDGcode= " << projPDG << " on nucleu 127 << " PDGcode= " << projPDG << " on nucleus Z= " << Z 108 << " A= " << A << " N= " << A - Z << 128 << " A= " << A << " N= " << N 109 << G4endl; 129 << G4endl; 110 130 111 G4double mass1 = G4NucleiProperties::GetNucl << 131 G4ParticleDefinition * theDef = 0; 112 G4LorentzVector lv0 = aTrack.Get4Momentum(); << 113 G4double etot = mass1 + lv0.e(); << 114 << 115 // select final state << 116 const G4ParticleDefinition* theSecondary = << 117 fXSection->SampleSecondaryType(part, Z, A) << 118 G4int pdg = theSecondary->GetPDGEncoding(); << 119 << 120 // omega(782) and f2(1270) << 121 G4bool isShortLived = (pdg == 223 || pdg == << 122 << 123 // atomic number of the recoil nucleus << 124 if (projPDG == -211) { --Z; } << 125 else if (projPDG == 211) { ++Z; } << 126 else if (projPDG == -321) { --Z; } << 127 else if (projPDG == 321) { ++Z; } << 128 else if (projPDG == 130) { << 129 if (theSecondary->GetPDGCharge() > 0.0) { << 130 else { ++Z; } << 131 } else { << 132 // not ready for other projectile << 133 return &theParticleChange; << 134 } << 135 132 136 // recoil nucleus << 133 G4double mass2 = G4NucleiProperties::GetNuclearMass((G4double)A, (G4double)Z); 137 const G4ParticleDefinition* theRecoil = null << 134 G4LorentzVector lv1 = aParticle->Get4Momentum(); 138 if (Z == 0 && A == 1) { theRecoil = G4Neutro << 135 G4LorentzVector lv0(0.0,0.0,0.0,mass2); 139 else if (Z == 1 && A == 1) { theRecoil = G4P << 140 else if (Z == 1 && A == 2) { theRecoil = G4D << 141 else if (Z == 1 && A == 3) { theRecoil = G4T << 142 else if (Z == 2 && A == 3) { theRecoil = G4H << 143 else if (Z == 2 && A == 4) { theRecoil = G4A << 144 else if (nist->GetIsotopeAbundance(Z, A) > 0 << 145 theRecoil = G4ParticleTable::GetParticleTa << 146 ->GetIonTable()->GetIon(Z, A, 0.0); << 147 } << 148 136 149 // check if there is enough energy for the f << 137 G4LorentzVector lv = lv0 + lv1; 150 // and sample mass of produced state << 138 G4ThreeVector bst = lv.boostVector(); 151 const G4double mass0 = theSecondary->GetPDGM << 139 lv1.boost(-bst); 152 G4double mass3 = (nullptr == theRecoil) ? << 140 lv0.boost(-bst); 153 G4NucleiProperties::GetNuclearMass(A, Z) : << 141 154 G4double mass2 = mass0; << 142 // Sample final particles 155 if (isShortLived && << 143 G4bool theHyperon = false; 156 !SampleMass(mass2, theSecondary->GetPDGW << 144 G4ParticleDefinition* theRecoil = 0; 157 return &theParticleChange; << 145 G4ParticleDefinition* theSecondary = 0; >> 146 >> 147 if(theParticle == theProton) { >> 148 theSecondary = theNeutron; >> 149 Z++; >> 150 } else if(theParticle == theNeutron) { >> 151 theSecondary = theProton; >> 152 Z--; >> 153 } else if(theParticle == thePiPlus) { >> 154 theSecondary = thePiZero; >> 155 Z++; >> 156 } else if(theParticle == thePiMinus) { >> 157 theSecondary = thePiZero; >> 158 Z--; >> 159 } else if(theParticle == theKPlus) { >> 160 if(G4UniformRand()<0.5) theSecondary = theK0S; >> 161 else theSecondary = theK0L; >> 162 Z++; >> 163 } else if(theParticle == theKMinus) { >> 164 if(G4UniformRand()<0.5) theSecondary = theK0S; >> 165 else theSecondary = theK0L; >> 166 Z--; >> 167 } else if(theParticle == theK0S || theParticle == theK0L) { >> 168 if(G4UniformRand()*A < G4double(Z)) { >> 169 theSecondary = theKPlus; >> 170 Z--; >> 171 } else { >> 172 theSecondary = theKMinus; >> 173 Z++; >> 174 } >> 175 } else if(theParticle == theANeutron) { >> 176 theSecondary = theAProton; >> 177 Z++; >> 178 } else if(theParticle == theAProton) { >> 179 theSecondary = theANeutron; >> 180 Z--; >> 181 } else if(theParticle == theL) { >> 182 G4double x = G4UniformRand(); >> 183 if(G4UniformRand()*A < G4double(Z)) { >> 184 if(x < 0.2) { >> 185 theSecondary = theS0; >> 186 } else if (x < 0.4) { >> 187 theSecondary = theSPlus; >> 188 Z--; >> 189 } else if (x < 0.6) { >> 190 theSecondary = theProton; >> 191 theRecoil = theL; >> 192 theHyperon = true; >> 193 A--; >> 194 } else if (x < 0.8) { >> 195 theSecondary = theProton; >> 196 theRecoil = theS0; >> 197 theHyperon = true; >> 198 A--; >> 199 } else { >> 200 theSecondary = theNeutron; >> 201 theRecoil = theSPlus; >> 202 theHyperon = true; >> 203 A--; >> 204 } >> 205 } else { >> 206 if(x < 0.2) { >> 207 theSecondary = theS0; >> 208 } else if (x < 0.4) { >> 209 theSecondary = theSMinus; >> 210 Z++; >> 211 } else if (x < 0.6) { >> 212 theSecondary = theNeutron; >> 213 theRecoil = theL; >> 214 A--; >> 215 theHyperon = true; >> 216 } else if (x < 0.8) { >> 217 theSecondary = theNeutron; >> 218 theRecoil = theS0; >> 219 theHyperon = true; >> 220 A--; >> 221 } else { >> 222 theSecondary = theProton; >> 223 theRecoil = theSMinus; >> 224 theHyperon = true; >> 225 A--; >> 226 } >> 227 } 158 } 228 } 159 229 160 // not possible kinematically << 230 if (Z == 1 && A == 2) theDef = theD; 161 if (etot <= mass2 + mass3) { << 231 else if (Z == 1 && A == 3) theDef = theT; >> 232 else if (Z == 2 && A == 3) theDef = theHe3; >> 233 else if (Z == 2 && A == 4) theDef = theA; >> 234 else { >> 235 theDef = >> 236 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0); >> 237 } >> 238 if(!theSecondary) { return &theParticleChange; } >> 239 >> 240 G4double m11 = theSecondary->GetPDGMass(); >> 241 G4double m21 = theDef->GetPDGMass(); >> 242 if(theRecoil) { m21 += theRecoil->GetPDGMass(); } >> 243 else { theRecoil = theDef; } >> 244 >> 245 G4double etot = lv0.e() + lv1.e(); >> 246 >> 247 // kinematiacally impossible >> 248 if(etot < m11 + m21) { >> 249 theParticleChange.SetEnergyChange(ekin); >> 250 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 162 return &theParticleChange; 251 return &theParticleChange; 163 } 252 } 164 253 165 // sample kinematics << 254 G4ThreeVector p1 = lv1.vect(); 166 G4LorentzVector lv1(0.0, 0.0, 0.0, mass1); << 255 G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot)); 167 G4LorentzVector lv = lv0 + lv1; << 256 // G4double e2 = etot - e1; 168 G4ThreeVector bst = lv.boostVector(); << 257 G4double ptot = std::sqrt(e1*e1 - m11*m11); 169 G4double ss = lv.mag2(); << 258 >> 259 G4double tmax = 4.0*ptot*ptot; >> 260 G4double g2 = GeV*GeV; >> 261 >> 262 G4double t = g2*SampleT(tmax/g2, A); >> 263 >> 264 if(verboseLevel>1) >> 265 G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax >> 266 << " ptot= " << ptot << G4endl; 170 267 171 // tmax = 4*momCMS^2 << 268 // Sampling in CM system 172 G4double e2 = ss + mass2*mass2 - mass3*mass3 << 269 G4double phi = G4UniformRand()*twopi; 173 G4double tmax = e2*e2/ss - 4*mass2*mass2; << 174 << 175 G4double t = SampleT(theSecondary, A, tmax); << 176 << 177 G4double phi = G4UniformRand()*CLHEP::twopi << 178 G4double cost = 1. - 2.0*t/tmax; 270 G4double cost = 1. - 2.0*t/tmax; >> 271 if(std::abs(cost) > 1.0) cost = 1.0; >> 272 G4double sint = std::sqrt((1.0-cost)*(1.0+cost)); 179 273 180 if (cost > 1.0) { cost = 1.0; } << 274 //if (verboseLevel > 1) 181 else if(cost < -1.0) { cost = -1.0; } << 275 // G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl; 182 276 183 G4double sint = std::sqrt((1.0-cost)*(1.0+co << 277 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); >> 278 v1 *= ptot; >> 279 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1); >> 280 G4LorentzVector nlv0 = lv0 + lv1 - nlv1; 184 281 185 if (verboseLevel>1) { << 282 nlv0.boost(bst); 186 G4cout << " t= " << t << " tmax(GeV^2)= " << 283 nlv1.boost(bst); 187 << " cos(t)=" << cost << " sin(t)=" << si << 188 } << 189 G4double momentumCMS = 0.5*std::sqrt(tmax); << 190 G4LorentzVector lv2(momentumCMS*sint*std::co << 191 momentumCMS*sint*std::sin(phi), << 192 momentumCMS*cost, << 193 std::sqrt(momentumCMS*momentumCMS + << 194 << 195 // kinematics in the final state, may be a w << 196 lv2.boost(bst); << 197 if (lv2.e() < mass2) { << 198 lv2.setE(mass2); << 199 } << 200 lv -= lv2; << 201 if (lv.e() < mass3) { << 202 lv.setE(mass3); << 203 } << 204 284 205 // prepare secondary particles << 206 theParticleChange.SetStatusChange(stopAndKil 285 theParticleChange.SetStatusChange(stopAndKill); 207 theParticleChange.SetEnergyChange(0.0); 286 theParticleChange.SetEnergyChange(0.0); >> 287 G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1); >> 288 theParticleChange.AddSecondary(aSec); 208 289 209 if (!isShortLived) { << 290 G4double erec = nlv0.e() - m21; 210 auto aSec = new G4DynamicParticle(theSecon << 211 theParticleChange.AddSecondary(aSec, secID << 212 } else { << 213 auto channel = theSecondary->GetDecayTable << 214 auto products = channel->DecayIt(mass2); << 215 G4ThreeVector bst1 = lv2.boostVector(); << 216 G4int N = products->entries(); << 217 for (G4int i=0; i<N; ++i) { << 218 auto p = (*products)[i]; << 219 auto lvp = p->Get4Momentum(); << 220 lvp.boost(bst1); << 221 p->Set4Momentum(lvp); << 222 theParticleChange.AddSecondary(p, secID) << 223 } << 224 delete products; << 225 } << 226 291 227 // recoil is a stable isotope << 292 //G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl; 228 if (nullptr != theRecoil) { << 293 229 auto aRec = new G4DynamicParticle(theRecoi << 294 if(theHyperon) { 230 theParticleChange.AddSecondary(aRec, secID << 295 theParticleChange.SetLocalEnergyDeposit(erec); >> 296 aSec = new G4DynamicParticle(); >> 297 aSec->SetDefinition(theRecoil); >> 298 aSec->SetKineticEnergy(0.0); >> 299 } else if(erec > lowEnergyRecoilLimit) { >> 300 aSec = new G4DynamicParticle(theRecoil, nlv0); >> 301 theParticleChange.AddSecondary(aSec); 231 } else { 302 } else { 232 // recoil is an unstable fragment << 303 if(erec < 0.0) erec = 0.0; 233 G4Fragment frag(A, Z, lv); << 304 theParticleChange.SetLocalEnergyDeposit(erec); 234 auto products = fHandler->BreakItUp(frag); << 235 for (auto & prod : *products) { << 236 auto dp = new G4DynamicParticle(prod->Ge << 237 theParticleChange.AddSecondary(dp, secID << 238 delete prod; << 239 } << 240 delete products; << 241 } 305 } 242 return &theParticleChange; 306 return &theParticleChange; 243 } 307 } 244 308 245 G4double G4ChargeExchange::SampleT(const G4Par << 309 G4double G4ChargeExchange::SampleT(G4double tmax, G4double A) 246 const G4int << 247 { 310 { 248 G4double aa, bb, cc, dd; 311 G4double aa, bb, cc, dd; 249 G4Pow* g4pow = G4Pow::GetInstance(); << 250 if (A <= 62.) { 312 if (A <= 62.) { 251 aa = g4pow->powZ(A, 1.63); << 313 aa = G4Pow::GetInstance()->powA(A, 1.63); 252 bb = 14.5*g4pow->powZ(A, 0.66); << 314 bb = 14.5*G4Pow::GetInstance()->powA(A, 0.66); 253 cc = 1.4*g4pow->powZ(A, 0.33); << 315 cc = 1.4*G4Pow::GetInstance()->powA(A, 0.33); 254 dd = 10.; 316 dd = 10.; 255 } else { 317 } else { 256 aa = g4pow->powZ(A, 1.33); << 318 aa = G4Pow::GetInstance()->powA(A, 1.33); 257 bb = 60.*g4pow->powZ(A, 0.33); << 319 bb = 60.*G4Pow::GetInstance()->powA(A, 0.33); 258 cc = 0.4*g4pow->powZ(A, 0.40); << 320 cc = 0.4*G4Pow::GetInstance()->powA(A, 0.40); 259 dd = 10.; 321 dd = 10.; 260 } 322 } 261 G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb; 323 G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb; 262 G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd; 324 G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd; 263 325 264 G4double t; 326 G4double t; 265 G4double y = bb; 327 G4double y = bb; 266 if(G4UniformRand()*(x1 + x2) < x2) y = dd; 328 if(G4UniformRand()*(x1 + x2) < x2) y = dd; 267 329 268 for (G4int i=0; i<maxN; ++i) { << 330 const G4int maxNumberOfLoops = 10000; >> 331 G4int loopCounter = 0; >> 332 do { 269 t = -G4Log(G4UniformRand())/y; 333 t = -G4Log(G4UniformRand())/y; 270 if (t <= tmax) { return t; } << 334 } while ( (t > tmax) && >> 335 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */ >> 336 if ( loopCounter >= maxNumberOfLoops ) { >> 337 t = 0.0; 271 } 338 } 272 return 0.0; << 273 } << 274 339 275 G4bool G4ChargeExchange::SampleMass(G4double& << 340 return t; 276 { << 277 // +- 4 width but above 2 pion mass << 278 const G4double e1 = std::max(M - 4*G, emin); << 279 const G4double e2 = std::min(M + 4*G, elim) << 280 if (e2 <= 0.0) { return false; } << 281 const G4double M2 = M*M; << 282 const G4double MG2 = M2*G*G; << 283 << 284 // sampling Breit-Wigner function << 285 for (G4int i=0; i<maxN; ++i) { << 286 G4double e = e1 + e2*G4UniformRand(); << 287 G4double x = e*e - M2; << 288 G4double y = MG2/(x*x + MG2); << 289 if (y >= G4UniformRand()) { << 290 M = e; << 291 return true; << 292 } << 293 } << 294 return false; << 295 } 341 } >> 342 296 343