source: trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiSigmaMinusInelastic.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

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