source: trunk/source/processes/hadronic/models/rpg/src/G4RPGStrangeProduction.cc@ 1347

Last change on this file since 1347 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 15.5 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: G4RPGStrangeProduction.cc,v 1.1 2007/07/18 21:04:21 dennis Exp $
27// GEANT4 tag $Name: geant4-09-03-ref-09 $
28//
29
30#include "G4RPGStrangeProduction.hh"
31// #include "G4AntiProton.hh"
32// #include "G4AntiNeutron.hh"
33#include "Randomize.hh"
34#include <iostream>
35#include "G4HadReentrentException.hh"
36#include <signal.h>
37
38
39G4RPGStrangeProduction::G4RPGStrangeProduction()
40 : G4RPGReaction() {}
41
42
43G4bool G4RPGStrangeProduction::
44ReactionStage(const G4HadProjectile* /*originalIncident*/,
45 G4ReactionProduct& modifiedOriginal,
46 G4bool& incidentHasChanged,
47 const G4DynamicParticle* originalTarget,
48 G4ReactionProduct& targetParticle,
49 G4bool& targetHasChanged,
50 const G4Nucleus& /*targetNucleus*/,
51 G4ReactionProduct& currentParticle,
52 G4FastVector<G4ReactionProduct,256>& vec,
53 G4int& vecLen,
54 G4bool /*leadFlag*/,
55 G4ReactionProduct& /*leadingStrangeParticle*/)
56{
57 // Derived from H. Fesefeldt's original FORTRAN code STPAIR
58 //
59 // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
60 // K+ Y0, K0 Y+, K0 Y-
61 // For antibaryon induced reactions half of the cross sections KB YB
62 // pairs are produced. Charge is not conserved, no experimental data available
63 // for exclusive reactions, therefore some average behaviour assumed.
64 // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
65 //
66
67 if( vecLen == 0 )return true;
68 //
69 // the following protects against annihilation processes
70 //
71 if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return true;
72
73 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
74 const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
75 G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
76 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
77 targetMass*targetMass +
78 2.0*targetMass*etOriginal ); // GeV
79 G4double currentMass = currentParticle.GetMass()/GeV;
80 G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
81 if( availableEnergy <= 1.0 )return true;
82
83 G4ParticleDefinition *aProton = G4Proton::Proton();
84 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
85 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
86 G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
87 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
88 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
89 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
90 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
91 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
92 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
93 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
94 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
95 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
96 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
97 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
98 G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
99
100 const G4double protonMass = aProton->GetPDGMass()/GeV;
101 const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
102 //
103 // determine the center of mass energy bin
104 //
105 const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
106
107 G4int ibin, i3, i4;
108 G4double avk, avy, avn, ran;
109 G4int i = 1;
110 while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
111 if( i == 12 )
112 ibin = 11;
113 else
114 ibin = i;
115 //
116 // the fortran code chooses a random replacement of produced kaons
117 // but does not take into account charge conservation
118 //
119 if( vecLen == 1 ) // we know that vecLen > 0
120 {
121 i3 = 0;
122 i4 = 1; // note that we will be adding a new secondary particle in this case only
123 }
124 else // otherwise 0 <= i3,i4 < vecLen
125 {
126 G4double ran = G4UniformRand();
127 while( ran == 1.0 )ran = G4UniformRand();
128 i4 = i3 = G4int( vecLen*ran );
129 while( i3 == i4 )
130 {
131 ran = G4UniformRand();
132 while( ran == 1.0 )ran = G4UniformRand();
133 i4 = G4int( vecLen*ran );
134 }
135 }
136
137 //
138 // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
139 //
140 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
141 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
142 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
143 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
144 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
145 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
146
147 avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
148 /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
149 avk = std::exp(avk);
150
151 avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
152 /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
153 avy = std::exp(avy);
154
155 avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
156 /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
157 avn = std::exp(avn);
158
159 if( avk+avy+avn <= 0.0 )return true;
160
161 if( currentMass < protonMass )avy /= 2.0;
162 if( targetMass < protonMass )avy = 0.0;
163 avy += avk+avn;
164 avk += avn;
165 ran = G4UniformRand();
166 if( ran < avn )
167 {
168 if( availableEnergy < 2.0 )return true;
169 if( vecLen == 1 ) // add a new secondary
170 {
171 G4ReactionProduct *p1 = new G4ReactionProduct;
172 if( G4UniformRand() < 0.5 )
173 {
174 vec[0]->SetDefinition( aNeutron );
175 p1->SetDefinition( anAntiNeutron );
176 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
177 vec[0]->SetMayBeKilled(false);
178 p1->SetMayBeKilled(false);
179 }
180 else
181 {
182 vec[0]->SetDefinition( aProton );
183 p1->SetDefinition( anAntiProton );
184 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
185 vec[0]->SetMayBeKilled(false);
186 p1->SetMayBeKilled(false);
187 }
188 vec.SetElement( vecLen++, p1 );
189 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
190 }
191 else
192 { // replace two secondaries
193 if( G4UniformRand() < 0.5 )
194 {
195 vec[i3]->SetDefinition( aNeutron );
196 vec[i4]->SetDefinition( anAntiNeutron );
197 vec[i3]->SetMayBeKilled(false);
198 vec[i4]->SetMayBeKilled(false);
199 }
200 else
201 {
202 vec[i3]->SetDefinition( aProton );
203 vec[i4]->SetDefinition( anAntiProton );
204 vec[i3]->SetMayBeKilled(false);
205 vec[i4]->SetMayBeKilled(false);
206 }
207 }
208 }
209 else if( ran < avk )
210 {
211 if( availableEnergy < 1.0 )return true;
212
213 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
214 0.6875, 0.7500, 0.8750, 1.000 };
215 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
216 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
217 ran = G4UniformRand();
218 i = 0;
219 while( (i<9) && (ran>=kkb[i]) )++i;
220 if( i == 9 )return true;
221 //
222 // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
223 // charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 -
224 //
225 switch( ipakkb1[i] )
226 {
227 case 10:
228 vec[i3]->SetDefinition( aKaonPlus );
229 vec[i3]->SetMayBeKilled(false);
230 break;
231 case 11:
232 vec[i3]->SetDefinition( aKaonZS );
233 vec[i3]->SetMayBeKilled(false);
234 break;
235 case 12:
236 vec[i3]->SetDefinition( aKaonZL );
237 vec[i3]->SetMayBeKilled(false);
238 break;
239 }
240
241 if( vecLen == 1 ) // add a secondary
242 {
243 G4ReactionProduct *p1 = new G4ReactionProduct;
244 switch( ipakkb2[i] )
245 {
246 case 11:
247 p1->SetDefinition( aKaonZS );
248 p1->SetMayBeKilled(false);
249 break;
250 case 12:
251 p1->SetDefinition( aKaonZL );
252 p1->SetMayBeKilled(false);
253 break;
254 case 13:
255 p1->SetDefinition( aKaonMinus );
256 p1->SetMayBeKilled(false);
257 break;
258 }
259 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
260 vec.SetElement( vecLen++, p1 );
261
262 }
263 else // replace
264 {
265 switch( ipakkb2[i] )
266 {
267 case 11:
268 vec[i4]->SetDefinition( aKaonZS );
269 vec[i4]->SetMayBeKilled(false);
270 break;
271 case 12:
272 vec[i4]->SetDefinition( aKaonZL );
273 vec[i4]->SetMayBeKilled(false);
274 break;
275 case 13:
276 vec[i4]->SetDefinition( aKaonMinus );
277 vec[i4]->SetMayBeKilled(false);
278 break;
279 }
280 }
281 }
282 else if( ran < avy )
283 {
284 if( availableEnergy < 1.6 )return true;
285
286 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
287 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
288 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
289 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
290 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
291 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
292 ran = G4UniformRand();
293 i = 0;
294 while( (i<12) && (ran>ky[i]) )++i;
295 if( i == 12 )return true;
296 if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
297 {
298 // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
299 // 0 + 0 0 0 0 + + + 0 + 0
300 //
301 // 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
302 // 0 + 0 0 0 0 - + - 0 - 0
303 switch( ipaky1[i] )
304 {
305 case 18:
306 targetParticle.SetDefinition( aLambda );
307 break;
308 case 20:
309 targetParticle.SetDefinition( aSigmaPlus );
310 break;
311 case 21:
312 targetParticle.SetDefinition( aSigmaZero );
313 break;
314 case 22:
315 targetParticle.SetDefinition( aSigmaMinus );
316 break;
317 }
318 targetHasChanged = true;
319 switch( ipaky2[i] )
320 {
321 case 10:
322 vec[i3]->SetDefinition( aKaonPlus );
323 vec[i3]->SetMayBeKilled(false);
324 break;
325 case 11:
326 vec[i3]->SetDefinition( aKaonZS );
327 vec[i3]->SetMayBeKilled(false);
328 break;
329 case 12:
330 vec[i3]->SetDefinition( aKaonZL );
331 vec[i3]->SetMayBeKilled(false);
332 break;
333 }
334 }
335 else // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
336 {
337 // ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
338 // 24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
339 if( (currentParticle.GetDefinition() == anAntiProton) ||
340 (currentParticle.GetDefinition() == anAntiNeutron) ||
341 (currentParticle.GetDefinition() == anAntiLambda) ||
342 (currentMass > sigmaMinusMass) )
343 {
344 switch( ipakyb1[i] )
345 {
346 case 19:
347 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
348 break;
349 case 23:
350 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
351 break;
352 case 24:
353 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
354 break;
355 case 25:
356 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
357 break;
358 }
359 incidentHasChanged = true;
360 switch( ipakyb2[i] )
361 {
362 case 11:
363 vec[i3]->SetDefinition( aKaonZS );
364 vec[i3]->SetMayBeKilled(false);
365 break;
366 case 12:
367 vec[i3]->SetDefinition( aKaonZL );
368 vec[i3]->SetMayBeKilled(false);
369 break;
370 case 13:
371 vec[i3]->SetDefinition( aKaonMinus );
372 vec[i3]->SetMayBeKilled(false);
373 break;
374 }
375 }
376 else
377 {
378 switch( ipaky1[i] )
379 {
380 case 18:
381 currentParticle.SetDefinitionAndUpdateE( aLambda );
382 break;
383 case 20:
384 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
385 break;
386 case 21:
387 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
388 break;
389 case 22:
390 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
391 break;
392 }
393 incidentHasChanged = true;
394 switch( ipaky2[i] )
395 {
396 case 10:
397 vec[i3]->SetDefinition( aKaonPlus );
398 vec[i3]->SetMayBeKilled(false);
399 break;
400 case 11:
401 vec[i3]->SetDefinition( aKaonZS );
402 vec[i3]->SetMayBeKilled(false);
403 break;
404 case 12:
405 vec[i3]->SetDefinition( aKaonZL );
406 vec[i3]->SetMayBeKilled(false);
407 break;
408 }
409 }
410 }
411 }
412 else return true;
413
414 //
415 // check the available energy
416 // if there is not enough energy for kkb/ky pair production
417 // then reduce the number of secondary particles
418 // NOTE:
419 // the number of secondaries may have been changed
420 // the incident and/or target particles may have changed
421 // charge conservation is ignored (as well as strangness conservation)
422 //
423 currentMass = currentParticle.GetMass()/GeV;
424 targetMass = targetParticle.GetMass()/GeV;
425
426 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
427 for( i=0; i<vecLen; ++i )
428 {
429 energyCheck -= vec[i]->GetMass()/GeV;
430 if( energyCheck < 0.0 ) // chop off the secondary List
431 {
432 vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
433 G4int j;
434 for(j=i; j<vecLen; j++) delete vec[j];
435 break;
436 }
437 }
438
439 return true;
440}
441
442
443 /* end of file */
Note: See TracBrowser for help on using the repository browser.