source: trunk/source/processes/hadronic/models/high_energy/src/G4HEKaonMinusInelastic.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 22.5 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//
27// $Id: G4HEKaonMinusInelastic.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
39// like nuclear reactions, nuclear fission without any cascading and all
40// processes for 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 "G4HEKaonMinusInelastic.hh"
46
47G4HadFinalState * G4HEKaonMinusInelastic::
48ApplyYourself( 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 A = targetNucleus.GetN();
54 const G4double Z = targetNucleus.GetZ();
55 G4HEVector incidentParticle(aParticle);
56
57 G4double atomicNumber = Z;
58 G4double atomicWeight = A;
59
60 G4int incidentCode = incidentParticle.getCode();
61 G4double incidentMass = incidentParticle.getMass();
62 G4double incidentTotalEnergy = incidentParticle.getEnergy();
63 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
64 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
65
66 if(incidentKineticEnergy < 1.)
67 {
68 G4cout << "GHEKaonMinusInelastic: incident energy < 1 GeV" << G4endl;
69 }
70 if(verboseLevel > 1)
71 {
72 G4cout << "G4HEKaonMinusInelastic::ApplyYourself" << G4endl;
73 G4cout << "incident particle " << incidentParticle.getName()
74 << "mass " << incidentMass
75 << "kinetic energy " << incidentKineticEnergy
76 << G4endl;
77 G4cout << "target material with (A,Z) = ("
78 << atomicWeight << "," << atomicNumber << ")" << G4endl;
79 }
80
81 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
82 atomicWeight, atomicNumber);
83 if(verboseLevel > 1)
84 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
85
86 incidentKineticEnergy -= inelasticity;
87
88 G4double excitationEnergyGNP = 0.;
89 G4double excitationEnergyDTA = 0.;
90
91 G4double excitation = NuclearExcitation(incidentKineticEnergy,
92 atomicWeight, atomicNumber,
93 excitationEnergyGNP,
94 excitationEnergyDTA);
95 if(verboseLevel > 1)
96 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
97 << excitationEnergyDTA << G4endl;
98
99
100 incidentKineticEnergy -= excitation;
101 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
102 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
103 *(incidentTotalEnergy+incidentMass));
104
105
106 G4HEVector targetParticle;
107 if(G4UniformRand() < atomicNumber/atomicWeight)
108 {
109 targetParticle.setDefinition("Proton");
110 }
111 else
112 {
113 targetParticle.setDefinition("Neutron");
114 }
115
116 G4double targetMass = targetParticle.getMass();
117 G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass
118 + 2.0*targetMass*incidentTotalEnergy);
119 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
120
121 // this was the meaning of inElastic in the
122 // original Gheisha stand-alone version.
123// G4bool inElastic = InElasticCrossSectionInFirstInt
124// (availableEnergy, incidentCode, incidentTotalMomentum);
125 // by unknown reasons, it has been replaced
126 // to the following code in Geant???
127 G4bool inElastic = true;
128// if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;
129
130
131 vecLength = 0;
132
133 if(verboseLevel > 1)
134 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
135 << incidentCode << G4endl;
136
137 G4bool successful = false;
138
139 if(inElastic || (!inElastic && atomicWeight < 1.5))
140 {
141 FirstIntInCasKaonMinus(inElastic, availableEnergy, pv, vecLength,
142 incidentParticle, targetParticle);
143
144 if(verboseLevel > 1)
145 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
146
147
148 if ((vecLength > 0) && (availableEnergy > 1.))
149 StrangeParticlePairProduction( availableEnergy, centerOfMassEnergy,
150 pv, vecLength,
151 incidentParticle, targetParticle);
152 HighEnergyCascading( successful, pv, vecLength,
153 excitationEnergyGNP, excitationEnergyDTA,
154 incidentParticle, targetParticle,
155 atomicWeight, atomicNumber);
156 if (!successful)
157 HighEnergyClusterProduction( successful, pv, vecLength,
158 excitationEnergyGNP, excitationEnergyDTA,
159 incidentParticle, targetParticle,
160 atomicWeight, atomicNumber);
161 if (!successful)
162 MediumEnergyCascading( successful, pv, vecLength,
163 excitationEnergyGNP, excitationEnergyDTA,
164 incidentParticle, targetParticle,
165 atomicWeight, atomicNumber);
166
167 if (!successful)
168 MediumEnergyClusterProduction( successful, pv, vecLength,
169 excitationEnergyGNP, excitationEnergyDTA,
170 incidentParticle, targetParticle,
171 atomicWeight, atomicNumber);
172 if (!successful)
173 QuasiElasticScattering( successful, pv, vecLength,
174 excitationEnergyGNP, excitationEnergyDTA,
175 incidentParticle, targetParticle,
176 atomicWeight, atomicNumber);
177 }
178 if (!successful)
179 {
180 ElasticScattering( successful, pv, vecLength,
181 incidentParticle,
182 atomicWeight, atomicNumber);
183 }
184
185 if (!successful)
186 {
187 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" << G4endl;
188 }
189 FillParticleChange(pv, vecLength);
190 delete [] pv;
191 theParticleChange.SetStatusChange(stopAndKill);
192 return & theParticleChange;
193 }
194
195void
196 G4HEKaonMinusInelastic::FirstIntInCasKaonMinus( G4bool &inElastic,
197 const G4double availableEnergy,
198 G4HEVector pv[],
199 G4int &vecLen,
200 G4HEVector incidentParticle,
201 G4HEVector targetParticle)
202
203// Kaon- undergoes interaction with nucleon within a nucleus. Check if it is
204// energetically possible to produce pions/kaons. In not, assume nuclear excitation
205// occurs and input particle is degraded in energy. No other particles are produced.
206// If reaction is possible, find the correct number of pions/protons/neutrons
207// produced using an interpolation to multiplicity data. Replace some pions or
208// protons/neutrons by kaons or strange baryons according to the average
209// multiplicity per inelastic reaction.
210
211 {
212 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp
213 static const G4double expxl = -expxu; // lower bound for arg. of exp
214
215 static const G4double protb = 0.7;
216 static const G4double neutb = 0.7;
217 static const G4double c = 1.25;
218
219 static const G4int numMul = 1200;
220 static const G4int numSec = 60;
221
222 G4int neutronCode = Neutron.getCode();
223 G4int protonCode = Proton.getCode();
224 G4int kaonMinusCode = KaonMinus.getCode();
225 G4int kaonZeroCode = KaonZero.getCode();
226 G4int antiKaonZeroCode = AntiKaonZero.getCode();
227
228 G4int targetCode = targetParticle.getCode();
229// G4double incidentMass = incidentParticle.getMass();
230// G4double incidentEnergy = incidentParticle.getEnergy();
231 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
232
233 static G4bool first = true;
234 static G4double protmul[numMul], protnorm[numSec]; // proton constants
235 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
236
237// misc. local variables
238// np = number of pi+, nm = number of pi-, nz = number of pi0
239
240 G4int i, counter, nt, np, nm, nz;
241
242 if( first )
243 { // compute normalization constants, this will only be done once
244 first = false;
245 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
246 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
247 counter = -1;
248 for( np=0; np<(numSec/3); np++ )
249 {
250 for( nm=Imax(0,np-1); nm<=np+1; nm++ )
251 {
252 for( nz=0; nz<numSec/3; nz++ )
253 {
254 if( ++counter < numMul )
255 {
256 nt = np+nm+nz;
257 if( (nt>0) && (nt<=numSec) )
258 {
259 protmul[counter] =
260 pmltpc(np,nm,nz,nt,protb,c) ;
261 protnorm[nt-1] += protmul[counter];
262 }
263 }
264 }
265 }
266 }
267 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
268 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
269 counter = -1;
270 for( np=0; np<numSec/3; np++ )
271 {
272 for( nm=np; nm<=(np+2); nm++ )
273 {
274 for( nz=0; nz<numSec/3; nz++ )
275 {
276 if( ++counter < numMul )
277 {
278 nt = np+nm+nz;
279 if( (nt>0) && (nt<=numSec) )
280 {
281 neutmul[counter] =
282 pmltpc(np,nm,nz,nt,neutb,c);
283 neutnorm[nt-1] += neutmul[counter];
284 }
285 }
286 }
287 }
288 }
289 for( i=0; i<numSec; i++ )
290 {
291 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
292 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
293 }
294 } // end of initialization
295
296
297 // initialize the first two places
298 // the same as beam and target
299 pv[0] = incidentParticle;
300 pv[1] = targetParticle;
301 vecLen = 2;
302
303 if (!inElastic || (availableEnergy <= PionPlus.getMass()))
304 return;
305
306
307// inelastic scattering
308
309 np = 0, nm = 0, nz = 0;
310 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
311 G4int iplab = G4int( incidentTotalMomentum*5.);
312 if( (iplab < 10) && (G4UniformRand() < cech[iplab]))
313 {
314 G4int iplab = Imin(19, G4int( incidentTotalMomentum*5.));
315 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
316 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
317 if( G4UniformRand() < cnk0[iplab] )
318 {
319 if( targetCode == protonCode )
320 {
321 pv[0] = AntiKaonZero;
322 pv[1] = Neutron;
323 return;
324 }
325 else
326 {
327 return;
328 }
329 }
330 G4double ran = G4UniformRand();
331 if( targetCode == protonCode ) // target is a proton
332 {
333 if( ran < 0.25 )
334 {
335 pv[0] = PionMinus;
336 pv[1] = SigmaPlus;
337 }
338 else if (ran < 0.50)
339 {
340 pv[0] = PionZero;
341 pv[1] = SigmaZero;
342 }
343 else if (ran < 0.75)
344 {
345 pv[0] = PionPlus;
346 pv[1] = SigmaMinus;
347 }
348 else
349 {
350 pv[0] = PionZero;
351 pv[1] = Lambda;
352 }
353 }
354 else
355 { // target is a neutron
356 if( ran < 0.25 )
357 {
358 }
359 else if (ran < 0.50)
360 {
361 pv[0] = PionMinus;
362 pv[1] = SigmaZero;
363 }
364 else if (ran < 0.75)
365 {
366 pv[0] = PionZero;
367 pv[1] = SigmaMinus;
368 }
369 else
370 {
371 pv[0] = PionMinus;
372 pv[1] = Lambda;
373 }
374 }
375 return;
376 }
377 else
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( Amin( expxu, Amax( 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=Imax(0,np-1); nm<=np+1; nm++ )
407 {
408 for (nz=0; nz<numSec/3; nz++) {
409 if (++counter < numMul) {
410 nt = np+nm+nz;
411 if ( (nt>0) && (nt<=numSec) ) {
412 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
413 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
414 if (std::fabs(dum) < 1.0) {
415 if( test >= 1.0e-10 )excs += dum*test;
416 } else {
417 excs += dum*test;
418 }
419
420 if (ran < excs) goto outOfLoop; //----------------------->
421 }
422 }
423 }
424 }
425 }
426
427 // 3 previous loops continued to the end
428 inElastic = false; // quasi-elastic scattering
429 return;
430 }
431 else
432 { // target must be a neutron
433 counter = -1;
434 for( np=0; np<numSec/3; np++ )
435 {
436 for( nm=np; nm<=(np+2); nm++ )
437 {
438 for (nz=0; nz<numSec/3; nz++) {
439 if (++counter < numMul) {
440 nt = np+nm+nz;
441 if ( (nt>=1) && (nt<=numSec) ) {
442 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
443 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
444 if (std::fabs(dum) < 1.0) {
445 if( test >= 1.0e-10 )excs += dum*test;
446 } else {
447 excs += dum*test;
448 }
449 if (ran < excs) goto outOfLoop; // -------------------------->
450 }
451 }
452 }
453 }
454 }
455 // 3 previous loops continued to the end
456 inElastic = false; // quasi-elastic scattering.
457 return;
458 }
459 }
460 outOfLoop: // <---------------------------------------------
461
462 if( targetCode == protonCode)
463 {
464 if( np == (1+nm))
465 {
466 pv[1] = Neutron;
467 }
468 else if (np == nm)
469 {
470 if( G4UniformRand() < 0.75)
471 {
472 }
473 else
474 {
475 pv[0] = AntiKaonZero;
476 pv[1] = Neutron;
477 }
478 }
479 else
480 {
481 pv[0] = AntiKaonZero;
482 }
483 }
484 else
485 {
486 if( np == (nm-1))
487 {
488 if( G4UniformRand() < 0.5)
489 {
490 pv[1] = Proton;
491 }
492 else
493 {
494 pv[0] = AntiKaonZero;
495 }
496 }
497 else if ( np == nm)
498 {
499 }
500 else
501 {
502 pv[0] = AntiKaonZero;
503 pv[1] = Proton;
504 }
505 }
506
507 if( G4UniformRand() < 0.5 )
508 {
509 if( ( (pv[0].getCode() == kaonMinusCode)
510 && (pv[1].getCode() == neutronCode) )
511 || ( (pv[0].getCode() == kaonZeroCode)
512 && (pv[1].getCode() == protonCode) )
513 || ( (pv[0].getCode() == antiKaonZeroCode)
514 && (pv[1].getCode() == protonCode) ) )
515 {
516 G4double ran = G4UniformRand();
517 if( pv[1].getCode() == protonCode)
518 {
519 if(ran < 0.68)
520 {
521 pv[0] = PionPlus;
522 pv[1] = Lambda;
523 }
524 else if (ran < 0.84)
525 {
526 pv[0] = PionZero;
527 pv[1] = SigmaPlus;
528 }
529 else
530 {
531 pv[0] = PionPlus;
532 pv[1] = SigmaZero;
533 }
534 }
535 else
536 {
537 if(ran < 0.68)
538 {
539 pv[0] = PionMinus;
540 pv[1] = Lambda;
541 }
542 else if (ran < 0.84)
543 {
544 pv[0] = PionMinus;
545 pv[1] = SigmaZero;
546 }
547 else
548 {
549 pv[0] = PionZero;
550 pv[1] = SigmaMinus;
551 }
552 }
553 }
554 else
555 {
556 G4double ran = G4UniformRand();
557 if (ran < 0.67)
558 {
559 pv[0] = PionZero;
560 pv[1] = Lambda;
561 }
562 else if (ran < 0.78)
563 {
564 pv[0] = PionMinus;
565 pv[1] = SigmaPlus;
566 }
567 else if (ran < 0.89)
568 {
569 pv[0] = PionZero;
570 pv[1] = SigmaZero;
571 }
572 else
573 {
574 pv[0] = PionPlus;
575 pv[1] = SigmaMinus;
576 }
577 }
578 }
579
580
581 nt = np + nm + nz;
582 while ( nt > 0)
583 {
584 G4double ran = G4UniformRand();
585 if ( ran < (G4double)np/nt)
586 {
587 if( np > 0 )
588 { pv[vecLen++] = PionPlus;
589 np--;
590 }
591 }
592 else if ( ran < (G4double)(np+nm)/nt)
593 {
594 if( nm > 0 )
595 {
596 pv[vecLen++] = PionMinus;
597 nm--;
598 }
599 }
600 else
601 {
602 if( nz > 0 )
603 {
604 pv[vecLen++] = PionZero;
605 nz--;
606 }
607 }
608 nt = np + nm + nz;
609 }
610 if (verboseLevel > 1)
611 {
612 G4cout << "Particles produced: " ;
613 G4cout << pv[0].getName() << " " ;
614 G4cout << pv[1].getName() << " " ;
615 for (i=2; i < vecLen; i++)
616 {
617 G4cout << pv[i].getName() << " " ;
618 }
619 G4cout << G4endl;
620 }
621 return;
622 }
623
624
625
626
627
628
629
630
631
Note: See TracBrowser for help on using the repository browser.