G4HEAntiLambdaInelastic Class Reference

#include <G4HEAntiLambdaInelastic.hh>

Inheritance diagram for G4HEAntiLambdaInelastic:

G4HEInelastic G4HadronicInteraction

Public Member Functions

 G4HEAntiLambdaInelastic (const G4String &name="G4HEAntiLambdaInelastic")
 ~G4HEAntiLambdaInelastic ()
virtual void ModelDescription (std::ostream &) const
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4int GetNumberOfSecondaries ()
void FirstIntInCasAntiLambda (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)

Data Fields

G4int vecLength

Detailed Description

Definition at line 51 of file G4HEAntiLambdaInelastic.hh.


Constructor & Destructor Documentation

G4HEAntiLambdaInelastic::G4HEAntiLambdaInelastic ( const G4String name = "G4HEAntiLambdaInelastic"  ) 

Definition at line 43 of file G4HEAntiLambdaInelastic.cc.

References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.

00044  : G4HEInelastic(name)
00045 {
00046   vecLength = 0;
00047   theMinEnergy = 20*GeV;
00048   theMaxEnergy = 10*TeV;
00049   MAXPART      = 2048;
00050   verboseLevel = 0;
00051   G4cout << "WARNING: model G4HEAntiLambdaInelastic is being deprecated and will\n"
00052          << "disappear in Geant4 version 10.0"  << G4endl;  
00053 }

G4HEAntiLambdaInelastic::~G4HEAntiLambdaInelastic (  )  [inline]

Definition at line 56 of file G4HEAntiLambdaInelastic.hh.

00056 {};


Member Function Documentation

G4HadFinalState * G4HEAntiLambdaInelastic::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus targetNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 71 of file G4HEAntiLambdaInelastic.cc.

References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasAntiLambda(), G4cout, G4endl, G4UniformRand, G4Nucleus::GetA_asInt(), G4HEVector::getCode(), G4HEVector::getEnergy(), G4HEVector::getMass(), G4HEVector::getName(), G4Nucleus::GetZ_asInt(), G4HEInelastic::HighEnergyCascading(), G4HEInelastic::HighEnergyClusterProduction(), G4HEInelastic::MAXPART, G4HEInelastic::MediumEnergyCascading(), G4HEInelastic::MediumEnergyClusterProduction(), G4HEInelastic::NuclearExcitation(), G4HEInelastic::NuclearInelasticity(), G4HEInelastic::QuasiElasticScattering(), G4HEVector::setDefinition(), G4HadFinalState::SetStatusChange(), stopAndKill, G4HEInelastic::StrangeParticlePairProduction(), G4HadronicInteraction::theParticleChange, vecLength, and G4HEInelastic::verboseLevel.

Referenced by G4HEAntiSigmaZeroInelastic::ApplyYourself().

00073 {
00074   G4HEVector* pv = new G4HEVector[MAXPART];
00075   const G4HadProjectile *aParticle = &aTrack;
00076   const G4double atomicWeight = targetNucleus.GetA_asInt();
00077   const G4double atomicNumber = targetNucleus.GetZ_asInt();
00078   G4HEVector incidentParticle(aParticle);
00079 
00080   G4int incidentCode = incidentParticle.getCode();
00081   G4double incidentMass = incidentParticle.getMass();
00082   G4double incidentTotalEnergy = incidentParticle.getEnergy();
00083 
00084   G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00085 
00086   if (incidentKineticEnergy < 1.)
00087     G4cout << "GHEAntiLambdaInelastic: incident energy < 1 GeV" << G4endl;
00088 
00089   if (verboseLevel > 1) {
00090     G4cout << "G4HEAntiLambdaInelastic::ApplyYourself" << G4endl;
00091     G4cout << "incident particle " << incidentParticle.getName()
00092            << "mass "              << incidentMass
00093            << "kinetic energy "    << incidentKineticEnergy
00094            << G4endl;
00095     G4cout << "target material with (A,Z) = (" 
00096            << atomicWeight << "," << atomicNumber << ")" << G4endl;
00097   }
00098 
00099   G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
00100                                               atomicWeight, atomicNumber);
00101   if (verboseLevel > 1)
00102     G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00103     
00104   incidentKineticEnergy -= inelasticity;
00105     
00106   G4double excitationEnergyGNP = 0.;
00107   G4double excitationEnergyDTA = 0.; 
00108 
00109   G4double excitation = NuclearExcitation(incidentKineticEnergy,
00110                                           atomicWeight, atomicNumber,
00111                                           excitationEnergyGNP,
00112                                           excitationEnergyDTA);
00113   if (verboseLevel > 1)
00114     G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
00115            << excitationEnergyDTA << G4endl;             
00116 
00117   incidentKineticEnergy -= excitation;
00118   incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00119   // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
00120   //                        *(incidentTotalEnergy+incidentMass));
00121   // DHW 19 May 2011: variable set but not used
00122 
00123   G4HEVector targetParticle;
00124   if (G4UniformRand() < atomicNumber/atomicWeight) { 
00125     targetParticle.setDefinition("Proton");
00126   } else { 
00127     targetParticle.setDefinition("Neutron");
00128   }
00129 
00130   G4double targetMass = targetParticle.getMass();
00131   G4double centerOfMassEnergy = 
00132        std::sqrt( incidentMass*incidentMass + targetMass*targetMass
00133        + 2.0*targetMass*incidentTotalEnergy);
00134   G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00135 
00136   G4bool inElastic = true;
00137   vecLength = 0;           
00138         
00139   if (verboseLevel > 1)
00140     G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00141            << incidentCode << G4endl;
00142 
00143   G4bool successful = false; 
00144     
00145   FirstIntInCasAntiLambda(inElastic, availableEnergy, pv, vecLength,
00146                           incidentParticle, targetParticle, atomicWeight);
00147 
00148   if (verboseLevel > 1)
00149     G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;  
00150 
00151   if ((vecLength > 0) && (availableEnergy > 1.)) 
00152     StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00153                                   pv, vecLength,
00154                                   incidentParticle, targetParticle);
00155   HighEnergyCascading(successful, pv, vecLength,
00156                       excitationEnergyGNP, excitationEnergyDTA,
00157                       incidentParticle, targetParticle,
00158                       atomicWeight, atomicNumber);
00159   if (!successful)
00160     HighEnergyClusterProduction(successful, pv, vecLength,
00161                                 excitationEnergyGNP, excitationEnergyDTA,
00162                                 incidentParticle, targetParticle,
00163                                 atomicWeight, atomicNumber);
00164   if (!successful) 
00165     MediumEnergyCascading(successful, pv, vecLength, 
00166                           excitationEnergyGNP, excitationEnergyDTA, 
00167                           incidentParticle, targetParticle,
00168                           atomicWeight, atomicNumber);
00169 
00170   if (!successful)
00171     MediumEnergyClusterProduction(successful, pv, vecLength,
00172                                   excitationEnergyGNP, excitationEnergyDTA,       
00173                                   incidentParticle, targetParticle,
00174                                   atomicWeight, atomicNumber);
00175   if (!successful)
00176     QuasiElasticScattering(successful, pv, vecLength,
00177                            excitationEnergyGNP, excitationEnergyDTA,
00178                            incidentParticle, targetParticle, 
00179                            atomicWeight, atomicNumber);
00180   if (!successful)
00181     ElasticScattering(successful, pv, vecLength,
00182                       incidentParticle,    
00183                       atomicWeight, atomicNumber);
00184 
00185   if (!successful)
00186     G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 
00187            << G4endl;
00188 
00189   FillParticleChange(pv, vecLength);
00190   delete [] pv;
00191   theParticleChange.SetStatusChange(stopAndKill);
00192   return & theParticleChange;
00193 }

void G4HEAntiLambdaInelastic::FirstIntInCasAntiLambda ( G4bool inElastic,
const G4double  availableEnergy,
G4HEVector  pv[],
G4int vecLen,
const G4HEVector incidentParticle,
const G4HEVector targetParticle,
const G4double  atomicWeight 
)

Definition at line 197 of file G4HEAntiLambdaInelastic.cc.

References G4HEInelastic::AntiLambda, G4HEInelastic::AntiSigmaMinus, G4HEInelastic::AntiSigmaPlus, G4HEInelastic::AntiSigmaZero, G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), CLHEP::detail::n, G4HEInelastic::Neutron, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, and G4HEInelastic::verboseLevel.

Referenced by ApplyYourself().

00213 {
00214   static const G4double expxu = 82.;     // upper bound for arg. of exp
00215   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00216 
00217   static const G4double protb = 0.7;
00218   static const G4double neutb = 0.7;
00219   static const G4double     c = 1.25;
00220 
00221   static const G4int   numMul   = 1200;
00222   static const G4int   numMulAn = 400;
00223   static const G4int   numSec   = 60;
00224 
00225   G4int protonCode = Proton.getCode();
00226 
00227   G4int targetCode = targetParticle.getCode();
00228   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00229 
00230   static G4bool first = true;
00231   static G4double protmul[numMul],  protnorm[numSec];   // proton constants
00232   static G4double protmulAn[numMulAn],protnormAn[numSec]; 
00233   static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
00234   static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00235 
00236   //  misc. local variables
00237   //  npos = number of pi+,  nneg = number of pi-,  nzero = number of pi0
00238 
00239   G4int i, counter, nt, npos, nneg, nzero;
00240 
00241    if( first ) 
00242      {                 // compute normalization constants, this will only be done once
00243        first = false;
00244        for( i=0; i<numMul  ; i++ ) protmul[i]  = 0.0;
00245        for( i=0; i<numSec  ; i++ ) protnorm[i] = 0.0;
00246        counter = -1;
00247        for( npos=0; npos<(numSec/3); npos++ ) 
00248           {
00249             for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ ) 
00250                {
00251                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00252                     {
00253                       if( ++counter < numMul ) 
00254                         {
00255                           nt = npos+nneg+nzero;
00256                           if( (nt>0) && (nt<=numSec) ) 
00257                             {
00258                               protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00259                               protnorm[nt-1] += protmul[counter];
00260                             }
00261                         }
00262                     }
00263                }
00264           }
00265        for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
00266        for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00267        counter = -1;
00268        for( npos=0; npos<numSec/3; npos++ ) 
00269           {
00270             for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ ) 
00271                {
00272                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00273                     {
00274                       if( ++counter < numMul ) 
00275                         {
00276                           nt = npos+nneg+nzero;
00277                           if( (nt>0) && (nt<=numSec) ) 
00278                             {
00279                                neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00280                                neutnorm[nt-1] += neutmul[counter];
00281                             }
00282                         }
00283                     }
00284                }
00285           }
00286        for( i=0; i<numSec; i++ ) 
00287           {
00288             if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00289             if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00290           }
00291                                                                    // annihilation
00292        for( i=0; i<numMulAn  ; i++ ) protmulAn[i]  = 0.0;
00293        for( i=0; i<numSec    ; i++ ) protnormAn[i] = 0.0;
00294        counter = -1;
00295        for( npos=1; npos<(numSec/3); npos++ ) 
00296           {
00297             nneg = std::max(0,npos-1); 
00298             for( nzero=0; nzero<numSec/3; nzero++ ) 
00299                {
00300                  if( ++counter < numMulAn ) 
00301                    {
00302                      nt = npos+nneg+nzero;
00303                      if( (nt>1) && (nt<=numSec) ) 
00304                        {
00305                          protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00306                          protnormAn[nt-1] += protmulAn[counter];
00307                        }
00308                    }
00309                }
00310           }
00311        for( i=0; i<numMulAn; i++ ) neutmulAn[i]  = 0.0;
00312        for( i=0; i<numSec;   i++ ) neutnormAn[i] = 0.0;
00313        counter = -1;
00314        for( npos=0; npos<numSec/3; npos++ ) 
00315           {
00316             nneg = npos; 
00317             for( nzero=0; nzero<numSec/3; nzero++ ) 
00318                {
00319                  if( ++counter < numMulAn ) 
00320                    {
00321                      nt = npos+nneg+nzero;
00322                      if( (nt>1) && (nt<=numSec) ) 
00323                        {
00324                           neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00325                           neutnormAn[nt-1] += neutmulAn[counter];
00326                        }
00327                    }
00328                }
00329           }
00330        for( i=0; i<numSec; i++ ) 
00331           {
00332             if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00333             if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00334           }
00335      }                                          // end of initialization
00336 
00337          
00338                                               // initialize the first two places
00339                                               // the same as beam and target                                    
00340    pv[0] = incidentParticle;
00341    pv[1] = targetParticle;
00342    vecLen = 2;
00343 
00344    if( !inElastic ) 
00345      {                                        // some two-body reactions 
00346        G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00347 
00348        G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
00349        if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 
00350          {           
00351            G4double ran = G4UniformRand();
00352 
00353            if ( targetCode == protonCode)
00354              {
00355                if(ran < 0.2)
00356                  {
00357                    pv[0] = AntiSigmaZero;
00358                  }
00359                else if (ran < 0.4)
00360                  {
00361                    pv[0] = AntiSigmaMinus;
00362                    pv[1] = Neutron;
00363                  }
00364                else if (ran < 0.6)
00365                  {
00366                    pv[0] = Proton;
00367                    pv[1] = AntiLambda;
00368                  }
00369                else if (ran < 0.8)
00370                  {
00371                    pv[0] = Proton;
00372                    pv[1] = AntiSigmaZero;
00373                  }
00374                else
00375                  {
00376                    pv[0] = Neutron;
00377                    pv[1] = AntiSigmaMinus;
00378                  }     
00379              }
00380            else
00381              {
00382                if (ran < 0.2)
00383                  {
00384                    pv[0] = AntiSigmaZero;
00385                  }
00386                else if (ran < 0.4)
00387                  {
00388                    pv[0] = AntiSigmaPlus;
00389                    pv[1] = Proton;
00390                  }
00391                else if (ran < 0.6)
00392                  {
00393                    pv[0] = Neutron;
00394                    pv[1] = AntiLambda;
00395                  }
00396                else if (ran < 0.8)
00397                  {
00398                    pv[0] = Neutron;
00399                    pv[1] = AntiSigmaZero;
00400                  }
00401                else
00402                  {
00403                    pv[0] = Proton;
00404                    pv[1] = AntiSigmaPlus;
00405                  } 
00406              }  
00407          }   
00408        return;
00409      }
00410    else if (availableEnergy <= PionPlus.getMass())
00411        return;
00412 
00413                                                   //   inelastic scattering
00414 
00415   npos = 0; nneg = 0; nzero = 0;
00416   G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88, 
00417                      0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40, 
00418                      0.39, 0.36, 0.33, 0.10, 0.01};
00419   G4int            iplab =      G4int( incidentTotalMomentum*10.);
00420   if ( iplab >  9) iplab = 10 + G4int( (incidentTotalMomentum  -1.)*5. );          
00421   if ( iplab > 14) iplab = 15 + G4int(  incidentTotalMomentum  -2.     );
00422   if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.); 
00423                    iplab = std::min(24, iplab);
00424 
00425   if (G4UniformRand() > anhl[iplab]) {  // non- annihilation channels
00426 
00427     // number of total particles vs. centre of mass Energy - 2*proton mass
00428     G4double aleab = std::log(availableEnergy);
00429     G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00430                         + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00431    
00432     // normalization constant for kno-distribution.
00433     // calculate first the sum of all constants, check for numerical problems.   
00434     G4double test, dum, anpn = 0.0;
00435 
00436     for (nt = 1; nt <= numSec; nt++) {
00437       test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00438       dum = pi*nt/(2.0*n*n);
00439       if (std::fabs(dum) < 1.0) { 
00440         if (test >= 1.0e-10) anpn += dum*test;
00441       } else { 
00442         anpn += dum*test;
00443       }
00444     }
00445    
00446     G4double ran = G4UniformRand();
00447     G4double excs = 0.0;
00448     if (targetCode == protonCode) {
00449       counter = -1;
00450       for (npos = 0; npos < numSec/3; npos++) {
00451         for (nneg = std::max(0,npos-2); nneg <= (npos+1); nneg++) {
00452           for (nzero = 0; nzero < numSec/3; nzero++) {
00453             if (++counter < numMul) {
00454               nt = npos+nneg+nzero;
00455               if ((nt > 0) && (nt <= numSec) ) {
00456                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00457                 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00458 
00459                 if (std::fabs(dum) < 1.0) { 
00460                   if (test >= 1.0e-10) excs += dum*test;
00461                 } else { 
00462                   excs += dum*test;
00463                 }
00464 
00465                 if (ran < excs) goto outOfLoop;      //----------------------->
00466               }   
00467             }    
00468           }     
00469         }                                                                                  
00470       } 
00471                             // 3 previous loops continued to the end
00472       inElastic = false;    // quasi-elastic scattering   
00473       return;
00474     } else {   // target must be a neutron
00475       counter = -1;
00476                for( npos=0; npos<numSec/3; npos++ ) 
00477                   {
00478                     for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ ) 
00479                        {
00480                          for( nzero=0; nzero<numSec/3; nzero++ ) 
00481                             {
00482                               if( ++counter < numMul ) 
00483                                 {
00484                                   nt = npos+nneg+nzero;
00485                                   if( (nt>0) && (nt<=numSec) ) 
00486                                     {
00487                                       test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00488                                       dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00489                                       if (std::fabs(dum) < 1.0) { 
00490                                         if( test >= 1.0e-10 )excs += dum*test;
00491                                       } else { 
00492                                         excs += dum*test;
00493                                       }
00494 
00495                                       if (ran < excs) goto outOfLoop;       // -------------------------->
00496                                     }
00497                                 }
00498                             }
00499                        }
00500                   }
00501                                                       // 3 previous loops continued to the end
00502                inElastic = false;                     // quasi-elastic scattering.
00503                return;
00504              }
00505           
00506        outOfLoop:           //  <------------------------------------------------------------------------   
00507     
00508        ran = G4UniformRand();
00509 
00510        if( targetCode == protonCode)
00511          {
00512            if( npos == nneg)
00513              {
00514                if (ran < 0.40)
00515                  {
00516                  }
00517                else if (ran < 0.8)
00518                  {
00519                    pv[0] = AntiSigmaZero;
00520                  }
00521                else
00522                  {
00523                    pv[0] = AntiSigmaMinus;
00524                    pv[1] = Neutron;
00525                  } 
00526              }
00527            else if (npos == (nneg+1))
00528              {
00529                if( ran < 0.25)
00530                  {
00531                    pv[1] = Neutron;
00532                  }
00533                else if (ran < 0.5)
00534                  {
00535                    pv[0] = AntiSigmaZero;
00536                    pv[1] = Neutron;
00537                  }
00538                else
00539                  {
00540                    pv[0] = AntiSigmaPlus;
00541                  }
00542              }
00543            else if (npos == (nneg-1))
00544              {  
00545                pv[0] = AntiSigmaMinus;
00546              }
00547            else      
00548              {
00549                pv[0] = AntiSigmaPlus;
00550                pv[1] = Neutron;
00551              } 
00552          }  
00553        else
00554          {
00555            if( npos == nneg)
00556              {
00557                if (ran < 0.4)
00558                  {
00559                  }
00560                else if(ran < 0.8)
00561                  {
00562                    pv[0] = AntiSigmaZero;
00563                  }
00564                else
00565                  {
00566                    pv[0] = AntiSigmaPlus;
00567                    pv[1] = Proton;
00568                  }
00569              } 
00570            else if ( npos == (nneg-1))
00571              {
00572                if (ran < 0.5)
00573                  {
00574                    pv[0] = AntiSigmaMinus;
00575                  }
00576                else if (ran < 0.75)
00577                  {
00578                    pv[1] = Proton;
00579                  }
00580                else
00581                  {
00582                    pv[0] = AntiSigmaZero;
00583                    pv[1] = Proton;
00584                  } 
00585              }
00586            else if (npos == (nneg+1))
00587              {
00588                pv[0] = AntiSigmaPlus;
00589              }
00590            else 
00591              {
00592                pv[0] = AntiSigmaMinus;
00593                pv[1] = Proton;
00594              }
00595          }      
00596      
00597      }
00598    else                                                               // annihilation
00599      {  
00600        if ( availableEnergy > 2. * PionPlus.getMass() )
00601          {
00602 
00603            G4double aleab = std::log(availableEnergy);
00604            G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00605                             + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00606    
00607                       //   normalization constant for kno-distribution.
00608                       //   calculate first the sum of all constants, check for numerical problems.   
00609            G4double test, dum, anpn = 0.0;
00610 
00611            for (nt=2; nt<=numSec; nt++) {
00612              test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00613              dum = pi*nt/(2.0*n*n);
00614 
00615              if (std::fabs(dum) < 1.0) { 
00616                if( test >= 1.0e-10 )anpn += dum*test;
00617              } else { 
00618                anpn += dum*test;
00619              }
00620            }
00621    
00622            G4double ran = G4UniformRand();
00623            G4double excs = 0.0;
00624            if (targetCode == protonCode) {
00625              counter = -1;
00626              for (npos=1; npos<numSec/3; npos++) {
00627                nneg = npos-1; 
00628                for( nzero=0; nzero<numSec/3; nzero++ ) 
00629                  {
00630                    if( ++counter < numMulAn ) 
00631                      {
00632                        nt = npos+nneg+nzero;
00633                        if( (nt>1) && (nt<=numSec) ) {
00634                          test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00635                          dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00636 
00637                          if (std::fabs(dum) < 1.0) { 
00638                            if( test >= 1.0e-10 )excs += dum*test;
00639                          } else { 
00640                            excs += dum*test;
00641                          }
00642 
00643                          if (ran < excs) goto outOfLoopAn;      //----------------------->
00644                        }
00645                      }    
00646                  }   
00647              }                                                                                 
00648                                           // 3 previous loops continued to the end
00649              inElastic = false;           // quasi-elastic scattering   
00650              return;
00651        
00652            } else {                   // target must be a neutron
00653              counter = -1;
00654              for (npos=0; npos<numSec/3; npos++) { 
00655                nneg = npos; 
00656                for( nzero=0; nzero<numSec/3; nzero++ ) {
00657                  if (++counter < numMulAn) {
00658                    nt = npos+nneg+nzero;
00659                    if( (nt>1) && (nt<=numSec) ) {
00660                      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00661                      dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00662 
00663                      if (std::fabs(dum) < 1.0) { 
00664                        if( test >= 1.0e-10 )excs += dum*test;
00665                      } else { 
00666                        excs += dum*test;
00667                      }
00668 
00669                      if (ran < excs) goto outOfLoopAn;       // -------------------------->
00670                    }
00671                  }
00672                }
00673              }
00674 
00675              inElastic = false;              // quasi-elastic scattering.
00676              return;
00677            }
00678        outOfLoopAn:           //  <---------------------------------------------------------   
00679        vecLen = 0;
00680          }
00681      }
00682 
00683    nt = npos + nneg + nzero;
00684    while ( nt > 0)
00685        {
00686          G4double ran = G4UniformRand();
00687          if ( ran < (G4double)npos/nt)
00688             { 
00689               if( npos > 0 ) 
00690                 { pv[vecLen++] = PionPlus;
00691                   npos--;
00692                 }
00693             }
00694          else if ( ran < (G4double)(npos+nneg)/nt)
00695             {   
00696               if( nneg > 0 )
00697                 { 
00698                   pv[vecLen++] = PionMinus;
00699                   nneg--;
00700                 }
00701             }
00702          else
00703             {
00704               if( nzero > 0 )
00705                 { 
00706                   pv[vecLen++] = PionZero;
00707                   nzero--;
00708                 }
00709             }
00710          nt = npos + nneg + nzero;
00711        } 
00712    if (verboseLevel > 1)
00713       {
00714         G4cout << "Particles produced: " ;
00715         G4cout << pv[0].getName() << " " ;
00716         G4cout << pv[1].getName() << " " ;
00717         for (i=2; i < vecLen; i++)   
00718             { 
00719               G4cout << pv[i].getName() << " " ;
00720             }
00721          G4cout << G4endl;
00722       }
00723    return;
00724  }

G4int G4HEAntiLambdaInelastic::GetNumberOfSecondaries (  )  [inline]

Definition at line 65 of file G4HEAntiLambdaInelastic.hh.

References vecLength.

Referenced by G4HEAntiSigmaZeroInelastic::ApplyYourself().

00065 {return vecLength;}

void G4HEAntiLambdaInelastic::ModelDescription ( std::ostream &   )  const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 56 of file G4HEAntiLambdaInelastic.cc.

00057 {
00058   outFile << "G4HEAntiLambdaInelastic is one of the High Energy\n"
00059           << "Parameterized (HEP) models used to implement inelastic\n"
00060           << "anti-Lambda scattering from nuclei.  It is a re-engineered\n"
00061           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00062           << "initial collision products into backward- and forward-going\n"
00063           << "clusters which are then decayed into final state hadrons.\n"
00064           << "The model does not conserve energy on an event-by-event\n"
00065           << "basis.  It may be applied to anti-Lambdas with initial energies\n"
00066           << "above 20 GeV.\n";
00067 }


Field Documentation

G4int G4HEAntiLambdaInelastic::vecLength

Definition at line 60 of file G4HEAntiLambdaInelastic.hh.

Referenced by ApplyYourself(), G4HEAntiLambdaInelastic(), and GetNumberOfSecondaries().


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:10 2013 for Geant4 by  doxygen 1.4.7