Geant4  10.00.p01
G4VScatteringCollision.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // @hpw@ misses the sampling of two breit wigner in a corelated fashion,
27 // @hpw@ to be usefull for resonance resonance scattering.
28 
29 #include <typeinfo>
30 
31 #include "globals.hh"
32 #include "G4SystemOfUnits.hh"
34 #include "G4KineticTrack.hh"
35 #include "G4VCrossSectionSource.hh"
36 #include "G4Proton.hh"
37 #include "G4Neutron.hh"
38 #include "G4XNNElastic.hh"
39 #include "G4AngularDistribution.hh"
40 #include "G4ThreeVector.hh"
41 #include "G4LorentzVector.hh"
42 #include "G4LorentzRotation.hh"
43 #include "G4KineticTrackVector.hh"
44 #include "Randomize.hh"
45 #include "G4PionPlus.hh"
46 
48 {
50 }
51 
52 
54 {
56 }
57 
58 
60  const G4KineticTrack& trk2) const
61 {
62  const G4VAngularDistribution* angDistribution = GetAngularDistribution();
63  G4LorentzVector p = trk1.Get4Momentum() + trk2.Get4Momentum();
64  G4double sqrtS = p.m();
65  G4double S = sqrtS * sqrtS;
66 
67  std::vector<const G4ParticleDefinition*> OutputDefinitions = GetOutgoingParticles();
68  if (OutputDefinitions.size() != 2)
69  throw G4HadronicException(__FILE__, __LINE__, "G4VScatteringCollision: Too many output particles!");
70 
71  if (OutputDefinitions[0]->IsShortLived() && OutputDefinitions[1]->IsShortLived())
72  {
73  if(getenv("G4KCDEBUG")) G4cerr << "two shortlived for Type = "<<typeid(*this).name()<<G4endl;
74  // throw G4HadronicException(__FILE__, __LINE__, "G4VScatteringCollision: can't handle two shortlived particles!"); // @hpw@
75  }
76 
77  G4double outm1 = OutputDefinitions[0]->GetPDGMass();
78  G4double outm2 = OutputDefinitions[1]->GetPDGMass();
79 
80  if (OutputDefinitions[0]->IsShortLived())
81  {
82  outm1 = SampleResonanceMass(outm1,
83  OutputDefinitions[0]->GetPDGWidth(),
84  G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass(),
85  sqrtS-(G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass()));
86 
87  }
88  if (OutputDefinitions[1]->IsShortLived())
89  {
90  outm2 = SampleResonanceMass(outm2, OutputDefinitions[1]->GetPDGWidth(),
91  G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass(),
92  sqrtS-outm1);
93  }
94 
95  // Angles of outgoing particles
96  G4double cosTheta = angDistribution->CosTheta(S, trk1.GetActualMass(), trk2.GetActualMass());
97  G4double phi = angDistribution->Phi();
98 
99  // Unit vector of three-momentum
100  G4LorentzRotation fromCMSFrame(p.boostVector());
101  G4LorentzRotation toCMSFrame(fromCMSFrame.inverse());
102  G4LorentzVector TempPtr = toCMSFrame*trk1.Get4Momentum();
103  G4LorentzRotation toZ;
104  toZ.rotateZ(-1*TempPtr.phi());
105  toZ.rotateY(-1*TempPtr.theta());
106  G4LorentzRotation toCMS(toZ.inverse());
107 
108  G4ThreeVector pFinal1(std::sin(std::acos(cosTheta))*std::cos(phi), std::sin(std::acos(cosTheta))*std::sin(phi), cosTheta);
109 
110  // Three momentum in cm system
111  G4double pCM = std::sqrt( (S-(outm1+outm2)*(outm1+outm2)) * (S-(outm1-outm2)*(outm1-outm2)) /(4.*S));
112  pFinal1 = pFinal1 * pCM;
113  G4ThreeVector pFinal2 = -pFinal1;
114 
115  G4double eFinal1 = std::sqrt(pFinal1.mag2() + outm1*outm1);
116  G4double eFinal2 = std::sqrt(pFinal2.mag2() + outm2*outm2);
117 
118  G4LorentzVector p4Final1(pFinal1, eFinal1);
119  G4LorentzVector p4Final2(pFinal2, eFinal2);
120  p4Final1 = toCMS*p4Final1;
121  p4Final2 = toCMS*p4Final2;
122 
123 
124  // Lorentz transformation
125  G4LorentzRotation toLabFrame(p.boostVector());
126  p4Final1 *= toLabFrame;
127  p4Final2 *= toLabFrame;
128 
129  // Final tracks are copies of incoming ones, with modified 4-momenta
130 
131  G4double chargeBalance = OutputDefinitions[0]->GetPDGCharge()+OutputDefinitions[1]->GetPDGCharge();
132  chargeBalance-= trk1.GetDefinition()->GetPDGCharge();
133  chargeBalance-= trk2.GetDefinition()->GetPDGCharge();
134  if(std::abs(chargeBalance) >.1)
135  {
136  G4cout << "Charges in "<<typeid(*this).name()<<G4endl;
137  G4cout << OutputDefinitions[0]->GetPDGCharge()<<" "<<OutputDefinitions[0]->GetParticleName()
138  << OutputDefinitions[1]->GetPDGCharge()<<" "<<OutputDefinitions[1]->GetParticleName()
139  << trk1.GetDefinition()->GetPDGCharge()<<" "<<trk1.GetDefinition()->GetParticleName()
140  << trk2.GetDefinition()->GetPDGCharge()<<" "<<trk2.GetDefinition()->GetParticleName()<<G4endl;
141  }
142  G4KineticTrack* final1 = new G4KineticTrack(const_cast<G4ParticleDefinition *>(OutputDefinitions[0]), 0.0,
143  trk1.GetPosition(), p4Final1);
144  G4KineticTrack* final2 = new G4KineticTrack(const_cast<G4ParticleDefinition *>(OutputDefinitions[1]), 0.0,
145  trk2.GetPosition(), p4Final2);
146 
147  G4KineticTrackVector* finalTracks = new G4KineticTrackVector;
148 
149  finalTracks->push_back(final1);
150  finalTracks->push_back(final2);
151 
152  return finalTracks;
153 }
154 
155 
156 
157 double G4VScatteringCollision::SampleResonanceMass(const double poleMass,
158  const double gamma,
159  const double aMinMass,
160  const double maxMass) const
161 {
162  // Chooses a mass randomly between minMass and maxMass
163  // according to a Breit-Wigner function with constant
164  // width gamma and pole poleMass
165 
166  G4double minMass = aMinMass;
167  if (minMass > maxMass) G4cerr << "##################### SampleResonanceMass: particle out of mass range" << G4endl;
168  if(minMass > maxMass) minMass -= G4PionPlus::PionPlus()->GetPDGMass();
169  if(minMass > maxMass) minMass = 0;
170 
171  if (gamma < 1E-10*GeV)
172  return std::max(minMass,std::min(maxMass, poleMass));
173  else {
174  double fmin = BrWigInt0(minMass, gamma, poleMass);
175  double fmax = BrWigInt0(maxMass, gamma, poleMass);
176  double f = fmin + (fmax-fmin)*G4UniformRand();
177  return BrWigInv(f, gamma, poleMass);
178  }
179 }
180 
183 }
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
const G4ThreeVector & GetPosition() const
G4double GetActualMass() const
virtual G4double CosTheta(G4double s, G4double m1, G4double m2) const =0
const G4String & GetParticleName() const
void establish_G4MT_TLS_G4VCollision()
G4ParticleDefinition * GetDefinition() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
virtual G4double Phi() const
virtual const G4VAngularDistribution * GetAngularDistribution() const
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static const double GeV
Definition: G4SIunits.hh:196
G4VAngularDistribution * theAngularDistribution
double SampleResonanceMass(const double poleMass, const double width, const double minMass, const double maxMass) const
virtual const std::vector< const G4ParticleDefinition * > & GetOutgoingParticles() const =0
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4KineticTrackVector * FinalState(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double BrWigInt0(const double x, const double gamma, const double m0) const
#define G4endl
Definition: G4ios.hh:61
double BrWigInv(const double x, const double gamma, const double m0) const
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
const G4LorentzVector & Get4Momentum() const
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector