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 #include <numeric>
00040
00041 #include "G4FermiPhaseSpaceDecay.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4HadronicException.hh"
00044
00045 G4FermiPhaseSpaceDecay::G4FermiPhaseSpaceDecay()
00046 {
00047 g4pow = G4Pow::GetInstance();
00048 }
00049
00050 G4FermiPhaseSpaceDecay::~G4FermiPhaseSpaceDecay()
00051 {
00052 }
00053
00054 std::vector<G4LorentzVector*> *
00055 G4FermiPhaseSpaceDecay::KopylovNBodyDecay(const G4double M,
00056 const std::vector<G4double>& mr) const
00057
00058 {
00059 size_t N = mr.size();
00060
00061 std::vector<G4LorentzVector*>* P = new std::vector<G4LorentzVector*>(N, 0);
00062
00063 G4double mtot = std::accumulate( mr.begin(), mr.end(), 0.0);
00064 G4double mu = mtot;
00065 G4double PFragMagCM = 0.0;
00066 G4double Mass = M;
00067 G4double T = Mass-mtot;
00068 G4LorentzVector PFragCM(0.0,0.0,0.0,0.0);
00069 G4LorentzVector PRestCM(0.0,0.0,0.0,0.0);
00070 G4LorentzVector PRestLab(0.0,0.0,0.0,Mass);
00071
00072 for (size_t k = N-1; k>0; --k)
00073 {
00074 mu -= mr[k];
00075 if (k>1) { T *= BetaKopylov(k); }
00076 else { T = 0.0; }
00077
00078 G4double RestMass = mu + T;
00079
00080 PFragMagCM = PtwoBody(Mass,mr[k],RestMass);
00081
00082
00083 G4ThreeVector RandVector(IsotropicVector(PFragMagCM));
00084
00085 PFragCM.setVect(RandVector);
00086 PFragCM.setE(std::sqrt(PFragMagCM*PFragMagCM + mr[k]*mr[k]));
00087
00088 PRestCM.setVect(-RandVector);
00089 PRestCM.setE(std::sqrt(PFragMagCM*PFragMagCM + RestMass*RestMass));
00090
00091
00092 G4ThreeVector BoostV = PRestLab.boostVector();
00093
00094 PFragCM.boost(BoostV);
00095 PRestCM.boost(BoostV);
00096 PRestLab = PRestCM;
00097
00098 (*P)[k] = new G4LorentzVector(PFragCM);
00099
00100 Mass = RestMass;
00101 }
00102
00103 (*P)[0] = new G4LorentzVector(PRestLab);
00104
00105 return P;
00106 }
00107
00108
00109 std::vector<G4LorentzVector*> *
00110 G4FermiPhaseSpaceDecay::NBodyDecay(G4double M, const std::vector<G4double>& mr) const
00111 {
00112
00113 size_t N = mr.size();
00114 size_t i, j;
00115
00116 G4double mtot = std::accumulate( mr.begin(), mr.end(), 0.0);
00117 G4double Emax = M - mtot + mr[0];
00118 G4double Emin = 0.0;
00119 G4double Wmax = 1.0;
00120 for (i = 1; i < N; i++)
00121 {
00122 Emax += mr[i];
00123 Emin += mr[i-1];
00124 Wmax *= PtwoBody(Emax, Emin, mr[i]);
00125 }
00126
00127 G4int ntries = 0;
00128 G4double weight = 1.0;
00129 std::vector<G4double> p(N, 0.0);
00130 std::vector<G4double> r(N,0.0);
00131 std::vector<G4double> vm(N, 0.0);
00132 r[N-1] = 1.0;
00133
00134 do
00135 {
00136
00137 for (i = 1; i < N-1; i++) { r[i] = G4UniformRand(); }
00138 std::sort(r.begin(),r.end(), std::less<G4double>());
00139
00140
00141 std::partial_sum(mr.begin(), mr.end(), vm.begin());
00142 std::transform(r.begin(), r.end(), r.begin(),
00143 std::bind2nd(std::multiplies<G4double>(), M-mtot));
00144 std::transform(r.begin(), r.end(), vm.begin(), vm.begin(), std::plus<G4double>());
00145
00146
00147 weight = 1.0;
00148 for (j = 0; j < N-1; j++)
00149 {
00150 p[j] = PtwoBody(vm[j+1],vm[j],mr[j+1]);
00151 weight *= p[j];
00152 }
00153 p[N-1] = PtwoBody(vm[N-2],mr[N-2],mr[N-1]);
00154
00155
00156 if (ntries++ > 1000000)
00157 {
00158 throw G4HadronicException(__FILE__, __LINE__, "Failed to decay");
00159 }
00160 }
00161 while ( weight < G4UniformRand()*Wmax );
00162
00163 std::vector<G4LorentzVector*> * P = new std::vector<G4LorentzVector*>(N, 0);
00164
00165 G4ThreeVector a3P = IsotropicVector(p[0]);
00166
00167 (*P)[0] = new G4LorentzVector( a3P, std::sqrt(a3P.mag2()+mr[0]*mr[0]) );
00168 (*P)[1] = new G4LorentzVector(-a3P, std::sqrt(a3P.mag2()+mr[1]*mr[1]) );
00169 for (i = 2; i < N; i++)
00170 {
00171 a3P = IsotropicVector(p[i-1]);
00172 (*P)[i] = new G4LorentzVector(a3P, std::sqrt(a3P.mag2() + mr[i]*mr[i]));
00173 G4ThreeVector Beta = -((*P)[i]->boostVector());
00174
00175 for (j = 0; j < i; j++)
00176 {
00177 (*P)[j]->boost(Beta);
00178 }
00179 }
00180
00181 return P;
00182 }
00183
00184 std::vector<G4LorentzVector*> *
00185 G4FermiPhaseSpaceDecay::TwoBodyDecay(G4double M,
00186 const std::vector<G4double>& mass) const
00187 {
00188 G4double m0 = mass.front();
00189 G4double m1 = mass.back();
00190 G4double mom = PtwoBody(M,m0,m1);
00191 G4ThreeVector p = IsotropicVector(mom);
00192
00193 G4LorentzVector * P41 = new G4LorentzVector;
00194 P41->setVect(p);
00195 P41->setE(std::sqrt(mom*mom + m0*m0));
00196
00197 G4LorentzVector * P42 = new G4LorentzVector;
00198 P42->setVect(-p);
00199 P42->setE(std::sqrt(mom*mom + m1*m1));
00200
00201 std::vector<G4LorentzVector*> * result = new std::vector<G4LorentzVector*>;
00202 result->push_back(P41);
00203 result->push_back(P42);
00204 return result;
00205 }
00206
00207 void
00208 G4FermiPhaseSpaceDecay::DumpProblem(G4double E, G4double P1, G4double P2,
00209 G4double P) const
00210 {
00211 G4cout << "G4FermiPhaseSpaceDecay: problem of decay of M(GeV)= " << E/GeV
00212 << " on M1(GeV)= " << P1/GeV << " and M2(GeV)= " << P2/GeV
00213 << " P(MeV)= " << P/MeV << " < 0" << G4endl;
00214
00215 if(P < -CLHEP::eV) {
00216 throw G4HadronicException(__FILE__, __LINE__,"Error in decay kinematics");
00217 }
00218 }
00219
00220