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 //
0027 // -------------------------------------------------------------------
0028 //
0029 // 21.03.17 V. Grichine based on G4hBremsstrahlungModel
0030 //
0031 // Class Description:
0032 //
0033 // Implementation of energy loss for LDMPhoton emission by hadrons
0034 //
0035
0036 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0037
0038 #include "G4LDMBremModel.hh"
0039
0040 #include "TestParameters.hh"
0041
0042 #include "G4LDMPhoton.hh"
0043 #include "G4Log.hh"
0044 #include "G4ParticleChangeForLoss.hh"
0045 #include "G4PhysicalConstants.hh"
0046 #include "G4SystemOfUnits.hh"
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 using namespace std;
0051
0052 G4LDMBremModel::G4LDMBremModel(const G4ParticleDefinition* p, const G4String& nam)
0053   : G4MuBremsstrahlungModel(p, nam)
0054 {
0055   fEpsilon = TestParameters::GetPointer()->GetAlphaFactor();
0056   theLDMPhoton = G4LDMPhoton::LDMPhoton();
0057   fLDMPhotonMass = theLDMPhoton->GetPDGMass();
0058   minThreshold = 1.2 * fLDMPhotonMass;
0059 }
0060
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0062
0063 G4LDMBremModel::~G4LDMBremModel() {}
0064
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066
0067 G4double G4LDMBremModel::ComputeDEDXPerVolume(const G4Material*, const G4ParticleDefinition*,
0068                                               G4double, G4double)
0069 {
0070   return 0.0;
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 G4double G4LDMBremModel::ComputeDMicroscopicCrossSection(G4double tkin, G4double Z,
0076                                                          G4double gammaEnergy)
0077 //  differential cross section
0078 {
0079   G4double dxsection = 0.;
0080
0081   if (gammaEnergy > tkin || tkin < minThreshold) return dxsection;
0082   /*
0083   G4cout << "G4LDMBremModel m= " << mass
0084          << "  " << particle->GetParticleName()
0085          << "  Egamma(GeV)= " << gammaEnergy/GeV
0086          << "  Ekin(GeV)= " << tkin/GeV << G4endl;
0087   */
0088   G4double E = tkin + mass;
0089   G4double v = gammaEnergy / E;
0090   G4double delta = 0.5 * mass * mass * v / (E - gammaEnergy);
0091   G4double rab0 = delta * sqrte;
0092
0093   G4int iz = std::max(1, std::min(G4lrint(Z), 99));
0094
0095   G4double z13 = 1.0 / nist->GetZ13(iz);
0096   G4double dn = mass * nist->GetA27(iz) / (70. * MeV);
0097
0098   G4double b = btf;
0099   if (1 == iz) b = bh;
0100
0101   // nucleus contribution logarithm
0102   G4double rab1 = b * z13;
0103   G4double fn =
0104     G4Log(rab1 / (dn * (electron_mass_c2 + rab0 * rab1)) * (mass + delta * (dn * sqrte - 2.)));
0105   if (fn < 0.) fn = 0.;
0106
0107   G4double x = 1.0 - v;
0108
0109   if (particle->GetPDGSpin() != 0) {
0110     x += 0.75 * v * v;
0111   }
0112
0113   dxsection = coeff * x * Z * Z * fn / gammaEnergy;
0114   return dxsection;
0115 }
0116
0117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0118
0119 G4double G4LDMBremModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
0120                                                     G4double kineticEnergy, G4double Z, G4double,
0121                                                     G4double cutEnergy, G4double maxEnergy)
0122 {
0123   G4double cross = 0.0;
0124
0125   if (kineticEnergy <= lowestKinEnergy) return cross;
0126
0127   G4double tmax = std::min(maxEnergy, kineticEnergy);
0128   G4double cut = std::min(cutEnergy, kineticEnergy);
0129
0130   cut = std::max(cut, minThreshold);
0131   if (cut >= tmax) return cross;
0132
0133   cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
0134
0135   if (tmax < kineticEnergy) {
0136     cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
0137   }
0138   cross *= fEpsilon * fEpsilon;
0139
0140   return cross;
0141 }
0142
0143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0144
0145 void G4LDMBremModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
0146                                        const G4MaterialCutsCouple* couple,
0147                                        const G4DynamicParticle* dp, G4double minEnergy,
0148                                        G4double maxEnergy)
0149 {
0150   G4double kineticEnergy = dp->GetKineticEnergy();
0151   // check against insufficient energy
0152   G4double tmax = std::min(kineticEnergy, maxEnergy);
0153   G4double tmin = std::min(kineticEnergy, minEnergy);
0154   tmin = std::max(tmin, minThreshold);
0155   if (tmin >= tmax) return;
0156
0157   // ===== sampling of energy transfer ======
0158
0159   G4ParticleMomentum partDirection = dp->GetMomentumDirection();
0160
0161   // select randomly one element constituing the material
0162   const G4Element* anElement = SelectRandomAtom(couple, particle, kineticEnergy);
0163   G4double Z = anElement->GetZ();
0164
0165   G4double totalEnergy = kineticEnergy + mass;
0166   G4double totalMomentum = sqrt(kineticEnergy * (kineticEnergy + 2.0 * mass));
0167
0168   G4double func1 = tmin * ComputeDMicroscopicCrossSection(kineticEnergy, Z, tmin);
0169
0170   G4double lnepksi, epksi;
0171   G4double func2;
0172
0173   G4double xmin = G4Log(tmin / MeV);
0174   G4double xmax = G4Log(tmax / tmin);
0175
0176   do {
0177     lnepksi = xmin + G4UniformRand() * xmax;
0178     epksi = MeV * G4Exp(lnepksi);
0179     func2 = epksi * ComputeDMicroscopicCrossSection(kineticEnergy, Z, epksi);
0180
0181     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
0182   } while (func2 < func1 * G4UniformRand());
0183
0184   G4double gEnergy = std::max(epksi, fLDMPhotonMass);
0185   G4double gMomentum = std::sqrt((epksi - fLDMPhotonMass) * (epksi + fLDMPhotonMass));
0186
0187   // ===== sample angle =====
0188
0189   G4double gam = totalEnergy / mass;
0190   G4double rmax = gam * std::min(1.0, totalEnergy / gEnergy - 1.0);
0191   G4double rmax2 = rmax * rmax;
0192   G4double x = G4UniformRand() * rmax2 / (1.0 + rmax2);
0193
0194   G4double theta = std::sqrt(x / (1.0 - x)) / gam;
0195   G4double sint = std::sin(theta);
0196   G4double phi = twopi * G4UniformRand();
0197   G4double dirx = sint * cos(phi), diry = sint * sin(phi), dirz = cos(theta);
0198
0199   G4ThreeVector gDirection(dirx, diry, dirz);
0200   gDirection.rotateUz(partDirection);
0201
0202   partDirection *= totalMomentum;
0203   partDirection -= gMomentum * gDirection;
0204   partDirection = partDirection.unit();
0205
0206   // primary change
0207
0208   kineticEnergy -= gEnergy;
0209
0210   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
0211   fParticleChange->SetProposedMomentumDirection(partDirection);
0212
0213   // save secondary
0214   G4DynamicParticle* aLDMPhoton = new G4DynamicParticle(theLDMPhoton, gDirection, gEnergy);
0215   vdp->push_back(aLDMPhoton);
0216 }
0217
0218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......