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 #include "G4LEpp.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036 #include "G4ios.hh"
00037
00038
00039 #include "G4LEppData.hh"
00040
00041 G4LEpp::G4LEpp():G4HadronicInteraction("G4LEpp")
00042 {
00043
00044
00045
00046
00047 SetCoulombEffects(0);
00048
00049 SetMinEnergy(0.);
00050 SetMaxEnergy(5.*GeV);
00051 }
00052
00053 G4LEpp::~G4LEpp()
00054 {
00055
00056 }
00057
00058
00059 void
00060 G4LEpp::SetCoulombEffects(G4int State)
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 }
00079
00080
00081 G4HadFinalState*
00082 G4LEpp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
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
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
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
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
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
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
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
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
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
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
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
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
00322 theParticleChange.AddSecondary(targetParticle);
00323 return &theParticleChange;
00324 }
00325
00326