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 "G4eeTo3PiModel.hh"
00050 #include "Randomize.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4PionPlus.hh"
00053 #include "G4PionMinus.hh"
00054 #include "G4PionZero.hh"
00055 #include "G4DynamicParticle.hh"
00056 #include "G4PhysicsVector.hh"
00057 #include "G4PhysicsLinearVector.hh"
00058 #include "G4eeCrossSections.hh"
00059 #include "G4RandomDirection.hh"
00060 #include <complex>
00061
00062
00063
00064 using namespace std;
00065
00066 G4eeTo3PiModel::G4eeTo3PiModel(G4eeCrossSections* cr):
00067 cross(cr)
00068 {
00069 massPi = G4PionPlus::PionPlus()->GetPDGMass();
00070 massPi0 = G4PionZero::PionZero()->GetPDGMass();
00071 massOm = 782.62*MeV;
00072 massPhi = 1019.46*MeV;
00073 gcash = 0.0;
00074 gmax = 1.0;
00075 }
00076
00077
00078
00079 G4eeTo3PiModel::~G4eeTo3PiModel()
00080 {}
00081
00082
00083
00084 G4double G4eeTo3PiModel::ThresholdEnergy() const
00085 {
00086 return std::max(LowEnergy(),2.0*massPi + massPi0);
00087 }
00088
00089
00090
00091 G4double G4eeTo3PiModel::PeakEnergy() const
00092 {
00093 G4double e = massOm;
00094 if(HighEnergy() > massPhi) e = massPhi;
00095 return e;
00096 }
00097
00098
00099
00100 G4double G4eeTo3PiModel::ComputeCrossSection(G4double e) const
00101 {
00102 G4double ee = std::min(HighEnergy(),e);
00103 return cross->CrossSection3pi(ee);
00104 }
00105
00106
00107
00108 G4PhysicsVector* G4eeTo3PiModel::PhysicsVector(G4double emin,
00109 G4double emax) const
00110 {
00111 G4double tmin = std::max(emin, ThresholdEnergy());
00112 G4double tmax = std::max(tmin, emax);
00113 G4int nbins = (G4int)((tmax - tmin)/(1.*MeV));
00114 G4PhysicsVector* v = new G4PhysicsLinearVector(emin,emax,nbins);
00115 v->SetSpline(true);
00116 return v;
00117 }
00118
00119
00120
00121 void G4eeTo3PiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
00122 G4double e, const G4ThreeVector& direction)
00123 {
00124 if(e < ThresholdEnergy()) return;
00125
00126 G4double x0 = massPi0/e;
00127 G4double x1 = massPi/e;
00128
00129 G4LorentzVector w0, w1, w2;
00130 G4ThreeVector dir0, dir1;
00131 G4double e0, p0, e2, p, gg, m01, m02, m12;
00132
00133
00134 G4double edel = 0.5*e*(1.0 + x0*x0 - 4.0*x1*x1) - massPi0;
00135
00136 do {
00137
00138 e0 = edel*G4UniformRand() + massPi0;
00139 p0 = sqrt(e0 - massPi0*massPi0);
00140 dir0 = G4RandomDirection();
00141 w0 = G4LorentzVector(p0*dir0.x(),p0*dir0.y(),p0*dir0.z(),e0);
00142
00143
00144 w1 = G4LorentzVector(-p0*dir0.x(),-p0*dir0.y(),-p0*dir0.z(),e-e0);
00145 G4ThreeVector bst = w1.boostVector();
00146 e2 = 0.25*w1.m2();
00147
00148
00149 p = sqrt(e2 - massPi*massPi);
00150 dir1 = G4RandomDirection();
00151 w2 = G4LorentzVector(p*dir1.x(),p*dir1.y(),p*dir1.z(),sqrt(e2));
00152 w2.boost(bst);
00153 G4double px2 = w2.x();
00154 G4double py2 = w2.y();
00155 G4double pz2 = w2.z();
00156
00157
00158 w1 -= w2;
00159 G4double px1 = w1.x();
00160 G4double py1 = w1.y();
00161 G4double pz1 = w1.z();
00162
00163 m01 = w0*w1;
00164 m02 = w0*w2;
00165 m12 = w1*w2;
00166
00167 G4double px = py1*pz2 - py2*pz1;
00168 G4double py = pz1*px2 - pz2*px1;
00169 G4double pz = px1*py2 - px2*py1;
00170
00171 gg = (px*px + py*py + pz*pz)*
00172 norm( 1.0/cross->DpRho(m01) + 1.0/cross->DpRho(m02)
00173 + 1.0/cross->DpRho(m12) );
00174 if(gg > gmax) {
00175 G4cout << "G4eeTo3PiModel::SampleSecondaries WARNING matrix element g= "
00176 << gg << " > " << gmax << " (majoranta)" << G4endl;
00177 }
00178 if(gg > gcash) gcash = gg;
00179
00180 } while( gmax*G4UniformRand() > gg );
00181
00182 w0.rotateUz(direction);
00183 w1.rotateUz(direction);
00184 w2.rotateUz(direction);
00185
00186
00187 G4DynamicParticle* dp0 =
00188 new G4DynamicParticle(G4PionZero::PionZero(), w0);
00189 G4DynamicParticle* dp1 =
00190 new G4DynamicParticle(G4PionPlus::PionPlus(), w1);
00191 G4DynamicParticle* dp2 =
00192 new G4DynamicParticle(G4PionMinus::PionMinus(), w2);
00193 newp->push_back(dp0);
00194 newp->push_back(dp1);
00195 newp->push_back(dp2);
00196 }
00197
00198
00199