source: trunk/source/processes/electromagnetic/highenergy/src/G4eeTo3PiModel.cc@ 1066

Last change on this file since 1066 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 6.1 KB
RevLine 
[968]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//
[1055]26// $Id: G4eeTo3PiModel.cc,v 1.2 2009/02/20 16:38:33 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
[968]28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class header file
32//
33//
34// File name: G4eeTo3PiModel
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 25.10.2003
39//
40// Modifications:
41//
42//
43// -------------------------------------------------------------------
44//
45
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50#include "G4eeTo3PiModel.hh"
51#include "Randomize.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
64using namespace std;
65
66G4eeTo3PiModel::G4eeTo3PiModel(G4eeCrossSections* cr):
67 cross(cr)
68{
69 massPi = G4PionPlus::PionPlus()->GetPDGMass();
70 massPi0 = G4PionZero::PionZero()->GetPDGMass();
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
79G4eeTo3PiModel::~G4eeTo3PiModel()
80{
81 G4cout << "### G4eeTo3PiModel::~G4eeTo3PiModel: gmax= "
82 << gmax << " gcash= " << gcash << G4endl;
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
[1055]87G4double G4eeTo3PiModel::ThresholdEnergy() const
88{
89 return std::max(LowEnergy(),2.0*massPi + massPi0);
90}
91
92//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93
94G4double G4eeTo3PiModel::PeakEnergy() const
95{
96 G4double e = massOm;
97 if(HighEnergy() > massPhi) e = massPhi;
98 return e;
99}
100
101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102
103G4double G4eeTo3PiModel::ComputeCrossSection(G4double e) const
104{
105 G4double ee = std::min(HighEnergy(),e);
106 return cross->CrossSection3pi(ee);
107}
108
109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110
[968]111G4PhysicsVector* G4eeTo3PiModel::PhysicsVector(G4double emin,
112 G4double emax) const
113{
114 G4double tmin = std::max(emin, ThresholdEnergy());
115 G4double tmax = std::max(tmin, emax);
116 G4int nbins = (G4int)((tmax - tmin)/(1.*MeV));
117 G4PhysicsVector* v = new G4PhysicsLinearVector(emin,emax,nbins);
118 v->SetSpline(true);
119 return v;
120}
121
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123
124void G4eeTo3PiModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
125 G4double e, const G4ThreeVector& direction)
126{
127 if(e < ThresholdEnergy()) return;
128
129 G4double x0 = massPi0/e;
130 G4double x1 = massPi/e;
131
132 G4LorentzVector w0, w1, w2;
133 G4ThreeVector dir0, dir1, mom, mom1, mom2;
134 G4double e0, p0, e2, p, g, m01, m02, m12;
135
136 // max pi0 energy
137 G4double edel = 0.5*e*(1.0 + x0*x0 - 4.0*x1*x1) - massPi0;
138
139 do {
140 // pi0 sample
141 e0 = edel*G4UniformRand() + massPi0;
142 p0 = sqrt(e0 - massPi0*massPi0);
143 dir0 = G4RandomDirection();
144 w0 = G4LorentzVector(p0*dir0.x(),p0*dir0.y(),p0*dir0.z(),e0);
145
146 // pi+pi- pair
147 w1 = G4LorentzVector(-p0*dir0.x(),-p0*dir0.y(),-p0*dir0.z(),e-e0);
148 G4ThreeVector bst = w1.boostVector();
149 e2 = 0.25*w1.m2();
150
151 // pi+
152 p = sqrt(e2 - massPi*massPi);
153 dir1 = G4RandomDirection();
154 w2 = G4LorentzVector(p*dir1.x(),p*dir1.y(),p*dir1.z(),sqrt(e2));
155 w2.boost(bst);
156 mom2 = w2.vect();
157
158 // pi-
159 w1 -= w2;
160 mom1 = w2.vect();
161
162 m01 = w0*w1;
163 m02 = w0*w2;
164 m12 = w1*w2;
165
166 mom = mom1*mom2;
167 g = mom.mag2()*norm( 1.0/cross->DpRho(m01) + 1.0/cross->DpRho(m02)
168 + 1.0/cross->DpRho(m12) );
169 if(g > gmax) {
170 G4cout << "G4eeTo3PiModel::SampleSecondaries WARNING matrix element g= "
171 << g << " > " << gmax << " (majoranta)" << G4endl;
172 }
173 if(g > gcash) gcash = g;
174
175 } while( gmax*G4UniformRand() > g );
176
177 w0.rotateUz(direction);
178 w1.rotateUz(direction);
179 w2.rotateUz(direction);
180
181 // create G4DynamicParticle objects
182 G4DynamicParticle* dp0 =
183 new G4DynamicParticle(G4PionZero::PionZero(), w0);
184 G4DynamicParticle* dp1 =
185 new G4DynamicParticle(G4PionPlus::PionPlus(), w1);
186 G4DynamicParticle* dp2 =
187 new G4DynamicParticle(G4PionMinus::PionMinus(), w2);
188 newp->push_back(dp0);
189 newp->push_back(dp1);
190 newp->push_back(dp2);
191}
192
193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
194
Note: See TracBrowser for help on using the repository browser.