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