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

Last change on this file since 1228 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 41.1 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: G4RPGFragmentation.cc,v 1.6 2008/06/09 18:13:16 dennis Exp $
27// GEANT4 tag $Name: geant4-09-03 $
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 // Based on 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 G4int itemp = G4int( G4UniformRand()*vecLen );
118 G4ReactionProduct pTemp = *vec[itemp];
119 *vec[itemp] = *vec[i];
120 *vec[i] = pTemp;
121 }
122
123 if (currentMass == 0.0 && targetMass == 0.0) {
124 // Target and projectile have annihilated. Replace them with the first
125 // two secondaries in the list. Current particle KE is maintained.
126
127 G4double ek = currentParticle.GetKineticEnergy();
128 G4ThreeVector m = currentParticle.GetMomentum();
129 currentParticle = *vec[0];
130 currentParticle.SetSide(1);
131 targetParticle = *vec[1];
132 targetParticle.SetSide(-1);
133 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
134 G4ReactionProduct *temp = vec[vecLen-1];
135 delete temp;
136 temp = vec[vecLen-2];
137 delete temp;
138 vecLen -= 2;
139 currentMass = currentParticle.GetMass()/GeV;
140 targetMass = targetParticle.GetMass()/GeV;
141 incidentHasChanged = true;
142 targetHasChanged = true;
143 currentParticle.SetKineticEnergy( ek );
144 currentParticle.SetMomentum( m );
145 veryForward = true;
146 }
147 const G4double atomicWeight = targetNucleus.GetN();
148 const G4double atomicNumber = targetNucleus.GetZ();
149 const G4double protonMass = aProton->GetPDGMass();
150
151 if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
152 && G4UniformRand() >= 0.7) {
153 G4ReactionProduct temp = currentParticle;
154 currentParticle = targetParticle;
155 targetParticle = temp;
156 incidentHasChanged = true;
157 targetHasChanged = true;
158 currentMass = currentParticle.GetMass()/GeV;
159 targetMass = targetParticle.GetMass()/GeV;
160 }
161 const G4double afc = std::min( 0.75,
162 0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
163 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
164
165 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
166 G4double forwardEnergy = freeEnergy/2.;
167 G4int forwardCount = 1; // number of particles in forward hemisphere
168
169 G4double backwardEnergy = freeEnergy/2.;
170 G4int backwardCount = 1; // number of particles in backward hemisphere
171
172 if(veryForward) {
173 if(currentParticle.GetSide()==-1)
174 {
175 forwardEnergy += currentMass;
176 forwardCount --;
177 backwardEnergy -= currentMass;
178 backwardCount ++;
179 }
180 if(targetParticle.GetSide()!=-1)
181 {
182 backwardEnergy += targetMass;
183 backwardCount --;
184 forwardEnergy -= targetMass;
185 forwardCount ++;
186 }
187 }
188
189 for (i=0; i<vecLen; ++i) {
190 if( vec[i]->GetSide() == -1 )
191 {
192 ++backwardCount;
193 backwardEnergy -= vec[i]->GetMass()/GeV;
194 } else {
195 ++forwardCount;
196 forwardEnergy -= vec[i]->GetMass()/GeV;
197 }
198 }
199
200 // Check that sum of forward particle masses does not exceed forwardEnergy,
201 // and similarly for backward particles. If so, move particles from one
202 // hemisphere to another.
203
204 if (backwardEnergy < 0.0) {
205 for (i = 0; i < vecLen; ++i) {
206 if (vec[i]->GetSide() == -1) {
207 backwardEnergy += vec[i]->GetMass()/GeV;
208 --backwardCount;
209 vec[i]->SetSide(1);
210 forwardEnergy -= vec[i]->GetMass()/GeV;
211 ++forwardCount;
212 if (backwardEnergy > 0.0) break;
213 }
214 }
215 }
216
217 if (forwardEnergy < 0.0) {
218 for (i = 0; i < vecLen; ++i) {
219 if (vec[i]->GetSide() == 1) {
220 forwardEnergy += vec[i]->GetMass()/GeV;
221 --forwardCount;
222 vec[i]->SetSide(-1);
223 backwardEnergy -= vec[i]->GetMass()/GeV;
224 ++backwardCount;
225 if (forwardEnergy > 0.0) break;
226 }
227 }
228 }
229
230 // Special cases for reactions near threshold
231
232 // 1. There is only one secondary
233 if (forwardEnergy > 0.0 && backwardEnergy < 0.0) {
234 forwardEnergy += backwardEnergy;
235 backwardEnergy = 0;
236 }
237
238 // 2. Nuclear binding energy is large
239 if (forwardEnergy + backwardEnergy < 0.0) return false;
240
241
242 // forwardEnergy now represents the total energy in the forward reaction
243 // hemisphere which is available for kinetic energy and particle creation.
244 // Similarly for backwardEnergy.
245
246 // Add particles from the intranuclear cascade.
247 // nuclearExcitationCount = number of new secondaries produced by nuclear
248 // excitation
249 // extraCount = number of nucleons within these new secondaries
250 //
251 // Note: eventually have to make sure that enough nucleons are available
252 // in the case of small target nuclei
253
254 G4double xtarg;
255 if( centerofmassEnergy < (2.0+G4UniformRand()) )
256 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
257 else
258 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
259 if( xtarg <= 0.0 )xtarg = 0.01;
260 G4int nuclearExcitationCount = G4Poisson( xtarg );
261 // To do: try reducing nuclearExcitationCount with increasing energy
262 // to simulate cut-off of cascade
263 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
264 G4int extraNucleonCount = 0;
265 G4double extraNucleonMass = 0.0;
266
267 if (nuclearExcitationCount > 0) {
268 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
269 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
270 G4int momentumBin = 0;
271 while( (momentumBin < 6) &&
272 (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
273 ++momentumBin;
274 momentumBin = std::min( 5, momentumBin );
275
276 // NOTE: in GENXPT, these new particles were given negative codes
277 // here I use NewlyAdded = true instead
278 //
279
280 for (i = 0; i < nuclearExcitationCount; ++i) {
281 G4ReactionProduct * pVec = new G4ReactionProduct();
282 if (G4UniformRand() < nucsup[momentumBin]) {
283
284 if (G4UniformRand() > 1.0-atomicNumber/atomicWeight)
285 pVec->SetDefinition( aProton );
286 else
287 pVec->SetDefinition( aNeutron );
288
289 // nucleon comes from nucleus -
290 // do not subtract its mass from backward energy
291 pVec->SetSide( -2 ); // -2 means backside nucleon
292 ++extraNucleonCount;
293 extraNucleonMass += pVec->GetMass()/GeV;
294 // To do: Remove chosen nucleon from target nucleus
295 pVec->SetNewlyAdded( true );
296 vec.SetElement( vecLen++, pVec );
297 ++backwardCount;
298
299 } else {
300
301 G4double ran = G4UniformRand();
302 if( ran < 0.3181 )
303 pVec->SetDefinition( aPiPlus );
304 else if( ran < 0.6819 )
305 pVec->SetDefinition( aPiZero );
306 else
307 pVec->SetDefinition( aPiMinus );
308
309 if (backwardEnergy > pVec->GetMass()/GeV) {
310 backwardEnergy -= pVec->GetMass()/GeV; // pion mass comes from free energy
311 ++backwardCount;
312 pVec->SetSide( -1 ); // backside particle, but not a nucleon
313 pVec->SetNewlyAdded( true );
314 vec.SetElement( vecLen++, pVec );
315 }
316
317 // To do: Change proton to neutron (or vice versa) in target nucleus depending
318 // on pion charge
319 }
320 }
321 }
322
323 // Define initial state vectors for Lorentz transformations
324 // The pseudoParticles have non-standard masses, hence the "pseudo"
325
326 G4ReactionProduct pseudoParticle[8];
327 for (i = 0; i < 8; ++i) pseudoParticle[i].SetZero();
328
329 pseudoParticle[0].SetMass( mOriginal*GeV );
330 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
331 pseudoParticle[0].SetTotalEnergy(
332 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
333
334 pseudoParticle[1].SetMass(protonMass); // this could be targetMass
335 pseudoParticle[1].SetTotalEnergy(protonMass);
336
337 pseudoParticle[3].SetMass(protonMass*(1+extraNucleonCount) );
338 pseudoParticle[3].SetTotalEnergy(protonMass*(1+extraNucleonCount) );
339
340 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
341 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
342
343 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
344 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
345
346 // Main loop for 4-momentum generation
347 // See Pitha-report (Aachen) for a detailed description of the method
348
349 G4double aspar, pt, et, x, pp, pp1, wgt;
350 G4int innerCounter, outerCounter;
351 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
352
353 G4double forwardKinetic = 0.0;
354 G4double backwardKinetic = 0.0;
355
356 // Process the secondary particles in reverse order
357 // The incident particle is done after the secondaries
358 // Nucleons, including the target, in the backward hemisphere are also
359 // done later
360
361 G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere
362 G4double totalEnergy, kineticEnergy, vecMass;
363 G4double phi;
364
365 for (i = vecLen-1; i >= 0; --i) {
366
367 if (vec[i]->GetNewlyAdded()) { // added from intranuclear cascade
368 if (vec[i]->GetSide() == -2) { // its a nucleon
369 if (backwardNucleonCount < 18) {
370 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
371 for(G4int i=0; i<vecLen; i++) delete vec[i];
372 vecLen = 0;
373 throw G4HadReentrentException(__FILE__, __LINE__,
374 "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon");
375 }
376 vec[i]->SetSide(-3);
377 ++backwardNucleonCount;
378 continue; // Don't generate momenta for the first 17 backward
379 // cascade nucleons. This gets done by the cluster
380 // model later on.
381 }
382 }
383 }
384
385 // Set pt and phi values, they are changed somewhat in the iteration loop
386 // Set mass parameter for lambda fragmentation model
387
388 vecMass = vec[i]->GetMass()/GeV;
389 G4double ran = -std::log(1.0-G4UniformRand())/3.5;
390
391 if (vec[i]->GetSide() == -2) { // backward nucleon
392 aspar = 0.20;
393 pt = std::sqrt( std::pow( ran, 1.2 ) );
394
395 } else { // not a backward nucleon
396 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
397 aspar = 0.75;
398 pt = std::sqrt( std::pow( ran, 1.7 ) );
399 } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
400 aspar = 0.70;
401 pt = std::sqrt( std::pow( ran, 1.7 ) );
402 } else { // vec[i] must be a baryon or ion
403 aspar = 0.65;
404 pt = std::sqrt( std::pow( ran, 1.5 ) );
405 }
406 }
407
408 pt = std::max( 0.001, pt );
409 phi = G4UniformRand()*twopi;
410 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
411 if (vec[i]->GetSide() > 0)
412 et = pseudoParticle[0].GetTotalEnergy()/GeV;
413 else
414 et = pseudoParticle[1].GetTotalEnergy()/GeV;
415
416 //
417 // Start of outer iteration loop
418 //
419 outerCounter = 0;
420 eliminateThisParticle = true;
421 resetEnergies = true;
422 dndl[0] = 0.0;
423
424 while (++outerCounter < 3) {
425
426 FragmentationIntegral(pt, et, aspar, vecMass);
427 innerCounter = 0;
428 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
429
430 // Start of inner iteration loop
431
432 while (++innerCounter < 7) {
433
434 ran = G4UniformRand()*dndl[19];
435 l = 1;
436 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
437 x = (G4double(l-1) + G4UniformRand())/19.;
438 if (vec[i]->GetSide() < 0) x *= -1.;
439 vec[i]->SetMomentum( x*et*GeV ); // set the z-momentum
440 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
441 vec[i]->SetTotalEnergy( totalEnergy*GeV );
442 kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
443
444 if (vec[i]->GetSide() > 0) { // forward side
445 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy ) {
446 // Leave at least 5% of the forward free energy for the projectile
447
448 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
449 forwardKinetic += kineticEnergy;
450 outerCounter = 2; // leave outer loop
451 eliminateThisParticle = false; // don't eliminate this particle
452 resetEnergies = false;
453 break; // leave inner loop
454 }
455 if( innerCounter > 5 )break; // leave inner loop
456 if( backwardEnergy >= vecMass ) // switch sides
457 {
458 vec[i]->SetSide(-1);
459 forwardEnergy += vecMass;
460 backwardEnergy -= vecMass;
461 ++backwardCount;
462 }
463 } else { // backward side
464 // if (extraNucleonCount > 19) x = 0.999;
465 // G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
466 // DHW: I think above lines were meant to be as follows:
467 G4double xxx = 0.999;
468 if (extraNucleonCount < 20) xxx = 0.95+0.05*extraNucleonCount/20.0;
469
470 if ((backwardKinetic+kineticEnergy) < xxx*backwardEnergy) {
471 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
472 backwardKinetic += kineticEnergy;
473 outerCounter = 2; // leave outer loop
474 eliminateThisParticle = false; // don't eliminate this particle
475 resetEnergies = false;
476 break; // leave inner loop
477 }
478 if (innerCounter > 5) break; // leave inner loop
479 if (forwardEnergy >= vecMass) { // switch sides
480 vec[i]->SetSide(1);
481 forwardEnergy -= vecMass;
482 backwardEnergy += vecMass;
483 backwardCount--;
484 }
485 }
486 G4ThreeVector momentum = vec[i]->GetMomentum();
487 vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
488 pt *= 0.9;
489 dndl[19] *= 0.9;
490 } // closes inner loop
491
492 // If we get here, the inner loop has been done 6 times.
493 // If necessary, reduce energies of the previously done particles if
494 // they are lighter than protons or are in the forward hemisphere.
495 // Then continue with outer loop.
496
497 if (resetEnergies)
498 ReduceEnergiesOfSecondaries(i+1, forwardKinetic, backwardKinetic,
499 vec, vecLen,
500 pseudoParticle[4], pseudoParticle[5],
501 pt);
502
503 } // closes outer loop
504
505 if (eliminateThisParticle && vec[i]->GetMayBeKilled()) {
506 // not enough energy, eliminate this particle
507
508 if (vec[i]->GetSide() > 0) {
509 --forwardCount;
510 forwardEnergy += vecMass;
511 } else {
512 --backwardCount;
513 if (vec[i]->GetSide() == -2) {
514 --extraNucleonCount;
515 extraNucleonMass -= vecMass;
516 } else {
517 backwardEnergy += vecMass;
518 }
519 }
520
521 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
522 G4ReactionProduct* temp = vec[vecLen-1];
523 delete temp;
524 // To do: modify target nucleus according to particle eliminated
525
526 if( --vecLen == 0 ){
527 G4cout << " FALSE RETURN DUE TO ENERGY BALANCE " << G4endl;
528 return false;
529 } // all the secondaries have been eliminated
530 }
531 } // closes main loop over secondaries
532
533 // Re-balance forward and backward energy if possible and if necessary
534
535 G4double forwardKEDiff = forwardEnergy - forwardKinetic;
536 G4double backwardKEDiff = backwardEnergy - backwardKinetic;
537
538 if (forwardKEDiff < 0.0 || backwardKEDiff < 0.0) {
539 ReduceEnergiesOfSecondaries(0, forwardKinetic, backwardKinetic,
540 vec, vecLen,
541 pseudoParticle[4], pseudoParticle[5],
542 pt);
543
544 forwardKEDiff = forwardEnergy - forwardKinetic;
545 backwardKEDiff = backwardEnergy - backwardKinetic;
546 if (backwardKEDiff < 0.0) {
547 if (forwardKEDiff + backwardKEDiff > 0.0) {
548 backwardEnergy = backwardKinetic;
549 forwardEnergy += backwardKEDiff;
550 forwardKEDiff = forwardEnergy - forwardKinetic;
551 backwardKEDiff = 0.0;
552 } else {
553 G4cout << " False return due to insufficient backward energy " << G4endl;
554 return false;
555 }
556 }
557
558 if (forwardKEDiff < 0.0) {
559 if (forwardKEDiff + backwardKEDiff > 0.0) {
560 forwardEnergy = forwardKinetic;
561 backwardEnergy += forwardKEDiff;
562 backwardKEDiff = backwardEnergy - backwardKinetic;
563 forwardKEDiff = 0.0;
564 } else {
565 G4cout << " False return due to insufficient forward energy " << G4endl;
566 return false;
567 }
568 }
569 }
570
571 // Generate momentum for incident (current) particle, which was placed
572 // in the forward hemisphere.
573 // Set mass parameter for lambda fragmentation model.
574 // Set pt and phi values, which are changed somewhat in the iteration loop
575
576 G4double ran = -std::log(1.0-G4UniformRand());
577 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
578 aspar = 0.60;
579 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
580 } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
581 aspar = 0.50;
582 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
583 } else {
584 aspar = 0.40;
585 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
586 }
587
588 phi = G4UniformRand()*twopi;
589 currentParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
590 et = pseudoParticle[0].GetTotalEnergy()/GeV;
591 dndl[0] = 0.0;
592 vecMass = currentParticle.GetMass()/GeV;
593
594 FragmentationIntegral(pt, et, aspar, vecMass);
595
596 ran = G4UniformRand()*dndl[19];
597 l = 1;
598 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
599 x = (G4double(l-1) + G4UniformRand())/19.;
600 currentParticle.SetMomentum( x*et*GeV ); // set the z-momentum
601
602 if (forwardEnergy < forwardKinetic) {
603 totalEnergy = vecMass + 0.04*std::fabs(normal());
604 G4cout << " Not enough forward energy: forwardEnergy = "
605 << forwardEnergy << " forwardKinetic = "
606 << forwardKinetic << " total energy left = "
607 << backwardKEDiff + forwardKEDiff << G4endl;
608 } else {
609 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
610 forwardKinetic = forwardEnergy;
611 }
612 currentParticle.SetTotalEnergy( totalEnergy*GeV );
613 pp = std::sqrt(std::abs( totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
614 pp1 = currentParticle.GetMomentum().mag();
615
616 if (pp1 < 1.0e-6*GeV) {
617 G4ThreeVector iso = Isotropic(pp);
618 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
619 } else {
620 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
621 }
622 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
623
624 // Current particle now finished
625
626 // Begin target particle
627
628 if (backwardNucleonCount < 18) {
629 targetParticle.SetSide(-3);
630 ++backwardNucleonCount;
631
632 } else {
633 // Set pt and phi values, they are changed somewhat in the iteration loop
634 // Set mass parameter for lambda fragmentation model
635
636 vecMass = targetParticle.GetMass()/GeV;
637 ran = -std::log(1.0-G4UniformRand());
638 aspar = 0.40;
639 pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
640 phi = G4UniformRand()*twopi;
641 targetParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
642 et = pseudoParticle[1].GetTotalEnergy()/GeV;
643 outerCounter = 0;
644 innerCounter = 0;
645 G4bool marginalEnergy = true;
646 dndl[0] = 0.0;
647 G4double xxx = 0.999;
648 if( extraNucleonCount < 20 ) xxx = 0.95+0.05*extraNucleonCount/20.0;
649 G4ThreeVector momentum;
650
651 while (++outerCounter < 4) {
652 FragmentationIntegral(pt, et, aspar, vecMass);
653
654 for (innerCounter = 0; innerCounter < 6; innerCounter++) {
655 ran = G4UniformRand()*dndl[19];
656 l = 1;
657 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
658 x = -(G4double(l-1) + G4UniformRand())/19.;
659 targetParticle.SetMomentum( x*et*GeV ); // set the z-momentum
660 totalEnergy = std::sqrt(x*et*x*et + pt*pt + vecMass*vecMass);
661 targetParticle.SetTotalEnergy( totalEnergy*GeV );
662
663 if ((backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy) {
664 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
665 backwardKinetic += totalEnergy - vecMass;
666 outerCounter = 3; // leave outer loop
667 marginalEnergy = false;
668 break; // leave inner loop
669 }
670 momentum = targetParticle.GetMomentum();
671 targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
672 pt *= 0.9;
673 dndl[19] *= 0.9;
674 }
675 }
676
677 if (marginalEnergy) {
678 G4cout << " Extra backward kinetic energy = "
679 << 0.999*backwardEnergy - backwardKinetic << G4endl;
680 totalEnergy = vecMass + 0.999*backwardEnergy - backwardKinetic;
681 targetParticle.SetTotalEnergy(totalEnergy*GeV);
682 pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
683 targetParticle.SetMomentum(momentum.x()/0.9, momentum.y()/0.9);
684 pp1 = targetParticle.GetMomentum().mag();
685 targetParticle.SetMomentum(targetParticle.GetMomentum() * pp/pp1 );
686 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
687 backwardKinetic = 0.999*backwardEnergy;
688 }
689
690 }
691
692 if (backwardEnergy < backwardKinetic)
693 G4cout << " Backward Edif = " << backwardEnergy - backwardKinetic << G4endl;
694 if (forwardEnergy != forwardKinetic)
695 G4cout << " Forward Edif = " << forwardEnergy - forwardKinetic << G4endl;
696
697 // Target particle finished.
698 // Now produce backward nucleons with a cluster model
699 // ps[2] = CM frame of projectile + target
700 // ps[3] = sum of projectile + nucleon cluster in lab frame
701 // ps[6] = proj + cluster 4-vector boosted into CM frame of proj + targ
702 // with secondaries, current and target particles subtracted
703 // = total 4-momentum of backward nucleon cluster
704
705 pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
706 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
707 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
708
709 if (backwardNucleonCount == 1) {
710 // Target particle is the only backward nucleon. Give it the remainder
711 // of the backward kinetic energy.
712
713 G4double ekin =
714 std::min(backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV);
715
716 if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
717 vecMass = targetParticle.GetMass()/GeV;
718 totalEnergy = ekin + vecMass;
719 targetParticle.SetTotalEnergy( totalEnergy*GeV );
720 pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
721 pp1 = pseudoParticle[6].GetMomentum().mag();
722 if (pp1 < 1.0e-6*GeV) {
723 G4ThreeVector iso = Isotropic(pp);
724 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
725 } else {
726 targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1));
727 }
728 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
729
730 } else if (backwardNucleonCount > 1) {
731 // Share remaining energy with up to 17 backward nucleons
732
733 G4int tempCount = 5;
734 if (backwardNucleonCount < 5) tempCount = backwardNucleonCount;
735 tempCount -= 2;
736
737 G4double clusterMass = 0.;
738 if (targetParticle.GetSide() == -3)
739 clusterMass = targetParticle.GetMass()/GeV;
740 for (i = 0; i < vecLen; ++i)
741 if (vec[i]->GetSide() == -3) clusterMass += vec[i]->GetMass()/GeV;
742 clusterMass += backwardEnergy - backwardKinetic;
743
744 totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
745 pseudoParticle[6].SetMass(clusterMass*GeV);
746
747 pp = std::sqrt(std::abs(totalEnergy*totalEnergy -
748 clusterMass*clusterMass) )*GeV;
749 pp1 = pseudoParticle[6].GetMomentum().mag();
750 if (pp1 < 1.0e-6*GeV) {
751 G4ThreeVector iso = Isotropic(pp);
752 pseudoParticle[6].SetMomentum(iso.x(), iso.y(), iso.z());
753 } else {
754 pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp/pp1));
755 }
756
757 std::vector<G4ReactionProduct*> tempList; // ptrs to backward nucleons
758 if (targetParticle.GetSide() == -3) tempList.push_back(&targetParticle);
759 for (i = 0; i < vecLen; ++i)
760 if (vec[i]->GetSide() == -3) tempList.push_back(vec[i]);
761
762 constantCrossSection = true;
763
764 if (tempList.size() > 1) {
765 G4int n_entry = 0;
766 wgt = GenerateNBodyEventT(pseudoParticle[6].GetMass(),
767 constantCrossSection, tempList);
768
769 if (targetParticle.GetSide() == -3) {
770 targetParticle = *tempList[0];
771 targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
772 n_entry++;
773 }
774
775 for (i = 0; i < vecLen; ++i) {
776 if (vec[i]->GetSide() == -3) {
777 *vec[i] = *tempList[n_entry];
778 vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
779 n_entry++;
780 }
781 }
782 }
783 } else return false;
784
785 if (vecLen == 0) return false; // all the secondaries have been eliminated
786
787 // Lorentz transformation to lab frame
788
789 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
790 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
791 for (i = 0; i < vecLen; ++i) vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
792
793 // Set leading strange particle flag.
794 // leadFlag will be true if original particle and incident particle are
795 // both strange, in which case the incident particle becomes the leading
796 // particle.
797 // leadFlag will also be true if the target particle is strange, but the
798 // incident particle is not, in which case the target particle becomes the
799 // leading particle.
800
801 G4bool leadingStrangeParticleHasChanged = true;
802 if (leadFlag)
803 {
804 if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
805 leadingStrangeParticleHasChanged = false;
806 if (leadingStrangeParticleHasChanged &&
807 (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()) )
808 leadingStrangeParticleHasChanged = false;
809 if( leadingStrangeParticleHasChanged )
810 {
811 for( i=0; i<vecLen; i++ )
812 {
813 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
814 {
815 leadingStrangeParticleHasChanged = false;
816 break;
817 }
818 }
819 }
820 if( leadingStrangeParticleHasChanged )
821 {
822 G4bool leadTest =
823 (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
824 leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
825 G4bool targetTest =
826 (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
827 targetParticle.GetDefinition()->GetParticleSubType() == "pi");
828
829 // following modified by JLC 22-Oct-97
830
831 if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
832 {
833 targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
834 targetHasChanged = true;
835 }
836 else
837 {
838 currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
839 incidentHasChanged = false;
840 }
841 }
842 } // end of if( leadFlag )
843
844 // Get number of final state nucleons and nucleons remaining in
845 // target nucleus
846
847 std::pair<G4int, G4int> finalStateNucleons =
848 GetFinalStateNucleons(originalTarget, vec, vecLen);
849
850 G4int protonsInFinalState = finalStateNucleons.first;
851 G4int neutronsInFinalState = finalStateNucleons.second;
852
853 G4int numberofFinalStateNucleons =
854 protonsInFinalState + neutronsInFinalState;
855
856 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
857 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
858 originalIncident->GetDefinition()->GetPDGMass() <
859 G4Lambda::Lambda()->GetPDGMass())
860 numberofFinalStateNucleons++;
861
862 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
863
864 G4int PinNucleus = std::max(0,
865 G4int(targetNucleus.GetZ()) - protonsInFinalState);
866 G4int NinNucleus = std::max(0,
867 G4int(targetNucleus.GetN()-targetNucleus.GetZ()) - neutronsInFinalState);
868
869 pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
870 pseudoParticle[3].SetMass( mOriginal*GeV );
871 pseudoParticle[3].SetTotalEnergy(
872 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
873
874 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
875 G4int diff = 0;
876 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
877 if(numberofFinalStateNucleons == 1) diff = 0;
878 pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
879 pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff) );
880 pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff) );
881
882 G4double theoreticalKinetic =
883 pseudoParticle[3].GetTotalEnergy() + pseudoParticle[4].GetTotalEnergy() -
884 currentParticle.GetMass() - targetParticle.GetMass();
885 for (i = 0; i < vecLen; ++i) theoreticalKinetic -= vec[i]->GetMass();
886
887 G4double simulatedKinetic =
888 currentParticle.GetKineticEnergy() + targetParticle.GetKineticEnergy();
889 for (i = 0; i < vecLen; ++i)
890 simulatedKinetic += vec[i]->GetKineticEnergy();
891
892 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
893 pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
894 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
895
896 pseudoParticle[7].SetZero();
897 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
898 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
899 for (i = 0; i < vecLen; ++i)
900 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
901
902 /*
903 // This code does not appear to do anything to the energy or angle spectra
904 if( vecLen <= 16 && vecLen > 0 )
905 {
906 // must create a new set of ReactionProducts here because GenerateNBody will
907 // modify the momenta for the particles, and we don't want to do this
908 //
909 G4ReactionProduct tempR[130];
910 tempR[0] = currentParticle;
911 tempR[1] = targetParticle;
912 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
913 G4FastVector<G4ReactionProduct,256> tempV1;
914 tempV1.Initialize( vecLen+2 );
915 G4int tempLen1 = 0;
916 for( i=0; i<vecLen+2; ++i )tempV1.SetElement( tempLen1++, &tempR[i] );
917 constantCrossSection = true;
918
919 wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() +
920 pseudoParticle[4].GetTotalEnergy(),
921 constantCrossSection, tempV1, tempLen1);
922 if (wgt == -1) {
923 G4double Qvalue = 0;
924 for (i = 0; i < tempLen1; i++) Qvalue += tempV1[i]->GetMass();
925 wgt = GenerateNBodyEvent(Qvalue,
926 constantCrossSection, tempV1, tempLen1);
927 }
928 if(wgt>-.5)
929 {
930 theoreticalKinetic = 0.0;
931 for( i=0; i<tempLen1; ++i )
932 {
933 pseudoParticle[6].Lorentz( *tempV1[i], pseudoParticle[4] );
934 theoreticalKinetic += pseudoParticle[6].GetKineticEnergy();
935 }
936 }
937 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
938 }
939 */
940
941 //
942 // Make sure that the kinetic energies are correct
943 //
944
945 if (simulatedKinetic != 0.0) {
946 wgt = theoreticalKinetic/simulatedKinetic;
947 theoreticalKinetic = currentParticle.GetKineticEnergy() * wgt;
948 simulatedKinetic = theoreticalKinetic;
949 currentParticle.SetKineticEnergy(theoreticalKinetic);
950 pp = currentParticle.GetTotalMomentum();
951 pp1 = currentParticle.GetMomentum().mag();
952 if (pp1 < 1.0e-6*GeV) {
953 G4ThreeVector iso = Isotropic(pp);
954 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
955 } else {
956 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp/pp1));
957 }
958
959 theoreticalKinetic = targetParticle.GetKineticEnergy() * wgt;
960 targetParticle.SetKineticEnergy(theoreticalKinetic);
961 simulatedKinetic += theoreticalKinetic;
962 pp = targetParticle.GetTotalMomentum();
963 pp1 = targetParticle.GetMomentum().mag();
964
965 if (pp1 < 1.0e-6*GeV) {
966 G4ThreeVector iso = Isotropic(pp);
967 targetParticle.SetMomentum(iso.x(), iso.y(), iso.z() );
968 } else {
969 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp/pp1) );
970 }
971
972 for (i = 0; i < vecLen; ++i ) {
973 theoreticalKinetic = vec[i]->GetKineticEnergy() * wgt;
974 simulatedKinetic += theoreticalKinetic;
975 vec[i]->SetKineticEnergy(theoreticalKinetic);
976 pp = vec[i]->GetTotalMomentum();
977 pp1 = vec[i]->GetMomentum().mag();
978 if( pp1 < 1.0e-6*GeV ) {
979 G4ThreeVector iso = Isotropic(pp);
980 vec[i]->SetMomentum(iso.x(), iso.y(), iso.z() );
981 } else {
982 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp/pp1) );
983 }
984 }
985 }
986
987 // Rotate(numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
988 // modifiedOriginal, originalIncident, targetNucleus,
989 // currentParticle, targetParticle, vec, vecLen );
990
991 // Add black track particles
992 // the total number of particles produced is restricted to 198
993 // this may have influence on very high energies
994
995 if( atomicWeight >= 1.5 )
996 {
997 // npnb is number of proton/neutron black track particles
998 // ndta is the number of deuterons, tritons, and alphas produced
999 // epnb is the kinetic energy available for proton/neutron black track
1000 // particles
1001 // edta is the kinetic energy available for deuteron/triton/alpha particles
1002
1003 G4int npnb = 0;
1004 G4int ndta = 0;
1005
1006 G4double epnb, edta;
1007 if (veryForward) {
1008 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
1009 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
1010 } else {
1011 epnb = targetNucleus.GetPNBlackTrackEnergy();
1012 edta = targetNucleus.GetDTABlackTrackEnergy();
1013 }
1014 /*
1015 G4ReactionProduct* fudge = new G4ReactionProduct();
1016 fudge->SetDefinition( aProton );
1017 G4double TT = epnb + edta;
1018 G4double MM = fudge->GetMass()/GeV;
1019 fudge->SetTotalEnergy(MM*GeV + TT*GeV);
1020 G4double pzz = std::sqrt(TT*(TT + 2.*MM));
1021 fudge->SetMomentum(0.0, 0.0, pzz*GeV);
1022 vec.SetElement(vecLen++, fudge);
1023 // G4cout << " Fudge = " << vec[vecLen-1]->GetKineticEnergy()/GeV
1024 << G4endl;
1025 */
1026
1027 const G4double pnCutOff = 0.001;
1028 const G4double dtaCutOff = 0.001;
1029 // const G4double kineticMinimum = 1.e-6;
1030 // const G4double kineticFactor = -0.010;
1031 // G4double sprob = 0.0; // sprob = probability of self-absorption in
1032 // heavy molecules
1033 // Not currently used (DHW 9 June 2008) const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
1034 // if (ekIncident >= 5.0) sprob = std::min(1.0, 0.6*std::log(ekIncident-4.0));
1035 if (epnb > pnCutOff)
1036 {
1037 npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
1038 if (numberofFinalStateNucleons + npnb > atomicWeight)
1039 npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
1040 npnb = std::min( npnb, 127-vecLen );
1041 }
1042 if( edta >= dtaCutOff )
1043 {
1044 ndta = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta));
1045 ndta = std::min( ndta, 127-vecLen );
1046 }
1047 if (npnb == 0 && ndta == 0) npnb = 1;
1048
1049 AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
1050 PinNucleus, NinNucleus, targetNucleus,
1051 vec, vecLen);
1052 }
1053
1054 // if( centerofmassEnergy <= (4.0+G4UniformRand()) )
1055 // MomentumCheck( modifiedOriginal, currentParticle, targetParticle,
1056 // vec, vecLen );
1057 //
1058 // calculate time delay for nuclear reactions
1059 //
1060
1061 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
1062 currentParticle.SetTOF(
1063 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
1064 else
1065 currentParticle.SetTOF( 1.0 );
1066 return true;
1067
1068}
1069
1070
1071void G4RPGFragmentation::
1072ReduceEnergiesOfSecondaries(G4int startingIndex,
1073 G4double& forwardKinetic,
1074 G4double& backwardKinetic,
1075 G4FastVector<G4ReactionProduct,256>& vec,
1076 G4int& vecLen,
1077 G4ReactionProduct& forwardPseudoParticle,
1078 G4ReactionProduct& backwardPseudoParticle,
1079 G4double& pt)
1080{
1081 // Reduce energies and pt of secondaries in order to maintain
1082 // energy conservation
1083
1084 G4double totalEnergy;
1085 G4double pp;
1086 G4double pp1;
1087 G4double px;
1088 G4double py;
1089 G4double mass;
1090 G4ReactionProduct* pVec;
1091 G4int i;
1092
1093 forwardKinetic = 0.0;
1094 backwardKinetic = 0.0;
1095 forwardPseudoParticle.SetZero();
1096 backwardPseudoParticle.SetZero();
1097
1098 for (i = startingIndex; i < vecLen; i++) {
1099 pVec = vec[i];
1100 if (pVec->GetSide() != -3) {
1101 mass = pVec->GetMass();
1102 totalEnergy = 0.95*pVec->GetTotalEnergy() + 0.05*mass;
1103 pVec->SetTotalEnergy(totalEnergy);
1104 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - mass*mass ) );
1105 pp1 = pVec->GetMomentum().mag();
1106 if (pp1 < 1.0e-6*GeV) {
1107 G4ThreeVector iso = Isotropic(pp);
1108 pVec->SetMomentum( iso.x(), iso.y(), iso.z() );
1109 } else {
1110 pVec->SetMomentum(pVec->GetMomentum() * (pp/pp1) );
1111 }
1112
1113 px = pVec->GetMomentum().x();
1114 py = pVec->GetMomentum().y();
1115 pt = std::max(1.0, std::sqrt( px*px + py*py ) )/GeV;
1116 if (pVec->GetSide() > 0) {
1117 forwardKinetic += pVec->GetKineticEnergy()/GeV;
1118 forwardPseudoParticle = forwardPseudoParticle + (*pVec);
1119 } else {
1120 backwardKinetic += pVec->GetKineticEnergy()/GeV;
1121 backwardPseudoParticle = backwardPseudoParticle + (*pVec);
1122 }
1123 }
1124 }
1125}
1126
1127 /* end of file */
Note: See TracBrowser for help on using the repository browser.