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