G4INCLDeltaDecayChannel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // INCL++ intra-nuclear cascade model
00027 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
00028 // Davide Mancusi, CEA
00029 // Alain Boudard, CEA
00030 // Sylvie Leray, CEA
00031 // Joseph Cugnon, University of Liege
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     //      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   }
00193 }

Generated on Mon May 27 17:48:34 2013 for Geant4 by  doxygen 1.4.7