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
00029
00030
00031 #include "G4NeutronHPFinalState.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "G4ParticleTable.hh"
00035 #include "G4Gamma.hh"
00036 #include "G4Neutron.hh"
00037
00038 void G4NeutronHPFinalState::adjust_final_state ( G4LorentzVector init_4p_lab )
00039 {
00040
00041 G4double minimum_energy = 1*keV;
00042
00043 if ( adjustResult != true ) return;
00044
00045 G4int nSecondaries = theResult.GetNumberOfSecondaries();
00046
00047 G4int sum_Z = 0;
00048 G4int sum_A = 0;
00049 G4int max_SecZ = 0;
00050 G4int max_SecA = 0;
00051 G4int imaxA = -1;
00052 for ( int i = 0 ; i < nSecondaries ; i++ )
00053 {
00054 sum_Z += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber();
00055 max_SecZ = std::max ( max_SecZ , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() );
00056 sum_A += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass();
00057 max_SecA = std::max ( max_SecA , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() );
00058 if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA ) imaxA = i;
00059 }
00060
00061 G4ParticleDefinition* resi_pd = NULL;
00062 G4bool needOneMoreSec = false;
00063 G4ParticleDefinition* oneMoreSec_pd = NULL;
00064 if ( (int)(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) == 0 )
00065 {
00066
00067 resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition();
00068 }
00069 else
00070 {
00071 if ( max_SecA > int(theBaseA + 1 - sum_A) )
00072 {
00073
00074 resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition();
00075 needOneMoreSec = true;
00076 }
00077 else
00078 {
00079
00080 resi_pd = G4ParticleTable::GetParticleTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
00081 }
00082
00083 if ( needOneMoreSec )
00084 {
00085 if ( int(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) > 0 )
00086 {
00087 if ( int(theBaseA + 1 - sum_A) > 1 ) G4cout << "More than one neutron is required for the balance of baryon number!" << G4endl;
00088 oneMoreSec_pd = G4Neutron::Neutron();
00089 }
00090 else
00091 {
00092 oneMoreSec_pd = G4ParticleTable::GetParticleTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
00093 }
00094 }
00095
00096 if ( resi_pd == NULL )
00097 {
00098
00099 if ( (int)(theNDLDataZ - sum_Z) == 0 && (int)(theNDLDataA + 1 - sum_A) == 0 )
00100 {
00101 G4int dif_Z = ( int ) ( theNDLDataZ - theBaseZ );
00102 G4int dif_A = ( int ) ( theNDLDataA - theBaseA );
00103 resi_pd = G4ParticleTable::GetParticleTable()->GetIon ( max_SecZ - dif_Z , max_SecA - dif_A , 0.0 );
00104 for ( int i = 0 ; i < nSecondaries ; i++ )
00105 {
00106 if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() == max_SecZ
00107 && theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA )
00108 {
00109 G4ThreeVector p = theResult.GetSecondary( i )->GetParticle()->GetMomentum();
00110 p = p * resi_pd->GetPDGMass()/ G4ParticleTable::GetParticleTable()->GetIon ( max_SecZ , max_SecA , 0.0 )->GetPDGMass();
00111 theResult.GetSecondary( i )->GetParticle()->SetDefinition( resi_pd );
00112 theResult.GetSecondary( i )->GetParticle()->SetMomentum( p );
00113 }
00114 }
00115 }
00116 }
00117 }
00118
00119
00120 G4LorentzVector secs_4p_lab( 0.0 );
00121
00122 G4int n_sec = theResult.GetNumberOfSecondaries();
00123 G4double fast = 0;
00124 G4double slow = 1;
00125 G4int ifast = 0;
00126 G4int islow = 0;
00127 G4int ires = -1;
00128
00129 for ( G4int i = 0 ; i < n_sec ; i++ )
00130 {
00131
00132
00133
00134
00135
00136
00137
00138 secs_4p_lab += theResult.GetSecondary( i )->GetParticle()->Get4Momentum();
00139
00140 G4double beta = 0;
00141 if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() != G4Gamma::Gamma() )
00142 {
00143 beta = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().beta();
00144 }
00145 else
00146 {
00147 beta = 1;
00148 }
00149
00150 if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() == resi_pd ) ires = i;
00151
00152 if ( slow > beta && beta != 0 )
00153 {
00154 slow = beta;
00155 islow = i;
00156 }
00157
00158 if ( fast <= beta )
00159 {
00160 if ( fast != 1 )
00161 {
00162 fast = beta;
00163 ifast = i;
00164 }
00165 else
00166 {
00167
00168 G4double e = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e();
00169 if ( e > theResult.GetSecondary( ifast )->GetParticle()->Get4Momentum().e() )
00170 {
00171
00172 ifast = i;
00173 }
00174 }
00175 }
00176 }
00177
00178
00179 G4LorentzVector dif_4p = init_4p_lab - secs_4p_lab;
00180
00181
00182
00183
00184
00185 G4LorentzVector p4(0);
00186 if ( ires == -1 )
00187 {
00188
00189 ires = nSecondaries;
00190 nSecondaries += 1;
00191
00192 G4DynamicParticle* res = new G4DynamicParticle ( resi_pd , dif_4p.v() );
00193 theResult.AddSecondary ( res );
00194
00195 p4 = res->Get4Momentum();
00196 if ( slow > p4.beta() )
00197 {
00198 slow = p4.beta();
00199 islow = ires;
00200 }
00201 dif_4p = init_4p_lab - ( secs_4p_lab + p4 );
00202 }
00203
00204 if ( needOneMoreSec )
00205 {
00206 nSecondaries += 1;
00207 G4DynamicParticle* one = new G4DynamicParticle ( oneMoreSec_pd , dif_4p.v() );
00208 theResult.AddSecondary ( one );
00209 p4 = one->Get4Momentum();
00210 if ( slow > p4.beta() )
00211 {
00212 slow = p4.beta();
00213 islow = nSecondaries-1;
00214 }
00215 dif_4p = init_4p_lab - ( secs_4p_lab + p4 );
00216 }
00217
00218
00219
00220 if ( dif_4p.v().mag() < std::abs( dif_4p.e() ) )
00221 {
00222
00223
00224
00225 if ( minimum_energy < dif_4p.v().mag() && dif_4p.v().mag() < 1*MeV )
00226 {
00227
00228 nSecondaries += 1;
00229 theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , dif_4p.v() ) );
00230
00231 }
00232 else
00233 {
00234
00235 }
00236
00237 }
00238 else
00239 {
00240
00241
00242
00243
00244
00245 p4 = theResult.GetSecondary( ires )->GetParticle()->Get4Momentum();
00246 theResult.GetSecondary( ires )->GetParticle()->SetMomentum( p4.v() + dif_4p.v() );
00247 dif_4p = init_4p_lab - ( secs_4p_lab - p4 + theResult.GetSecondary( ires )->GetParticle()->Get4Momentum() );
00248
00249
00250
00251 }
00252
00253 G4double dif_e = dif_4p.e() - ( dif_4p.v() ).mag();
00254
00255
00256 if ( dif_e > 0 )
00257 {
00258
00259
00260
00261 nSecondaries += 2;
00262 G4double e1 = ( dif_4p.e() -dif_4p.v().mag() ) / 2;
00263
00264 if ( minimum_energy < e1 )
00265 {
00266 G4double costh = 2.*G4UniformRand()-1.;
00267 G4double phi = twopi*G4UniformRand();
00268 G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi),
00269 std::sin(std::acos(costh))*std::sin(phi),
00270 costh);
00271 theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );
00272 theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );
00273 }
00274 else
00275 {
00276
00277 }
00278
00279 }
00280 else
00281 {
00282
00283
00284 G4double ke0 = theResult.GetSecondary( ifast )->GetParticle()->GetKineticEnergy();
00285 G4ThreeVector p0 = theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();
00286 G4ThreeVector dir = ( theResult.GetSecondary( ifast )->GetParticle()->GetMomentum() ).unit();
00287
00288
00289
00290 if ( ke0 + dif_e > 0 )
00291 {
00292 theResult.GetSecondary( ifast )->GetParticle()->SetKineticEnergy( ke0 + dif_e );
00293 G4ThreeVector dp = p0 - theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();
00294
00295 G4ThreeVector p = theResult.GetSecondary( islow )->GetParticle()->GetMomentum();
00296
00297 theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p + dp );
00298 }
00299 else
00300 {
00301
00302 }
00303
00304 }
00305
00306 }