#include <G4INCLDeltaProductionChannel.hh>
Inheritance diagram for G4INCL::DeltaProductionChannel:
Public Member Functions | |
DeltaProductionChannel (Particle *, Particle *, Nucleus *) | |
virtual | ~DeltaProductionChannel () |
FinalState * | getFinalState () |
Definition at line 46 of file G4INCLDeltaProductionChannel.hh.
Definition at line 46 of file G4INCLDeltaProductionChannel.cc.
00049 :theNucleus(n), particle1(p1), particle2(p2) 00050 {}
G4INCL::DeltaProductionChannel::~DeltaProductionChannel | ( | ) | [virtual] |
FinalState * G4INCL::DeltaProductionChannel::getFinalState | ( | ) | [virtual] |
Implements G4INCL::IChannel.
Definition at line 84 of file G4INCLDeltaProductionChannel.cc.
References G4INCL::FinalState::addModifiedParticle(), G4INCL::DeltaMinus, G4INCL::DeltaPlus, G4INCL::DeltaPlusPlus, G4INCL::DeltaZero, G4INCL::ParticleTable::effectiveNucleonMass, G4INCL::ParticleTable::effectiveNucleonMass2, G4INCL::ParticleTable::getIsospin(), G4INCL::Particle::getMomentum(), G4INCL::Particle::getType(), G4INCL::ThreeVector::getX(), G4INCL::ThreeVector::getY(), G4INCL::ThreeVector::getZ(), G4INCL::Particle::isDelta(), G4INCL::ThreeVector::mag(), G4INCL::KinematicsUtils::momentumInCM(), G4INCL::KinematicsUtils::momentumInLab(), G4INCL::Neutron, G4INCL::ThreeVector::perp2(), G4INCL::Proton, G4INCL::Particle::setEnergy(), G4INCL::Particle::setHelicity(), G4INCL::Particle::setMass(), G4INCL::Particle::setMomentum(), G4INCL::Particle::setType(), G4INCL::Random::shoot(), G4INCL::Math::sign(), G4INCL::KinematicsUtils::totalEnergyInCM(), and G4INCL::Math::twoPi.
00084 { 00093 // 100 IF (K4.NE.1) GO TO 101 // ThA K4 = 2 by default 00094 // ParticleType p1TypeOld = particle1->getType(); 00095 // ParticleType p2TypeOld = particle2->getType(); 00096 G4double ecm = KinematicsUtils::totalEnergyInCM(particle1, particle2); 00097 00098 const G4int isospin = ParticleTable::getIsospin(particle1->getType()) + 00099 ParticleTable::getIsospin(particle2->getType()); 00100 00101 // Calculate the outcome of the channel: 00102 G4double pin = particle1->getMomentum().mag(); 00103 G4double rndm = 0.0, b = 0.0; 00104 00105 G4double xmdel = sampleDeltaMass(ecm); 00106 // deltaProduction103: // This label is not used 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) { // pn case 00114 rndm = Random::shoot(); 00115 if (rndm < 0.5) index2=1; 00116 } 00117 00118 // G4double x=0.001*0.5*ecm*std::sqrt(ecm*ecm-4.*ParticleTable::effectiveNucleonMass2) 00119 // / ParticleTable::effectiveNucleonMass; 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 // delta production: correction of the angular distribution 02/09/02 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 // end of correction angular distribution of delta production 00166 G4double e3 = std::sqrt(xp1*xp1+xp2*xp2+xp3*xp3 00167 +ParticleTable::effectiveNucleonMass2); 00168 // if(k4.ne.0) go to 161 00169 00170 // long-lived delta 00171 G4int m1 = 0; 00172 G4int m2 = 0; 00173 if (index != 1) { 00174 ThreeVector mom(xp1, xp2, xp3); 00175 particle1->setMomentum(mom); 00176 // e1=ecm-eout1 00177 m1=1; 00178 } else { 00179 ThreeVector mom(-xp1, -xp2, -xp3); 00180 particle1->setMomentum(mom); 00181 // e1=ecm-eout1 00182 m1=1; 00183 } 00184 00185 particle1->setEnergy(ecm - e3); 00186 particle2->setEnergy(e3); 00187 particle2->setMomentum(-particle1->getMomentum()); 00188 00189 // SYMMETRIZATION OF CHARGES IN pn -> N DELTA 00190 // THE TEST ON "INDEX" ABOVE SYMETRIZES THE EXCITATION OF ONE 00191 // OF THE NUCLEONS WITH RESPECT TO THE DELTA EXCITATION 00192 // (SEE NOTE 16/10/97) 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 }