G4LEAntiNeutronInelastic.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 // Hadronic Process: AntiNeutron Inelastic Process
00029 // J.L. Chuma, TRIUMF, 18-Feb-1997
00030 // J.P.Wellisch: 23-Apr-97: Added theNucleus.SetParameters call
00031 // J.P. Wellisch: 23-Apr-97: nm = np+1; in line 392
00032 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
00033 //
00034  
00035 #include <iostream>
00036 
00037 #include "G4LEAntiNeutronInelastic.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "Randomize.hh"
00041 
00042 G4LEAntiNeutronInelastic::G4LEAntiNeutronInelastic(const G4String& name)
00043  :G4InelasticInteraction(name)
00044 {
00045   SetMinEnergy(0.0);
00046   SetMaxEnergy(25.*GeV);
00047   G4cout << "WARNING: model G4LEAntiNeutronInelastic is being deprecated and will\n"
00048          << "disappear in Geant4 version 10.0"  << G4endl;
00049 }
00050 
00051 
00052 void G4LEAntiNeutronInelastic::ModelDescription(std::ostream& outFile) const
00053 {
00054   outFile << "G4LEAntiNeutronInelastic is one of the Low Energy\n"
00055           << "Parameterized (LEP) models used to implement inelastic\n"
00056           << "anti-neutron scattering from nuclei.  It is a re-engineered\n"
00057           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00058           << "initial collision products into backward- and forward-going\n"
00059           << "clusters which are then decayed into final state hadrons.\n"
00060           << "The model does not conserve energy on an event-by-event\n"
00061           << "basis.  It may be applied to anti-neutrons with initial\n"
00062           << "energies between 0 and 25 GeV.\n";
00063 }
00064 
00065 
00066 G4HadFinalState*
00067 G4LEAntiNeutronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00068                                         G4Nucleus& targetNucleus)
00069 { 
00070   const G4HadProjectile* originalIncident = &aTrack;
00071 
00072   // create the target particle
00073   G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00074     
00075   if (verboseLevel > 1) {
00076     const G4Material *targetMaterial = aTrack.GetMaterial();
00077     G4cout << "G4LEAntiNeutronInelastic::ApplyYourself called" << G4endl;
00078     G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00079     G4cout << "target material = " << targetMaterial->GetName() << ", ";
00080     G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00081            << G4endl;
00082   }
00083 
00084   // Fermi motion and evaporation
00085   // As of Geant3, the Fermi energy calculation had not been Done
00086   G4double ek = originalIncident->GetKineticEnergy()/MeV;
00087   G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00088   G4ReactionProduct modifiedOriginal;
00089   modifiedOriginal = *originalIncident;
00090     
00091   G4double tkin = targetNucleus.Cinema(ek);
00092   ek += tkin;
00093   modifiedOriginal.SetKineticEnergy(ek*MeV);
00094   G4double et = ek + amas;
00095   G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00096   G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00097   if (pp > 0.0) {
00098     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00099     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00100   }
00101 
00102   // calculate black track energies
00103   tkin = targetNucleus.EvaporationEffects(ek);
00104   ek -= tkin;
00105   modifiedOriginal.SetKineticEnergy(ek*MeV);
00106   et = ek + amas;
00107   p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00108   pp = modifiedOriginal.GetMomentum().mag()/MeV;
00109   if (pp > 0.0) {
00110     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00111     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00112   }
00113     
00114   G4ReactionProduct currentParticle = modifiedOriginal;
00115   G4ReactionProduct targetParticle;
00116   targetParticle = *originalTarget;
00117   currentParticle.SetSide(1); // incident always goes in forward hemisphere
00118   targetParticle.SetSide(-1);  // target always goes in backward hemisphere
00119   G4bool incidentHasChanged = false;
00120   G4bool targetHasChanged = false;
00121   G4bool quasiElastic = false;
00122   G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
00123   G4int vecLen = 0;
00124   vec.Initialize(0);
00125     
00126   const G4double cutOff = 0.1*MeV;
00127   const G4double anni = std::min(1.3*currentParticle.GetTotalMomentum()/GeV, 0.4);
00128     
00129   if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
00130       (G4UniformRand() > anni) )
00131     Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00132             incidentHasChanged, targetHasChanged, quasiElastic);
00133   else
00134     quasiElastic = true;
00135     
00136   CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00137                    modifiedOriginal, targetNucleus, currentParticle,
00138                    targetParticle, incidentHasChanged, targetHasChanged,
00139                    quasiElastic);
00140     
00141   SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00142 
00143   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00144 
00145   delete originalTarget;
00146   return &theParticleChange;
00147 }
00148 
00149 
00150 void G4LEAntiNeutronInelastic::Cascade(
00151    G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00152    G4int& vecLen,
00153    const G4HadProjectile *originalIncident,
00154    G4ReactionProduct &currentParticle,
00155    G4ReactionProduct &targetParticle,
00156    G4bool &incidentHasChanged,
00157    G4bool &targetHasChanged,
00158    G4bool &quasiElastic)
00159 {
00160   // Derived from original FORTRAN code CASNB by H. Fesefeldt (13-Sep-1987)
00161   //
00162   // AntiNeutron undergoes interaction with nucleon within a nucleus.  Check if
00163   // it is energetically possible to produce pions/kaons.  In not, assume
00164   // nuclear excitation occurs and input particle is degraded in energy.  No
00165   // other particles are produced.  If reaction is possible, find the correct
00166   // number of pions/protons/neutrons produced using an interpolation to
00167   // multiplicity data.  Replace some pions or protons/neutrons by kaons or
00168   // strange baryons according to the average  multiplicity per inelastic
00169   // reaction.
00170 
00171   const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00172   const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00173   const G4double pOriginal = originalIncident->GetTotalMomentum()/MeV;
00174   const G4double targetMass = targetParticle.GetMass()/MeV;
00175   G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00176                                         targetMass*targetMass +
00177                                         2.0*targetMass*etOriginal );
00178   G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00179     
00180   static G4bool first = true;
00181   const G4int numMul = 1200;
00182   const G4int numMulA = 400;
00183   const G4int numSec = 60;
00184   static G4double protmul[numMul], protnorm[numSec]; // proton constants
00185   static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
00186   static G4double protmulA[numMulA], protnormA[numSec]; // proton constants
00187   static G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
00188 
00189   // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
00190   G4int counter;
00191   G4int nt=0;
00192   G4int npos = 0, nneg = 0, nzero = 0;
00193   G4double test;
00194   const G4double c = 1.25;    
00195   const G4double b[] = { 0.70, 0.70 };
00196 
00197   if (first) {      // compute normalization constants (only be done once)
00198     first = false;
00199     G4int i;
00200     for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00201     for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00202     counter = -1;
00203     for (npos = 0; npos < (numSec/3); ++npos) {
00204       for (nneg = std::max(0,npos-2); nneg <= npos; ++nneg) {
00205         for (nzero = 0; nzero < numSec/3; ++nzero) {
00206           if (++counter < numMul) {
00207             nt = npos+nneg+nzero;
00208             if (nt > 0 && nt <= numSec) {
00209               protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00210               protnorm[nt-1] += protmul[counter];
00211             }
00212           }
00213         }
00214       }
00215     }
00216     for (i = 0; i < numMul; ++i) neutmul[i] = 0.0;
00217     for (i = 0; i < numSec; ++i) neutnorm[i] = 0.0;
00218     counter = -1;
00219     for (npos = 0; npos < numSec/3; ++npos) {
00220       for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
00221         for (nzero = 0; nzero < numSec/3; ++nzero) {
00222           if (++counter < numMul) {
00223             nt = npos+nneg+nzero;
00224             if ( (nt > 0) && (nt <= numSec) ) {
00225               neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00226               neutnorm[nt-1] += neutmul[counter];
00227             }
00228           }
00229         }
00230       }
00231     }
00232     for (i = 0; i < numSec; ++i) {
00233       if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
00234       if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
00235     }
00236 
00237     // do the same for annihilation channels
00238     for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
00239     for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
00240     counter = -1;
00241     for (npos = 1; npos < (numSec/3); ++npos) {
00242       nneg = npos-1;
00243       for (nzero = 0; nzero < numSec/3; ++nzero) {
00244         if (++counter < numMulA) {
00245           nt = npos+nneg+nzero;
00246           if (nt > 1 && nt <= numSec) {
00247             protmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00248             protnormA[nt-1] += protmulA[counter];
00249           }
00250         }
00251       }
00252     }
00253     for (i = 0; i < numMulA; ++i) neutmulA[i] = 0.0;
00254     for (i = 0; i < numSec; ++i) neutnormA[i] = 0.0;
00255     counter = -1;
00256     for (npos = 0; npos < numSec/3; ++npos) {
00257       nneg = npos;
00258       for (nzero = 0; nzero < numSec/3; ++nzero) {
00259         if (++counter < numMulA) {
00260           nt = npos+nneg+nzero;
00261           if (nt > 1 && nt <= numSec) {
00262             neutmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00263             neutnormA[nt-1] += neutmulA[counter];
00264           }
00265         }
00266       }
00267     }
00268     for (i = 0; i < numSec; ++i) {
00269       if (protnormA[i] > 0.0) protnormA[i] = 1.0/protnormA[i];
00270       if (neutnormA[i] > 0.0) neutnormA[i] = 1.0/neutnormA[i];
00271     }
00272   }   // end of initialization
00273 
00274   const G4double expxu = 82.;           // upper bound for arg. of exp
00275   const G4double expxl = -expxu;        // lower bound for arg. of exp
00276   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00277   G4ParticleDefinition *aProton = G4Proton::Proton();
00278   G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
00279   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00280     
00281   // energetically possible to produce pion(s)  -->  inelastic scattering
00282   //                                   otherwise quasi-elastic scattering
00283     
00284   const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
00285                            0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
00286                            0.39,0.36,0.33,0.10,0.01};
00287   G4int iplab = G4int( pOriginal/GeV*10.0 );
00288   if( iplab >  9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
00289   if( iplab > 14 )iplab = G4int(  pOriginal/GeV- 2.0 ) + 15;
00290   if( iplab > 22 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
00291   if( iplab > 24 )iplab = 24;
00292   if (G4UniformRand() > anhl[iplab]) {
00293     if (availableEnergy <= aPiPlus->GetPDGMass()/MeV) {
00294       quasiElastic = true;
00295       return;
00296     }
00297     G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
00298     const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
00299     G4double w0, wp, wt, wm;
00300     if ((availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) ) {
00301       // suppress high multiplicity events at low momentum
00302       // only one pion will be produced
00303       npos = nneg = nzero = 0;
00304       if (targetParticle.GetDefinition() == aProton) {
00305         test = std::exp(std::min(expxu,
00306                         std::max( expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
00307         w0 = test;
00308         wp = test;
00309         if (G4UniformRand() < w0/(w0+wp) ) nzero = 1;
00310         else npos = 1;
00311       } else {
00312         // target is a neutron
00313         test = std::exp(std::min(expxu,
00314                         std::max( expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
00315         w0 = test;
00316         wp = test;
00317         test = std::exp(std::min(expxu,
00318                         std::max( expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
00319         wm = test;
00320         wt = w0+wp+wm;
00321         wp += w0;
00322         G4double ran = G4UniformRand();
00323         if( ran < w0/wt )
00324           nzero = 1;
00325         else if( ran < wp/wt )
00326           npos = 1;
00327         else
00328           nneg = 1;
00329       }
00330     } else {
00331       // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
00332       G4double n, anpn;
00333       GetNormalizationConstant( availableEnergy, n, anpn );
00334       G4double ran = G4UniformRand();
00335       G4double dum, excs = 0.0;
00336       if (targetParticle.GetDefinition() == aProton) {
00337         counter = -1;
00338         for (npos = 0; npos < numSec/3 && ran>=excs; ++npos) {
00339           for (nneg = std::max(0,npos-2); nneg <= npos && ran >= excs; ++nneg) {
00340             for (nzero = 0; nzero < numSec/3 && ran >= excs; ++nzero) {
00341               if (++counter < numMul) {
00342                 nt = npos + nneg + nzero;
00343                 if (nt > 0) {
00344                   test = std::exp(std::min(expxu,
00345                                   std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00346                   dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00347                   if( std::fabs(dum) < 1.0 ) {
00348                     if( test >= 1.0e-10 )excs += dum*test;
00349                   }
00350                     else
00351                       excs += dum*test;
00352                 }
00353               }
00354             }
00355           }
00356         }
00357         if (ran >= excs) {
00358           // 3 previous loops continued to the end
00359           quasiElastic = true;
00360           return;
00361         }
00362         npos--; nneg--; nzero--;
00363 
00364       } else {
00365         // target must be a neutron
00366         counter = -1;
00367         for (npos = 0; npos < numSec/3 && ran >= excs; ++npos) {
00368           for (nneg = std::max(0,npos-1); nneg <= (npos+1) && ran >= excs; ++nneg) {
00369             for (nzero = 0; nzero < numSec/3 && ran>=excs; ++nzero) {
00370               if (++counter < numMul) {
00371                 nt = npos+nneg+nzero;
00372                 if ((nt >= 1) && (nt <= numSec) ) {
00373                   test = std::exp(std::min(expxu,
00374                                   std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00375                   dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00376                   if (std::fabs(dum) < 1.0) {
00377                     if (test >= 1.0e-10) excs += dum*test;
00378                   }
00379                   else
00380                     excs += dum*test;
00381                 }
00382               }
00383             }
00384           }
00385         }
00386         if (ran >= excs) {
00387           // 3 previous loops continued to the end
00388           quasiElastic = true;
00389           return;
00390         }
00391         npos--; nneg--; nzero--;
00392       }
00393     }
00394 
00395     if (targetParticle.GetDefinition() == aProton) {
00396       switch (npos-nneg)
00397       {
00398        case 1:
00399          if( G4UniformRand() < 0.5 )
00400          {
00401            currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00402            incidentHasChanged = true;
00403          }
00404          else
00405          {
00406            targetParticle.SetDefinitionAndUpdateE( aNeutron );
00407            targetHasChanged = true;
00408          }
00409          break;
00410        case 2:
00411          currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00412          targetParticle.SetDefinitionAndUpdateE( aNeutron );
00413          incidentHasChanged = true;
00414          targetHasChanged = true;
00415          break;
00416        default:
00417          break;
00418       }
00419     } else {
00420       // target must be a neutron
00421       switch (npos-nneg)
00422       {
00423         case 0:
00424           if (G4UniformRand() < 0.33) {
00425              currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00426              targetParticle.SetDefinitionAndUpdateE( aProton );
00427              incidentHasChanged = true;
00428              targetHasChanged = true;
00429           }
00430           break;
00431         case 1:
00432           currentParticle.SetDefinitionAndUpdateE( anAntiProton );
00433           incidentHasChanged = true;
00434           break;
00435         default:
00436           targetParticle.SetDefinitionAndUpdateE( aProton );
00437           targetHasChanged = true;
00438           break;
00439       }
00440     }
00441   } else {
00442     // random number <= anhl[iplab]
00443     if (centerofmassEnergy <= 2*aPiPlus->GetPDGMass()/MeV) {
00444       quasiElastic = true;
00445       return;
00446     }
00447 
00448     // annihilation channels
00449     G4double n, anpn;
00450     GetNormalizationConstant(-centerofmassEnergy, n, anpn);
00451     G4double ran = G4UniformRand();
00452     G4double dum, excs = 0.0;
00453     if (targetParticle.GetDefinition() == aProton) {
00454       counter = -1;
00455       for (npos = 1; (npos < numSec/3) && (ran>=excs); ++npos) {
00456         nneg = npos-1;
00457         for (nzero = 0; (nzero < numSec/3) && (ran>=excs); ++nzero) {
00458           if (++counter < numMulA) {
00459             nt = npos+nneg+nzero;
00460             if (nt > 1 && nt <= numSec) {
00461               test = std::exp(std::min(expxu,
00462                               std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00463               dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
00464               if (std::fabs(dum) < 1.0)
00465                 {
00466                   if( test >= 1.0e-10 )excs += dum*test;
00467                 }
00468                 else
00469                   excs += dum*test;
00470               }
00471             }
00472           }
00473         }
00474     } else {
00475       // target must be a neutron
00476       counter = -1;
00477       for (npos = 0; (npos < numSec/3) && (ran >= excs); ++npos) {
00478         nneg = npos;
00479         for (nzero = 0; (nzero < numSec/3) && (ran >= excs); ++nzero) {
00480           if (++counter < numMulA) {
00481             nt = npos+nneg+nzero;
00482             if ((nt > 1) && (nt <= numSec) ) {
00483               test = std::exp(std::min(expxu,
00484                               std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00485               dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
00486               if (std::fabs(dum) < 1.0) {
00487                 if (test >= 1.0e-10 )excs += dum*test;
00488               }
00489                 else
00490                   excs += dum*test;
00491             }
00492           }
00493         }
00494       }
00495     }
00496     if (ran >= excs) {
00497       // 3 previous loops continued to the end
00498       quasiElastic = true;
00499       return;
00500     }
00501     npos--; nzero--;
00502     currentParticle.SetMass( 0.0 );
00503     targetParticle.SetMass( 0.0 );
00504   }
00505 
00506   while(npos+nneg+nzero < 3) nzero++;
00507   SetUpPions(npos, nneg, nzero, vec, vecLen);
00508   return;
00509 }

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