| 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 | //
|
|---|
| 27 | // $Id: G4HEAntiProtonInelastic.cc,v 1.14 2008/03/17 20:49:17 dennis Exp $
|
|---|
| 28 | // GEANT4 tag $Name: geant4-09-02 $
|
|---|
| 29 | //
|
|---|
| 30 | //
|
|---|
| 31 |
|
|---|
| 32 | #include "globals.hh"
|
|---|
| 33 | #include "G4ios.hh"
|
|---|
| 34 |
|
|---|
| 35 | //
|
|---|
| 36 | // G4 Process: Gheisha High Energy Collision model.
|
|---|
| 37 | // This includes the high energy cascading model, the two-body-resonance model
|
|---|
| 38 | // and the low energy two-body model. Not included are the low energy stuff like
|
|---|
| 39 | // nuclear reactions, nuclear fission without any cascading and all processes for
|
|---|
| 40 | // particles at rest.
|
|---|
| 41 | // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
|
|---|
| 42 | // H. Fesefeldt, RWTH-Aachen, 23-October-1996
|
|---|
| 43 | // Last modified: 29-July-1998
|
|---|
| 44 |
|
|---|
| 45 | #include "G4HEAntiProtonInelastic.hh"
|
|---|
| 46 |
|
|---|
| 47 | G4HadFinalState * G4HEAntiProtonInelastic::
|
|---|
| 48 | ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &targetNucleus )
|
|---|
| 49 | {
|
|---|
| 50 | G4HEVector * pv = new G4HEVector[MAXPART];
|
|---|
| 51 | const G4HadProjectile *aParticle = &aTrack;
|
|---|
| 52 | // G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
|
|---|
| 53 | const G4double atomicWeight = targetNucleus.GetN();
|
|---|
| 54 | const G4double atomicNumber = targetNucleus.GetZ();
|
|---|
| 55 | G4HEVector incidentParticle(aParticle);
|
|---|
| 56 |
|
|---|
| 57 | G4int incidentCode = incidentParticle.getCode();
|
|---|
| 58 | G4double incidentMass = incidentParticle.getMass();
|
|---|
| 59 | G4double incidentTotalEnergy = incidentParticle.getEnergy();
|
|---|
| 60 | G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
|
|---|
| 61 | G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
|
|---|
| 62 |
|
|---|
| 63 | if(incidentKineticEnergy < 1.)
|
|---|
| 64 | {
|
|---|
| 65 | G4cout << "GHEAntiProtonInelastic: incident energy < 1 GeV" << G4endl;
|
|---|
| 66 | }
|
|---|
| 67 | if(verboseLevel > 1)
|
|---|
| 68 | {
|
|---|
| 69 | G4cout << "G4HEAntiProtonInelastic::ApplyYourself" << G4endl;
|
|---|
| 70 | G4cout << "incident particle " << incidentParticle.getName()
|
|---|
| 71 | << "mass " << incidentMass
|
|---|
| 72 | << "kinetic energy " << incidentKineticEnergy
|
|---|
| 73 | << G4endl;
|
|---|
| 74 | G4cout << "target material with (A,Z) = ("
|
|---|
| 75 | << atomicWeight << "," << atomicNumber << ")" << G4endl;
|
|---|
| 76 | }
|
|---|
| 77 |
|
|---|
| 78 | G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
|
|---|
| 79 | atomicWeight, atomicNumber);
|
|---|
| 80 | if(verboseLevel > 1)
|
|---|
| 81 | G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
|
|---|
| 82 |
|
|---|
| 83 |
|
|---|
| 84 | incidentKineticEnergy -= inelasticity;
|
|---|
| 85 |
|
|---|
| 86 | G4double excitationEnergyGNP = 0.;
|
|---|
| 87 | G4double excitationEnergyDTA = 0.;
|
|---|
| 88 |
|
|---|
| 89 | G4double excitation = NuclearExcitation(incidentKineticEnergy,
|
|---|
| 90 | atomicWeight, atomicNumber,
|
|---|
| 91 | excitationEnergyGNP,
|
|---|
| 92 | excitationEnergyDTA);
|
|---|
| 93 | if(verboseLevel > 1)
|
|---|
| 94 | G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
|
|---|
| 95 | << excitationEnergyDTA << G4endl;
|
|---|
| 96 |
|
|---|
| 97 |
|
|---|
| 98 | incidentKineticEnergy -= excitation;
|
|---|
| 99 | incidentTotalEnergy = incidentKineticEnergy + incidentMass;
|
|---|
| 100 | incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
|
|---|
| 101 | *(incidentTotalEnergy+incidentMass));
|
|---|
| 102 |
|
|---|
| 103 | G4HEVector targetParticle;
|
|---|
| 104 | if(G4UniformRand() < atomicNumber/atomicWeight)
|
|---|
| 105 | {
|
|---|
| 106 | targetParticle.setDefinition("Proton");
|
|---|
| 107 | }
|
|---|
| 108 | else
|
|---|
| 109 | {
|
|---|
| 110 | targetParticle.setDefinition("Neutron");
|
|---|
| 111 | }
|
|---|
| 112 |
|
|---|
| 113 | G4double targetMass = targetParticle.getMass();
|
|---|
| 114 | G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass
|
|---|
| 115 | + 2.0*targetMass*incidentTotalEnergy);
|
|---|
| 116 | G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
|
|---|
| 117 |
|
|---|
| 118 | // this was the meaning of inElastic in the
|
|---|
| 119 | // original Gheisha stand-alone version.
|
|---|
| 120 | // G4bool inElastic = InElasticCrossSectionInFirstInt
|
|---|
| 121 | // (availableEnergy, incidentCode, incidentTotalMomentum);
|
|---|
| 122 | // by unknown reasons, it has been replaced
|
|---|
| 123 | // to the following code in Geant???
|
|---|
| 124 | G4bool inElastic = true;
|
|---|
| 125 | // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;
|
|---|
| 126 |
|
|---|
| 127 | vecLength = 0;
|
|---|
| 128 |
|
|---|
| 129 | if(verboseLevel > 1)
|
|---|
| 130 | G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
|
|---|
| 131 | << incidentCode << G4endl;
|
|---|
| 132 |
|
|---|
| 133 | G4bool successful = false;
|
|---|
| 134 |
|
|---|
| 135 | if(inElastic || (!inElastic && atomicWeight < 1.5))
|
|---|
| 136 | {
|
|---|
| 137 | FirstIntInCasAntiProton(inElastic, availableEnergy, pv, vecLength,
|
|---|
| 138 | incidentParticle, targetParticle, atomicWeight);
|
|---|
| 139 |
|
|---|
| 140 | if(verboseLevel > 1)
|
|---|
| 141 | G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
|
|---|
| 142 |
|
|---|
| 143 |
|
|---|
| 144 | if ((vecLength > 0) && (availableEnergy > 1.))
|
|---|
| 145 | StrangeParticlePairProduction( availableEnergy, centerOfMassEnergy,
|
|---|
| 146 | pv, vecLength,
|
|---|
| 147 | incidentParticle, targetParticle);
|
|---|
| 148 | HighEnergyCascading( successful, pv, vecLength,
|
|---|
| 149 | excitationEnergyGNP, excitationEnergyDTA,
|
|---|
| 150 | incidentParticle, targetParticle,
|
|---|
| 151 | atomicWeight, atomicNumber);
|
|---|
| 152 | if (!successful)
|
|---|
| 153 | HighEnergyClusterProduction( successful, pv, vecLength,
|
|---|
| 154 | excitationEnergyGNP, excitationEnergyDTA,
|
|---|
| 155 | incidentParticle, targetParticle,
|
|---|
| 156 | atomicWeight, atomicNumber);
|
|---|
| 157 | if (!successful)
|
|---|
| 158 | MediumEnergyCascading( successful, pv, vecLength,
|
|---|
| 159 | excitationEnergyGNP, excitationEnergyDTA,
|
|---|
| 160 | incidentParticle, targetParticle,
|
|---|
| 161 | atomicWeight, atomicNumber);
|
|---|
| 162 |
|
|---|
| 163 | if (!successful)
|
|---|
| 164 | MediumEnergyClusterProduction( successful, pv, vecLength,
|
|---|
| 165 | excitationEnergyGNP, excitationEnergyDTA,
|
|---|
| 166 | incidentParticle, targetParticle,
|
|---|
| 167 | atomicWeight, atomicNumber);
|
|---|
| 168 | if (!successful)
|
|---|
| 169 | QuasiElasticScattering( successful, pv, vecLength,
|
|---|
| 170 | excitationEnergyGNP, excitationEnergyDTA,
|
|---|
| 171 | incidentParticle, targetParticle,
|
|---|
| 172 | atomicWeight, atomicNumber);
|
|---|
| 173 | }
|
|---|
| 174 | if (!successful)
|
|---|
| 175 | {
|
|---|
| 176 | ElasticScattering( successful, pv, vecLength,
|
|---|
| 177 | incidentParticle,
|
|---|
| 178 | atomicWeight, atomicNumber);
|
|---|
| 179 | }
|
|---|
| 180 |
|
|---|
| 181 | if (!successful)
|
|---|
| 182 | {
|
|---|
| 183 | G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" << G4endl;
|
|---|
| 184 | }
|
|---|
| 185 | FillParticleChange(pv, vecLength);
|
|---|
| 186 | delete [] pv;
|
|---|
| 187 | theParticleChange.SetStatusChange(stopAndKill);
|
|---|
| 188 | return & theParticleChange;
|
|---|
| 189 | }
|
|---|
| 190 |
|
|---|
| 191 | void
|
|---|
| 192 | G4HEAntiProtonInelastic::FirstIntInCasAntiProton( G4bool &inElastic,
|
|---|
| 193 | const G4double availableEnergy,
|
|---|
| 194 | G4HEVector pv[],
|
|---|
| 195 | G4int &vecLen,
|
|---|
| 196 | G4HEVector incidentParticle,
|
|---|
| 197 | G4HEVector targetParticle,
|
|---|
| 198 | const G4double atomicWeight)
|
|---|
| 199 |
|
|---|
| 200 | // AntiProton undergoes interaction with nucleon within a nucleus. Check if it is
|
|---|
| 201 | // energetically possible to produce pions/kaons. In not, assume nuclear excitation
|
|---|
| 202 | // occurs and input particle is degraded in energy. No other particles are produced.
|
|---|
| 203 | // If reaction is possible, find the correct number of pions/protons/neutrons
|
|---|
| 204 | // produced using an interpolation to multiplicity data. Replace some pions or
|
|---|
| 205 | // protons/neutrons by kaons or strange baryons according to the average
|
|---|
| 206 | // multiplicity per inelastic reaction.
|
|---|
| 207 |
|
|---|
| 208 | {
|
|---|
| 209 | static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp
|
|---|
| 210 | static const G4double expxl = -expxu; // lower bound for arg. of exp
|
|---|
| 211 |
|
|---|
| 212 | static const G4double protb = 0.7;
|
|---|
| 213 | static const G4double neutb = 0.7;
|
|---|
| 214 | static const G4double c = 1.25;
|
|---|
| 215 |
|
|---|
| 216 | static const G4int numMul = 1200;
|
|---|
| 217 | static const G4int numMulAn = 400;
|
|---|
| 218 | static const G4int numSec = 60;
|
|---|
| 219 |
|
|---|
| 220 | G4int neutronCode = Neutron.getCode();
|
|---|
| 221 | G4int protonCode = Proton.getCode();
|
|---|
| 222 |
|
|---|
| 223 | G4int targetCode = targetParticle.getCode();
|
|---|
| 224 | // G4double incidentMass = incidentParticle.getMass();
|
|---|
| 225 | // G4double incidentEnergy = incidentParticle.getEnergy();
|
|---|
| 226 | G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
|
|---|
| 227 |
|
|---|
| 228 | static G4bool first = true;
|
|---|
| 229 | static G4double protmul[numMul], protnorm[numSec]; // proton constants
|
|---|
| 230 | static G4double protmulAn[numMulAn],protnormAn[numSec];
|
|---|
| 231 | static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
|
|---|
| 232 | static G4double neutmulAn[numMulAn],neutnormAn[numSec];
|
|---|
| 233 |
|
|---|
| 234 | // misc. local variables
|
|---|
| 235 | // np = number of pi+, nm = number of pi-, nz = number of pi0
|
|---|
| 236 |
|
|---|
| 237 | G4int i, counter, nt, np, nm, nz;
|
|---|
| 238 |
|
|---|
| 239 | if( first )
|
|---|
| 240 | { // compute normalization constants, this will only be done once
|
|---|
| 241 | first = false;
|
|---|
| 242 | for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
|
|---|
| 243 | for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
|
|---|
| 244 | counter = -1;
|
|---|
| 245 | for( np=0; np<(numSec/3); np++ )
|
|---|
| 246 | {
|
|---|
| 247 | for( nm=Imax(0,np-1); nm<=(np+1); nm++ )
|
|---|
| 248 | {
|
|---|
| 249 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 250 | {
|
|---|
| 251 | if( ++counter < numMul )
|
|---|
| 252 | {
|
|---|
| 253 | nt = np+nm+nz;
|
|---|
| 254 | if( (nt>0) && (nt<=numSec) )
|
|---|
| 255 | {
|
|---|
| 256 | protmul[counter] = pmltpc(np,nm,nz,nt,protb,c);
|
|---|
| 257 | protnorm[nt-1] += protmul[counter];
|
|---|
| 258 | }
|
|---|
| 259 | }
|
|---|
| 260 | }
|
|---|
| 261 | }
|
|---|
| 262 | }
|
|---|
| 263 | for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
|
|---|
| 264 |
|
|---|
| 265 | for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
|
|---|
| 266 | counter = -1;
|
|---|
| 267 | for( np=0; np<numSec/3; np++ )
|
|---|
| 268 | {
|
|---|
| 269 | for( nm=np; nm<=(np+2); nm++ )
|
|---|
| 270 | {
|
|---|
| 271 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 272 | {
|
|---|
| 273 | if( ++counter < numMul )
|
|---|
| 274 | {
|
|---|
| 275 | nt = np+nm+nz;
|
|---|
| 276 | if( (nt>0) && (nt<=numSec) )
|
|---|
| 277 | {
|
|---|
| 278 | neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c);
|
|---|
| 279 | neutnorm[nt-1] += neutmul[counter];
|
|---|
| 280 | }
|
|---|
| 281 | }
|
|---|
| 282 | }
|
|---|
| 283 | }
|
|---|
| 284 | }
|
|---|
| 285 | for( i=0; i<numSec; i++ )
|
|---|
| 286 | {
|
|---|
| 287 | if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
|
|---|
| 288 | if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
|
|---|
| 289 | }
|
|---|
| 290 | // annihilation
|
|---|
| 291 | for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
|
|---|
| 292 | for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
|
|---|
| 293 | counter = -1;
|
|---|
| 294 | for( np=1; np<(numSec/3); np++ )
|
|---|
| 295 | {
|
|---|
| 296 | nm = np;
|
|---|
| 297 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 298 | {
|
|---|
| 299 | if( ++counter < numMulAn )
|
|---|
| 300 | {
|
|---|
| 301 | nt = np+nm+nz;
|
|---|
| 302 | if( (nt>0) && (nt<=numSec) )
|
|---|
| 303 | {
|
|---|
| 304 | protmulAn[counter] = pmltpc(np,nm,nz,nt,protb,c);
|
|---|
| 305 | protnormAn[nt-1] += protmulAn[counter];
|
|---|
| 306 | }
|
|---|
| 307 | }
|
|---|
| 308 | }
|
|---|
| 309 | }
|
|---|
| 310 | for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
|
|---|
| 311 | for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
|
|---|
| 312 | counter = -1;
|
|---|
| 313 | for( np=1; np<numSec/3; np++ )
|
|---|
| 314 | {
|
|---|
| 315 | nm = np+1;
|
|---|
| 316 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 317 | {
|
|---|
| 318 | if( ++counter < numMulAn )
|
|---|
| 319 | {
|
|---|
| 320 | nt = np+nm+nz;
|
|---|
| 321 | if( (nt>0) && (nt<=numSec) )
|
|---|
| 322 | {
|
|---|
| 323 | neutmulAn[counter] = pmltpc(np,nm,nz,nt,neutb,c);
|
|---|
| 324 | neutnormAn[nt-1] += neutmulAn[counter];
|
|---|
| 325 | }
|
|---|
| 326 | }
|
|---|
| 327 | }
|
|---|
| 328 | }
|
|---|
| 329 | for( i=0; i<numSec; i++ )
|
|---|
| 330 | {
|
|---|
| 331 | if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
|
|---|
| 332 | if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
|
|---|
| 333 | }
|
|---|
| 334 | } // end of initialization
|
|---|
| 335 |
|
|---|
| 336 |
|
|---|
| 337 | // initialize the first two places
|
|---|
| 338 | // the same as beam and target
|
|---|
| 339 | pv[0] = incidentParticle;
|
|---|
| 340 | pv[1] = targetParticle;
|
|---|
| 341 | vecLen = 2;
|
|---|
| 342 |
|
|---|
| 343 | if( !inElastic )
|
|---|
| 344 | { // pb p --> nb n
|
|---|
| 345 | if( targetCode == protonCode )
|
|---|
| 346 | {
|
|---|
| 347 | G4double cech[] = {0.14, 0.170, 0.180, 0.180, 0.180, 0.170, 0.170, 0.160, 0.155, 0.145,
|
|---|
| 348 | 0.11, 0.082, 0.065, 0.050, 0.041, 0.035, 0.028, 0.024, 0.010, 0.000};
|
|---|
| 349 |
|
|---|
| 350 | G4int iplab = G4int( incidentTotalMomentum*10.);
|
|---|
| 351 | if (iplab > 9) iplab = Imin(19, G4int( incidentTotalMomentum) + 9);
|
|---|
| 352 | if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
|
|---|
| 353 | { // charge exchange pi+ n -> pi0 p
|
|---|
| 354 | pv[0] = AntiNeutron;
|
|---|
| 355 | pv[1] = Neutron;
|
|---|
| 356 | }
|
|---|
| 357 | }
|
|---|
| 358 | return;
|
|---|
| 359 | }
|
|---|
| 360 | else if (availableEnergy <= PionPlus.getMass())
|
|---|
| 361 | return;
|
|---|
| 362 |
|
|---|
| 363 | // inelastic scattering
|
|---|
| 364 |
|
|---|
| 365 | np = 0; nm = 0; nz = 0;
|
|---|
| 366 | G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.90,
|
|---|
| 367 | 0.60, 0.52, 0.47, 0.44, 0.41, 0.39, 0.37, 0.35, 0.34, 0.24,
|
|---|
| 368 | 0.19, 0.15, 0.12, 0.10, 0.09, 0.07, 0.06, 0.05, 0.00};
|
|---|
| 369 | G4int iplab = G4int( incidentTotalMomentum*10.);
|
|---|
| 370 | if ( iplab > 9) iplab = 9 + G4int( incidentTotalMomentum);
|
|---|
| 371 | if ( iplab > 18) iplab = 18 + G4int( incidentTotalMomentum*10.);
|
|---|
| 372 | iplab = Imin(28, iplab);
|
|---|
| 373 |
|
|---|
| 374 | if ( G4UniformRand() > anhl[iplab] )
|
|---|
| 375 | {
|
|---|
| 376 |
|
|---|
| 377 | G4double eab = availableEnergy;
|
|---|
| 378 | G4int ieab = G4int( eab*5.0 );
|
|---|
| 379 |
|
|---|
| 380 | G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
|
|---|
| 381 | if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
|
|---|
| 382 | {
|
|---|
| 383 | // suppress high multiplicity events at low momentum
|
|---|
| 384 | // only one additional pion will be produced
|
|---|
| 385 | G4double w0, wp, wm, wt, ran;
|
|---|
| 386 | if( targetCode == neutronCode ) // target is a neutron
|
|---|
| 387 | {
|
|---|
| 388 | w0 = - sqr(1.+neutb)/(2.*c*c);
|
|---|
| 389 | w0 = std::exp(w0);
|
|---|
| 390 | wm = - sqr(-1.+neutb)/(2.*c*c);
|
|---|
| 391 | wm = std::exp(wm);
|
|---|
| 392 | if( G4UniformRand() < w0/(w0+wm) )
|
|---|
| 393 | { np = 0; nm = 0; nz = 1; }
|
|---|
| 394 | else
|
|---|
| 395 | { np = 0; nm = 1; nz = 0; }
|
|---|
| 396 | }
|
|---|
| 397 | else
|
|---|
| 398 | { // target is a proton
|
|---|
| 399 | w0 = -sqr(1.+protb)/(2.*c*c);
|
|---|
| 400 | w0 = std::exp(w0);
|
|---|
| 401 | wp = w0;
|
|---|
| 402 | wm = -sqr(-1.+protb)/(2.*c*c);
|
|---|
| 403 | wm = std::exp(wm);
|
|---|
| 404 | wt = w0+wp+wm;
|
|---|
| 405 | wp = w0+wp;
|
|---|
| 406 | ran = G4UniformRand();
|
|---|
| 407 | if( ran < w0/wt)
|
|---|
| 408 | { np = 0; nm = 0; nz = 1; }
|
|---|
| 409 | else if( ran < wp/wt)
|
|---|
| 410 | { np = 1; nm = 0; nz = 0; }
|
|---|
| 411 | else
|
|---|
| 412 | { np = 0; nm = 1; nz = 0; }
|
|---|
| 413 | }
|
|---|
| 414 | }
|
|---|
| 415 | else
|
|---|
| 416 | {
|
|---|
| 417 | // number of total particles vs. centre of mass Energy - 2*proton mass
|
|---|
| 418 |
|
|---|
| 419 | G4double aleab = std::log(availableEnergy);
|
|---|
| 420 | G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
|
|---|
| 421 | + aleab*(0.117712+0.0136912*aleab))) - 2.0;
|
|---|
| 422 |
|
|---|
| 423 | // normalization constant for kno-distribution.
|
|---|
| 424 | // calculate first the sum of all constants, check for numerical problems.
|
|---|
| 425 | G4double test, dum, anpn = 0.0;
|
|---|
| 426 |
|
|---|
| 427 | for (nt=1; nt<=numSec; nt++) {
|
|---|
| 428 | test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
|
|---|
| 429 | dum = pi*nt/(2.0*n*n);
|
|---|
| 430 | if (std::fabs(dum) < 1.0) {
|
|---|
| 431 | if( test >= 1.0e-10 )anpn += dum*test;
|
|---|
| 432 | } else {
|
|---|
| 433 | anpn += dum*test;
|
|---|
| 434 | }
|
|---|
| 435 | }
|
|---|
| 436 |
|
|---|
| 437 | G4double ran = G4UniformRand();
|
|---|
| 438 | G4double excs = 0.0;
|
|---|
| 439 | if( targetCode == protonCode )
|
|---|
| 440 | {
|
|---|
| 441 | counter = -1;
|
|---|
| 442 | for( np=0; np<numSec/3; np++ )
|
|---|
| 443 | {
|
|---|
| 444 | for( nm=Imax(0,np-1); nm<=(np+1); nm++ )
|
|---|
| 445 | {
|
|---|
| 446 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 447 | {
|
|---|
| 448 | if( ++counter < numMul )
|
|---|
| 449 | {
|
|---|
| 450 | nt = np+nm+nz;
|
|---|
| 451 | if ( (nt>0) && (nt<=numSec) ) {
|
|---|
| 452 | test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
|
|---|
| 453 | dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
|
|---|
| 454 | if (std::fabs(dum) < 1.0) {
|
|---|
| 455 | if( test >= 1.0e-10 )excs += dum*test;
|
|---|
| 456 | } else {
|
|---|
| 457 | excs += dum*test;
|
|---|
| 458 | }
|
|---|
| 459 |
|
|---|
| 460 | if (ran < excs) goto outOfLoop; //----------------------->
|
|---|
| 461 | }
|
|---|
| 462 | }
|
|---|
| 463 | }
|
|---|
| 464 | }
|
|---|
| 465 | }
|
|---|
| 466 |
|
|---|
| 467 | // 3 previous loops continued to the end
|
|---|
| 468 | inElastic = false; // quasi-elastic scattering
|
|---|
| 469 | return;
|
|---|
| 470 | }
|
|---|
| 471 | else
|
|---|
| 472 | { // target must be a neutron
|
|---|
| 473 | counter = -1;
|
|---|
| 474 | for( np=0; np<numSec/3; np++ )
|
|---|
| 475 | {
|
|---|
| 476 | for( nm=np; nm<=(np+2); nm++ )
|
|---|
| 477 | {
|
|---|
| 478 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 479 | {
|
|---|
| 480 | if( ++counter < numMul )
|
|---|
| 481 | {
|
|---|
| 482 | nt = np+nm+nz;
|
|---|
| 483 | if ( (nt>=1) && (nt<=numSec) ) {
|
|---|
| 484 | test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
|
|---|
| 485 | dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
|
|---|
| 486 | if (std::fabs(dum) < 1.0) {
|
|---|
| 487 | if( test >= 1.0e-10 )excs += dum*test;
|
|---|
| 488 | } else {
|
|---|
| 489 | excs += dum*test;
|
|---|
| 490 | }
|
|---|
| 491 |
|
|---|
| 492 | if (ran < excs) goto outOfLoop; // -------------------------->
|
|---|
| 493 | }
|
|---|
| 494 | }
|
|---|
| 495 | }
|
|---|
| 496 | }
|
|---|
| 497 | }
|
|---|
| 498 | // 3 previous loops continued to the end
|
|---|
| 499 | inElastic = false; // quasi-elastic scattering.
|
|---|
| 500 | return;
|
|---|
| 501 | }
|
|---|
| 502 | }
|
|---|
| 503 | outOfLoop: // <------------------------------------------------------------------------
|
|---|
| 504 |
|
|---|
| 505 | if( targetCode == neutronCode)
|
|---|
| 506 | {
|
|---|
| 507 | if( np == nm)
|
|---|
| 508 | {
|
|---|
| 509 | }
|
|---|
| 510 | else if (np == (nm-1))
|
|---|
| 511 | {
|
|---|
| 512 | if( G4UniformRand() < 0.5)
|
|---|
| 513 | {
|
|---|
| 514 | pv[1] = Proton;
|
|---|
| 515 | }
|
|---|
| 516 | else
|
|---|
| 517 | {
|
|---|
| 518 | pv[0] = AntiNeutron;
|
|---|
| 519 | }
|
|---|
| 520 | }
|
|---|
| 521 | else
|
|---|
| 522 | {
|
|---|
| 523 | pv[0] = AntiNeutron;
|
|---|
| 524 | pv[1] = Proton;
|
|---|
| 525 | }
|
|---|
| 526 | }
|
|---|
| 527 | else
|
|---|
| 528 | {
|
|---|
| 529 | if( np == nm)
|
|---|
| 530 | {
|
|---|
| 531 | if( G4UniformRand() < 0.25)
|
|---|
| 532 | {
|
|---|
| 533 | pv[0] = AntiNeutron;
|
|---|
| 534 | pv[1] = Neutron;
|
|---|
| 535 | }
|
|---|
| 536 | else
|
|---|
| 537 | {
|
|---|
| 538 | }
|
|---|
| 539 | }
|
|---|
| 540 | else if ( np == (1+nm))
|
|---|
| 541 | {
|
|---|
| 542 | pv[1] = Neutron;
|
|---|
| 543 | }
|
|---|
| 544 | else
|
|---|
| 545 | {
|
|---|
| 546 | pv[0] = AntiNeutron;
|
|---|
| 547 | }
|
|---|
| 548 | }
|
|---|
| 549 |
|
|---|
| 550 | }
|
|---|
| 551 | else // annihilation
|
|---|
| 552 | {
|
|---|
| 553 | if ( availableEnergy > 2. * PionPlus.getMass() )
|
|---|
| 554 | {
|
|---|
| 555 |
|
|---|
| 556 | G4double aleab = std::log(availableEnergy);
|
|---|
| 557 | G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
|
|---|
| 558 | + aleab*(0.117712+0.0136912*aleab))) - 2.0;
|
|---|
| 559 |
|
|---|
| 560 | // normalization constant for kno-distribution.
|
|---|
| 561 | // calculate first the sum of all constants, check for numerical problems.
|
|---|
| 562 | G4double test, dum, anpn = 0.0;
|
|---|
| 563 |
|
|---|
| 564 | for (nt=2; nt<=numSec; nt++) {
|
|---|
| 565 | test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
|
|---|
| 566 | dum = pi*nt/(2.0*n*n);
|
|---|
| 567 | if (std::fabs(dum) < 1.0) {
|
|---|
| 568 | if( test >= 1.0e-10 )anpn += dum*test;
|
|---|
| 569 | } else {
|
|---|
| 570 | anpn += dum*test;
|
|---|
| 571 | }
|
|---|
| 572 | }
|
|---|
| 573 |
|
|---|
| 574 | G4double ran = G4UniformRand();
|
|---|
| 575 | G4double excs = 0.0;
|
|---|
| 576 | if( targetCode == protonCode )
|
|---|
| 577 | {
|
|---|
| 578 | counter = -1;
|
|---|
| 579 | for( np=1; np<numSec/3; np++ )
|
|---|
| 580 | {
|
|---|
| 581 | nm = np;
|
|---|
| 582 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 583 | {
|
|---|
| 584 | if( ++counter < numMulAn )
|
|---|
| 585 | {
|
|---|
| 586 | nt = np+nm+nz;
|
|---|
| 587 | if ( (nt>0) && (nt<=numSec) ) {
|
|---|
| 588 | test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
|
|---|
| 589 | dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
|
|---|
| 590 | if (std::fabs(dum) < 1.0) {
|
|---|
| 591 | if( test >= 1.0e-10 )excs += dum*test;
|
|---|
| 592 | } else {
|
|---|
| 593 | excs += dum*test;
|
|---|
| 594 | }
|
|---|
| 595 |
|
|---|
| 596 | if (ran < excs) goto outOfLoopAn; //----------------------->
|
|---|
| 597 | }
|
|---|
| 598 | }
|
|---|
| 599 | }
|
|---|
| 600 | }
|
|---|
| 601 | // 3 previous loops continued to the end
|
|---|
| 602 | inElastic = false; // quasi-elastic scattering
|
|---|
| 603 | return;
|
|---|
| 604 | }
|
|---|
| 605 | else
|
|---|
| 606 | { // target must be a neutron
|
|---|
| 607 | counter = -1;
|
|---|
| 608 | for( np=1; np<numSec/3; np++ )
|
|---|
| 609 | {
|
|---|
| 610 | nm = np+1;
|
|---|
| 611 | for( nz=0; nz<numSec/3; nz++ )
|
|---|
| 612 | {
|
|---|
| 613 | if( ++counter < numMulAn )
|
|---|
| 614 | {
|
|---|
| 615 | nt = np+nm+nz;
|
|---|
| 616 | if ( (nt>=1) && (nt<=numSec) ) {
|
|---|
| 617 | test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
|
|---|
| 618 | dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
|
|---|
| 619 | if (std::fabs(dum) < 1.0) {
|
|---|
| 620 | if( test >= 1.0e-10 )excs += dum*test;
|
|---|
| 621 | } else {
|
|---|
| 622 | excs += dum*test;
|
|---|
| 623 | }
|
|---|
| 624 |
|
|---|
| 625 | if (ran < excs) goto outOfLoopAn; // -------------------------->
|
|---|
| 626 | }
|
|---|
| 627 | }
|
|---|
| 628 | }
|
|---|
| 629 | }
|
|---|
| 630 | inElastic = false; // quasi-elastic scattering.
|
|---|
| 631 | return;
|
|---|
| 632 | }
|
|---|
| 633 | outOfLoopAn: // <------------------------------------------------------------------
|
|---|
| 634 | vecLen = 0;
|
|---|
| 635 | }
|
|---|
| 636 | }
|
|---|
| 637 |
|
|---|
| 638 | nt = np + nm + nz;
|
|---|
| 639 | while ( nt > 0)
|
|---|
| 640 | {
|
|---|
| 641 | G4double ran = G4UniformRand();
|
|---|
| 642 | if ( ran < (G4double)np/nt)
|
|---|
| 643 | {
|
|---|
| 644 | if( np > 0 )
|
|---|
| 645 | { pv[vecLen++] = PionPlus;
|
|---|
| 646 | np--;
|
|---|
| 647 | }
|
|---|
| 648 | }
|
|---|
| 649 | else if ( ran < (G4double)(np+nm)/nt)
|
|---|
| 650 | {
|
|---|
| 651 | if( nm > 0 )
|
|---|
| 652 | {
|
|---|
| 653 | pv[vecLen++] = PionMinus;
|
|---|
| 654 | nm--;
|
|---|
| 655 | }
|
|---|
| 656 | }
|
|---|
| 657 | else
|
|---|
| 658 | {
|
|---|
| 659 | if( nz > 0 )
|
|---|
| 660 | {
|
|---|
| 661 | pv[vecLen++] = PionZero;
|
|---|
| 662 | nz--;
|
|---|
| 663 | }
|
|---|
| 664 | }
|
|---|
| 665 | nt = np + nm + nz;
|
|---|
| 666 | }
|
|---|
| 667 | if (verboseLevel > 1)
|
|---|
| 668 | {
|
|---|
| 669 | G4cout << "Particles produced: " ;
|
|---|
| 670 | G4cout << pv[0].getName() << " " ;
|
|---|
| 671 | G4cout << pv[1].getName() << " " ;
|
|---|
| 672 | for (i=2; i < vecLen; i++)
|
|---|
| 673 | {
|
|---|
| 674 | G4cout << pv[i].getName() << " " ;
|
|---|
| 675 | }
|
|---|
| 676 | G4cout << G4endl;
|
|---|
| 677 | }
|
|---|
| 678 | return;
|
|---|
| 679 | }
|
|---|
| 680 |
|
|---|
| 681 |
|
|---|
| 682 |
|
|---|
| 683 |
|
|---|
| 684 |
|
|---|
| 685 |
|
|---|
| 686 |
|
|---|
| 687 |
|
|---|
| 688 |
|
|---|