Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/event/src/G4SPSAngDistribution.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /event/src/G4SPSAngDistribution.cc (Version 11.3.0) and /event/src/G4SPSAngDistribution.cc (Version 9.6)


  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 // G4SPSAngDistribution class implementation   <<  26 ///////////////////////////////////////////////////////////////////////////////
                                                   >>  27 //
                                                   >>  28 // MODULE:        G4SPSAngDistribution.cc
                                                   >>  29 //
                                                   >>  30 // Version:      1.0
                                                   >>  31 // Date:         5/02/04
                                                   >>  32 // Author:       Fan Lei 
                                                   >>  33 // Organisation: QinetiQ ltd.
                                                   >>  34 // Customer:     ESA/ESTEC
                                                   >>  35 //
                                                   >>  36 ///////////////////////////////////////////////////////////////////////////////
                                                   >>  37 //
                                                   >>  38 //
                                                   >>  39 // CHANGE HISTORY
                                                   >>  40 // --------------
                                                   >>  41 //
                                                   >>  42 //
                                                   >>  43 // Version 1.0, 05/02/2004, Fan Lei, Created.
                                                   >>  44 //    Based on the G4GeneralParticleSource class in Geant4 v6.0
                                                   >>  45 //
                                                   >>  46 ///////////////////////////////////////////////////////////////////////////////
 27 //                                                 47 //
 28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004  << 
 29 // Customer: ESA/ESTEC                         << 
 30 // Revisions: Andrea Dotti, SLAC               << 
 31 // ------------------------------------------- << 
 32                                                    48 
 33 #include "G4SPSAngDistribution.hh"                 49 #include "G4SPSAngDistribution.hh"
 34                                                    50 
 35 #include "Randomize.hh"                            51 #include "Randomize.hh"
 36 #include "G4PhysicalConstants.hh"                  52 #include "G4PhysicalConstants.hh"
 37                                                    53 
 38 G4SPSAngDistribution::G4SPSAngDistribution()   <<  54 G4SPSAngDistribution::G4SPSAngDistribution()
                                                   >>  55   : Theta(0.), Phi(0.)
 39 {                                                  56 {
 40   // Angular distribution Variables                57   // Angular distribution Variables
 41   G4ThreeVector zero;                              58   G4ThreeVector zero;
 42   particle_momentum_direction = G4ParticleMome     59   particle_momentum_direction = G4ParticleMomentum(0,0,-1);
 43                                                    60 
 44   AngDistType = "planar";                          61   AngDistType = "planar"; 
 45   AngRef1 = CLHEP::HepXHat;                        62   AngRef1 = CLHEP::HepXHat;
 46   AngRef2 = CLHEP::HepYHat;                        63   AngRef2 = CLHEP::HepYHat;
 47   AngRef3 = CLHEP::HepZHat;                        64   AngRef3 = CLHEP::HepZHat;
 48   MinTheta = 0.;                                   65   MinTheta = 0.;
 49   MaxTheta = pi;                                   66   MaxTheta = pi;
 50   MinPhi = 0.;                                     67   MinPhi = 0.;
 51   MaxPhi = twopi;                                  68   MaxPhi = twopi;
 52   DR = 0.;                                         69   DR = 0.;
 53   DX = 0.;                                         70   DX = 0.;
 54   DY = 0.;                                         71   DY = 0.;
 55   FocusPoint = G4ThreeVector(0., 0., 0.);          72   FocusPoint = G4ThreeVector(0., 0., 0.);
 56   UserDistType = "NULL";                           73   UserDistType = "NULL";
 57   UserWRTSurface = true;                           74   UserWRTSurface = true;
 58   UserAngRef = false;                              75   UserAngRef = false;
 59   IPDFThetaExist = false;                          76   IPDFThetaExist = false;
 60   IPDFPhiExist = false;                            77   IPDFPhiExist = false;
 61   verbosityLevel = 0;                          <<  78   verbosityLevel = 0 ;
 62                                                << 
 63   G4MUTEXINIT(mutex);                          << 
 64 }                                                  79 }
 65                                                    80 
 66 G4SPSAngDistribution::~G4SPSAngDistribution()      81 G4SPSAngDistribution::~G4SPSAngDistribution()
 67 {                                              <<  82 {}
 68     G4MUTEXDESTROY(mutex);                     << 
 69 }                                              << 
 70                                                    83 
 71 void G4SPSAngDistribution::SetAngDistType(cons <<  84 //
                                                   >>  85 void G4SPSAngDistribution::SetAngDistType(G4String atype)
 72 {                                                  86 {
 73   G4AutoLock l(&mutex);                        << 
 74   if(atype != "iso" && atype != "cos" && atype     87   if(atype != "iso" && atype != "cos" && atype != "user" && atype != "planar"
 75      && atype != "beam1d" && atype != "beam2d"     88      && atype != "beam1d" && atype != "beam2d"  && atype != "focused")
 76   {                                            <<  89     G4cout << "Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user" << G4endl;
 77     G4cout << "Error, distribution must be iso << 
 78            << G4endl;                          << 
 79   }                                            << 
 80   else                                             90   else
 81   {                                            << 
 82     AngDistType = atype;                           91     AngDistType = atype;
 83   }                                            <<  92   if (AngDistType == "cos") MaxTheta = pi/2. ;
 84   if (AngDistType == "cos")  { MaxTheta = pi/2 <<  93   if (AngDistType == "user") {
 85   if (AngDistType == "user")                   <<  94     UDefThetaH = IPDFThetaH = ZeroPhysVector ;
 86   {                                            <<  95     IPDFThetaExist = false ;
 87     UDefThetaH = IPDFThetaH = ZeroPhysVector;  <<  96     UDefPhiH = IPDFPhiH = ZeroPhysVector ;
 88     IPDFThetaExist = false;                    <<  97     IPDFPhiExist = false ;
 89     UDefPhiH = IPDFPhiH = ZeroPhysVector;      << 
 90     IPDFPhiExist = false;                      << 
 91   }                                                98   }
 92 }                                                  99 }
 93                                                   100 
 94 void G4SPSAngDistribution::DefineAngRefAxes(co << 101 void G4SPSAngDistribution::DefineAngRefAxes(G4String refname, G4ThreeVector ref)
 95                                             co << 
 96 {                                                 102 {
 97   G4AutoLock l(&mutex);                        << 103   if(refname == "angref1")
 98   if (refname == "angref1")                    << 
 99     AngRef1 = ref.unit(); // x'                   104     AngRef1 = ref.unit(); // x'
100   else if (refname == "angref2")               << 105   else if(refname == "angref2")
101     AngRef2 = ref.unit(); // vector in x'y' pl    106     AngRef2 = ref.unit(); // vector in x'y' plane
102                                                   107 
103   // User defines x' (AngRef1) and a vector in    108   // User defines x' (AngRef1) and a vector in the x'y'
104   // plane (AngRef2). Then, AngRef1 x AngRef2     109   // plane (AngRef2). Then, AngRef1 x AngRef2 = AngRef3
105   // the z' vector. Then, AngRef3 x AngRef1 =     110   // the z' vector. Then, AngRef3 x AngRef1 = AngRef2
106   // which will now be y'.                        111   // which will now be y'.
107                                                   112 
108   AngRef3 = AngRef1.cross(AngRef2); // z'         113   AngRef3 = AngRef1.cross(AngRef2); // z'
109   AngRef2 = AngRef3.cross(AngRef1); // y'         114   AngRef2 = AngRef3.cross(AngRef1); // y'
110   UserAngRef = true ;                             115   UserAngRef = true ;
111   if(verbosityLevel == 2)                         116   if(verbosityLevel == 2)
112   {                                            << 117     {
113     G4cout << "Angular distribution rotation a << 118       G4cout << "Angular distribution rotation axes " << AngRef1 << " " << AngRef2 << " " << AngRef3 << G4endl;
114            << " " << AngRef2 << " " << AngRef3 << 119     }
115   }                                            << 
116 }                                                 120 }
117                                                   121 
118 void G4SPSAngDistribution::SetMinTheta(G4doubl    122 void G4SPSAngDistribution::SetMinTheta(G4double mint)
119 {                                                 123 {
120   G4AutoLock l(&mutex);                        << 
121   MinTheta = mint;                                124   MinTheta = mint;
122 }                                                 125 }
123                                                   126 
124 void G4SPSAngDistribution::SetMinPhi(G4double     127 void G4SPSAngDistribution::SetMinPhi(G4double minp)
125 {                                                 128 {
126   G4AutoLock l(&mutex);                        << 
127   MinPhi = minp;                                  129   MinPhi = minp;
128 }                                                 130 }
129                                                   131 
130 void G4SPSAngDistribution::SetMaxTheta(G4doubl    132 void G4SPSAngDistribution::SetMaxTheta(G4double maxt)
131 {                                                 133 {
132   G4AutoLock l(&mutex);                        << 
133   MaxTheta = maxt;                                134   MaxTheta = maxt;
134 }                                                 135 }
135                                                   136 
136 void G4SPSAngDistribution::SetMaxPhi(G4double     137 void G4SPSAngDistribution::SetMaxPhi(G4double maxp)
137 {                                                 138 {
138   G4AutoLock l(&mutex);                        << 
139   MaxPhi = maxp;                                  139   MaxPhi = maxp;
140 }                                                 140 }
141                                                   141 
142 void G4SPSAngDistribution::SetBeamSigmaInAngR(    142 void G4SPSAngDistribution::SetBeamSigmaInAngR(G4double r)
143 {                                                 143 {
144   G4AutoLock l(&mutex);                        << 
145   DR = r;                                         144   DR = r;
146 }                                                 145 }
147                                                   146 
148 void G4SPSAngDistribution::SetBeamSigmaInAngX(    147 void G4SPSAngDistribution::SetBeamSigmaInAngX(G4double r)
149 {                                                 148 {
150   G4AutoLock l(&mutex);                        << 
151   DX = r;                                         149   DX = r;
152 }                                                 150 }
153                                                   151 
154 void G4SPSAngDistribution::SetBeamSigmaInAngY(    152 void G4SPSAngDistribution::SetBeamSigmaInAngY(G4double r)
155 {                                                 153 {
156   G4AutoLock l(&mutex);                        << 
157   DY = r;                                         154   DY = r;
158 }                                                 155 }
159                                                   156 
160 void G4SPSAngDistribution::                    << 157 void G4SPSAngDistribution::UserDefAngTheta(G4ThreeVector input)
161 SetParticleMomentumDirection(const G4ParticleM << 
162 {                                                 158 {
163   G4AutoLock l(&mutex);                        << 
164   particle_momentum_direction = aMomentumDirec << 
165 }                                              << 
166                                                << 
167 void G4SPSAngDistribution::SetPosDistribution( << 
168 {                                              << 
169   G4AutoLock l(&mutex);                        << 
170   posDist = a;                                 << 
171 }                                              << 
172                                                << 
173 void G4SPSAngDistribution::SetBiasRndm(G4SPSRa << 
174 {                                              << 
175   G4AutoLock l(&mutex);                        << 
176   angRndm = a;                                 << 
177 }                                              << 
178                                                << 
179 void G4SPSAngDistribution::SetVerbosity(G4int  << 
180 {                                              << 
181   G4AutoLock l(&mutex);                        << 
182   verbosityLevel = a;                          << 
183 }                                              << 
184                                                << 
185 void G4SPSAngDistribution::UserDefAngTheta(con << 
186 {                                              << 
187   G4AutoLock l(&mutex);                        << 
188   if(UserDistType == "NULL") UserDistType = "t    159   if(UserDistType == "NULL") UserDistType = "theta";
189   if(UserDistType == "phi") UserDistType = "bo    160   if(UserDistType == "phi") UserDistType = "both";  
190   G4double thi, val;                              161   G4double thi, val;
191   thi = input.x();                                162   thi = input.x();
192   val = input.y();                                163   val = input.y();
193   if(verbosityLevel >= 1) G4cout << "In UserDe << 164   if(verbosityLevel >= 1)
                                                   >> 165     G4cout << "In UserDefAngTheta" << G4endl;
194   UDefThetaH.InsertValues(thi, val);              166   UDefThetaH.InsertValues(thi, val);
195 }                                                 167 }
196                                                   168 
197 G4String G4SPSAngDistribution::GetDistType()   << 169 void G4SPSAngDistribution::UserDefAngPhi(G4ThreeVector input)
198 {                                              << 
199   G4AutoLock l(&mutex);                        << 
200   return AngDistType;                          << 
201 }                                              << 
202                                                << 
203 G4double G4SPSAngDistribution::GetMinTheta()   << 
204 {                                                 170 {
205   G4AutoLock l(&mutex);                        << 
206   return MinTheta;                             << 
207 }                                              << 
208                                                << 
209 G4double G4SPSAngDistribution::GetMaxTheta()   << 
210 {                                              << 
211   G4AutoLock l(&mutex);                        << 
212   return MaxTheta;                             << 
213 }                                              << 
214                                                << 
215 G4double G4SPSAngDistribution::GetMinPhi()     << 
216 {                                              << 
217   G4AutoLock l(&mutex);                        << 
218   return MinPhi;                               << 
219 }                                              << 
220                                                << 
221 G4double G4SPSAngDistribution::GetMaxPhi()     << 
222 {                                              << 
223   G4AutoLock l(&mutex);                        << 
224   return MaxPhi;                               << 
225 }                                              << 
226                                                << 
227 G4ThreeVector G4SPSAngDistribution::GetDirecti << 
228 {                                              << 
229   G4AutoLock l(&mutex);                        << 
230   return particle_momentum_direction;          << 
231 }                                              << 
232                                                << 
233 void G4SPSAngDistribution::UserDefAngPhi(const << 
234 {                                              << 
235   G4AutoLock l(&mutex);                        << 
236   if(UserDistType == "NULL") UserDistType = "p    171   if(UserDistType == "NULL") UserDistType = "phi";
237   if(UserDistType == "theta") UserDistType = "    172   if(UserDistType == "theta") UserDistType = "both";  
238   G4double phhi, val;                             173   G4double phhi, val;
239   phhi = input.x();                               174   phhi = input.x();
240   val = input.y();                                175   val = input.y();
241   if(verbosityLevel >= 1) G4cout << "In UserDe << 176   if(verbosityLevel >= 1)
                                                   >> 177     G4cout << "In UserDefAngPhi" << G4endl;
242   UDefPhiH.InsertValues(phhi, val);               178   UDefPhiH.InsertValues(phhi, val); 
243 }                                                 179 }
244                                                   180 
245 void G4SPSAngDistribution::SetFocusPoint(const << 181 void G4SPSAngDistribution::SetFocusPoint(G4ThreeVector input)
246 {                                                 182 {
247   G4AutoLock l(&mutex);                        << 
248   FocusPoint = input;                             183   FocusPoint = input;
249 }                                                 184 }
250                                                   185 
251 void G4SPSAngDistribution::SetUserWRTSurface(G    186 void G4SPSAngDistribution::SetUserWRTSurface(G4bool wrtSurf)
252 {                                                 187 {
253   G4AutoLock l(&mutex);                        << 188   // This is only applied in user mode?
254                                                << 
255   // if UserWRTSurface = true then the user wa    189   // if UserWRTSurface = true then the user wants momenta with respect
256   // to the surface normals.                      190   // to the surface normals.
257   // When doing this theta has to be 0-90 only    191   // When doing this theta has to be 0-90 only otherwise there will be
258   // errors, which currently are flagged anywh    192   // errors, which currently are flagged anywhere.
259   //                                           << 
260   UserWRTSurface = wrtSurf;                       193   UserWRTSurface = wrtSurf;
261 }                                                 194 }
262                                                   195 
263 void G4SPSAngDistribution::SetUseUserAngAxis(G    196 void G4SPSAngDistribution::SetUseUserAngAxis(G4bool userang)
264 {                                                 197 {
265   G4AutoLock l(&mutex);                        << 
266                                                << 
267   // if UserAngRef = true  the angular distrib    198   // if UserAngRef = true  the angular distribution is defined wrt 
268   // the user defined coordinates              << 199   // the user defined co-ordinates
269   //                                           << 
270   UserAngRef = userang;                           200   UserAngRef = userang;
271 }                                                 201 }
272                                                   202 
273 void G4SPSAngDistribution::GenerateBeamFlux(G4 << 203 void G4SPSAngDistribution::GenerateBeamFlux()
274 {                                                 204 {
275   G4double theta, phi;                            205   G4double theta, phi;
276   G4double px, py, pz;                            206   G4double px, py, pz;
277   if (AngDistType == "beam1d")                 << 207   if (AngDistType == "beam1d") 
278   {                                            << 
279     theta = G4RandGauss::shoot(0.0,DR);        << 
280     phi = twopi * G4UniformRand();             << 
281   }                                            << 
282   else                                         << 
283   {                                            << 
284     px = G4RandGauss::shoot(0.0,DX);           << 
285     py = G4RandGauss::shoot(0.0,DY);           << 
286     theta = std::sqrt (px*px + py*py);         << 
287     if (theta != 0.)                           << 
288     {                                             208     { 
289       phi = std::acos(px/theta);               << 209       theta = G4RandGauss::shoot(0.0,DR);
290       if ( py < 0.) phi = -phi;                << 210       phi = twopi * G4UniformRand();
291     }                                             211     }
292     else                                       << 212   else 
293     {                                          << 213     { 
294       phi = 0.0;                               << 214       px = G4RandGauss::shoot(0.0,DX);
                                                   >> 215       py = G4RandGauss::shoot(0.0,DY);
                                                   >> 216       theta = std::sqrt (px*px + py*py);
                                                   >> 217       if (theta != 0.) { 
                                                   >> 218   phi = std::acos(px/theta);
                                                   >> 219   if ( py < 0.) phi = -phi;
                                                   >> 220       }
                                                   >> 221       else
                                                   >> 222   { 
                                                   >> 223     phi = 0.0;
                                                   >> 224   }
295     }                                             225     }
296   }                                            << 
297   px = -std::sin(theta) * std::cos(phi);          226   px = -std::sin(theta) * std::cos(phi);
298   py = -std::sin(theta) * std::sin(phi);          227   py = -std::sin(theta) * std::sin(phi);
299   pz = -std::cos(theta);                          228   pz = -std::cos(theta);
300   G4double finx, finy,  finz;                  << 229   G4double finx, finy,  finz ;
301   finx=px, finy=py, finz=pz;                   << 230   finx = px, finy =py, finz =pz;
302   if (UserAngRef)                              << 231   if (UserAngRef){
303   {                                            << 
304     // Apply Angular Rotation Matrix              232     // Apply Angular Rotation Matrix
305     // x * AngRef1, y * AngRef2 and z * AngRef    233     // x * AngRef1, y * AngRef2 and z * AngRef3
306     finx = (px * AngRef1.x()) + (py * AngRef2.    234     finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
307     finy = (px * AngRef1.y()) + (py * AngRef2.    235     finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
308     finz = (px * AngRef1.z()) + (py * AngRef2.    236     finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
309     G4double ResMag = std::sqrt((finx*finx) +     237     G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
310     finx = finx/ResMag;                           238     finx = finx/ResMag;
311     finy = finy/ResMag;                           239     finy = finy/ResMag;
312     finz = finz/ResMag;                           240     finz = finz/ResMag;
313   }                                               241   }
314   mom.setX(finx);                              << 242   particle_momentum_direction.setX(finx);
315   mom.setY(finy);                              << 243   particle_momentum_direction.setY(finy);
316   mom.setZ(finz);                              << 244   particle_momentum_direction.setZ(finz);
317                                                << 
318   // particle_momentum_direction now holds uni << 
319                                                   245 
                                                   >> 246   // particle_momentum_direction now holds unit momentum vector.
320   if(verbosityLevel >= 1)                         247   if(verbosityLevel >= 1)
321   {                                            << 248     G4cout << "Generating beam vector: " << particle_momentum_direction << G4endl;
322     G4cout << "Generating beam vector: " << mo << 
323   }                                            << 
324 }                                                 249 }
325                                                   250 
326 void G4SPSAngDistribution::GenerateFocusedFlux << 251 void G4SPSAngDistribution::GenerateFocusedFlux()
327 {                                                 252 {
328   mom = (FocusPoint - posDist->GetParticlePos( << 253   particle_momentum_direction = (FocusPoint - posDist->particle_position).unit();
329                                                << 254   // 
330   // particle_momentum_direction now holds uni    255   // particle_momentum_direction now holds unit momentum vector.
331                                                << 
332   if(verbosityLevel >= 1)                         256   if(verbosityLevel >= 1)
333   {                                            << 257     G4cout << "Generating focused vector: " << particle_momentum_direction << G4endl;
334     G4cout << "Generating focused vector: " << << 
335   }                                            << 
336 }                                                 258 }
337                                                   259 
338 void G4SPSAngDistribution::GenerateIsotropicFl << 260 void G4SPSAngDistribution::GenerateIsotropicFlux()
339 {                                                 261 {
340   // generates isotropic flux.                    262   // generates isotropic flux.
341   // No vectors are needed.                       263   // No vectors are needed.
342                                                << 
343   G4double rndm, rndm2;                           264   G4double rndm, rndm2;
344   G4double px, py, pz;                            265   G4double px, py, pz;
345                                                   266 
                                                   >> 267   //
346   G4double sintheta, sinphi,costheta,cosphi;      268   G4double sintheta, sinphi,costheta,cosphi;
347   rndm = angRndm->GenRandTheta();                 269   rndm = angRndm->GenRandTheta();
348   costheta = std::cos(MinTheta) - rndm * (std: << 270   costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
349                                         - std: << 
350   sintheta = std::sqrt(1. - costheta*costheta)    271   sintheta = std::sqrt(1. - costheta*costheta);
351                                                   272   
352   rndm2 = angRndm->GenRandPhi();                  273   rndm2 = angRndm->GenRandPhi();
353   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;       274   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 
354   sinphi = std::sin(Phi);                         275   sinphi = std::sin(Phi);
355   cosphi = std::cos(Phi);                         276   cosphi = std::cos(Phi);
356                                                   277 
357   px = -sintheta * cosphi;                        278   px = -sintheta * cosphi;
358   py = -sintheta * sinphi;                        279   py = -sintheta * sinphi;
359   pz = -costheta;                                 280   pz = -costheta;
360                                                   281 
361   // For volume and point source use mother or << 282   // for volume and ponit source use mother or user defined co-ordinates
362   // for plane and surface source user surface << 283   // for plane and surface source user surface-normal or userdefined co-ordinates
363   // coordinates                               << 
364   //                                              284   //
365   G4double finx, finy, finz;                      285   G4double finx, finy, finz;
366   if (posDist->GetSourcePosType() == "Point"   << 286   if (posDist->SourcePosType == "Point" || posDist->SourcePosType == "Volume") {
367    || posDist->GetSourcePosType() == "Volume") << 287     if (UserAngRef){
368   {                                            << 
369     if (UserAngRef)                            << 
370     {                                          << 
371       // Apply Rotation Matrix                    288       // Apply Rotation Matrix
372       // x * AngRef1, y * AngRef2 and z * AngR    289       // x * AngRef1, y * AngRef2 and z * AngRef3
373       finx = (px * AngRef1.x()) + (py * AngRef    290       finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
374       finy = (px * AngRef1.y()) + (py * AngRef    291       finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
375       finz = (px * AngRef1.z()) + (py * AngRef    292       finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
376     }                                          << 293     } else {
377     else                                       << 
378     {                                          << 
379       finx = px;                                  294       finx = px;
380       finy = py;                                  295       finy = py;
381       finz = pz;                                  296       finz = pz;
382     }                                             297     }
383   }                                            << 298   } else {    // for plane and surface source   
384   else                                         << 299     if (UserAngRef){
385   {    // for plane and surface source         << 
386     if (UserAngRef)                            << 
387     {                                          << 
388       // Apply Rotation Matrix                    300       // Apply Rotation Matrix
389       // x * AngRef1, y * AngRef2 and z * AngR    301       // x * AngRef1, y * AngRef2 and z * AngRef3
390       finx = (px * AngRef1.x()) + (py * AngRef    302       finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
391       finy = (px * AngRef1.y()) + (py * AngRef    303       finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
392       finz = (px * AngRef1.z()) + (py * AngRef    304       finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
393     }                                          << 305     } else {
394     else                                       << 306       finx = (px*posDist->SideRefVec1.x()) + (py*posDist->SideRefVec2.x()) + (pz*posDist->SideRefVec3.x());
395     {                                          << 307       finy = (px*posDist->SideRefVec1.y()) + (py*posDist->SideRefVec2.y()) + (pz*posDist->SideRefVec3.y());
396       finx = (px*posDist->GetSideRefVec1().x() << 308       finz = (px*posDist->SideRefVec1.z()) + (py*posDist->SideRefVec2.z()) + (pz*posDist->SideRefVec3.z());
397            + (py*posDist->GetSideRefVec2().x() << 
398            + (pz*posDist->GetSideRefVec3().x() << 
399       finy = (px*posDist->GetSideRefVec1().y() << 
400            + (py*posDist->GetSideRefVec2().y() << 
401            + (pz*posDist->GetSideRefVec3().y() << 
402       finz = (px*posDist->GetSideRefVec1().z() << 
403            + (py*posDist->GetSideRefVec2().z() << 
404            + (pz*posDist->GetSideRefVec3().z() << 
405     }                                             309     }
406   }                                               310   }
407   G4double ResMag = std::sqrt((finx*finx) + (f    311   G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
408   finx = finx/ResMag;                             312   finx = finx/ResMag;
409   finy = finy/ResMag;                             313   finy = finy/ResMag;
410   finz = finz/ResMag;                             314   finz = finz/ResMag;
411                                                   315 
412   mom.setX(finx);                              << 316   particle_momentum_direction.setX(finx);
413   mom.setY(finy);                              << 317   particle_momentum_direction.setY(finy);
414   mom.setZ(finz);                              << 318   particle_momentum_direction.setZ(finz);
415                                                   319 
416   // particle_momentum_direction now holds uni    320   // particle_momentum_direction now holds unit momentum vector.
417                                                << 
418   if(verbosityLevel >= 1)                         321   if(verbosityLevel >= 1)
419   {                                            << 322     G4cout << "Generating isotropic vector: " << particle_momentum_direction << G4endl;
420     G4cout << "Generating isotropic vector: "  << 
421   }                                            << 
422 }                                                 323 }
423                                                   324 
424 void G4SPSAngDistribution::GenerateCosineLawFl << 325 void G4SPSAngDistribution::GenerateCosineLawFlux()
425 {                                                 326 {
426   // Method to generate flux distributed with     327   // Method to generate flux distributed with a cosine law
427                                                << 
428   G4double px, py, pz;                            328   G4double px, py, pz;
429   G4double rndm, rndm2;                           329   G4double rndm, rndm2;
430                                                << 330   //
431   G4double sintheta, sinphi,costheta,cosphi;      331   G4double sintheta, sinphi,costheta,cosphi;
432   rndm = angRndm->GenRandTheta();                 332   rndm = angRndm->GenRandTheta();
433   sintheta = std::sqrt( rndm * (std::sin(MaxTh << 333   sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta) - std::sin(MinTheta)*std::sin(MinTheta) ) 
434                               - std::sin(MinTh << 334        +std::sin(MinTheta)*std::sin(MinTheta) );
435                       + std::sin(MinTheta)*std << 
436   costheta = std::sqrt(1. -sintheta*sintheta);    335   costheta = std::sqrt(1. -sintheta*sintheta);
437                                                   336   
438   rndm2 = angRndm->GenRandPhi();                  337   rndm2 = angRndm->GenRandPhi();
439   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;       338   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 
440   sinphi = std::sin(Phi);                         339   sinphi = std::sin(Phi);
441   cosphi = std::cos(Phi);                         340   cosphi = std::cos(Phi);
442                                                   341 
443   px = -sintheta * cosphi;                        342   px = -sintheta * cosphi;
444   py = -sintheta * sinphi;                        343   py = -sintheta * sinphi;
445   pz = -costheta;                                 344   pz = -costheta;
446                                                   345 
447   // for volume and point source use mother or << 346   // for volume and ponit source use mother or user defined co-ordinates
448   // for plane and surface source user surface << 347   // for plane and surface source user surface-normal or userdefined co-ordinates
449   // coordinates                               << 
450   //                                              348   //
451   G4double finx, finy, finz;                      349   G4double finx, finy, finz;
452   if (posDist->GetSourcePosType() == "Point"   << 350   if (posDist->SourcePosType == "Point" || posDist->SourcePosType == "Volume") {
453    || posDist->GetSourcePosType() == "Volume") << 351     if (UserAngRef){
454   {                                            << 
455     if (UserAngRef)                            << 
456     {                                          << 
457       // Apply Rotation Matrix                    352       // Apply Rotation Matrix
458       finx = (px * AngRef1.x()) + (py * AngRef    353       finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
459       finy = (px * AngRef1.y()) + (py * AngRef    354       finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
460       finz = (px * AngRef1.z()) + (py * AngRef    355       finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
461     }                                          << 356     } else {
462     else                                       << 
463     {                                          << 
464       finx = px;                                  357       finx = px;
465       finy = py;                                  358       finy = py;
466       finz = pz;                                  359       finz = pz;
467     }                                             360     }
468   }                                            << 361   } else {    // for plane and surface source   
469   else                                         << 362     if (UserAngRef){
470   {    // for plane and surface source         << 
471     if (UserAngRef)                            << 
472     {                                          << 
473       // Apply Rotation Matrix                    363       // Apply Rotation Matrix
474       finx = (px * AngRef1.x()) + (py * AngRef    364       finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
475       finy = (px * AngRef1.y()) + (py * AngRef    365       finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
476       finz = (px * AngRef1.z()) + (py * AngRef    366       finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
477     }                                          << 367     } else {
478     else                                       << 368       finx = (px*posDist->SideRefVec1.x()) + (py*posDist->SideRefVec2.x()) + (pz*posDist->SideRefVec3.x());
479     {                                          << 369       finy = (px*posDist->SideRefVec1.y()) + (py*posDist->SideRefVec2.y()) + (pz*posDist->SideRefVec3.y());
480       finx = (px*posDist->GetSideRefVec1().x() << 370       finz = (px*posDist->SideRefVec1.z()) + (py*posDist->SideRefVec2.z()) + (pz*posDist->SideRefVec3.z());
481            + (py*posDist->GetSideRefVec2().x() << 
482            + (pz*posDist->GetSideRefVec3().x() << 
483       finy = (px*posDist->GetSideRefVec1().y() << 
484            + (py*posDist->GetSideRefVec2().y() << 
485            + (pz*posDist->GetSideRefVec3().y() << 
486       finz = (px*posDist->GetSideRefVec1().z() << 
487            + (py*posDist->GetSideRefVec2().z() << 
488            + (pz*posDist->GetSideRefVec3().z() << 
489     }                                             371     }
490   }                                               372   }
491   G4double ResMag = std::sqrt((finx*finx) + (f    373   G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
492   finx = finx/ResMag;                             374   finx = finx/ResMag;
493   finy = finy/ResMag;                             375   finy = finy/ResMag;
494   finz = finz/ResMag;                             376   finz = finz/ResMag;
495                                                   377 
496   mom.setX(finx);                              << 378   particle_momentum_direction.setX(finx);
497   mom.setY(finy);                              << 379   particle_momentum_direction.setY(finy);
498   mom.setZ(finz);                              << 380   particle_momentum_direction.setZ(finz);
499                                                   381 
500   // particle_momentum_direction now contains     382   // particle_momentum_direction now contains unit momentum vector.
501                                                << 
502   if(verbosityLevel >= 1)                         383   if(verbosityLevel >= 1)
503   {                                            << 384     {
504     G4cout << "Resultant cosine-law unit momen << 385       G4cout << "Resultant cosine-law unit momentum vector " << particle_momentum_direction << G4endl;
505   }                                            << 386     }
506 }                                                 387 }
507                                                   388 
508 void G4SPSAngDistribution::GeneratePlanarFlux( << 389 void G4SPSAngDistribution::GeneratePlanarFlux()
509 {                                                 390 {
510   // particle_momentum_direction now contains     391   // particle_momentum_direction now contains unit momentum vector.
511   // nothing need be done here as the m-direct    392   // nothing need be done here as the m-directions have been set directly
512   // under this option                            393   // under this option
513                                                << 
514   if(verbosityLevel >= 1)                         394   if(verbosityLevel >= 1)
515   {                                            << 395     {
516     G4cout << "Resultant Planar wave  momentum << 396       G4cout << "Resultant Planar wave  momentum vector " << particle_momentum_direction << G4endl;
517   }                                            << 397     }
518 }                                                 398 }
519                                                   399 
520 void G4SPSAngDistribution::GenerateUserDefFlux << 400 void G4SPSAngDistribution::GenerateUserDefFlux()
521 {                                                 401 {
522   G4double rndm, px, py, pz, pmag;                402   G4double rndm, px, py, pz, pmag;
523                                                   403 
524   if(UserDistType == "NULL")                      404   if(UserDistType == "NULL")
525   {                                            << 
526     G4cout << "Error: UserDistType undefined"     405     G4cout << "Error: UserDistType undefined" << G4endl;
527   }                                            << 406   else if(UserDistType == "theta") {
528   else if(UserDistType == "theta")             << 
529   {                                            << 
530     Theta = 10.;                                  407     Theta = 10.;
531     while(Theta > MaxTheta || Theta < MinTheta    408     while(Theta > MaxTheta || Theta < MinTheta)
532     {                                          << 
533       Theta = GenerateUserDefTheta();             409       Theta = GenerateUserDefTheta();
534     }                                          << 
535     Phi = 10.;                                    410     Phi = 10.;
536     while(Phi > MaxPhi || Phi < MinPhi)        << 411     while(Phi > MaxPhi || Phi < MinPhi) {
537     {                                          << 
538       rndm = angRndm->GenRandPhi();               412       rndm = angRndm->GenRandPhi();
539       Phi = twopi * rndm;                         413       Phi = twopi * rndm;
540     }                                             414     }
541   }                                               415   }
542   else if(UserDistType == "phi")               << 416   else if(UserDistType == "phi") {
543   {                                            << 
544     Theta = 10.;                                  417     Theta = 10.;
545     while(Theta > MaxTheta || Theta < MinTheta    418     while(Theta > MaxTheta || Theta < MinTheta)
546     {                                          << 419       {
547       rndm = angRndm->GenRandTheta();          << 420   rndm = angRndm->GenRandTheta();
548       Theta = std::acos(1. - (2. * rndm));     << 421   Theta = std::acos(1. - (2. * rndm));
549     }                                          << 422       }
550     Phi = 10.;                                    423     Phi = 10.;
551     while(Phi > MaxPhi || Phi < MinPhi)           424     while(Phi > MaxPhi || Phi < MinPhi)
552     {                                          << 
553       Phi = GenerateUserDefPhi();                 425       Phi = GenerateUserDefPhi();
554     }                                          << 
555   }                                               426   }
556   else if(UserDistType == "both")                 427   else if(UserDistType == "both")
557   {                                            << 
558     Theta = 10.;                               << 
559     while(Theta > MaxTheta || Theta < MinTheta << 
560     {                                             428     {
561       Theta = GenerateUserDefTheta();          << 429       Theta = 10.;
562     }                                          << 430       while(Theta > MaxTheta || Theta < MinTheta)      
563     Phi = 10.;                                 << 431   Theta = GenerateUserDefTheta();
564     while(Phi > MaxPhi || Phi < MinPhi)        << 432       Phi = 10.;
565     {                                          << 433       while(Phi > MaxPhi || Phi < MinPhi)
566       Phi = GenerateUserDefPhi();              << 434   Phi = GenerateUserDefPhi();
567     }                                             435     }
568   }                                            << 
569   px = -std::sin(Theta) * std::cos(Phi);          436   px = -std::sin(Theta) * std::cos(Phi);
570   py = -std::sin(Theta) * std::sin(Phi);          437   py = -std::sin(Theta) * std::sin(Phi);
571   pz = -std::cos(Theta);                          438   pz = -std::cos(Theta);
572                                                   439 
573   pmag = std::sqrt((px*px) + (py*py) + (pz*pz)    440   pmag = std::sqrt((px*px) + (py*py) + (pz*pz));
574                                                   441 
575   if(!UserWRTSurface)                          << 442   if(!UserWRTSurface) {
576   {                                            << 
577     G4double finx, finy, finz;                    443     G4double finx, finy, finz;      
578     if (UserAngRef)                            << 444     if (UserAngRef) {
579     {                                          << 
580       // Apply Rotation Matrix                    445       // Apply Rotation Matrix
581       // x * AngRef1, y * AngRef2 and z * AngR    446       // x * AngRef1, y * AngRef2 and z * AngRef3
582       finx = (px * AngRef1.x()) + (py * AngRef    447       finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
583       finy = (px * AngRef1.y()) + (py * AngRef    448       finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
584       finz = (px * AngRef1.z()) + (py * AngRef    449       finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
585     }                                          << 450     } else {  // use mother co-ordinates
586     else    // use mother coordinates          << 
587     {                                          << 
588       finx = px;                                  451       finx = px;
589       finy = py;                                  452       finy = py;
590       finz = pz;                                  453       finz = pz;
591     }                                             454     }
592     G4double ResMag = std::sqrt((finx*finx) +     455     G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
593     finx = finx/ResMag;                           456     finx = finx/ResMag;
594     finy = finy/ResMag;                           457     finy = finy/ResMag;
595     finz = finz/ResMag;                           458     finz = finz/ResMag;
596                                                   459     
597     mom.setX(finx);                            << 460     particle_momentum_direction.setX(finx);
598     mom.setY(finy);                            << 461     particle_momentum_direction.setY(finy);
599     mom.setZ(finz);                            << 462     particle_momentum_direction.setZ(finz);     
600   }                                               463   } 
601   else    // UserWRTSurface = true             << 464   else  {  // UserWRTSurface = true
602   {                                            << 
603     G4double pxh = px/pmag;                       465     G4double pxh = px/pmag;
604     G4double pyh = py/pmag;                       466     G4double pyh = py/pmag;
605     G4double pzh = pz/pmag;                       467     G4double pzh = pz/pmag;
606     if(verbosityLevel > 1)                     << 468     if(verbosityLevel > 1) {
607     {                                          << 469       G4cout <<"SideRefVecs " <<posDist->SideRefVec1<<posDist->SideRefVec2<<posDist->SideRefVec3<<G4endl;
608       G4cout << "SideRefVecs " << posDist->Get << 470       G4cout <<"Raw Unit vector "<<pxh<<","<<pyh<<","<<pzh<<G4endl;
609              << posDist->GetSideRefVec2() << p << 471     }
610              << G4endl;                        << 472     G4double resultx = (pxh*posDist->SideRefVec1.x()) + (pyh*posDist->SideRefVec2.x()) + 
611       G4cout << "Raw Unit vector " << pxh      << 473       (pzh*posDist->SideRefVec3.x());
612              << "," << pyh << "," << pzh << G4 << 
613     }                                          << 
614     G4double resultx = (pxh*posDist->GetSideRe << 
615                      + (pyh*posDist->GetSideRe << 
616                      + (pzh*posDist->GetSideRe << 
617                                                   474     
618     G4double resulty = (pxh*posDist->GetSideRe << 475     G4double resulty = (pxh*posDist->SideRefVec1.y()) + (pyh*posDist->SideRefVec2.y()) + 
619                      + (pyh*posDist->GetSideRe << 476       (pzh*posDist->SideRefVec3.y());
620                      + (pzh*posDist->GetSideRe << 
621                                                   477     
622     G4double resultz = (pxh*posDist->GetSideRe << 478     G4double resultz = (pxh*posDist->SideRefVec1.z()) + (pyh*posDist->SideRefVec2.z()) + 
623                      + (pyh*posDist->GetSideRe << 479       (pzh*posDist->SideRefVec3.z());
624                      + (pzh*posDist->GetSideRe << 
625                                                   480     
626     G4double ResMag = std::sqrt((resultx*resul << 481     G4double ResMag = std::sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz));
627                               + (resulty*resul << 
628                               + (resultz*resul << 
629     resultx = resultx/ResMag;                     482     resultx = resultx/ResMag;
630     resulty = resulty/ResMag;                     483     resulty = resulty/ResMag;
631     resultz = resultz/ResMag;                     484     resultz = resultz/ResMag;
632                                                   485     
633     mom.setX(resultx);                         << 486     particle_momentum_direction.setX(resultx);
634     mom.setY(resulty);                         << 487     particle_momentum_direction.setY(resulty);
635     mom.setZ(resultz);                         << 488     particle_momentum_direction.setZ(resultz);
636   }                                               489   }
637                                                   490   
638   // particle_momentum_direction now contains     491   // particle_momentum_direction now contains unit momentum vector.
639                                                << 
640   if(verbosityLevel > 0 )                         492   if(verbosityLevel > 0 )
641   {                                            << 493     {
642     G4cout << "Final User Defined momentum vec << 494       G4cout << "Final User Defined momentum vector " << particle_momentum_direction << G4endl;
643            << particle_momentum_direction << G << 495     }
644   }                                            << 
645 }                                                 496 }
646                                                   497 
647 G4double G4SPSAngDistribution::GenerateUserDef    498 G4double G4SPSAngDistribution::GenerateUserDefTheta()
648 {                                                 499 {
649   // Create cumulative histogram if not alread << 500   // Create cumulative histogram if not already done so. Then use RandFlat
650   // Then use RandFlat::shoot to generate the  << 501   //::shoot to generate the output Theta value.
651                                                << 
652   if(UserDistType == "NULL" || UserDistType ==    502   if(UserDistType == "NULL" || UserDistType == "phi")
653   {                                            << 
654     // No user defined theta distribution      << 
655     G4cout << "Error ***********************"  << 
656     G4cout << "UserDistType = " << UserDistTyp << 
657     return (0.);                               << 
658   }                                            << 
659                                                << 
660   // UserDistType = theta or both and so a the << 
661   // is defined. This should be integrated if  << 
662   G4AutoLock l(&mutex);                        << 
663   if(!IPDFThetaExist)                          << 
664   {                                            << 
665     // IPDF has not been created, so create it << 
666     //                                         << 
667     G4double bins[1024],vals[1024], sum;       << 
668     G4int ii;                                  << 
669     auto  maxbin = G4int(UDefThetaH.GetVectorL << 
670     bins[0] = UDefThetaH.GetLowEdgeEnergy(std: << 
671     vals[0] = UDefThetaH(std::size_t(0));      << 
672     sum = vals[0];                             << 
673     for(ii=1; ii<maxbin; ++ii)                 << 
674     {                                             503     {
675       bins[ii] = UDefThetaH.GetLowEdgeEnergy(s << 504       // No user defined theta distribution
676       vals[ii] = UDefThetaH(std::size_t(ii)) + << 505       G4cout << "Error ***********************" << G4endl;
677       sum = sum + UDefThetaH(std::size_t(ii)); << 506       G4cout << "UserDistType = " << UserDistType << G4endl;
                                                   >> 507       return (0.);
678     }                                             508     }
679     for(ii=0; ii<maxbin; ++ii)                 << 509   else
680     {                                             510     {
681       vals[ii] = vals[ii]/sum;                 << 511       // UserDistType = theta or both and so a theta distribution
682       IPDFThetaH.InsertValues(bins[ii], vals[i << 512       // is defined. This should be integrated if not already done.
                                                   >> 513       if(IPDFThetaExist == false)
                                                   >> 514   {
                                                   >> 515     // IPDF has not been created, so create it
                                                   >> 516     G4double bins[1024],vals[1024], sum;
                                                   >> 517     G4int ii;
                                                   >> 518     G4int maxbin = G4int(UDefThetaH.GetVectorLength());
                                                   >> 519     bins[0] = UDefThetaH.GetLowEdgeEnergy(size_t(0));
                                                   >> 520     vals[0] = UDefThetaH(size_t(0));
                                                   >> 521     sum = vals[0];
                                                   >> 522     for(ii=1;ii<maxbin;ii++)
                                                   >> 523       {
                                                   >> 524         bins[ii] = UDefThetaH.GetLowEdgeEnergy(size_t(ii));
                                                   >> 525         vals[ii] = UDefThetaH(size_t(ii)) + vals[ii-1];
                                                   >> 526         sum = sum + UDefThetaH(size_t(ii));
                                                   >> 527       }
                                                   >> 528     for(ii=0;ii<maxbin;ii++)
                                                   >> 529       {
                                                   >> 530         vals[ii] = vals[ii]/sum;
                                                   >> 531         IPDFThetaH.InsertValues(bins[ii], vals[ii]);
                                                   >> 532       }
                                                   >> 533     // Make IPDFThetaExist = true
                                                   >> 534     IPDFThetaExist = true;
                                                   >> 535   }
                                                   >> 536       // IPDF has been create so carry on
                                                   >> 537       G4double rndm = G4UniformRand();
                                                   >> 538       return(IPDFThetaH.GetEnergy(rndm));
683     }                                             539     }
684       IPDFThetaExist = true;                   << 
685   }                                            << 
686   l.unlock();                                  << 
687                                                << 
688   // IPDF has been created so carry on         << 
689   //                                           << 
690   G4double rndm = G4UniformRand();             << 
691   return IPDFThetaH.GetEnergy(rndm);           << 
692 }                                                 540 }
693                                                   541 
694 G4double G4SPSAngDistribution::GenerateUserDef    542 G4double G4SPSAngDistribution::GenerateUserDefPhi()
695 {                                                 543 {
696   // Create cumulative histogram if not alread << 544   // Create cumulative histogram if not already done so. Then use RandFlat
697   // Then use RandFlat::shoot to generate the  << 545   //::shoot to generate the output Theta value.
698                                                   546 
699   if(UserDistType == "NULL" || UserDistType ==    547   if(UserDistType == "NULL" || UserDistType == "theta")
700   {                                            << 
701     // No user defined phi distribution        << 
702     G4cout << "Error ***********************"  << 
703     G4cout << "UserDistType = " << UserDistTyp << 
704     return(0.);                                << 
705   }                                            << 
706                                                << 
707   // UserDistType = phi or both and so a phi d << 
708   // is defined. This should be integrated if  << 
709   G4AutoLock l(&mutex);                        << 
710   if(!IPDFPhiExist)                            << 
711   {                                            << 
712     // IPDF has not been created, so create it << 
713     //                                         << 
714     G4double bins[1024],vals[1024], sum;       << 
715     G4int ii;                                  << 
716     auto  maxbin = G4int(UDefPhiH.GetVectorLen << 
717     bins[0] = UDefPhiH.GetLowEdgeEnergy(std::s << 
718     vals[0] = UDefPhiH(std::size_t(0));        << 
719     sum = vals[0];                             << 
720     for(ii=1; ii<maxbin; ++ii)                 << 
721     {                                             548     {
722       bins[ii] = UDefPhiH.GetLowEdgeEnergy(std << 549       // No user defined phi distribution
723       vals[ii] = UDefPhiH(std::size_t(ii)) + v << 550       G4cout << "Error ***********************" << G4endl;
724       sum = sum + UDefPhiH(std::size_t(ii));   << 551       G4cout << "UserDistType = " << UserDistType << G4endl;
                                                   >> 552       return(0.);
725     }                                             553     }
726     for(ii=0; ii<maxbin; ++ii)                 << 554   else
727     {                                             555     {
728       vals[ii] = vals[ii]/sum;                 << 556       // UserDistType = phi or both and so a phi distribution
729       IPDFPhiH.InsertValues(bins[ii], vals[ii] << 557       // is defined. This should be integrated if not already done.
                                                   >> 558       if(IPDFPhiExist == false)
                                                   >> 559   {
                                                   >> 560     // IPDF has not been created, so create it
                                                   >> 561     G4double bins[1024],vals[1024], sum;
                                                   >> 562     G4int ii;
                                                   >> 563     G4int maxbin = G4int(UDefPhiH.GetVectorLength());
                                                   >> 564     bins[0] = UDefPhiH.GetLowEdgeEnergy(size_t(0));
                                                   >> 565     vals[0] = UDefPhiH(size_t(0));
                                                   >> 566     sum = vals[0];
                                                   >> 567     for(ii=1;ii<maxbin;ii++)
                                                   >> 568       {
                                                   >> 569         bins[ii] = UDefPhiH.GetLowEdgeEnergy(size_t(ii));
                                                   >> 570         vals[ii] = UDefPhiH(size_t(ii)) + vals[ii-1];
                                                   >> 571         sum = sum + UDefPhiH(size_t(ii));
                                                   >> 572       }
                                                   >> 573 
                                                   >> 574     for(ii=0;ii<maxbin;ii++)
                                                   >> 575       {
                                                   >> 576         vals[ii] = vals[ii]/sum;
                                                   >> 577         IPDFPhiH.InsertValues(bins[ii], vals[ii]);
                                                   >> 578       }
                                                   >> 579     // Make IPDFPhiExist = true
                                                   >> 580     IPDFPhiExist = true;
                                                   >> 581   }
                                                   >> 582       // IPDF has been create so carry on
                                                   >> 583       G4double rndm = G4UniformRand();
                                                   >> 584       return(IPDFPhiH.GetEnergy(rndm));
730     }                                             585     }
731     IPDFPhiExist = true;                       << 
732   }                                            << 
733   l.unlock();                                  << 
734                                                << 
735   // IPDF has been create so carry on          << 
736   //                                           << 
737   G4double rndm = G4UniformRand();             << 
738   return IPDFPhiH.GetEnergy(rndm);             << 
739 }                                                 586 }
740                                                << 587 //
741 void G4SPSAngDistribution::ReSetHist(const G4S << 588 void G4SPSAngDistribution::ReSetHist(G4String atype)
742 {                                                 589 {
743   G4AutoLock l(&mutex);                        << 590   if (atype == "theta") {
744   if (atype == "theta")                        << 
745   {                                            << 
746     UDefThetaH = IPDFThetaH = ZeroPhysVector ;    591     UDefThetaH = IPDFThetaH = ZeroPhysVector ;
747     IPDFThetaExist = false ;                   << 592     IPDFThetaExist = false ;}
748   }                                            << 593   else if (atype == "phi"){    
749   else if (atype == "phi")                     << 
750   {                                            << 
751     UDefPhiH = IPDFPhiH = ZeroPhysVector ;        594     UDefPhiH = IPDFPhiH = ZeroPhysVector ;
752     IPDFPhiExist = false ;                     << 595     IPDFPhiExist = false ;} 
753   }                                            << 596   else {
754   else                                         << 
755   {                                            << 
756     G4cout << "Error, histtype not accepted "     597     G4cout << "Error, histtype not accepted " << G4endl;
757   }                                               598   }
758 }                                                 599 }
759                                                   600 
                                                   >> 601 
760 G4ParticleMomentum G4SPSAngDistribution::Gener    602 G4ParticleMomentum G4SPSAngDistribution::GenerateOne()
761 {                                                 603 {
762   // Local copy for thread safety              << 
763   //                                           << 
764   G4ParticleMomentum localM = particle_momentu << 
765                                                << 
766   // Angular stuff                                604   // Angular stuff
767   //                                           << 
768   if(AngDistType == "iso")                        605   if(AngDistType == "iso")
769     GenerateIsotropicFlux(localM);             << 606     GenerateIsotropicFlux();
770   else if(AngDistType == "cos")                   607   else if(AngDistType == "cos")
771     GenerateCosineLawFlux(localM);             << 608     GenerateCosineLawFlux();
772   else if(AngDistType == "planar")                609   else if(AngDistType == "planar")
773     GeneratePlanarFlux(localM);                << 610     GeneratePlanarFlux();
774   else if(AngDistType == "beam1d" || AngDistTy    611   else if(AngDistType == "beam1d" || AngDistType == "beam2d" )
775     GenerateBeamFlux(localM);                  << 612     GenerateBeamFlux();
776   else if(AngDistType == "user")                  613   else if(AngDistType == "user")
777     GenerateUserDefFlux(localM);               << 614     GenerateUserDefFlux();
778   else if(AngDistType == "focused")               615   else if(AngDistType == "focused")
779     GenerateFocusedFlux(localM);               << 616     GenerateFocusedFlux();
780   else                                            617   else
781     G4cout << "Error: AngDistType has unusual     618     G4cout << "Error: AngDistType has unusual value" << G4endl;
782   return localM;                               << 619   return particle_momentum_direction;
783 }                                                 620 }
                                                   >> 621 
                                                   >> 622 
                                                   >> 623 
                                                   >> 624 
                                                   >> 625 
                                                   >> 626 
                                                   >> 627 
                                                   >> 628 
                                                   >> 629 
784                                                   630