Geant4  10.00.p01
G4eeTo3PiModel.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 // $Id: G4eeTo3PiModel.cc 66241 2012-12-13 18:34:42Z gunter $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class header file
31 //
32 //
33 // File name: G4eeTo3PiModel
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 25.10.2003
38 //
39 // Modifications:
40 //
41 //
42 // -------------------------------------------------------------------
43 //
44 
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 #include "G4eeTo3PiModel.hh"
50 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4PionPlus.hh"
53 #include "G4PionMinus.hh"
54 #include "G4PionZero.hh"
55 #include "G4DynamicParticle.hh"
56 #include "G4PhysicsVector.hh"
57 #include "G4PhysicsLinearVector.hh"
58 #include "G4eeCrossSections.hh"
59 #include "G4RandomDirection.hh"
60 #include <complex>
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
64 using namespace std;
65 
67  cross(cr)
68 {
71  massOm = 782.62*MeV;
72  massPhi = 1019.46*MeV;
73  gcash = 0.0;
74  gmax = 1.0;
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 
80 {}
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
85 {
86  return std::max(LowEnergy(),2.0*massPi + massPi0);
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
92 {
93  G4double e = massOm;
94  if(HighEnergy() > massPhi) e = massPhi;
95  return e;
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 
101 {
102  G4double ee = std::min(HighEnergy(),e);
103  return cross->CrossSection3pi(ee);
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 
109  G4double emax) const
110 {
111  G4double tmin = std::max(emin, ThresholdEnergy());
112  G4double tmax = std::max(tmin, emax);
113  G4int nbins = (G4int)((tmax - tmin)/(1.*MeV));
114  G4PhysicsVector* v = new G4PhysicsLinearVector(emin,emax,nbins);
115  v->SetSpline(true);
116  return v;
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120 
121 void G4eeTo3PiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
122  G4double e, const G4ThreeVector& direction)
123 {
124  if(e < ThresholdEnergy()) return;
125 
126  G4double x0 = massPi0/e;
127  G4double x1 = massPi/e;
128 
129  G4LorentzVector w0, w1, w2;
130  G4ThreeVector dir0, dir1;
131  G4double e0, p0, e2, p, gg, m01, m02, m12;
132 
133  // max pi0 energy
134  G4double edel = 0.5*e*(1.0 + x0*x0 - 4.0*x1*x1) - massPi0;
135 
136  do {
137  // pi0 sample
138  e0 = edel*G4UniformRand() + massPi0;
139  p0 = sqrt(e0 - massPi0*massPi0);
140  dir0 = G4RandomDirection();
141  w0 = G4LorentzVector(p0*dir0.x(),p0*dir0.y(),p0*dir0.z(),e0);
142 
143  // pi+pi- pair
144  w1 = G4LorentzVector(-p0*dir0.x(),-p0*dir0.y(),-p0*dir0.z(),e-e0);
145  G4ThreeVector bst = w1.boostVector();
146  e2 = 0.25*w1.m2();
147 
148  // pi+
149  p = sqrt(e2 - massPi*massPi);
150  dir1 = G4RandomDirection();
151  w2 = G4LorentzVector(p*dir1.x(),p*dir1.y(),p*dir1.z(),sqrt(e2));
152  w2.boost(bst);
153  G4double px2 = w2.x();
154  G4double py2 = w2.y();
155  G4double pz2 = w2.z();
156 
157  // pi-
158  w1 -= w2;
159  G4double px1 = w1.x();
160  G4double py1 = w1.y();
161  G4double pz1 = w1.z();
162 
163  m01 = w0*w1;
164  m02 = w0*w2;
165  m12 = w1*w2;
166 
167  G4double px = py1*pz2 - py2*pz1;
168  G4double py = pz1*px2 - pz2*px1;
169  G4double pz = px1*py2 - px2*py1;
170 
171  gg = (px*px + py*py + pz*pz)*
172  norm( 1.0/cross->DpRho(m01) + 1.0/cross->DpRho(m02)
173  + 1.0/cross->DpRho(m12) );
174  if(gg > gmax) {
175  G4cout << "G4eeTo3PiModel::SampleSecondaries WARNING matrix element g= "
176  << gg << " > " << gmax << " (majoranta)" << G4endl;
177  }
178  if(gg > gcash) gcash = gg;
179 
180  } while( gmax*G4UniformRand() > gg );
181 
182  w0.rotateUz(direction);
183  w1.rotateUz(direction);
184  w2.rotateUz(direction);
185 
186  // create G4DynamicParticle objects
187  G4DynamicParticle* dp0 =
189  G4DynamicParticle* dp1 =
191  G4DynamicParticle* dp2 =
193  newp->push_back(dp0);
194  newp->push_back(dp1);
195  newp->push_back(dp2);
196 }
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199 
G4double LowEnergy() const
static const double MeV
Definition: G4SIunits.hh:193
CLHEP::Hep3Vector G4ThreeVector
static const G4double e2
G4ThreeVector G4RandomDirection()
virtual G4PhysicsVector * PhysicsVector(G4double, G4double) const
G4double CrossSection3pi(G4double)
int G4int
Definition: G4Types.hh:78
void SetSpline(G4bool)
std::complex< G4double > DpRho(G4double e)
virtual G4double PeakEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
virtual G4double ThresholdEnergy() const
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
G4double GetPDGMass() const
virtual G4double ComputeCrossSection(G4double) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
G4double HighEnergy() const
G4eeCrossSections * cross
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, G4double, const G4ThreeVector &)
double G4double
Definition: G4Types.hh:76
G4eeTo3PiModel(G4eeCrossSections *)
virtual ~G4eeTo3PiModel()
CLHEP::HepLorentzVector G4LorentzVector