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 // G4 Model: Charge and strangness exchange ba 26 // G4 Model: Charge and strangness exchange based on G4LightMedia model 27 // 28 May 2006 V.Ivanchenko 27 // 28 May 2006 V.Ivanchenko 28 // 28 // 29 // Modified: 29 // Modified: 30 // 07-Jun-06 V.Ivanchenko fix problem of rotat 30 // 07-Jun-06 V.Ivanchenko fix problem of rotation of final state 31 // 25-Jul-06 V.Ivanchenko add 19 MeV low energ 31 // 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 32 // 12-Jun-12 A.Ribon fix warnings of shadowed variables 33 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log 33 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log and G4Pow 34 // 34 // 35 35 36 #include "G4ChargeExchange.hh" 36 #include "G4ChargeExchange.hh" 37 #include "G4ChargeExchangeXS.hh" 37 #include "G4ChargeExchangeXS.hh" 38 #include "G4PhysicalConstants.hh" 38 #include "G4PhysicalConstants.hh" 39 #include "G4SystemOfUnits.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 "Randomize.hh" 43 #include "Randomize.hh" 44 #include "G4NucleiProperties.hh" 44 #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 45 53 #include "G4Exp.hh" 46 #include "G4Exp.hh" 54 #include "G4Log.hh" 47 #include "G4Log.hh" 55 #include "G4Pow.hh" 48 #include "G4Pow.hh" 56 49 57 #include "G4HadronicParameters.hh" 50 #include "G4HadronicParameters.hh" 58 #include "G4PhysicsModelCatalog.hh" 51 #include "G4PhysicsModelCatalog.hh" 59 52 60 namespace << 61 { << 62 constexpr G4int maxN = 1000; << 63 constexpr G4double emin = 2*136.9*CLHEP::MeV << 64 } << 65 53 66 G4ChargeExchange::G4ChargeExchange(G4ChargeExc 54 G4ChargeExchange::G4ChargeExchange(G4ChargeExchangeXS* ptr) 67 : G4HadronicInteraction("ChargeExchange"), 55 : G4HadronicInteraction("ChargeExchange"), 68 fXSection(ptr), fXSWeightFactor(1.0) << 56 fXSection(ptr) 69 { 57 { 70 lowEnergyLimit = 1.*CLHEP::MeV; 58 lowEnergyLimit = 1.*CLHEP::MeV; 71 secID = G4PhysicsModelCatalog::GetModelID( " 59 secID = G4PhysicsModelCatalog::GetModelID( "model_ChargeExchange" ); 72 nist = G4NistManager::Instance(); << 73 fHandler = new G4ExcitationHandler(); << 74 if (nullptr != fXSection) { << 75 fXSWeightFactor = 1.0/fXSection->GetCrossS << 76 } << 77 } << 78 << 79 G4ChargeExchange::~G4ChargeExchange() << 80 { << 81 delete fHandler; << 82 } 60 } 83 61 84 G4HadFinalState* G4ChargeExchange::ApplyYourse 62 G4HadFinalState* G4ChargeExchange::ApplyYourself( 85 const G4HadProjectile& aTrack, G4Nucleus& 63 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus) 86 { 64 { 87 theParticleChange.Clear(); 65 theParticleChange.Clear(); 88 auto part = aTrack.GetDefinition(); 66 auto part = aTrack.GetDefinition(); >> 67 G4int pdg = part->GetPDGEncoding(); 89 G4double ekin = aTrack.GetKineticEnergy(); 68 G4double ekin = aTrack.GetKineticEnergy(); 90 69 91 G4int A = targetNucleus.GetA_asInt(); 70 G4int A = targetNucleus.GetA_asInt(); 92 G4int Z = targetNucleus.GetZ_asInt(); 71 G4int Z = targetNucleus.GetZ_asInt(); 93 72 94 if (ekin <= lowEnergyLimit) { 73 if (ekin <= lowEnergyLimit) { 95 return &theParticleChange; 74 return &theParticleChange; 96 } 75 } 97 theParticleChange.SetWeightChange(fXSWeightF << 98 76 99 G4int projPDG = part->GetPDGEncoding(); 77 G4int projPDG = part->GetPDGEncoding(); 100 << 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) 78 if (verboseLevel > 1) 106 G4cout << "G4ChargeExchange for " << part- 79 G4cout << "G4ChargeExchange for " << part->GetParticleName() 107 << " PDGcode= " << projPDG << " on nucleu 80 << " PDGcode= " << projPDG << " on nucleus Z= " << Z 108 << " A= " << A << " N= " << A - Z 81 << " A= " << A << " N= " << A - Z 109 << G4endl; 82 << G4endl; 110 83 111 G4double mass1 = G4NucleiProperties::GetNucl 84 G4double mass1 = G4NucleiProperties::GetNuclearMass(A, Z); 112 G4LorentzVector lv0 = aTrack.Get4Momentum(); 85 G4LorentzVector lv0 = aTrack.Get4Momentum(); 113 G4double etot = mass1 + lv0.e(); 86 G4double etot = mass1 + lv0.e(); 114 87 115 // select final state 88 // select final state >> 89 const G4ParticleDefinition* theRecoil = nullptr; 116 const G4ParticleDefinition* theSecondary = 90 const G4ParticleDefinition* theSecondary = 117 fXSection->SampleSecondaryType(part, Z, A) 91 fXSection->SampleSecondaryType(part, Z, A); 118 G4int pdg = theSecondary->GetPDGEncoding(); << 119 << 120 // omega(782) and f2(1270) << 121 G4bool isShortLived = (pdg == 223 || pdg == << 122 92 123 // atomic number of the recoil nucleus 93 // atomic number of the recoil nucleus 124 if (projPDG == -211) { --Z; } << 94 if (pdg == -211) { --Z; } 125 else if (projPDG == 211) { ++Z; } << 95 else if (pdg == 211) { ++Z; } 126 else if (projPDG == -321) { --Z; } << 96 else if (pdg == -321) { --Z; } 127 else if (projPDG == 321) { ++Z; } << 97 else if (pdg == 321) { ++Z; } 128 else if (projPDG == 130) { << 98 else if (pdg == 130) { 129 if (theSecondary->GetPDGCharge() > 0.0) { 99 if (theSecondary->GetPDGCharge() > 0.0) { --Z; } 130 else { ++Z; } 100 else { ++Z; } 131 } else { 101 } else { 132 // not ready for other projectile 102 // not ready for other projectile 133 return &theParticleChange; 103 return &theParticleChange; 134 } 104 } 135 105 136 // recoil nucleus << 106 if (Z == 0 && A == 1) theRecoil = G4Neutron::Neutron(); 137 const G4ParticleDefinition* theRecoil = null << 107 else if (Z == 1 && A == 1) theRecoil = G4Proton::Proton(); 138 if (Z == 0 && A == 1) { theRecoil = G4Neutro << 108 else if (Z == 1 && A == 2) theRecoil = G4Deuteron::Deuteron(); 139 else if (Z == 1 && A == 1) { theRecoil = G4P << 109 else if (Z == 1 && A == 3) theRecoil = G4Triton::Triton(); 140 else if (Z == 1 && A == 2) { theRecoil = G4D << 110 else if (Z == 2 && A == 3) theRecoil = G4He3::He3(); 141 else if (Z == 1 && A == 3) { theRecoil = G4T << 111 else if (Z == 2 && A == 4) theRecoil = G4Alpha::Alpha(); 142 else if (Z == 2 && A == 3) { theRecoil = G4H << 112 else { 143 else if (Z == 2 && A == 4) { theRecoil = G4A << 144 else if (nist->GetIsotopeAbundance(Z, A) > 0 << 145 theRecoil = G4ParticleTable::GetParticleTa 113 theRecoil = G4ParticleTable::GetParticleTable() 146 ->GetIonTable()->GetIon(Z, A, 0.0); 114 ->GetIonTable()->GetIon(Z, A, 0.0); 147 } 115 } >> 116 if (nullptr == theRecoil) { return &theParticleChange; } 148 117 149 // check if there is enough energy for the f << 118 G4double mass2 = theSecondary->GetPDGMass(); 150 // and sample mass of produced state << 119 G4double mass3 = theRecoil->GetPDGMass(); 151 const G4double mass0 = theSecondary->GetPDGM << 152 G4double mass3 = (nullptr == theRecoil) ? << 153 G4NucleiProperties::GetNuclearMass(A, Z) : << 154 G4double mass2 = mass0; << 155 if (isShortLived && << 156 !SampleMass(mass2, theSecondary->GetPDGW << 157 return &theParticleChange; << 158 } << 159 120 160 // not possible kinematically << 161 if (etot <= mass2 + mass3) { 121 if (etot <= mass2 + mass3) { >> 122 // not possible kinematically 162 return &theParticleChange; 123 return &theParticleChange; 163 } 124 } 164 125 165 // sample kinematics 126 // sample kinematics 166 G4LorentzVector lv1(0.0, 0.0, 0.0, mass1); 127 G4LorentzVector lv1(0.0, 0.0, 0.0, mass1); 167 G4LorentzVector lv = lv0 + lv1; 128 G4LorentzVector lv = lv0 + lv1; 168 G4ThreeVector bst = lv.boostVector(); 129 G4ThreeVector bst = lv.boostVector(); 169 G4double ss = lv.mag2(); 130 G4double ss = lv.mag2(); 170 131 171 // tmax = 4*momCMS^2 132 // tmax = 4*momCMS^2 172 G4double e2 = ss + mass2*mass2 - mass3*mass3 133 G4double e2 = ss + mass2*mass2 - mass3*mass3; 173 G4double tmax = e2*e2/ss - 4*mass2*mass2; 134 G4double tmax = e2*e2/ss - 4*mass2*mass2; 174 135 175 G4double t = SampleT(theSecondary, A, tmax); 136 G4double t = SampleT(theSecondary, A, tmax); 176 137 177 G4double phi = G4UniformRand()*CLHEP::twopi 138 G4double phi = G4UniformRand()*CLHEP::twopi; 178 G4double cost = 1. - 2.0*t/tmax; 139 G4double cost = 1. - 2.0*t/tmax; 179 140 180 if (cost > 1.0) { cost = 1.0; } 141 if (cost > 1.0) { cost = 1.0; } 181 else if(cost < -1.0) { cost = -1.0; } 142 else if(cost < -1.0) { cost = -1.0; } 182 143 183 G4double sint = std::sqrt((1.0-cost)*(1.0+co 144 G4double sint = std::sqrt((1.0-cost)*(1.0+cost)); 184 145 185 if (verboseLevel>1) { 146 if (verboseLevel>1) { 186 G4cout << " t= " << t << " tmax(GeV^2)= " 147 G4cout << " t= " << t << " tmax(GeV^2)= " << tmax/(GeV*GeV) 187 << " cos(t)=" << cost << " sin(t)=" << si 148 << " cos(t)=" << cost << " sin(t)=" << sint << G4endl; 188 } 149 } 189 G4double momentumCMS = 0.5*std::sqrt(tmax); 150 G4double momentumCMS = 0.5*std::sqrt(tmax); 190 G4LorentzVector lv2(momentumCMS*sint*std::co 151 G4LorentzVector lv2(momentumCMS*sint*std::cos(phi), 191 momentumCMS*sint*std::sin(phi), 152 momentumCMS*sint*std::sin(phi), 192 momentumCMS*cost, 153 momentumCMS*cost, 193 std::sqrt(momentumCMS*momentumCMS + 154 std::sqrt(momentumCMS*momentumCMS + mass2*mass2)); 194 155 195 // kinematics in the final state, may be a w << 156 // kinematics in the final stae, may be a warning should be added if 196 lv2.boost(bst); 157 lv2.boost(bst); 197 if (lv2.e() < mass2) { 158 if (lv2.e() < mass2) { 198 lv2.setE(mass2); 159 lv2.setE(mass2); 199 } 160 } 200 lv -= lv2; 161 lv -= lv2; 201 if (lv.e() < mass3) { 162 if (lv.e() < mass3) { 202 lv.setE(mass3); 163 lv.setE(mass3); 203 } 164 } 204 165 205 // prepare secondary particles 166 // prepare secondary particles 206 theParticleChange.SetStatusChange(stopAndKil 167 theParticleChange.SetStatusChange(stopAndKill); 207 theParticleChange.SetEnergyChange(0.0); 168 theParticleChange.SetEnergyChange(0.0); 208 169 209 if (!isShortLived) { << 170 G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, lv2); 210 auto aSec = new G4DynamicParticle(theSecon << 171 theParticleChange.AddSecondary(aSec, secID); 211 theParticleChange.AddSecondary(aSec, secID << 172 212 } else { << 173 G4DynamicParticle * aRec = new G4DynamicParticle(theRecoil, lv); 213 auto channel = theSecondary->GetDecayTable << 174 theParticleChange.AddSecondary(aRec, secID); 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 << 227 // recoil is a stable isotope << 228 if (nullptr != theRecoil) { << 229 auto aRec = new G4DynamicParticle(theRecoi << 230 theParticleChange.AddSecondary(aRec, secID << 231 } else { << 232 // recoil is an unstable fragment << 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 } << 242 return &theParticleChange; 175 return &theParticleChange; 243 } 176 } 244 177 245 G4double G4ChargeExchange::SampleT(const G4Par 178 G4double G4ChargeExchange::SampleT(const G4ParticleDefinition*, 246 const G4int << 179 G4int A, G4double tmax) const 247 { 180 { 248 G4double aa, bb, cc, dd; 181 G4double aa, bb, cc, dd; 249 G4Pow* g4pow = G4Pow::GetInstance(); 182 G4Pow* g4pow = G4Pow::GetInstance(); 250 if (A <= 62.) { 183 if (A <= 62.) { 251 aa = g4pow->powZ(A, 1.63); 184 aa = g4pow->powZ(A, 1.63); 252 bb = 14.5*g4pow->powZ(A, 0.66); 185 bb = 14.5*g4pow->powZ(A, 0.66); 253 cc = 1.4*g4pow->powZ(A, 0.33); 186 cc = 1.4*g4pow->powZ(A, 0.33); 254 dd = 10.; 187 dd = 10.; 255 } else { 188 } else { 256 aa = g4pow->powZ(A, 1.33); 189 aa = g4pow->powZ(A, 1.33); 257 bb = 60.*g4pow->powZ(A, 0.33); 190 bb = 60.*g4pow->powZ(A, 0.33); 258 cc = 0.4*g4pow->powZ(A, 0.40); 191 cc = 0.4*g4pow->powZ(A, 0.40); 259 dd = 10.; 192 dd = 10.; 260 } 193 } 261 G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb; 194 G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb; 262 G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd; 195 G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd; 263 196 264 G4double t; 197 G4double t; 265 G4double y = bb; 198 G4double y = bb; 266 if(G4UniformRand()*(x1 + x2) < x2) y = dd; 199 if(G4UniformRand()*(x1 + x2) < x2) y = dd; 267 200 268 for (G4int i=0; i<maxN; ++i) { << 201 const G4int maxN = 10000; >> 202 G4int count = 0; >> 203 do { 269 t = -G4Log(G4UniformRand())/y; 204 t = -G4Log(G4UniformRand())/y; 270 if (t <= tmax) { return t; } << 205 } while ( (t > tmax) && ++count < maxN ); 271 } << 206 /* Loop checking, 10.08.2015, A.Ribon */ 272 return 0.0; << 207 if ( count >= maxN ) { 273 } << 208 t = 0.0; 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 } 209 } 294 return false; << 210 return t; 295 } 211 } 296 212