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