Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 91869 2015-08-07 15:21:02Z gcosmo $
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  G4double maxkinEnergy,
68  G4double binWidth)
69 : G4Vee2hadrons(cr,
70  0.41612*GeV, //threshold
71  maxkinEnergy,
72  binWidth)
73 {
74  G4cout << "####G4eeTo3PiModel####" << G4endl;
75 
76  massPi = G4PionPlus::PionPlus()->GetPDGMass();
77  massPi0 = G4PionZero::PionZero()->GetPDGMass();
78  massOm = 782.62*MeV;
79  massPhi = 1019.46*MeV;
80  gmax = 3.0e-8;
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
86 {}
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91 {
92  G4double e = massOm;
93  if(HighEnergy() > massPhi) { e = massPhi; }
94  return e;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
100 {
101  return cross->CrossSection3pi(e);
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 
106 void G4eeTo3PiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
107  G4double e, const G4ThreeVector& direction)
108 {
109 
110  G4double x0 = massPi0/e;
111  G4double x1 = massPi/e;
112 
113  G4LorentzVector w0(0.,0.,0.,0.), w1(0.,0.,0.,0.), w2(0.,0.,0.,0.);
114  G4ThreeVector dir0, dir1;
115  G4double e0, p0, e2, p, gg, m01, m02, m12;
116 
117  // max pi0 energy
118  G4double edel = 0.5*e*(1.0 + x0*x0 - 4.0*x1*x1) - massPi0;
119 
120  const G4int nmax = 200;
121  G4int nn = 0;
122  do {
123  ++nn;
124  // pi0 sample
125  e0 = edel*G4UniformRand() + massPi0;
126  p0 = sqrt(e0*e0 - massPi0*massPi0);
127  dir0 = G4RandomDirection();
128  w0 = G4LorentzVector(p0*dir0.x(),p0*dir0.y(),p0*dir0.z(),e0);
129 
130  // pi+pi- pair
131  w1 = G4LorentzVector(-p0*dir0.x(),-p0*dir0.y(),-p0*dir0.z(),e-e0);
132  G4ThreeVector bst = w1.boostVector();
133  e2 = 0.25*w1.m2();
134 
135  // pi+
136  p = sqrt(e2 - massPi*massPi);
137  dir1 = G4RandomDirection();
138  w2 = G4LorentzVector(p*dir1.x(),p*dir1.y(),p*dir1.z(),sqrt(e2));
139  // pi-
140  w1.set(-w2.px(), -w2.py(), -w2.pz(), w2.e());
141 
142  w1.boost(bst);
143  w2.boost(bst);
144 
145  G4double px2 = w2.x();
146  G4double py2 = w2.y();
147  G4double pz2 = w2.z();
148 
149  G4double px1 = w1.x();
150  G4double py1 = w1.y();
151  G4double pz1 = w1.z();
152 
153  m01 = w0*w1;
154  m02 = w0*w2;
155  m12 = w1*w2;
156 
157  G4double px = py1*pz2 - py2*pz1;
158  G4double py = pz1*px2 - pz2*px1;
159  G4double pz = px1*py2 - px2*py1;
160 
161  gg = (px*px + py*py + pz*pz)*
162  norm( 1.0/cross->DpRho(m01) + 1.0/cross->DpRho(m02)
163  + 1.0/cross->DpRho(m12) );
164 
165  if(gg > gmax) {
166  G4cout << "G4eeTo3PiModel::SampleSecondaries WARNING matrix element g= "
167  << gg << " > " << gmax << " (majoranta)" << G4endl;
168  gmax = gg;
169  }
170  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
171  } while( gmax*G4UniformRand() > gg || nn < nmax);
172 
173  w0.rotateUz(direction);
174  w1.rotateUz(direction);
175  w2.rotateUz(direction);
176 
177  // create G4DynamicParticle objects
178  G4DynamicParticle* dp0 =
180  G4DynamicParticle* dp1 =
182  G4DynamicParticle* dp2 =
184  newp->push_back(dp0);
185  newp->push_back(dp1);
186  newp->push_back(dp2);
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
190 
Hep3Vector boostVector() const
double x() const
G4eeTo3PiModel(G4eeCrossSections *, G4double, G4double)
const char * p
Definition: xmltok.h:285
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, G4double, const G4ThreeVector &) override
G4ThreeVector G4RandomDirection()
G4double CrossSection3pi(G4double)
int G4int
Definition: G4Types.hh:78
virtual G4double PeakEnergy() const override
double z() const
std::complex< G4double > DpRho(G4double e)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
virtual G4double ComputeCrossSection(G4double) const override
double py() const
HepLorentzVector & boost(double, double, double)
double px() const
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
const G4int nmax
HepLorentzVector & rotateUz(const Hep3Vector &)
static G4PionZero * PionZero()
Definition: G4PionZero.cc:108
void set(double x, double y, double z, double t)
G4double GetPDGMass() const
double m2() const
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
double y() const
static constexpr double GeV
Definition: G4SIunits.hh:217
double pz() const
G4eeCrossSections * cross
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double HighEnergy() const
double G4double
Definition: G4Types.hh:76
virtual ~G4eeTo3PiModel()
CLHEP::HepLorentzVector G4LorentzVector