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

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

update ti head
File size: 29.9 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//
[962]26// $Id: G4RPGTwoCluster.cc,v 1.5 2008/06/09 18:13:35 dennis Exp $
[1340]27// GEANT4 tag $Name: geant4-09-03-ref-09 $
[819]28//
29
30#include "G4RPGTwoCluster.hh"
31#include "Randomize.hh"
32#include "G4Poisson.hh"
33#include <iostream>
34#include "G4HadReentrentException.hh"
35#include <signal.h>
36
37
38G4RPGTwoCluster::G4RPGTwoCluster()
39 : G4RPGReaction() {}
40
41
42G4bool G4RPGTwoCluster::
43ReactionStage(const G4HadProjectile* originalIncident,
44 G4ReactionProduct& modifiedOriginal,
45 G4bool& incidentHasChanged,
46 const G4DynamicParticle* originalTarget,
47 G4ReactionProduct& targetParticle,
48 G4bool& targetHasChanged,
49 const G4Nucleus& targetNucleus,
50 G4ReactionProduct& currentParticle,
51 G4FastVector<G4ReactionProduct,256>& vec,
52 G4int& vecLen,
53 G4bool leadFlag,
54 G4ReactionProduct& leadingStrangeParticle)
55{
56 // Derived from H. Fesefeldt's FORTRAN code TWOCLU
57 //
58 // A simple two cluster model is used to generate x- and pt- values for
59 // incident, target, and all secondary particles.
60 // This should be sufficient for low energy interactions.
61 //
62
63 G4int i;
64 G4ParticleDefinition *aProton = G4Proton::Proton();
65 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
66 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
67 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
68 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
69 G4bool veryForward = false;
70
71 const G4double protonMass = aProton->GetPDGMass()/MeV;
72 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
73 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
74 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
75 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
76 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
77 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
78 targetMass*targetMass +
79 2.0*targetMass*etOriginal ); // GeV
80 G4double currentMass = currentParticle.GetMass()/GeV;
81 targetMass = targetParticle.GetMass()/GeV;
82
83 if( currentMass == 0.0 && targetMass == 0.0 )
84 {
85 G4double ek = currentParticle.GetKineticEnergy();
86 G4ThreeVector m = currentParticle.GetMomentum();
87 currentParticle = *vec[0];
88 targetParticle = *vec[1];
89 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
90 if(vecLen<2)
91 {
92 for(G4int i=0; i<vecLen; i++) delete vec[i];
93 vecLen = 0;
94 throw G4HadReentrentException(__FILE__, __LINE__,
95 "G4RPGTwoCluster::ReactionStage : Negative number of particles");
96 }
97 delete vec[vecLen-1];
98 delete vec[vecLen-2];
99 vecLen -= 2;
100 currentMass = currentParticle.GetMass()/GeV;
101 targetMass = targetParticle.GetMass()/GeV;
102 incidentHasChanged = true;
103 targetHasChanged = true;
104 currentParticle.SetKineticEnergy( ek );
105 currentParticle.SetMomentum( m );
106 veryForward = true;
107 }
108
109 const G4double atomicWeight = targetNucleus.GetN();
110 const G4double atomicNumber = targetNucleus.GetZ();
111 //
112 // particles have been distributed in forward and backward hemispheres
113 // in center of mass system of the hadron nucleon interaction
114 //
115
116 // Incident particle always in forward hemisphere
117
118 G4int forwardCount = 1; // number of particles in forward hemisphere
119 currentParticle.SetSide( 1 );
120 G4double forwardMass = currentParticle.GetMass()/GeV;
121 G4double cMass = forwardMass;
122
123 // Target particle always in backward hemisphere
124
125 G4int backwardCount = 1; // number of particles in backward hemisphere
126 targetParticle.SetSide( -1 );
127 G4double backwardMass = targetParticle.GetMass()/GeV;
128 G4double bMass = backwardMass;
129
[962]130 // G4int backwardNucleonCount = 1; // number of nucleons in backward hemisphere
[819]131
132 for( i=0; i<vecLen; ++i )
133 {
134 if( vec[i]->GetSide() < 0 )vec[i]->SetSide( -1 ); // to take care of
135 // case where vec has been preprocessed by GenerateXandPt
136 // and some of them have been set to -2 or -3
137 if( vec[i]->GetSide() == -1 )
138 {
139 ++backwardCount;
140 backwardMass += vec[i]->GetMass()/GeV;
141 }
142 else
143 {
144 ++forwardCount;
145 forwardMass += vec[i]->GetMass()/GeV;
146 }
147 }
148
149 //
150 // Add nucleons and some pions from intra-nuclear cascade
151 //
152
153 G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
154 if(term1 < 0) term1 = 0.0001; // making sure xtarg<0;
155 const G4double afc = 0.312 + 0.2 * std::log(term1);
156 G4double xtarg;
157 if( centerofmassEnergy < 2.0+G4UniformRand() ) // added +2 below, JLC 4Jul97
158 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
159 else
160 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
161 if( xtarg <= 0.0 )xtarg = 0.01;
162 G4int nuclearExcitationCount = G4Poisson( xtarg );
163
164 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
[962]165 // G4int extraNucleonCount = 0;
166 // G4double extraMass = 0.0;
167 // G4double extraNucleonMass = 0.0;
[819]168 if( nuclearExcitationCount > 0 )
169 {
170 G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );
171 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
172 //
173 // NOTE: in TWOCLU, these new particles were given negative codes
174 // here we use NewlyAdded = true instead
175 //
176 for( i=0; i<nuclearExcitationCount; ++i )
177 {
178 G4ReactionProduct* pVec = new G4ReactionProduct();
179 if( G4UniformRand() < nucsup[momentumBin] ) // add proton or neutron
180 {
181 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
182 pVec->SetDefinition( aProton );
183 else
184 pVec->SetDefinition( aNeutron );
[962]185 // Not used ++backwardNucleonCount;
186 // Not used ++extraNucleonCount;
187 // Not used extraNucleonMass += pVec->GetMass()/GeV;
[819]188 }
189 else
190 { // add a pion
191 G4double ran = G4UniformRand();
192 if( ran < 0.3181 )
193 pVec->SetDefinition( aPiPlus );
194 else if( ran < 0.6819 )
195 pVec->SetDefinition( aPiZero );
196 else
197 pVec->SetDefinition( aPiMinus );
[962]198
199 // DHW: add following two lines to correct energy balance
200 // ++backwardCount;
201 // backwardMass += pVec->GetMass()/GeV;
[819]202 }
203 pVec->SetSide( -2 ); // backside particle
[962]204 // Not used extraMass += pVec->GetMass()/GeV;
[819]205 pVec->SetNewlyAdded( true );
206 vec.SetElement( vecLen++, pVec );
207 }
208 }
209
[962]210 // Masses of particles added from cascade not included in energy balance.
211 // That's correct for nucleons from the intra-nuclear cascade but not for
212 // pions from the cascade.
213
[819]214 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
215 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
216 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
217 G4bool secondaryDeleted;
218 G4double pMass;
219
220 while( eAvailable <= 0.0 ) // must eliminate a particle
221 {
222 secondaryDeleted = false;
223 for( i=(vecLen-1); i>=0; --i )
224 {
225 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
226 {
227 pMass = vec[i]->GetMass()/GeV;
228 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
229 --forwardCount;
230 forwardEnergy += pMass;
231 forwardMass -= pMass;
232 secondaryDeleted = true;
233 break;
234 }
235 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
236 {
237 pMass = vec[i]->GetMass()/GeV;
238 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
239 --backwardCount;
240 backwardEnergy += pMass;
241 backwardMass -= pMass;
242 secondaryDeleted = true;
243 break;
244 }
245 } // breaks go down to here
246
247 if( secondaryDeleted )
248 {
249 delete vec[vecLen-1];
250 --vecLen;
251 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
252 }
253 else
254 {
255 if( vecLen == 0 ) return false; // all secondaries have been eliminated
256 if( targetParticle.GetSide() == -1 )
257 {
258 pMass = targetParticle.GetMass()/GeV;
259 targetParticle = *vec[0];
260 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
261 --backwardCount;
262 backwardEnergy += pMass;
263 backwardMass -= pMass;
264 secondaryDeleted = true;
265 }
266 else if( targetParticle.GetSide() == 1 )
267 {
268 pMass = targetParticle.GetMass()/GeV;
269 targetParticle = *vec[0];
270 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
271 --forwardCount;
272 forwardEnergy += pMass;
273 forwardMass -= pMass;
274 secondaryDeleted = true;
275 }
276
277 if( secondaryDeleted )
278 {
279 delete vec[vecLen-1];
280 --vecLen;
281 }
282 else
283 {
284 if( currentParticle.GetSide() == -1 )
285 {
286 pMass = currentParticle.GetMass()/GeV;
287 currentParticle = *vec[0];
288 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
289 --backwardCount;
290 backwardEnergy += pMass;
291 backwardMass -= pMass;
292 secondaryDeleted = true;
293 }
294 else if( currentParticle.GetSide() == 1 )
295 {
296 pMass = currentParticle.GetMass()/GeV;
297 currentParticle = *vec[0];
298 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
299 --forwardCount;
300 forwardEnergy += pMass;
301 forwardMass -= pMass;
302 secondaryDeleted = true;
303 }
304 if( secondaryDeleted )
305 {
306 delete vec[vecLen-1];
307 --vecLen;
308 }
309 else break;
310
311 } // secondary not deleted
312 } // secondary not deleted
313
314 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
315 } // while
316
317 //
318 // This is the start of the TwoCluster function
319 // Choose multi-particle resonance masses by sampling
320 // P(M) = gc[g(M-M0)]**(c-1) *exp[-(g(M-M0))**c]
321 // for M > M0
322 //
323 // Use for the forward and backward clusters, but not
324 // the cascade cluster
325
326 const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
327 const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
328 G4int ntc = 0;
329
330 if (forwardCount < 1 || backwardCount < 1) return false; // array bounds protection
331
332 G4double rmc = forwardMass;
333 if (forwardCount > 1) {
334 ntc = std::min(3,forwardCount-2);
335 rmc += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
336 }
337 G4double rmd = backwardMass;
338 if( backwardCount > 1 ) {
339 ntc = std::min(3,backwardCount-2);
340 rmd += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
341 }
342
343 while( rmc+rmd > centerofmassEnergy )
344 {
345 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
346 {
347 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
348 rmc *= temp;
349 rmd *= temp;
350 }
351 else
352 {
353 rmc = 0.1*forwardMass + 0.9*rmc;
354 rmd = 0.1*backwardMass + 0.9*rmd;
355 }
356 }
357
358 G4ReactionProduct pseudoParticle[8];
359 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
360
361 pseudoParticle[1].SetMass( mOriginal*GeV );
362 pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
363 pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
364
365 pseudoParticle[2].SetMass( protonMass*MeV );
366 pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
367 pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
368 //
369 // transform into center of mass system
370 //
371 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
372 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
373 pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
374
375 // Calculate cm momentum for forward and backward masses
376 // W = sqrt(pf*pf + rmc*rmc) + sqrt(pf*pf + rmd*rmd)
377 // Solve for pf
378
379 const G4double pfMin = 0.0001;
380 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
381 pf *= pf;
382 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
383 pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
384 //
385 // set final state masses and energies in centre of mass system
386 //
387 pseudoParticle[3].SetMass( rmc*GeV );
388 pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
389
390 pseudoParticle[4].SetMass( rmd*GeV );
391 pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
392
393 //
394 // Get cm scattering angle by sampling t from tmin to tmax
395 //
396 const G4double bMin = 0.01;
397 const G4double b1 = 4.0;
398 const G4double b2 = 1.6;
399 G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
400 G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
401 G4double factor = 1.0 - std::exp(-dtb);
402 G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
403
404 costheta = std::max(-1.0, std::min(1.0, costheta));
405 G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
406 G4double phi = G4UniformRand() * twopi;
407 //
408 // calculate final state momenta in centre of mass system
409 //
410 pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
411 pf*sintheta*std::sin(phi)*GeV,
412 pf*costheta*GeV );
413 pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());
414
415 // Backward cluster of nucleons and pions from intra-nuclear cascade
416 // Set up in lab system and transform to cms
417
418 G4double pp, pp1;
419 if( nuclearExcitationCount > 0 )
420 {
421 const G4double ga = 1.2;
422 G4double ekit1 = 0.04;
423 G4double ekit2 = 0.6; // Max KE of cascade particle
424 if( ekOriginal <= 5.0 )
425 {
426 ekit1 *= ekOriginal*ekOriginal/25.0;
427 ekit2 *= ekOriginal*ekOriginal/25.0;
428 }
429 G4double scale = std::pow(ekit2/ekit1, 1.0-ga) - 1.0;
430 for( i=0; i<vecLen; ++i )
431 {
432 if( vec[i]->GetSide() == -2 )
433 {
434 G4double kineticE = ekit1*std::pow((1.0 + G4UniformRand()*scale), 1.0/(1.0-ga) );
435 vec[i]->SetKineticEnergy( kineticE*GeV );
436 G4double vMass = vec[i]->GetMass()/MeV;
437 G4double totalE = kineticE*GeV + vMass;
438 pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
439 G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
440 G4double sint = std::sqrt(1.0-cost*cost);
441 phi = twopi*G4UniformRand();
442 vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
443 pp*sint*std::sin(phi)*MeV,
444 pp*cost*MeV );
445 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
446 }
447 }
448 }
449
450 //
451 // Fragmentation of forward and backward clusters
452 //
453
454 currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
455 currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
456
457 targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
458 targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
459
460 pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
461 pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
462 pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
463
464 pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
465 pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
466 pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
467
468 G4double wgt;
469 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
470 if( forwardCount > 1 ) // tempV will contain the forward particles
471 {
472 G4FastVector<G4ReactionProduct,256> tempV;
473 tempV.Initialize( forwardCount );
474 G4bool constantCrossSection = true;
475 G4int tempLen = 0;
476 if( currentParticle.GetSide() == 1 )
477 tempV.SetElement( tempLen++, &currentParticle );
478 if( targetParticle.GetSide() == 1 )
479 tempV.SetElement( tempLen++, &targetParticle );
480 for( i=0; i<vecLen; ++i )
481 {
482 if( vec[i]->GetSide() == 1 )
483 {
484 if( tempLen < 18 )
485 tempV.SetElement( tempLen++, vec[i] );
486 else
487 {
488 vec[i]->SetSide( -1 );
489 continue;
490 }
491 }
492 }
493 if( tempLen >= 2 )
494 {
495 wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
496 constantCrossSection, tempV, tempLen );
497 if( currentParticle.GetSide() == 1 )
498 currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
499 if( targetParticle.GetSide() == 1 )
500 targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
501 for( i=0; i<vecLen; ++i )
502 {
503 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
504 }
505 }
506 }
507 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
508 if( backwardCount > 1 ) // tempV will contain the backward particles,
509 { // but not those created from the intranuclear cascade
510 G4FastVector<G4ReactionProduct,256> tempV;
511 tempV.Initialize( backwardCount );
512 G4bool constantCrossSection = true;
513 G4int tempLen = 0;
514 if( currentParticle.GetSide() == -1 )
515 tempV.SetElement( tempLen++, &currentParticle );
516 if( targetParticle.GetSide() == -1 )
517 tempV.SetElement( tempLen++, &targetParticle );
518 for( i=0; i<vecLen; ++i )
519 {
520 if( vec[i]->GetSide() == -1 )
521 {
522 if( tempLen < 18 )
523 tempV.SetElement( tempLen++, vec[i] );
524 else
525 {
526 vec[i]->SetSide( -2 );
527 vec[i]->SetKineticEnergy( 0.0 );
528 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
529 continue;
530 }
531 }
532 }
533 if( tempLen >= 2 )
534 {
535 wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
536 constantCrossSection, tempV, tempLen );
537 if( currentParticle.GetSide() == -1 )
538 currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
539 if( targetParticle.GetSide() == -1 )
540 targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
541 for( i=0; i<vecLen; ++i )
542 {
543 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
544 }
545 }
546 }
547
548 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
549 //
550 // Lorentz transformation in lab system
551 //
552 currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
553 targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
554 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
555
556 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
557 //
558 // sometimes the leading strange particle is lost, set it back
559 //
560 G4bool dum = true;
561 if( leadFlag )
562 {
563 // leadFlag will be true
564 // iff original particle is strange AND if incident particle is strange
565 // leadFlag is set to the incident particle
566 // or
567 // target particle is strange leadFlag is set to the target particle
568
569 if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
570 dum = false;
571 else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
572 dum = false;
573 else
574 {
575 for( i=0; i<vecLen; ++i )
576 {
577 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
578 {
579 dum = false;
580 break;
581 }
582 }
583 }
584 if( dum )
585 {
586 G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
587 G4double ekin;
588 if( ((leadMass < protonMass) && (targetParticle.GetMass()/MeV < protonMass)) ||
589 ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
590 {
591 ekin = targetParticle.GetKineticEnergy()/GeV;
592 pp1 = targetParticle.GetMomentum().mag()/MeV; // old momentum
593 targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
594 targetParticle.SetKineticEnergy( ekin*GeV );
595 pp = targetParticle.GetTotalMomentum()/MeV; // new momentum
596 if( pp1 < 1.0e-3 ) {
597 G4ThreeVector iso = Isotropic(pp);
598 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
599 } else {
600 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
601 }
602 targetHasChanged = true;
603 }
604 else
605 {
606 ekin = currentParticle.GetKineticEnergy()/GeV;
607 pp1 = currentParticle.GetMomentum().mag()/MeV;
608 currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
609 currentParticle.SetKineticEnergy( ekin*GeV );
610 pp = currentParticle.GetTotalMomentum()/MeV;
611 if( pp1 < 1.0e-3 ) {
612 G4ThreeVector iso = Isotropic(pp);
613 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
614 } else {
615 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
616 }
617 incidentHasChanged = true;
618 }
619 }
620 } // end of if( leadFlag )
621
622 // Get number of final state nucleons and nucleons remaining in
623 // target nucleus
624
625 std::pair<G4int, G4int> finalStateNucleons =
626 GetFinalStateNucleons(originalTarget, vec, vecLen);
627
628 G4int protonsInFinalState = finalStateNucleons.first;
629 G4int neutronsInFinalState = finalStateNucleons.second;
630
631 G4int numberofFinalStateNucleons =
632 protonsInFinalState + neutronsInFinalState;
633
634 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
635 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
636 originalIncident->GetDefinition()->GetPDGMass() <
637 G4Lambda::Lambda()->GetPDGMass())
638 numberofFinalStateNucleons++;
639
640 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
641
642 G4int PinNucleus = std::max(0,
643 G4int(targetNucleus.GetZ()) - protonsInFinalState);
644 G4int NinNucleus = std::max(0,
645 G4int(targetNucleus.GetN()-targetNucleus.GetZ()) - neutronsInFinalState);
646 //
647 // for various reasons, the energy balance is not sufficient,
648 // check that, energy balance, angle of final system, etc.
649 //
650 pseudoParticle[4].SetMass( mOriginal*GeV );
651 pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
652 pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
653
654 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
655 G4int diff = 0;
656 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
657 if(numberofFinalStateNucleons == 1) diff = 0;
658 pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
659 pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
660 pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
661
662 G4double theoreticalKinetic =
663 pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
664
665 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
666 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
667 pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
668
669 if( vecLen < 16 )
670 {
671 G4ReactionProduct tempR[130];
672 tempR[0] = currentParticle;
673 tempR[1] = targetParticle;
674 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
675
676 G4FastVector<G4ReactionProduct,256> tempV;
677 tempV.Initialize( vecLen+2 );
678 G4bool constantCrossSection = true;
679 G4int tempLen = 0;
680 for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
681
682 if( tempLen >= 2 )
683 {
684 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
685 wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
686 pseudoParticle[5].GetTotalEnergy()/MeV,
687 constantCrossSection, tempV, tempLen );
688 if (wgt == -1) {
689 G4double Qvalue = 0;
690 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
691 wgt = GenerateNBodyEvent( Qvalue/MeV,
692 constantCrossSection, tempV, tempLen );
693 }
694 theoreticalKinetic = 0.0;
695 for( i=0; i<vecLen+2; ++i )
696 {
697 pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
698 pseudoParticle[7].SetMass( tempV[i]->GetMass() );
699 pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
700 pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
701 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
702 }
703 }
704 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
705 }
706 else
707 {
708 theoreticalKinetic -=
709 ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
710 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
711 }
712 G4double simulatedKinetic =
713 currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
714 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
715
716 // make sure that kinetic energies are correct
717 // the backward nucleon cluster is not produced within proper kinematics!!!
718
719 if( simulatedKinetic != 0.0 )
720 {
721 wgt = (theoreticalKinetic)/simulatedKinetic;
722 currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
723 pp = currentParticle.GetTotalMomentum()/MeV;
724 pp1 = currentParticle.GetMomentum().mag()/MeV;
725 if( pp1 < 0.001*MeV ) {
726 G4ThreeVector iso = Isotropic(pp);
727 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
728 } else {
729 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
730 }
731
732 targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
733 pp = targetParticle.GetTotalMomentum()/MeV;
734 pp1 = targetParticle.GetMomentum().mag()/MeV;
735 if( pp1 < 0.001*MeV ) {
736 G4ThreeVector iso = Isotropic(pp);
737 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
738 } else {
739 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
740 }
741
742 for( i=0; i<vecLen; ++i )
743 {
744 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
745 pp = vec[i]->GetTotalMomentum()/MeV;
746 pp1 = vec[i]->GetMomentum().mag()/MeV;
747 if( pp1 < 0.001 ) {
748 G4ThreeVector iso = Isotropic(pp);
749 vec[i]->SetMomentum( iso.x(), iso.y(), iso.z() );
750 } else {
751 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
752 }
753 }
754 }
755 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
756
757 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
758 modifiedOriginal, originalIncident, targetNucleus,
759 currentParticle, targetParticle, vec, vecLen );
760
761 // Add black track particles
762 // the total number of particles produced is restricted to 198
763 // this may have influence on very high energies
764
765 if( atomicWeight >= 1.5 )
766 {
767 // npnb is number of proton/neutron black track particles
768 // ndta is the number of deuterons, tritons, and alphas produced
769 // epnb is the kinetic energy available for proton/neutron black track
770 // particles
771 // edta is the kinetic energy available for deuteron/triton/alpha
772 // particles
773
774 G4int npnb = 0;
775 G4int ndta = 0;
776
777 G4double epnb, edta;
778 if (veryForward) {
779 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
780 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
781 } else {
782 epnb = targetNucleus.GetPNBlackTrackEnergy();
783 edta = targetNucleus.GetDTABlackTrackEnergy();
784 }
785
786 const G4double pnCutOff = 0.001; // GeV
787 const G4double dtaCutOff = 0.001; // GeV
[962]788 // const G4double kineticMinimum = 1.e-6;
789 // const G4double kineticFactor = -0.005;
[819]790
[962]791 // G4double sprob = 0.0; // sprob = probability of self-absorption in
[819]792 // heavy molecules
[962]793 // Not currently used (DHW 9 June 2008) const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
794 // if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
[819]795
796 if( epnb >= pnCutOff )
797 {
798 npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
799 if( numberofFinalStateNucleons + npnb > atomicWeight )
800 npnb = G4int(atomicWeight - numberofFinalStateNucleons);
801 npnb = std::min( npnb, 127-vecLen );
802 }
803 if( edta >= dtaCutOff )
804 {
805 ndta = G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
806 ndta = std::min( ndta, 127-vecLen );
807 }
808 if (npnb == 0 && ndta == 0) npnb = 1;
809
810 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
811
[962]812 AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
[819]813 PinNucleus, NinNucleus, targetNucleus,
814 vec, vecLen );
815 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
816 }
817
818 //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
819 // MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
820 //
821 // calculate time delay for nuclear reactions
822 //
823 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
824 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
825 else
826 currentParticle.SetTOF( 1.0 );
827
828 return true;
829}
830
831 /* end of file */
Note: See TracBrowser for help on using the repository browser.