| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file eventgenerator/particleGun/src/Prima <<
27 /// \brief Implementation of the PrimaryGenera <<
28 // 26 //
29 // << 27 // $Id: PrimaryGeneratorAction2.cc,v 1.3 2010-07-16 07:37:48 maire Exp $
30 // << 28 // GEANT4 tag $Name: not supported by cvs2svn $
31 //....oooOO0OOooo........oooOO0OOooo........oo << 29 //
32 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 32
34 #include "PrimaryGeneratorAction2.hh" 33 #include "PrimaryGeneratorAction2.hh"
35 <<
36 #include "PrimaryGeneratorAction.hh" 34 #include "PrimaryGeneratorAction.hh"
37 35
38 #include "G4Event.hh" 36 #include "G4Event.hh"
39 #include "G4ParticleDefinition.hh" <<
40 #include "G4ParticleGun.hh" 37 #include "G4ParticleGun.hh"
41 #include "G4ParticleTable.hh" 38 #include "G4ParticleTable.hh"
42 #include "G4PhysicalConstants.hh" << 39 #include "G4ParticleDefinition.hh"
43 #include "G4SystemOfUnits.hh" <<
44 #include "Randomize.hh" 40 #include "Randomize.hh"
45 41
46 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 43
48 PrimaryGeneratorAction2::PrimaryGeneratorActio << 44 PrimaryGeneratorAction2::PrimaryGeneratorAction2(G4ParticleGun* gun)
49 { << 45 : particleGun(gun)
>> 46 {
50 // energy distribution 47 // energy distribution
51 // 48 //
52 InitFunction(); 49 InitFunction();
53 } 50 }
54 51
55 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 53
57 void PrimaryGeneratorAction2::GeneratePrimarie << 54 PrimaryGeneratorAction2::~PrimaryGeneratorAction2()
58 { << 55 { }
59 // uniform solid angle <<
60 G4double cosAlpha = 1. - 2 * G4UniformRand() <<
61 G4double sinAlpha = std::sqrt(1. - cosAlpha <<
62 G4double psi = twopi * G4UniformRand(); // <<
63 G4ThreeVector dir(sinAlpha * std::cos(psi), <<
64 <<
65 fParticleGun->SetParticleMomentumDirection(d <<
66 56
67 // set energy from a tabulated distribution << 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 // <<
69 // G4double energy = RejectAccept(); <<
70 G4double energy = InverseCumul(); <<
71 fParticleGun->SetParticleEnergy(energy); <<
72 58
73 // create vertex << 59 void PrimaryGeneratorAction2::GeneratePrimaries(G4Event* anEvent)
74 // << 60 {
75 fParticleGun->GeneratePrimaryVertex(anEvent) << 61 //cosAlpha uniform in [cos(0), cos(pi)]
>> 62 G4double cosAlpha = 1. - 2*G4UniformRand();
>> 63 G4double sinAlpha = std::sqrt(1. - cosAlpha*cosAlpha);
>> 64 G4double psi = twopi*G4UniformRand(); //psi uniform in [0, 2*pi]
>> 65 G4ThreeVector dir(sinAlpha*std::cos(psi),sinAlpha*std::sin(psi),cosAlpha);
>> 66
>> 67 particleGun->SetParticleMomentumDirection(dir);
>> 68
>> 69 //set energy from a tabulated distribution
>> 70 //
>> 71 //G4double energy = RejectAccept();
>> 72 G4double energy = InverseCumul();
>> 73 particleGun->SetParticleEnergy(energy);
>> 74
>> 75 //create vertex
>> 76 //
>> 77 particleGun->GeneratePrimaryVertex(anEvent);
76 } 78 }
77 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 80
79 void PrimaryGeneratorAction2::InitFunction() 81 void PrimaryGeneratorAction2::InitFunction()
80 { 82 {
81 // tabulated function << 83 // tabulated function
82 // Y is assumed positive, linear per segment << 84 // f is assumed positive, linear per segment, continuous
83 // 85 //
84 fNPoints = 16; << 86 nPoints = 16;
85 const G4double xx[] = {37 * keV, 39 * keV, << 87 const G4double xx[] =
86 71 * keV, 75 * keV, << 88 { 37*keV, 39*keV, 45*keV, 51*keV, 57*keV, 69*keV, 71*keV, 75*keV,
87 125 * keV, 145 * keV, << 89 83*keV, 91*keV, 97*keV, 107*keV, 125*keV, 145*keV, 159*keV, 160*keV };
88 << 90
89 const G4double yy[] = {0.000, 0.077, 0.380 << 91 const G4double ff[] =
90 17.644, 18.518, 17.77 << 92 { 0.000, 0.077, 0.380, 2.044, 5.535, 15.077, 12.443, 14.766,
91 << 93 17.644, 18.518, 17.772, 14.776, 8.372, 3.217, 0.194, 0.000 };
92 // copy arrays in std::vector and compute fM << 94
93 // << 95 //copy arrays in std::vector and compute fMax
94 fX.resize(fNPoints); << 96 //
95 fY.resize(fNPoints); << 97 x.resize(nPoints); f.resize(nPoints);
96 fYmax = 0.; << 98 fMax = 0.;
97 for (G4int j = 0; j < fNPoints; j++) { << 99 for (G4int j=0; j<nPoints; j++) {
98 fX[j] = xx[j]; << 100 x[j] = xx[j]; f[j] = ff[j];
99 fY[j] = yy[j]; << 101 if (fMax < f[j]) fMax = f[j];
100 if (fYmax < fY[j]) fYmax = fY[j]; <<
101 }; 102 };
102 << 103
103 // compute slopes << 104 //compute slopes
104 // 105 //
105 fSlp.resize(fNPoints); << 106 a.resize(nPoints);
106 for (G4int j = 0; j < fNPoints - 1; j++) { << 107 for (G4int j=0; j<nPoints-1; j++) {
107 fSlp[j] = (fY[j + 1] - fY[j]) / (fX[j + 1] << 108 a[j] = (f[j+1] - f[j])/(x[j+1] - x[j]);
108 }; 109 };
109 << 110
110 // compute cumulative function << 111 //compute cumulative function
111 // 112 //
112 fYC.resize(fNPoints); << 113 Fc.resize(nPoints);
113 fYC[0] = 0.; << 114 Fc[0] = 0.;
114 for (G4int j = 1; j < fNPoints; j++) { << 115 for (G4int j=1; j<nPoints; j++) {
115 fYC[j] = fYC[j - 1] + 0.5 * (fY[j] + fY[j << 116 Fc[j] = Fc[j-1] + 0.5*(f[j] + f[j-1])*(x[j] - x[j-1]);
116 }; << 117 };
117 } 118 }
118 119
119 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 121
121 G4double PrimaryGeneratorAction2::RejectAccept 122 G4double PrimaryGeneratorAction2::RejectAccept()
122 { 123 {
123 // tabulated function << 124 // tabulated function
124 // Y is assumed positive, linear per segment << 125 // f is assumed positive, linear per segment, continuous
125 // (see Particle Data Group: pdg.lbl.gov --> << 126 //
126 // << 127 G4double x_rndm = 0., y_rndm = 0., f_inter = -1.;
127 G4double Xrndm = 0., Yrndm = 0., Yinter = -1 << 128
128 << 129 while (y_rndm > f_inter) {
129 while (Yrndm > Yinter) { << 130 //choose a point randomly
130 // choose a point randomly << 131 x_rndm = x[0] + G4UniformRand()*(x[nPoints-1] - x[0]);
131 Xrndm = fX[0] + G4UniformRand() * (fX[fNPo << 132 y_rndm = G4UniformRand()*fMax;
132 Yrndm = G4UniformRand() * fYmax; << 133 //find bin
133 // find bin << 134 G4int j = nPoints-2;
134 G4int j = fNPoints - 2; << 135 while ((x[j] > x_rndm) && (j > 0)) j--;
135 while ((fX[j] > Xrndm) && (j > 0)) << 136 //compute f(x_rndm) by linear interpolation
136 j--; << 137 f_inter = f[j] + a[j]*(x_rndm - x[j]);
137 // compute Y(x_rndm) by linear interpolati <<
138 Yinter = fY[j] + fSlp[j] * (Xrndm - fX[j]) <<
139 }; 138 };
140 return Xrndm; << 139 return x_rndm;
141 } 140 }
142 141
143 //....oooOO0OOooo........oooOO0OOooo........oo 142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 143
145 G4double PrimaryGeneratorAction2::InverseCumul 144 G4double PrimaryGeneratorAction2::InverseCumul()
146 { 145 {
147 // tabulated function 146 // tabulated function
148 // Y is assumed positive, linear per segment << 147 // f is assumed positive, linear per segment, continuous
149 // --> cumulative function is second order p 148 // --> cumulative function is second order polynomial
150 // (see Particle Data Group: pdg.lbl.gov --> << 149
151 << 150 //choose y randomly
152 // choose y randomly << 151 G4double y_rndm = G4UniformRand()*Fc[nPoints-1];
153 G4double Yrndm = G4UniformRand() * fYC[fNPoi << 152 //find bin
154 // find bin << 153 G4int j = nPoints-2;
155 G4int j = fNPoints - 2; << 154 while ((Fc[j] > y_rndm) && (j > 0)) j--;
156 while ((fYC[j] > Yrndm) && (j > 0)) << 155 //y_rndm --> x_rndm : Fc(x) is second order polynomial
157 j--; << 156 G4double x_rndm = x[j];
158 // y_rndm --> x_rndm : fYC(x) is second ord << 157 G4double aa = a[j];
159 G4double Xrndm = fX[j]; << 158 if (aa != 0.) {
160 G4double a = fSlp[j]; << 159 G4double b = f[j]/aa, c = 2*(y_rndm - Fc[j])/aa;
161 if (a != 0.) { << 160 G4double delta = b*b + c;
162 G4double b = fY[j] / a, c = 2 * (Yrndm - f << 161 G4int sign = 1; if (aa < 0.) sign = -1;
163 G4double delta = b * b + c; << 162 x_rndm += sign*std::sqrt(delta) - b;
164 G4int sign = 1; << 163 } else if (f[j] > 0.) {
165 if (a < 0.) sign = -1; << 164 x_rndm += (y_rndm - Fc[j])/f[j];
166 Xrndm += sign * std::sqrt(delta) - b; <<
167 } <<
168 else if (fY[j] > 0.) { <<
169 Xrndm += (Yrndm - fYC[j]) / fY[j]; <<
170 }; 165 };
171 return Xrndm; << 166 return x_rndm;
172 } 167 }
173 168
174 //....oooOO0OOooo........oooOO0OOooo........oo 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
175 170