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 "G4INCLDeltaProductionChannel.hh"
00038 #include "G4INCLKinematicsUtils.hh"
00039 #include "G4INCLBinaryCollisionAvatar.hh"
00040 #include "G4INCLRandom.hh"
00041 #include "G4INCLGlobals.hh"
00042 #include "G4INCLLogger.hh"
00043
00044 namespace G4INCL {
00045
00046 DeltaProductionChannel::DeltaProductionChannel(Particle *p1,
00047 Particle *p2,
00048 Nucleus *n)
00049 :theNucleus(n), particle1(p1), particle2(p2)
00050 {}
00051
00052 DeltaProductionChannel::~DeltaProductionChannel() {}
00053
00054 G4double DeltaProductionChannel::sampleDeltaMass(G4double ecm) {
00055 const G4double ramass = 0.0;
00056 const G4int maxTries = 100000;
00057 G4int nTries = 0;
00058 deltaProd101: G4double rndm = Random::shoot();
00059 nTries++;
00060 G4double y = std::tan(Math::pi*(rndm-0.5));
00061 G4double x = 1232.+0.5*130.*y+ramass;
00062 if (x < ParticleTable::effectiveDeltaDecayThreshold && (nTries < maxTries))
00063 goto deltaProd101;
00064 if (ecm < x + ParticleTable::effectiveNucleonMass + 1.0 && (nTries < maxTries)) goto deltaProd101;
00065
00066
00067
00068 y=ecm*ecm;
00069 G4double q2=(y-1.157776E6)*(y-6.4E5)/y/4.0;
00070 G4double q3=std::pow(std::sqrt(q2), 3.);
00071 G4double f3max=q3/(q3+5.832E6);
00072 y=x*x;
00073 q2=(y-1.157776E6)*(y-6.4E5)/y/4.0;
00074 q3=std::pow(std::sqrt(q2), 3.);
00075 G4double f3=q3/(q3+5.832E6);
00076 rndm = Random::shoot();
00077 if (rndm > f3/f3max && (nTries < maxTries)) goto deltaProd101;
00078 if(nTries >= maxTries) {
00079 WARN("DeltaProductionChannel::sampleDeltaMass loop was stopped because maximum number of tries was reached. Delta mass " << x << " MeV with CM energy " << ecm << " MeV may be unphysical." << std::endl);
00080 }
00081 return x;
00082 }
00083
00084 FinalState* DeltaProductionChannel::getFinalState() {
00093
00094
00095
00096 G4double ecm = KinematicsUtils::totalEnergyInCM(particle1, particle2);
00097
00098 const G4int isospin = ParticleTable::getIsospin(particle1->getType()) +
00099 ParticleTable::getIsospin(particle2->getType());
00100
00101
00102 G4double pin = particle1->getMomentum().mag();
00103 G4double rndm = 0.0, b = 0.0;
00104
00105 G4double xmdel = sampleDeltaMass(ecm);
00106
00107 G4double pnorm = KinematicsUtils::momentumInCM(ecm, ParticleTable::effectiveNucleonMass, xmdel);
00108 if (pnorm <= 0.0) pnorm=0.000001;
00109 G4int index=0;
00110 G4int index2=0;
00111 rndm = Random::shoot();
00112 if (rndm < 0.5) index=1;
00113 if (isospin == 0) {
00114 rndm = Random::shoot();
00115 if (rndm < 0.5) index2=1;
00116 }
00117
00118
00119
00120 G4double x = 0.001 * KinematicsUtils::momentumInLab(ecm*ecm, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass);
00121 if(x < 1.4) {
00122 b=(5.287/(1.+std::exp((1.3-x)/0.05)))*1.e-6;
00123 } else {
00124 b=(4.65+0.706*(x-1.4))*1.e-6;
00125 }
00126 G4double xkh = 2.*b*pin*pnorm;
00127 rndm = Random::shoot();
00128 G4double ctet=1.0+std::log(1.-rndm*(1.-std::exp(-2.*xkh)))/xkh;
00129 if(std::abs(ctet) > 1.0) ctet = Math::sign(ctet);
00130 G4double stet = std::sqrt(1.-ctet*ctet);
00131
00132 rndm = Random::shoot();
00133 G4double fi = Math::twoPi*rndm;
00134 G4double cfi = std::cos(fi);
00135 G4double sfi = std::sin(fi);
00136
00137
00138 G4double xx = particle1->getMomentum().perp2();
00139 G4double zz = std::pow(particle1->getMomentum().getZ(), 2);
00140 G4double xp1, xp2, xp3;
00141 if (xx >= zz*1.e-8) {
00142 G4double yn = std::sqrt(xx);
00143 G4double zn = yn*pin;
00144 G4double ex[3], ey[3], ez[3];
00145 G4double p1 = particle1->getMomentum().getX();
00146 G4double p2 = particle1->getMomentum().getY();
00147 G4double p3 = particle1->getMomentum().getZ();
00148 ez[0] = p1/pin;
00149 ez[1] = p2/pin;
00150 ez[2] = p3/pin;
00151 ex[0] = p2/yn;
00152 ex[1] = -p1/yn;
00153 ex[2] = 0.0;
00154 ey[0] = p1*p3/zn;
00155 ey[1] = p2*p3/zn;
00156 ey[2] = -xx/zn;
00157 xp1 = (ex[0]*cfi*stet+ey[0]*sfi*stet+ez[0]*ctet)*pnorm;
00158 xp2 = (ex[1]*cfi*stet+ey[1]*sfi*stet+ez[1]*ctet)*pnorm;
00159 xp3 = (ex[2]*cfi*stet+ey[2]*sfi*stet+ez[2]*ctet)*pnorm;
00160 }else {
00161 xp1=pnorm*stet*cfi;
00162 xp2=pnorm*stet*sfi;
00163 xp3=pnorm*ctet;
00164 }
00165
00166 G4double e3 = std::sqrt(xp1*xp1+xp2*xp2+xp3*xp3
00167 +ParticleTable::effectiveNucleonMass2);
00168
00169
00170
00171 G4int m1 = 0;
00172 G4int m2 = 0;
00173 if (index != 1) {
00174 ThreeVector mom(xp1, xp2, xp3);
00175 particle1->setMomentum(mom);
00176
00177 m1=1;
00178 } else {
00179 ThreeVector mom(-xp1, -xp2, -xp3);
00180 particle1->setMomentum(mom);
00181
00182 m1=1;
00183 }
00184
00185 particle1->setEnergy(ecm - e3);
00186 particle2->setEnergy(e3);
00187 particle2->setMomentum(-particle1->getMomentum());
00188
00189
00190
00191
00192
00193 G4int is1 = ParticleTable::getIsospin(particle1->getType());
00194 G4int is2 = ParticleTable::getIsospin(particle2->getType());
00195 if (isospin == 0) {
00196 if(index2 == 1) {
00197 G4int isi=is1;
00198 is1=is2;
00199 is2=isi;
00200 }
00201 particle1->setHelicity(0.0);
00202 } else {
00203 rndm = Random::shoot();
00204 if (rndm >= 0.25) {
00205 is1=3*is1*m1-(1-m1)*is1;
00206 is2=3*is2*m2-(1-m2)*is2;
00207 }
00208 particle1->setHelicity(ctet*ctet);
00209 }
00210
00211 if(is1 == ParticleTable::getIsospin(Proton) && m1 == 0) {
00212 particle1->setType(Proton);
00213 } else if(is1 == ParticleTable::getIsospin(Neutron) && m1 == 0) {
00214 particle1->setType(Neutron);
00215 } else if(is1 == ParticleTable::getIsospin(DeltaMinus) && m1 == 1) {
00216 particle1->setType(DeltaMinus);
00217 } else if(is1 == ParticleTable::getIsospin(DeltaZero) && m1 == 1) {
00218 particle1->setType(DeltaZero);
00219 } else if(is1 == ParticleTable::getIsospin(DeltaPlus) && m1 == 1) {
00220 particle1->setType(DeltaPlus);
00221 } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus) && m1 == 1) {
00222 particle1->setType(DeltaPlusPlus);
00223 }
00224
00225 if(is2 == ParticleTable::getIsospin(Proton) && m2 == 0) {
00226 particle2->setType(Proton);
00227 } else if(is2 == ParticleTable::getIsospin(Neutron) && m2 == 0) {
00228 particle2->setType(Neutron);
00229 } else if(is2 == ParticleTable::getIsospin(DeltaMinus) && m2 == 1) {
00230 particle2->setType(DeltaMinus);
00231 } else if(is2 == ParticleTable::getIsospin(DeltaZero) && m2 == 1) {
00232 particle2->setType(DeltaZero);
00233 } else if(is2 == ParticleTable::getIsospin(DeltaPlus) && m2 == 1) {
00234 particle2->setType(DeltaPlus);
00235 } else if(is2 == ParticleTable::getIsospin(DeltaPlusPlus) && m2 == 1) {
00236 particle2->setType(DeltaPlusPlus);
00237 }
00238
00239 if(particle1->isDelta()) particle1->setMass(xmdel);
00240 if(particle2->isDelta()) particle2->setMass(xmdel);
00241
00242 FinalState *fs = new FinalState;
00243 fs->addModifiedParticle(particle1);
00244 fs->addModifiedParticle(particle2);
00245 return fs;
00246 }
00247 }