G4RPGFragmentation.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028  
00029 #include <iostream>
00030 #include <signal.h>
00031 
00032 #include "G4RPGFragmentation.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4AntiProton.hh"
00036 #include "G4AntiNeutron.hh"
00037 #include "Randomize.hh"
00038 #include "G4Poisson.hh"
00039 #include "G4HadReentrentException.hh"
00040 
00041 G4RPGFragmentation::G4RPGFragmentation()
00042  : G4RPGReaction()
00043 {
00044   for (G4int i = 0; i < 20; i++) dndl[i] = 0.0;
00045 }
00046 
00047 
00048 void G4RPGFragmentation::
00049 FragmentationIntegral(G4double pt, G4double et, G4double parMass, G4double secMass)
00050 {
00051   pt = std::max( 0.001, pt );
00052   G4double dx = 1./(19.*pt);
00053   G4double x;
00054   G4double term1;
00055   G4double term2;
00056 
00057   for (G4int i = 1; i < 20; i++) {
00058     x = (G4double(i) - 0.5)*dx;
00059     term1 = 1. + parMass*parMass*x*x;
00060     term2 = pt*x*et*pt*x*et + pt*pt + secMass*secMass;
00061     dndl[i] = dx / std::sqrt( term1*term1*term1*term2 )
00062             + dndl[i-1];
00063   }
00064 }
00065 
00066 
00067 G4bool G4RPGFragmentation::
00068 ReactionStage(const G4HadProjectile* originalIncident,
00069               G4ReactionProduct& modifiedOriginal,
00070               G4bool& incidentHasChanged,
00071               const G4DynamicParticle* originalTarget,
00072               G4ReactionProduct& targetParticle,
00073               G4bool& targetHasChanged,
00074               const G4Nucleus& targetNucleus,
00075               G4ReactionProduct& currentParticle,
00076               G4FastVector<G4ReactionProduct,256>& vec,
00077               G4int& vecLen,
00078               G4bool leadFlag,
00079               G4ReactionProduct& leadingStrangeParticle) 
00080 {
00081   // 
00082   // Based on H. Fesefeldt's original FORTRAN code GENXPT
00083   //
00084   // Generation of x- and pT- values for incident, target, and all secondary 
00085   // particles using a simple single variable description E D3S/DP3= F(Q) 
00086   // with Q^2 = (M*X)^2 + PT^2.  Final state kinematics are produced by an 
00087   // FF-type iterative cascade method
00088   //
00089   // Internal units are GeV
00090   //
00091     
00092   // Protection in case no secondary has been created. In that case use  
00093   // two-body scattering
00094   //
00095   if (vecLen == 0) return false;
00096 
00097   G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus();
00098   G4ParticleDefinition* aProton = G4Proton::Proton();
00099   G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
00100   G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus();
00101   G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
00102 
00103   G4int i, l;
00104   G4bool veryForward = false;
00105     
00106   const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
00107   const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
00108   const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
00109   const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
00110   G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
00111   G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00112                                       targetMass*targetMass +
00113                                       2.0*targetMass*etOriginal );  // GeV
00114   G4double currentMass = currentParticle.GetMass()/GeV;
00115   targetMass = targetParticle.GetMass()/GeV;
00116 
00117   // Randomize the order of the secondary particles.
00118   // Note that the current and target particles are not affected.
00119 
00120   for (i=0; i<vecLen; ++i) {
00121     G4int itemp = G4int( G4UniformRand()*vecLen );
00122     G4ReactionProduct pTemp = *vec[itemp];
00123     *vec[itemp] = *vec[i];
00124     *vec[i] = pTemp;
00125   }
00126 
00127   if (currentMass == 0.0 && targetMass == 0.0) {
00128     // Target and projectile have annihilated.  Replace them with the first 
00129     // two secondaries in the list.  Current particle KE is maintained.
00130  
00131     G4double ek = currentParticle.GetKineticEnergy();
00132     G4ThreeVector mom = currentParticle.GetMomentum();
00133     currentParticle = *vec[0];
00134     currentParticle.SetSide(1);
00135     targetParticle = *vec[1];
00136     targetParticle.SetSide(-1);
00137     for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
00138     G4ReactionProduct *temp = vec[vecLen-1];
00139     delete temp;
00140     temp = vec[vecLen-2];
00141     delete temp;
00142     vecLen -= 2;
00143     currentMass = currentParticle.GetMass()/GeV;
00144     targetMass = targetParticle.GetMass()/GeV;
00145     incidentHasChanged = true;
00146     targetHasChanged = true;
00147     currentParticle.SetKineticEnergy( ek );
00148     currentParticle.SetMomentum(mom);
00149     veryForward = true;
00150   }
00151   const G4double atomicWeight = targetNucleus.GetA_asInt();
00152   const G4double atomicNumber = targetNucleus.GetZ_asInt();
00153   const G4double protonMass = aProton->GetPDGMass();
00154 
00155   if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
00156       && G4UniformRand() >= 0.7) {
00157     G4ReactionProduct temp = currentParticle;
00158     currentParticle = targetParticle;
00159     targetParticle = temp;
00160     incidentHasChanged = true;
00161     targetHasChanged = true;
00162     currentMass = currentParticle.GetMass()/GeV;
00163     targetMass = targetParticle.GetMass()/GeV;
00164   }
00165   const G4double afc = std::min( 0.75,
00166         0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
00167         std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
00168     
00169   G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
00170   G4double forwardEnergy = freeEnergy/2.;
00171   G4int forwardCount = 1;         // number of particles in forward hemisphere
00172     
00173   G4double backwardEnergy = freeEnergy/2.;
00174   G4int backwardCount = 1;        // number of particles in backward hemisphere
00175 
00176   if(veryForward) {
00177     if(currentParticle.GetSide()==-1)
00178     {
00179       forwardEnergy += currentMass;
00180       forwardCount --;
00181       backwardEnergy -= currentMass;
00182       backwardCount ++;
00183     }
00184     if(targetParticle.GetSide()!=-1)
00185     {
00186       backwardEnergy += targetMass;
00187       backwardCount --;
00188       forwardEnergy -= targetMass;
00189       forwardCount ++;
00190     }
00191   }
00192 
00193   for (i=0; i<vecLen; ++i) {
00194     if( vec[i]->GetSide() == -1 )
00195     {
00196       ++backwardCount;
00197       backwardEnergy -= vec[i]->GetMass()/GeV;
00198     } else {
00199       ++forwardCount;
00200       forwardEnergy -= vec[i]->GetMass()/GeV;
00201     }
00202   }
00203 
00204   // Check that sum of forward particle masses does not exceed forwardEnergy, 
00205   // and similarly for backward particles.  If so, move particles from one
00206   // hemisphere to another.
00207 
00208   if (backwardEnergy < 0.0) {
00209     for (i = 0; i < vecLen; ++i) {
00210       if (vec[i]->GetSide() == -1) { 
00211         backwardEnergy += vec[i]->GetMass()/GeV;
00212         --backwardCount;
00213         vec[i]->SetSide(1);
00214         forwardEnergy -= vec[i]->GetMass()/GeV;
00215         ++forwardCount;
00216         if (backwardEnergy > 0.0) break;
00217       }
00218     }
00219   }
00220 
00221   if (forwardEnergy < 0.0) {
00222     for (i = 0; i < vecLen; ++i) {
00223       if (vec[i]->GetSide() == 1) { 
00224         forwardEnergy += vec[i]->GetMass()/GeV;
00225         --forwardCount;
00226         vec[i]->SetSide(-1);
00227         backwardEnergy -= vec[i]->GetMass()/GeV;
00228         ++backwardCount;
00229         if (forwardEnergy > 0.0) break;
00230       }
00231     }
00232   }
00233  
00234   // Special cases for reactions near threshold
00235 
00236   // 1. There is only one secondary 
00237   if (forwardEnergy > 0.0 && backwardEnergy < 0.0) {
00238     forwardEnergy += backwardEnergy;
00239     backwardEnergy = 0;
00240   }
00241 
00242   // 2. Nuclear binding energy is large 
00243   if (forwardEnergy + backwardEnergy < 0.0) return false;
00244 
00245 
00246   // forwardEnergy now represents the total energy in the forward reaction 
00247   // hemisphere which is available for kinetic energy and particle creation.
00248   // Similarly for backwardEnergy.
00249 
00250   // Add particles from the intranuclear cascade.
00251   // nuclearExcitationCount = number of new secondaries produced by nuclear 
00252   // excitation
00253   // extraCount = number of nucleons within these new secondaries
00254   //
00255   // Note: eventually have to make sure that enough nucleons are available 
00256   // in the case of small target nuclei
00257 
00258   G4double xtarg;
00259   if( centerofmassEnergy < (2.0+G4UniformRand()) )
00260     xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
00261   else
00262     xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
00263   if( xtarg <= 0.0 )xtarg = 0.01;
00264   G4int nuclearExcitationCount = G4Poisson( xtarg );
00265   // To do: try reducing nuclearExcitationCount with increasing energy 
00266   //        to simulate cut-off of cascade
00267   if(atomicWeight<1.0001) nuclearExcitationCount = 0;
00268   G4int extraNucleonCount = 0;
00269   G4double extraNucleonMass = 0.0;
00270 
00271   if (nuclearExcitationCount > 0) {
00272     const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
00273     const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
00274     G4int momentumBin = 0;
00275     while( (momentumBin < 6) &&
00276            (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
00277       ++momentumBin;
00278     momentumBin = std::min( 5, momentumBin );
00279 
00280     //  NOTE: in GENXPT, these new particles were given negative codes
00281     //        here I use  NewlyAdded = true  instead
00282     //
00283 
00284     for (i = 0; i < nuclearExcitationCount; ++i) {
00285       G4ReactionProduct * pVec = new G4ReactionProduct();
00286       if (G4UniformRand() < nucsup[momentumBin]) {
00287 
00288         if (G4UniformRand() > 1.0-atomicNumber/atomicWeight)
00289           pVec->SetDefinition( aProton );
00290         else
00291           pVec->SetDefinition( aNeutron );
00292 
00293         // nucleon comes from nucleus - 
00294         // do not subtract its mass from backward energy
00295         pVec->SetSide( -2 );                // -2 means backside nucleon
00296         ++extraNucleonCount;
00297         extraNucleonMass += pVec->GetMass()/GeV;
00298         // To do: Remove chosen nucleon from target nucleus
00299         pVec->SetNewlyAdded( true );
00300         vec.SetElement( vecLen++, pVec );    
00301         ++backwardCount;
00302 
00303       } else {
00304 
00305         G4double ran = G4UniformRand();
00306         if( ran < 0.3181 )
00307           pVec->SetDefinition( aPiPlus );
00308         else if( ran < 0.6819 )
00309           pVec->SetDefinition( aPiZero );
00310         else
00311           pVec->SetDefinition( aPiMinus );
00312 
00313         if (backwardEnergy > pVec->GetMass()/GeV) {
00314           backwardEnergy -= pVec->GetMass()/GeV;    // pion mass comes from free energy
00315           ++backwardCount;
00316           pVec->SetSide( -1 );         // backside particle, but not a nucleon
00317           pVec->SetNewlyAdded( true );
00318           vec.SetElement( vecLen++, pVec );
00319         }
00320 
00321         // To do: Change proton to neutron (or vice versa) in target nucleus depending 
00322         //        on pion charge
00323       }
00324     }
00325   }
00326 
00327   // Define initial state vectors for Lorentz transformations
00328   // The pseudoParticles have non-standard masses, hence the "pseudo"
00329 
00330   G4ReactionProduct pseudoParticle[8];
00331   for (i = 0; i < 8; ++i) pseudoParticle[i].SetZero();
00332     
00333   pseudoParticle[0].SetMass( mOriginal*GeV );
00334   pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00335   pseudoParticle[0].SetTotalEnergy(
00336    std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
00337 
00338   pseudoParticle[1].SetMass(protonMass);     // this could be targetMass
00339   pseudoParticle[1].SetTotalEnergy(protonMass);
00340     
00341   pseudoParticle[3].SetMass(protonMass*(1+extraNucleonCount) );
00342   pseudoParticle[3].SetTotalEnergy(protonMass*(1+extraNucleonCount) );
00343     
00344   pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
00345   pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
00346     
00347   pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
00348   pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
00349     
00350   // Main loop for 4-momentum generation
00351   // See Pitha-report (Aachen) for a detailed description of the method
00352 
00353   G4double aspar, pt, et, x, pp, pp1, wgt;
00354   G4int    innerCounter, outerCounter;
00355   G4bool   eliminateThisParticle, resetEnergies, constantCrossSection;
00356     
00357   G4double forwardKinetic = 0.0;
00358   G4double backwardKinetic = 0.0;
00359 
00360   // Process the secondary particles in reverse order
00361   // The incident particle is done after the secondaries
00362   // Nucleons, including the target, in the backward hemisphere are also 
00363   // done later
00364 
00365   G4int backwardNucleonCount = 0;  // number of nucleons in backward hemisphere
00366   G4double totalEnergy, kineticEnergy, vecMass;
00367   G4double phi;
00368 
00369   for (i = vecLen-1; i >= 0; --i) {
00370 
00371     if (vec[i]->GetNewlyAdded()) {         // added from intranuclear cascade
00372       if (vec[i]->GetSide() == -2) {       // its a nucleon
00373         if (backwardNucleonCount < 18) {
00374           if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00375             for (G4int j = 0; j < vecLen; j++) delete vec[j];
00376             vecLen = 0;
00377             throw G4HadReentrentException(__FILE__, __LINE__,
00378             "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon");
00379           }
00380           vec[i]->SetSide(-3);
00381           ++backwardNucleonCount;
00382           continue;     // Don't generate momenta for the first 17 backward 
00383                         // cascade nucleons. This gets done by the cluster 
00384                         // model later on.
00385         }
00386       }
00387     }
00388 
00389     // Set pt and phi values, they are changed somewhat in the iteration loop
00390     // Set mass parameter for lambda fragmentation model
00391 
00392     vecMass = vec[i]->GetMass()/GeV;
00393     G4double ran = -std::log(1.0-G4UniformRand())/3.5;
00394 
00395     if (vec[i]->GetSide() == -2) {  // backward nucleon
00396       aspar = 0.20;
00397       pt = std::sqrt( std::pow( ran, 1.2 ) );
00398 
00399     } else {                        // not a backward nucleon
00400       if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00401         aspar = 0.75;
00402         pt = std::sqrt( std::pow( ran, 1.7 ) );
00403       } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
00404         aspar = 0.70;
00405         pt = std::sqrt( std::pow( ran, 1.7 ) );
00406       } else {                        // vec[i] must be a baryon or ion 
00407         aspar = 0.65;
00408         pt = std::sqrt( std::pow( ran, 1.5 ) );
00409       }
00410     }
00411 
00412     pt = std::max( 0.001, pt );
00413     phi = G4UniformRand()*twopi;
00414     vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00415     if (vec[i]->GetSide() > 0)
00416       et = pseudoParticle[0].GetTotalEnergy()/GeV;
00417     else
00418       et = pseudoParticle[1].GetTotalEnergy()/GeV;
00419 
00420     //
00421     // Start of outer iteration loop
00422     //
00423     outerCounter = 0;
00424     eliminateThisParticle = true;
00425     resetEnergies = true;
00426     dndl[0] = 0.0;
00427 
00428     while (++outerCounter < 3) {
00429 
00430       FragmentationIntegral(pt, et, aspar, vecMass);
00431       innerCounter = 0;
00432       vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00433 
00434       // Start of inner iteration loop
00435 
00436       while (++innerCounter < 7) {
00437 
00438         ran = G4UniformRand()*dndl[19];
00439         l = 1;
00440         while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
00441         x = (G4double(l-1) + G4UniformRand())/19.;
00442         if (vec[i]->GetSide() < 0) x *= -1.;
00443         vec[i]->SetMomentum( x*et*GeV );              // set the z-momentum
00444         totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
00445         vec[i]->SetTotalEnergy( totalEnergy*GeV );
00446         kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
00447 
00448         if (vec[i]->GetSide() > 0) {                  // forward side
00449           if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy ) {
00450             // Leave at least 5% of the forward free energy for the projectile
00451 
00452             pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
00453             forwardKinetic += kineticEnergy;
00454             outerCounter = 2;                  // leave outer loop
00455             eliminateThisParticle = false;     // don't eliminate this particle
00456             resetEnergies = false;
00457             break;                             // leave inner loop
00458           }
00459           if( innerCounter > 5 )break;         // leave inner loop
00460           if( backwardEnergy >= vecMass )      // switch sides
00461           {
00462             vec[i]->SetSide(-1);
00463             forwardEnergy += vecMass;
00464             backwardEnergy -= vecMass;
00465             ++backwardCount;
00466           }
00467         } else {                               // backward side
00468           //  if (extraNucleonCount > 19) x = 0.999;
00469           //  G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
00470           //  DHW: I think above lines were meant to be as follows:
00471           G4double xxx = 0.999;
00472           if (extraNucleonCount < 20) xxx = 0.95+0.05*extraNucleonCount/20.0;
00473 
00474           if ((backwardKinetic+kineticEnergy) < xxx*backwardEnergy) {
00475             pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
00476             backwardKinetic += kineticEnergy;
00477             outerCounter = 2;                  // leave outer loop
00478             eliminateThisParticle = false;     // don't eliminate this particle
00479             resetEnergies = false;
00480             break;                             // leave inner loop
00481           }
00482           if (innerCounter > 5) break;         // leave inner loop
00483           if (forwardEnergy >= vecMass) {      // switch sides
00484             vec[i]->SetSide(1);
00485             forwardEnergy -= vecMass;
00486             backwardEnergy += vecMass;
00487             backwardCount--;
00488           }
00489         }
00490         G4ThreeVector momentum = vec[i]->GetMomentum();
00491         vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
00492         pt *= 0.9;
00493         dndl[19] *= 0.9;
00494       }                      // closes inner loop
00495 
00496       // If we get here, the inner loop has been done 6 times.
00497       // If necessary, reduce energies of the previously done particles if 
00498       // they are lighter than protons or are in the forward hemisphere.
00499       // Then continue with outer loop.
00500 
00501       if (resetEnergies)
00502         ReduceEnergiesOfSecondaries(i+1, forwardKinetic, backwardKinetic,
00503                                     vec, vecLen,
00504                                     pseudoParticle[4], pseudoParticle[5],
00505                                     pt);
00506 
00507     } // closes outer loop
00508               
00509     if (eliminateThisParticle && vec[i]->GetMayBeKilled()) {
00510       // not enough energy, eliminate this particle
00511 
00512       if (vec[i]->GetSide() > 0) {
00513         --forwardCount;
00514         forwardEnergy += vecMass;
00515       } else {
00516         --backwardCount;
00517         if (vec[i]->GetSide() == -2) {
00518           --extraNucleonCount;
00519           extraNucleonMass -= vecMass;
00520         } else {
00521           backwardEnergy += vecMass;
00522         }
00523       }
00524 
00525       for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];    // shift up
00526       G4ReactionProduct* temp = vec[vecLen-1];
00527       delete temp;
00528       // To do: modify target nucleus according to particle eliminated  
00529 
00530       if( --vecLen == 0 ){
00531         G4cout << " FALSE RETURN DUE TO ENERGY BALANCE " << G4endl;
00532         return false;
00533       }  // all the secondaries have been eliminated
00534     }
00535   } // closes main loop over secondaries
00536 
00537   // Re-balance forward and backward energy if possible and if necessary
00538 
00539   G4double forwardKEDiff = forwardEnergy - forwardKinetic;
00540   G4double backwardKEDiff = backwardEnergy - backwardKinetic;
00541 
00542   if (forwardKEDiff < 0.0 || backwardKEDiff < 0.0) {
00543     ReduceEnergiesOfSecondaries(0, forwardKinetic, backwardKinetic,
00544                                 vec, vecLen,
00545                                 pseudoParticle[4], pseudoParticle[5],
00546                                 pt);
00547   
00548     forwardKEDiff = forwardEnergy - forwardKinetic;
00549     backwardKEDiff = backwardEnergy - backwardKinetic;
00550     if (backwardKEDiff < 0.0) {
00551       if (forwardKEDiff + backwardKEDiff > 0.0) {
00552         backwardEnergy = backwardKinetic;
00553         forwardEnergy += backwardKEDiff;
00554         forwardKEDiff = forwardEnergy - forwardKinetic;
00555         backwardKEDiff = 0.0;
00556       } else {
00557         G4cout << " False return due to insufficient backward energy " << G4endl; 
00558         return false;
00559       }
00560     }
00561       
00562     if (forwardKEDiff < 0.0) {
00563       if (forwardKEDiff + backwardKEDiff > 0.0) {
00564         forwardEnergy = forwardKinetic;
00565         backwardEnergy += forwardKEDiff;
00566         backwardKEDiff = backwardEnergy - backwardKinetic;
00567         forwardKEDiff = 0.0;
00568       } else {
00569         G4cout << " False return due to insufficient forward energy " << G4endl; 
00570         return false;
00571       }
00572     }
00573   }
00574 
00575   // Generate momentum for incident (current) particle, which was placed 
00576   // in the forward hemisphere. 
00577   // Set mass parameter for lambda fragmentation model.
00578   // Set pt and phi values, which are changed somewhat in the iteration loop
00579 
00580   G4double ran = -std::log(1.0-G4UniformRand());
00581   if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
00582     aspar = 0.60;
00583     pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
00584   } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
00585     aspar = 0.50;
00586     pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
00587   } else {
00588     aspar = 0.40;
00589     pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
00590   }
00591 
00592   phi = G4UniformRand()*twopi;
00593   currentParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
00594   et = pseudoParticle[0].GetTotalEnergy()/GeV;
00595   dndl[0] = 0.0;
00596   vecMass = currentParticle.GetMass()/GeV;
00597 
00598   FragmentationIntegral(pt, et, aspar, vecMass);
00599 
00600   ran = G4UniformRand()*dndl[19];
00601   l = 1;
00602   while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
00603   x = (G4double(l-1) + G4UniformRand())/19.;
00604   currentParticle.SetMomentum( x*et*GeV );             // set the z-momentum
00605 
00606   if (forwardEnergy < forwardKinetic) {
00607     totalEnergy = vecMass + 0.04*std::fabs(normal());
00608     G4cout << " Not enough forward energy: forwardEnergy = " 
00609            << forwardEnergy << " forwardKinetic = "  
00610            << forwardKinetic << " total energy left = " 
00611            << backwardKEDiff + forwardKEDiff << G4endl;
00612   } else {
00613     totalEnergy = vecMass + forwardEnergy - forwardKinetic;
00614     forwardKinetic = forwardEnergy;
00615   }
00616   currentParticle.SetTotalEnergy( totalEnergy*GeV );
00617   pp = std::sqrt(std::abs( totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
00618   pp1 = currentParticle.GetMomentum().mag();
00619 
00620   if (pp1 < 1.0e-6*GeV) {
00621     G4ThreeVector iso = Isotropic(pp);
00622     currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00623   } else {
00624     currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00625   }
00626   pseudoParticle[4] = pseudoParticle[4] + currentParticle;
00627 
00628   // Current particle now finished
00629 
00630   // Begin target particle
00631 
00632   if (backwardNucleonCount < 18) {
00633     targetParticle.SetSide(-3);
00634     ++backwardNucleonCount;
00635 
00636   } else {
00637     // Set pt and phi values, they are changed somewhat in the iteration loop
00638     // Set mass parameter for lambda fragmentation model
00639 
00640     vecMass = targetParticle.GetMass()/GeV;
00641     ran = -std::log(1.0-G4UniformRand());
00642     aspar = 0.40;
00643     pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
00644     phi = G4UniformRand()*twopi;
00645     targetParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
00646     et = pseudoParticle[1].GetTotalEnergy()/GeV;
00647     outerCounter = 0;
00648     innerCounter = 0;
00649     G4bool marginalEnergy = true;
00650     dndl[0] = 0.0;
00651     G4double xxx = 0.999;
00652     if( extraNucleonCount < 20 ) xxx = 0.95+0.05*extraNucleonCount/20.0;
00653     G4ThreeVector momentum;
00654 
00655     while (++outerCounter < 4) {
00656       FragmentationIntegral(pt, et, aspar, vecMass);
00657 
00658       for (innerCounter = 0; innerCounter < 6; innerCounter++) {
00659         ran = G4UniformRand()*dndl[19];
00660         l = 1;
00661         while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
00662         x = -(G4double(l-1) + G4UniformRand())/19.;
00663         targetParticle.SetMomentum( x*et*GeV );        // set the z-momentum
00664         totalEnergy = std::sqrt(x*et*x*et + pt*pt + vecMass*vecMass);
00665         targetParticle.SetTotalEnergy( totalEnergy*GeV );
00666 
00667         if ((backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy) {
00668           pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00669           backwardKinetic += totalEnergy - vecMass;
00670           outerCounter = 3;                 // leave outer loop
00671           marginalEnergy = false;
00672           break;                            // leave inner loop
00673         }
00674         momentum = targetParticle.GetMomentum();
00675         targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
00676         pt *= 0.9;
00677         dndl[19] *= 0.9;
00678       }
00679     }
00680 
00681     if (marginalEnergy) {
00682       G4cout << " Extra backward kinetic energy = "
00683              << 0.999*backwardEnergy - backwardKinetic << G4endl;
00684       totalEnergy = vecMass + 0.999*backwardEnergy - backwardKinetic;
00685       targetParticle.SetTotalEnergy(totalEnergy*GeV);
00686       pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
00687       targetParticle.SetMomentum(momentum.x()/0.9, momentum.y()/0.9);
00688       pp1 = targetParticle.GetMomentum().mag();
00689       targetParticle.SetMomentum(targetParticle.GetMomentum() * pp/pp1 );
00690       pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00691       backwardKinetic = 0.999*backwardEnergy;
00692     }
00693 
00694   }
00695 
00696   if (backwardEnergy < backwardKinetic) 
00697     G4cout << " Backward Edif = " << backwardEnergy - backwardKinetic << G4endl;
00698   if (forwardEnergy != forwardKinetic) 
00699     G4cout << " Forward Edif = " << forwardEnergy - forwardKinetic << G4endl;
00700 
00701   // Target particle finished.
00702   // Now produce backward nucleons with a cluster model
00703   // ps[2] = CM frame of projectile + target
00704   // ps[3] = sum of projectile + nucleon cluster in lab frame
00705   // ps[6] = proj + cluster 4-vector boosted into CM frame of proj + targ
00706   //         with secondaries, current and target particles subtracted
00707   //       = total 4-momentum of backward nucleon cluster  
00708 
00709   pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
00710   pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
00711   pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
00712 
00713   if (backwardNucleonCount == 1) {   
00714     // Target particle is the only backward nucleon.  Give it the remainder
00715     // of the backward kinetic energy. 
00716 
00717     G4double ekin =
00718       std::min(backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV);
00719 
00720     if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
00721     vecMass = targetParticle.GetMass()/GeV;
00722     totalEnergy = ekin + vecMass;
00723     targetParticle.SetTotalEnergy( totalEnergy*GeV );
00724     pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
00725     pp1 = pseudoParticle[6].GetMomentum().mag();
00726     if (pp1 < 1.0e-6*GeV) { 
00727       G4ThreeVector iso = Isotropic(pp);
00728       targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00729     } else {
00730       targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1));
00731     }
00732     pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00733 
00734   } else if (backwardNucleonCount > 1) {
00735     // Share remaining energy with up to 17 backward nucleons
00736 
00737     G4int tempCount = 5;
00738     if (backwardNucleonCount < 5) tempCount = backwardNucleonCount;
00739     tempCount -= 2;
00740 
00741     G4double clusterMass = 0.;
00742     if (targetParticle.GetSide() == -3) 
00743       clusterMass = targetParticle.GetMass()/GeV;
00744     for (i = 0; i < vecLen; ++i)
00745       if (vec[i]->GetSide() == -3) clusterMass += vec[i]->GetMass()/GeV;
00746     clusterMass += backwardEnergy - backwardKinetic;
00747 
00748     totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
00749     pseudoParticle[6].SetMass(clusterMass*GeV);
00750 
00751     pp = std::sqrt(std::abs(totalEnergy*totalEnergy - 
00752                             clusterMass*clusterMass) )*GeV;
00753     pp1 = pseudoParticle[6].GetMomentum().mag();
00754     if (pp1 < 1.0e-6*GeV) {
00755       G4ThreeVector iso = Isotropic(pp);
00756       pseudoParticle[6].SetMomentum(iso.x(), iso.y(), iso.z());
00757     } else {
00758       pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp/pp1));
00759     }
00760 
00761     std::vector<G4ReactionProduct*> tempList;  // ptrs to backward nucleons
00762     if (targetParticle.GetSide() == -3) tempList.push_back(&targetParticle);
00763     for (i = 0; i < vecLen; ++i)
00764       if (vec[i]->GetSide() == -3) tempList.push_back(vec[i]);
00765 
00766     constantCrossSection = true;
00767 
00768     if (tempList.size() > 1) {
00769       G4int n_entry = 0;
00770       wgt = GenerateNBodyEventT(pseudoParticle[6].GetMass(), 
00771                                 constantCrossSection, tempList);
00772 
00773       if (targetParticle.GetSide() == -3) {
00774         targetParticle = *tempList[0];
00775         targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
00776         n_entry++;
00777       }
00778 
00779       for (i = 0; i < vecLen; ++i) {
00780         if (vec[i]->GetSide() == -3) {
00781           *vec[i] = *tempList[n_entry];
00782           vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
00783           n_entry++;
00784         }
00785       }
00786     }
00787   } else return false;
00788 
00789   if (vecLen == 0) return false;  // all the secondaries have been eliminated
00790 
00791   // Lorentz transformation to lab frame
00792 
00793   currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
00794   targetParticle.Lorentz( targetParticle, pseudoParticle[1] );    
00795   for (i = 0; i < vecLen; ++i) vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
00796 
00797   // Set leading strange particle flag.
00798   // leadFlag will be true if original particle and incident particle are
00799   // both strange, in which case the incident particle becomes the leading 
00800   // particle. 
00801   // leadFlag will also be true if the target particle is strange, but the
00802   // incident particle is not, in which case the target particle becomes the 
00803   // leading particle.
00804 
00805   G4bool leadingStrangeParticleHasChanged = true;
00806   if (leadFlag)
00807   {
00808     if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
00809       leadingStrangeParticleHasChanged = false;
00810     if (leadingStrangeParticleHasChanged &&
00811         (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()) )
00812       leadingStrangeParticleHasChanged = false;
00813     if( leadingStrangeParticleHasChanged )
00814     {
00815       for( i=0; i<vecLen; i++ )
00816       {
00817         if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
00818         {
00819           leadingStrangeParticleHasChanged = false;
00820           break;
00821         }
00822       }
00823     }
00824     if( leadingStrangeParticleHasChanged )
00825     {
00826       G4bool leadTest = 
00827         (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
00828          leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
00829       G4bool targetTest =
00830         (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
00831          targetParticle.GetDefinition()->GetParticleSubType() == "pi");
00832         
00833       // following modified by JLC 22-Oct-97
00834           
00835       if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
00836       {
00837         targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
00838         targetHasChanged = true;
00839       }
00840       else
00841       {
00842         currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
00843         incidentHasChanged = false;
00844       }
00845     }
00846   }   // end of if( leadFlag )
00847 
00848   // Get number of final state nucleons and nucleons remaining in
00849   // target nucleus
00850     
00851   std::pair<G4int, G4int> finalStateNucleons = 
00852     GetFinalStateNucleons(originalTarget, vec, vecLen);
00853 
00854   G4int protonsInFinalState = finalStateNucleons.first;
00855   G4int neutronsInFinalState = finalStateNucleons.second;
00856 
00857   G4int numberofFinalStateNucleons = 
00858     protonsInFinalState + neutronsInFinalState;
00859 
00860   if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00861       targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00862       originalIncident->GetDefinition()->GetPDGMass() < 
00863                                  G4Lambda::Lambda()->GetPDGMass())
00864     numberofFinalStateNucleons++;
00865 
00866   numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
00867 
00868   G4int PinNucleus = std::max(0, 
00869     G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
00870   G4int NinNucleus = std::max(0,
00871     G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
00872 
00873   pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00874   pseudoParticle[3].SetMass( mOriginal*GeV );
00875   pseudoParticle[3].SetTotalEnergy(
00876    std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
00877     
00878   G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
00879   G4int diff = 0;
00880   if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() )  diff = 1;
00881   if(numberofFinalStateNucleons == 1) diff = 0;
00882   pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
00883   pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff) );
00884   pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff) );
00885     
00886   G4double theoreticalKinetic = 
00887     pseudoParticle[3].GetTotalEnergy() + pseudoParticle[4].GetTotalEnergy() -
00888     currentParticle.GetMass() - targetParticle.GetMass();
00889   for (i = 0; i < vecLen; ++i) theoreticalKinetic -= vec[i]->GetMass();
00890     
00891   G4double simulatedKinetic = 
00892     currentParticle.GetKineticEnergy() + targetParticle.GetKineticEnergy();
00893   for (i = 0; i < vecLen; ++i) 
00894     simulatedKinetic += vec[i]->GetKineticEnergy();
00895     
00896   pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
00897   pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
00898   pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
00899     
00900   pseudoParticle[7].SetZero();
00901   pseudoParticle[7] = pseudoParticle[7] + currentParticle;
00902   pseudoParticle[7] = pseudoParticle[7] + targetParticle;
00903   for (i = 0; i < vecLen; ++i)
00904     pseudoParticle[7] = pseudoParticle[7] + *vec[i];
00905 
00906     /*
00907     // This code does not appear to do anything to the energy or angle spectra
00908     if( vecLen <= 16 && vecLen > 0 )
00909     {
00910       // must create a new set of ReactionProducts here because GenerateNBody will
00911       // modify the momenta for the particles, and we don't want to do this
00912       //
00913       G4ReactionProduct tempR[130];
00914       tempR[0] = currentParticle;
00915       tempR[1] = targetParticle;
00916       for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
00917       G4FastVector<G4ReactionProduct,256> tempV1;
00918       tempV1.Initialize( vecLen+2 );
00919       G4int tempLen1 = 0;
00920       for( i=0; i<vecLen+2; ++i )tempV1.SetElement( tempLen1++, &tempR[i] );
00921       constantCrossSection = true;
00922 
00923       wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() +
00924                                pseudoParticle[4].GetTotalEnergy(),
00925                                constantCrossSection, tempV1, tempLen1);
00926       if (wgt == -1) {
00927         G4double Qvalue = 0;
00928         for (i = 0; i < tempLen1; i++) Qvalue += tempV1[i]->GetMass();
00929         wgt = GenerateNBodyEvent(Qvalue,
00930                                  constantCrossSection, tempV1, tempLen1);
00931       }
00932       if(wgt>-.5)
00933       {
00934         theoreticalKinetic = 0.0;
00935         for( i=0; i<tempLen1; ++i )
00936         {
00937           pseudoParticle[6].Lorentz( *tempV1[i], pseudoParticle[4] );
00938           theoreticalKinetic += pseudoParticle[6].GetKineticEnergy();
00939         }
00940       }
00941       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
00942     }
00943     */
00944 
00945   //
00946   // Make sure that the kinetic energies are correct
00947   //
00948 
00949   if (simulatedKinetic != 0.0) {
00950     wgt = theoreticalKinetic/simulatedKinetic;
00951     theoreticalKinetic = currentParticle.GetKineticEnergy() * wgt;
00952     simulatedKinetic = theoreticalKinetic;
00953     currentParticle.SetKineticEnergy(theoreticalKinetic);
00954     pp = currentParticle.GetTotalMomentum();
00955     pp1 = currentParticle.GetMomentum().mag();
00956     if (pp1 < 1.0e-6*GeV) {
00957       G4ThreeVector iso = Isotropic(pp);
00958       currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00959     } else {
00960       currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp/pp1));
00961     }
00962 
00963     theoreticalKinetic = targetParticle.GetKineticEnergy() * wgt;
00964     targetParticle.SetKineticEnergy(theoreticalKinetic);
00965     simulatedKinetic += theoreticalKinetic;
00966     pp = targetParticle.GetTotalMomentum();
00967     pp1 = targetParticle.GetMomentum().mag();
00968 
00969     if (pp1 < 1.0e-6*GeV) {
00970       G4ThreeVector iso = Isotropic(pp);
00971       targetParticle.SetMomentum(iso.x(), iso.y(), iso.z() );
00972     } else {
00973       targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp/pp1) );
00974     }
00975 
00976     for (i = 0; i < vecLen; ++i ) {
00977       theoreticalKinetic = vec[i]->GetKineticEnergy() * wgt;
00978       simulatedKinetic += theoreticalKinetic;
00979       vec[i]->SetKineticEnergy(theoreticalKinetic);
00980       pp = vec[i]->GetTotalMomentum();
00981       pp1 = vec[i]->GetMomentum().mag();
00982       if( pp1 < 1.0e-6*GeV ) {
00983         G4ThreeVector iso = Isotropic(pp);
00984         vec[i]->SetMomentum(iso.x(), iso.y(), iso.z() );
00985       } else {
00986         vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp/pp1) );
00987       }
00988     }
00989   }
00990 
00991   //    Rotate(numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
00992   //           modifiedOriginal, originalIncident, targetNucleus,
00993   //           currentParticle, targetParticle, vec, vecLen );
00994 
00995   // Add black track particles
00996   // the total number of particles produced is restricted to 198
00997   // this may have influence on very high energies
00998 
00999   if( atomicWeight >= 1.5 )
01000   {
01001     // npnb is number of proton/neutron black track particles
01002     // ndta is the number of deuterons, tritons, and alphas produced
01003     // epnb is the kinetic energy available for proton/neutron black track 
01004     //    particles
01005     // edta is the kinetic energy available for deuteron/triton/alpha particles
01006 
01007     G4int npnb = 0;
01008     G4int ndta = 0;
01009       
01010     G4double epnb, edta;
01011     if (veryForward) {
01012       epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
01013       edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
01014     } else {
01015       epnb = targetNucleus.GetPNBlackTrackEnergy();
01016       edta = targetNucleus.GetDTABlackTrackEnergy();
01017     }
01018       /*
01019       G4ReactionProduct* fudge = new G4ReactionProduct();
01020       fudge->SetDefinition( aProton );
01021       G4double TT = epnb + edta;
01022       G4double MM = fudge->GetMass()/GeV;
01023       fudge->SetTotalEnergy(MM*GeV + TT*GeV);
01024       G4double pzz = std::sqrt(TT*(TT + 2.*MM));
01025       fudge->SetMomentum(0.0, 0.0, pzz*GeV);
01026       vec.SetElement(vecLen++, fudge);    
01027       // G4cout << " Fudge = " << vec[vecLen-1]->GetKineticEnergy()/GeV 
01028                 << G4endl;
01029       */
01030 
01031     const G4double pnCutOff = 0.001;
01032     const G4double dtaCutOff = 0.001;
01033     //    const G4double kineticMinimum = 1.e-6;
01034     //    const G4double kineticFactor = -0.010;
01035     //    G4double sprob = 0.0;  // sprob = probability of self-absorption in 
01036                            // heavy molecules
01037     //  Not currently used (DHW 9 June 2008)  const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
01038     //    if (ekIncident >= 5.0) sprob = std::min(1.0, 0.6*std::log(ekIncident-4.0));
01039     if (epnb > pnCutOff)
01040     {
01041       npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
01042       if (numberofFinalStateNucleons + npnb > atomicWeight)
01043         npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
01044       npnb = std::min( npnb, 127-vecLen );
01045     }
01046     if( edta >= dtaCutOff )
01047     {
01048       ndta = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta));
01049       ndta = std::min( ndta, 127-vecLen );
01050     }
01051     if (npnb == 0 && ndta == 0) npnb = 1;
01052 
01053     AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
01054                            PinNucleus, NinNucleus, targetNucleus,
01055                            vec, vecLen);
01056   }
01057 
01058   //  if( centerofmassEnergy <= (4.0+G4UniformRand()) )
01059   //    MomentumCheck( modifiedOriginal, currentParticle, targetParticle, 
01060   //                     vec, vecLen );
01061   //
01062   //  calculate time delay for nuclear reactions
01063   //
01064 
01065   if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
01066     currentParticle.SetTOF( 
01067          1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
01068   else
01069     currentParticle.SetTOF( 1.0 );
01070   return true;
01071 
01072 }
01073 
01074 
01075 void G4RPGFragmentation::
01076 ReduceEnergiesOfSecondaries(G4int startingIndex,
01077                             G4double& forwardKinetic,
01078                             G4double& backwardKinetic,
01079                             G4FastVector<G4ReactionProduct,256>& vec,
01080                             G4int& vecLen,
01081                             G4ReactionProduct& forwardPseudoParticle,
01082                             G4ReactionProduct& backwardPseudoParticle,
01083                             G4double& pt)
01084 {
01085   // Reduce energies and pt of secondaries in order to maintain 
01086   // energy conservation
01087 
01088   G4double totalEnergy;
01089   G4double pp;
01090   G4double pp1;
01091   G4double px;
01092   G4double py;
01093   G4double mass;
01094   G4ReactionProduct* pVec; 
01095   G4int i;
01096 
01097   forwardKinetic = 0.0;
01098   backwardKinetic = 0.0;
01099   forwardPseudoParticle.SetZero();
01100   backwardPseudoParticle.SetZero();
01101 
01102   for (i = startingIndex; i < vecLen; i++) {
01103     pVec = vec[i];
01104     if (pVec->GetSide() != -3) {
01105       mass = pVec->GetMass();
01106       totalEnergy = 0.95*pVec->GetTotalEnergy() + 0.05*mass;
01107       pVec->SetTotalEnergy(totalEnergy);
01108       pp = std::sqrt( std::abs( totalEnergy*totalEnergy - mass*mass ) );
01109       pp1 = pVec->GetMomentum().mag();
01110       if (pp1 < 1.0e-6*GeV) {
01111         G4ThreeVector iso = Isotropic(pp);
01112         pVec->SetMomentum( iso.x(), iso.y(), iso.z() );
01113       } else {
01114         pVec->SetMomentum(pVec->GetMomentum() * (pp/pp1) );
01115       }
01116 
01117       px = pVec->GetMomentum().x();
01118       py = pVec->GetMomentum().y();
01119       pt = std::max(1.0, std::sqrt( px*px + py*py ) )/GeV;
01120       if (pVec->GetSide() > 0) {
01121         forwardKinetic += pVec->GetKineticEnergy()/GeV;
01122         forwardPseudoParticle = forwardPseudoParticle + (*pVec);
01123       } else {
01124         backwardKinetic += pVec->GetKineticEnergy()/GeV;
01125         backwardPseudoParticle = backwardPseudoParticle + (*pVec);
01126       }
01127     }
01128   }
01129 }
01130 
01131  /* end of file */

Generated on Mon May 27 17:49:44 2013 for Geant4 by  doxygen 1.4.7