G4MuBremsstrahlungModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4MuBremsstrahlungModel
00034 //
00035 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
00036 //
00037 // Creation date: 24.06.2002
00038 //
00039 // Modifications:
00040 //
00041 // 04-12-02 Change G4DynamicParticle constructor in PostStepDoIt (V.Ivanchenko)
00042 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
00043 // 24-01-03 Fix for compounds (V.Ivanchenko)
00044 // 27-01-03 Make models region aware (V.Ivanchenko)
00045 // 13-02-03 Add name (V.Ivanchenko)
00046 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
00047 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
00048 // 03-08-05 Angular correlations according to PRM (V.Ivanchenko)
00049 // 13-02-06 add ComputeCrossSectionPerAtom (mma)
00050 // 21-03-06 Fix problem of initialisation in case when cuts are not defined (VI)
00051 // 07-11-07 Improve sampling of final state (A.Bogdanov)
00052 // 28-02-08 Use precomputed Z^1/3 and Log(A) (V.Ivanchenko)
00053 //
00054 
00055 //
00056 // Class Description:
00057 //
00058 //
00059 // -------------------------------------------------------------------
00060 //
00061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00063 
00064 #include "G4MuBremsstrahlungModel.hh"
00065 #include "G4PhysicalConstants.hh"
00066 #include "G4SystemOfUnits.hh"
00067 #include "G4Gamma.hh"
00068 #include "G4MuonMinus.hh"
00069 #include "G4MuonPlus.hh"
00070 #include "Randomize.hh"
00071 #include "G4Material.hh"
00072 #include "G4Element.hh"
00073 #include "G4ElementVector.hh"
00074 #include "G4ProductionCutsTable.hh"
00075 #include "G4ParticleChangeForLoss.hh"
00076 
00077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00079 
00080 using namespace std;
00081 
00082 G4MuBremsstrahlungModel::G4MuBremsstrahlungModel(const G4ParticleDefinition* p,
00083                                                  const G4String& nam)
00084   : G4VEmModel(nam),
00085     particle(0),
00086     sqrte(sqrt(exp(1.))),
00087     bh(202.4),
00088     bh1(446.),
00089     btf(183.),
00090     btf1(1429.),
00091     fParticleChange(0),
00092     lowestKinEnergy(1.0*GeV),
00093     minThreshold(1.0*keV)
00094 {
00095   theGamma = G4Gamma::Gamma();
00096   nist = G4NistManager::Instance();
00097 
00098   mass = rmass = cc = coeff = 1.0;
00099 
00100   fDN[0] = 0.0;
00101   for(G4int i=1; i<93; ++i) {
00102     G4double dn = 1.54*nist->GetA27(i);
00103     fDN[i] = dn;
00104     if(1 < i) {
00105       fDN[i] /= std::pow(dn, 1./G4double(i));
00106     }
00107   }
00108 
00109   if(p) { SetParticle(p); }
00110 }
00111 
00112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00113 
00114 G4MuBremsstrahlungModel::~G4MuBremsstrahlungModel()
00115 {
00116   size_t n = partialSumSigma.size();
00117   if(n > 0) {
00118     for(size_t i=0; i<n; i++) {
00119       delete partialSumSigma[i];
00120     }
00121   }
00122 }
00123 
00124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00125 
00126 G4double G4MuBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
00127                                                const G4MaterialCutsCouple*)
00128 {
00129   return minThreshold;
00130 }
00131 
00132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00133 
00134 void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
00135                                          const G4DataVector& cuts)
00136 {
00137   if(p) { SetParticle(p); }
00138 
00139   // partial cross section is computed for fixed energy
00140   G4double fixedEnergy = 0.5*HighEnergyLimit();
00141 
00142   const G4ProductionCutsTable* theCoupleTable=
00143         G4ProductionCutsTable::GetProductionCutsTable();
00144   if(theCoupleTable) {
00145     G4int numOfCouples = theCoupleTable->GetTableSize();
00146 
00147     G4int nn = partialSumSigma.size();
00148     G4int nc = cuts.size();
00149 
00150     // do we need to perform initialisation?
00151     if(nn == numOfCouples) { return; }
00152 
00153     // clear old data    
00154     if(nn > 0) {
00155       for (G4int ii=0; ii<nn; ii++){
00156         G4DataVector* a = partialSumSigma[ii];
00157         if ( a ) { delete a; }
00158       } 
00159       partialSumSigma.clear();
00160     }
00161     // fill new data
00162     if (numOfCouples>0) {
00163       for (G4int i=0; i<numOfCouples; i++) {
00164         G4double cute = DBL_MAX;
00165 
00166         // protection for usage with extrapolator
00167         if(i < nc) { cute = cuts[i]; }
00168 
00169         const G4MaterialCutsCouple* couple = 
00170           theCoupleTable->GetMaterialCutsCouple(i);
00171         const G4Material* material = couple->GetMaterial();
00172         G4DataVector* dv = ComputePartialSumSigma(material,fixedEnergy,cute);
00173         partialSumSigma.push_back(dv);
00174       }
00175     }
00176   }
00177 
00178   // define pointer to G4ParticleChange
00179   if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
00180 }
00181 
00182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00183 
00184 G4double G4MuBremsstrahlungModel::ComputeDEDXPerVolume(
00185                                               const G4Material* material,
00186                                               const G4ParticleDefinition*,
00187                                                     G4double kineticEnergy,
00188                                                     G4double cutEnergy)
00189 {
00190   G4double dedx = 0.0;
00191   if (kineticEnergy <= lowestKinEnergy) return dedx;
00192 
00193   G4double tmax = kineticEnergy;
00194   G4double cut  = std::min(cutEnergy,tmax);
00195   if(cut < minThreshold) cut = minThreshold;
00196 
00197   const G4ElementVector* theElementVector = material->GetElementVector();
00198   const G4double* theAtomicNumDensityVector =
00199     material->GetAtomicNumDensityVector();
00200 
00201   //  loop for elements in the material
00202   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
00203 
00204     G4double loss = 
00205       ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
00206 
00207     dedx += loss*theAtomicNumDensityVector[i];
00208   }
00209   //  G4cout << "BR e= " << kineticEnergy << "  dedx= " << dedx << G4endl;
00210   if(dedx < 0.) dedx = 0.;
00211   return dedx;
00212 }
00213 
00214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00215 
00216 G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
00217                                                    G4double tkin, G4double cut)
00218 {
00219   G4double totalEnergy = mass + tkin;
00220   G4double ak1 = 0.05;
00221   G4int    k2=5;
00222   G4double xgi[]={0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
00223   G4double wgi[]={0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
00224   G4double loss = 0.;
00225 
00226   G4double vcut = cut/totalEnergy;
00227   G4double vmax = tkin/totalEnergy;
00228 
00229   G4double aaa = 0.;
00230   G4double bbb = vcut;
00231   if(vcut>vmax) bbb=vmax ;
00232   G4int kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
00233   if (kkk < 1) { kkk = 1; }
00234   G4double hhh=(bbb-aaa)/float(kkk) ;
00235 
00236   G4double aa = aaa;
00237   for(G4int l=0; l<kkk; l++)
00238   {
00239     for(G4int i=0; i<6; i++)
00240     {
00241       G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
00242       loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
00243     }
00244     aa += hhh;
00245   }
00246 
00247   loss *=hhh*totalEnergy ;
00248 
00249   return loss;
00250 }
00251 
00252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00253 
00254 G4double G4MuBremsstrahlungModel::ComputeMicroscopicCrossSection(
00255                                            G4double tkin,
00256                                            G4double Z,
00257                                            G4double cut)
00258 {
00259   G4double totalEnergy = tkin + mass;
00260   G4double ak1 = 2.3;
00261   G4int    k2  = 4;
00262   G4double xgi[]={0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
00263   G4double wgi[]={0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
00264   G4double cross = 0.;
00265 
00266   if(cut >= tkin) return cross;
00267 
00268   G4double vcut = cut/totalEnergy;
00269   G4double vmax = tkin/totalEnergy;
00270 
00271   G4double aaa = log(vcut);
00272   G4double bbb = log(vmax);
00273   G4int    kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
00274   G4double hhh = (bbb-aaa)/G4double(kkk);
00275 
00276   G4double aa = aaa;
00277 
00278   for(G4int l=0; l<kkk; l++)
00279   {
00280     for(G4int i=0; i<6; i++)
00281     {
00282       G4double ep = exp(aa + xgi[i]*hhh)*totalEnergy;
00283       cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
00284     }
00285     aa += hhh;
00286   }
00287 
00288   cross *=hhh;
00289 
00290   //G4cout << "BR e= " << tkin<< "  cross= " << cross/barn << G4endl;
00291 
00292   return cross;
00293 }
00294 
00295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00296 
00297 G4double G4MuBremsstrahlungModel::ComputeDMicroscopicCrossSection(
00298                                            G4double tkin,
00299                                            G4double Z,
00300                                            G4double gammaEnergy)
00301 //  differential cross section
00302 {
00303   G4double dxsection = 0.;
00304 
00305   if( gammaEnergy > tkin) return dxsection ;
00306 
00307   G4double E = tkin + mass ;
00308   G4double v = gammaEnergy/E ;
00309   G4double delta = 0.5*mass*mass*v/(E-gammaEnergy) ;
00310   G4double rab0=delta*sqrte ;
00311 
00312   G4int iz = G4int(Z);
00313   if(iz < 1) iz = 1;
00314   else if(iz > 92) iz = 92;
00315 
00316   G4double z13 = 1.0/nist->GetZ13(iz);
00317   G4double dnstar = fDN[iz];
00318 
00319   G4double b,b1;
00320 
00321   if(1 == iz) {
00322     b  = bh;
00323     b1 = bh1;
00324   } else {
00325     b  = btf;
00326     b1 = btf1;
00327   }
00328 
00329   // nucleus contribution logarithm
00330   G4double rab1=b*z13;
00331   G4double fn=log(rab1/(dnstar*(electron_mass_c2+rab0*rab1))*
00332               (mass+delta*(dnstar*sqrte-2.))) ;
00333   if(fn <0.) fn = 0. ;
00334   // electron contribution logarithm
00335   G4double epmax1=E/(1.+0.5*mass*rmass/E) ;
00336   G4double fe=0.;
00337   if(gammaEnergy<epmax1)
00338   {
00339     G4double rab2=b1*z13*z13 ;
00340     fe=log(rab2*mass/((1.+delta*rmass/(electron_mass_c2*sqrte))*
00341                               (electron_mass_c2+rab0*rab2))) ;
00342     if(fe<0.) fe=0. ;
00343   }
00344 
00345   dxsection = coeff*(1.-v*(1. - 0.75*v))*Z*(fn*Z + fe)/gammaEnergy;
00346 
00347   return dxsection;
00348 }
00349 
00350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00351 
00352 G4double G4MuBremsstrahlungModel::ComputeCrossSectionPerAtom(
00353                                            const G4ParticleDefinition*,
00354                                                  G4double kineticEnergy,
00355                                                  G4double Z, G4double,
00356                                                  G4double cutEnergy,
00357                                                  G4double maxEnergy)
00358 {
00359   G4double cross = 0.0;
00360   if (kineticEnergy <= lowestKinEnergy) return cross;
00361   G4double tmax = std::min(maxEnergy, kineticEnergy);
00362   G4double cut  = std::min(cutEnergy, kineticEnergy);
00363   if(cut < minThreshold) cut = minThreshold;
00364   if (cut >= tmax) return cross;
00365 
00366   cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
00367   if(tmax < kineticEnergy) {
00368     cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
00369   }
00370   return cross;
00371 }
00372 
00373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00374 
00375 G4DataVector* G4MuBremsstrahlungModel::ComputePartialSumSigma(
00376                                        const G4Material* material,
00377                                        G4double kineticEnergy,
00378                                        G4double cut)
00379 
00380 // Build the table of cross section per element. 
00381 // The table is built for material 
00382 // This table is used to select randomly an element in the material.
00383 {
00384   G4int nElements = material->GetNumberOfElements();
00385   const G4ElementVector* theElementVector = material->GetElementVector();
00386   const G4double* theAtomNumDensityVector = 
00387     material->GetAtomicNumDensityVector();
00388 
00389   G4DataVector* dv = new G4DataVector();
00390 
00391   G4double cross = 0.0;
00392 
00393   for (G4int i=0; i<nElements; i++ ) {
00394     cross += theAtomNumDensityVector[i] 
00395       * ComputeMicroscopicCrossSection(kineticEnergy, 
00396                                        (*theElementVector)[i]->GetZ(), cut);
00397     dv->push_back(cross);
00398   }
00399   return dv;
00400 }
00401 
00402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00403 
00404 void G4MuBremsstrahlungModel::SampleSecondaries(
00405                               std::vector<G4DynamicParticle*>* vdp,
00406                               const G4MaterialCutsCouple* couple,
00407                               const G4DynamicParticle* dp,
00408                               G4double minEnergy,
00409                               G4double maxEnergy)
00410 {
00411   G4double kineticEnergy = dp->GetKineticEnergy();
00412   // check against insufficient energy
00413   G4double tmax = std::min(kineticEnergy, maxEnergy);
00414   G4double tmin = std::min(kineticEnergy, minEnergy);
00415   if(tmin < minThreshold) tmin = minThreshold;
00416   if(tmin >= tmax) return;
00417 
00418   // ===== sampling of energy transfer ======
00419 
00420   G4ParticleMomentum partDirection = dp->GetMomentumDirection();
00421 
00422   // select randomly one element constituing the material
00423   const G4Element* anElement = SelectRandomAtom(couple);
00424   G4double Z = anElement->GetZ();
00425 
00426   G4double totalEnergy   = kineticEnergy + mass;
00427   G4double totalMomentum = sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
00428 
00429   G4double func1 = tmin*
00430     ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
00431 
00432   G4double lnepksi, epksi;
00433   G4double func2;
00434 
00435   G4double xmin = log(tmin/MeV);
00436   G4double xmax = log(kineticEnergy/tmin);
00437 
00438   do {
00439     lnepksi = xmin + G4UniformRand()*xmax;
00440     epksi   = MeV*exp(lnepksi);
00441     func2   = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
00442 
00443   } while(func2 < func1*G4UniformRand());
00444 
00445   G4double gEnergy = epksi;
00446 
00447   // ===== sample angle =====
00448 
00449   G4double gam  = totalEnergy/mass;
00450   G4double rmax = gam*std::min(1.0, totalEnergy/gEnergy - 1.0);
00451   G4double rmax2= rmax*rmax;
00452   G4double x = G4UniformRand()*rmax2/(1.0 + rmax2);
00453 
00454   G4double theta = sqrt(x/(1.0 - x))/gam;
00455   G4double sint  = sin(theta);
00456   G4double phi   = twopi * G4UniformRand() ;
00457   G4double dirx  = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
00458 
00459   G4ThreeVector gDirection(dirx, diry, dirz);
00460   gDirection.rotateUz(partDirection);
00461 
00462   partDirection *= totalMomentum;
00463   partDirection -= gEnergy*gDirection;
00464   partDirection = partDirection.unit();
00465 
00466   // primary change
00467   kineticEnergy -= gEnergy;
00468   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00469   fParticleChange->SetProposedMomentumDirection(partDirection);
00470 
00471   // save secondary
00472   G4DynamicParticle* aGamma = 
00473     new G4DynamicParticle(theGamma,gDirection,gEnergy);
00474   vdp->push_back(aGamma);
00475 }
00476 
00477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00478 
00479 const G4Element* G4MuBremsstrahlungModel::SelectRandomAtom(
00480            const G4MaterialCutsCouple* couple) const
00481 {
00482   // select randomly 1 element within the material
00483 
00484   const G4Material* material = couple->GetMaterial();
00485   G4int nElements = material->GetNumberOfElements();
00486   const G4ElementVector* theElementVector = material->GetElementVector();
00487   if(1 == nElements) { return (*theElementVector)[0]; }
00488   else if(1 > nElements) { return 0; }
00489 
00490   G4DataVector* dv = partialSumSigma[couple->GetIndex()];
00491   G4double rval = G4UniformRand()*((*dv)[nElements-1]);
00492   for (G4int i=0; i<nElements; i++) {
00493     if (rval <= (*dv)[i]) { return (*theElementVector)[i]; }
00494   }
00495   return (*theElementVector)[nElements-1];
00496 }
00497 
00498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Generated on Mon May 27 17:48:53 2013 for Geant4 by  doxygen 1.4.7