dmparticle/src/G4LDMBremModel.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file G4LDMBremModel.cc
0027 /// \brief Implementation of the G4LDMBremModel class
0028
0029 // -------------------------------------------------------------------
0030 //
0031 // 21.03.17 V. Grichine based on G4hBremsstrahlungModel
0032 //
0033 // Class Description:
0034 //
0035 // Implementation of energy loss for LDMPhoton emission by hadrons
0036 //
0037
0038 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0039
0040 #include "G4LDMBremModel.hh"
0041
0042 #include "TestParameters.hh"
0043
0044 #include "G4LDMPhoton.hh"
0045 #include "G4Log.hh"
0046 #include "G4ParticleChangeForLoss.hh"
0047 #include "G4PhysicalConstants.hh"
0048 #include "G4SystemOfUnits.hh"
0049
0050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0051
0052 using namespace std;
0053
0054 G4LDMBremModel::G4LDMBremModel(const G4ParticleDefinition* p, const G4String& nam)
0055   : G4MuBremsstrahlungModel(p, nam)
0056 {
0057   fEpsilon = TestParameters::GetPointer()->GetAlphaFactor();
0058   theLDMPhoton = G4LDMPhoton::LDMPhoton();
0059   fLDMPhotonMass = theLDMPhoton->GetPDGMass();
0060   minThreshold = 1.2 * fLDMPhotonMass;
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 G4LDMBremModel::~G4LDMBremModel() {}
0066
0067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0068
0069 G4double G4LDMBremModel::ComputeDEDXPerVolume(const G4Material*, const G4ParticleDefinition*,
0070                                               G4double, G4double)
0071 {
0072   return 0.0;
0073 }
0074
0075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0076
0077 G4double G4LDMBremModel::ComputeDMicroscopicCrossSection(G4double tkin, G4double Z,
0078                                                          G4double gammaEnergy)
0079 //  differential cross section
0080 {
0081   G4double dxsection = 0.;
0082
0083   if (gammaEnergy > tkin || tkin < minThreshold) return dxsection;
0084   /*
0085   G4cout << "G4LDMBremModel m= " << mass
0086          << "  " << particle->GetParticleName()
0087          << "  Egamma(GeV)= " << gammaEnergy/GeV
0088          << "  Ekin(GeV)= " << tkin/GeV << G4endl;
0089   */
0090   G4double E = tkin + mass;
0091   G4double v = gammaEnergy / E;
0092   G4double delta = 0.5 * mass * mass * v / (E - gammaEnergy);
0093   G4double rab0 = delta * sqrte;
0094
0095   G4int iz = std::max(1, std::min(G4lrint(Z), 99));
0096
0097   G4double z13 = 1.0 / nist->GetZ13(iz);
0098   G4double dn = mass * nist->GetA27(iz) / (70. * MeV);
0099
0100   G4double b = btf;
0101   if (1 == iz) b = bh;
0102
0103   // nucleus contribution logarithm
0104   G4double rab1 = b * z13;
0105   G4double fn =
0106     G4Log(rab1 / (dn * (electron_mass_c2 + rab0 * rab1)) * (mass + delta * (dn * sqrte - 2.)));
0107   if (fn < 0.) fn = 0.;
0108
0109   G4double x = 1.0 - v;
0110
0111   if (particle->GetPDGSpin() != 0) {
0112     x += 0.75 * v * v;
0113   }
0114
0115   dxsection = coeff * x * Z * Z * fn / gammaEnergy;
0116   return dxsection;
0117 }
0118
0119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0120
0121 G4double G4LDMBremModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
0122                                                     G4double kineticEnergy, G4double Z, G4double,
0123                                                     G4double cutEnergy, G4double maxEnergy)
0124 {
0125   G4double cross = 0.0;
0126
0127   if (kineticEnergy <= lowestKinEnergy) return cross;
0128
0129   G4double tmax = std::min(maxEnergy, kineticEnergy);
0130   G4double cut = std::min(cutEnergy, kineticEnergy);
0131
0132   cut = std::max(cut, minThreshold);
0133   if (cut >= tmax) return cross;
0134
0135   cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
0136
0137   if (tmax < kineticEnergy) {
0138     cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
0139   }
0140   cross *= fEpsilon * fEpsilon;
0141
0142   return cross;
0143 }
0144
0145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0146
0147 void G4LDMBremModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
0148                                        const G4MaterialCutsCouple* couple,
0149                                        const G4DynamicParticle* dp, G4double minEnergy,
0150                                        G4double maxEnergy)
0151 {
0152   G4double kineticEnergy = dp->GetKineticEnergy();
0153   // check against insufficient energy
0154   G4double tmax = std::min(kineticEnergy, maxEnergy);
0155   G4double tmin = std::min(kineticEnergy, minEnergy);
0156   tmin = std::max(tmin, minThreshold);
0157   if (tmin >= tmax) return;
0158
0159   // ===== sampling of energy transfer ======
0160
0161   G4ParticleMomentum partDirection = dp->GetMomentumDirection();
0162
0163   // select randomly one element constituing the material
0164   const G4Element* anElement = SelectRandomAtom(couple, particle, kineticEnergy);
0165   G4double Z = anElement->GetZ();
0166
0167   G4double totalEnergy = kineticEnergy + mass;
0168   G4double totalMomentum = sqrt(kineticEnergy * (kineticEnergy + 2.0 * mass));
0169
0170   G4double func1 = tmin * ComputeDMicroscopicCrossSection(kineticEnergy, Z, tmin);
0171
0172   G4double lnepksi, epksi;
0173   G4double func2;
0174
0175   G4double xmin = G4Log(tmin / MeV);
0176   G4double xmax = G4Log(tmax / tmin);
0177
0178   do {
0179     lnepksi = xmin + G4UniformRand() * xmax;
0180     epksi = MeV * G4Exp(lnepksi);
0181     func2 = epksi * ComputeDMicroscopicCrossSection(kineticEnergy, Z, epksi);
0182
0183     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
0184   } while (func2 < func1 * G4UniformRand());
0185
0186   G4double gEnergy = std::max(epksi, fLDMPhotonMass);
0187   G4double gMomentum = std::sqrt((epksi - fLDMPhotonMass) * (epksi + fLDMPhotonMass));
0188
0189   // ===== sample angle =====
0190
0191   G4double gam = totalEnergy / mass;
0192   G4double rmax = gam * std::min(1.0, totalEnergy / gEnergy - 1.0);
0193   G4double rmax2 = rmax * rmax;
0194   G4double x = G4UniformRand() * rmax2 / (1.0 + rmax2);
0195
0196   G4double theta = std::sqrt(x / (1.0 - x)) / gam;
0197   G4double sint = std::sin(theta);
0198   G4double phi = twopi * G4UniformRand();
0199   G4double dirx = sint * cos(phi), diry = sint * sin(phi), dirz = cos(theta);
0200
0201   G4ThreeVector gDirection(dirx, diry, dirz);
0202   gDirection.rotateUz(partDirection);
0203
0204   partDirection *= totalMomentum;
0205   partDirection -= gMomentum * gDirection;
0206   partDirection = partDirection.unit();
0207
0208   // primary change
0209
0210   kineticEnergy -= gEnergy;
0211
0212   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
0213   fParticleChange->SetProposedMomentumDirection(partDirection);
0214
0215   // save secondary
0216   G4DynamicParticle* aLDMPhoton = new G4DynamicParticle(theLDMPhoton, gDirection, gEnergy);
0217   vdp->push_back(aLDMPhoton);
0218 }
0219
0220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......