G4RPGPionSuppression.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 "G4RPGPionSuppression.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "Randomize.hh"
00035 #include "G4HadReentrentException.hh"
00036 
00037 G4RPGPionSuppression::G4RPGPionSuppression()
00038   : G4RPGReaction() {}
00039 
00040 
00041 G4bool G4RPGPionSuppression::
00042 ReactionStage(const G4HadProjectile* /*originalIncident*/,
00043               G4ReactionProduct& modifiedOriginal,
00044               G4bool& incidentHasChanged,
00045               const G4DynamicParticle* /*originalTarget*/,
00046               G4ReactionProduct& targetParticle,
00047               G4bool& targetHasChanged,
00048               const G4Nucleus& targetNucleus,
00049               G4ReactionProduct& currentParticle,
00050               G4FastVector<G4ReactionProduct,256>& vec,
00051               G4int& vecLen,
00052               G4bool /*leadFlag*/,
00053               G4ReactionProduct& /*leadingStrangeParticle*/)
00054 {
00055   // This code was originally in the fortran code TWOCLU
00056   //
00057   // Suppress charged pions, for various reasons
00058   //
00059   G4double mOriginal = modifiedOriginal.GetMass()/GeV;
00060   G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
00061   G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
00062   G4double cmEnergy = std::sqrt( mOriginal*mOriginal + targetMass*targetMass +
00063                            2.0*targetMass*etOriginal ); 
00064   G4double eAvailable = cmEnergy - mOriginal - targetMass;
00065   for (G4int i = 0; i < vecLen; i++) eAvailable -= vec[i]->GetMass()/GeV;
00066 
00067   const G4double atomicWeight = targetNucleus.GetA_asInt();
00068   const G4double atomicNumber = targetNucleus.GetZ_asInt();
00069   const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/GeV;
00070     
00071   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00072   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00073   G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
00074   G4ParticleDefinition *aProton = G4Proton::Proton();
00075   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00076   G4double piMass = aPiPlus->GetPDGMass()/GeV;
00077   G4double nucleonMass = aNeutron->GetPDGMass()/GeV;
00078     
00079   const G4bool antiTest =
00080     modifiedOriginal.GetDefinition() != G4AntiProton::AntiProton() &&
00081     modifiedOriginal.GetDefinition() != G4AntiNeutron::AntiNeutron() &&
00082     modifiedOriginal.GetDefinition() != G4AntiLambda::AntiLambda() &&
00083     modifiedOriginal.GetDefinition() != G4AntiSigmaPlus::AntiSigmaPlus() &&
00084     modifiedOriginal.GetDefinition() != G4AntiSigmaMinus::AntiSigmaMinus() &&
00085     modifiedOriginal.GetDefinition() != G4AntiXiZero::AntiXiZero() &&
00086     modifiedOriginal.GetDefinition() != G4AntiXiMinus::AntiXiMinus() &&
00087     modifiedOriginal.GetDefinition() != G4AntiOmegaMinus::AntiOmegaMinus();
00088 
00089   if( antiTest && (
00090         currentParticle.GetDefinition() == aPiPlus ||
00091         currentParticle.GetDefinition() == aPiZero ||
00092         currentParticle.GetDefinition() == aPiMinus ) &&
00093       ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
00094       ( G4UniformRand() <= atomicWeight/300.0 ) )
00095   {
00096     if (eAvailable > nucleonMass - piMass) {
00097       if( G4UniformRand() > atomicNumber/atomicWeight )
00098         currentParticle.SetDefinitionAndUpdateE( aNeutron );
00099       else
00100         currentParticle.SetDefinitionAndUpdateE( aProton );
00101       incidentHasChanged = true;
00102     }
00103   }
00104 
00105   if( antiTest && (
00106         targetParticle.GetDefinition() == aPiPlus ||
00107         targetParticle.GetDefinition() == aPiZero ||
00108         targetParticle.GetDefinition() == aPiMinus ) &&
00109       ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
00110       ( G4UniformRand() <= atomicWeight/300.0 ) )
00111   {
00112     if (eAvailable > nucleonMass - piMass) {
00113       if( G4UniformRand() > atomicNumber/atomicWeight )
00114         targetParticle.SetDefinitionAndUpdateE( aNeutron );
00115       else
00116         targetParticle.SetDefinitionAndUpdateE( aProton );
00117       targetHasChanged = true;
00118     }
00119   }
00120 
00121   for( G4int i=0; i<vecLen; ++i )
00122   {
00123     if( antiTest && (
00124           vec[i]->GetDefinition() == aPiPlus ||
00125           vec[i]->GetDefinition() == aPiZero ||
00126           vec[i]->GetDefinition() == aPiMinus ) &&
00127         ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
00128         ( G4UniformRand() <= atomicWeight/300.0 ) )
00129     {
00130       if (eAvailable > nucleonMass - piMass) {
00131         if( G4UniformRand() > atomicNumber/atomicWeight )
00132           vec[i]->SetDefinitionAndUpdateE( aNeutron );
00133         else
00134           vec[i]->SetDefinitionAndUpdateE( aProton );
00135       }
00136     }
00137   }
00138 
00139   return true;
00140 }
00141 
00142  
00143  /* end of file */

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