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

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

import all except CVS
File size: 46.5 KB
RevLine 
[819]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: G4RPGFragmentation.cc,v 1.3 2007/12/06 01:13:14 dennis Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29
30#include "G4RPGFragmentation.hh"
31#include "G4AntiProton.hh"
32#include "G4AntiNeutron.hh"
33#include "Randomize.hh"
34#include "G4Poisson.hh"
35#include <iostream>
36#include "G4HadReentrentException.hh"
37#include <signal.h>
38
39
40G4RPGFragmentation::G4RPGFragmentation()
41 : G4RPGReaction() {}
42
43
44void G4RPGFragmentation::
45FragmentationIntegral(G4double pt, G4double et, G4double parMass, G4double secMass)
46{
47 pt = std::max( 0.001, pt );
48 G4double dx = 1./(19.*pt);
49 G4double x;
50 G4double term1;
51 G4double term2;
52
53 for (G4int i = 1; i < 20; i++) {
54 x = (G4double(i) - 0.5)*dx;
55 term1 = 1. + parMass*parMass*x*x;
56 term2 = pt*x*et*pt*x*et + pt*pt + secMass*secMass;
57 dndl[i] = dx / std::sqrt( term1*term1*term1*term2 )
58 + dndl[i-1];
59 }
60}
61
62
63G4bool G4RPGFragmentation::
64ReactionStage(const G4HadProjectile* originalIncident,
65 G4ReactionProduct& modifiedOriginal,
66 G4bool& incidentHasChanged,
67 const G4DynamicParticle* originalTarget,
68 G4ReactionProduct& targetParticle,
69 G4bool& targetHasChanged,
70 const G4Nucleus& targetNucleus,
71 G4ReactionProduct& currentParticle,
72 G4FastVector<G4ReactionProduct,256>& vec,
73 G4int& vecLen,
74 G4bool leadFlag,
75 G4ReactionProduct& leadingStrangeParticle)
76{
77 //
78 // Derived from H. Fesefeldt's original FORTRAN code GENXPT
79 //
80 // Generation of x- and pT- values for incident, target, and all secondary
81 // particles using a simple single variable description E D3S/DP3= F(Q)
82 // with Q^2 = (M*X)^2 + PT^2. Final state kinematics are produced by an
83 // FF-type iterative cascade method
84 //
85 // Internal units are GeV
86 //
87
88 // Protection in case no secondary has been created. In that case use
89 // two-body scattering
90 //
91 if(vecLen == 0) return false;
92
93 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
94 G4ParticleDefinition *aProton = G4Proton::Proton();
95 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
96 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
97 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
98
99 G4int i, l;
100 G4bool veryForward = false;
101
102 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
103 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
104 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
105 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
106 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
107 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
108 targetMass*targetMass +
109 2.0*targetMass*etOriginal ); // GeV
110 G4double currentMass = currentParticle.GetMass()/GeV;
111 targetMass = targetParticle.GetMass()/GeV;
112 //
113 // randomize the order of the secondary particles
114 // note that the current and target particles are not affected
115 //
116 for( i=0; i<vecLen; ++i )
117 {
118 G4int itemp = G4int( G4UniformRand()*vecLen );
119 G4ReactionProduct pTemp = *vec[itemp];
120 *vec[itemp] = *vec[i];
121 *vec[i] = pTemp;
122 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
123 }
124
125 if( currentMass == 0.0 && targetMass == 0.0 )
126 {
127 // Target and projectile have annihilated. Replace them with the first
128 // two secondaries in the list. Current particle KE is maintained.
129
130 G4double ek = currentParticle.GetKineticEnergy();
131 G4ThreeVector m = currentParticle.GetMomentum();
132 currentParticle = *vec[0];
133 targetParticle = *vec[1];
134 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
135 G4ReactionProduct *temp = vec[vecLen-1];
136 delete temp;
137 temp = vec[vecLen-2];
138 delete temp;
139 vecLen -= 2;
140 currentMass = currentParticle.GetMass()/GeV;
141 targetMass = targetParticle.GetMass()/GeV;
142 incidentHasChanged = true;
143 targetHasChanged = true;
144 currentParticle.SetKineticEnergy( ek );
145 currentParticle.SetMomentum( m );
146 veryForward = true;
147 }
148 const G4double atomicWeight = targetNucleus.GetN();
149 const G4double atomicNumber = targetNucleus.GetZ();
150 const G4double protonMass = aProton->GetPDGMass()/MeV;
151
152 if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
153 && G4UniformRand() >= 0.7) {
154 G4ReactionProduct temp = currentParticle;
155 currentParticle = targetParticle;
156 targetParticle = temp;
157 incidentHasChanged = true;
158 targetHasChanged = true;
159 currentMass = currentParticle.GetMass()/GeV;
160 targetMass = targetParticle.GetMass()/GeV;
161 }
162 const G4double afc = std::min( 0.75,
163 0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
164 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
165
166 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
167 G4double forwardEnergy = freeEnergy/2.;
168 G4int forwardCount = 1; // number of particles in forward hemisphere
169
170 G4double backwardEnergy = freeEnergy/2.;
171 G4int backwardCount = 1; // number of particles in backward hemisphere
172
173 if(veryForward)
174 {
175 if(currentParticle.GetSide()==-1)
176 {
177 forwardEnergy += currentMass;
178 forwardCount --;
179 backwardEnergy -= currentMass;
180 backwardCount ++;
181 }
182 if(targetParticle.GetSide()!=-1)
183 {
184 backwardEnergy += targetMass;
185 backwardCount --;
186 forwardEnergy -= targetMass;
187 forwardCount ++;
188 }
189 }
190
191 for( i=0; i<vecLen; ++i )
192 {
193 if( vec[i]->GetSide() == -1 )
194 {
195 ++backwardCount;
196 backwardEnergy -= vec[i]->GetMass()/GeV;
197 } else {
198 ++forwardCount;
199 forwardEnergy -= vec[i]->GetMass()/GeV;
200 }
201 }
202 //
203 // Add particles from intranuclear cascade.
204 // nuclearExcitationCount = number of new secondaries produced by nuclear excitation
205 // extraCount = number of nucleons within these new secondaries
206 //
207 G4double xtarg;
208 if( centerofmassEnergy < (2.0+G4UniformRand()) )
209 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
210 else
211 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
212 if( xtarg <= 0.0 )xtarg = 0.01;
213 G4int nuclearExcitationCount = G4Poisson( xtarg );
214 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
215 G4int extraNucleonCount = 0;
216 G4double extraNucleonMass = 0.0;
217 if( nuclearExcitationCount > 0 )
218 {
219 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
220 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
221 G4int momentumBin = 0;
222 while( (momentumBin < 6) &&
223 (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
224 ++momentumBin;
225 momentumBin = std::min( 5, momentumBin );
226 //
227 // NOTE: in GENXPT, these new particles were given negative codes
228 // here I use NewlyAdded = true instead
229 //
230
231 for( i=0; i<nuclearExcitationCount; ++i )
232 {
233 G4ReactionProduct * pVec = new G4ReactionProduct();
234 if( G4UniformRand() < nucsup[momentumBin] )
235 {
236 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
237 pVec->SetDefinition( aProton );
238 else
239 pVec->SetDefinition( aNeutron );
240 pVec->SetSide( -2 ); // -2 means backside nucleon
241 ++extraNucleonCount;
242 backwardEnergy += pVec->GetMass()/GeV;
243 extraNucleonMass += pVec->GetMass()/GeV;
244 }
245 else
246 {
247 G4double ran = G4UniformRand();
248 if( ran < 0.3181 )
249 pVec->SetDefinition( aPiPlus );
250 else if( ran < 0.6819 )
251 pVec->SetDefinition( aPiZero );
252 else
253 pVec->SetDefinition( aPiMinus );
254 pVec->SetSide( -1 ); // backside particle, but not a nucleon
255 }
256 pVec->SetNewlyAdded( true ); // true is the same as IPA(i)<0
257 vec.SetElement( vecLen++, pVec );
258 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
259 backwardEnergy -= pVec->GetMass()/GeV;
260 ++backwardCount;
261 }
262 }
263 //
264 // assume conservation of kinetic energy in forward & backward hemispheres
265 //
266 G4int is, iskip;
267 while( forwardEnergy <= 0.0 ) // must eliminate a particle from the forward side
268 {
269 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
270 iskip = G4int(G4UniformRand()*forwardCount) + 1; // 1 <= iskip <= forwardCount
271 is = 0;
272 G4int forwardParticlesLeft = 0;
273 for( i=(vecLen-1); i>=0; --i )
274 {
275 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
276 {
277 forwardParticlesLeft = 1;
278 if( ++is == iskip )
279 {
280 forwardEnergy += vec[i]->GetMass()/GeV;
281 for( G4int j=i; j<(vecLen-1); j++ )*vec[j] = *vec[j+1]; // shift up
282 --forwardCount;
283 delete vec[vecLen-1];
284 if( --vecLen == 0 )return false; // all the secondaries have been eliminated
285 break; // --+
286 } // |
287 } // |
288 } // break goes down to here
289 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
290 if( forwardParticlesLeft == 0 )
291 {
292 G4int iremove = -1;
293 for (G4int i = 0; i < vecLen; i++) {
294 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
295 iremove = i;
296 break;
297 }
298 }
299 if (iremove == -1) {
300 for (G4int i = 0; i < vecLen; i++) {
301 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
302 iremove = i;
303 break;
304 }
305 }
306 }
307 if (iremove == -1) iremove = 0;
308
309 forwardEnergy += vec[iremove]->GetMass()/GeV;
310 if (vec[iremove]->GetSide() > 0) --forwardCount;
311
312 for (G4int i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
313 delete vec[vecLen-1];
314 vecLen--;
315 if (vecLen == 0) return false; // all secondaries have been eliminated
316 break;
317 }
318 } // while
319
320 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
321 while( backwardEnergy <= 0.0 ) // must eliminate a particle from the backward side
322 {
323 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
324 iskip = G4int(G4UniformRand()*backwardCount) + 1; // 1 <= iskip <= backwardCount
325 is = 0;
326 G4int backwardParticlesLeft = 0;
327 for( i=(vecLen-1); i>=0; --i )
328 {
329 if( vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled())
330 {
331 backwardParticlesLeft = 1;
332 if( ++is == iskip ) // eliminate the i'th particle
333 {
334 if( vec[i]->GetSide() == -2 )
335 {
336 --extraNucleonCount;
337 extraNucleonMass -= vec[i]->GetMass()/GeV;
338 backwardEnergy -= vec[i]->GetTotalEnergy()/GeV;
339 }
340 backwardEnergy += vec[i]->GetTotalEnergy()/GeV;
341 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
342 --backwardCount;
343 delete vec[vecLen-1];
344 if( --vecLen == 0 )return false; // all the secondaries have been eliminated
345 break;
346 }
347 }
348 }
349 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
350 if( backwardParticlesLeft == 0 )
351 {
352 G4int iremove = -1;
353 for (G4int i = 0; i < vecLen; i++) {
354 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
355 iremove = i;
356 break;
357 }
358 }
359 if (iremove == -1) {
360 for (G4int i = 0; i < vecLen; i++) {
361 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
362 iremove = i;
363 break;
364 }
365 }
366 }
367 if (iremove == -1) iremove = 0;
368
369 backwardEnergy += vec[iremove]->GetMass()/GeV;
370 if (vec[iremove]->GetSide() > 0) --backwardCount;
371
372 for (G4int i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
373 delete vec[vecLen-1];
374 vecLen--;
375 if (vecLen == 0) return false; // all secondaries have been eliminated
376 break;
377 }
378 } // while
379
380 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
381 //
382 // define initial state vectors for Lorentz transformations
383 // the pseudoParticles have non-standard masses, hence the "pseudo"
384 //
385 G4ReactionProduct pseudoParticle[8];
386 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
387
388 pseudoParticle[0].SetMass( mOriginal*GeV );
389 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
390 pseudoParticle[0].SetTotalEnergy(
391 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
392
393 pseudoParticle[1].SetMass( protonMass*MeV ); // this could be targetMass
394 pseudoParticle[1].SetTotalEnergy( protonMass*MeV );
395
396 pseudoParticle[3].SetMass( protonMass*(1+extraNucleonCount)*MeV );
397 pseudoParticle[3].SetTotalEnergy( protonMass*(1+extraNucleonCount)*MeV );
398
399 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
400 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
401
402 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
403 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
404
405 //
406 // main loop for 4-momentum generation
407 // see Pitha-report (Aachen) for a detailed description of the method
408 //
409 G4double aspar, pt, et, x, pp, pp1, wgt;
410 G4int innerCounter, outerCounter;
411 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
412
413 G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
414 //
415 // process the secondary particles in reverse order
416 // the incident particle is Done after the secondaries
417 // nucleons, including the target, in the backward hemisphere are also Done later
418 //
419 G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere
420 G4double totalEnergy, kineticEnergy, vecMass;
421
422 for( i=(vecLen-1); i>=0; --i )
423 {
424 G4double phi = G4UniformRand()*twopi;
425 if( vec[i]->GetNewlyAdded() ) // added from intranuclear cascade
426 {
427 if( vec[i]->GetSide() == -2 ) // is a nucleon
428 {
429 if( backwardNucleonCount < 18 )
430 {
431 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
432 for(G4int i=0; i<vecLen; i++) delete vec[i];
433 vecLen = 0;
434 throw G4HadReentrentException(__FILE__, __LINE__,
435 "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon");
436 }
437 vec[i]->SetSide( -3 );
438 ++backwardNucleonCount;
439 continue;
440 }
441 }
442 }
443 //
444 // set pt and phi values, they are changed somewhat in the iteration loop
445 // set mass parameter for lambda fragmentation model
446 //
447 vecMass = vec[i]->GetMass()/GeV;
448 G4double ran = -std::log(1.0-G4UniformRand())/3.5;
449 if( vec[i]->GetSide() == -2 ) // backward nucleon
450 {
451 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon" ||
452 vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
453 aspar = 0.75;
454 pt = std::sqrt( std::pow( ran, 1.7 ) );
455 } else { // vec[i] must be a proton, neutron,
456 aspar = 0.20; // lambda, sigma, xsi, or ion
457 pt = std::sqrt( std::pow( ran, 1.2 ) );
458 }
459
460 } else { // not a backward nucleon
461
462 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
463 aspar = 0.75;
464 pt = std::sqrt( std::pow( ran, 1.7 ) );
465 } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
466 aspar = 0.70;
467 pt = std::sqrt( std::pow( ran, 1.7 ) );
468 } else { // vec[i] must be a proton, neutron,
469 aspar = 0.65; // lambda, sigma, xsi, or ion
470 pt = std::sqrt( std::pow( ran, 1.5 ) );
471 }
472 }
473 pt = std::max( 0.001, pt );
474 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
475 if( vec[i]->GetSide() > 0 )
476 et = pseudoParticle[0].GetTotalEnergy()/GeV;
477 else
478 et = pseudoParticle[1].GetTotalEnergy()/GeV;
479
480 //
481 // start of outer iteration loop
482 //
483 outerCounter = 0;
484 eliminateThisParticle = true;
485 resetEnergies = true;
486 dndl[0] = 0.0;
487
488 while( ++outerCounter < 3 )
489 {
490 FragmentationIntegral(pt, et, aspar, vecMass);
491
492 innerCounter = 0;
493 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
494 //
495 // start of inner iteration loop
496 //
497 while( ++innerCounter < 7 )
498 {
499 ran = G4UniformRand()*dndl[19];
500 l = 1;
501 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
502 x = (G4double(l-1) + G4UniformRand())/19.;
503 if( vec[i]->GetSide() < 0 )x *= -1.;
504 vec[i]->SetMomentum( x*et*GeV ); // set the z-momentum
505 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
506 vec[i]->SetTotalEnergy( totalEnergy*GeV );
507 kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
508 if( vec[i]->GetSide() > 0 ) // forward side
509 {
510 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy )
511 {
512 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
513 forwardKinetic += kineticEnergy;
514 outerCounter = 2; // leave outer loop
515 eliminateThisParticle = false; // don't eliminate this particle
516 resetEnergies = false;
517 break; // leave inner loop
518 }
519 if( innerCounter > 5 )break; // leave inner loop
520 if( backwardEnergy >= vecMass ) // switch sides
521 {
522 vec[i]->SetSide( -1 );
523 forwardEnergy += vecMass;
524 backwardEnergy -= vecMass;
525 ++backwardCount;
526 }
527 } else { // backward side
528 if( extraNucleonCount > 19 ) // commented out to duplicate ?bug? in GENXPT
529 x = 0.999;
530 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
531 if( (backwardKinetic+kineticEnergy) < xxx*backwardEnergy )
532 {
533 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
534 backwardKinetic += kineticEnergy;
535 outerCounter = 2; // leave outer loop
536 eliminateThisParticle = false; // don't eliminate this particle
537 resetEnergies = false;
538 break; // leave inner loop
539 }
540 if( innerCounter > 5 )break; // leave inner loop
541 if( forwardEnergy >= vecMass ) // switch sides
542 {
543 vec[i]->SetSide( 1 );
544 forwardEnergy -= vecMass;
545 backwardEnergy += vecMass;
546 backwardCount--;
547 }
548 }
549 G4ThreeVector momentum = vec[i]->GetMomentum();
550 vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
551 pt *= 0.9;
552 dndl[19] *= 0.9;
553 } // closes inner loop
554 if( resetEnergies ) {
555 // If we get to here, the inner loop has been done 6 times.
556 // Reset the kinetic energies of previously done particles, if
557 // they are lighter than protons and in the forward hemisphere,
558 // then continue with outer loop.
559 //
560 forwardKinetic = 0.0;
561 backwardKinetic = 0.0;
562 pseudoParticle[4].SetZero();
563 pseudoParticle[5].SetZero();
564 for( l=i+1; l<vecLen; ++l ) {
565 if (vec[l]->GetSide() > 0 ||
566 vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
567 vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
568
569 G4double tempMass = vec[l]->GetMass()/MeV;
570 totalEnergy = 0.95*vec[l]->GetTotalEnergy()/MeV + 0.05*tempMass;
571 totalEnergy = std::max( tempMass, totalEnergy );
572 vec[l]->SetTotalEnergy( totalEnergy*MeV );
573 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) );
574 pp1 = vec[l]->GetMomentum().mag()/MeV;
575 if( pp1 < 1.0e-6*GeV ) {
576 G4ThreeVector iso = Isotropic(pp);
577 vec[l]->SetMomentum( iso.x(), iso.y(), iso.z() );
578 } else {
579 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
580 }
581 G4double px = vec[l]->GetMomentum().x()/MeV;
582 G4double py = vec[l]->GetMomentum().y()/MeV;
583 pt = std::max( 1.0, std::sqrt( px*px + py*py ) )/GeV;
584 if( vec[l]->GetSide() > 0 )
585 {
586 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
587 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
588 } else {
589 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
590 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
591 }
592 } // if pi, K or forward
593 } // for l
594 } // if resetEnergies
595 } // closes outer loop
596
597 if( eliminateThisParticle && vec[i]->GetMayBeKilled()) // not enough energy, eliminate this particle
598 {
599 if( vec[i]->GetSide() > 0 )
600 {
601 --forwardCount;
602 forwardEnergy += vecMass;
603 } else {
604 if( vec[i]->GetSide() == -2 )
605 {
606 --extraNucleonCount;
607 extraNucleonMass -= vecMass;
608 backwardEnergy -= vecMass;
609 }
610 --backwardCount;
611 backwardEnergy += vecMass;
612 }
613 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
614 G4ReactionProduct *temp = vec[vecLen-1];
615 delete temp;
616 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
617 if( --vecLen == 0 )return false; // all the secondaries have been eliminated
618 }
619 } // closes main for loop
620
621 //
622 // for the incident particle: it was placed in the forward hemisphere
623 // set pt and phi values, they are changed somewhat in the iteration loop
624 // set mass parameter for lambda fragmentation model
625 //
626 G4double phi = G4UniformRand()*twopi;
627 G4double ran = -std::log(1.0-G4UniformRand());
628 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
629 aspar = 0.60;
630 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
631 } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
632 aspar = 0.50;
633 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
634 } else {
635 aspar = 0.40;
636 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
637 }
638
639 currentParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
640 et = pseudoParticle[0].GetTotalEnergy()/GeV;
641 dndl[0] = 0.0;
642 vecMass = currentParticle.GetMass()/GeV;
643
644 FragmentationIntegral(pt, et, aspar, vecMass);
645
646 ran = G4UniformRand()*dndl[19];
647 l = 1;
648 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
649 x = (G4double(l-1) + G4UniformRand())/19.;
650 currentParticle.SetMomentum( x*et*GeV ); // set the z-momentum
651 if( forwardEnergy < forwardKinetic )
652 totalEnergy = vecMass + 0.04*std::fabs(normal());
653 else
654 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
655 currentParticle.SetTotalEnergy( totalEnergy*GeV );
656 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
657 pp1 = currentParticle.GetMomentum().mag()/MeV;
658
659 if( pp1 < 1.0e-6*GeV ) {
660 G4ThreeVector iso = Isotropic(pp);
661 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
662 } else {
663 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
664 }
665 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
666
667 //
668 // Current particle now finished
669 //
670 // Begin target particle
671 //
672
673 if( backwardNucleonCount < 18 )
674 {
675 targetParticle.SetSide( -3 );
676 ++backwardNucleonCount;
677 }
678 else
679 {
680 // set pt and phi values, they are changed somewhat in the iteration loop
681 // set mass parameter for lambda fragmentation model
682 //
683 vecMass = targetParticle.GetMass()/GeV;
684 ran = -std::log(1.0-G4UniformRand());
685 aspar = 0.40;
686 pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
687 targetParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
688 et = pseudoParticle[1].GetTotalEnergy()/GeV;
689 outerCounter = 0;
690 eliminateThisParticle = true; // should never eliminate the target particle
691 resetEnergies = true;
692 dndl[0] = 0.0;
693
694 while( ++outerCounter < 3 ) // start of outer iteration loop
695 {
696 FragmentationIntegral(pt, et, aspar, vecMass);
697
698 innerCounter = 0;
699 while( ++innerCounter < 7 ) // start of inner iteration loop
700 {
701 ran = G4UniformRand()*dndl[19];
702 l = 1;
703 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
704 x = (G4double(l-1) + G4UniformRand())/19.;
705 if( targetParticle.GetSide() < 0 )x *= -1.;
706 targetParticle.SetMomentum( x*et*GeV ); // set the z-momentum
707 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
708 targetParticle.SetTotalEnergy( totalEnergy*GeV );
709 if( targetParticle.GetSide() < 0 )
710 {
711 if( extraNucleonCount > 19 )x=0.999;
712 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
713 if( (backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy )
714 {
715 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
716 backwardKinetic += totalEnergy - vecMass;
717 // pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
718 // pseudoParticle[6].SetMomentum( 0.0 ); // set z-momentum
719 outerCounter = 2; // leave outer loop
720 eliminateThisParticle = false; // don't eliminate this particle
721 resetEnergies = false;
722 break; // leave inner loop
723 }
724 if( innerCounter > 5 )break; // leave inner loop
725 if( forwardEnergy >= vecMass ) // switch sides
726 {
727 targetParticle.SetSide( 1 );
728 forwardEnergy -= vecMass;
729 backwardEnergy += vecMass;
730 --backwardCount;
731 }
732 G4ThreeVector momentum = targetParticle.GetMomentum();
733 targetParticle.SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
734 pt *= 0.9;
735 dndl[19] *= 0.9;
736 }
737 else // target has gone to forward side
738 {
739 if( forwardEnergy < forwardKinetic )
740 totalEnergy = vecMass + 0.04*std::fabs(normal());
741 else
742 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
743 targetParticle.SetTotalEnergy( totalEnergy*GeV );
744 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
745 pp1 = targetParticle.GetMomentum().mag()/MeV;
746 if( pp1 < 1.0e-6*GeV ) {
747 G4ThreeVector iso = Isotropic(pp);
748 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
749 } else {
750 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
751 }
752
753 pseudoParticle[4] = pseudoParticle[4] + targetParticle;
754 outerCounter = 2; // leave outer loop
755 eliminateThisParticle = false; // don't eliminate this particle
756 resetEnergies = false;
757 break; // leave inner loop
758 }
759 } // closes inner loop
760
761 if( resetEnergies ) {
762 // If we get to here, the inner loop has been done 6 times.
763 // Reset the kinetic energies of previously done particles,
764 // if they are lighter than protons and in the forward hemisphere,
765 // then continue with outer loop.
766
767 forwardKinetic = backwardKinetic = 0.0;
768 pseudoParticle[4].SetZero();
769 pseudoParticle[5].SetZero();
770 for( l=0; l<vecLen; ++l ) {
771 if (vec[l]->GetSide() > 0 ||
772 vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
773 vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
774 G4double tempMass = vec[l]->GetMass()/GeV;
775 totalEnergy =
776 std::max( tempMass, 0.95*vec[l]->GetTotalEnergy()/GeV + 0.05*tempMass );
777 vec[l]->SetTotalEnergy( totalEnergy*GeV );
778 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) )*GeV;
779 pp1 = vec[l]->GetMomentum().mag()/MeV;
780 if( pp1 < 1.0e-6*GeV ) {
781 G4ThreeVector iso = Isotropic(pp);
782 vec[l]->SetMomentum( iso.x(), iso.y(), iso.z() );
783 } else {
784 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
785 }
786 pt = std::max( 0.001*GeV, std::sqrt( sqr(vec[l]->GetMomentum().x()/MeV) +
787 sqr(vec[l]->GetMomentum().y()/MeV) ) )/GeV;
788 if( vec[l]->GetSide() > 0)
789 {
790 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
791 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
792 } else {
793 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
794 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
795 }
796 } // if pi, K or forward
797 } // for l
798 } // if (resetEnergies)
799 } // closes outer loop
800
801// if( eliminateThisParticle ) // not enough energy, eliminate target
802// {
803// G4cerr << "Warning: eliminating target particle" << G4endl;
804// exit( EXIT_FAILURE );
805// }
806 }
807 //
808 // Target particle finished.
809 //
810 // Now produce backward nucleons with a cluster model
811 //
812 pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
813 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
814 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
815 if( backwardNucleonCount == 1 ) // target particle is the only backward nucleon
816 {
817 G4double ekin =
818 std::min( backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV );
819
820 if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
821 vecMass = targetParticle.GetMass()/GeV;
822 totalEnergy = ekin+vecMass;
823 targetParticle.SetTotalEnergy( totalEnergy*GeV );
824 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
825 pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
826 if( pp1 < 1.0e-6*GeV ) {
827 G4ThreeVector iso = Isotropic(pp);
828 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
829 } else {
830 targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1) );
831 }
832 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
833 }
834 else if (backwardNucleonCount > 1)
835 {
836 const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
837 const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
838
839 G4int tempCount = 5;
840 if (backwardNucleonCount < 5) tempCount = backwardNucleonCount;
841 tempCount -= 2;
842
843 G4double rmb = 0.;
844 if( targetParticle.GetSide() == -3 ) rmb += targetParticle.GetMass()/GeV;
845 for( i=0; i<vecLen; ++i )
846 {
847 if( vec[i]->GetSide() == -3 ) rmb += vec[i]->GetMass()/GeV;
848 }
849 rmb += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[tempCount]) / gpar[tempCount];
850 totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
851 vecMass = std::min( rmb, totalEnergy );
852 pseudoParticle[6].SetMass( vecMass*GeV );
853 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
854 pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
855 if( pp1 < 1.0e-6*GeV ) {
856 G4ThreeVector iso = Isotropic(pp);
857 pseudoParticle[6].SetMomentum( iso.x(), iso.y(), iso.z() );
858 } else {
859 pseudoParticle[6].SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
860 }
861 G4FastVector<G4ReactionProduct,256> tempV; // tempV contains the backward nucleons
862 tempV.Initialize( backwardNucleonCount );
863 G4int tempLen = 0;
864 if( targetParticle.GetSide() == -3 )tempV.SetElement( tempLen++, &targetParticle );
865 for( i=0; i<vecLen; ++i )
866 {
867 if( vec[i]->GetSide() == -3 )tempV.SetElement( tempLen++, vec[i] );
868 }
869 if( tempLen != backwardNucleonCount )
870 {
871 G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
872 G4cerr << "tempLen = " << tempLen;
873 G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
874 G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
875 G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
876 for( i=0; i<vecLen; ++i )
877 G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
878 exit( EXIT_FAILURE );
879 }
880 constantCrossSection = true;
881 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
882 if( tempLen >= 2 )
883 {
884 wgt = GenerateNBodyEvent(
885 pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen );
886 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
887 if( targetParticle.GetSide() == -3 )
888 {
889 targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
890 // tempV contains the real stuff
891 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
892 }
893 for( i=0; i<vecLen; ++i )
894 {
895 if( vec[i]->GetSide() == -3 )
896 {
897 vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
898 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
899 }
900 }
901 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
902 }
903 }
904 else return false;
905
906 //
907 // Lorentz transformation in lab system
908 //
909 if( vecLen == 0 )return false; // all the secondaries have been eliminated
910 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
911
912 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
913 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
914 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[1] );
915
916 // leadFlag will be true if original particle and incident particle are
917 // both strange, in which case the incident particle becomes the leading
918 // particle.
919 // leadFlag will also be true if the target particle is strange, but the
920 // incident particle is not, in which case the target particle becomes the
921 // leading particle.
922
923 G4bool leadingStrangeParticleHasChanged = true;
924 if( leadFlag )
925 {
926 if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
927 leadingStrangeParticleHasChanged = false;
928 if( leadingStrangeParticleHasChanged &&
929 ( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() ) )
930 leadingStrangeParticleHasChanged = false;
931 if( leadingStrangeParticleHasChanged )
932 {
933 for( i=0; i<vecLen; i++ )
934 {
935 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
936 {
937 leadingStrangeParticleHasChanged = false;
938 break;
939 }
940 }
941 }
942 if( leadingStrangeParticleHasChanged )
943 {
944 G4bool leadTest =
945 (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
946 leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
947 G4bool targetTest =
948 (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
949 targetParticle.GetDefinition()->GetParticleSubType() == "pi");
950
951 // following modified by JLC 22-Oct-97
952
953 if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
954 {
955 targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
956 targetHasChanged = true;
957 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
958 }
959 else
960 {
961 currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
962 incidentHasChanged = false;
963 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
964 }
965 }
966 } // end of if( leadFlag )
967
968 // Get number of final state nucleons and nucleons remaining in
969 // target nucleus
970
971 std::pair<G4int, G4int> finalStateNucleons =
972 GetFinalStateNucleons(originalTarget, vec, vecLen);
973
974 G4int protonsInFinalState = finalStateNucleons.first;
975 G4int neutronsInFinalState = finalStateNucleons.second;
976
977 G4int numberofFinalStateNucleons =
978 protonsInFinalState + neutronsInFinalState;
979
980 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
981 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
982 originalIncident->GetDefinition()->GetPDGMass() <
983 G4Lambda::Lambda()->GetPDGMass())
984 numberofFinalStateNucleons++;
985
986 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
987
988 G4int PinNucleus = std::max(0,
989 G4int(targetNucleus.GetZ()) - protonsInFinalState);
990 G4int NinNucleus = std::max(0,
991 G4int(targetNucleus.GetN()-targetNucleus.GetZ()) - neutronsInFinalState);
992
993 pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
994 pseudoParticle[3].SetMass( mOriginal*GeV );
995 pseudoParticle[3].SetTotalEnergy(
996 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
997
998 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
999 G4int diff = 0;
1000 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
1001 if(numberofFinalStateNucleons == 1) diff = 0;
1002 pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
1003 pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV );
1004 pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV );
1005
1006 G4double theoreticalKinetic =
1007 pseudoParticle[3].GetTotalEnergy()/MeV +
1008 pseudoParticle[4].GetTotalEnergy()/MeV -
1009 currentParticle.GetMass()/MeV -
1010 targetParticle.GetMass()/MeV;
1011
1012 G4double simulatedKinetic =
1013 currentParticle.GetKineticEnergy()/MeV + targetParticle.GetKineticEnergy()/MeV;
1014
1015 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
1016 pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
1017 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
1018
1019 pseudoParticle[7].SetZero();
1020 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
1021 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
1022
1023 for( i=0; i<vecLen; ++i )
1024 {
1025 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
1026 simulatedKinetic += vec[i]->GetKineticEnergy()/MeV;
1027 theoreticalKinetic -= vec[i]->GetMass()/MeV;
1028 }
1029
1030 if( vecLen <= 16 && vecLen > 0 )
1031 {
1032 // must create a new set of ReactionProducts here because GenerateNBody will
1033 // modify the momenta for the particles, and we don't want to do this
1034 //
1035 G4ReactionProduct tempR[130];
1036 tempR[0] = currentParticle;
1037 tempR[1] = targetParticle;
1038 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
1039 G4FastVector<G4ReactionProduct,256> tempV;
1040 tempV.Initialize( vecLen+2 );
1041 G4int tempLen = 0;
1042 for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
1043 constantCrossSection = true;
1044
1045 wgt = GenerateNBodyEvent( pseudoParticle[3].GetTotalEnergy()/MeV+
1046 pseudoParticle[4].GetTotalEnergy()/MeV,
1047 constantCrossSection, tempV, tempLen );
1048 if (wgt == -1) {
1049 G4double Qvalue = 0;
1050 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
1051 wgt = GenerateNBodyEvent( Qvalue/MeV,
1052 constantCrossSection, tempV, tempLen );
1053 }
1054 if(wgt>-.5)
1055 {
1056 theoreticalKinetic = 0.0;
1057 for( i=0; i<tempLen; ++i )
1058 {
1059 pseudoParticle[6].Lorentz( *tempV[i], pseudoParticle[4] );
1060 theoreticalKinetic += pseudoParticle[6].GetKineticEnergy()/MeV;
1061 }
1062 }
1063 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1064 }
1065 //
1066 // Make sure, that the kinetic energies are correct
1067 //
1068 if( simulatedKinetic != 0.0 )
1069 {
1070 wgt = (theoreticalKinetic)/simulatedKinetic;
1071 theoreticalKinetic = currentParticle.GetKineticEnergy()/MeV * wgt;
1072 simulatedKinetic = theoreticalKinetic;
1073 currentParticle.SetKineticEnergy( theoreticalKinetic*MeV );
1074 pp = currentParticle.GetTotalMomentum()/MeV;
1075 pp1 = currentParticle.GetMomentum().mag()/MeV;
1076 if( pp1 < 1.0e-6*GeV ) {
1077 G4ThreeVector iso = Isotropic(pp);
1078 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
1079 } else {
1080 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
1081 }
1082 theoreticalKinetic = targetParticle.GetKineticEnergy()/MeV * wgt;
1083 targetParticle.SetKineticEnergy( theoreticalKinetic*MeV );
1084 simulatedKinetic += theoreticalKinetic;
1085 pp = targetParticle.GetTotalMomentum()/MeV;
1086 pp1 = targetParticle.GetMomentum().mag()/MeV;
1087 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1088 if( pp1 < 1.0e-6*GeV ) {
1089 G4ThreeVector iso = Isotropic(pp);
1090 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
1091 } else {
1092 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
1093 }
1094
1095 for( i=0; i<vecLen; ++i ) {
1096 theoreticalKinetic = vec[i]->GetKineticEnergy()/MeV * wgt;
1097 simulatedKinetic += theoreticalKinetic;
1098 vec[i]->SetKineticEnergy( theoreticalKinetic*MeV );
1099 pp = vec[i]->GetTotalMomentum()/MeV;
1100 pp1 = vec[i]->GetMomentum().mag()/MeV;
1101 if( pp1 < 1.0e-6*GeV ) {
1102 G4ThreeVector iso = Isotropic(pp);
1103 vec[i]->SetMomentum( iso.x(), iso.y(), iso.z() );
1104 } else {
1105 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
1106 }
1107 }
1108 }
1109 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1110
1111 Rotate( numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
1112 modifiedOriginal, originalIncident, targetNucleus,
1113 currentParticle, targetParticle, vec, vecLen );
1114 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1115 //
1116 // add black track particles
1117 // the total number of particles produced is restricted to 198
1118 // this may have influence on very high energies
1119 //
1120 if( atomicWeight >= 1.5 )
1121 {
1122 // npnb is number of proton/neutron black track particles
1123 // ndta is the number of deuterons, tritons, and alphas produced
1124 // epnb is the kinetic energy available for proton/neutron black track particles
1125 // edta is the kinetic energy available for deuteron/triton/alpha particles
1126 //
1127 G4int npnb = 0;
1128 G4int ndta = 0;
1129
1130 G4double epnb, edta;
1131 if (veryForward) {
1132 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
1133 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
1134 } else {
1135 epnb = targetNucleus.GetPNBlackTrackEnergy();
1136 edta = targetNucleus.GetDTABlackTrackEnergy();
1137 }
1138
1139 const G4double pnCutOff = 0.001;
1140 const G4double dtaCutOff = 0.001;
1141 const G4double kineticMinimum = 1.e-6;
1142 const G4double kineticFactor = -0.010;
1143 G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules
1144 const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
1145 if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
1146 if( epnb >= pnCutOff )
1147 {
1148 npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
1149 if( numberofFinalStateNucleons + npnb > atomicWeight )
1150 npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
1151 npnb = std::min( npnb, 127-vecLen );
1152 }
1153 if( edta >= dtaCutOff )
1154 {
1155 ndta = G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
1156 ndta = std::min( ndta, 127-vecLen );
1157 }
1158 if (npnb == 0 && ndta == 0) npnb = 1;
1159
1160 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
1161
1162 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
1163 kineticFactor, modifiedOriginal,
1164 PinNucleus, NinNucleus, targetNucleus,
1165 vec, vecLen);
1166 }
1167 // if( centerofmassEnergy <= (4.0+G4UniformRand()) )
1168 // MomentumCheck( modifiedOriginal, currentParticle, targetParticle,
1169 // vec, vecLen );
1170 //
1171 // calculate time delay for nuclear reactions
1172 //
1173 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
1174 currentParticle.SetTOF(
1175 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
1176 else
1177 currentParticle.SetTOF( 1.0 );
1178 return true;
1179
1180}
1181
1182 /* end of file */
Note: See TracBrowser for help on using the repository browser.