00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "G4HadLeadBias.hh"
00029 #include "G4Gamma.hh"
00030 #include "G4PionZero.hh"
00031 #include "Randomize.hh"
00032 #include "G4HadFinalState.hh"
00033
00034 G4HadFinalState * G4HadLeadBias::Bias(G4HadFinalState * result)
00035 {
00036
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
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
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
00083
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
00160 result->AddSecondaries(buffer);
00161
00162
00163 return result;
00164 }