source: trunk/source/processes/hadronic/models/low_energy/src/G4LEAntiLambdaInelastic.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: 22.1 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: G4LEAntiLambdaInelastic.cc,v 1.11 2006/06/29 20:44:41 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30 // Hadronic Process: AntiLambda Inelastic Process
31 // J.L. Chuma, TRIUMF, 19-Feb-1997
32 // Last modified: 27-Mar-1997
33 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
34
35#include "G4LEAntiLambdaInelastic.hh"
36#include "Randomize.hh"
37
38 G4HadFinalState *
39 G4LEAntiLambdaInelastic::ApplyYourself( const G4HadProjectile &aTrack,
40 G4Nucleus &targetNucleus )
41 {
42 const G4HadProjectile *originalIncident = &aTrack;
43 //
44 // create the target particle
45 //
46 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
47
48 if( verboseLevel > 1 )
49 {
50 const G4Material *targetMaterial = aTrack.GetMaterial();
51 G4cout << "G4LEAntiLambdaInelastic::ApplyYourself called" << G4endl;
52 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
53 G4cout << "target material = " << targetMaterial->GetName() << ", ";
54 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
55 << G4endl;
56 }
57 //
58 // Fermi motion and evaporation
59 // As of Geant3, the Fermi energy calculation had not been Done
60 //
61 G4double ek = originalIncident->GetKineticEnergy()/MeV;
62 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
63 G4ReactionProduct modifiedOriginal;
64 modifiedOriginal = *originalIncident;
65
66 G4double tkin = targetNucleus.Cinema( ek );
67 ek += tkin;
68 modifiedOriginal.SetKineticEnergy( ek*MeV );
69 G4double et = ek + amas;
70 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
71 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
72 if( pp > 0.0 )
73 {
74 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
75 modifiedOriginal.SetMomentum( momentum * (p/pp) );
76 }
77 //
78 // calculate black track energies
79 //
80 tkin = targetNucleus.EvaporationEffects( ek );
81 ek -= tkin;
82 modifiedOriginal.SetKineticEnergy( ek*MeV );
83 et = ek + amas;
84 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
85 pp = modifiedOriginal.GetMomentum().mag()/MeV;
86 if( pp > 0.0 )
87 {
88 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
89 modifiedOriginal.SetMomentum( momentum * (p/pp) );
90 }
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 const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
106 if( (originalIncident->GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
107 Cascade( vec, vecLen,
108 originalIncident, currentParticle, targetParticle,
109 incidentHasChanged, targetHasChanged, quasiElastic );
110
111 CalculateMomenta( vec, vecLen,
112 originalIncident, originalTarget, modifiedOriginal,
113 targetNucleus, currentParticle, targetParticle,
114 incidentHasChanged, targetHasChanged, quasiElastic );
115
116 SetUpChange( vec, vecLen,
117 currentParticle, targetParticle,
118 incidentHasChanged );
119
120 delete originalTarget;
121 return &theParticleChange;
122 }
123
124 void
125 G4LEAntiLambdaInelastic::Cascade(
126 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
127 G4int &vecLen,
128 const G4HadProjectile *originalIncident,
129 G4ReactionProduct &currentParticle,
130 G4ReactionProduct &targetParticle,
131 G4bool &incidentHasChanged,
132 G4bool &targetHasChanged,
133 G4bool &quasiElastic )
134 {
135 // derived from original FORTRAN code CASAL0 by H. Fesefeldt (13-Sep-1987)
136 //
137 // AntiLambda undergoes interaction with nucleon within a nucleus. Check if it is
138 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
139 // occurs and input particle is degraded in energy. No other particles are produced.
140 // If reaction is possible, find the correct number of pions/protons/neutrons
141 // produced using an interpolation to multiplicity data. Replace some pions or
142 // protons/neutrons by kaons or strange baryons according to the average
143 // multiplicity per Inelastic reaction.
144 //
145 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
146 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
147 const G4double targetMass = targetParticle.GetMass()/MeV;
148 const G4double pOriginal = originalIncident->GetTotalMomentum()/GeV;
149 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
150 targetMass*targetMass +
151 2.0*targetMass*etOriginal );
152 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
153
154 static G4bool first = true;
155 const G4int numMul = 1200;
156 const G4int numMulA = 400;
157 const G4int numSec = 60;
158 static G4double protmul[numMul], protnorm[numSec]; // proton constants
159 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
160 static G4double protmulA[numMulA], protnormA[numSec]; // proton constants
161 static G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
162 // np = number of pi+, nm = number of pi-, nz = number of pi0
163 G4int nt=0, np=0, nm=0, nz=0;
164 G4double test;
165 const G4double c = 1.25;
166 const G4double b[] = { 0.7, 0.7 };
167 if( first ) // compute normalization constants, this will only be Done once
168 {
169 first = false;
170 G4int i;
171 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
172 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
173 G4int counter = -1;
174 for( np=0; np<(numSec/3); ++np )
175 {
176 for( nm=std::max(0,np-2); nm<=(np+1); ++nm )
177 {
178 for( nz=0; nz<numSec/3; ++nz )
179 {
180 if( ++counter < numMul )
181 {
182 nt = np+nm+nz;
183 if( nt>0 && nt<=numSec )
184 {
185 protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
186 protnorm[nt-1] += protmul[counter];
187 }
188 }
189 }
190 }
191 }
192 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
193 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
194 counter = -1;
195 for( np=0; np<numSec/3; ++np )
196 {
197 for( nm=std::max(0,np-1); nm<=(np+2); ++nm )
198 {
199 for( nz=0; nz<numSec/3; ++nz )
200 {
201 if( ++counter < numMul )
202 {
203 nt = np+nm+nz;
204 if( nt>0 && nt<=numSec )
205 {
206 neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
207 neutnorm[nt-1] += neutmul[counter];
208 }
209 }
210 }
211 }
212 }
213 for( i=0; i<numSec; ++i )
214 {
215 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
216 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
217 }
218 //
219 // do the same for annihilation channels
220 //
221 for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
222 for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
223 counter = -1;
224 for( np=1; np<(numSec/3); ++np )
225 {
226 nm = np-1;
227 for( nz=0; nz<numSec/3; ++nz )
228 {
229 if( ++counter < numMulA )
230 {
231 nt = np+nm+nz;
232 if( nt>1 && nt<=numSec )
233 {
234 protmulA[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
235 protnormA[nt-1] += protmulA[counter];
236 }
237 }
238 }
239 }
240 for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
241 for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
242 counter = -1;
243 for( np=0; np<numSec/3; ++np )
244 {
245 nm = np;
246 for( nz=0; nz<numSec/3; ++nz )
247 {
248 if( ++counter < numMulA )
249 {
250 nt = np+nm+nz;
251 if( nt>1 && nt<=numSec )
252 {
253 neutmulA[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
254 neutnormA[nt-1] += neutmulA[counter];
255 }
256 }
257 }
258 }
259 for( i=0; i<numSec; ++i )
260 {
261 if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
262 if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
263 }
264 } // end of initialization
265 const G4double expxu = 82.; // upper bound for arg. of exp
266 const G4double expxl = -expxu; // lower bound for arg. of exp
267
268 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
269 G4ParticleDefinition *aProton = G4Proton::Proton();
270 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
271 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
272 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
273 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
274 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
275 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
276 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
277 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
278 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
279 const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
280 0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
281 0.39,0.36,0.33,0.10,0.01};
282 G4int iplab = G4int( pOriginal*10.0 );
283 if( iplab > 9 )iplab = G4int( (pOriginal- 1.0)*5.0 ) + 10;
284 if( iplab > 14 )iplab = G4int( pOriginal- 2.0 ) + 15;
285 if( iplab > 22 )iplab = G4int( (pOriginal-10.0)/10.0 ) + 23;
286 if( iplab > 24 )iplab = 24;
287 if( G4UniformRand() > anhl[iplab] )
288 {
289 if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
290 { // not energetically possible to produce pion(s)
291 quasiElastic = true;
292 return;
293 }
294 G4double n, anpn;
295 GetNormalizationConstant( availableEnergy, n, anpn );
296 G4double ran = G4UniformRand();
297 G4double dum, excs = 0.0;
298 if( targetParticle.GetDefinition() == aProton )
299 {
300 G4int counter = -1;
301 for( np=0; np<numSec/3 && ran>=excs; ++np )
302 {
303 for( nm=std::max(0,np-2); nm<=(np+1) && ran>=excs; ++nm )
304 {
305 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
306 {
307 if( ++counter < numMul )
308 {
309 nt = np+nm+nz;
310 if( nt>0 && nt<=numSec )
311 {
312 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
313 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
314 if( std::fabs(dum) < 1.0 )
315 {
316 if( test >= 1.0e-10 )excs += dum*test;
317 }
318 else
319 excs += dum*test;
320 }
321 }
322 }
323 }
324 }
325 if( ran >= excs ) // 3 previous loops continued to the end
326 {
327 quasiElastic = true;
328 return;
329 }
330 np--; nm--; nz--;
331 G4int ncht = std::min( 4, std::max( 1, np-nm+2 ) );
332 switch( ncht )
333 {
334 case 1:
335 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
336 incidentHasChanged = true;
337 break;
338 case 2:
339 if( G4UniformRand() < 0.5 )
340 {
341 if( G4UniformRand() < 0.5 )
342 {
343 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
344 incidentHasChanged = true;
345 }
346 else
347 {
348 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
349 incidentHasChanged = true;
350 targetParticle.SetDefinitionAndUpdateE( aNeutron );
351 targetHasChanged = true;
352 }
353 }
354 else
355 {
356 if( G4UniformRand() >= 0.5 )
357 {
358 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
359 incidentHasChanged = true;
360 targetParticle.SetDefinitionAndUpdateE( aNeutron );
361 targetHasChanged = true;
362 }
363 }
364 break;
365 case 3:
366 if( G4UniformRand() < 0.5 )
367 {
368 if( G4UniformRand() < 0.5 )
369 {
370 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
371 incidentHasChanged = true;
372 targetParticle.SetDefinitionAndUpdateE( aNeutron );
373 targetHasChanged = true;
374 }
375 else
376 {
377 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
378 incidentHasChanged = true;
379 }
380 }
381 else
382 {
383 if( G4UniformRand() < 0.5 )
384 {
385 targetParticle.SetDefinitionAndUpdateE( aNeutron );
386 targetHasChanged = true;
387 }
388 else
389 {
390 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
391 incidentHasChanged = true;
392 }
393 }
394 break;
395 case 4:
396 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
397 incidentHasChanged = true;
398 targetParticle.SetDefinitionAndUpdateE( aNeutron );
399 targetHasChanged = true;
400 break;
401 }
402 }
403 else // target must be a neutron
404 {
405 G4int counter = -1;
406 for( np=0; np<numSec/3 && ran>=excs; ++np )
407 {
408 for( nm=std::max(0,np-1); nm<=(np+2) && ran>=excs; ++nm )
409 {
410 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
411 {
412 if( ++counter < numMul )
413 {
414 nt = np+nm+nz;
415 if( nt>0 && nt<=numSec )
416 {
417 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
418 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
419 if( std::fabs(dum) < 1.0 )
420 {
421 if( test >= 1.0e-10 )excs += dum*test;
422 }
423 else
424 excs += dum*test;
425 }
426 }
427 }
428 }
429 }
430 if( ran >= excs ) // 3 previous loops continued to the end
431 {
432 quasiElastic = true;
433 return;
434 }
435 np--; nm--; nz--;
436 G4int ncht = std::min( 4, std::max( 1, np-nm+3 ) );
437 switch( ncht )
438 {
439 case 1:
440 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
441 incidentHasChanged = true;
442 targetParticle.SetDefinitionAndUpdateE( aProton );
443 targetHasChanged = true;
444 break;
445 case 2:
446 if( G4UniformRand() < 0.5 )
447 {
448 if( G4UniformRand() < 0.5 )
449 {
450 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
451 incidentHasChanged = true;
452 targetParticle.SetDefinitionAndUpdateE( aProton );
453 targetHasChanged = true;
454 }
455 else
456 {
457 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
458 incidentHasChanged = true;
459 }
460 }
461 else
462 {
463 if( G4UniformRand() < 0.5 )
464 {
465 targetParticle.SetDefinitionAndUpdateE( aProton );
466 targetHasChanged = true;
467 }
468 else
469 {
470 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
471 incidentHasChanged = true;
472 }
473 }
474 break;
475 case 3:
476 if( G4UniformRand() < 0.5 )
477 {
478 if( G4UniformRand() < 0.5 )
479 {
480 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
481 incidentHasChanged = true;
482 }
483 else
484 {
485 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
486 incidentHasChanged = true;
487 targetParticle.SetDefinitionAndUpdateE( aProton );
488 targetHasChanged = true;
489 }
490 }
491 else
492 {
493 if( G4UniformRand() >= 0.5 )
494 {
495 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
496 incidentHasChanged = true;
497 targetParticle.SetDefinitionAndUpdateE( aProton );
498 targetHasChanged = true;
499 }
500 }
501 break;
502 default:
503 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
504 incidentHasChanged = true;
505 break;
506 }
507 }
508 }
509 else // random number <= anhl[iplab]
510 {
511 if( centerofmassEnergy <= aPiPlus->GetPDGMass()/MeV+aKaonPlus->GetPDGMass()/MeV )
512 {
513 quasiElastic = true;
514 return;
515 }
516 //
517 // annihilation channels
518 //
519 G4double n, anpn;
520 GetNormalizationConstant( -centerofmassEnergy, n, anpn );
521 G4double ran = G4UniformRand();
522 G4double dum, excs = 0.0;
523 if( targetParticle.GetDefinition() == aProton )
524 {
525 G4int counter = -1;
526 for( np=1; np<numSec/3 && ran>=excs; ++np )
527 {
528 nm = np-1;
529 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
530 {
531 if( ++counter < numMulA )
532 {
533 nt = np+nm+nz;
534 if( nt>1 && nt<=numSec )
535 {
536 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
537 dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
538 if( std::fabs(dum) < 1.0 )
539 {
540 if( test >= 1.0e-10 )excs += dum*test;
541 }
542 else
543 excs += dum*test;
544 }
545 }
546 }
547 }
548 }
549 else // target must be a neutron
550 {
551 G4int counter = -1;
552 for( np=0; np<numSec/3 && ran>=excs; ++np )
553 {
554 nm = np;
555 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
556 {
557 if( ++counter < numMulA )
558 {
559 nt = np+nm+nz;
560 if( nt>1 && nt<=numSec )
561 {
562 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
563 dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
564 if( std::fabs(dum) < 1.0 )
565 {
566 if( test >= 1.0e-10 )excs += dum*test;
567 }
568 else
569 excs += dum*test;
570 }
571 }
572 }
573 }
574 }
575 if( ran >= excs ) // 3 previous loops continued to the end
576 {
577 quasiElastic = true;
578 return;
579 }
580 np--; nz--;
581 currentParticle.SetMass( 0.0 );
582 targetParticle.SetMass( 0.0 );
583 }
584 SetUpPions( np, nm, nz, vec, vecLen );
585 if( currentParticle.GetMass() == 0.0 )
586 {
587 if( nz == 0 )
588 {
589 if( nm > 0 )
590 {
591 for( G4int i=0; i<vecLen; ++i )
592 {
593 if( vec[i]->GetDefinition() == aPiMinus )
594 {
595 vec[i]->SetDefinitionAndUpdateE( aKaonMinus );
596 break;
597 }
598 }
599 }
600 }
601 else // nz > 0
602 {
603 if( nm == 0 )
604 {
605 for( G4int i=0; i<vecLen; ++i )
606 {
607 if( vec[i]->GetDefinition() == aPiZero )
608 {
609 vec[i]->SetDefinitionAndUpdateE( aKaonZL );
610 break;
611 }
612 }
613 }
614 else // nm > 0
615 {
616 if( G4UniformRand() < 0.5 )
617 {
618 if( nm > 0 )
619 {
620 for( G4int i=0; i<vecLen; ++i )
621 {
622 if( vec[i]->GetDefinition() == aPiMinus )
623 {
624 vec[i]->SetDefinitionAndUpdateE( aKaonMinus );
625 break;
626 }
627 }
628 }
629 }
630 else // random number >= 0.5
631 {
632 for( G4int i=0; i<vecLen; ++i )
633 {
634 if( vec[i]->GetDefinition() == aPiZero )
635 {
636 vec[i]->SetDefinitionAndUpdateE( aKaonZL );
637 break;
638 }
639 }
640 }
641 }
642 }
643 }
644 return;
645 }
646
647 /* end of file */
648
Note: See TracBrowser for help on using the repository browser.