G4INCL::DeltaDecayChannel Class Reference

#include <G4INCLDeltaDecayChannel.hh>

Inheritance diagram for G4INCL::DeltaDecayChannel:

G4INCL::IChannel

Public Member Functions

 DeltaDecayChannel (Nucleus *n, Particle *, ThreeVector const)
virtual ~DeltaDecayChannel ()
FinalStategetFinalState ()

Static Public Member Functions

static G4double computeDecayTime (Particle *p)

Detailed Description

Definition at line 47 of file G4INCLDeltaDecayChannel.hh.


Constructor & Destructor Documentation

G4INCL::DeltaDecayChannel::DeltaDecayChannel ( Nucleus n,
Particle ,
ThreeVector  const 
)

Definition at line 45 of file G4INCLDeltaDecayChannel.cc.

00046     :theParticle(p), theNucleus(n), incidentDirection(dir)
00047   { }

G4INCL::DeltaDecayChannel::~DeltaDecayChannel (  )  [virtual]

Definition at line 49 of file G4INCLDeltaDecayChannel.cc.

00049 {}


Member Function Documentation

G4double G4INCL::DeltaDecayChannel::computeDecayTime ( Particle p  )  [static]

Definition at line 51 of file G4INCLDeltaDecayChannel.cc.

References G4INCL::ParticleTable::effectiveNucleonMass, G4INCL::ParticleTable::effectivePionMass, G4INCL::Particle::getEnergy(), G4INCL::Particle::getMass(), G4INCL::PhysicalConstants::hc, G4INCL::KinematicsUtils::momentumInCM(), and G4INCL::Random::shoot().

Referenced by G4INCL::StandardPropagationModel::generateDecays().

00051                                                           {
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   }

FinalState * G4INCL::DeltaDecayChannel::getFinalState (  )  [virtual]

Implements G4INCL::IChannel.

Definition at line 74 of file G4INCLDeltaDecayChannel.cc.

References G4INCL::FinalState::addCreatedParticle(), G4INCL::FinalState::addModifiedParticle(), G4INCL::Particle::adjustEnergyFromMomentum(), G4INCL::DeltaMinus, G4INCL::DeltaPlus, G4INCL::DeltaPlusPlus, G4INCL::DeltaZero, FATAL, G4INCL::ParticleTable::getINCLMass(), G4INCL::Particle::getMass(), G4INCL::Particle::getPosition(), G4INCL::Particle::getType(), G4INCL::ThreeVector::getX(), G4INCL::ThreeVector::getY(), G4INCL::ThreeVector::getZ(), G4INCL::ThreeVector::mag(), G4INCL::KinematicsUtils::momentumInCM(), G4INCL::Neutron, G4INCL::ThreeVector::perp(), G4INCL::PiMinus, G4INCL::PiPlus, G4INCL::PiZero, G4INCL::Proton, G4INCL::Particle::setHelicity(), G4INCL::Particle::setMomentum(), G4INCL::Particle::setType(), G4INCL::Random::shoot(), and G4INCL::UnknownParticle.

00074                                                {
00075     //      SUBROUTINE DECAY2(P1,P2,P3,WP,ij,
00076     //     s       X1,X2,hel,B1,B2,B3)
00077 
00078     // This routine describes the anisotropic decay of a particle of mass
00079     // xi into 2 particles of masses x1,x2.
00080     // The anisotropy is supposed to follow a 1+3*hel*(cos(theta))**2
00081     // law with respect to the direction of the incoming particle.
00082     // In the input, p1,p2,p3 is the momentum of particle xi.
00083     // In the output, p1,p2,p3 is the momentum of particle x1 , while
00084     // q1,q2,q3 is the momentum of particle x2.
00085 
00086     //  COMMON/bl12/QQ1(200),QQ2(200),QQ3(200),QQ4(200),
00087     // s            YY1(200),YY2(200),YY3(200),YM(200),IPI(200)
00088     //   common/hazard/ial,IY1,IY2,IY3,IY4,IY5,IY6,IY7,IY8,IY9,IY10,
00089     // s               IY11,IY12,IY13,IY14,IY15,IY16,IY17,IY18,IY19
00090 
00091     // DATA IY8,IY9,IY10/82345,92345,45681/                           
00092     // PCM(E,A,C)=0.5*SQRT((E**2-(A+C)**2)*(E**2-(A-C)**2))/E            P-N20800
00093     // XI=YM(ij)
00094 
00095     // XE=WP                                                             P-N20810
00096     // B1=P1/XE                                                          P-N20820
00097     // B2=P2/XE                                                          P-N20830
00098     // B3=P3/XE                                                          
00099     // XQ=PCM(XI,X1,X2)                                                  
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     //      call loren(q1,q2,q3,b1,b2,b3,wq)
00183     //      call loren(p1,p2,p3,b1,b2,b3,wp)
00184     //      qq1(ij)=q1
00185     //      qq2(ij)=q2
00186     //      qq3(ij)=q3
00187     //      qq4(ij)=wq
00188     //      ym(ij)=xi
00189     //      RETURN                                                            P-N21120
00190     //      END                                                               P-N21130
00191     return fs;
00192   }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:54:06 2013 for Geant4 by  doxygen 1.4.7