Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 35 36 #include "globals.hh" 37 38 /* 39 * G4INCLReflectionAvatar.cc 40 * 41 * \date Jun 8, 2009 42 * \author Pekka Kaitaniemi 43 */ 44 45 #include "G4INCLSurfaceAvatar.hh" 46 #include "G4INCLRandom.hh" 47 #include "G4INCLReflectionChannel.hh" 48 #include "G4INCLTransmissionChannel.hh" 49 #include "G4INCLClustering.hh" 50 #include <sstream> 51 #include <string> 52 53 namespace G4INCL { 54 55 SurfaceAvatar::SurfaceAvatar(G4INCL::Particl 56 :IAvatar(time), theParticle(aParticle), th 57 particlePIn(0.), 58 particlePOut(0.), 59 particleTOut(0.), 60 TMinusV(0.), 61 TMinusV2(0.), 62 particleMass(0.), 63 sinIncidentAngle(0.), 64 cosIncidentAngle(0.), 65 sinRefractionAngle(0.), 66 cosRefractionAngle(0.), 67 refractionIndexRatio(0.), 68 internalReflection(false) 69 { 70 setType(SurfaceAvatarType); 71 } 72 73 SurfaceAvatar::~SurfaceAvatar() 74 { 75 } 76 77 G4INCL::IChannel* SurfaceAvatar::getChannel( 78 if(theParticle->isTargetSpectator()) { 79 INCL_DEBUG("Particle " << theParticle->g 80 return new ReflectionChannel(theNucleus, 81 } 82 83 // We forbid transmission of resonances be 84 // delta particle below Tf can lead to neg 85 // CDPP assumes that particles stay in the 86 if(theParticle->isResonance()) { 87 const G4double theFermiEnergy = theNucle 88 if(theParticle->getKineticEnergy()<theFe 89 INCL_DEBUG("Particle " << theParticle- 90 << " Tf=" << theFermiEnergy << ", 91 return new ReflectionChannel(theNucleu 92 } 93 } 94 95 // Don't try to make a cluster if the lead 96 const G4double transmissionProbability = g 97 const G4double TOut = TMinusV; 98 const G4double kOut = particlePOut; 99 const G4double cosR = cosRefractionAngle; 100 101 INCL_DEBUG("Transmission probability for p 102 /* Don't attempt to construct clusters whe 103 * trying to escape during a nucleus-nucle 104 * this is that projectile spectators will 105 * projectile remnant, and trying to clust 106 * G4double counting. Moreover, applying t 107 * projectile spectators makes the code *r 108 * large. 109 */ 110 if(theParticle->isNucleonorLambda() 111 && (!theParticle->isProjectileSpectato 112 && transmissionProbability>1.E-4) { 113 Cluster *candidateCluster = 0; 114 115 candidateCluster = Clustering::getCluste 116 if(candidateCluster != 0 && 117 Clustering::clusterCanEscape(theNucl 118 119 INCL_DEBUG("Cluster algorithm succeede 120 121 // Check if the cluster can penetrate 122 const G4double clusterTransmissionProb 123 const G4double x = Random::shoot(); 124 125 INCL_DEBUG("Transmission probability f 126 127 if (x <= clusterTransmissionProbabilit 128 theNucleus->getStore()->getBook().in 129 INCL_DEBUG("Cluster " << candidateCl 130 return new TransmissionChannel(theNu 131 } else { 132 INCL_DEBUG("Cluster " << candidateCl 133 delete candidateCluster; 134 } 135 } else { 136 delete candidateCluster; 137 } 138 } 139 140 // If we haven't transmitted a cluster (ma 141 // disabled or maybe we just can't produce 142 143 // Always transmit projectile spectators i 144 // transmission is energetically allowed 145 if(theParticle->isProjectileSpectator() && 146 INCL_DEBUG("Particle " << theParticle->g 147 return new TransmissionChannel(theNucleu 148 } 149 150 // Transmit or reflect depending on the tr 151 const G4double x = Random::shoot(); 152 153 if(x <= transmissionProbability) { // Tran 154 INCL_DEBUG("Particle " << theParticle->g 155 if(theParticle->isKaon()) theNucleus->se 156 if(theNucleus->getStore()->getConfig()-> 157 return new TransmissionChannel(theNucl 158 } else { 159 return new TransmissionChannel(theNucl 160 } 161 } else { // Reflection 162 INCL_DEBUG("Particle " << theParticle->g 163 return new ReflectionChannel(theNucleus, 164 } 165 } 166 167 void SurfaceAvatar::fillFinalState(FinalStat 168 getChannel()->fillFinalState(fs); 169 } 170 171 void SurfaceAvatar::preInteraction() {} 172 173 void SurfaceAvatar::postInteraction(FinalSta 174 ParticleList const &outgoing = fs->getOutg 175 if(!outgoing.empty()) { // Transmission 176 // assert(outgoing.size()==1); 177 Particle *out = outgoing.front(); 178 out->rpCorrelate(); 179 if(out->isCluster()) { 180 Cluster *clusterOut = dynamic_cast<Clu 181 ParticleList const &components = clust 182 for(ParticleIter i=components.begin(), 183 if(!(*i)->isTargetSpectator()) 184 theNucleus->getStore()->getBook(). 185 } 186 out->setBiasCollisionVector(components 187 } else if(!theParticle->isTargetSpectato 188 // assert(out==theParticle); 189 theNucleus->getStore()->getBook().decr 190 } 191 } 192 } 193 194 std::string SurfaceAvatar::dump() const { 195 std::stringstream ss; 196 ss << "(avatar " << theTime << " 'reflecti 197 << "(list " << '\n' 198 << theParticle->dump() 199 << "))" << '\n'; 200 return ss.str(); 201 } 202 203 G4double SurfaceAvatar::getTransmissionProba 204 205 particleMass = particle->getMass(); 206 const G4double V = particle->getPotentialE 207 208 // Correction to the particle kinetic ener 209 const G4int theA = theNucleus->getA(); 210 const G4int theZ = theNucleus->getZ(); 211 const G4int theS = theNucleus->getS(); 212 const G4double correction = particle->getE 213 particleTOut = particle->getKineticEnergy( 214 215 if (particleTOut <= V) // No transmission 216 return 0.0; 217 218 TMinusV = particleTOut-V; 219 TMinusV2 = TMinusV*TMinusV; 220 221 // Momenta in and out 222 const G4double particlePIn2 = particle->g 223 const G4double particlePOut2 = 2.*particle 224 particlePIn = std::sqrt(particlePIn2); 225 particlePOut = std::sqrt(particlePOut2); 226 227 if (0. > V) // Automatic transmission for 228 return 1.0; 229 230 // Compute the transmission probability 231 G4double theTransmissionProbability; 232 if(theNucleus->getStore()->getConfig()->ge 233 // Use the formula with refraction 234 initializeRefractionVariables(particle); 235 236 if(internalReflection) 237 return 0.; // total internal reflectio 238 239 // Intermediate variables for calculatio 240 const G4double x = refractionIndexRatio* 241 const G4double y = (x - cosRefractionAng 242 243 theTransmissionProbability = 1. - y*y; 244 } else { 245 // Use the formula without refraction 246 // Intermediate variable for calculation 247 const G4double y = particlePIn+particleP 248 249 // The transmission probability for a po 250 theTransmissionProbability = 4.*particle 251 } 252 253 // For neutral and negative particles, no 254 // Also, no Coulomb if the particle takes 255 const G4int particleZ = particle->getZ(); 256 if (particleZ <= 0 || particleZ >= theZ) 257 return theTransmissionProbability; 258 259 // Nominal Coulomb barrier 260 const G4double theTransmissionBarrier = th 261 if (TMinusV >= theTransmissionBarrier) // 262 return theTransmissionProbability; 263 264 // Coulomb-penetration factor 265 const G4double px = std::sqrt(TMinusV/theT 266 const G4double logCoulombTransmission = 267 particleZ*(theZ-particleZ)/137.03*std::s 268 *(Math::arcCos(px)-px*std::sqrt(1.-px*px 269 INCL_DEBUG("Coulomb barrier, logCoulombTra 270 if (logCoulombTransmission > 35.) // Trans 271 return 0.; 272 theTransmissionProbability *= std::exp(-2. 273 274 return theTransmissionProbability; 275 } 276 277 void SurfaceAvatar::initializeRefractionVari 278 cosIncidentAngle = particle->getCosRPAngle 279 if(cosIncidentAngle>1.) 280 cosIncidentAngle=1.; 281 sinIncidentAngle = std::sqrt(1. - cosIncid 282 refractionIndexRatio = particlePIn/particl 283 const G4double sinCandidate = refractionIn 284 internalReflection = (std::fabs(sinCandida 285 if(internalReflection) { 286 sinRefractionAngle = 1.; 287 cosRefractionAngle = 0.; 288 } else { 289 sinRefractionAngle = sinCandidate; 290 cosRefractionAngle = std::sqrt(1. - sinR 291 } 292 INCL_DEBUG("Refraction parameters initiali 293 << " cosIncidentAngle=" << cosIncid 294 << " sinIncidentAngle=" << sinIncid 295 << " cosRefractionAngle=" << cosRef 296 << " sinRefractionAngle=" << sinRef 297 << " refractionIndexRatio=" << refr 298 << " internalReflection=" << intern 299 } 300 } 301