Geant4 Cross Reference |
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 "G4INCLDeltaProductionChannel.hh" 39 #include "G4INCLKinematicsUtils.hh" 40 #include "G4INCLBinaryCollisionAvatar.hh" 41 #include "G4INCLRandom.hh" 42 #include "G4INCLGlobals.hh" 43 #include "G4INCLLogger.hh" 44 45 namespace G4INCL { 46 47 const G4int DeltaProductionChannel::maxTries = 100000; 48 49 DeltaProductionChannel::DeltaProductionChannel(Particle *p1, 50 Particle *p2) 51 : particle1(p1), particle2(p2) 52 {} 53 54 DeltaProductionChannel::~DeltaProductionChannel() {} 55 56 G4double DeltaProductionChannel::sampleDeltaMass(G4double ecm) { 57 const G4double maxDeltaMass = ecm - ParticleTable::effectiveNucleonMass - 1.0; 58 const G4double maxDeltaMassRndm = std::atan((maxDeltaMass-ParticleTable::effectiveDeltaMass)*2./ParticleTable::effectiveDeltaWidth); 59 const G4double deltaMassRndmRange = maxDeltaMassRndm - ParticleTable::minDeltaMassRndm; 60 // assert(deltaMassRndmRange>0.); 61 62 G4double y=ecm*ecm; 63 G4double q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2 64 G4double q3=std::pow(std::sqrt(q2), 3.); 65 const G4double f3max=q3/(q3+5.832E6); // 5.832E6 = 180^3 66 G4double x; 67 68 G4int nTries = 0; 69 G4bool success = false; 70 while(!success) { /* Loop checking, 10.07.2015, D.Mancusi */ 71 if(++nTries >= maxTries) { 72 INCL_WARN("DeltaProductionChannel::sampleDeltaMass loop was stopped because maximum number of tries was reached. Minimum delta mass " 73 << ParticleTable::minDeltaMass << " MeV with CM energy " << ecm << " MeV may be unphysical." << '\n'); 74 return ParticleTable::minDeltaMass; 75 } 76 77 G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange; 78 y = std::tan(rndm); 79 x = ParticleTable::effectiveDeltaMass + 0.5*ParticleTable::effectiveDeltaWidth*y; 80 // assert(x>=ParticleTable::minDeltaMass && ecm >= x + ParticleTable::effectiveNucleonMass + 1.0); 81 82 // generation of the delta mass with the penetration factor 83 // (see prc56(1997)2431) 84 y=x*x; 85 q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2 86 q3=std::pow(std::sqrt(q2), 3.); 87 const G4double f3=q3/(q3+5.832E6); // 5.832E6 = 180^3 88 rndm = Random::shoot(); 89 if (rndm*f3max < f3) 90 success = true; 91 } 92 return x; 93 } 94 95 void DeltaProductionChannel::fillFinalState(FinalState *fs) { 96 /** 97 * Delta production 98 * 99 * The production is not isotropic in this version it has the same 100 * exp(b*t) structure as the nn elastic scattering (formula 2.3 of 101 * j.cugnon et al, nucl phys a352(1981)505) parametrization of b 102 * taken from ref. prc56(1997)2431 103 */ 104 // 100 IF (K4.NE.1) GO TO 101 // ThA K4 = 2 by default 105 // ParticleType p1TypeOld = particle1->getType(); 106 // ParticleType p2TypeOld = particle2->getType(); 107 G4double ecm = KinematicsUtils::totalEnergyInCM(particle1, particle2); 108 109 const G4int isospin = ParticleTable::getIsospin(particle1->getType()) + 110 ParticleTable::getIsospin(particle2->getType()); 111 112 // Calculate the outcome of the channel: 113 const ThreeVector &particle1Momentum = particle1->getMomentum(); 114 G4double pin = particle1Momentum.mag(); 115 G4double rndm = 0.0, b = 0.0; 116 117 G4double xmdel = sampleDeltaMass(ecm); 118 // deltaProduction103: // This label is not used 119 G4double pnorm = KinematicsUtils::momentumInCM(ecm, ParticleTable::effectiveNucleonMass, xmdel); 120 if (pnorm <= 0.0) pnorm=0.000001; 121 G4int index=0; 122 G4int index2=0; 123 rndm = Random::shoot(); 124 if (rndm < 0.5) index=1; 125 if (isospin == 0) { // pn case 126 rndm = Random::shoot(); 127 if (rndm < 0.5) index2=1; 128 } 129 130 // G4double x=0.001*0.5*ecm*std::sqrt(ecm*ecm-4.*ParticleTable::effectiveNucleonMass2) 131 // / ParticleTable::effectiveNucleonMass; 132 G4double x = 0.001 * KinematicsUtils::momentumInLab(ecm*ecm, ParticleTable::effectiveNucleonMass, ParticleTable::effectiveNucleonMass); 133 if(x < 1.4) { 134 b=(5.287/(1.+std::exp((1.3-x)/0.05)))*1.e-6; 135 } else { 136 b=(4.65+0.706*(x-1.4))*1.e-6; 137 } 138 G4double xkh = 2.*b*pin*pnorm; 139 rndm = Random::shoot(); 140 G4double ctet=1.0+std::log(1.-rndm*(1.-std::exp(-2.*xkh)))/xkh; 141 if(std::abs(ctet) > 1.0) ctet = Math::sign(ctet); 142 G4double stet = std::sqrt(1.-ctet*ctet); 143 144 rndm = Random::shoot(); 145 G4double fi = Math::twoPi*rndm; 146 G4double cfi = std::cos(fi); 147 G4double sfi = std::sin(fi); 148 // delta production: correction of the angular distribution 02/09/02 149 150 G4double xx = particle1Momentum.perp2(); 151 const G4double particle1MomentumZ = particle1Momentum.getZ(); 152 G4double zz = std::pow(particle1MomentumZ, 2); 153 G4double xp1, xp2, xp3; 154 if (xx >= zz*1.e-8) { 155 G4double yn = std::sqrt(xx); 156 G4double zn = yn*pin; 157 G4double ex[3], ey[3], ez[3]; 158 G4double p1 = particle1Momentum.getX(); 159 G4double p2 = particle1Momentum.getY(); 160 G4double p3 = particle1MomentumZ; 161 ez[0] = p1/pin; 162 ez[1] = p2/pin; 163 ez[2] = p3/pin; 164 ex[0] = p2/yn; 165 ex[1] = -p1/yn; 166 ex[2] = 0.0; 167 ey[0] = p1*p3/zn; 168 ey[1] = p2*p3/zn; 169 ey[2] = -xx/zn; 170 xp1 = (ex[0]*cfi*stet+ey[0]*sfi*stet+ez[0]*ctet)*pnorm; 171 xp2 = (ex[1]*cfi*stet+ey[1]*sfi*stet+ez[1]*ctet)*pnorm; 172 xp3 = (ex[2]*cfi*stet+ey[2]*sfi*stet+ez[2]*ctet)*pnorm; 173 }else { 174 xp1=pnorm*stet*cfi; 175 xp2=pnorm*stet*sfi; 176 xp3=pnorm*ctet; 177 } 178 // end of correction angular distribution of delta production 179 G4double e3 = std::sqrt(xp1*xp1+xp2*xp2+xp3*xp3 180 +ParticleTable::effectiveNucleonMass2); 181 // if(k4.ne.0) go to 161 182 183 // long-lived delta 184 if (index != 1) { 185 ThreeVector mom(xp1, xp2, xp3); 186 particle1->setMomentum(mom); 187 // e1=ecm-eout1 188 } else { 189 ThreeVector mom(-xp1, -xp2, -xp3); 190 particle1->setMomentum(mom); 191 // e1=ecm-eout1 192 } 193 194 particle1->setEnergy(ecm - e3); 195 particle2->setEnergy(e3); 196 particle2->setMomentum(-particle1->getMomentum()); 197 198 // SYMMETRIZATION OF CHARGES IN pn -> N DELTA 199 // THE TEST ON "INDEX" ABOVE SYMETRIZES THE EXCITATION OF ONE 200 // OF THE NUCLEONS WITH RESPECT TO THE DELTA EXCITATION 201 // (SEE NOTE 16/10/97) 202 G4int is1 = ParticleTable::getIsospin(particle1->getType()); 203 G4int is2 = ParticleTable::getIsospin(particle2->getType()); 204 if (isospin == 0) { 205 if(index2 == 1) { 206 G4int isi=is1; 207 is1=is2; 208 is2=isi; 209 } 210 particle1->setHelicity(0.0); 211 } else { 212 rndm = Random::shoot(); 213 if (rndm >= 0.25) { 214 is1=3*is1; 215 is2=-is2; 216 } 217 particle1->setHelicity(ctet*ctet); 218 } 219 220 if(is1 == ParticleTable::getIsospin(DeltaMinus)) { 221 particle1->setType(DeltaMinus); 222 } else if(is1 == ParticleTable::getIsospin(DeltaZero)) { 223 particle1->setType(DeltaZero); 224 } else if(is1 == ParticleTable::getIsospin(DeltaPlus)) { 225 particle1->setType(DeltaPlus); 226 } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus)) { 227 particle1->setType(DeltaPlusPlus); 228 } 229 230 if(is2 == ParticleTable::getIsospin(Proton)) { 231 particle2->setType(Proton); 232 } else if(is2 == ParticleTable::getIsospin(Neutron)) { 233 particle2->setType(Neutron); 234 } 235 236 if(particle1->isDelta()) particle1->setMass(xmdel); 237 if(particle2->isDelta()) particle2->setMass(xmdel); 238 239 fs->addModifiedParticle(particle1); 240 fs->addModifiedParticle(particle2); 241 } 242 } 243