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 "G4INCLDeltaDecayChannel.hh"
00038 #include "G4INCLKinematicsUtils.hh"
00039 #include "G4INCLBinaryCollisionAvatar.hh"
00040 #include "G4INCLRandom.hh"
00041 #include "G4INCLGlobals.hh"
00042
00043 namespace G4INCL {
00044
00045 DeltaDecayChannel::DeltaDecayChannel(Nucleus *n, Particle *p, ThreeVector const dir)
00046 :theParticle(p), theNucleus(n), incidentDirection(dir)
00047 { }
00048
00049 DeltaDecayChannel::~DeltaDecayChannel() {}
00050
00051 G4double DeltaDecayChannel::computeDecayTime(Particle *p) {
00052 const G4double m = p->getMass();
00053 const G4double g0 = 115.0;
00054 G4double gg = g0;
00055 if(m > 1500.0) gg = 200.0;
00056 const G4double geff = p->getEnergy()/m;
00057 const G4double qqq = KinematicsUtils::momentumInCM(m, ParticleTable::effectiveNucleonMass, ParticleTable::effectivePionMass);
00058 const G4double psf = std::pow(qqq, 3)/(std::pow(qqq, 3) + 5832000.0);
00059 const G4double tdel = -G4INCL::PhysicalConstants::hc/(gg*psf)*std::log(Random::shoot())*geff;
00060 return tdel;
00061 }
00062
00063 void DeltaDecayChannel::sampleAngles(G4double *ctet_par, G4double *stet_par, G4double *phi_par) {
00064 const G4double hel = theParticle->getHelicity();
00065 do {
00066 (*ctet_par) = -1.0 + 2.0*Random::shoot();
00067 if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
00068 } while(Random::shoot() > ((1.0 + 3.0 * hel * (*ctet_par) * (*ctet_par))
00069 /(1.0 + 3.0 * hel)));
00070 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
00071 (*phi_par) = Math::twoPi * Random::shoot();
00072 }
00073
00074 FinalState* DeltaDecayChannel::getFinalState() {
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101 const G4double deltaMass = theParticle->getMass();
00102
00103 G4double fi, ctet, stet;
00104 sampleAngles(&ctet, &stet, &fi);
00105
00106 G4double cfi = std::cos(fi);
00107 G4double sfi = std::sin(fi);
00108 G4double beta = incidentDirection.mag();
00109
00110 G4double q1, q2, q3;
00111 G4double sal=0.0;
00112 if (beta >= 1.0e-10)
00113 sal = incidentDirection.perp()/beta;
00114 if (sal >= 1.0e-6) {
00115 G4double b1 = incidentDirection.getX();
00116 G4double b2 = incidentDirection.getY();
00117 G4double b3 = incidentDirection.getZ();
00118 G4double cal = b3/beta;
00119 G4double t1 = ctet+cal*stet*sfi/sal;
00120 G4double t2 = stet/sal;
00121 q1=(b1*t1+b2*t2*cfi)/beta;
00122 q2=(b2*t1-b1*t2*cfi)/beta;
00123 q3=(b3*t1/beta-t2*sfi);
00124 } else {
00125 q1 = stet*cfi;
00126 q2 = stet*sfi;
00127 q3 = ctet;
00128 }
00129 theParticle->setHelicity(0.0);
00130
00131 ParticleType pionType;
00132 switch(theParticle->getType()) {
00133 case DeltaPlusPlus:
00134 theParticle->setType(Proton);
00135 pionType = PiPlus;
00136 break;
00137 case DeltaPlus:
00138 if(Random::shoot() < 1.0/3.0) {
00139 theParticle->setType(Neutron);
00140 pionType = PiPlus;
00141 } else {
00142 theParticle->setType(Proton);
00143 pionType = PiZero;
00144 }
00145 break;
00146 case DeltaZero:
00147 if(Random::shoot() < 1.0/3.0) {
00148 theParticle->setType(Proton);
00149 pionType = PiMinus;
00150 } else {
00151 theParticle->setType(Neutron);
00152 pionType = PiZero;
00153 }
00154 break;
00155 case DeltaMinus:
00156 theParticle->setType(Neutron);
00157 pionType = PiMinus;
00158 break;
00159 default:
00160 FATAL("Unrecognized delta type; type=" << theParticle->getType() << std::endl);
00161 pionType = UnknownParticle;
00162 break;
00163 }
00164
00165 G4double xq = KinematicsUtils::momentumInCM(deltaMass,
00166 theParticle->getMass(),
00167 ParticleTable::getINCLMass(pionType));
00168
00169 q1 *= xq;
00170 q2 *= xq;
00171 q3 *= xq;
00172
00173 ThreeVector pionMomentum(q1, q2, q3);
00174 ThreeVector pionPosition(theParticle->getPosition());
00175 Particle *pion = new Particle(pionType, pionMomentum, pionPosition);
00176 theParticle->setMomentum(-pionMomentum);
00177 theParticle->adjustEnergyFromMomentum();
00178
00179 FinalState *fs = new FinalState;
00180 fs->addModifiedParticle(theParticle);
00181 fs->addCreatedParticle(pion);
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 return fs;
00192 }
00193 }