- Timestamp:
- Nov 5, 2010, 3:45:55 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/cascade/cascade/src/G4BigBanger.cc
r1337 r1340 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BigBanger.cc,v 1. 30 2010/06/25 09:43:58 gunterExp $27 // Geant4 tag: $Name: geant4-09-04-beta-01$26 // $Id: G4BigBanger.cc,v 1.40 2010/09/28 20:15:00 mkelsey Exp $ 27 // Geant4 tag: $Name: hadr-casc-V09-03-85 $ 28 28 // 29 29 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly … … 34 34 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide() 35 35 // 20100517 M. Kelsey -- Inherit from common base class, clean up code 36 // 20100628 M. Kelsey -- Use old "bindingEnergy" fn as wrapper, add balance 37 // checking after bang. 38 // 20100630 M. Kelsey -- Just do simple boost for target, instead of using 39 // G4LorentzConverter with dummy bullet. 40 // 20100701 M. Kelsey -- Re-throw momentum list, not just angles! 41 // 20100714 M. Kelsey -- Move conservation checking to base class 42 // 20100726 M. Kelsey -- Move std::vector<> buffer to .hh file 43 // 20100923 M. Kelsey -- Migrate to integer A and Z 36 44 37 45 #include "G4BigBanger.hh" … … 41 49 #include "G4InuclSpecialFunctions.hh" 42 50 #include "G4ParticleLargerEkin.hh" 43 #include "G4LorentzConvertor.hh"44 #include "G4HadTmpUtil.hh"45 #include "G4NucleiProperties.hh"46 51 #include <algorithm> 47 52 … … 50 55 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator; 51 56 52 G4BigBanger::G4BigBanger() : G4 VCascadeCollider("G4BigBanger") {}57 G4BigBanger::G4BigBanger() : G4CascadeColliderBase("G4BigBanger") {} 53 58 54 59 void … … 56 61 G4CollisionOutput& output) { 57 62 58 if (verboseLevel > 3) { 59 G4cout << " >>> G4BigBanger::collide" << G4endl; 60 } 63 if (verboseLevel) G4cout << " >>> G4BigBanger::collide" << G4endl; 61 64 62 65 // primitive explosion model A -> nucleons to prevent too exotic evaporation 63 64 const G4double small_ekin = 1.0e-6;65 66 G4LorentzVector totscm;67 G4LorentzVector totlab;68 66 69 67 G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target); … … 73 71 } 74 72 75 G4double A = nuclei_target->getA(); 76 G4double Z = nuclei_target->getZ(); 73 G4int A = nuclei_target->getA(); 74 G4int Z = nuclei_target->getZ(); 75 77 76 G4LorentzVector PEX = nuclei_target->getMomentum(); 78 77 G4double EEXS = nuclei_target->getExitationEnergy(); 79 G4InuclElementaryParticle dummy(small_ekin, 1); 80 G4LorentzConvertor toTheNucleiSystemRestFrame; 81 82 toTheNucleiSystemRestFrame.setBullet(dummy); 83 toTheNucleiSystemRestFrame.setTarget(nuclei_target); 84 toTheNucleiSystemRestFrame.toTheTargetRestFrame(); 85 86 G4double etot = (EEXS - G4NucleiProperties::GetBindingEnergy(G4lrint(A), G4lrint(Z) ) ) * MeV/GeV; // To Bertini units 78 79 G4ThreeVector toTheLabFrame = PEX.boostVector(); // From rest to lab 80 81 // This "should" be difference between E-target and sum of m(nucleons) 82 G4double etot = (EEXS - bindingEnergy(A,Z)) * MeV/GeV; // To Bertini units 87 83 if (etot < 0.0) etot = 0.0; 88 84 … … 90 86 G4cout << " BigBanger: target " << G4endl; 91 87 nuclei_target->printParticle(); 92 G4cout << " BigBanger: a " << A << " z " << Z << " eexs " << EEXS << " etot " << 93 etot << " nm " << nuclei_target->getMass() << G4endl; 94 } 95 88 G4cout << " etot " << etot << G4endl; 89 } 90 91 if (verboseLevel > 3) { 92 G4LorentzVector PEXrest = PEX; 93 PEXrest.boost(-toTheLabFrame); 94 G4cout << " target rest frame: px " << PEXrest.px() << " py " 95 << PEXrest.py() << " pz " << PEXrest.pz() << " E " << PEXrest.e() 96 << G4endl; 97 } 98 96 99 generateBangInSCM(etot, A, Z); 97 100 … … 102 105 } 103 106 104 if(!particles.empty()) { // convert back to Lab 105 particleIterator ipart; 106 107 for(ipart = particles.begin(); ipart != particles.end(); ipart++) { 108 if (verboseLevel > 2) { 109 totscm += ipart->getMomentum(); 110 } 111 G4LorentzVector mom = 112 toTheNucleiSystemRestFrame.backToTheLab(ipart->getMomentum()); 113 ipart->setMomentum(mom); 114 115 if (verboseLevel > 2) { 116 mom = ipart->getMomentum(); 117 totlab += mom; 118 } 119 } 120 121 std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin()); 122 123 if (verboseLevel > 2) { 124 G4cout << " In SCM: total outgoing momentum " << G4endl 125 << " E " << totscm.e() << " px " << totscm.x() 126 << " py " << totscm.y() << " pz " << totscm.z() << G4endl; 127 G4cout << " In Lab: mom cons " << G4endl 128 << " E " << PEX.e() + 0.001 * EEXS - totlab.e() 129 << " px " << PEX.x() - totlab.x() 130 << " py " << PEX.y() - totlab.y() 131 << " pz " << PEX.z() - totlab.z() << G4endl; 132 } 107 if (particles.empty()) { // No bang! Don't know why... 108 G4cerr << " >>> G4BigBanger unable to process fragment " 109 << nuclei_target->getDefinition()->GetParticleName() << G4endl; 110 111 // FIXME: This will violate baryon number, momentum, energy, etc. 112 return; 113 } 114 115 // convert back to Lab 116 G4LorentzVector totscm; 117 G4LorentzVector totlab; 118 119 if (verboseLevel > 2) G4cout << " BigBanger: boosting to lab" << G4endl; 120 121 particleIterator ipart; 122 for(ipart = particles.begin(); ipart != particles.end(); ipart++) { 123 G4LorentzVector mom = ipart->getMomentum(); 124 if (verboseLevel > 2) totscm += mom; 125 126 mom.boost(toTheLabFrame); 127 if (verboseLevel > 2) totlab += mom; 128 129 ipart->setMomentum(mom); 130 if (verboseLevel > 2) ipart->printParticle(); 131 } 132 133 std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin()); 134 135 validateOutput(0, target, particles); // Checks <vector> directly 136 137 if (verboseLevel > 2) { 138 G4cout << " In SCM: total outgoing momentum " << G4endl 139 << " E " << totscm.e() << " px " << totscm.x() 140 << " py " << totscm.y() << " pz " << totscm.z() << G4endl; 141 G4cout << " In Lab: mom cons " << G4endl 142 << " E " << PEX.e() - totlab.e() // PEX now includes EEXS 143 << " px " << PEX.x() - totlab.x() 144 << " py " << PEX.y() - totlab.y() 145 << " pz " << PEX.z() - totlab.z() << G4endl; 133 146 } 134 147 135 148 output.addOutgoingParticles(particles); 136 137 return;138 149 } 139 150 140 void G4BigBanger::generateBangInSCM(G4double etot, G4 double a, G4doublez) {151 void G4BigBanger::generateBangInSCM(G4double etot, G4int a, G4int z) { 141 152 if (verboseLevel > 3) { 142 153 G4cout << " >>> G4BigBanger::generateBangInSCM" << G4endl; 143 154 } 144 155 145 // Proton and neutron masses146 const G4double mp = G4InuclElementaryParticle::getParticleMass(1);147 const G4double mn = G4InuclElementaryParticle::getParticleMass(2);148 149 156 const G4double ang_cut = 0.9999; 150 157 const G4int itry_max = 1000; 151 158 152 G4int ia = int(a + 0.1); 153 G4int iz = int(z + 0.1); 154 155 if (verboseLevel > 2) { 156 G4cout << " ia " << ia << " iz " << iz << G4endl; 159 if (verboseLevel > 2) { 160 G4cout << " a " << a << " z " << z << G4endl; 157 161 } 158 162 159 163 particles.clear(); // Reset output vector before filling 160 164 161 if(ia == 1) { 162 // abnormal situation 163 G4double m = iz > 0 ? mp : mn; 164 G4double pmod = std::sqrt((etot + 2.0 * m) * etot); 165 G4LorentzVector mom = generateWithRandomAngles(pmod, m); 166 167 G4int knd = iz > 0 ? 1 : 2; 168 169 particles.push_back(G4InuclElementaryParticle(mom, knd, 8)); // modelId included 170 165 if (a == 1) { // Special -- bare nucleon doesn't really "explode" 166 G4int knd = (z>0) ? 1 : 2; 167 particles.push_back(G4InuclElementaryParticle(knd)); // zero momentum 171 168 return; 172 } ;169 } 173 170 174 generateMomentumModules(etot, a, z); 171 // NOTE: If distribution fails, need to regenerate magnitudes and angles! 172 //*** generateMomentumModules(etot, a, z); 173 174 scm_momentums.reserve(a); 175 G4LorentzVector tot_mom; 176 175 177 G4bool bad = true; 176 178 G4int itry = 0; 177 178 179 while(bad && itry < itry_max) { 179 180 itry++; 180 s td::vector<G4LorentzVector> scm_momentums;181 G4LorentzVector tot_mom; 182 183 if (ia == 2) {184 // FIXME: This isn't actually a correct four-vector, wrong mass!181 scm_momentums.clear(); 182 183 generateMomentumModules(etot, a, z); 184 if (a == 2) { 185 // This is only a three-vector, not a four-vector 185 186 G4LorentzVector mom = generateWithRandomAngles(momModules[0]); 186 187 tot_mom += mom;188 189 187 scm_momentums.push_back(mom); 190 191 G4LorentzVector mom1 = -mom; 192 193 scm_momentums.push_back(mom1); 188 scm_momentums.push_back(-mom); // Only safe since three-vector! 194 189 bad = false; 195 190 } else { 196 for(G4int i = 0; i < ia - 2; i++) { 197 // FIXME: This isn't actually a correct four-vector, wrong mass! 191 tot_mom *= 0.; // Easy way to reset accumulator 192 193 for(G4int i = 0; i < a-2; i++) { // All but last two are thrown 194 // This is only a three-vector, not a four-vector 198 195 G4LorentzVector mom = generateWithRandomAngles(momModules[i]); 199 196 scm_momentums.push_back(mom); 200 197 tot_mom += mom; 201 202 scm_momentums.push_back(mom);203 198 }; 204 199 205 200 // handle last two 206 201 G4double tot_mod = tot_mom.rho(); 207 G4double ct = -0.5 * (tot_mod * tot_mod + momModules[ia - 2] * momModules[ia - 2] - 208 momModules[ia - 1] * momModules[ia - 1]) / tot_mod / momModules[ia - 2]; 209 210 if (verboseLevel > 2) { 211 G4cout << " ct last " << ct << G4endl; 212 } 202 G4double ct = -0.5*(tot_mod*tot_mod + momModules[a-2]*momModules[a-2] 203 - momModules[a-1]*momModules[a-1]) / tot_mod 204 / momModules[a-2]; 205 206 if (verboseLevel > 2) G4cout << " ct last " << ct << G4endl; 213 207 214 208 if(std::fabs(ct) < ang_cut) { 215 // FIXME: These are not actually four-vectors, just three-momenta 216 G4LorentzVector mom2 = generateWithFixedTheta(ct, momModules[ia - 2]); 217 // rotate to the normal system 209 // This is only a three-vector, not a four-vector 210 G4LorentzVector mom2 = generateWithFixedTheta(ct, momModules[a - 2]); 211 212 // rotate to the normal system 218 213 G4LorentzVector apr = tot_mom/tot_mod; 219 214 G4double a_tr = std::sqrt(apr.x()*apr.x() + apr.y()*apr.y()); … … 224 219 225 220 scm_momentums.push_back(mom); 226 // and the last one 227 G4LorentzVector mom1= - mom - tot_mom; 221 222 // and the last one (again, not actually a four-vector!) 223 G4LorentzVector mom1 = -mom - tot_mom; 224 228 225 scm_momentums.push_back(mom1); 229 226 bad = false; 230 } ;231 } ;232 if(!bad) {233 for(G4int i = 0; i < ia; i++) { 234 G4int knd = i < iz ? 1 : 2; 235 236 particles.push_back(G4InuclElementaryParticle(scm_momentums[i], knd, 8));237 };227 } // if (abs(ct) < ang_cut) 228 } // (a > 2) 229 } // while (bad && itry<itry_max) 230 231 if (!bad) { 232 for(G4int i = 0; i < a; i++) { 233 G4int knd = i < z ? 1 : 2; 234 particles.push_back(G4InuclElementaryParticle(scm_momentums[i], knd, 8)); 238 235 }; 239 }; 240 if (verboseLevel > 2) { 241 if(itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl; 236 }; 237 238 if (verboseLevel > 2) { 239 if (itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl; 242 240 } 243 241 … … 245 243 } 246 244 247 void G4BigBanger::generateMomentumModules(G4double etot, G4 double a, G4doublez) {245 void G4BigBanger::generateMomentumModules(G4double etot, G4int a, G4int z) { 248 246 if (verboseLevel > 3) { 249 247 G4cout << " >>> G4BigBanger::generateMomentumModules" << G4endl; … … 256 254 momModules.clear(); // Reset buffer for filling 257 255 258 G4int ia = G4int(a + 0.1);259 G4int iz = G4int(z + 0.1);260 261 256 G4double xtot = 0.0; 262 G4double promax = maxProbability(a); 263 264 G4int i; 265 for(i = 0; i < ia; i++) { 266 G4double x = generateX(ia, a, promax); 267 268 if (verboseLevel > 2) { 269 G4cout << " i " << i << " x " << x << G4endl; 257 258 if (a > 2) { // For "large" nuclei, energy is distributed 259 G4double promax = maxProbability(a); 260 261 for(G4int i = 0; i < a; i++) { 262 G4double x = generateX(a, promax); 263 264 if (verboseLevel > 2) { 265 G4cout << " i " << i << " x " << x << G4endl; 266 } 267 momModules.push_back(x); 268 xtot += x; 270 269 } 271 momModules.push_back(x); 272 xtot += x; 273 }; 274 for(i = 0; i < ia; i++) { 275 G4double m = i < iz ? mp : mn; 270 } else { // Two-body case is special, must be 50% 271 xtot = 1.; 272 momModules.push_back(0.5); 273 momModules.push_back(0.5); 274 } 275 276 for(G4int i = 0; i < a; i++) { 277 G4double m = i < z ? mp : mn; 276 278 277 279 momModules[i] = momModules[i] * etot / xtot; … … 286 288 } 287 289 288 G4double G4BigBanger::xProbability(G4double x, G4int ia) const { 289 290 291 if (verboseLevel > 3) { 292 G4cout << " >>> G4BigBanger::xProbability" << G4endl; 293 } 294 295 G4int ihalf = ia / 2; 290 G4double G4BigBanger::xProbability(G4double x, G4int a) const { 291 if (verboseLevel > 3) G4cout << " >>> G4BigBanger::xProbability" << G4endl; 292 296 293 G4double ekpr = 0.0; 297 294 … … 299 296 ekpr = x * x; 300 297 301 if (2 * ihalf == ia) { // even A302 ekpr *= std::sqrt(1.0 - x) * std::pow((1.0 - x), G4int(G4double(3 * ia - 6) / 2.0));298 if (a%2 == 0) { // even A 299 ekpr *= std::sqrt(1.0 - x) * std::pow((1.0 - x), (3*a-6)/2); 303 300 } 304 301 else { 305 ekpr *= std::pow((1.0 - x), G4int(G4double(3 * ia - 5) / 2.0));302 ekpr *= std::pow((1.0 - x), (3*a-5)/2); 306 303 }; 307 304 }; … … 310 307 } 311 308 312 G4double G4BigBanger::maxProbability(G4double a) const { 313 309 G4double G4BigBanger::maxProbability(G4int a) const { 314 310 if (verboseLevel > 3) { 315 311 G4cout << " >>> G4BigBanger::maxProbability" << G4endl; 316 312 } 317 313 318 return xProbability(1.0 / (a - 1.0) / 1.5, G4int(a + 0.1)); 319 } 320 321 G4double G4BigBanger::generateX(G4int ia, 322 G4double a, 323 G4double promax) const { 324 325 if (verboseLevel > 3) { 326 G4cout << " >>> G4BigBanger::generateX" << G4endl; 327 } 314 return xProbability(2./3./(a-1.0), a); 315 } 316 317 G4double G4BigBanger::generateX(G4int a, G4double promax) const { 318 if (verboseLevel > 3) G4cout << " >>> G4BigBanger::generateX" << G4endl; 328 319 329 320 const G4int itry_max = 1000; … … 335 326 x = inuclRndm(); 336 327 337 if(xProbability(x, ia) >= promax * inuclRndm()) return x;328 if(xProbability(x, a) >= promax * inuclRndm()) return x; 338 329 }; 339 330 if (verboseLevel > 2) {
Note: See TracChangeset
for help on using the changeset viewer.