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 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00037 #include "G4INCLElasticChannel.hh"
00038 #include "G4INCLRandom.hh"
00039 #include "G4INCLKinematicsUtils.hh"
00040 #include "G4INCLParticleTable.hh"
00041 #include "G4INCLCrossSections.hh"
00042 #include "G4INCLGlobals.hh"
00043
00044 namespace G4INCL {
00045
00046 ElasticChannel::ElasticChannel(Nucleus *n, Particle *p1, Particle *p2)
00047 :theNucleus(n), particle1(p1), particle2(p2)
00048 {
00049 }
00050
00051 ElasticChannel::~ElasticChannel()
00052 {
00053 }
00054
00055 FinalState* ElasticChannel::getFinalState()
00056 {
00057 ParticleType p1TypeOld = particle1->getType();
00058 ParticleType p2TypeOld = particle2->getType();
00059
00060
00061
00062
00063 const G4double s = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
00064 const G4double pl = KinematicsUtils::momentumInLab(s, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
00065
00066 const G4int isospin = ParticleTable::getIsospin(particle1->getType()) +
00067 ParticleTable::getIsospin(particle2->getType());
00068
00069
00070 G4double psq = particle1->getMomentum().mag2();
00071 G4double pnorm = std::sqrt(psq);
00072 G4double b = CrossSections::calculateNNDiffCrossSection(pl, isospin);
00073 G4double btmax = 4.0 * psq * b;
00074 G4double z = std::exp(-btmax);
00075 G4double ranres = Random::shoot();
00076 G4double y = 1.0 - ranres * (1.0 - z);
00077 G4double T = std::log(y)/b;
00078 G4int iexpi = 0;
00079 G4double apt = 1.0;
00080
00081
00082 if((particle1->getType() == Proton && particle2->getType() == Neutron) ||
00083 (particle1->getType() == Neutron && particle2->getType() == Proton)) {
00084 if(pl > 800.0) {
00085 const G4double x = 0.001 * pl;
00086 apt = (800.0/pl)*(800.0/pl);
00087 G4double cpt = std::max(6.23 * std::exp(-1.79*x), 0.3);
00088 G4double alphac = 100.0 * 1.0e-6;
00089 G4double aaa = (1 + apt) * (1 - std::exp(-btmax))/b;
00090 G4double argu = psq * alphac;
00091
00092 if(argu >= 8) {
00093 argu = 0.0;
00094 } else {
00095 argu = std::exp(-4.0 * argu);
00096 }
00097
00098 G4double aac = cpt * (1.0 - argu)/alphac;
00099 G4double fracpn = aaa/(aac + aaa);
00100 if(Random::shoot() > fracpn) {
00101 z = std::exp(-4.0 * psq *alphac);
00102 iexpi = 1;
00103 y = 1.0 - ranres*(1.0 - z);
00104 T = std::log(y)/alphac;
00105 }
00106 }
00107 }
00108
00109 G4double ctet = 1.0 + 0.5*T/psq;
00110 if(std::abs(ctet) > 1.0) ctet = Math::sign(ctet);
00111 G4double stet = std::sqrt(1.0 - ctet*ctet);
00112 G4double rndm = Random::shoot();
00113
00114 G4double fi = Math::twoPi * rndm;
00115 G4double cfi = std::cos(fi);
00116 G4double sfi = std::sin(fi);
00117
00118 G4double xx = particle1->getMomentum().perp2();
00119 G4double zz = std::pow(particle1->getMomentum().getZ(), 2);
00120
00121 if(xx >= (zz * 1.0e-8)) {
00122 ThreeVector p = particle1->getMomentum();
00123 G4double yn = std::sqrt(xx);
00124 G4double zn = yn * pnorm;
00125 G4double ex[3], ey[3], ez[3];
00126 ez[0] = p.getX() / pnorm;
00127 ez[1] = p.getY() / pnorm;
00128 ez[2] = p.getZ() / pnorm;
00129
00130
00131 ex[0] = p.getY() / yn;
00132 ex[1] = -p.getX() / yn;
00133 ex[2] = 0.0;
00134
00135 ey[0] = p.getX() * p.getZ() / zn;
00136 ey[1] = p.getY() * p.getZ() / zn;
00137 ey[2] = -xx/zn;
00138
00139 G4double pX = (ex[0]*cfi*stet + ey[0]*sfi*stet + ez[0]*ctet) * pnorm;
00140 G4double pY = (ex[1]*cfi*stet + ey[1]*sfi*stet + ez[1]*ctet) * pnorm;
00141 G4double pZ = (ex[2]*cfi*stet + ey[2]*sfi*stet + ez[2]*ctet) * pnorm;
00142
00143 ThreeVector p1momentum = ThreeVector(pX, pY, pZ);
00144 particle1->setMomentum(p1momentum);
00145 particle2->setMomentum(-p1momentum);
00146 } else {
00147 G4double momZ = particle1->getMomentum().getZ();
00148 G4double pX = momZ * cfi * stet;
00149 G4double pY = momZ * sfi * stet;
00150 G4double pZ = momZ * ctet;
00151
00152 ThreeVector p1momentum(pX, pY, pZ);
00153 particle1->setMomentum(p1momentum);
00154 particle2->setMomentum(-p1momentum);
00155 }
00156
00157
00158
00159 if((particle1->getType() == Proton && particle2->getType() == Neutron) ||
00160 (particle1->getType() == Neutron && particle2->getType() == Proton)) {
00161 rndm = Random::shoot();
00162 apt = 1.0;
00163 if(pl > 800.0) {
00164 apt = std::pow(800.0/pl, 2);
00165 }
00166 if(iexpi == 1 || rndm > 1.0/(1.0 + apt)) {
00167 particle1->setType(p2TypeOld);
00168 particle2->setType(p1TypeOld);
00169 }
00170 }
00171
00172
00173
00174
00175 FinalState *fs = new FinalState();
00176 fs->addModifiedParticle(particle1);
00177 fs->addModifiedParticle(particle2);
00178
00179 return fs;
00180
00181 }
00182
00183 }