source: trunk/source/processes/hadronic/models/low_energy/src/G4LEAntiKaonZeroInelastic.cc@ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 18.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: G4LEAntiKaonZeroInelastic.cc,v 1.10 2006/06/29 20:44:39 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30 // Hadronic Process: Low Energy KaonZeroLong Inelastic Process
31 // J.L. Chuma, TRIUMF, 11-Feb-1997
32 // Last modified: 27-Mar-1997
33 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
34
35#include "G4LEAntiKaonZeroInelastic.hh"
36#include "Randomize.hh"
37#include "G4HadReentrentException.hh"
38
39 G4HadFinalState *
40 G4LEAntiKaonZeroInelastic::ApplyYourself( const G4HadProjectile &aTrack,
41 G4Nucleus &targetNucleus )
42 {
43 const G4HadProjectile *originalIncident = &aTrack;
44 //
45 // create the target particle
46 //
47 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
48
49 if( verboseLevel > 1 )
50 {
51 const G4Material *targetMaterial = aTrack.GetMaterial();
52 G4cout << "G4LEAntiKaonZeroInelastic::ApplyYourself called" << G4endl;
53 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
54 G4cout << "target material = " << targetMaterial->GetName() << ", ";
55 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
56 << G4endl;
57 }
58 //
59 // Fermi motion and evaporation
60 // As of Geant3, the Fermi energy calculation had not been Done
61 //
62 G4double ek = originalIncident->GetKineticEnergy()/MeV;
63 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
64 G4ReactionProduct modifiedOriginal;
65 modifiedOriginal = *originalIncident;
66
67 G4double tkin = targetNucleus.Cinema( ek );
68 ek += tkin;
69 modifiedOriginal.SetKineticEnergy( ek*MeV );
70 G4double et = ek + amas;
71 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
72 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
73 if( pp > 0.0 )
74 {
75 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
76 modifiedOriginal.SetMomentum( momentum * (p/pp) );
77 }
78 //
79 // calculate black track energies
80 //
81 tkin = targetNucleus.EvaporationEffects( ek );
82 ek -= tkin;
83 modifiedOriginal.SetKineticEnergy( ek*MeV );
84 et = ek + amas;
85 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
86 pp = modifiedOriginal.GetMomentum().mag()/MeV;
87 if( pp > 0.0 )
88 {
89 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
90 modifiedOriginal.SetMomentum( momentum * (p/pp) );
91 }
92 G4ReactionProduct currentParticle = modifiedOriginal;
93 G4ReactionProduct targetParticle;
94 targetParticle = *originalTarget;
95 currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
96 targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
97 G4bool incidentHasChanged = false;
98 G4bool targetHasChanged = false;
99 G4bool quasiElastic = false;
100 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
101 G4int vecLen = 0;
102 vec.Initialize( 0 );
103
104 const G4double cutOff = 0.1;
105 if( currentParticle.GetKineticEnergy()/MeV > cutOff )
106 Cascade( vec, vecLen,
107 originalIncident, currentParticle, targetParticle,
108 incidentHasChanged, targetHasChanged, quasiElastic );
109
110 try
111 {
112 CalculateMomenta( vec, vecLen,
113 originalIncident, originalTarget, modifiedOriginal,
114 targetNucleus, currentParticle, targetParticle,
115 incidentHasChanged, targetHasChanged, quasiElastic );
116 }
117 catch(G4HadReentrentException aR)
118 {
119 aR.Report(G4cout);
120 throw G4HadReentrentException(__FILE__, __LINE__, "Bailing out");
121 }
122 SetUpChange( vec, vecLen,
123 currentParticle, targetParticle,
124 incidentHasChanged );
125
126 delete originalTarget;
127 return &theParticleChange;
128 }
129
130 void
131 G4LEAntiKaonZeroInelastic::Cascade(
132 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
133 G4int& vecLen,
134 const G4HadProjectile *originalIncident,
135 G4ReactionProduct &currentParticle,
136 G4ReactionProduct &targetParticle,
137 G4bool &incidentHasChanged,
138 G4bool &targetHasChanged,
139 G4bool &quasiElastic )
140 {
141 // derived from original FORTRAN code CASK0B by H. Fesefeldt (13-Sep-1987)
142 //
143 // K0Long undergoes interaction with nucleon within a nucleus. Check if it is
144 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
145 // occurs and input particle is degraded in energy. No other particles are produced.
146 // If reaction is possible, find the correct number of pions/protons/neutrons
147 // produced using an interpolation to multiplicity data. Replace some pions or
148 // protons/neutrons by kaons or strange baryons according to the average
149 // multiplicity per Inelastic reaction.
150 //
151 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
152 const G4double etOriginal = originalIncident->Get4Momentum().e()/MeV;
153 const G4double pOriginal = originalIncident->Get4Momentum().vect().mag()/MeV;
154 const G4double targetMass = targetParticle.GetMass()/MeV;
155 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
156 targetMass*targetMass +
157 2.0*targetMass*etOriginal );
158 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
159
160 static G4bool first = true;
161 const G4int numMul = 1200;
162 const G4int numSec = 60;
163 static G4double protmul[numMul], protnorm[numSec]; // proton constants
164 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
165 // np = number of pi+, nm = number of pi-, nz = number of pi0
166 G4int counter, nt=0, np=0, nm=0, nz=0;
167 const G4double c = 1.25;
168 const G4double b[] = { 0.7, 0.7 };
169 if( first ) // compute normalization constants, this will only be Done once
170 {
171 first = false;
172 G4int i;
173 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
174 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
175 counter = -1;
176 for( np=0; np<numSec/3; ++np )
177 {
178 for( nm=std::max(0,np-2); nm<=np; ++nm )
179 {
180 for( nz=0; nz<numSec/3; ++nz )
181 {
182 if( ++counter < numMul )
183 {
184 nt = np+nm+nz;
185 if( nt>0 && nt<=numSec )
186 {
187 protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
188 protnorm[nt-1] += protmul[counter];
189 }
190 }
191 }
192 }
193 }
194 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
195 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
196 counter = -1;
197 for( np=0; np<(numSec/3); ++np )
198 {
199 for( nm=std::max(0,np-1); nm<=(np+1); ++nm )
200 {
201 for( nz=0; nz<numSec/3; ++nz )
202 {
203 if( ++counter < numMul )
204 {
205 nt = np+nm+nz;
206 if( nt>0 && nt<=numSec )
207 {
208 neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
209 neutnorm[nt-1] += neutmul[counter];
210 }
211 }
212 }
213 }
214 }
215 for( i=0; i<numSec; ++i )
216 {
217 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
218 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
219 }
220 } // end of initialization
221
222 const G4double expxu = 82.; // upper bound for arg. of exp
223 const G4double expxl = -expxu; // lower bound for arg. of exp
224 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
225 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
226 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
227 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
228 G4ParticleDefinition *aProton = G4Proton::Proton();
229 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
230 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
231 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
232 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
233 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
234 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
235 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
236 const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
237 G4int iplab = G4int(std::min( 9.0, 5.0*pOriginal*MeV/GeV ));
238 if( (pOriginal*MeV/GeV <= 2.0) && (G4UniformRand() < cech[iplab]) )
239 {
240 np = nm = nz = nt = 0;
241 iplab = G4int(std::min( 19.0, pOriginal*MeV/GeV*10.0 ));
242 const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
243 0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
244 if( G4UniformRand() > cnk0[iplab] )
245 {
246 G4double ran = G4UniformRand();
247 if( ran < 0.25 ) // k0Long n --> pi- s+
248 {
249 if( targetParticle.GetDefinition() == aNeutron )
250 {
251 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
252 targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
253 incidentHasChanged = true;
254 targetHasChanged = true;
255 }
256 }
257 else if( ran < 0.50 ) // k0Long p --> pi+ s0 or k0Long n --> pi0 s0
258 {
259 if( targetParticle.GetDefinition() == aNeutron )
260 currentParticle.SetDefinitionAndUpdateE( aPiZero );
261 else
262 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
263 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
264 incidentHasChanged = true;
265 targetHasChanged = true;
266 }
267 else if( ran < 0.75 ) // k0Long n --> pi+ s-
268 {
269 if( targetParticle.GetDefinition() == aNeutron )
270 {
271 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
272 targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
273 incidentHasChanged = true;
274 targetHasChanged = true;
275 }
276 }
277 else // k0Long p --> pi+ L or k0Long n --> pi0 L
278 {
279 if( targetParticle.GetDefinition() == aNeutron )
280 currentParticle.SetDefinitionAndUpdateE( aPiZero );
281 else
282 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
283 targetParticle.SetDefinitionAndUpdateE( aLambda );
284 incidentHasChanged = true;
285 targetHasChanged = true;
286 }
287 }
288 else // ran <= cnk0
289 {
290 quasiElastic = true;
291 if( targetParticle.GetDefinition() == aNeutron )
292 {
293 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
294 targetParticle.SetDefinitionAndUpdateE( aProton );
295 incidentHasChanged = true;
296 targetHasChanged = true;
297 }
298 }
299 }
300 else // (pOriginal > 2.0*GeV) || (random number >= cech[iplab])
301 {
302 if( availableEnergy < aPiPlus->GetPDGMass()/MeV )
303 {
304 quasiElastic = true;
305 return;
306 }
307 G4double n, anpn;
308 GetNormalizationConstant( availableEnergy, n, anpn );
309 G4double ran = G4UniformRand();
310 G4double dum, test, excs = 0.0;
311 if( targetParticle.GetDefinition() == aProton )
312 {
313 counter = -1;
314 for( np=0; (np<numSec/3) && (ran>=excs); ++np )
315 {
316 for( nm=std::max(0,np-2); nm<=np && ran>=excs; ++nm )
317 {
318 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
319 {
320 if( ++counter < numMul )
321 {
322 nt = np+nm+nz;
323 if( nt>0 && nt<=numSec )
324 {
325 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
326 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
327 if( std::fabs(dum) < 1.0 )
328 {
329 if( test >= 1.0e-10 )excs += dum*test;
330 }
331 else
332 excs += dum*test;
333 }
334 }
335 }
336 }
337 }
338 if( ran >= excs ) // 3 previous loops continued to the end
339 {
340 quasiElastic = true;
341 return;
342 }
343 np--; nm--; nz--;
344 switch( np-nm )
345 {
346 case 1:
347 if( G4UniformRand() < 0.5 )
348 {
349 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
350 incidentHasChanged = true;
351 }
352 else
353 {
354 targetParticle.SetDefinitionAndUpdateE( aNeutron );
355 targetHasChanged = true;
356 }
357 case 0:
358 break;
359 default:
360 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
361 targetParticle.SetDefinitionAndUpdateE( aNeutron );
362 incidentHasChanged = true;
363 targetHasChanged = true;
364 break;
365 }
366 }
367 else // target must be a neutron
368 {
369 counter = -1;
370 for( np=0; np<numSec/3 && ran>=excs; ++np )
371 {
372 for( nm=std::max(0,np-1); nm<=(np+1) && ran>=excs; ++nm )
373 {
374 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
375 {
376 if( ++counter < numMul )
377 {
378 nt = np+nm+nz;
379 if( nt>0 && nt<=numSec )
380 {
381 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
382 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
383 if( std::fabs(dum) < 1.0 )
384 {
385 if( test >= 1.0e-10 )excs += dum*test;
386 }
387 else
388 excs += dum*test;
389 }
390 }
391 }
392 }
393 }
394 if( ran >= excs ) // 3 previous loops continued to the end
395 {
396 quasiElastic = true;
397 return;
398 }
399 np--; nm--; nz--;
400 switch( np-nm )
401 {
402 case 0:
403 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
404 targetParticle.SetDefinitionAndUpdateE( aProton );
405 incidentHasChanged = true;
406 targetHasChanged = true;
407 break;
408 case 1:
409 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
410 incidentHasChanged = true;
411 break;
412 default:
413 targetParticle.SetDefinitionAndUpdateE( aProton );
414 targetHasChanged = true;
415 break;
416 }
417 }
418 if( G4UniformRand() >= 0.5 )
419 {
420 if( currentParticle.GetDefinition() == aKaonMinus &&
421 targetParticle.GetDefinition() == aNeutron )
422 {
423 ran = G4UniformRand();
424 if( ran < 0.68 )
425 {
426 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
427 targetParticle.SetDefinitionAndUpdateE( aLambda );
428 }
429 else if( ran < 0.84 )
430 {
431 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
432 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
433 }
434 else
435 {
436 currentParticle.SetDefinitionAndUpdateE( aPiZero );
437 targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
438 }
439 }
440 else if( (currentParticle.GetDefinition() == aKaonZS ||
441 currentParticle.GetDefinition() == aKaonZL ) &&
442 targetParticle.GetDefinition() == aProton )
443 {
444 ran = G4UniformRand();
445 if( ran < 0.68 )
446 {
447 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
448 targetParticle.SetDefinitionAndUpdateE( aLambda );
449 }
450 else if( ran < 0.84 )
451 {
452 currentParticle.SetDefinitionAndUpdateE( aPiZero );
453 targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
454 }
455 else
456 {
457 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
458 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
459 }
460 }
461 else
462 {
463 ran = G4UniformRand();
464 if( ran < 0.67 )
465 {
466 currentParticle.SetDefinitionAndUpdateE( aPiZero );
467 targetParticle.SetDefinitionAndUpdateE( aLambda );
468 }
469 else if( ran < 0.78 )
470 {
471 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
472 targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
473 }
474 else if( ran < 0.89 )
475 {
476 currentParticle.SetDefinitionAndUpdateE( aPiZero );
477 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
478 }
479 else
480 {
481 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
482 targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
483 }
484 }
485 incidentHasChanged = true;
486 targetHasChanged = true;
487 }
488 }
489 if( currentParticle.GetDefinition() == aKaonZL )
490 {
491 if( G4UniformRand() >= 0.5 )
492 {
493 currentParticle.SetDefinitionAndUpdateE( aKaonZS );
494 incidentHasChanged = true;
495 }
496 }
497 if( targetParticle.GetDefinition() == aKaonZL )
498 {
499 if( G4UniformRand() >= 0.5 )
500 {
501 targetParticle.SetDefinitionAndUpdateE( aKaonZS );
502 targetHasChanged = true;
503 }
504 }
505 SetUpPions( np, nm, nz, vec, vecLen );
506 return;
507 }
508
509 /* end of file */
510
Note: See TracBrowser for help on using the repository browser.