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
00034
00035
00036
00037
00038 #include "G4MuonVDNuclearModel.hh"
00039
00040 #include "Randomize.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4CascadeInterface.hh"
00044 #include "G4TheoFSGenerator.hh"
00045 #include "G4GeneratorPrecompoundInterface.hh"
00046 #include "G4ExcitationHandler.hh"
00047 #include "G4PreCompoundModel.hh"
00048 #include "G4LundStringFragmentation.hh"
00049 #include "G4ExcitedStringDecay.hh"
00050 #include "G4FTFModel.hh"
00051
00052 G4MuonVDNuclearModel::G4MuonVDNuclearModel()
00053 : G4HadronicInteraction("G4MuonVDNuclearModel")
00054 {
00055 SetMinEnergy(0.0);
00056 SetMaxEnergy(1*PeV);
00057 CutFixed = 0.2*GeV;
00058 NBIN = 1000;
00059
00060 for (G4int k = 0; k < 5; k++) {
00061 for (G4int j = 0; j < 8; j++) {
00062 for (G4int i = 0; i < 1001; i++) {
00063 proba[k][j][i] = 0.0;
00064 ya[i] = 0.0;
00065 }
00066 }
00067 }
00068
00069 MakeSamplingTable();
00070
00071
00072 ftfp = new G4TheoFSGenerator();
00073 precoInterface = new G4GeneratorPrecompoundInterface();
00074 theHandler = new G4ExcitationHandler();
00075 preEquilib = new G4PreCompoundModel(theHandler);
00076 precoInterface->SetDeExcitation(preEquilib);
00077 ftfp->SetTransport(precoInterface);
00078 theFragmentation = new G4LundStringFragmentation();
00079 theStringDecay = new G4ExcitedStringDecay(theFragmentation);
00080 theStringModel = new G4FTFModel;
00081 theStringModel->SetFragmentationModel(theStringDecay);
00082 ftfp->SetHighEnergyGenerator(theStringModel);
00083
00084
00085 bert = new G4CascadeInterface();
00086 }
00087
00088
00089 G4MuonVDNuclearModel::~G4MuonVDNuclearModel()
00090 {
00091 delete ftfp;
00092 delete preEquilib;
00093 delete theFragmentation;
00094 delete theStringDecay;
00095 delete theStringModel;
00096 delete bert;
00097 }
00098
00099
00100 G4HadFinalState*
00101 G4MuonVDNuclearModel::ApplyYourself(const G4HadProjectile& aTrack,
00102 G4Nucleus& targetNucleus)
00103 {
00104 theParticleChange.Clear();
00105
00106
00107 G4double epmax = aTrack.GetTotalEnergy() - 0.5*proton_mass_c2;
00108 if (epmax <= CutFixed) {
00109 theParticleChange.SetStatusChange(isAlive);
00110 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00111 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00112 return &theParticleChange;
00113 }
00114
00115
00116 G4DynamicParticle* transferredPhoton = CalculateEMVertex(aTrack, targetNucleus);
00117
00118
00119 CalculateHadronicVertex(transferredPhoton, targetNucleus);
00120 return &theParticleChange;
00121 }
00122
00123
00124 G4DynamicParticle*
00125 G4MuonVDNuclearModel::CalculateEMVertex(const G4HadProjectile& aTrack,
00126 G4Nucleus& targetNucleus)
00127 {
00128
00129 G4double KineticEnergy = aTrack.GetKineticEnergy();
00130 G4double TotalEnergy = aTrack.GetTotalEnergy();
00131 G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass();
00132 G4double lnZ = std::log(G4double(targetNucleus.GetZ_asInt() ) );
00133
00134 G4double epmin = CutFixed;
00135 G4double epmax = TotalEnergy - 0.5*proton_mass_c2;
00136 G4double m0 = 0.2*GeV;
00137
00138 G4double delmin = 1.e10;
00139 G4double del;
00140 G4int izz = 0;
00141 G4int itt = 0;
00142 G4int NBINminus1 = NBIN - 1;
00143
00144 G4int nzdat = 5;
00145 G4double zdat[] = {1.,4.,13.,29.,92.};
00146 for (G4int iz = 0; iz < nzdat; iz++) {
00147 del = std::abs(lnZ-std::log(zdat[iz]));
00148 if (del < delmin) {
00149 delmin = del;
00150 izz = iz;
00151 }
00152 }
00153
00154 G4int ntdat = 8;
00155 G4double tdat[] = {1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
00156 delmin = 1.e10;
00157 for (G4int it = 0; it < ntdat; it++) {
00158 del = std::abs(std::log(KineticEnergy)-std::log(tdat[it]) );
00159 if (del < delmin) {
00160 delmin = del;
00161 itt = it;
00162 }
00163 }
00164
00165
00166 G4double r = G4UniformRand();
00167
00168 G4int iy = -1;
00169 do {
00170 iy += 1 ;
00171 } while (((proba[izz][itt][iy]) < r)&&(iy < NBINminus1)) ;
00172
00173
00174
00175 G4double y;
00176 if (iy < NBIN)
00177 y = ya[iy] + G4UniformRand() * (ya[iy+1] - ya[iy]);
00178 else
00179 y = ya[iy];
00180
00181 G4double x = std::exp(y);
00182 G4double ep = epmin*std::exp(x*std::log(epmax/epmin) );
00183
00184
00185 G4double yy = ep/TotalEnergy;
00186 G4double tmin = Mass*Mass*yy*yy/(1.-yy);
00187 G4double tmax = 2.*proton_mass_c2*ep;
00188 G4double t1;
00189 G4double t2;
00190 if (m0 < ep) {
00191 t1 = m0*m0;
00192 t2 = ep*ep;
00193 } else {
00194 t1 = ep*ep;
00195 t2 = m0*m0;
00196 }
00197
00198 G4double w1 = tmax*t1;
00199 G4double w2 = tmax+t1;
00200 G4double w3 = tmax*(tmin+t1)/(tmin*w2);
00201 G4double y1 = 1.-yy;
00202 G4double y2 = 0.5*yy*yy;
00203 G4double y3 = y1+y2;
00204
00205 G4double t;
00206 G4double rej;
00207
00208
00209 G4int ntry = 0;
00210 do
00211 {
00212 ntry += 1;
00213 t = w1/(w2*std::exp(G4UniformRand()*std::log(w3))-tmax);
00214 rej = (1.-t/tmax)*(y1*(1.-tmin/t)+y2)/(y3*(1.-t/t2));
00215 } while (G4UniformRand() > rej) ;
00216
00217
00218 G4double sinth2 =
00219 0.5*(t-tmin)/(2.*(TotalEnergy*(TotalEnergy-ep)-Mass*Mass)-tmin);
00220 G4double theta = std::acos(1. - 2.*sinth2);
00221
00222 G4double phi = twopi*G4UniformRand();
00223 G4double sinth = std::sin(theta);
00224 G4double dirx = sinth*std::cos(phi);
00225 G4double diry = sinth*std::sin(phi);
00226 G4double dirz = std::cos(theta);
00227 G4ThreeVector finalDirection(dirx,diry,dirz);
00228 G4ThreeVector ParticleDirection(aTrack.Get4Momentum().vect().unit() );
00229 finalDirection.rotateUz(ParticleDirection);
00230
00231 G4double NewKinEnergy = KineticEnergy - ep;
00232 G4double finalMomentum = std::sqrt(NewKinEnergy*(NewKinEnergy+2.*Mass) );
00233 G4double Ef = NewKinEnergy + Mass;
00234 G4double initMomentum = std::sqrt(KineticEnergy*(TotalEnergy+Mass) );
00235
00236
00237 theParticleChange.SetStatusChange(isAlive);
00238 theParticleChange.SetEnergyChange(NewKinEnergy);
00239 theParticleChange.SetMomentumChange(finalDirection);
00240
00241
00242 G4LorentzVector primaryMomentum(initMomentum*ParticleDirection, TotalEnergy);
00243 G4LorentzVector fsMomentum(finalMomentum*finalDirection, Ef);
00244 G4LorentzVector momentumTransfer = primaryMomentum - fsMomentum;
00245
00246 G4DynamicParticle* gamma =
00247 new G4DynamicParticle(G4Gamma::Gamma(), momentumTransfer);
00248
00249 return gamma;
00250 }
00251
00252
00253 void
00254 G4MuonVDNuclearModel::CalculateHadronicVertex(G4DynamicParticle* incident,
00255 G4Nucleus& target)
00256 {
00257 G4HadFinalState* hfs = 0;
00258 G4double gammaE = incident->GetTotalEnergy();
00259
00260 if (gammaE < 10*GeV) {
00261 G4HadProjectile projectile(*incident);
00262 hfs = bert->ApplyYourself(projectile, target);
00263 } else {
00264
00265 G4double piMass = G4PionZero::PionZero()->GetPDGMass();
00266 G4double piKE = incident->GetTotalEnergy() - piMass;
00267 G4double piMom = std::sqrt(piKE*(piKE + 2*piMass) );
00268 G4ThreeVector piMomentum(incident->GetMomentumDirection() );
00269 piMomentum *= piMom;
00270 G4DynamicParticle theHadron(G4PionZero::PionZero(), piMomentum);
00271 G4HadProjectile projectile(theHadron);
00272 hfs = ftfp->ApplyYourself(projectile, target);
00273 }
00274
00275 delete incident;
00276
00277
00278 theParticleChange.AddSecondaries(hfs);
00279 }
00280
00281
00282 void G4MuonVDNuclearModel::MakeSamplingTable()
00283 {
00284 G4double adat[] = {1.01,9.01,26.98,63.55,238.03};
00285 G4double zdat[] = {1.,4.,13.,29.,92.};
00286 G4int nzdat = 5;
00287
00288 G4double tdat[] = {1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
00289 G4int ntdat = 8;
00290
00291 G4int nbin;
00292 G4double KineticEnergy;
00293 G4double TotalEnergy;
00294 G4double Maxep;
00295 G4double CrossSection;
00296
00297 G4double c;
00298 G4double y;
00299 G4double ymin,ymax;
00300 G4double dy,yy;
00301 G4double dx,x;
00302 G4double ep;
00303
00304 G4double AtomicNumber;
00305 G4double AtomicWeight;
00306
00307 for (G4int iz = 0; iz < nzdat; iz++) {
00308 AtomicNumber = zdat[iz];
00309 AtomicWeight = adat[iz]*(g/mole);
00310
00311 for (G4int it = 0; it < ntdat; it++) {
00312 KineticEnergy = tdat[it];
00313 TotalEnergy = KineticEnergy + G4MuonMinus::MuonMinus()->GetPDGMass();
00314 Maxep = TotalEnergy - 0.5*proton_mass_c2;
00315
00316 CrossSection = 0.0;
00317
00318
00319
00320 c = std::log(Maxep/CutFixed);
00321 ymin = -5.0;
00322 ymax = 0.0;
00323 dy = (ymax-ymin)/NBIN;
00324
00325 nbin=-1;
00326
00327 y = ymin - 0.5*dy;
00328 yy = ymin - dy;
00329 for (G4int i = 0; i < NBIN; i++) {
00330 y += dy;
00331 x = std::exp(y);
00332 yy += dy;
00333 dx = std::exp(yy+dy)-std::exp(yy);
00334
00335 ep = CutFixed*std::exp(c*x);
00336
00337 CrossSection +=
00338 ep*dx*muNucXS.ComputeDDMicroscopicCrossSection(KineticEnergy,
00339 AtomicNumber,
00340 AtomicWeight, ep);
00341 if (nbin < NBIN) {
00342 nbin += 1;
00343 ya[nbin] = y;
00344 proba[iz][it][nbin] = CrossSection;
00345 }
00346 }
00347 ya[NBIN] = 0.;
00348
00349 if (CrossSection > 0.0) {
00350 for (G4int ib = 0; ib <= nbin; ib++) proba[iz][it][ib] /= CrossSection;
00351 }
00352 }
00353 }
00354
00355
00356
00357
00358 }
00359