Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 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 Helsinki Institute of Physics, Finland 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 #include "G4INCLPhaseSpaceGenerator.hh" 39 #include "G4INCLPhaseSpaceKopylov.hh" 40 #include "G4INCLPhaseSpaceRauboldLynch.hh" 41 42 namespace G4INCL { 43 44 namespace { 45 G4ThreadLocal IPhaseSpaceGenerator *thePhaseSpaceGenerator; 46 47 G4ThreadLocal Particle *biasMe; 48 49 /** \brief Actually perform the biasing 50 * 51 * \param particles list of particles to bias 52 * \param pInVec momentum of the particle to be biased before the collision 53 * \param slope the parameter \f$B\f$ in \f$\exp(B\cdot t)\f$ 54 */ 55 void bias(ParticleList &particles, const ThreeVector &pInVec, const G4double slope) { 56 const G4double pIn = pInVec.mag(); 57 const ThreeVector collisionAxis = pInVec/pIn; 58 const ThreeVector pMomVec = biasMe->getMomentum(); 59 const G4double pMom = pMomVec.mag(); 60 if(pMom ==0.) return; 61 const G4double pMomCosAng = pMomVec.dot(collisionAxis)/pMom; 62 const G4double pMomAng = Math::arcCos(pMomCosAng); // Angle between the original axis of the dominant particle and is new one after generate 63 64 // compute the target angle for the biasing 65 // it is drawn from a exp(Bt) distribution 66 const G4double cosAngSlope = 2e-6 * slope * pIn * pMom; 67 const G4double cosAng = 1. + std::log(1. - Random::shoot()*(1.-std::exp(-2.*cosAngSlope)))/cosAngSlope; 68 const G4double ang = Math::arcCos(cosAng); 69 70 // compute the rotation angle 71 const G4double rotationAngle = ang - pMomAng; 72 73 // generate the rotation axis; it is perpendicular to collisionAxis and 74 // pMomVec 75 ThreeVector rotationAxis; 76 if(pMomAng>1E-10) { 77 rotationAxis = collisionAxis.vector(pMomVec); 78 const G4double axisLength = rotationAxis.mag(); 79 const G4double oneOverLength = 1./axisLength; 80 rotationAxis *= oneOverLength; 81 } else { 82 // need to jump through some hoops if collisionAxis is nearly aligned 83 // with pMomVec 84 rotationAxis = collisionAxis.anyOrthogonal(); 85 } 86 87 // apply the rotation 88 particles.rotateMomentum(rotationAngle, rotationAxis); 89 } 90 91 } 92 93 namespace PhaseSpaceGenerator { 94 void generate(const G4double sqrtS, ParticleList &particles) { 95 return thePhaseSpaceGenerator->generate(sqrtS, particles); 96 } 97 98 void generateBiased(const G4double sqrtS, ParticleList &particles, const size_t index, const G4double slope) { 99 // assert(index<particles.size()); 100 // store the incoming momentum of particle[index]; it will be used to 101 // compute t when biasing 102 biasMe = particles[index]; 103 const ThreeVector pInVec = biasMe->getMomentum(); 104 generate(sqrtS, particles); 105 // Extremely rare event try to bias with vector null 106 if(pInVec.mag() != 0.) bias(particles, pInVec, slope); 107 } 108 109 void setPhaseSpaceGenerator(IPhaseSpaceGenerator *g) { 110 thePhaseSpaceGenerator = g; 111 } 112 113 IPhaseSpaceGenerator *getPhaseSpaceGenerator() { 114 return thePhaseSpaceGenerator; 115 } 116 117 void deletePhaseSpaceGenerator() { 118 delete thePhaseSpaceGenerator; 119 thePhaseSpaceGenerator = NULL; 120 } 121 122 void initialize(Config const * const theConfig) { 123 PhaseSpaceGeneratorType psg = theConfig->getPhaseSpaceGeneratorType(); 124 if(psg==RauboldLynchType) 125 setPhaseSpaceGenerator(new PhaseSpaceRauboldLynch); 126 else if(psg==KopylovType) 127 setPhaseSpaceGenerator(new PhaseSpaceKopylov); 128 else 129 setPhaseSpaceGenerator(NULL); 130 } 131 } 132 } 133