G4LEAntiXiZeroInelastic.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: AntiXiZero Inelastic Process
00029 // J.L. Chuma, TRIUMF, 20-Feb-1997
00030 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
00031 //
00032 // NOTE:  The FORTRAN version of the cascade, CASAXO, simply called the
00033 //        routine for the XiZero particle.  Hence, the ApplyYourself function
00034 //        below is just a copy of the ApplyYourself from the XiZero particle.
00035  
00036 #include "G4LEAntiXiZeroInelastic.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "Randomize.hh"
00040 
00041 void G4LEAntiXiZeroInelastic::ModelDescription(std::ostream& outFile) const
00042 {
00043   outFile << "G4LEAntiXiZeroInelastic is one of the Low Energy Parameterized\n"
00044           << "(LEP) models used to implement inelastic antiXi0 scattering\n"
00045           << "from nuclei.  It is a re-engineered version of the GHEISHA\n"
00046           << "code of H. Fesefeldt.  It divides the initial collision\n"
00047           << "products into backward- and forward-going clusters which are\n"
00048           << "then decayed into final state hadrons.  The model does not\n"
00049           << "conserve energy on an event-by-event basis.  It may be\n"
00050           << "applied to antiXi0 with initial energies between 0 and 25\n"
00051           << "GeV.\n";
00052 }
00053 
00054 G4HadFinalState*
00055 G4LEAntiXiZeroInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00056                                        G4Nucleus& targetNucleus)
00057 { 
00058   const G4HadProjectile *originalIncident = &aTrack;
00059 
00060   // create the target particle
00061   G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00062     
00063   if (verboseLevel > 1) {
00064     const G4Material *targetMaterial = aTrack.GetMaterial();
00065     G4cout << "G4LEAntiXiZeroInelastic::ApplyYourself called" << G4endl;
00066     G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00067     G4cout << "target material = " << targetMaterial->GetName() << ", ";
00068     G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00069            << G4endl;
00070   }
00071 
00072   // Fermi motion and evaporation
00073   // As of Geant3, the Fermi energy calculation had not been Done
00074   G4double ek = originalIncident->GetKineticEnergy()/MeV;
00075   G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00076   G4ReactionProduct modifiedOriginal;
00077   modifiedOriginal = *originalIncident;
00078     
00079   G4double tkin = targetNucleus.Cinema( ek );
00080   ek += tkin;
00081   modifiedOriginal.SetKineticEnergy( ek*MeV );
00082   G4double et = ek + amas;
00083   G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00084   G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00085   if (pp > 0.0) {
00086     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00087     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00088   }
00089 
00090   // calculate black track energies
00091   tkin = targetNucleus.EvaporationEffects( ek );
00092   ek -= tkin;
00093   modifiedOriginal.SetKineticEnergy( ek*MeV );
00094   et = ek + amas;
00095   p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00096   pp = modifiedOriginal.GetMomentum().mag()/MeV;
00097   if (pp > 0.0) {
00098     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00099     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00100   }
00101   G4ReactionProduct currentParticle = modifiedOriginal;
00102   G4ReactionProduct targetParticle;
00103   targetParticle = *originalTarget;
00104   currentParticle.SetSide(1); // incident always goes in forward hemisphere
00105   targetParticle.SetSide(-1);  // target always goes in backward hemisphere
00106   G4bool incidentHasChanged = false;
00107   G4bool targetHasChanged = false;
00108   G4bool quasiElastic = false;
00109   G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
00110   G4int vecLen = 0;
00111   vec.Initialize(0);
00112     
00113   const G4double cutOff = 0.1;
00114   const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
00115   if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
00116       (G4UniformRand() > anni) )
00117     Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00118             incidentHasChanged, targetHasChanged, quasiElastic);
00119     
00120   CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00121                    modifiedOriginal, targetNucleus, currentParticle,
00122                    targetParticle, incidentHasChanged, targetHasChanged,
00123                    quasiElastic);
00124     
00125   SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00126 
00127   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00128 
00129   delete originalTarget;
00130   return &theParticleChange;
00131 }
00132 
00133  
00134 void G4LEAntiXiZeroInelastic::Cascade(
00135    G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00136    G4int& vecLen,
00137    const G4HadProjectile *originalIncident,
00138    G4ReactionProduct &currentParticle,
00139    G4ReactionProduct &targetParticle,
00140    G4bool &incidentHasChanged,
00141    G4bool &targetHasChanged,
00142    G4bool &quasiElastic)
00143 {
00144     // derived from original FORTRAN code CASAX0 by H. Fesefeldt (20-Jan-1989)
00145     // which is just a copy of CASX0 (cascade for Xi0)
00146     //
00147     // AntiXiZero undergoes interaction with nucleon within a nucleus.  Check if it is
00148     // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
00149     // occurs and input particle is degraded in energy. No other particles are produced.
00150     // If reaction is possible, find the correct number of pions/protons/neutrons
00151     // produced using an interpolation to multiplicity data.  Replace some pions or
00152     // protons/neutrons by kaons or strange baryons according to the average
00153     // multiplicity per inelastic reaction.
00154     //
00155     const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00156     const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00157     const G4double targetMass = targetParticle.GetMass()/MeV;
00158     G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00159                                         targetMass*targetMass +
00160                                         2.0*targetMass*etOriginal );
00161     G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00162     if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV )
00163     {
00164       quasiElastic = true;
00165       return;
00166     }
00167     static G4bool first = true;
00168     const G4int numMul = 1200;
00169     const G4int numSec = 60;
00170     static G4double protmul[numMul], protnorm[numSec]; // proton constants
00171     static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
00172     // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
00173     G4int counter, nt=0, npos=0, nneg=0, nzero=0;
00174     G4double test;
00175     const G4double c = 1.25;    
00176     const G4double b[] = { 0.7, 0.7 };
00177     if( first )       // compute normalization constants, this will only be Done once
00178     {
00179       first = false;
00180       G4int i;
00181       for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00182       for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00183       counter = -1;
00184       for( npos=0; npos<(numSec/3); ++npos )
00185       {
00186         for( nneg=std::max(0,npos-2); nneg<=(npos+1); ++nneg )
00187         {
00188           for( nzero=0; nzero<numSec/3; ++nzero )
00189           {
00190             if( ++counter < numMul )
00191             {
00192               nt = npos+nneg+nzero;
00193               if( nt>0 && nt<=numSec )
00194               {
00195                 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00196                 protnorm[nt-1] += protmul[counter];
00197               }
00198             }
00199           }
00200         }
00201       }
00202       for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00203       for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00204       counter = -1;
00205       for( npos=0; npos<numSec/3; ++npos )
00206       {
00207         for( nneg=std::max(0,npos-1); nneg<=(npos+2); ++nneg )
00208         {
00209           for( nzero=0; nzero<numSec/3; ++nzero )
00210           {
00211             if( ++counter < numMul )
00212             {
00213               nt = npos+nneg+nzero;
00214               if( nt>0 && nt<=numSec )
00215               {
00216                 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00217                 neutnorm[nt-1] += neutmul[counter];
00218               }
00219             }
00220           }
00221         }
00222       }
00223       for( i=0; i<numSec; ++i )
00224       {
00225         if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00226         if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00227       }
00228     }   // end of initialization
00229     
00230     const G4double expxu = 82.;           // upper bound for arg. of exp
00231     const G4double expxl = -expxu;        // lower bound for arg. of exp
00232     G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00233     G4ParticleDefinition *aProton = G4Proton::Proton();
00234     G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00235     G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00236     G4ParticleDefinition *aXiMinus = G4XiMinus::XiMinus();
00237     //
00238     // energetically possible to produce pion(s)  -->  inelastic scattering
00239     //
00240     G4double n, anpn;
00241     GetNormalizationConstant( availableEnergy, n, anpn );
00242     G4double ran = G4UniformRand();
00243     G4double dum, excs = 0.0;
00244     if( targetParticle.GetDefinition() == aProton )
00245     {
00246       counter = -1;
00247       for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00248       {
00249         for( nneg=std::max(0,npos-2); nneg<=(npos+1) && ran>=excs; ++nneg )
00250         {
00251           for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00252           {
00253             if( ++counter < numMul )
00254             {
00255               nt = npos+nneg+nzero;
00256               if( nt>0 && nt<=numSec )
00257               {
00258                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00259                 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00260                 if( std::fabs(dum) < 1.0 )
00261                 {
00262                   if( test >= 1.0e-10 )excs += dum*test;
00263                 }
00264                 else
00265                   excs += dum*test;
00266               }
00267             }
00268           }
00269         }
00270       }
00271       if( ran >= excs )  // 3 previous loops continued to the end
00272       {
00273         quasiElastic = true;
00274         return;
00275       }
00276       npos--; nneg--; nzero--;
00277       //
00278       // number of secondary mesons determined by kno distribution
00279       // check for total charge of final state mesons to determine
00280       // the kind of baryons to be produced, taking into account
00281       // charge and strangeness conservation
00282       //
00283       if( npos < nneg+1 )
00284       {
00285         if( npos != nneg )   // charge mismatch
00286         {
00287           currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00288           incidentHasChanged = true;
00289           //
00290           // correct the strangeness by replacing a pi- by a kaon-
00291           //
00292           vec.Initialize( 1 );
00293           G4ReactionProduct *p = new G4ReactionProduct;
00294           p->SetDefinition( aKaonMinus );
00295           (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00296           vec.SetElement( vecLen++, p );
00297           --nneg;
00298         }
00299       }
00300       else if( npos == nneg+1 )
00301       {
00302         if( G4UniformRand() < 0.5 )
00303         {
00304           targetParticle.SetDefinitionAndUpdateE( aNeutron );
00305           targetHasChanged = true;
00306         }
00307         else
00308         {
00309           currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00310           incidentHasChanged = true;
00311         }
00312       }
00313       else
00314       {
00315         currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00316         incidentHasChanged = true;
00317         targetParticle.SetDefinitionAndUpdateE( aNeutron );
00318         targetHasChanged = true;
00319       }
00320     }
00321     else  // target must be a neutron
00322     {
00323       counter = -1;
00324       for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00325       {
00326         for( nneg=std::max(0,npos-1); nneg<=(npos+2) && ran>=excs; ++nneg )
00327         {
00328           for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00329           {
00330             if( ++counter < numMul )
00331             {
00332               nt = npos+nneg+nzero;
00333               if( nt>0 && nt<=numSec )
00334               {
00335                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00336                 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00337                 if( std::fabs(dum) < 1.0 )
00338                 {
00339                   if( test >= 1.0e-10 )excs += dum*test;
00340                 }
00341                 else
00342                   excs += dum*test;
00343               }
00344             }
00345           }
00346         }
00347       }
00348       if( ran >= excs )  // 3 previous loops continued to the end
00349       {
00350         quasiElastic = true;
00351         return;
00352       }
00353       npos--; nneg--; nzero--;
00354       if( npos < nneg )
00355       {
00356         if( npos+1 == nneg )
00357         {
00358           targetParticle.SetDefinitionAndUpdateE( aProton );
00359           targetHasChanged = true;
00360         }
00361         else   // charge mismatch
00362         {
00363           currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00364           incidentHasChanged = true;
00365           targetParticle.SetDefinitionAndUpdateE( aProton );
00366           targetHasChanged = true;
00367           //
00368           // correct the strangeness by replacing a pi- by a kaon-
00369           //
00370           vec.Initialize( 1 );
00371           G4ReactionProduct *p = new G4ReactionProduct;
00372           p->SetDefinition( aKaonMinus );
00373           (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00374           vec.SetElement( vecLen++, p );
00375           --nneg;
00376         }
00377       }
00378       else if( npos == nneg )
00379       {
00380         if( G4UniformRand() >= 0.5 )
00381         {
00382           currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00383           incidentHasChanged = true;
00384           targetParticle.SetDefinitionAndUpdateE( aProton );
00385           targetHasChanged = true;
00386         }
00387       }
00388       else
00389       {
00390         currentParticle.SetDefinitionAndUpdateE( aXiMinus );
00391         incidentHasChanged = true;
00392       }        
00393     }
00394     SetUpPions( npos, nneg, nzero, vec, vecLen );
00395     return;
00396 }
00397 
00398  /* end of file */
00399  

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