ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eeTo3PiModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eeTo3PiModel.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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 //
32 // File name: G4eeTo3PiModel
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 25.10.2003
37 //
38 // Modifications:
39 //
40 //
41 // -------------------------------------------------------------------
42 //
43 
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
48 #include "G4eeTo3PiModel.hh"
49 #include "Randomize.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4PionPlus.hh"
52 #include "G4PionMinus.hh"
53 #include "G4PionZero.hh"
54 #include "G4DynamicParticle.hh"
55 #include "G4PhysicsVector.hh"
56 #include "G4PhysicsLinearVector.hh"
57 #include "G4eeCrossSections.hh"
58 #include "G4RandomDirection.hh"
59 #include <complex>
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 
63 using namespace std;
64 
66  G4double maxkinEnergy,
67  G4double binWidth)
68 : G4Vee2hadrons(cr,
69  0.41612*GeV, //threshold
70  maxkinEnergy,
71  binWidth)
72 {
73  G4cout << "####G4eeTo3PiModel####" << G4endl;
74 
77  massOm = 782.62*MeV;
78  massPhi = 1019.46*MeV;
79  gmax = 3.0e-8;
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
85 {}
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
90 {
91  G4double e = massOm;
92  if(HighEnergy() > massPhi) { e = massPhi; }
93  return e;
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
99 {
100  return cross->CrossSection3pi(e);
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 
105 void G4eeTo3PiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
106  G4double e, const G4ThreeVector& direction)
107 {
108 
109  G4double x0 = massPi0/e;
110  G4double x1 = massPi/e;
111 
112  G4LorentzVector w0(0.,0.,0.,0.), w1(0.,0.,0.,0.), w2(0.,0.,0.,0.);
113  G4ThreeVector dir0, dir1;
114  G4double e0, p0, e2, p, gg, m01, m02, m12;
115 
116  // max pi0 energy
117  G4double edel = 0.5*e*(1.0 + x0*x0 - 4.0*x1*x1) - massPi0;
118 
119  const G4int nmax = 200;
120  G4int nn = 0;
121  do {
122  ++nn;
123  // pi0 sample
124  e0 = edel*G4UniformRand() + massPi0;
125  p0 = sqrt(e0*e0 - massPi0*massPi0);
126  dir0 = G4RandomDirection();
127  w0 = G4LorentzVector(p0*dir0.x(),p0*dir0.y(),p0*dir0.z(),e0);
128 
129  // pi+pi- pair
130  w1 = G4LorentzVector(-p0*dir0.x(),-p0*dir0.y(),-p0*dir0.z(),e-e0);
131  G4ThreeVector bst = w1.boostVector();
132  e2 = 0.25*w1.m2();
133 
134  // pi+
135  p = sqrt(e2 - massPi*massPi);
136  dir1 = G4RandomDirection();
137  w2 = G4LorentzVector(p*dir1.x(),p*dir1.y(),p*dir1.z(),sqrt(e2));
138  // pi-
139  w1.set(-w2.px(), -w2.py(), -w2.pz(), w2.e());
140 
141  w1.boost(bst);
142  w2.boost(bst);
143 
144  G4double px2 = w2.x();
145  G4double py2 = w2.y();
146  G4double pz2 = w2.z();
147 
148  G4double px1 = w1.x();
149  G4double py1 = w1.y();
150  G4double pz1 = w1.z();
151 
152  m01 = w0*w1;
153  m02 = w0*w2;
154  m12 = w1*w2;
155 
156  G4double px = py1*pz2 - py2*pz1;
157  G4double py = pz1*px2 - pz2*px1;
158  G4double pz = px1*py2 - px2*py1;
159 
160  gg = (px*px + py*py + pz*pz)*
161  norm( 1.0/cross->DpRho(m01) + 1.0/cross->DpRho(m02)
162  + 1.0/cross->DpRho(m12) );
163 
164  if(gg > gmax) {
165  G4cout << "G4eeTo3PiModel::SampleSecondaries WARNING matrix element g= "
166  << gg << " > " << gmax << " (majoranta)" << G4endl;
167  gmax = gg;
168  }
169  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
170  } while( gmax*G4UniformRand() > gg || nn < nmax);
171 
172  w0.rotateUz(direction);
173  w1.rotateUz(direction);
174  w2.rotateUz(direction);
175 
176  // create G4DynamicParticle objects
177  G4DynamicParticle* dp0 =
179  G4DynamicParticle* dp1 =
181  G4DynamicParticle* dp2 =
183  newp->push_back(dp0);
184  newp->push_back(dp1);
185  newp->push_back(dp2);
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189