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 7.0)


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