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 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland 31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1 35 35 36 #include "globals.hh" 36 #include "globals.hh" 37 37 38 /* 38 /* 39 * Particle.cc 39 * Particle.cc 40 * 40 * 41 * \date Jun 5, 2009 41 * \date Jun 5, 2009 42 * \author Pekka Kaitaniemi 42 * \author Pekka Kaitaniemi 43 */ 43 */ 44 44 45 #include "G4INCLParticle.hh" 45 #include "G4INCLParticle.hh" 46 #include "G4INCLParticleTable.hh" 46 #include "G4INCLParticleTable.hh" 47 47 48 namespace G4INCL { 48 namespace G4INCL { 49 49 50 #ifdef INCLXX_IN_GEANT4_MODE << 50 G4ThreadLocal long Particle::nextID = 1; 51 std::vector<G4double> Particle::INCLBiasVe << 52 #else << 53 G4ThreadLocal std::vector<G4double> Partic << 54 //G4VectorCache<G4double> Particle::INCLBias << 55 #endif << 56 G4ThreadLocal long Particle::nextID = 1; << 57 G4ThreadLocal G4int Particle::nextBiasedColl << 58 51 59 Particle::Particle() 52 Particle::Particle() 60 : theZ(0), theA(0), theS(0), << 53 : theZ(0), theA(0), 61 theParticipantType(TargetSpectator), 54 theParticipantType(TargetSpectator), 62 theType(UnknownParticle), 55 theType(UnknownParticle), 63 theEnergy(0.0), 56 theEnergy(0.0), 64 thePropagationEnergy(&theEnergy), 57 thePropagationEnergy(&theEnergy), 65 theFrozenEnergy(theEnergy), 58 theFrozenEnergy(theEnergy), 66 theMomentum(ThreeVector(0.,0.,0.)), 59 theMomentum(ThreeVector(0.,0.,0.)), 67 thePropagationMomentum(&theMomentum), 60 thePropagationMomentum(&theMomentum), 68 theFrozenMomentum(theMomentum), 61 theFrozenMomentum(theMomentum), 69 thePosition(ThreeVector(0.,0.,0.)), 62 thePosition(ThreeVector(0.,0.,0.)), 70 nCollisions(0), 63 nCollisions(0), 71 nDecays(0), 64 nDecays(0), 72 thePotentialEnergy(0.0), 65 thePotentialEnergy(0.0), 73 rpCorrelated(false), 66 rpCorrelated(false), 74 uncorrelatedMomentum(0.), 67 uncorrelatedMomentum(0.), 75 theParticleBias(1.), << 76 theNKaon(0), << 77 #ifdef INCLXX_IN_GEANT4_MODE << 78 theParentResonancePDGCode(0), << 79 theParentResonanceID(0), << 80 #endif << 81 theHelicity(0.0), 68 theHelicity(0.0), 82 emissionTime(0.0), 69 emissionTime(0.0), 83 outOfWell(false), 70 outOfWell(false), 84 theMass(0.) 71 theMass(0.) 85 { 72 { 86 ID = nextID; 73 ID = nextID; 87 nextID++; 74 nextID++; 88 } 75 } 89 76 90 Particle::Particle(ParticleType t, G4double 77 Particle::Particle(ParticleType t, G4double energy, 91 ThreeVector const &momentum, ThreeVector 78 ThreeVector const &momentum, ThreeVector const &position) 92 : theEnergy(energy), 79 : theEnergy(energy), 93 thePropagationEnergy(&theEnergy), 80 thePropagationEnergy(&theEnergy), 94 theFrozenEnergy(theEnergy), 81 theFrozenEnergy(theEnergy), 95 theMomentum(momentum), 82 theMomentum(momentum), 96 thePropagationMomentum(&theMomentum), 83 thePropagationMomentum(&theMomentum), 97 theFrozenMomentum(theMomentum), 84 theFrozenMomentum(theMomentum), 98 thePosition(position), 85 thePosition(position), 99 nCollisions(0), nDecays(0), 86 nCollisions(0), nDecays(0), 100 thePotentialEnergy(0.), << 87 thePotentialEnergy(0.), 101 rpCorrelated(false), << 88 rpCorrelated(false), 102 uncorrelatedMomentum(theMomentum.mag()), << 89 uncorrelatedMomentum(theMomentum.mag()), 103 theParticleBias(1.), << 90 theHelicity(0.0), 104 theNKaon(0), << 91 emissionTime(0.0), outOfWell(false) 105 #ifdef INCLXX_IN_GEANT4_MODE << 106 theParentResonancePDGCode(0), << 107 theParentResonanceID(0), << 108 #endif << 109 theHelicity(0.0), << 110 emissionTime(0.0), outOfWell(false) << 111 { 92 { 112 theParticipantType = TargetSpectator; 93 theParticipantType = TargetSpectator; 113 ID = nextID; 94 ID = nextID; 114 nextID++; 95 nextID++; 115 if(theEnergy <= 0.0) { 96 if(theEnergy <= 0.0) { 116 INCL_WARN("Particle with energy " << the 97 INCL_WARN("Particle with energy " << theEnergy << " created." << '\n'); 117 } 98 } 118 setType(t); 99 setType(t); 119 setMass(getInvariantMass()); 100 setMass(getInvariantMass()); 120 } 101 } 121 102 122 Particle::Particle(ParticleType t, 103 Particle::Particle(ParticleType t, 123 ThreeVector const &momentum, ThreeVector 104 ThreeVector const &momentum, ThreeVector const &position) 124 : thePropagationEnergy(&theEnergy), 105 : thePropagationEnergy(&theEnergy), 125 theMomentum(momentum), 106 theMomentum(momentum), 126 thePropagationMomentum(&theMomentum), 107 thePropagationMomentum(&theMomentum), 127 theFrozenMomentum(theMomentum), 108 theFrozenMomentum(theMomentum), 128 thePosition(position), 109 thePosition(position), 129 nCollisions(0), nDecays(0), 110 nCollisions(0), nDecays(0), 130 thePotentialEnergy(0.), << 111 thePotentialEnergy(0.), 131 rpCorrelated(false), << 112 rpCorrelated(false), 132 uncorrelatedMomentum(theMomentum.mag()), << 113 uncorrelatedMomentum(theMomentum.mag()), 133 theParticleBias(1.), << 114 theHelicity(0.0), 134 theNKaon(0), << 115 emissionTime(0.0), outOfWell(false) 135 #ifdef INCLXX_IN_GEANT4_MODE << 136 theParentResonancePDGCode(0), << 137 theParentResonanceID(0), << 138 #endif << 139 theHelicity(0.0), << 140 emissionTime(0.0), outOfWell(false) << 141 { 116 { 142 theParticipantType = TargetSpectator; 117 theParticipantType = TargetSpectator; 143 ID = nextID; 118 ID = nextID; 144 nextID++; 119 nextID++; 145 setType(t); 120 setType(t); 146 if( isResonance() ) { 121 if( isResonance() ) { 147 INCL_ERROR("Cannot create resonance with 122 INCL_ERROR("Cannot create resonance without specifying its momentum four-vector." << '\n'); 148 } 123 } 149 G4double energy = std::sqrt(theMomentum.ma 124 G4double energy = std::sqrt(theMomentum.mag2() + theMass*theMass); 150 theEnergy = energy; 125 theEnergy = energy; 151 theFrozenEnergy = theEnergy; 126 theFrozenEnergy = theEnergy; 152 } 127 } 153 128 154 const ThreeVector &Particle::adjustMomentumF 129 const ThreeVector &Particle::adjustMomentumFromEnergy() { 155 const G4double p2 = theMomentum.mag2(); 130 const G4double p2 = theMomentum.mag2(); 156 G4double newp2 = theEnergy*theEnergy - the 131 G4double newp2 = theEnergy*theEnergy - theMass*theMass; 157 if( newp2<0.0 ) { 132 if( newp2<0.0 ) { 158 INCL_ERROR("Particle has E^2 < m^2." << 133 INCL_ERROR("Particle has E^2 < m^2." << '\n' << print()); 159 newp2 = 0.0; 134 newp2 = 0.0; 160 theEnergy = theMass; 135 theEnergy = theMass; 161 } 136 } 162 137 163 theMomentum *= std::sqrt(newp2/p2); 138 theMomentum *= std::sqrt(newp2/p2); 164 return theMomentum; 139 return theMomentum; 165 } 140 } 166 141 167 G4double Particle::adjustEnergyFromMomentum( 142 G4double Particle::adjustEnergyFromMomentum() { 168 theEnergy = std::sqrt(theMomentum.mag2() + 143 theEnergy = std::sqrt(theMomentum.mag2() + theMass*theMass); 169 return theEnergy; 144 return theEnergy; 170 } 145 } 171 146 172 void ParticleList::rotatePositionAndMomentum 147 void ParticleList::rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) const { 173 for(const_iterator i=begin(), e=end(); i!= 148 for(const_iterator i=begin(), e=end(); i!=e; ++i) { 174 (*i)->rotatePositionAndMomentum(angle, a 149 (*i)->rotatePositionAndMomentum(angle, axis); 175 } 150 } 176 } 151 } 177 152 178 void ParticleList::rotatePosition(const G4do 153 void ParticleList::rotatePosition(const G4double angle, const ThreeVector &axis) const { 179 for(const_iterator i=begin(), e=end(); i!= 154 for(const_iterator i=begin(), e=end(); i!=e; ++i) { 180 (*i)->rotatePosition(angle, axis); 155 (*i)->rotatePosition(angle, axis); 181 } 156 } 182 } 157 } 183 158 184 void ParticleList::rotateMomentum(const G4do 159 void ParticleList::rotateMomentum(const G4double angle, const ThreeVector &axis) const { 185 for(const_iterator i=begin(), e=end(); i!= 160 for(const_iterator i=begin(), e=end(); i!=e; ++i) { 186 (*i)->rotateMomentum(angle, axis); 161 (*i)->rotateMomentum(angle, axis); 187 } 162 } 188 } 163 } 189 164 190 void ParticleList::boost(const ThreeVector & 165 void ParticleList::boost(const ThreeVector &b) const { 191 for(const_iterator i=begin(), e=end(); i!= 166 for(const_iterator i=begin(), e=end(); i!=e; ++i) { 192 (*i)->boost(b); 167 (*i)->boost(b); 193 } 168 } 194 } << 195 << 196 G4double ParticleList::getParticleListBias() << 197 if(G4int((*this).size())==0) return 1.; << 198 std::vector<G4int> MergedVector; << 199 for(ParticleIter i = (*this).begin(), e = << 200 MergedVector = Particle::MergeVectorBi << 201 } << 202 return Particle::getBiasFromVector(std::mo << 203 } << 204 << 205 std::vector<G4int> ParticleList::getParticle << 206 std::vector<G4int> MergedVector; << 207 if(G4int((*this).size())==0) return Merged << 208 for(ParticleIter i = (*this).begin(), e = << 209 MergedVector = Particle::MergeVectorBi << 210 } << 211 return MergedVector; << 212 } << 213 << 214 void Particle::FillINCLBiasVector(G4double n << 215 // assert(G4int(Particle::INCLBiasVector.size( << 216 //assert(G4int(Particle::INCLBiasVector.Si << 217 // assert(std::fabs(newBias - 1.) > 1E-6); << 218 Particle::INCLBiasVector.push_back(newBias); << 219 //Particle::INCLBiasVector.Push_back(newBias << 220 Particle::nextBiasedCollisionID++; << 221 } << 222 << 223 G4double Particle::getBiasFromVector(std::ve << 224 if(VectorBias.empty()) return 1.; << 225 << 226 G4double ParticleBias = 1.; << 227 << 228 for(G4int i=0; i<G4int(VectorBias.size()); << 229 ParticleBias *= Particle::INCLBiasVect << 230 } << 231 << 232 return ParticleBias; << 233 } << 234 << 235 std::vector<G4int> Particle::MergeVectorBias << 236 std::vector<G4int> MergedVectorBias; << 237 std::vector<G4int> VectorBias1 = p1->getBi << 238 std::vector<G4int> VectorBias2 = p2->getBi << 239 G4int i = 0; << 240 G4int j = 0; << 241 if(VectorBias1.size()==0 && VectorBias2.si << 242 else if(VectorBias1.size()==0) return Vect << 243 else if(VectorBias2.size()==0) return Vect << 244 << 245 while(i < G4int(VectorBias1.size()) || j < << 246 if(VectorBias1[i]==VectorBias2[j]){ << 247 MergedVectorBias.push_back(VectorB << 248 i++; << 249 j++; << 250 if(i == G4int(VectorBias1.size())) << 251 for(;j<G4int(VectorBias2.size( << 252 } << 253 else if(j == G4int(VectorBias2.siz << 254 for(;i<G4int(VectorBias1.size( << 255 } << 256 } else if(VectorBias1[i]<VectorBias2[j << 257 MergedVectorBias.push_back(VectorB << 258 i++; << 259 if(i == G4int(VectorBias1.size())) << 260 for(;j<G4int(VectorBias2.size( << 261 } << 262 } << 263 else { << 264 MergedVectorBias.push_back(VectorB << 265 j++; << 266 if(j == G4int(VectorBias2.size())) << 267 for(;i<G4int(VectorBias1.size( << 268 } << 269 } << 270 } << 271 return MergedVectorBias; << 272 } << 273 << 274 std::vector<G4int> Particle::MergeVectorBias << 275 std::vector<G4int> MergedVectorBias; << 276 std::vector<G4int> VectorBias = p2->getBia << 277 G4int i = 0; << 278 G4int j = 0; << 279 if(p1.size()==0 && VectorBias.size()==0) r << 280 else if(p1.size()==0) return VectorBias; << 281 else if(VectorBias.size()==0) return p1; << 282 << 283 while(i < G4int(p1.size()) || j < G4int(Ve << 284 if(p1[i]==VectorBias[j]){ << 285 MergedVectorBias.push_back(p1[i]); << 286 i++; << 287 j++; << 288 if(i == G4int(p1.size())){ << 289 for(;j<G4int(VectorBias.size() << 290 } << 291 else if(j == G4int(VectorBias.size << 292 for(;i<G4int(p1.size());i++) M << 293 } << 294 } else if(p1[i]<VectorBias[j]){ << 295 MergedVectorBias.push_back(p1[i]); << 296 i++; << 297 if(i == G4int(p1.size())){ << 298 for(;j<G4int(VectorBias.size() << 299 } << 300 } << 301 else { << 302 MergedVectorBias.push_back(VectorB << 303 j++; << 304 if(j == G4int(VectorBias.size())){ << 305 for(;i<G4int(p1.size());i++) M << 306 } << 307 } << 308 } << 309 return MergedVectorBias; << 310 } << 311 << 312 G4double Particle::getTotalBias() { << 313 G4double TotalBias = 1.; << 314 for(G4int i=0; i<G4int(INCLBiasVector.si << 315 return TotalBias; << 316 } << 317 << 318 void Particle::setINCLBiasVector(std::vector << 319 Particle::INCLBiasVector = std::move(New << 320 } 169 } 321 } 170 } 322 171