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 "G4ParticleDefinition.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4DecayProducts.hh"
00044 #include "G4VDecayChannel.hh"
00045 #include "G4NeutronBetaDecayChannel.hh"
00046 #include "Randomize.hh"
00047 #include "G4RotationMatrix.hh"
00048 #include "G4LorentzVector.hh"
00049 #include "G4LorentzRotation.hh"
00050
00051 G4NeutronBetaDecayChannel::G4NeutronBetaDecayChannel()
00052 :G4VDecayChannel(),
00053 aENuCorr(-0.102)
00054 {
00055 }
00056
00057 G4NeutronBetaDecayChannel::G4NeutronBetaDecayChannel(
00058 const G4String& theParentName,
00059 G4double theBR)
00060 :G4VDecayChannel("Neutron Decay"),
00061 aENuCorr(-0.102)
00062 {
00063
00064 if (theParentName == "neutron") {
00065 SetBR(theBR);
00066 SetParent("neutron");
00067 SetNumberOfDaughters(3);
00068 SetDaughter(0, "e-");
00069 SetDaughter(1, "anti_nu_e");
00070 SetDaughter(2, "proton");
00071 } else if (theParentName == "anti_neutron") {
00072 SetBR(theBR);
00073 SetParent("anti_neutron");
00074 SetNumberOfDaughters(3);
00075 SetDaughter(0, "e+");
00076 SetDaughter(1, "nu_e");
00077 SetDaughter(2, "anti_proton");
00078 } else {
00079 #ifdef G4VERBOSE
00080 if (GetVerboseLevel()>0) {
00081 G4cout << "G4NeutronBetaDecayChannel:: constructor :";
00082 G4cout << " parent particle is not neutron but ";
00083 G4cout << theParentName << G4endl;
00084 }
00085 #endif
00086 }
00087 }
00088
00089 G4NeutronBetaDecayChannel::~G4NeutronBetaDecayChannel()
00090 {
00091 }
00092
00093 G4NeutronBetaDecayChannel::G4NeutronBetaDecayChannel(const G4NeutronBetaDecayChannel &right)
00094 : G4VDecayChannel(right),
00095 aENuCorr(-0.102)
00096 {
00097 }
00098
00099
00100 G4NeutronBetaDecayChannel & G4NeutronBetaDecayChannel::operator=(const G4NeutronBetaDecayChannel & right)
00101 {
00102 if (this != &right) {
00103 kinematics_name = right.kinematics_name;
00104 verboseLevel = right.verboseLevel;
00105 rbranch = right.rbranch;
00106
00107
00108 parent_name = new G4String(*right.parent_name);
00109
00110
00111 ClearDaughtersName();
00112
00113
00114 numberOfDaughters = right.numberOfDaughters;
00115 if ( numberOfDaughters >0 ) {
00116 if (daughters_name !=0) ClearDaughtersName();
00117 daughters_name = new G4String*[numberOfDaughters];
00118
00119 for (G4int index=0; index < numberOfDaughters; index++) {
00120 daughters_name[index] = new G4String(*right.daughters_name[index]);
00121 }
00122 }
00123 }
00124 return *this;
00125 }
00126
00127 G4DecayProducts *G4NeutronBetaDecayChannel::DecayIt(G4double)
00128 {
00129
00130
00131
00132
00133 #ifdef G4VERBOSE
00134 if (GetVerboseLevel()>1) G4cout << "G4NeutronBetaDecayChannel::DecayIt ";
00135 #endif
00136
00137 if (parent == 0) FillParent();
00138 if (daughters == 0) FillDaughters();
00139
00140
00141 G4double parentmass = parent->GetPDGMass();
00142
00143
00144 G4double daughtermass[3];
00145 G4double sumofdaughtermass = 0.0;
00146 for (G4int index=0; index<3; index++){
00147 daughtermass[index] = daughters[index]->GetPDGMass();
00148 sumofdaughtermass += daughtermass[index];
00149 }
00150 G4double xmax = parentmass-sumofdaughtermass;
00151
00152
00153 G4ThreeVector dummy;
00154 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00155
00156
00157 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00158 delete parentparticle;
00159
00160
00161 G4double daughtermomentum[3];
00162
00163
00164 G4double x;
00165 G4double p;
00166 G4double dm = daughtermass[0];
00167 G4double w;
00168 G4double r;
00169 G4double r0;
00170 do {
00171 x = xmax*G4UniformRand();
00172 p = std::sqrt(x*(x+2.0*dm));
00173 w = 1.0-2.0*G4UniformRand();
00174 r = p*(x+dm)*(xmax-x)*(xmax-x)*(1.0+aENuCorr*p/(x+dm)*w);
00175 r0 = G4UniformRand()*(xmax+dm)*(xmax+dm)*xmax*xmax*(1.0+aENuCorr);
00176 } while (r < r0);
00177
00178
00179
00180
00181 G4double costheta = 2.*G4UniformRand()-1.0;
00182 G4double theta = std::acos(costheta)*rad;
00183 G4double phi = twopi*G4UniformRand()*rad;
00184 G4RotationMatrix rm;
00185 rm.rotateY(theta);
00186 rm.rotateZ(phi);
00187
00188
00189 daughtermomentum[0] = p;
00190 G4ThreeVector direction0(0.0, 0.0, 1.0);
00191 direction0 = rm * direction0;
00192 G4DynamicParticle * daughterparticle0
00193 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
00194 products->PushProducts(daughterparticle0);
00195
00196
00197 G4double eNu;
00198 eNu = (parentmass-daughtermass[2])*(parentmass+daughtermass[2])+(dm*dm)-2.*parentmass*(x+dm);
00199 eNu /= 2.*(parentmass+p*w-(x+dm));
00200 G4double cosn = w;
00201 G4double sinn = std::sqrt((1.0-cosn)*(1.0+cosn));
00202
00203 G4ThreeVector direction1(sinn, 0.0, cosn);
00204 direction1 = rm * direction1;
00205 G4DynamicParticle * daughterparticle1
00206 = new G4DynamicParticle( daughters[1], direction1*eNu);
00207 products->PushProducts(daughterparticle1);
00208
00209
00210 G4double eP;
00211 eP = parentmass-eNu-(x+dm)-daughtermass[2];
00212 G4double pPx = -eNu*sinn;
00213 G4double pPz = -p-eNu*cosn;
00214 G4double pP = std::sqrt(eP*(eP+2.*daughtermass[2]));
00215 G4ThreeVector direction2(pPx/pP, 0.0, pPz/pP);
00216 G4DynamicParticle * daughterparticle2
00217 = new G4DynamicParticle( daughters[2], direction2);
00218 products->PushProducts(daughterparticle2);
00219
00220
00221
00222 #ifdef G4VERBOSE
00223 if (GetVerboseLevel()>1) {
00224 G4cout << "G4NeutronBetaDecayChannel::DecayIt ";
00225 G4cout << " create decay products in rest frame " <<G4endl;
00226 products->DumpInfo();
00227 }
00228 #endif
00229 return products;
00230 }
00231
00232
00233
00234
00235
00236