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 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4DecayProducts.hh"
00041 #include "G4VDecayChannel.hh"
00042 #include "G4DalitzDecayChannel.hh"
00043 #include "G4PhaseSpaceDecayChannel.hh"
00044 #include "Randomize.hh"
00045 #include "G4LorentzVector.hh"
00046 #include "G4LorentzRotation.hh"
00047
00048 G4DalitzDecayChannel::G4DalitzDecayChannel()
00049 :G4VDecayChannel()
00050 {
00051 }
00052
00053 G4DalitzDecayChannel::G4DalitzDecayChannel(
00054 const G4String& theParentName,
00055 G4double theBR,
00056 const G4String& theLeptonName,
00057 const G4String& theAntiLeptonName)
00058 :G4VDecayChannel("Dalitz Decay",1)
00059 {
00060
00061 SetParent(theParentName);
00062 SetBR(theBR);
00063 SetNumberOfDaughters(3);
00064 G4String gammaName = "gamma";
00065 SetDaughter(idGamma, gammaName);
00066 SetDaughter(idLepton, theLeptonName);
00067 SetDaughter(idAntiLepton, theAntiLeptonName);
00068 }
00069
00070 G4DalitzDecayChannel::~G4DalitzDecayChannel()
00071 {
00072 }
00073
00074 G4DalitzDecayChannel::G4DalitzDecayChannel(const G4DalitzDecayChannel &right):
00075 G4VDecayChannel(right)
00076 {
00077 }
00078
00079 G4DalitzDecayChannel & G4DalitzDecayChannel::operator=(const G4DalitzDecayChannel & right)
00080 {
00081 if (this != &right) {
00082 kinematics_name = right.kinematics_name;
00083 verboseLevel = right.verboseLevel;
00084 rbranch = right.rbranch;
00085
00086
00087 parent_name = new G4String(*right.parent_name);
00088
00089
00090 ClearDaughtersName();
00091
00092
00093 numberOfDaughters = right.numberOfDaughters;
00094 if ( numberOfDaughters >0 ) {
00095 if (daughters_name !=0) ClearDaughtersName();
00096 daughters_name = new G4String*[numberOfDaughters];
00097
00098 for (G4int index=0; index < numberOfDaughters; index++) {
00099 daughters_name[index] = new G4String(*right.daughters_name[index]);
00100 }
00101 }
00102 }
00103 return *this;
00104 }
00105
00106 G4DecayProducts *G4DalitzDecayChannel::DecayIt(G4double)
00107 {
00108 #ifdef G4VERBOSE
00109 if (GetVerboseLevel()>1) G4cout << "G4DalitzDecayChannel::DecayIt ";
00110 #endif
00111 if (parent == 0) FillParent();
00112 if (daughters == 0) FillDaughters();
00113
00114
00115 G4double parentmass = parent->GetPDGMass();
00116
00117
00118 G4ThreeVector dummy;
00119 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00120
00121
00122 G4double leptonmass = daughters[idLepton]->GetPDGMass();
00123
00124
00125 G4double xmin = 2.0*std::log(2.0*leptonmass);
00126 G4double xmax = 2.0*std::log(parentmass);
00127 G4double wmax = 1.5;
00128 G4double x, w, ww, w1, w2, w3, t;
00129 do {
00130 x = G4UniformRand()*(xmax-xmin) + xmin;
00131 w = G4UniformRand()*wmax;
00132 t = std::exp(x);
00133 w1 = (1.0-4.0*leptonmass*leptonmass/t);
00134 if ( w1 > 0.0) {
00135 w2 = ( 1.0 + 2.0*leptonmass*leptonmass/t );
00136 w3 = ( 1.0 - t/parentmass/parentmass );
00137 w3 = w3 * w3 * w3;
00138 ww = w3 * w2 * std::sqrt(w1);
00139 } else {
00140 ww = 0.0;
00141 }
00142 } while (w > ww);
00143
00144
00145 G4double Pgamma =
00146 G4PhaseSpaceDecayChannel::Pmx(parentmass, 0.0, std::sqrt(t));
00147 G4double costheta = 2.*G4UniformRand()-1.0;
00148 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00149 G4double phi = twopi*G4UniformRand()*rad;
00150 G4ThreeVector gdirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00151
00152
00153 G4DynamicParticle * gammaparticle
00154 = new G4DynamicParticle(daughters[idGamma] , gdirection, Pgamma);
00155
00156
00157 G4double beta = Pgamma/(parentmass-Pgamma);
00158
00159
00160 G4double Plepton =
00161 G4PhaseSpaceDecayChannel::Pmx(std::sqrt(t),leptonmass, leptonmass);
00162 G4double Elepton = std::sqrt(Plepton*Plepton + leptonmass*leptonmass );
00163 costheta = 2.*G4UniformRand()-1.0;
00164 sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00165 phi = twopi*G4UniformRand()*rad;
00166 G4ThreeVector ldirection(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00167
00168 G4DynamicParticle * leptonparticle
00169 = new G4DynamicParticle(daughters[idLepton] ,
00170 ldirection, Elepton-leptonmass );
00171 G4DynamicParticle * antileptonparticle
00172 = new G4DynamicParticle(daughters[idAntiLepton] ,
00173 -1.0*ldirection, Elepton-leptonmass );
00174
00175 G4LorentzVector p4 = leptonparticle->Get4Momentum();
00176 p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
00177 leptonparticle->Set4Momentum(p4);
00178 p4 = antileptonparticle->Get4Momentum();
00179 p4.boost( -1.0*gdirection.x()*beta, -1.0*gdirection.y()*beta, -1.0*gdirection.z()*beta);
00180 antileptonparticle->Set4Momentum(p4);
00181
00182
00183 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00184 delete parentparticle;
00185 products->PushProducts(gammaparticle);
00186 products->PushProducts(leptonparticle);
00187 products->PushProducts(antileptonparticle);
00188
00189 #ifdef G4VERBOSE
00190 if (GetVerboseLevel()>1) {
00191 G4cout << "G4DalitzDecayChannel::DecayIt ";
00192 G4cout << " create decay products in rest frame " <<G4endl;
00193 products->DumpInfo();
00194 }
00195 #endif
00196 return products;
00197 }
00198
00199
00200
00201
00202