source: trunk/source/processes/hadronic/models/high_energy/src/G4HEPionMinusInelastic.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: 19.6 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: G4HEPionMinusInelastic.cc,v 1.15 2008/03/17 20:49:17 dennis Exp $
28// GEANT4 tag $Name: geant4-09-02 $
29
30// 11-OCT-2007 F.W. Jones: fixed incorrect Imax (should be Imin) in
31// sampling of charge exchange.
32
33
34#include "globals.hh"
35#include "G4ios.hh"
36
37//
38// G4 Process: Gheisha High Energy Collision model.
39// This includes the high energy cascading model, the two-body-resonance model
40// and the low energy two-body model. Not included are the low energy stuff like
41// nuclear reactions, nuclear fission without any cascading and all processes for
42// particles at rest.
43// First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
44// H. Fesefeldt, RWTH-Aachen, 23-October-1996
45// Last modified: 29-July-1998
46
47#include "G4HEPionMinusInelastic.hh"
48
49G4HadFinalState * G4HEPionMinusInelastic::
50ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &targetNucleus )
51 {
52 G4HEVector * pv = new G4HEVector[MAXPART];
53 const G4HadProjectile *aParticle = &aTrack;
54// G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
55 const G4double A = targetNucleus.GetN();
56 const G4double Z = targetNucleus.GetZ();
57 G4HEVector incidentParticle(aParticle);
58
59 G4double atomicNumber = Z;
60 G4double atomicWeight = A;
61
62 G4int incidentCode = incidentParticle.getCode();
63 G4double incidentMass = incidentParticle.getMass();
64 G4double incidentTotalEnergy = incidentParticle.getEnergy();
65 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
66 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
67
68 if(incidentKineticEnergy < 1.)
69 {
70 G4cout << "GHEPionMinusInelastic: incident energy < 1 GeV" << G4endl ;
71 }
72 if(verboseLevel > 1)
73 {
74 G4cout << "G4HEPionMinusInelastic::ApplyYourself" << G4endl;
75 G4cout << "incident particle " << incidentParticle.getName()
76 << "mass " << incidentMass
77 << "kinetic energy " << incidentKineticEnergy
78 << G4endl;
79 G4cout << "target material with (A,Z) = ("
80 << atomicWeight << "," << atomicNumber << ")" << G4endl;
81 }
82
83 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
84 atomicWeight, atomicNumber);
85 if(verboseLevel > 1)
86 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
87
88 incidentKineticEnergy -= inelasticity;
89
90 G4double excitationEnergyGNP = 0.;
91 G4double excitationEnergyDTA = 0.;
92
93 G4double excitation = NuclearExcitation(incidentKineticEnergy,
94 atomicWeight, atomicNumber,
95 excitationEnergyGNP,
96 excitationEnergyDTA);
97 if(verboseLevel > 1)
98 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
99 << excitationEnergyDTA << G4endl;
100
101
102 incidentKineticEnergy -= excitation;
103 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
104 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
105 *(incidentTotalEnergy+incidentMass));
106
107 G4HEVector targetParticle;
108
109 if(G4UniformRand() < atomicNumber/atomicWeight)
110 {
111 targetParticle.setDefinition("Proton");
112 }
113 else
114 {
115 targetParticle.setDefinition("Neutron");
116 }
117
118 G4double targetMass = targetParticle.getMass();
119 G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass
120 + 2.0*targetMass*incidentTotalEnergy);
121 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
122
123 // this was the meaning of inElastic in the
124 // original Gheisha stand-alone version.
125// G4bool inElastic = InElasticCrossSectionInFirstInt
126// (availableEnergy, incidentCode, incidentTotalMomentum);
127 // by unknown reasons, it has been replaced
128 // to the following code in Geant???
129 G4bool inElastic = true;
130// if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;
131
132 vecLength = 0;
133
134 if(verboseLevel > 1)
135 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
136 << incidentCode << G4endl;
137
138 G4bool successful = false;
139
140 if(inElastic || (!inElastic && atomicWeight < 1.5))
141 {
142 FirstIntInCasPionMinus(inElastic, availableEnergy, pv, vecLength,
143 incidentParticle, targetParticle);
144
145 if(verboseLevel > 1)
146 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
147
148
149 if ((vecLength > 0) && (availableEnergy > 1.))
150 StrangeParticlePairProduction( availableEnergy, centerOfMassEnergy,
151 pv, vecLength,
152 incidentParticle, targetParticle);
153 HighEnergyCascading( successful, pv, vecLength,
154 excitationEnergyGNP, excitationEnergyDTA,
155 incidentParticle, targetParticle,
156 atomicWeight, atomicNumber);
157 if (!successful)
158 HighEnergyClusterProduction( successful, pv, vecLength,
159 excitationEnergyGNP, excitationEnergyDTA,
160 incidentParticle, targetParticle,
161 atomicWeight, atomicNumber);
162 if (!successful)
163 MediumEnergyCascading( successful, pv, vecLength,
164 excitationEnergyGNP, excitationEnergyDTA,
165 incidentParticle, targetParticle,
166 atomicWeight, atomicNumber);
167
168 if (!successful)
169 MediumEnergyClusterProduction( successful, pv, vecLength,
170 excitationEnergyGNP, excitationEnergyDTA,
171 incidentParticle, targetParticle,
172 atomicWeight, atomicNumber);
173 if (!successful)
174 QuasiElasticScattering( successful, pv, vecLength,
175 excitationEnergyGNP, excitationEnergyDTA,
176 incidentParticle, targetParticle,
177 atomicWeight, atomicNumber);
178 }
179 if (!successful)
180 {
181 ElasticScattering( successful, pv, vecLength,
182 incidentParticle,
183 atomicWeight, atomicNumber);
184 }
185
186 if (!successful)
187 {
188 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" << G4endl;
189 }
190 FillParticleChange(pv, vecLength);
191 delete [] pv;
192 theParticleChange.SetStatusChange(stopAndKill);
193 return & theParticleChange;
194 }
195
196void
197G4HEPionMinusInelastic::FirstIntInCasPionMinus( G4bool &inElastic,
198 const G4double availableEnergy,
199 G4HEVector pv[],
200 G4int &vecLen,
201 G4HEVector incidentParticle,
202 G4HEVector targetParticle)
203
204// Pion- undergoes interaction with nucleon within a nucleus. Check if it is
205// energetically possible to produce pions/kaons. In not, assume nuclear excitation
206// occurs and input particle is degraded in energy. No other particles are produced.
207// If reaction is possible, find the correct number of pions/protons/neutrons
208// produced using an interpolation to multiplicity data. Replace some pions or
209// protons/neutrons by kaons or strange baryons according to the average
210// multiplicity per inelastic reaction.
211
212 {
213 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp
214 static const G4double expxl = -expxu; // lower bound for arg. of exp
215
216 static const G4double protb = 0.7;
217 static const G4double neutb = 0.7;
218 static const G4double c = 1.25;
219
220 static const G4int numMul = 1200;
221 static const G4int numSec = 60;
222
223// G4int neutronCode = Neutron.getCode();
224 G4int protonCode = Proton.getCode();
225
226 G4int targetCode = targetParticle.getCode();
227// G4double incidentMass = incidentParticle.getMass();
228// G4double incidentEnergy = incidentParticle.getEnergy();
229 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
230
231 static G4bool first = true;
232 static G4double protmul[numMul], protnorm[numSec]; // proton constants
233 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
234
235// misc. local variables
236// np = number of pi+, nm = number of pi-, nz = number of pi0
237
238 G4int i, counter, nt, np, nm, nz;
239
240 if( first )
241 { // compute normalization constants, this will only be done once
242 first = false;
243 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
244 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
245 counter = -1;
246 for( np=0; np<(numSec/3); np++ )
247 {
248 for( nm=Imax(0,np-1); nm<=(np+1); nm++ )
249 {
250 for( nz=0; nz<numSec/3; nz++ )
251 {
252 if( ++counter < numMul )
253 {
254 nt = np+nm+nz;
255 if( (nt>0) && (nt<=numSec) )
256 {
257 protmul[counter] =
258 pmltpc(np,nm,nz,nt,protb,c) ;
259 protnorm[nt-1] += protmul[counter];
260 }
261 }
262 }
263 }
264 }
265 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
266 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
267 counter = -1;
268 for( np=0; np<numSec/3; np++ )
269 {
270 for( nm=np; nm<=(np+2); nm++ )
271 {
272 for( nz=0; nz<numSec/3; nz++ )
273 {
274 if( ++counter < numMul )
275 {
276 nt = np+nm+nz;
277 if( (nt>0) && (nt<=numSec) )
278 {
279 neutmul[counter] =
280 pmltpc(np,nm,nz,nt,neutb,c);
281 neutnorm[nt-1] += neutmul[counter];
282 }
283 }
284 }
285 }
286 }
287 for( i=0; i<numSec; i++ )
288 {
289 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
290 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
291 }
292 } // end of initialization
293
294
295 // initialize the first two places
296 // the same as beam and target
297 pv[0] = incidentParticle;
298 pv[1] = targetParticle;
299 vecLen = 2;
300
301 if (!inElastic || (availableEnergy <= PionPlus.getMass()))
302 return;
303
304
305 // inelastic scattering
306
307 np = 0, nm = 0, nz = 0;
308 G4double eab = availableEnergy;
309 G4int ieab = G4int( eab*5.0 );
310
311 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
312 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
313 {
314 // suppress high multiplicity events at low momentum
315 // only one additional pion will be produced
316 G4double cech[] = {1., 0.95, 0.79, 0.32, 0.19, 0.16, 0.14, 0.12, 0.10, 0.08};
317 G4int iplab = Imin(9, G4int( incidentTotalMomentum*5.));
318 if( G4UniformRand() < cech[iplab] )
319 {
320 if( targetCode == protonCode )
321 {
322 pv[0] = PionZero;
323 pv[1] = Neutron;
324 return;
325 }
326 else
327 {
328 return;
329 }
330 }
331
332 G4double w0, wp, wm, wt, ran;
333 if( targetCode == protonCode ) // target is a proton
334 {
335 w0 = - sqr(1.+protb)/(2.*c*c);
336 wp = w0 = std::exp(w0);
337 wm = - sqr(-1.+protb)/(2.*c*c);
338 wm = std::exp(wm);
339 wp *= 10.;
340 wt = w0+wp+wm;
341 wp = w0+wp;
342 ran = G4UniformRand();
343 if( ran < w0/wt )
344 { np = 0; nm = 0; nz = 1; }
345 else if ( ran < wp/wt )
346 { np = 1; nm = 0; nz = 0; }
347 else
348 { np = 0; nm = 1; nz = 0; }
349 }
350 else
351 { // target is a neutron
352 w0 = -sqr(1.+neutb)/(2.*c*c);
353 wm = -sqr(-1.+neutb)/(2.*c*c);
354 w0 = std::exp(w0);
355 wm = std::exp(wm);
356 if( G4UniformRand() < w0/(w0+wm) )
357 { np = 0; nm = 0; nz = 1; }
358 else
359 { np = 0; nm = 1; nz = 0; }
360 }
361 }
362 else
363 {
364 // number of total particles vs. centre of mass Energy - 2*proton mass
365
366 G4double aleab = std::log(availableEnergy);
367 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
368 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
369
370 // normalization constant for kno-distribution.
371 // calculate first the sum of all constants, check for numerical problems.
372 G4double test, dum, anpn = 0.0;
373
374 for (nt=1; nt<=numSec; nt++) {
375 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
376 dum = pi*nt/(2.0*n*n);
377 if (std::fabs(dum) < 1.0) {
378 if( test >= 1.0e-10 )anpn += dum*test;
379 } else {
380 anpn += dum*test;
381 }
382 }
383
384 G4double ran = G4UniformRand();
385 G4double excs = 0.0;
386 if( targetCode == protonCode )
387 {
388 counter = -1;
389 for (np=0; np<numSec/3; np++) {
390 for (nm=Imax(0,np-1); nm<=(np+1); nm++) {
391 for (nz=0; nz<numSec/3; nz++) {
392 if (++counter < numMul) {
393 nt = np+nm+nz;
394 if ( (nt>0) && (nt<=numSec) ) {
395 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
396 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
397 if (std::fabs(dum) < 1.0) {
398 if( test >= 1.0e-10 )excs += dum*test;
399 } else {
400 excs += dum*test;
401 }
402 if (ran < excs) goto outOfLoop; //----------------------->
403 }
404 }
405 }
406 }
407 }
408
409 // 3 previous loops continued to the end
410 inElastic = false; // quasi-elastic scattering
411 return;
412 }
413 else
414 { // target must be a neutron
415 counter = -1;
416 for (np=0; np<numSec/3; np++) {
417 for (nm=np; nm<=(np+2); nm++) {
418 for (nz=0; nz<numSec/3; nz++) {
419 if (++counter < numMul) {
420 nt = np+nm+nz;
421 if ( (nt>=1) && (nt<=numSec) ) {
422 test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
423 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
424 if (std::fabs(dum) < 1.0) {
425 if( test >= 1.0e-10 )excs += dum*test;
426 } else {
427 excs += dum*test;
428 }
429 if (ran < excs) goto outOfLoop; // ----------------->
430 }
431 }
432 }
433 }
434 }
435 // 3 previous loops continued to the end
436 inElastic = false; // quasi-elastic scattering.
437 return;
438 }
439 }
440 outOfLoop: // <-----------------------------------------------
441
442 if( targetCode == protonCode)
443 {
444 if( np == (1+nm))
445 {
446 pv[1] = Neutron;
447 }
448 else if (np == nm)
449 {
450 if( G4UniformRand() < 0.75)
451 {
452 }
453 else
454 {
455 pv[0] = PionZero;
456 pv[1] = Neutron;
457 }
458 }
459 else
460 {
461 pv[0] = PionZero;
462 }
463 }
464 else
465 {
466 if( np == (nm-1))
467 {
468 if( G4UniformRand() < 0.5)
469 {
470 pv[1] = Proton;
471 }
472 else
473 {
474 pv[0] = PionZero;
475 }
476 }
477 else if ( np == nm)
478 {
479 }
480 else
481 {
482 pv[0] = PionZero;
483 pv[1] = Proton;
484 }
485 }
486
487
488 nt = np + nm + nz;
489 while ( nt > 0)
490 {
491 G4double ran = G4UniformRand();
492 if ( ran < (G4double)np/nt)
493 {
494 if( np > 0 )
495 { pv[vecLen++] = PionPlus;
496 np--;
497 }
498 }
499 else if ( ran < (G4double)(np+nm)/nt)
500 {
501 if( nm > 0 )
502 {
503 pv[vecLen++] = PionMinus;
504 nm--;
505 }
506 }
507 else
508 {
509 if( nz > 0 )
510 {
511 pv[vecLen++] = PionZero;
512 nz--;
513 }
514 }
515 nt = np + nm + nz;
516 }
517 if (verboseLevel > 1)
518 {
519 G4cout << "Particles produced: " ;
520 G4cout << pv[0].getName() << " " ;
521 G4cout << pv[1].getName() << " " ;
522 for (i=2; i < vecLen; i++)
523 {
524 G4cout << pv[i].getName() << " " ;
525 }
526 G4cout << G4endl;
527 }
528 return;
529 }
530
531
532
533
534
535
536
537
538
Note: See TracBrowser for help on using the repository browser.