Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLPionResonanceDecayChannel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // INCL++ intra-nuclear cascade model
 27 // Alain Boudard, CEA-Saclay, France
 28 // Joseph Cugnon, University of Liege, Belgium
 29 // Jean-Christophe David, CEA-Saclay, France
 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
 31 // Sylvie Leray, CEA-Saclay, France
 32 // Davide Mancusi, CEA-Saclay, France
 33 //
 34 #define INCLXX_IN_GEANT4_MODE 1
 35 
 36 #include "globals.hh"
 37 
 38 #include "G4INCLPionResonanceDecayChannel.hh"
 39 #include "G4INCLKinematicsUtils.hh"
 40 #include "G4INCLBinaryCollisionAvatar.hh"
 41 #include "G4INCLRandom.hh"
 42 #include "G4INCLGlobals.hh"
 43 #include "G4INCLPhaseSpaceGenerator.hh"
 44 // #include <cassert>
 45 
 46 
 47 namespace G4INCL {
 48   
 49   const G4double PionResonanceDecayChannel::angularSlope = 0.; // Slope to be defined, if needed.
 50   
 51   PionResonanceDecayChannel::PionResonanceDecayChannel(Particle *p, ThreeVector const &dir)
 52     :theParticle(p), incidentDirection(dir)
 53   { }
 54   
 55   PionResonanceDecayChannel::~PionResonanceDecayChannel() {}
 56   
 57 //  Decay during the intranuclear cascade for Omega only
 58   G4double PionResonanceDecayChannel::computeDecayTime(Particle *p) {
 59   const G4double m = p->getMass();
 60   const G4double geff = p->getEnergy()/m;
 61 //    const G4double geta = 1.31e-3;
 62     const G4double gomega = 8.49;
 63     G4double gg=0.;
 64     switch (p->getType()) {
 65 /*      case Eta:
 66         gg=geta;
 67         break;*/
 68       case Omega:
 69         gg=gomega;
 70         break;
 71       default:
 72         INCL_FATAL("Unrecognized pion resonance type; type=" << p->getType() << '\n');
 73         break;
 74     }
 75     const G4double tpires = -G4INCL::PhysicalConstants::hc/gg*std::log(Random::shoot())*geff;
 76     return tpires;
 77   }
 78   
 79   void PionResonanceDecayChannel::sampleAngles(G4double *ctet_par, G4double *stet_par, G4double *phi_par) {
 80 
 81     (*ctet_par) = -1.0 + 2.0*Random::shoot();
 82     if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
 83     (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
 84     (*phi_par) = Math::twoPi * Random::shoot();
 85   }
 86   
 87   void PionResonanceDecayChannel::fillFinalState(FinalState *fs) {
 88     
 89     ParticleType createdType;
 90     ParticleType pionType1=Neutron; // to avoid forgetting pionType definition when 3 particles are emitted
 91     ParticleType pionType2=Neutron; 
 92 
 93     const G4double sqrtS = theParticle->getMass();
 94     G4int nbpart = 3; // number of emitted particles
 95     G4double drnd=Random::shoot();
 96     switch (theParticle->getType()) {
 97       case Eta:
 98         if (drnd < 0.3972) { // renormalized to the only four decays taken into account here
 99 //    2 photons
100           nbpart=2;
101           theParticle->setType(Photon);
102           createdType = Photon;
103         }
104         else if (drnd < 0.7265) {
105 //    3 pi0
106           theParticle->setType(PiZero);
107           pionType1 = PiZero;
108           pionType2 = PiZero;
109         }
110         else if (drnd < 0.9575) {
111 //    pi+ pi- pi0
112         theParticle->setType(PiZero);
113         pionType1 = PiPlus;
114         pionType2 = PiMinus;
115       }
116         else {
117 //    pi+ pi- photon
118           theParticle->setType(Photon);
119           pionType1 = PiPlus;
120           pionType2 = PiMinus;
121         }
122       break;
123       case Omega:
124         if (drnd < 0.9009) { // renormalized to the only three decays taken into account here
125 //    pi+ pi- pi0
126           theParticle->setType(PiZero);
127           pionType1 = PiPlus;
128           pionType2 = PiMinus;
129         }
130         else if (drnd < 0.9845) {
131 //    pi0 photon
132           nbpart=2;
133           theParticle->setType(PiZero);
134           createdType = Photon;
135         }
136         else {
137 //    pi+ pi-
138           nbpart=2;
139           theParticle->setType(PiPlus);
140           createdType = PiMinus;
141         }
142         break;
143       default:
144         INCL_FATAL("Unrecognized pion resonance type; type=" << theParticle->getType() << '\n');
145         break;
146     }
147 
148     if (nbpart == 2) {      
149       G4double fi, ctet, stet;
150       sampleAngles(&ctet, &stet, &fi);
151       
152       G4double cfi = std::cos(fi);
153       G4double sfi = std::sin(fi);
154       G4double beta = incidentDirection.mag();
155       
156       G4double q1, q2, q3;
157       G4double sal=0.0;
158       if (beta >= 1.0e-10)
159         sal = incidentDirection.perp()/beta;
160       if (sal >= 1.0e-6) {
161         G4double b1 = incidentDirection.getX();
162         G4double b2 = incidentDirection.getY();
163         G4double b3 = incidentDirection.getZ();
164         G4double cal = b3/beta;
165         G4double t1 = ctet+cal*stet*sfi/sal;
166         G4double t2 = stet/sal;
167         q1=(b1*t1+b2*t2*cfi)/beta;
168         q2=(b2*t1-b1*t2*cfi)/beta;
169         q3=(b3*t1/beta-t2*sfi);
170       } else {
171         q1 = stet*cfi;
172         q2 = stet*sfi;
173         q3 = ctet;
174       }
175             
176       G4double xq = KinematicsUtils::momentumInCM(sqrtS,
177                             theParticle->getMass(),
178                             ParticleTable::getINCLMass(createdType));
179       q1 *= xq;
180       q2 *= xq;
181       q3 *= xq;
182       
183       ThreeVector createdMomentum(q1, q2, q3);
184       ThreeVector createdPosition(theParticle->getPosition());
185       Particle *createdParticle = new Particle(createdType, createdMomentum, createdPosition);
186       theParticle->setMomentum(-createdMomentum);
187       theParticle->adjustEnergyFromMomentum();
188       
189       fs->addModifiedParticle(theParticle);
190       fs->addCreatedParticle(createdParticle);
191       
192     }
193     else if (nbpart == 3) {
194 // assert(pionType1!=Neutron && pionType2!=Neutron);
195       ParticleList list;
196       list.push_back(theParticle);
197       const ThreeVector &rposdecay = theParticle->getPosition();
198       const ThreeVector zero;
199       Particle *Pion1 = new Particle(pionType1,zero,rposdecay);
200       Particle *Pion2 = new Particle(pionType2,zero,rposdecay);
201       list.push_back(Pion1);
202       list.push_back(Pion2);
203       
204       fs->addModifiedParticle(theParticle);
205       fs->addCreatedParticle(Pion1);
206       fs->addCreatedParticle(Pion2);
207 
208       //      PhaseSpaceGenerator::generateBiased(sqrtS, list, 0, angularSlope); Biasing?
209       PhaseSpaceGenerator::generate(sqrtS, list);
210     }
211 
212   }
213 }
214