#include <G4LEpp.hh>
Inheritance diagram for G4LEpp:
Public Member Functions | |
G4LEpp () | |
~G4LEpp () | |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
void | SetCoulombEffects (G4int State) |
Definition at line 50 of file G4LEpp.hh.
G4LEpp::G4LEpp | ( | ) |
Definition at line 41 of file G4LEpp.cc.
References SetCoulombEffects(), G4HadronicInteraction::SetMaxEnergy(), and G4HadronicInteraction::SetMinEnergy().
00041 :G4HadronicInteraction("G4LEpp") 00042 { 00043 // theParticleChange.SetNumberOfSecondaries(1); 00044 // SetMinEnergy(10.*MeV); 00045 // SetMaxEnergy(1200.*MeV); 00046 00047 SetCoulombEffects(0); 00048 00049 SetMinEnergy(0.); 00050 SetMaxEnergy(5.*GeV); 00051 }
G4LEpp::~G4LEpp | ( | ) |
G4HadFinalState * G4LEpp::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 82 of file G4LEpp.cc.
References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4DynamicParticle::GetDefinition(), G4HadProjectile::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4DynamicParticle::GetMomentum(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4DynamicParticle::GetTotalEnergy(), G4HadProjectile::GetTotalEnergy(), G4DynamicParticle::GetTotalMomentum(), G4HadProjectile::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4INCL::Math::pi, G4Proton::Proton(), G4Nucleus::ReturnTargetParticle(), G4DynamicParticle::SetDefinition(), G4HadFinalState::SetEnergyChange(), G4DynamicParticle::SetMomentum(), G4HadFinalState::SetMomentumChange(), G4HadronicInteraction::theParticleChange, and G4HadronicInteraction::verboseLevel.
00083 { 00084 theParticleChange.Clear(); 00085 const G4HadProjectile* aParticle = &aTrack; 00086 00087 G4double P = aParticle->GetTotalMomentum(); 00088 G4double Px = aParticle->Get4Momentum().x(); 00089 G4double Py = aParticle->Get4Momentum().y(); 00090 G4double Pz = aParticle->Get4Momentum().z(); 00091 G4double ek = aParticle->GetKineticEnergy(); 00092 G4ThreeVector theInitial = aParticle->Get4Momentum().vect(); 00093 00094 if (verboseLevel > 1) { 00095 G4double E = aParticle->GetTotalEnergy(); 00096 G4double E0 = aParticle->GetDefinition()->GetPDGMass(); 00097 G4double Q = aParticle->GetDefinition()->GetPDGCharge(); 00098 G4int A = targetNucleus.GetA_asInt(); 00099 G4int Z = targetNucleus.GetZ_asInt(); 00100 G4cout << "G4LEpp:ApplyYourself: incident particle: " 00101 << aParticle->GetDefinition()->GetParticleName() << G4endl; 00102 G4cout << "P = " << P/GeV << " GeV/c" 00103 << ", Px = " << Px/GeV << " GeV/c" 00104 << ", Py = " << Py/GeV << " GeV/c" 00105 << ", Pz = " << Pz/GeV << " GeV/c" << G4endl; 00106 G4cout << "E = " << E/GeV << " GeV" 00107 << ", kinetic energy = " << ek/GeV << " GeV" 00108 << ", mass = " << E0/GeV << " GeV" 00109 << ", charge = " << Q << G4endl; 00110 G4cout << "G4LEpp:ApplyYourself: material:" << G4endl; 00111 G4cout << "A = " << A 00112 << ", Z = " << Z 00113 << ", atomic mass " 00114 << G4Proton::Proton()->GetPDGMass()/GeV << "GeV" 00115 << G4endl; 00116 // 00117 // GHEISHA ADD operation to get total energy, mass, charge 00118 // 00119 E += proton_mass_c2; 00120 G4double E02 = E*E - P*P; 00121 E0 = std::sqrt(std::fabs(E02)); 00122 if (E02 < 0)E0 *= -1; 00123 Q += Z; 00124 G4cout << "G4LEpp:ApplyYourself: total:" << G4endl; 00125 G4cout << "E = " << E/GeV << " GeV" 00126 << ", mass = " << E0/GeV << " GeV" 00127 << ", charge = " << Q << G4endl; 00128 } 00129 00130 // Find energy bin 00131 00132 G4int je1 = 0; 00133 G4int je2 = NENERGY - 1; 00134 ek = ek/GeV; 00135 do { 00136 G4int midBin = (je1 + je2)/2; 00137 if (ek < elab[midBin]) 00138 je2 = midBin; 00139 else 00140 je1 = midBin; 00141 } while (je2 - je1 > 1); 00142 G4double delab = elab[je2] - elab[je1]; 00143 00144 // Sample the angle 00145 00146 G4float sample = G4UniformRand(); 00147 G4int ke1 = 0; 00148 G4int ke2 = NANGLE - 1; 00149 G4double dsig = sig[je2][0] - sig[je1][0]; 00150 G4double rc = dsig/delab; 00151 G4double b = sig[je1][0] - rc*elab[je1]; 00152 G4double sigint1 = rc*ek + b; 00153 G4double sigint2 = 0.; 00154 00155 if (verboseLevel > 1) G4cout << "sample=" << sample << G4endl 00156 << ke1 << " " << ke2 << " " 00157 << sigint1 << " " << sigint2 << G4endl; 00158 00159 do { 00160 G4int midBin = (ke1 + ke2)/2; 00161 dsig = sig[je2][midBin] - sig[je1][midBin]; 00162 rc = dsig/delab; 00163 b = sig[je1][midBin] - rc*elab[je1]; 00164 G4double sigint = rc*ek + b; 00165 if (sample < sigint) { 00166 ke2 = midBin; 00167 sigint2 = sigint; 00168 } 00169 else { 00170 ke1 = midBin; 00171 sigint1 = sigint; 00172 } 00173 if (verboseLevel > 1)G4cout << ke1 << " " << ke2 << " " 00174 << sigint1 << " " << sigint2 << G4endl; 00175 } while (ke2 - ke1 > 1); 00176 00177 dsig = sigint2 - sigint1; 00178 rc = 1./dsig; 00179 b = ke1 - rc*sigint1; 00180 G4double kint = rc*sample + b; 00181 G4double theta = (0.5 + kint)*pi/180.; 00182 if (theta < 0.) { theta = 0.; } 00183 00184 if (verboseLevel > 1) { 00185 G4cout << " energy bin " << je1 << " energy=" << elab[je1] << G4endl; 00186 G4cout << " angle bin " << kint << " angle=" << theta/degree << G4endl; 00187 } 00188 00189 // Get the target particle 00190 G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle(); 00191 00192 G4double E1 = aParticle->GetTotalEnergy(); 00193 G4double M1 = aParticle->GetDefinition()->GetPDGMass(); 00194 G4double E2 = targetParticle->GetTotalEnergy(); 00195 G4double M2 = targetParticle->GetDefinition()->GetPDGMass(); 00196 G4double totalEnergy = E1 + E2; 00197 G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P); 00198 00199 // Transform into centre of mass system 00200 00201 G4double px = (M2/pseudoMass)*Px; 00202 G4double py = (M2/pseudoMass)*Py; 00203 G4double pz = (M2/pseudoMass)*Pz; 00204 G4double p = std::sqrt(px*px + py*py + pz*pz); 00205 00206 if (verboseLevel > 1) { 00207 G4cout << " E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl; 00208 G4cout << " E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl; 00209 G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " " 00210 << pz/GeV << " " << p/GeV << G4endl; 00211 } 00212 00213 // First scatter w.r.t. Z axis 00214 G4double phi = G4UniformRand()*twopi; 00215 G4double pxnew = p*std::sin(theta)*std::cos(phi); 00216 G4double pynew = p*std::sin(theta)*std::sin(phi); 00217 G4double pznew = p*std::cos(theta); 00218 00219 // Rotate according to the direction of the incident particle 00220 if (px*px + py*py > 0) { 00221 G4double cost, sint, ph, cosp, sinp; 00222 cost = pz/p; 00223 sint = (std::sqrt(std::fabs((1-cost)*(1+cost))) + std::sqrt(px*px+py*py)/p)/2; 00224 py < 0 ? ph = 3*halfpi : ph = halfpi; 00225 if (std::fabs(px) > 0.000001*GeV) ph = std::atan2(py,px); 00226 cosp = std::cos(ph); 00227 sinp = std::sin(ph); 00228 px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew); 00229 py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew); 00230 pz = (-sint*pxnew + cost*pznew); 00231 } 00232 else { 00233 px = pxnew; 00234 py = pynew; 00235 pz = pznew; 00236 } 00237 00238 if (verboseLevel > 1) { 00239 G4cout << " AFTER SCATTER..." << G4endl; 00240 G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " " 00241 << pz/GeV << " " << p/GeV << G4endl; 00242 } 00243 00244 // Transform to lab system 00245 00246 G4double E1pM2 = E1 + M2; 00247 G4double betaCM = P/E1pM2; 00248 G4double betaCMx = Px/E1pM2; 00249 G4double betaCMy = Py/E1pM2; 00250 G4double betaCMz = Pz/E1pM2; 00251 G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P); 00252 00253 if (verboseLevel > 1) { 00254 G4cout << " betaCM " << betaCMx << " " << betaCMy << " " 00255 << betaCMz << " " << betaCM << G4endl; 00256 G4cout << " gammaCM " << gammaCM << G4endl; 00257 } 00258 00259 // Now following GLOREN... 00260 00261 G4double BETA[5], PA[5], PB[5]; 00262 BETA[1] = -betaCMx; 00263 BETA[2] = -betaCMy; 00264 BETA[3] = -betaCMz; 00265 BETA[4] = gammaCM; 00266 00267 //The incident particle... 00268 00269 PA[1] = px; 00270 PA[2] = py; 00271 PA[3] = pz; 00272 PA[4] = std::sqrt(M1*M1 + p*p); 00273 00274 G4double BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3]; 00275 G4double BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4]; 00276 00277 PB[1] = PA[1] + BPGAM * BETA[1]; 00278 PB[2] = PA[2] + BPGAM * BETA[2]; 00279 PB[3] = PA[3] + BPGAM * BETA[3]; 00280 PB[4] = (PA[4] - BETPA) * BETA[4]; 00281 00282 G4DynamicParticle* newP = new G4DynamicParticle; 00283 newP->SetDefinition(const_cast<G4ParticleDefinition *>(aParticle->GetDefinition()) ); 00284 newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3])); 00285 00286 //The target particle... 00287 00288 PA[1] = -px; 00289 PA[2] = -py; 00290 PA[3] = -pz; 00291 PA[4] = std::sqrt(M2*M2 + p*p); 00292 00293 BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3]; 00294 BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4]; 00295 00296 PB[1] = PA[1] + BPGAM * BETA[1]; 00297 PB[2] = PA[2] + BPGAM * BETA[2]; 00298 PB[3] = PA[3] + BPGAM * BETA[3]; 00299 PB[4] = (PA[4] - BETPA) * BETA[4]; 00300 00301 targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3])); 00302 00303 if (verboseLevel > 1) { 00304 G4cout << " particle 1 momentum in LAB " 00305 << newP->GetMomentum()/GeV 00306 << " " << newP->GetTotalMomentum()/GeV << G4endl; 00307 G4cout << " particle 2 momentum in LAB " 00308 << targetParticle->GetMomentum()/GeV 00309 << " " << targetParticle->GetTotalMomentum()/GeV << G4endl; 00310 G4cout << " TOTAL momentum in LAB " 00311 << (newP->GetMomentum()+targetParticle->GetMomentum())/GeV 00312 << " " 00313 << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV 00314 << G4endl; 00315 } 00316 00317 theParticleChange.SetMomentumChange( newP->GetMomentumDirection()); 00318 theParticleChange.SetEnergyChange(newP->GetKineticEnergy()); 00319 delete newP; 00320 00321 // Recoil particle 00322 theParticleChange.AddSecondary(targetParticle); 00323 return &theParticleChange; 00324 }
void G4LEpp::SetCoulombEffects | ( | G4int | State | ) |
Definition at line 60 of file G4LEpp.cc.
References G4HadronicInteraction::SetMaxEnergy().
Referenced by G4LEpp().
00061 { 00062 if (State) { 00063 for(G4int i=0; i<NANGLE; i++) 00064 { 00065 sig[i] = SigCoul[i]; 00066 } 00067 elab = ElabCoul; 00068 SetMaxEnergy(1.2*GeV); 00069 } 00070 else { 00071 for(G4int i=0; i<NANGLE; i++) 00072 { 00073 sig[i] = Sig[i]; 00074 } 00075 elab = Elab; 00076 SetMaxEnergy(5.*GeV); 00077 } 00078 }