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 #include <typeinfo>
00030
00031 #include "globals.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4VScatteringCollision.hh"
00034 #include "G4KineticTrack.hh"
00035 #include "G4VCrossSectionSource.hh"
00036 #include "G4Proton.hh"
00037 #include "G4Neutron.hh"
00038 #include "G4XNNElastic.hh"
00039 #include "G4AngularDistribution.hh"
00040 #include "G4ThreeVector.hh"
00041 #include "G4LorentzVector.hh"
00042 #include "G4LorentzRotation.hh"
00043 #include "G4KineticTrackVector.hh"
00044 #include "Randomize.hh"
00045 #include "G4PionPlus.hh"
00046
00047 G4VScatteringCollision::G4VScatteringCollision()
00048 {
00049 theAngularDistribution = new G4AngularDistribution(true);
00050 }
00051
00052
00053 G4VScatteringCollision::~G4VScatteringCollision()
00054 {
00055 delete theAngularDistribution;
00056 }
00057
00058
00059 G4KineticTrackVector* G4VScatteringCollision::FinalState(const G4KineticTrack& trk1,
00060 const G4KineticTrack& trk2) const
00061 {
00062 const G4VAngularDistribution* angDistribution = GetAngularDistribution();
00063 G4LorentzVector p = trk1.Get4Momentum() + trk2.Get4Momentum();
00064 G4double sqrtS = p.m();
00065 G4double S = sqrtS * sqrtS;
00066
00067 std::vector<const G4ParticleDefinition*> OutputDefinitions = GetOutgoingParticles();
00068 if (OutputDefinitions.size() != 2)
00069 throw G4HadronicException(__FILE__, __LINE__, "G4VScatteringCollision: Too many output particles!");
00070
00071 if (OutputDefinitions[0]->IsShortLived() && OutputDefinitions[1]->IsShortLived())
00072 {
00073 if(getenv("G4KCDEBUG")) G4cerr << "two shortlived for Type = "<<typeid(*this).name()<<G4endl;
00074
00075 }
00076
00077 G4double outm1 = OutputDefinitions[0]->GetPDGMass();
00078 G4double outm2 = OutputDefinitions[1]->GetPDGMass();
00079
00080 if (OutputDefinitions[0]->IsShortLived())
00081 {
00082 outm1 = SampleResonanceMass(outm1,
00083 OutputDefinitions[0]->GetPDGWidth(),
00084 G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass(),
00085 sqrtS-(G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass()));
00086
00087 }
00088 if (OutputDefinitions[1]->IsShortLived())
00089 {
00090 outm2 = SampleResonanceMass(outm2, OutputDefinitions[1]->GetPDGWidth(),
00091 G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass(),
00092 sqrtS-outm1);
00093 }
00094
00095
00096 G4double cosTheta = angDistribution->CosTheta(S, trk1.GetActualMass(), trk2.GetActualMass());
00097 G4double phi = angDistribution->Phi();
00098
00099
00100 G4LorentzRotation fromCMSFrame(p.boostVector());
00101 G4LorentzRotation toCMSFrame(fromCMSFrame.inverse());
00102 G4LorentzVector TempPtr = toCMSFrame*trk1.Get4Momentum();
00103 G4LorentzRotation toZ;
00104 toZ.rotateZ(-1*TempPtr.phi());
00105 toZ.rotateY(-1*TempPtr.theta());
00106 G4LorentzRotation toCMS(toZ.inverse());
00107
00108 G4ThreeVector pFinal1(std::sin(std::acos(cosTheta))*std::cos(phi), std::sin(std::acos(cosTheta))*std::sin(phi), cosTheta);
00109
00110
00111 G4double pCM = std::sqrt( (S-(outm1+outm2)*(outm1+outm2)) * (S-(outm1-outm2)*(outm1-outm2)) /(4.*S));
00112 pFinal1 = pFinal1 * pCM;
00113 G4ThreeVector pFinal2 = -pFinal1;
00114
00115 G4double eFinal1 = std::sqrt(pFinal1.mag2() + outm1*outm1);
00116 G4double eFinal2 = std::sqrt(pFinal2.mag2() + outm2*outm2);
00117
00118 G4LorentzVector p4Final1(pFinal1, eFinal1);
00119 G4LorentzVector p4Final2(pFinal2, eFinal2);
00120 p4Final1 = toCMS*p4Final1;
00121 p4Final2 = toCMS*p4Final2;
00122
00123
00124
00125 G4LorentzRotation toLabFrame(p.boostVector());
00126 p4Final1 *= toLabFrame;
00127 p4Final2 *= toLabFrame;
00128
00129
00130
00131 G4double chargeBalance = OutputDefinitions[0]->GetPDGCharge()+OutputDefinitions[1]->GetPDGCharge();
00132 chargeBalance-= trk1.GetDefinition()->GetPDGCharge();
00133 chargeBalance-= trk2.GetDefinition()->GetPDGCharge();
00134 if(std::abs(chargeBalance) >.1)
00135 {
00136 G4cout << "Charges in "<<typeid(*this).name()<<G4endl;
00137 G4cout << OutputDefinitions[0]->GetPDGCharge()<<" "<<OutputDefinitions[0]->GetParticleName()
00138 << OutputDefinitions[1]->GetPDGCharge()<<" "<<OutputDefinitions[1]->GetParticleName()
00139 << trk1.GetDefinition()->GetPDGCharge()<<" "<<trk1.GetDefinition()->GetParticleName()
00140 << trk2.GetDefinition()->GetPDGCharge()<<" "<<trk2.GetDefinition()->GetParticleName()<<G4endl;
00141 }
00142 G4KineticTrack* final1 = new G4KineticTrack(const_cast<G4ParticleDefinition *>(OutputDefinitions[0]), 0.0,
00143 trk1.GetPosition(), p4Final1);
00144 G4KineticTrack* final2 = new G4KineticTrack(const_cast<G4ParticleDefinition *>(OutputDefinitions[1]), 0.0,
00145 trk2.GetPosition(), p4Final2);
00146
00147 G4KineticTrackVector* finalTracks = new G4KineticTrackVector;
00148
00149 finalTracks->push_back(final1);
00150 finalTracks->push_back(final2);
00151
00152 return finalTracks;
00153 }
00154
00155
00156
00157 double G4VScatteringCollision::SampleResonanceMass(const double poleMass,
00158 const double gamma,
00159 const double aMinMass,
00160 const double maxMass) const
00161 {
00162
00163
00164
00165
00166 G4double minMass = aMinMass;
00167 if (minMass > maxMass) G4cerr << "##################### SampleResonanceMass: particle out of mass range" << G4endl;
00168 if(minMass > maxMass) minMass -= G4PionPlus::PionPlus()->GetPDGMass();
00169 if(minMass > maxMass) minMass = 0;
00170
00171 if (gamma < 1E-10*GeV)
00172 return std::max(minMass,std::min(maxMass, poleMass));
00173 else {
00174 double fmin = BrWigInt0(minMass, gamma, poleMass);
00175 double fmax = BrWigInt0(maxMass, gamma, poleMass);
00176 double f = fmin + (fmax-fmin)*G4UniformRand();
00177 return BrWigInv(f, gamma, poleMass);
00178 }
00179 }