G4INCL::DeltaProductionChannel Class Reference

#include <G4INCLDeltaProductionChannel.hh>

Inheritance diagram for G4INCL::DeltaProductionChannel:

G4INCL::IChannel

Public Member Functions

 DeltaProductionChannel (Particle *, Particle *, Nucleus *)
virtual ~DeltaProductionChannel ()
FinalStategetFinalState ()

Detailed Description

Definition at line 46 of file G4INCLDeltaProductionChannel.hh.


Constructor & Destructor Documentation

G4INCL::DeltaProductionChannel::DeltaProductionChannel ( Particle ,
Particle ,
Nucleus  
)

Definition at line 46 of file G4INCLDeltaProductionChannel.cc.

00049     :theNucleus(n), particle1(p1), particle2(p2)
00050   {}

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

Definition at line 52 of file G4INCLDeltaProductionChannel.cc.

00052 {}


Member Function Documentation

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   }


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