#include <G4HadLeadBias.hh>
Inheritance diagram for G4HadLeadBias:
Public Member Functions | |
virtual G4HadFinalState * | Bias (G4HadFinalState *result) |
virtual | ~G4HadLeadBias () |
Definition at line 36 of file G4HadLeadBias.hh.
virtual G4HadLeadBias::~G4HadLeadBias | ( | ) | [inline, virtual] |
G4HadFinalState * G4HadLeadBias::Bias | ( | G4HadFinalState * | result | ) | [virtual] |
Implements G4VLeadingParticleBiasing.
Definition at line 34 of file G4HadLeadBias.cc.
References G4HadFinalState::AddSecondaries(), buffer, G4HadFinalState::ClearSecondaries(), G4UniformRand, G4Gamma::Gamma(), G4ParticleDefinition::GetBaryonNumber(), G4DynamicParticle::GetDefinition(), G4HadFinalState::GetEnergyChange(), G4DynamicParticle::GetKineticEnergy(), G4ParticleDefinition::GetLeptonNumber(), G4HadFinalState::GetNumberOfSecondaries(), G4HadSecondary::GetParticle(), G4HadFinalState::GetSecondary(), G4HadFinalState::GetStatusChange(), isAlive, and G4PionZero::PionZero().
00035 { 00036 // G4cerr << "bias enter"<<G4endl; 00037 G4int nMeson(0), nBaryon(0), npi0(0), ngamma(0), nLepton(0); 00038 G4int i(0); 00039 G4int maxE = -1; 00040 G4double emax = 0; 00041 if(result->GetStatusChange()==isAlive) 00042 { 00043 emax = result->GetEnergyChange(); 00044 } 00045 //G4cout << "max energy "<<G4endl; 00046 for(i=0;i<result->GetNumberOfSecondaries();i++) 00047 { 00048 if(result->GetSecondary(i)->GetParticle()->GetKineticEnergy()>emax) 00049 { 00050 maxE = i; 00051 emax = result->GetSecondary(i)->GetParticle()->GetKineticEnergy(); 00052 } 00053 } 00054 //G4cout <<"loop1"<<G4endl; 00055 for(i=0; i<result->GetNumberOfSecondaries(); i++) 00056 { 00057 const G4DynamicParticle* aSecTrack = result->GetSecondary(i)->GetParticle(); 00058 if(i==maxE) 00059 { 00060 } 00061 else if(aSecTrack->GetDefinition()->GetBaryonNumber()!=0) 00062 { 00063 nBaryon++; 00064 } 00065 else if(aSecTrack->GetDefinition()->GetLeptonNumber()!=0) 00066 { 00067 nLepton++; 00068 } 00069 else if(aSecTrack->GetDefinition()==G4Gamma::Gamma()) 00070 { 00071 ngamma++; 00072 } 00073 else if(aSecTrack->GetDefinition()==G4PionZero::PionZero()) 00074 { 00075 npi0++; 00076 } 00077 else 00078 { 00079 nMeson++; 00080 } 00081 } 00082 //G4cout << "BiasDebug 1 = "<<result->GetNumberOfSecondaries()<<" " 00083 // <<nMeson<<" "<< nBaryon<<" "<< npi0<<" "<< ngamma<<" "<< nLepton<<G4endl; 00084 G4double mesonWeight = nMeson; 00085 G4double baryonWeight = nBaryon; 00086 G4double gammaWeight = ngamma; 00087 G4double npi0Weight = npi0; 00088 G4double leptonWeight = nLepton; 00089 G4int randomMeson = static_cast<G4int>((nMeson+1)*G4UniformRand()); 00090 G4int randomBaryon = static_cast<G4int>((nBaryon+1)*G4UniformRand()); 00091 G4int randomGamma = static_cast<G4int>((ngamma+1)*G4UniformRand()); 00092 G4int randomPi0 = static_cast<G4int>((npi0+1)*G4UniformRand()); 00093 G4int randomLepton = static_cast<G4int>((nLepton+1)*G4UniformRand()); 00094 00095 std::vector<G4HadSecondary> buffer; 00096 G4int cMeson(0), cBaryon(0), cpi0(0), cgamma(0), cLepton(0); 00097 for(i=0; i<result->GetNumberOfSecondaries(); i++) 00098 { 00099 G4bool aCatch = false; 00100 G4double weight = 1; 00101 const G4HadSecondary* aSecTrack = result->GetSecondary(i); 00102 G4ParticleDefinition* aSecDef = aSecTrack->GetParticle()->GetDefinition(); 00103 if(i==maxE) 00104 { 00105 aCatch = true; 00106 weight = 1; 00107 } 00108 else if(aSecDef->GetBaryonNumber()!=0) 00109 { 00110 if(++cBaryon==randomBaryon) 00111 { 00112 aCatch = true; 00113 weight = baryonWeight; 00114 } 00115 } 00116 else if(aSecDef->GetLeptonNumber()!=0) 00117 { 00118 if(++cLepton==randomLepton) 00119 { 00120 aCatch = true; 00121 weight = leptonWeight; 00122 } 00123 } 00124 else if(aSecDef==G4Gamma::Gamma()) 00125 { 00126 if(++cgamma==randomGamma) 00127 { 00128 aCatch = true; 00129 weight = gammaWeight; 00130 } 00131 } 00132 else if(aSecDef==G4PionZero::PionZero()) 00133 { 00134 if(++cpi0==randomPi0) 00135 { 00136 aCatch = true; 00137 weight = npi0Weight; 00138 } 00139 } 00140 else 00141 { 00142 if(++cMeson==randomMeson) 00143 { 00144 aCatch = true; 00145 weight = mesonWeight; 00146 } 00147 } 00148 if(aCatch) 00149 { 00150 buffer.push_back(*aSecTrack); 00151 buffer.back().SetWeight(aSecTrack->GetWeight()*weight); 00152 } 00153 else 00154 { 00155 delete aSecTrack; 00156 } 00157 } 00158 result->ClearSecondaries(); 00159 // G4cerr << "pre"<<G4endl; 00160 result->AddSecondaries(buffer); 00161 // G4cerr << "bias exit"<<G4endl; 00162 00163 return result; 00164 }