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
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 #include "G4eeToPGammaModel.hh"
00050 #include "Randomize.hh"
00051 #include "G4PhysicalConstants.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4PionZero.hh"
00054 #include "G4Eta.hh"
00055 #include "G4Gamma.hh"
00056 #include "G4DynamicParticle.hh"
00057 #include "G4PhysicsVector.hh"
00058 #include "G4PhysicsLinearVector.hh"
00059 #include "G4eeCrossSections.hh"
00060
00061
00062
00063 using namespace std;
00064
00065 G4eeToPGammaModel::G4eeToPGammaModel(G4eeCrossSections* cr, const G4String& npart):
00066 cross(cr)
00067 {
00068 pi0 = G4PionZero::PionZero();
00069 if(npart == "pi0") {
00070 massR = 782.62*MeV;
00071 particle = pi0;
00072 } else {
00073 massR = 1019.46*MeV;
00074 particle = G4Eta::Eta();
00075 }
00076 massP = particle->GetPDGMass();
00077 }
00078
00079
00080
00081 G4eeToPGammaModel::~G4eeToPGammaModel()
00082 {}
00083
00084
00085
00086 G4double G4eeToPGammaModel::ThresholdEnergy() const
00087 {
00088 return LowEnergy();
00089 }
00090
00091
00092
00093 G4double G4eeToPGammaModel::PeakEnergy() const
00094 {
00095 return massR;
00096 }
00097
00098
00099
00100 G4double G4eeToPGammaModel::ComputeCrossSection(G4double e) const
00101 {
00102 G4double ee = std::min(HighEnergy(),e);
00103 G4double xs;
00104 if(particle == pi0) xs = cross->CrossSectionPi0G(ee);
00105 else xs = cross->CrossSectionEtaG(ee);
00106 return xs;
00107 }
00108
00109
00110
00111 G4PhysicsVector* G4eeToPGammaModel::PhysicsVector(G4double emin,
00112 G4double emax) const
00113 {
00114 G4double tmin = std::max(emin, ThresholdEnergy());
00115 G4double tmax = std::max(tmin, emax);
00116 G4int nbins = (G4int)((tmax - tmin)/(5.*MeV));
00117 G4PhysicsVector* v = new G4PhysicsLinearVector(emin,emax,nbins);
00118 v->SetSpline(true);
00119 return v;
00120 }
00121
00122
00123
00124 void G4eeToPGammaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
00125 G4double e, const G4ThreeVector& direction)
00126 {
00127 G4double egam = 0.5*e*(1.0 - massP*massP/(massR*massR));
00128 G4double tkin = e - egam - massP;
00129 if(tkin < 0.0) tkin = 0.0;
00130 G4double cost;
00131 do {
00132 cost = 2.0*G4UniformRand() - 1.0;
00133 } while( 2.0*G4UniformRand() > 1.0 + cost*cost );
00134
00135 G4double sint = sqrt(1.0 - cost*cost);
00136 G4double phi = twopi * G4UniformRand();
00137
00138 G4ThreeVector dir(sint*cos(phi),sint*sin(phi), cost);
00139 dir.rotateUz(direction);
00140
00141
00142 G4DynamicParticle* p1 =
00143 new G4DynamicParticle(particle,dir,tkin);
00144 G4DynamicParticle* p2 =
00145 new G4DynamicParticle(G4Gamma::Gamma(),-dir,egam);
00146 newp->push_back(p1);
00147 newp->push_back(p2);
00148 }
00149
00150
00151