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 #include "G4INCLPhaseSpaceRauboldLynch.hh" 38 #include "G4INCLPhaseSpaceRauboldLynch.hh" 39 #include "G4INCLRandom.hh" 39 #include "G4INCLRandom.hh" 40 #include "G4INCLKinematicsUtils.hh" 40 #include "G4INCLKinematicsUtils.hh" 41 #include <algorithm> 41 #include <algorithm> 42 #include <numeric> 42 #include <numeric> 43 #include <functional> 43 #include <functional> 44 44 45 namespace G4INCL { 45 namespace G4INCL { 46 46 47 const G4double PhaseSpaceRauboldLynch::wMaxM 47 const G4double PhaseSpaceRauboldLynch::wMaxMasslessX[PhaseSpaceRauboldLynch::wMaxNE] = { 48 0.01, 48 0.01, 49 0.017433288222, 49 0.017433288222, 50 0.0303919538231, 50 0.0303919538231, 51 0.0529831690628, 51 0.0529831690628, 52 0.0923670857187, 52 0.0923670857187, 53 0.161026202756, 53 0.161026202756, 54 0.280721620394, 54 0.280721620394, 55 0.489390091848, 55 0.489390091848, 56 0.853167852417, 56 0.853167852417, 57 1.48735210729, 57 1.48735210729, 58 2.5929437974, 58 2.5929437974, 59 4.52035365636, 59 4.52035365636, 60 7.88046281567, 60 7.88046281567, 61 13.7382379588, 61 13.7382379588, 62 23.9502661999, 62 23.9502661999, 63 41.7531893656, 63 41.7531893656, 64 72.7895384398, 64 72.7895384398, 65 126.896100317, 65 126.896100317, 66 221.221629107, 66 221.221629107, 67 385.662042116, 67 385.662042116, 68 672.33575365, 68 672.33575365, 69 1172.10229753, 69 1172.10229753, 70 2043.35971786, 70 2043.35971786, 71 3562.24789026, 71 3562.24789026, 72 6210.16941892, 72 6210.16941892, 73 10826.3673387, 73 10826.3673387, 74 18873.9182214, 74 18873.9182214, 75 32903.4456231, 75 32903.4456231, 76 57361.5251045, 76 57361.5251045, 77 100000.0 77 100000.0 78 }; 78 }; 79 79 80 const G4double PhaseSpaceRauboldLynch::wMaxM 80 const G4double PhaseSpaceRauboldLynch::wMaxMasslessY[PhaseSpaceRauboldLynch::wMaxNE] = { 81 -4.69180212056, 81 -4.69180212056, 82 -4.13600582212, 82 -4.13600582212, 83 -3.58020940928, 83 -3.58020940928, 84 -3.02441303018, 84 -3.02441303018, 85 -2.46861663471, 85 -2.46861663471, 86 -1.91282025644, 86 -1.91282025644, 87 -1.35702379603, 87 -1.35702379603, 88 -0.801227342882, 88 -0.801227342882, 89 -0.245430866403, 89 -0.245430866403, 90 0.310365269464, 90 0.310365269464, 91 0.866161720461, 91 0.866161720461, 92 1.42195839972, 92 1.42195839972, 93 1.97775459763, 93 1.97775459763, 94 2.5335510251, 94 2.5335510251, 95 3.08934734235, 95 3.08934734235, 96 3.64514378357, 96 3.64514378357, 97 4.20094013304, 97 4.20094013304, 98 4.75673663205, 98 4.75673663205, 99 5.3125329676, 99 5.3125329676, 100 5.86832950014, 100 5.86832950014, 101 6.42412601468, 101 6.42412601468, 102 6.97992197797, 102 6.97992197797, 103 7.53571856765, 103 7.53571856765, 104 8.09151503031, 104 8.09151503031, 105 8.64731144398, 105 8.64731144398, 106 9.20310774148, 106 9.20310774148, 107 9.7589041009, 107 9.7589041009, 108 10.314700625, 108 10.314700625, 109 10.8704971207, 109 10.8704971207, 110 11.4262935304 110 11.4262935304 111 }; 111 }; 112 112 113 const G4double PhaseSpaceRauboldLynch::wMaxC 113 const G4double PhaseSpaceRauboldLynch::wMaxCorrectionX[PhaseSpaceRauboldLynch::wMaxNE] = { 114 8.77396538453e-07, 114 8.77396538453e-07, 115 1.52959067398e-06, 115 1.52959067398e-06, 116 2.66657950812e-06, 116 2.66657950812e-06, 117 4.6487249132e-06, 117 4.6487249132e-06, 118 8.10425612766e-06, 118 8.10425612766e-06, 119 1.41283832898e-05, 119 1.41283832898e-05, 120 2.46304178003e-05, 120 2.46304178003e-05, 121 4.2938917254e-05, 121 4.2938917254e-05, 122 7.4856652043e-05, 122 7.4856652043e-05, 123 0.00013049975904, 123 0.00013049975904, 124 0.000227503991225, 124 0.000227503991225, 125 0.000396614265067, 125 0.000396614265067, 126 0.000691429079588, 126 0.000691429079588, 127 0.00120538824295, 127 0.00120538824295, 128 0.00210138806588, 128 0.00210138806588, 129 0.00366341038188, 129 0.00366341038188, 130 0.00638652890627, 130 0.00638652890627, 131 0.0111338199161, 131 0.0111338199161, 132 0.0194099091609, 132 0.0194099091609, 133 0.0338378540766, 133 0.0338378540766, 134 0.0589905062931, 134 0.0589905062931, 135 0.102839849857, 135 0.102839849857, 136 0.179283674326, 136 0.179283674326, 137 0.312550396803, 137 0.312550396803, 138 0.544878115136, 138 0.544878115136, 139 0.949901722703, 139 0.949901722703, 140 1.65599105145, 140 1.65599105145, 141 2.88693692929, 141 2.88693692929, 142 5.03288035671, 142 5.03288035671, 143 8.77396538453 143 8.77396538453 144 }; 144 }; 145 const G4double PhaseSpaceRauboldLynch::wMaxC 145 const G4double PhaseSpaceRauboldLynch::wMaxCorrectionY[PhaseSpaceRauboldLynch::wMaxNE] = { 146 7.28682764024, 146 7.28682764024, 147 7.0089298602, 147 7.0089298602, 148 6.7310326046, 148 6.7310326046, 149 6.45313436085, 149 6.45313436085, 150 6.17523713298, 150 6.17523713298, 151 5.89734080335, 151 5.89734080335, 152 5.61944580818, 152 5.61944580818, 153 5.34155326611, 153 5.34155326611, 154 5.06366530628, 154 5.06366530628, 155 4.78578514804, 155 4.78578514804, 156 4.50791742491, 156 4.50791742491, 157 4.23007259554, 157 4.23007259554, 158 3.97566018117, 158 3.97566018117, 159 3.72116551935, 159 3.72116551935, 160 3.44355099732, 160 3.44355099732, 161 3.1746329047, 161 3.1746329047, 162 2.89761210229, 162 2.89761210229, 163 2.62143335535, 163 2.62143335535, 164 2.34649707964, 164 2.34649707964, 165 2.07366126445, 165 2.07366126445, 166 1.8043078447, 166 1.8043078447, 167 1.54056992953, 167 1.54056992953, 168 1.27913996411, 168 1.27913996411, 169 0.981358535322, 169 0.981358535322, 170 0.73694629682, 170 0.73694629682, 171 0.587421056631, 171 0.587421056631, 172 0.417178268911, 172 0.417178268911, 173 0.280881750176, 173 0.280881750176, 174 0.187480311508, 174 0.187480311508, 175 0.116321254846 175 0.116321254846 176 }; 176 }; 177 177 178 const G4double PhaseSpaceRauboldLynch::wMaxI 178 const G4double PhaseSpaceRauboldLynch::wMaxInterpolationMargin = std::log(1.5); 179 179 180 PhaseSpaceRauboldLynch::PhaseSpaceRauboldLyn 180 PhaseSpaceRauboldLynch::PhaseSpaceRauboldLynch() : 181 nParticles(0), 181 nParticles(0), 182 sqrtS(0.), 182 sqrtS(0.), 183 availableEnergy(0.), 183 availableEnergy(0.), 184 maxGeneratedWeight(0.) 184 maxGeneratedWeight(0.) 185 { 185 { 186 std::vector<G4double> wMaxMasslessXV(wMaxM 186 std::vector<G4double> wMaxMasslessXV(wMaxMasslessX, wMaxMasslessX + wMaxNE); 187 std::vector<G4double> wMaxMasslessYV(wMaxM 187 std::vector<G4double> wMaxMasslessYV(wMaxMasslessY, wMaxMasslessY + wMaxNE); 188 wMaxMassless = new InterpolationTable(wMax 188 wMaxMassless = new InterpolationTable(wMaxMasslessXV, wMaxMasslessYV); 189 std::vector<G4double> wMaxCorrectionXV(wMa 189 std::vector<G4double> wMaxCorrectionXV(wMaxCorrectionX, wMaxCorrectionX + wMaxNE); 190 std::vector<G4double> wMaxCorrectionYV(wMa 190 std::vector<G4double> wMaxCorrectionYV(wMaxCorrectionY, wMaxCorrectionY + wMaxNE); 191 wMaxCorrection = new InterpolationTable(wM 191 wMaxCorrection = new InterpolationTable(wMaxCorrectionXV, wMaxCorrectionYV); 192 192 193 // initialize the precalculated coefficien 193 // initialize the precalculated coefficients 194 prelog[0] = 0.; 194 prelog[0] = 0.; 195 for(size_t i=1; i<wMaxNP; ++i) { 195 for(size_t i=1; i<wMaxNP; ++i) { 196 prelog[i] = -std::log(G4double(i)); 196 prelog[i] = -std::log(G4double(i)); 197 } 197 } 198 } 198 } 199 199 200 PhaseSpaceRauboldLynch::~PhaseSpaceRauboldLy 200 PhaseSpaceRauboldLynch::~PhaseSpaceRauboldLynch() { 201 delete wMaxMassless; 201 delete wMaxMassless; 202 delete wMaxCorrection; 202 delete wMaxCorrection; 203 } 203 } 204 204 205 void PhaseSpaceRauboldLynch::generate(const 205 void PhaseSpaceRauboldLynch::generate(const G4double sqs, ParticleList &particles) { 206 maxGeneratedWeight = 0.; 206 maxGeneratedWeight = 0.; 207 207 208 sqrtS = sqs; 208 sqrtS = sqs; 209 209 210 // initialize the structures containing th 210 // initialize the structures containing the particle masses 211 initialize(particles); 211 initialize(particles); 212 212 213 // initialize the maximum weight 213 // initialize the maximum weight 214 const G4double weightMax = computeMaximumW 214 const G4double weightMax = computeMaximumWeightParam(); 215 215 216 const G4int maxIter = 500; 216 const G4int maxIter = 500; 217 G4int iter = 0; 217 G4int iter = 0; 218 G4double weight, r; 218 G4double weight, r; 219 do { 219 do { 220 weight = computeWeight(); 220 weight = computeWeight(); 221 maxGeneratedWeight = std::max(weight, ma 221 maxGeneratedWeight = std::max(weight, maxGeneratedWeight); 222 // assert(weight<=weightMax); 222 // assert(weight<=weightMax); 223 r = Random::shoot(); 223 r = Random::shoot(); 224 } while(++iter<maxIter && r*weightMax>weig 224 } while(++iter<maxIter && r*weightMax>weight); /* Loop checking, 10.07.2015, D.Mancusi */ 225 225 226 #ifndef INCLXX_IN_GEANT4_MODE 226 #ifndef INCLXX_IN_GEANT4_MODE 227 if(iter>=maxIter) { 227 if(iter>=maxIter) { 228 INCL_WARN("Number of tries exceeded in P 228 INCL_WARN("Number of tries exceeded in PhaseSpaceRauboldLynch::generate()\n" 229 << "nParticles=" << nParticles 229 << "nParticles=" << nParticles << ", sqrtS=" << sqrtS << '\n'); 230 } 230 } 231 #endif 231 #endif 232 232 233 generateEvent(particles); 233 generateEvent(particles); 234 } 234 } 235 235 236 void PhaseSpaceRauboldLynch::initialize(Part 236 void PhaseSpaceRauboldLynch::initialize(ParticleList &particles) { 237 nParticles = particles.size(); 237 nParticles = particles.size(); 238 // assert(nParticles>2); 238 // assert(nParticles>2); 239 239 240 // masses and sum of masses 240 // masses and sum of masses 241 masses.resize(nParticles); 241 masses.resize(nParticles); 242 sumMasses.resize(nParticles); 242 sumMasses.resize(nParticles); 243 std::transform(particles.begin(), particle << 243 std::transform(particles.begin(), particles.end(), masses.begin(), std::mem_fun(&Particle::getMass)); 244 std::partial_sum(masses.begin(), masses.en 244 std::partial_sum(masses.begin(), masses.end(), sumMasses.begin()); 245 245 246 // sanity check 246 // sanity check 247 availableEnergy = sqrtS-sumMasses[nPartic 247 availableEnergy = sqrtS-sumMasses[nParticles-1]; 248 // assert(availableEnergy>-1.e-5); 248 // assert(availableEnergy>-1.e-5); 249 if(availableEnergy<0.) 249 if(availableEnergy<0.) 250 availableEnergy = 0.; 250 availableEnergy = 0.; 251 251 252 rnd.resize(nParticles); 252 rnd.resize(nParticles); 253 invariantMasses.resize(nParticles); 253 invariantMasses.resize(nParticles); 254 momentaCM.resize(nParticles-1); 254 momentaCM.resize(nParticles-1); 255 } 255 } 256 256 257 G4double PhaseSpaceRauboldLynch::computeMaxi 257 G4double PhaseSpaceRauboldLynch::computeMaximumWeightNaive() { 258 // compute the maximum weight 258 // compute the maximum weight 259 G4double eMMax = sqrtS + masses[0]; 259 G4double eMMax = sqrtS + masses[0]; 260 G4double eMMin = 0.; 260 G4double eMMin = 0.; 261 G4double wMax = 1.; 261 G4double wMax = 1.; 262 for(size_t i=1; i<nParticles; i++) { 262 for(size_t i=1; i<nParticles; i++) { 263 eMMin += masses[i-1]; 263 eMMin += masses[i-1]; 264 eMMax += masses[i]; 264 eMMax += masses[i]; 265 wMax *= KinematicsUtils::momentumInCM(eM 265 wMax *= KinematicsUtils::momentumInCM(eMMax, eMMin, masses[i]); 266 } 266 } 267 return wMax; 267 return wMax; 268 } 268 } 269 269 270 G4double PhaseSpaceRauboldLynch::computeMaxi 270 G4double PhaseSpaceRauboldLynch::computeMaximumWeightParam() { 271 #ifndef INCLXX_IN_GEANT4_MODE 271 #ifndef INCLXX_IN_GEANT4_MODE 272 if(nParticles>=wMaxNP) { 272 if(nParticles>=wMaxNP) { 273 INCL_WARN("The requested number of parti 273 INCL_WARN("The requested number of particles (" << nParticles << ") requires extrapolation the tables in PhaseSpaceRauboldLynch. YMMV." << '\n'); 274 } 274 } 275 275 276 if(availableEnergy>wMaxMasslessX[wMaxNE-1] 276 if(availableEnergy>wMaxMasslessX[wMaxNE-1] || availableEnergy<wMaxMasslessX[0] ) { 277 INCL_WARN("The requested available energ 277 INCL_WARN("The requested available energy (" << availableEnergy << " MeV) requires extrapolation the tables in PhaseSpaceRauboldLynch. YMMV." << '\n'); 278 } 278 } 279 #endif 279 #endif 280 280 281 // compute the maximum weight 281 // compute the maximum weight 282 const G4double logMassless = ((*wMaxMassle 282 const G4double logMassless = ((*wMaxMassless)(availableEnergy)+prelog[nParticles])*(nParticles-1); 283 const G4double reducedSqrtS = availableEne 283 const G4double reducedSqrtS = availableEnergy/sumMasses[nParticles-1]; 284 const G4double correction = (*wMaxCorrecti 284 const G4double correction = (*wMaxCorrection)(reducedSqrtS); 285 const G4double wMax = std::exp(correction* 285 const G4double wMax = std::exp(correction*G4double(nParticles-1) + logMassless + wMaxInterpolationMargin); 286 if(wMax>0.) 286 if(wMax>0.) 287 return wMax; 287 return wMax; 288 else 288 else 289 return computeMaximumWeightNaive(); 289 return computeMaximumWeightNaive(); 290 } 290 } 291 291 292 G4double PhaseSpaceRauboldLynch::computeWeig 292 G4double PhaseSpaceRauboldLynch::computeWeight() { 293 // Generate nParticles-2 sorted random num 293 // Generate nParticles-2 sorted random numbers, bracketed by 0 and 1 294 rnd[0] = 0.; 294 rnd[0] = 0.; 295 for(size_t i=1; i<nParticles-1; ++i) 295 for(size_t i=1; i<nParticles-1; ++i) 296 rnd[i] = Random::shoot(); 296 rnd[i] = Random::shoot(); 297 rnd[nParticles-1] = 1.; 297 rnd[nParticles-1] = 1.; 298 std::sort(rnd.begin()+1, rnd.begin()+nPart 298 std::sort(rnd.begin()+1, rnd.begin()+nParticles-1); 299 << 299 300 // invariant masses 300 // invariant masses 301 for(size_t i=0; i<nParticles; ++i) 301 for(size_t i=0; i<nParticles; ++i) 302 invariantMasses[i] = rnd[i]*availableEne 302 invariantMasses[i] = rnd[i]*availableEnergy + sumMasses[i]; 303 303 304 // compute the CM momenta and the weight f 304 // compute the CM momenta and the weight for the current event 305 G4double weight=KinematicsUtils::momentumI << 305 G4double weight=KinematicsUtils::momentumInCM(invariantMasses[1], invariantMasses[0], masses[1]);; 306 momentaCM[0] = weight; 306 momentaCM[0] = weight; 307 for(size_t i=1; i<nParticles-1; ++i) { 307 for(size_t i=1; i<nParticles-1; ++i) { 308 G4double momentumCM; << 308 const G4double momentumCM = KinematicsUtils::momentumInCM(invariantMasses[i+1], invariantMasses[i], masses[i+1]); 309 // assert(invariantMasses[i+1]-invariantMasses << 310 if(invariantMasses[i+1]-invariantMasses[ << 311 else momentumCM = KinematicsUtils::momen << 312 momentaCM[i] = momentumCM; 309 momentaCM[i] = momentumCM; 313 weight *= momentumCM; 310 weight *= momentumCM; 314 } 311 } >> 312 >> 313 // assert(weight==weight); 315 314 316 return weight; 315 return weight; 317 } 316 } 318 317 319 void PhaseSpaceRauboldLynch::generateEvent(P 318 void PhaseSpaceRauboldLynch::generateEvent(ParticleList &particles) { 320 Particle *p = particles[0]; 319 Particle *p = particles[0]; 321 ThreeVector momentum = Random::normVector( 320 ThreeVector momentum = Random::normVector(momentaCM[0]); 322 p->setMomentum(momentum); 321 p->setMomentum(momentum); 323 p->adjustEnergyFromMomentum(); 322 p->adjustEnergyFromMomentum(); 324 323 325 ThreeVector boostV; 324 ThreeVector boostV; 326 325 327 for(size_t i=1; i<nParticles; ++i) { 326 for(size_t i=1; i<nParticles; ++i) { 328 p = particles[i]; 327 p = particles[i]; 329 p->setMomentum(-momentum); 328 p->setMomentum(-momentum); 330 p->adjustEnergyFromMomentum(); 329 p->adjustEnergyFromMomentum(); 331 330 332 if(i==nParticles-1) 331 if(i==nParticles-1) 333 break; 332 break; 334 333 335 momentum = Random::normVector(momentaCM[ 334 momentum = Random::normVector(momentaCM[i]); 336 335 337 const G4double iM = invariantMasses[i]; 336 const G4double iM = invariantMasses[i]; 338 const G4double recoilE = std::sqrt(momen 337 const G4double recoilE = std::sqrt(momentum.mag2() + iM*iM); 339 boostV = -momentum/recoilE; 338 boostV = -momentum/recoilE; 340 for(size_t j=0; j<=i; ++j) 339 for(size_t j=0; j<=i; ++j) 341 particles[j]->boost(boostV); 340 particles[j]->boost(boostV); 342 } 341 } 343 } 342 } 344 343 345 G4double PhaseSpaceRauboldLynch::getMaxGener 344 G4double PhaseSpaceRauboldLynch::getMaxGeneratedWeight() const { 346 return maxGeneratedWeight; 345 return maxGeneratedWeight; 347 } 346 } 348 } 347 } 349 348