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 #include "G4PionRadiativeDecayChannel.hh"
00041
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "Randomize.hh"
00045 #include "G4DecayProducts.hh"
00046 #include "G4LorentzVector.hh"
00047
00048 G4PionRadiativeDecayChannel::G4PionRadiativeDecayChannel()
00049 : G4VDecayChannel(),
00050 beta(0.),cib(0.),csdp(0.),csdm(0.),cif(0.),cig(0.),
00051 xl(0.), yl(0.), xu(0.), yu(0.), d2wmax(0.)
00052 {
00053 }
00054
00055 G4PionRadiativeDecayChannel::
00056 G4PionRadiativeDecayChannel(const G4String& theParentName,
00057 G4double theBR)
00058 : G4VDecayChannel("Radiative Pion Decay",1)
00059 {
00060
00061 if (theParentName == "pi+") {
00062 SetBR(theBR);
00063 SetParent("pi+");
00064 SetNumberOfDaughters(3);
00065 SetDaughter(0, "e+");
00066 SetDaughter(1, "gamma");
00067 SetDaughter(2, "nu_e");
00068 } else if (theParentName == "pi-") {
00069 SetBR(theBR);
00070 SetParent("pi-");
00071 SetNumberOfDaughters(3);
00072 SetDaughter(0, "e-");
00073 SetDaughter(1, "gamma");
00074 SetDaughter(2, "anti_nu_e");
00075 } else {
00076 #ifdef G4VERBOSE
00077 if (GetVerboseLevel()>0) {
00078 G4cout << "G4RadiativePionDecayChannel:: constructor :";
00079 G4cout << " parent particle is not charged pion but ";
00080 G4cout << theParentName << G4endl;
00081 }
00082 #endif
00083 }
00084
00085 beta = 3.6612e-03;
00086
00087 cib = 1.16141e-03;
00088 csdp = 3.45055e-02;
00089 csdm = 5.14122e-03;
00090 cif = 4.63543e-05;
00091 cig = 1.78928e-05;
00092
00093 xl = 2.*0.1*MeV/139.57*MeV;
00094 yl = ((1.-xl) + std::sqrt((1-xl)*(1-xl)+4*beta*beta))/2.;
00095
00096 xu = 1. - (yl - std::sqrt(yl*yl-4.*beta*beta))/2.;
00097 yu = 1. + beta*beta;
00098
00099 d2wmax = D2W(xl,yl);
00100
00101 }
00102
00103 G4PionRadiativeDecayChannel::~G4PionRadiativeDecayChannel()
00104 {
00105 }
00106 G4PionRadiativeDecayChannel::G4PionRadiativeDecayChannel(const G4PionRadiativeDecayChannel &right)
00107 :G4VDecayChannel(right),
00108 beta(right.beta),cib(right.cib),csdp(right.csdp),
00109 csdm(right.csdm),cif(right.cif),cig(right.cig),
00110 xl(right.xl), yl(right.yl), xu(right.xu), yu(right.yu),
00111 d2wmax(right.d2wmax)
00112 {
00113 }
00114
00115 G4PionRadiativeDecayChannel & G4PionRadiativeDecayChannel::operator=(const G4PionRadiativeDecayChannel & right)
00116 {
00117 if (this != &right) {
00118 kinematics_name = right.kinematics_name;
00119 verboseLevel = right.verboseLevel;
00120 rbranch = right.rbranch;
00121
00122
00123 parent_name = new G4String(*right.parent_name);
00124
00125
00126 ClearDaughtersName();
00127
00128
00129 numberOfDaughters = right.numberOfDaughters;
00130 if ( numberOfDaughters >0 ) {
00131 if (daughters_name !=0) ClearDaughtersName();
00132 daughters_name = new G4String*[numberOfDaughters];
00133
00134 for (G4int index=0; index < numberOfDaughters; index++) {
00135 daughters_name[index] = new G4String(*right.daughters_name[index]);
00136 }
00137 }
00138 beta = right.beta;
00139 cib = right.cib;
00140 csdp = right.csdp;
00141 csdm = right.csdm;
00142 cif = right.cif;
00143 cig = right.cig;
00144 xl = right.xl;
00145 yl = right.yl;
00146 xu = right.xu;
00147 yu = right.yu;
00148 d2wmax = right.d2wmax;
00149 }
00150 return *this;
00151 }
00152
00153 G4DecayProducts *G4PionRadiativeDecayChannel::DecayIt(G4double)
00154 {
00155
00156 #ifdef G4VERBOSE
00157 if (GetVerboseLevel()>1)
00158 G4cout << "G4PionRadiativeDecayChannel::DecayIt ";
00159 #endif
00160
00161 if (parent == 0) FillParent();
00162 if (daughters == 0) FillDaughters();
00163
00164
00165 G4double parentmass = parent->GetPDGMass();
00166
00167 G4double EMPI = parentmass;
00168
00169
00170 G4double daughtermass[3];
00171 G4double sumofdaughtermass = 0.0;
00172 for (G4int index=0; index<3; index++){
00173 daughtermass[index] = daughters[index]->GetPDGMass();
00174 sumofdaughtermass += daughtermass[index];
00175 }
00176
00177 G4double EMASS = daughtermass[0];
00178
00179
00180 G4ThreeVector dummy;
00181 G4DynamicParticle * parentparticle =
00182 new G4DynamicParticle( parent, dummy, 0.0);
00183
00184 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00185 delete parentparticle;
00186
00187 G4double x, y, d2w;
00188
00189 do {
00190
00191 do {
00192
00193 x = xl + G4UniformRand()*(xu-xl);
00194 y = yl + G4UniformRand()*(yu-yl);
00195
00196 } while (x+y <= 1.);
00197
00198 d2w = D2W(x,y);
00199
00200 } while (d2w <= G4UniformRand()*d2wmax);
00201
00202
00203
00204
00205
00206 G4double cthetaGE = (y*(x-2.)+2.*(1.-x+beta*beta)) /
00207 (x*std::sqrt(y*y-4.*beta*beta));
00208
00209
00210
00211
00212 G4double G = x * EMPI/2.;
00213 G4double E = y * EMPI/2.;
00214
00215
00216
00217
00218 if (E < EMASS) E = EMASS;
00219
00220
00221 G4double daughtermomentum[2];
00222
00223 daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS);
00224
00225 G4double cthetaE = 2.*G4UniformRand()-1.;
00226 G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE);
00227
00228 G4double phiE = twopi*G4UniformRand()*rad;
00229 G4double cphiE = std::cos(phiE);
00230 G4double sphiE = std::sin(phiE);
00231
00232
00233
00234 G4double px = sthetaE*cphiE;
00235 G4double py = sthetaE*sphiE;
00236 G4double pz = cthetaE;
00237
00238 G4ThreeVector direction0(px,py,pz);
00239
00240 G4DynamicParticle * daughterparticle0
00241 = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
00242
00243 products->PushProducts(daughterparticle0);
00244
00245 daughtermomentum[1] = G;
00246
00247 G4double sthetaGE = std::sqrt(1.-cthetaGE*cthetaGE);
00248
00249 G4double phiGE = twopi*G4UniformRand()*rad;
00250 G4double cphiGE = std::cos(phiGE);
00251 G4double sphiGE = std::sin(phiGE);
00252
00253
00254
00255 px = sthetaGE*cphiGE;
00256 py = sthetaGE*sphiGE;
00257 pz = cthetaGE;
00258
00259 G4ThreeVector direction1(px,py,pz);
00260
00261 direction1.rotateUz(direction0);
00262
00263 G4DynamicParticle * daughterparticle1
00264 = new G4DynamicParticle( daughters[1], daughtermomentum[1]*direction1);
00265
00266 products->PushProducts(daughterparticle1);
00267
00268
00269 #ifdef G4VERBOSE
00270 if (GetVerboseLevel()>1) {
00271 G4cout << "G4PionRadiativeDecayChannel::DecayIt ";
00272 G4cout << " create decay products in rest frame " <<G4endl;
00273 products->DumpInfo();
00274 }
00275 #endif
00276
00277 return products;
00278
00279 }