Geant4  10.02.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 {
57 }
58 
59 
61  const G4KineticTrack& trk2) const
62 {
63  const G4VAngularDistribution* angDistribution = GetAngularDistribution();
64  G4LorentzVector p = trk1.Get4Momentum() + trk2.Get4Momentum();
65  G4double sqrtS = p.m();
66  G4double S = sqrtS * sqrtS;
67 
68  std::vector<const G4ParticleDefinition*> OutputDefinitions = GetOutgoingParticles();
69  if (OutputDefinitions.size() != 2)
70  throw G4HadronicException(__FILE__, __LINE__, "G4VScatteringCollision: Too many output particles!");
71 
72  if (OutputDefinitions[0]->IsShortLived() && OutputDefinitions[1]->IsShortLived())
73  {
74  if(getenv("G4KCDEBUG")) G4cerr << "two shortlived for Type = "<<typeid(*this).name()<<G4endl;
75  // throw G4HadronicException(__FILE__, __LINE__, "G4VScatteringCollision: can't handle two shortlived particles!"); // @hpw@
76  }
77 
78  G4double outm1 = OutputDefinitions[0]->GetPDGMass();
79  G4double outm2 = OutputDefinitions[1]->GetPDGMass();
80 
81  if (OutputDefinitions[0]->IsShortLived())
82  {
83  outm1 = SampleResonanceMass(outm1,
84  OutputDefinitions[0]->GetPDGWidth(),
85  G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass(),
86  sqrtS-(G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass()));
87 
88  }
89  if (OutputDefinitions[1]->IsShortLived())
90  {
91  outm2 = SampleResonanceMass(outm2, OutputDefinitions[1]->GetPDGWidth(),
92  G4Neutron::NeutronDefinition()->GetPDGMass()+G4PionPlus::PionPlus()->GetPDGMass(),
93  sqrtS-outm1);
94  }
95 
96  // Angles of outgoing particles
97  G4double cosTheta = angDistribution->CosTheta(S, trk1.GetActualMass(), trk2.GetActualMass());
98  G4double phi = angDistribution->Phi();
99 
100  // Unit vector of three-momentum
101  G4LorentzRotation fromCMSFrame(p.boostVector());
102  G4LorentzRotation toCMSFrame(fromCMSFrame.inverse());
103  G4LorentzVector TempPtr = toCMSFrame*trk1.Get4Momentum();
104  G4LorentzRotation toZ;
105  toZ.rotateZ(-1*TempPtr.phi());
106  toZ.rotateY(-1*TempPtr.theta());
107  G4LorentzRotation toCMS(toZ.inverse());
108 
109  G4ThreeVector pFinal1(std::sin(std::acos(cosTheta))*std::cos(phi), std::sin(std::acos(cosTheta))*std::sin(phi), cosTheta);
110 
111  // Three momentum in cm system
112  G4double pCM = std::sqrt( (S-(outm1+outm2)*(outm1+outm2)) * (S-(outm1-outm2)*(outm1-outm2)) /(4.*S));
113  pFinal1 = pFinal1 * pCM;
114  G4ThreeVector pFinal2 = -pFinal1;
115 
116  G4double eFinal1 = std::sqrt(pFinal1.mag2() + outm1*outm1);
117  G4double eFinal2 = std::sqrt(pFinal2.mag2() + outm2*outm2);
118 
119  G4LorentzVector p4Final1(pFinal1, eFinal1);
120  G4LorentzVector p4Final2(pFinal2, eFinal2);
121  p4Final1 = toCMS*p4Final1;
122  p4Final2 = toCMS*p4Final2;
123 
124 
125  // Lorentz transformation
126  G4LorentzRotation toLabFrame(p.boostVector());
127  p4Final1 *= toLabFrame;
128  p4Final2 *= toLabFrame;
129 
130  // Final tracks are copies of incoming ones, with modified 4-momenta
131 
132  G4double chargeBalance = OutputDefinitions[0]->GetPDGCharge()+OutputDefinitions[1]->GetPDGCharge();
133  chargeBalance-= trk1.GetDefinition()->GetPDGCharge();
134  chargeBalance-= trk2.GetDefinition()->GetPDGCharge();
135  if(std::abs(chargeBalance) >.1)
136  {
137  G4cout << "Charges in "<<typeid(*this).name()<<G4endl;
138  G4cout << OutputDefinitions[0]->GetPDGCharge()<<" "<<OutputDefinitions[0]->GetParticleName()
139  << OutputDefinitions[1]->GetPDGCharge()<<" "<<OutputDefinitions[1]->GetParticleName()
140  << trk1.GetDefinition()->GetPDGCharge()<<" "<<trk1.GetDefinition()->GetParticleName()
141  << trk2.GetDefinition()->GetPDGCharge()<<" "<<trk2.GetDefinition()->GetParticleName()<<G4endl;
142  }
143  G4KineticTrack* final1 = new G4KineticTrack(OutputDefinitions[0], 0.0, trk1.GetPosition(), p4Final1);
144  G4KineticTrack* final2 = new G4KineticTrack(OutputDefinitions[1], 0.0, trk2.GetPosition(), p4Final2);
145 
146  G4KineticTrackVector* finalTracks = new G4KineticTrackVector;
147 
148  finalTracks->push_back(final1);
149  finalTracks->push_back(final2);
150 
151  return finalTracks;
152 }
153 
154 
155 
156 double G4VScatteringCollision::SampleResonanceMass(const double poleMass,
157  const double gamma,
158  const double aMinMass,
159  const double maxMass) const
160 {
161  // Chooses a mass randomly between minMass and maxMass
162  // according to a Breit-Wigner function with constant
163  // width gamma and pole poleMass
164 
165  G4double minMass = aMinMass;
166  if (minMass > maxMass) G4cerr << "##################### SampleResonanceMass: particle out of mass range" << G4endl;
167  if(minMass > maxMass) minMass -= G4PionPlus::PionPlus()->GetPDGMass();
168  if(minMass > maxMass) minMass = 0;
169 
170  if (gamma < 1E-10*GeV)
171  return std::max(minMass,std::min(maxMass, poleMass));
172  else {
173  double fmin = BrWigInt0(minMass, gamma, poleMass);
174  double fmax = BrWigInt0(maxMass, gamma, poleMass);
175  double f = fmin + (fmax-fmin)*G4UniformRand();
176  return BrWigInv(f, gamma, poleMass);
177  }
178 }
179 
181 {
184 }
double S(double temp)
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()
#define G4UniformRand()
Definition: Randomize.hh:97
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:214
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
const G4ParticleDefinition * GetDefinition() const
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector