source: trunk/source/processes/hadronic/models/low_energy/src/G4LEPionPlusInelastic.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: 14.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: G4LEPionPlusInelastic.cc,v 1.15 2007/02/24 06:28:52 dennis Exp $
28// GEANT4 tag $Name: geant4-09-02 $
29//
30 // Hadronic Process: PionPlus Inelastic Process
31 // J.L. Chuma, TRIUMF, 19-Nov-1996
32 // Last modified: 27-Mar-1997
33 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
34 // fixing charge exchange - HPW Sep 2002.
35
36#include "G4LEPionPlusInelastic.hh"
37#include "Randomize.hh"
38
39 G4HadFinalState *
40 G4LEPionPlusInelastic::ApplyYourself( const G4HadProjectile &aTrack,
41 G4Nucleus &targetNucleus )
42 {
43 const G4HadProjectile *originalIncident = &aTrack;
44 if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
45 {
46 theParticleChange.SetStatusChange(isAlive);
47 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
48 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
49 return &theParticleChange;
50 }
51
52 // create the target particle
53
54 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
55// G4double targetMass = originalTarget->GetDefinition()->GetPDGMass();
56 G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
57
58 if( verboseLevel > 1 )
59 {
60 const G4Material *targetMaterial = aTrack.GetMaterial();
61 G4cout << "G4LEPionPlusInelastic::ApplyYourself called" << G4endl;
62 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
63 G4cout << "target material = " << targetMaterial->GetName() << ", ";
64 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
65 << G4endl;
66 }
67 G4ReactionProduct currentParticle(
68 const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition() ) );
69 currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
70 currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
71
72 // Fermi motion and evaporation
73 // As of Geant3, the Fermi energy calculation had not been Done
74
75 G4double ek = originalIncident->GetKineticEnergy();
76 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
77
78 G4double tkin = targetNucleus.Cinema( ek );
79 ek += tkin;
80 currentParticle.SetKineticEnergy( ek );
81 G4double et = ek + amas;
82 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
83 G4double pp = currentParticle.GetMomentum().mag();
84 if( pp > 0.0 )
85 {
86 G4ThreeVector momentum = currentParticle.GetMomentum();
87 currentParticle.SetMomentum( momentum * (p/pp) );
88 }
89
90 // calculate black track energies
91
92 tkin = targetNucleus.EvaporationEffects( ek );
93 ek -= tkin;
94 currentParticle.SetKineticEnergy( ek );
95 et = ek + amas;
96 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
97 pp = currentParticle.GetMomentum().mag();
98 if( pp > 0.0 )
99 {
100 G4ThreeVector momentum = currentParticle.GetMomentum();
101 currentParticle.SetMomentum( momentum * (p/pp) );
102 }
103
104 G4ReactionProduct modifiedOriginal = currentParticle;
105
106 currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
107 targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
108 G4bool incidentHasChanged = false;
109 G4bool targetHasChanged = false;
110 G4bool quasiElastic = false;
111 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
112 G4int vecLen = 0;
113 vec.Initialize( 0 );
114
115 const G4double cutOff = 0.1*MeV;
116 if( currentParticle.GetKineticEnergy() > cutOff )
117 Cascade( vec, vecLen,
118 originalIncident, currentParticle, targetParticle,
119 incidentHasChanged, targetHasChanged, quasiElastic );
120
121 CalculateMomenta( vec, vecLen,
122 originalIncident, originalTarget, modifiedOriginal,
123 targetNucleus, currentParticle, targetParticle,
124 incidentHasChanged, targetHasChanged, quasiElastic );
125
126 SetUpChange( vec, vecLen,
127 currentParticle, targetParticle,
128 incidentHasChanged );
129
130 delete originalTarget;
131 return &theParticleChange;
132 }
133
134 void
135 G4LEPionPlusInelastic::Cascade(
136 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
137 G4int& vecLen,
138 const G4HadProjectile *originalIncident,
139 G4ReactionProduct &currentParticle,
140 G4ReactionProduct &targetParticle,
141 G4bool &incidentHasChanged,
142 G4bool &targetHasChanged,
143 G4bool &quasiElastic )
144 {
145 // derived from original FORTRAN code CASPIP by H. Fesefeldt (18-Sep-1987)
146 //
147 // pi+ undergoes interaction with nucleon within nucleus.
148 // Check if energetically possible to produce pions/kaons.
149 // If not assume nuclear excitation occurs and input particle
150 // is degraded in energy. No other particles produced.
151 // If reaction is possible find correct number of pions/protons/neutrons
152 // produced using an interpolation to multiplicity data.
153 // Replace some pions or protons/neutrons by kaons or strange baryons
154 // according to average multiplicity per inelastic reactions.
155 //
156 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
157 const G4double etOriginal = originalIncident->GetTotalEnergy();
158 const G4double pOriginal = originalIncident->GetTotalMomentum();
159 const G4double targetMass = targetParticle.GetMass();
160 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
161 targetMass*targetMass +
162 2.0*targetMass*etOriginal );
163 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
164 static G4bool first = true;
165 const G4int numMul = 1200;
166 const G4int numSec = 60;
167 static G4double protmul[numMul], protnorm[numSec]; // proton constants
168 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
169 // np = number of pi+, nm = number of pi-, nz = number of pi0
170 G4int counter, nt=0, np=0, nm=0, nz=0;
171 const G4double c = 1.25;
172 const G4double b[] = { 0.70, 0.70 };
173 if( first ) { // compute normalization constants, this will only be Done once
174 first = false;
175 G4int i;
176 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
177 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
178 counter = -1;
179 for( np=0; np<(numSec/3); ++np ) {
180 for( nm=std::max(0,np-2); nm<=np; ++nm ) {
181 for( nz=0; nz<numSec/3; ++nz ) {
182 if( ++counter < numMul ) {
183 nt = np+nm+nz;
184 if( nt > 0 ) {
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 for( nm=std::max(0,np-1); nm<=(np+1); ++nm ) {
197 for( nz=0; nz<numSec/3; ++nz ) {
198 if( ++counter < numMul ) {
199 nt = np+nm+nz;
200 if( (nt>0) && (nt<=numSec) ) {
201 neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
202 neutnorm[nt-1] += neutmul[counter];
203 }
204 }
205 }
206 }
207 }
208 for( i=0; i<numSec; ++i ) {
209 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
210 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
211 }
212 } // end of initialization
213
214 const G4double expxu = 82.; // upper bound for arg. of exp
215 const G4double expxl = -expxu; // lower bound for arg. of exp
216 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
217 G4ParticleDefinition *aProton = G4Proton::Proton();
218 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
219 G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
220 const G4double supp[] = {0.,0.2,0.45,0.55,0.65,0.75,0.85,0.90,0.94,0.98};
221 G4double test, w0, wp, wt, wm;
222 if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
223 {
224 // suppress high multiplicity events at low momentum
225 // only one pion will be produced
226 // charge exchange reaction is included in inelastic cross section
227
228 const G4double cech[] = {1.,0.95,0.79,0.32,0.19,0.16,0.14,0.12,0.10,0.08};
229 G4int iplab = G4int(std::min( 9.0, pOriginal/GeV*5.0 ));
230 if( G4UniformRand() <= cech[iplab] )
231 {
232 if( targetParticle.GetDefinition() == aNeutron )
233 {
234 currentParticle.SetDefinitionAndUpdateE( aPiZero ); // charge exchange
235 targetParticle.SetDefinitionAndUpdateE( aProton );
236 incidentHasChanged = true;
237 targetHasChanged = true;
238 }
239 }
240
241 if( availableEnergy <= G4PionMinus::PionMinus()->GetPDGMass() )
242 {
243 quasiElastic = true;
244 return;
245 }
246
247 nm = np = nz = 0;
248 if( targetParticle.GetDefinition() == aProton ) {
249 test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[0])/(2.0*c*c) ) ) );
250 w0 = test;
251 wp = test;
252 if( G4UniformRand() < w0/(w0+wp) )
253 nz =1;
254 else
255 np = 1;
256 } else { // target is a neutron
257 test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[1])/(2.0*c*c) ) ) );
258 w0 = test;
259 wp = test;
260 test = std::exp( std::min( expxu, std::max( expxl, -sqr(-1.0+b[1])/(2.0*c*c) ) ) );
261 wm = test;
262 wt = w0+wp+wm;
263 wp = w0+wp;
264 G4double ran = G4UniformRand();
265 if( ran < w0/wt )
266 nz = 1;
267 else if( ran < wp/wt )
268 np = 1;
269 else
270 nm = 1;
271 }
272 } else {
273 if( availableEnergy <= G4PionMinus::PionMinus()->GetPDGMass() )
274 {
275 quasiElastic = true;
276 return;
277 }
278 G4double n, anpn;
279 GetNormalizationConstant( availableEnergy, n, anpn );
280 G4double ran = G4UniformRand();
281 G4double dum, excs = 0.0;
282 if( targetParticle.GetDefinition() == aProton ) {
283 counter = -1;
284 for( np=0; (np<numSec/3) && (ran>=excs); ++np ) {
285 for( nm=std::max(0,np-2); (nm<=np) && (ran>=excs); ++nm ) {
286 for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz ) {
287 if( ++counter < numMul ) {
288 nt = np+nm+nz;
289 if( nt > 0 ) {
290 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
291 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
292 if( std::fabs(dum) < 1.0 ) {
293 if( test >= 1.0e-10 )excs += dum*test;
294 } else {
295 excs += dum*test;
296 }
297 }
298 }
299 }
300 }
301 }
302 if( ran >= excs )
303 {
304 quasiElastic = true;
305 return; // 3 previous loops continued to the end
306 }
307 np--; nm--; nz--;
308 } else { // target must be a neutron
309 counter = -1;
310 for( np=0; (np<numSec/3) && (ran>=excs); ++np ) {
311 for( nm=std::max(0,np-1); (nm<=(np+1)) && (ran>=excs); ++nm ) {
312 for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz ) {
313 if( ++counter < numMul ) {
314 nt = np+nm+nz;
315 if( (nt>=1) && (nt<=numSec) ) {
316 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
317 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
318 if( std::fabs(dum) < 1.0 ) {
319 if( test >= 1.0e-10 )excs += dum*test;
320 } else {
321 excs += dum*test;
322 }
323 }
324 }
325 }
326 }
327 }
328 if( ran >= excs ) // 3 previous loops continued to the end
329 {
330 quasiElastic = true;
331 return; // 3 previous loops continued to the end
332 }
333 np--; nm--; nz--;
334 }
335 }
336 if( targetParticle.GetDefinition() == aProton ) {
337 switch( np-nm ) {
338 case 1:
339 if( G4UniformRand() < 0.5 ) {
340 currentParticle.SetDefinitionAndUpdateE( aPiZero );
341 incidentHasChanged = true;
342 } else {
343 targetParticle.SetDefinitionAndUpdateE( aNeutron );
344 targetHasChanged = true;
345 }
346 break;
347 case 2:
348 currentParticle.SetDefinitionAndUpdateE( aPiZero );
349 targetParticle.SetDefinitionAndUpdateE( aNeutron );
350 incidentHasChanged = true;
351 targetHasChanged = true;
352 break;
353 default:
354 break;
355 }
356 } else {
357 switch( np-nm ) {
358 case 0:
359 if( G4UniformRand() < 0.25 ) {
360 currentParticle.SetDefinitionAndUpdateE( aPiZero );
361 targetParticle.SetDefinitionAndUpdateE( aProton );
362 incidentHasChanged = true;
363 targetHasChanged = true;
364 }
365 break;
366 case 1:
367 currentParticle.SetDefinitionAndUpdateE( aPiZero );
368 incidentHasChanged = true;
369 break;
370 default:
371 targetParticle.SetDefinitionAndUpdateE( aProton );
372 targetHasChanged = true;
373 break;
374 }
375 }
376 SetUpPions( np, nm, nz, vec, vecLen );
377 return;
378 }
379
380 /* end of file */
381
Note: See TracBrowser for help on using the repository browser.