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