ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VScatteringCollision.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VScatteringCollision.cc
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(std::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 {
185 }