ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PrimaryGeneratorAction2.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PrimaryGeneratorAction2.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 //
28 //
29 //
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
35 #include "PrimaryGeneratorAction.hh"
36 
37 #include "G4Event.hh"
38 #include "G4ParticleGun.hh"
39 #include "G4ParticleTable.hh"
40 #include "G4ParticleDefinition.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh"
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 
48 : fParticleGun(gun)
49 {
50  // energy distribution
51  //
52  InitFunction();
53 }
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
58 { }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
63 {
64  //cosAlpha uniform in [cos(0), cos(pi)]
65  G4double cosAlpha = 1. - 2*G4UniformRand();
66  G4double sinAlpha = std::sqrt(1. - cosAlpha*cosAlpha);
67  G4double psi = twopi*G4UniformRand(); //psi uniform in [0, 2*pi]
68  G4ThreeVector dir(sinAlpha*std::cos(psi),sinAlpha*std::sin(psi),cosAlpha);
69 
71 
72  //set energy from a tabulated distribution
73  //
74  //G4double energy = RejectAccept();
77 
78  //create vertex
79  //
81 }
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 
85 {
86  // tabulated function
87  // Y is assumed positive, linear per segment, continuous
88  //
89  fNPoints = 16;
90  const G4double xx[] =
91  { 37*keV, 39*keV, 45*keV, 51*keV, 57*keV, 69*keV, 71*keV, 75*keV,
92  83*keV, 91*keV, 97*keV, 107*keV, 125*keV, 145*keV, 159*keV, 160*keV };
93 
94  const G4double yy[] =
95  { 0.000, 0.077, 0.380, 2.044, 5.535, 15.077, 12.443, 14.766,
96  17.644, 18.518, 17.772, 14.776, 8.372, 3.217, 0.194, 0.000 };
97 
98  //copy arrays in std::vector and compute fMax
99  //
100  fX.resize(fNPoints); fY.resize(fNPoints);
101  fYmax = 0.;
102  for (G4int j=0; j<fNPoints; j++) {
103  fX[j] = xx[j]; fY[j] = yy[j];
104  if (fYmax < fY[j]) fYmax = fY[j];
105  };
106 
107  //compute slopes
108  //
109  fSlp.resize(fNPoints);
110  for (G4int j=0; j<fNPoints-1; j++) {
111  fSlp[j] = (fY[j+1] - fY[j])/(fX[j+1] - fX[j]);
112  };
113 
114  //compute cumulative function
115  //
116  fYC.resize(fNPoints);
117  fYC[0] = 0.;
118  for (G4int j=1; j<fNPoints; j++) {
119  fYC[j] = fYC[j-1] + 0.5*(fY[j] + fY[j-1])*(fX[j] - fX[j-1]);
120  };
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 
126 {
127  // tabulated function
128  // Y is assumed positive, linear per segment, continuous
129  //
130  G4double Xrndm = 0., Yrndm = 0., Yinter = -1.;
131 
132  while (Yrndm > Yinter) {
133  //choose a point randomly
134  Xrndm = fX[0] + G4UniformRand()*(fX[fNPoints-1] - fX[0]);
135  Yrndm = G4UniformRand()*fYmax;
136  //find bin
137  G4int j = fNPoints-2;
138  while ((fX[j] > Xrndm) && (j > 0)) j--;
139  //compute Y(x_rndm) by linear interpolation
140  Yinter = fY[j] + fSlp[j]*(Xrndm - fX[j]);
141  };
142  return Xrndm;
143 }
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 
148 {
149  // tabulated function
150  // Y is assumed positive, linear per segment, continuous
151  // --> cumulative function is second order polynomial
152 
153  //choose y randomly
154  G4double Yrndm = G4UniformRand()*fYC[fNPoints-1];
155  //find bin
156  G4int j = fNPoints-2;
157  while ((fYC[j] > Yrndm) && (j > 0)) j--;
158  //y_rndm --> x_rndm : fYC(x) is second order polynomial
159  G4double Xrndm = fX[j];
160  G4double a = fSlp[j];
161  if (a != 0.) {
162  G4double b = fY[j]/a, c = 2*(Yrndm - fYC[j])/a;
163  G4double delta = b*b + c;
164  G4int sign = 1; if (a < 0.) sign = -1;
165  Xrndm += sign*std::sqrt(delta) - b;
166  } else if (fY[j] > 0.) {
167  Xrndm += (Yrndm - fYC[j])/fY[j];
168  };
169  return Xrndm;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......