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