00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "G4LEnp.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "Randomize.hh"
00038 #include "G4ios.hh"
00039
00040
00041 #include "G4LEnpData.hh"
00042 #include "Randomize.hh"
00043
00044
00045 G4LEnp::G4LEnp():G4HadronicInteraction("G4LEnp")
00046 {
00047
00048
00049
00050
00051 SetMinEnergy(0.);
00052 SetMaxEnergy(5.*GeV);
00053 }
00054
00055 G4LEnp::~G4LEnp()
00056 {
00057 theParticleChange.Clear();
00058 }
00059
00060 G4HadFinalState*
00061 G4LEnp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
00062 {
00063 theParticleChange.Clear();
00064 const G4HadProjectile* aParticle = &aTrack;
00065
00066 G4double P = aParticle->GetTotalMomentum();
00067 G4double Px = aParticle->Get4Momentum().x();
00068 G4double Py = aParticle->Get4Momentum().y();
00069 G4double Pz = aParticle->Get4Momentum().z();
00070 G4double ek = aParticle->GetKineticEnergy();
00071 G4ThreeVector theInitial = aParticle->Get4Momentum().vect();
00072
00073 if (verboseLevel > 1) {
00074 G4double E = aParticle->GetTotalEnergy();
00075 G4double E0 = aParticle->GetDefinition()->GetPDGMass();
00076 G4double Q = aParticle->GetDefinition()->GetPDGCharge();
00077 G4int A = targetNucleus.GetA_asInt();
00078 G4int Z = targetNucleus.GetZ_asInt();
00079 G4cout << "G4LEnp:ApplyYourself: incident particle: "
00080 << aParticle->GetDefinition()->GetParticleName() << G4endl;
00081 G4cout << "P = " << P/GeV << " GeV/c"
00082 << ", Px = " << Px/GeV << " GeV/c"
00083 << ", Py = " << Py/GeV << " GeV/c"
00084 << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
00085 G4cout << "E = " << E/GeV << " GeV"
00086 << ", kinetic energy = " << ek/GeV << " GeV"
00087 << ", mass = " << E0/GeV << " GeV"
00088 << ", charge = " << Q << G4endl;
00089 G4cout << "G4LEnp:ApplyYourself: material:" << G4endl;
00090 G4cout << "A = " << A
00091 << ", Z = " << Z
00092 << ", atomic mass "
00093 << G4Proton::Proton()->GetPDGMass()/GeV << "GeV"
00094 << G4endl;
00095
00096
00097
00098 E += proton_mass_c2;
00099 G4double E02 = E*E - P*P;
00100 E0 = std::sqrt(std::abs(E02));
00101 if (E02 < 0)E0 *= -1;
00102 Q += Z;
00103 G4cout << "G4LEnp:ApplyYourself: total:" << G4endl;
00104 G4cout << "E = " << E/GeV << " GeV"
00105 << ", mass = " << E0/GeV << " GeV"
00106 << ", charge = " << Q << G4endl;
00107 }
00108
00109
00110
00111 G4int je1 = 0;
00112 G4int je2 = NENERGY - 1;
00113 ek = ek/GeV;
00114 do {
00115 G4int midBin = (je1 + je2)/2;
00116 if (ek < elab[midBin])
00117 je2 = midBin;
00118 else
00119 je1 = midBin;
00120 } while (je2 - je1 > 1);
00121 G4double delab = elab[je2] - elab[je1];
00122
00123
00124
00125 G4float sample = G4UniformRand();
00126 G4int ke1 = 0;
00127 G4int ke2 = NANGLE - 1;
00128 G4double dsig = sig[je2][0] - sig[je1][0];
00129 G4double rc = dsig/delab;
00130 G4double b = sig[je1][0] - rc*elab[je1];
00131 G4double sigint1 = rc*ek + b;
00132 G4double sigint2 = 0.;
00133
00134 if (verboseLevel > 1) {
00135 G4cout << "sample=" << sample << G4endl
00136 << ke1 << " " << ke2 << " "
00137 << sigint1 << " " << sigint2 << G4endl;
00138 }
00139 do {
00140 G4int midBin = (ke1 + ke2)/2;
00141 dsig = sig[je2][midBin] - sig[je1][midBin];
00142 rc = dsig/delab;
00143 b = sig[je1][midBin] - rc*elab[je1];
00144 G4double sigint = rc*ek + b;
00145 if (sample < sigint) {
00146 ke2 = midBin;
00147 sigint2 = sigint;
00148 }
00149 else {
00150 ke1 = midBin;
00151 sigint1 = sigint;
00152 }
00153 if (verboseLevel > 1) {
00154 G4cout << ke1 << " " << ke2 << " "
00155 << sigint1 << " " << sigint2 << G4endl;
00156 }
00157 } while (ke2 - ke1 > 1);
00158
00159 dsig = sigint2 - sigint1;
00160 rc = 1./dsig;
00161 b = ke1 - rc*sigint1;
00162 G4double kint = rc*sample + b;
00163 G4double theta = (0.5 + kint)*pi/180.;
00164
00165 if (verboseLevel > 1) {
00166 G4cout << " energy bin " << je1 << " energy=" << elab[je1] << G4endl;
00167 G4cout << " angle bin " << kint << " angle=" << theta/degree << G4endl;
00168 }
00169
00170
00171
00172 G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
00173
00174 G4double E1 = aParticle->GetTotalEnergy();
00175 G4double M1 = aParticle->GetDefinition()->GetPDGMass();
00176 G4double E2 = targetParticle->GetTotalEnergy();
00177 G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
00178 G4double totalEnergy = E1 + E2;
00179 G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
00180
00181
00182
00183 G4double px = (M2/pseudoMass)*Px;
00184 G4double py = (M2/pseudoMass)*Py;
00185 G4double pz = (M2/pseudoMass)*Pz;
00186 G4double p = std::sqrt(px*px + py*py + pz*pz);
00187
00188 if (verboseLevel > 1) {
00189 G4cout << " E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
00190 G4cout << " E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
00191 G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
00192 << pz/GeV << " " << p/GeV << G4endl;
00193 }
00194
00195
00196 G4double phi = G4UniformRand()*twopi;
00197 G4double pxnew = p*std::sin(theta)*std::cos(phi);
00198 G4double pynew = p*std::sin(theta)*std::sin(phi);
00199 G4double pznew = p*std::cos(theta);
00200
00201
00202 if (px*px + py*py > 0) {
00203 G4double cost, sint, ph, cosp, sinp;
00204 cost = pz/p;
00205 sint = (std::sqrt(std::fabs((1-cost)*(1+cost))) + std::sqrt(px*px+py*py)/p)/2;
00206 py < 0 ? ph = 3*halfpi : ph = halfpi;
00207 if (std::abs(px) > 0.000001*GeV) ph = std::atan2(py,px);
00208 cosp = std::cos(ph);
00209 sinp = std::sin(ph);
00210 px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
00211 py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
00212 pz = (-sint*pxnew + cost*pznew);
00213 }
00214 else {
00215 px = pxnew;
00216 py = pynew;
00217 pz = pznew;
00218 }
00219
00220 if (verboseLevel > 1) {
00221 G4cout << " AFTER SCATTER..." << G4endl;
00222 G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
00223 << pz/GeV << " " << p/GeV << G4endl;
00224 }
00225
00226
00227
00228 G4double E1pM2 = E1 + M2;
00229 G4double betaCM = P/E1pM2;
00230 G4double betaCMx = Px/E1pM2;
00231 G4double betaCMy = Py/E1pM2;
00232 G4double betaCMz = Pz/E1pM2;
00233 G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
00234
00235 if (verboseLevel > 1) {
00236 G4cout << " betaCM " << betaCMx << " " << betaCMy << " "
00237 << betaCMz << " " << betaCM << G4endl;
00238 G4cout << " gammaCM " << gammaCM << G4endl;
00239 }
00240
00241
00242
00243 G4double BETA[5], PA[5], PB[5];
00244 BETA[1] = -betaCMx;
00245 BETA[2] = -betaCMy;
00246 BETA[3] = -betaCMz;
00247 BETA[4] = gammaCM;
00248
00249
00250
00251 PA[1] = px;
00252 PA[2] = py;
00253 PA[3] = pz;
00254 PA[4] = std::sqrt(M1*M1 + p*p);
00255
00256 G4double BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
00257 G4double BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
00258
00259 PB[1] = PA[1] + BPGAM * BETA[1];
00260 PB[2] = PA[2] + BPGAM * BETA[2];
00261 PB[3] = PA[3] + BPGAM * BETA[3];
00262 PB[4] = (PA[4] - BETPA) * BETA[4];
00263
00264 G4DynamicParticle* newP = new G4DynamicParticle;
00265 newP->SetDefinition(const_cast<G4ParticleDefinition *>(aParticle->GetDefinition()));
00266 newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
00267
00268
00269
00270 PA[1] = -px;
00271 PA[2] = -py;
00272 PA[3] = -pz;
00273 PA[4] = std::sqrt(M2*M2 + p*p);
00274
00275 BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
00276 BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
00277
00278 PB[1] = PA[1] + BPGAM * BETA[1];
00279 PB[2] = PA[2] + BPGAM * BETA[2];
00280 PB[3] = PA[3] + BPGAM * BETA[3];
00281 PB[4] = (PA[4] - BETPA) * BETA[4];
00282
00283 targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
00284
00285 if (verboseLevel > 1) {
00286 G4cout << " particle 1 momentum in LAB "
00287 << newP->GetMomentum()*(1./GeV)
00288 << " " << newP->GetTotalMomentum()/GeV << G4endl;
00289 G4cout << " particle 2 momentum in LAB "
00290 << targetParticle->GetMomentum()*(1./GeV)
00291 << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
00292 G4cout << " TOTAL momentum in LAB "
00293 << (newP->GetMomentum()+targetParticle->GetMomentum())*(1./GeV)
00294 << " "
00295 << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
00296 << G4endl;
00297 }
00298
00299 theParticleChange.SetMomentumChange(newP->GetMomentumDirection());
00300 theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
00301 delete newP;
00302 theParticleChange.AddSecondary(targetParticle);
00303
00304 return &theParticleChange;
00305 }
00306
00307