source: trunk/source/processes/hadronic/models/rpg/src/G4RPGPiPlusInelastic.cc@ 819

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

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