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 #include "G4RPGOmegaMinusInelastic.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "Randomize.hh"
00033
00034 G4HadFinalState*
00035 G4RPGOmegaMinusInelastic::ApplyYourself( const G4HadProjectile &aTrack,
00036 G4Nucleus &targetNucleus )
00037 {
00038 const G4HadProjectile *originalIncident = &aTrack;
00039 if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
00040 {
00041 theParticleChange.SetStatusChange(isAlive);
00042 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00043 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00044 return &theParticleChange;
00045 }
00046
00047
00048
00049 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00050
00051 G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
00052
00053 if( verboseLevel > 1 )
00054 {
00055 const G4Material *targetMaterial = aTrack.GetMaterial();
00056 G4cout << "G4RPGOmegaMinusInelastic::ApplyYourself called" << G4endl;
00057 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
00058 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00059 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00060 << G4endl;
00061 }
00062 G4ReactionProduct currentParticle( const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition() ));
00063 currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
00064 currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
00065
00066
00067
00068
00069 G4double ek = originalIncident->GetKineticEnergy();
00070 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00071
00072 G4double tkin = targetNucleus.Cinema( ek );
00073 ek += tkin;
00074 currentParticle.SetKineticEnergy( ek );
00075 G4double et = ek + amas;
00076 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00077 G4double pp = currentParticle.GetMomentum().mag();
00078 if( pp > 0.0 )
00079 {
00080 G4ThreeVector momentum = currentParticle.GetMomentum();
00081 currentParticle.SetMomentum( momentum * (p/pp) );
00082 }
00083
00084
00085
00086 tkin = targetNucleus.EvaporationEffects( ek );
00087 ek -= tkin;
00088 currentParticle.SetKineticEnergy( ek );
00089 et = ek + amas;
00090 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00091 pp = currentParticle.GetMomentum().mag();
00092 if( pp > 0.0 )
00093 {
00094 G4ThreeVector momentum = currentParticle.GetMomentum();
00095 currentParticle.SetMomentum( momentum * (p/pp) );
00096 }
00097
00098 G4ReactionProduct modifiedOriginal = currentParticle;
00099
00100 currentParticle.SetSide( 1 );
00101 targetParticle.SetSide( -1 );
00102 G4bool incidentHasChanged = false;
00103 G4bool targetHasChanged = false;
00104 G4bool quasiElastic = false;
00105 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00106 G4int vecLen = 0;
00107 vec.Initialize( 0 );
00108
00109 const G4double cutOff = 0.1*MeV;
00110 if( currentParticle.GetKineticEnergy() > cutOff )
00111 Cascade( vec, vecLen,
00112 originalIncident, currentParticle, targetParticle,
00113 incidentHasChanged, targetHasChanged, quasiElastic );
00114
00115 CalculateMomenta( vec, vecLen,
00116 originalIncident, originalTarget, modifiedOriginal,
00117 targetNucleus, currentParticle, targetParticle,
00118 incidentHasChanged, targetHasChanged, quasiElastic );
00119
00120 SetUpChange( vec, vecLen,
00121 currentParticle, targetParticle,
00122 incidentHasChanged );
00123
00124 delete originalTarget;
00125 return &theParticleChange;
00126 }
00127
00128
00129 void G4RPGOmegaMinusInelastic::Cascade(
00130 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00131 G4int& vecLen,
00132 const G4HadProjectile *originalIncident,
00133 G4ReactionProduct ¤tParticle,
00134 G4ReactionProduct &targetParticle,
00135 G4bool &incidentHasChanged,
00136 G4bool &targetHasChanged,
00137 G4bool &quasiElastic )
00138 {
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
00149 const G4double etOriginal = originalIncident->GetTotalEnergy();
00150
00151 const G4double targetMass = targetParticle.GetMass();
00152 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00153 targetMass*targetMass +
00154 2.0*targetMass*etOriginal );
00155 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00156 if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass() )
00157 {
00158 quasiElastic = true;
00159 return;
00160 }
00161 static G4bool first = true;
00162 const G4int numMul = 1200;
00163 const G4int numSec = 60;
00164 static G4double protmul[numMul], protnorm[numSec];
00165 static G4double neutmul[numMul], neutnorm[numSec];
00166
00167 G4int counter, nt=0, np=0, nneg=0, nz=0;
00168 G4double test;
00169 const G4double c = 1.25;
00170 const G4double b[] = { 0.70, 0.70 };
00171 if( first )
00172 {
00173 first = false;
00174 G4int i;
00175 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00176 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00177 counter = -1;
00178 for( np=0; np<(numSec/3); ++np )
00179 {
00180 for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
00181 {
00182 for( nz=0; nz<numSec/3; ++nz )
00183 {
00184 if( ++counter < numMul )
00185 {
00186 nt = np+nneg+nz;
00187 if( nt > 0 )
00188 {
00189 protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
00190 protnorm[nt-1] += protmul[counter];
00191 }
00192 }
00193 }
00194 }
00195 }
00196 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00197 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00198 counter = -1;
00199 for( np=0; np<numSec/3; ++np )
00200 {
00201 for( nneg=np; nneg<=(np+2); ++nneg )
00202 {
00203 for( nz=0; nz<numSec/3; ++nz )
00204 {
00205 if( ++counter < numMul )
00206 {
00207 nt = np+nneg+nz;
00208 if( (nt>0) && (nt<=numSec) )
00209 {
00210 neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
00211 neutnorm[nt-1] += neutmul[counter];
00212 }
00213 }
00214 }
00215 }
00216 }
00217 for( i=0; i<numSec; ++i )
00218 {
00219 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00220 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00221 }
00222 }
00223
00224 const G4double expxu = 82.;
00225 const G4double expxl = -expxu;
00226 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00227 G4ParticleDefinition *aProton = G4Proton::Proton();
00228 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00229 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00230 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00231 G4ParticleDefinition *aXiZero = G4XiZero::XiZero();
00232
00233
00234
00235 G4double n, anpn;
00236 GetNormalizationConstant( availableEnergy, n, anpn );
00237 G4double ran = G4UniformRand();
00238 G4double dum, excs = 0.0;
00239 if( targetParticle.GetDefinition() == aProton )
00240 {
00241 counter = -1;
00242 for( np=0; np<numSec/3 && ran>=excs; ++np )
00243 {
00244 for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
00245 {
00246 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00247 {
00248 if( ++counter < numMul )
00249 {
00250 nt = np+nneg+nz;
00251 if( nt > 0 )
00252 {
00253 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00254 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00255 if( std::fabs(dum) < 1.0 )
00256 {
00257 if( test >= 1.0e-10 )excs += dum*test;
00258 }
00259 else
00260 excs += dum*test;
00261 }
00262 }
00263 }
00264 }
00265 }
00266 if( ran >= excs )
00267 {
00268 quasiElastic = true;
00269 return;
00270 }
00271 np--; nneg--; nz--;
00272 }
00273 else
00274 {
00275 counter = -1;
00276 for( np=0; np<numSec/3 && ran>=excs; ++np )
00277 {
00278 for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
00279 {
00280 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00281 {
00282 if( ++counter < numMul )
00283 {
00284 nt = np+nneg+nz;
00285 if( (nt>=1) && (nt<=numSec) )
00286 {
00287 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00288 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00289 if( std::fabs(dum) < 1.0 )
00290 {
00291 if( test >= 1.0e-10 )excs += dum*test;
00292 }
00293 else
00294 excs += dum*test;
00295 }
00296 }
00297 }
00298 }
00299 }
00300 if( ran >= excs )
00301 {
00302 quasiElastic = true;
00303 return;
00304 }
00305 np--; nneg--; nz--;
00306 }
00307
00308
00309
00310
00311
00312 G4int nvefix = 0;
00313 if( targetParticle.GetDefinition() == aProton )
00314 {
00315 if( nneg > np )
00316 {
00317 if( nneg == np+1 )
00318 {
00319 currentParticle.SetDefinitionAndUpdateE( aXiZero );
00320 nvefix = 1;
00321 }
00322 else
00323 {
00324 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00325 nvefix = 2;
00326 }
00327 incidentHasChanged = true;
00328 }
00329 else if( nneg < np )
00330 {
00331 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00332 targetHasChanged = true;
00333 }
00334 }
00335 else
00336 {
00337 if( np+1 < nneg )
00338 {
00339 if( nneg == np+2 )
00340 {
00341 currentParticle.SetDefinitionAndUpdateE( aXiZero );
00342 incidentHasChanged = true;
00343 nvefix = 1;
00344 }
00345 else
00346 {
00347 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00348 incidentHasChanged = true;
00349 nvefix = 2;
00350 }
00351 targetParticle.SetDefinitionAndUpdateE( aProton );
00352 targetHasChanged = true;
00353 }
00354 else if( nneg == np+1 )
00355 {
00356 targetParticle.SetDefinitionAndUpdateE( aProton );
00357 targetHasChanged = true;
00358 }
00359 }
00360
00361 SetUpPions(np, nneg, nz, vec, vecLen);
00362 for (G4int i = 0; i < vecLen && nvefix > 0; ++i) {
00363 if (vec[i]->GetDefinition() == aPiMinus) {
00364 if( nvefix >= 1 )vec[i]->SetDefinitionAndUpdateE(aKaonMinus);
00365 --nvefix;
00366 }
00367 }
00368
00369 return;
00370 }
00371
00372
00373