Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SPSAngDistribution.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
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 //
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 //
47 //
48 
49 #include "G4SPSAngDistribution.hh"
50 
51 #include "Randomize.hh"
52 #include "G4PhysicalConstants.hh"
53 
55  : Theta(0.), Phi(0.)
56 {
57  // Angular distribution Variables
58  G4ThreeVector zero;
59  particle_momentum_direction = G4ParticleMomentum(0,0,-1);
60 
61  AngDistType = "planar";
62  AngRef1 = CLHEP::HepXHat;
63  AngRef2 = CLHEP::HepYHat;
64  AngRef3 = CLHEP::HepZHat;
65  MinTheta = 0.;
66  MaxTheta = pi;
67  MinPhi = 0.;
68  MaxPhi = twopi;
69  DR = 0.;
70  DX = 0.;
71  DY = 0.;
72  FocusPoint = G4ThreeVector(0., 0., 0.);
73  UserDistType = "NULL";
74  UserWRTSurface = true;
75  UserAngRef = false;
76  IPDFThetaExist = false;
77  IPDFPhiExist = false;
78  verbosityLevel = 0 ;
79 }
80 
82 {}
83 
84 //
86 {
87  if(atype != "iso" && atype != "cos" && atype != "user" && atype != "planar"
88  && atype != "beam1d" && atype != "beam2d" && atype != "focused")
89  G4cout << "Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user" << G4endl;
90  else
91  AngDistType = atype;
92  if (AngDistType == "cos") MaxTheta = pi/2. ;
93  if (AngDistType == "user") {
94  UDefThetaH = IPDFThetaH = ZeroPhysVector ;
95  IPDFThetaExist = false ;
96  UDefPhiH = IPDFPhiH = ZeroPhysVector ;
97  IPDFPhiExist = false ;
98  }
99 }
100 
102 {
103  if(refname == "angref1")
104  AngRef1 = ref.unit(); // x'
105  else if(refname == "angref2")
106  AngRef2 = ref.unit(); // vector in x'y' plane
107 
108  // User defines x' (AngRef1) and a vector in the x'y'
109  // plane (AngRef2). Then, AngRef1 x AngRef2 = AngRef3
110  // the z' vector. Then, AngRef3 x AngRef1 = AngRef2
111  // which will now be y'.
112 
113  AngRef3 = AngRef1.cross(AngRef2); // z'
114  AngRef2 = AngRef3.cross(AngRef1); // y'
115  UserAngRef = true ;
116  if(verbosityLevel == 2)
117  {
118  G4cout << "Angular distribution rotation axes " << AngRef1 << " " << AngRef2 << " " << AngRef3 << G4endl;
119  }
120 }
121 
123 {
124  MinTheta = mint;
125 }
126 
128 {
129  MinPhi = minp;
130 }
131 
133 {
134  MaxTheta = maxt;
135 }
136 
138 {
139  MaxPhi = maxp;
140 }
141 
143 {
144  DR = r;
145 }
146 
148 {
149  DX = r;
150 }
151 
153 {
154  DY = r;
155 }
156 
158 {
159  if(UserDistType == "NULL") UserDistType = "theta";
160  if(UserDistType == "phi") UserDistType = "both";
161  G4double thi, val;
162  thi = input.x();
163  val = input.y();
164  if(verbosityLevel >= 1)
165  G4cout << "In UserDefAngTheta" << G4endl;
166  UDefThetaH.InsertValues(thi, val);
167 }
168 
170 {
171  if(UserDistType == "NULL") UserDistType = "phi";
172  if(UserDistType == "theta") UserDistType = "both";
173  G4double phhi, val;
174  phhi = input.x();
175  val = input.y();
176  if(verbosityLevel >= 1)
177  G4cout << "In UserDefAngPhi" << G4endl;
178  UDefPhiH.InsertValues(phhi, val);
179 }
180 
182 {
183  FocusPoint = input;
184 }
185 
187 {
188  // This is only applied in user mode?
189  // if UserWRTSurface = true then the user wants momenta with respect
190  // to the surface normals.
191  // When doing this theta has to be 0-90 only otherwise there will be
192  // errors, which currently are flagged anywhere.
193  UserWRTSurface = wrtSurf;
194 }
195 
197 {
198  // if UserAngRef = true the angular distribution is defined wrt
199  // the user defined co-ordinates
200  UserAngRef = userang;
201 }
202 
203 void G4SPSAngDistribution::GenerateBeamFlux()
204 {
205  G4double theta, phi;
206  G4double px, py, pz;
207  if (AngDistType == "beam1d")
208  {
209  theta = G4RandGauss::shoot(0.0,DR);
210  phi = twopi * G4UniformRand();
211  }
212  else
213  {
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  }
225  }
226  px = -std::sin(theta) * std::cos(phi);
227  py = -std::sin(theta) * std::sin(phi);
228  pz = -std::cos(theta);
229  G4double finx, finy, finz ;
230  finx = px, finy =py, finz =pz;
231  if (UserAngRef){
232  // Apply Angular Rotation Matrix
233  // x * AngRef1, y * AngRef2 and z * AngRef3
234  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
235  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
236  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
237  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
238  finx = finx/ResMag;
239  finy = finy/ResMag;
240  finz = finz/ResMag;
241  }
242  particle_momentum_direction.setX(finx);
243  particle_momentum_direction.setY(finy);
244  particle_momentum_direction.setZ(finz);
245 
246  // particle_momentum_direction now holds unit momentum vector.
247  if(verbosityLevel >= 1)
248  G4cout << "Generating beam vector: " << particle_momentum_direction << G4endl;
249 }
250 
251 void G4SPSAngDistribution::GenerateFocusedFlux()
252 {
253  particle_momentum_direction = (FocusPoint - posDist->particle_position).unit();
254  //
255  // particle_momentum_direction now holds unit momentum vector.
256  if(verbosityLevel >= 1)
257  G4cout << "Generating focused vector: " << particle_momentum_direction << G4endl;
258 }
259 
260 void G4SPSAngDistribution::GenerateIsotropicFlux()
261 {
262  // generates isotropic flux.
263  // No vectors are needed.
264  G4double rndm, rndm2;
265  G4double px, py, pz;
266 
267  //
268  G4double sintheta, sinphi,costheta,cosphi;
269  rndm = angRndm->GenRandTheta();
270  costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
271  sintheta = std::sqrt(1. - costheta*costheta);
272 
273  rndm2 = angRndm->GenRandPhi();
274  Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
275  sinphi = std::sin(Phi);
276  cosphi = std::cos(Phi);
277 
278  px = -sintheta * cosphi;
279  py = -sintheta * sinphi;
280  pz = -costheta;
281 
282  // for volume and ponit source use mother or user defined co-ordinates
283  // for plane and surface source user surface-normal or userdefined co-ordinates
284  //
285  G4double finx, finy, finz;
286  if (posDist->SourcePosType == "Point" || posDist->SourcePosType == "Volume") {
287  if (UserAngRef){
288  // Apply Rotation Matrix
289  // x * AngRef1, y * AngRef2 and z * AngRef3
290  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
291  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
292  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
293  } else {
294  finx = px;
295  finy = py;
296  finz = pz;
297  }
298  } else { // for plane and surface source
299  if (UserAngRef){
300  // Apply Rotation Matrix
301  // x * AngRef1, y * AngRef2 and z * AngRef3
302  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
303  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
304  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
305  } else {
306  finx = (px*posDist->SideRefVec1.x()) + (py*posDist->SideRefVec2.x()) + (pz*posDist->SideRefVec3.x());
307  finy = (px*posDist->SideRefVec1.y()) + (py*posDist->SideRefVec2.y()) + (pz*posDist->SideRefVec3.y());
308  finz = (px*posDist->SideRefVec1.z()) + (py*posDist->SideRefVec2.z()) + (pz*posDist->SideRefVec3.z());
309  }
310  }
311  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
312  finx = finx/ResMag;
313  finy = finy/ResMag;
314  finz = finz/ResMag;
315 
316  particle_momentum_direction.setX(finx);
317  particle_momentum_direction.setY(finy);
318  particle_momentum_direction.setZ(finz);
319 
320  // particle_momentum_direction now holds unit momentum vector.
321  if(verbosityLevel >= 1)
322  G4cout << "Generating isotropic vector: " << particle_momentum_direction << G4endl;
323 }
324 
325 void G4SPSAngDistribution::GenerateCosineLawFlux()
326 {
327  // Method to generate flux distributed with a cosine law
328  G4double px, py, pz;
329  G4double rndm, rndm2;
330  //
331  G4double sintheta, sinphi,costheta,cosphi;
332  rndm = angRndm->GenRandTheta();
333  sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta) - std::sin(MinTheta)*std::sin(MinTheta) )
334  +std::sin(MinTheta)*std::sin(MinTheta) );
335  costheta = std::sqrt(1. -sintheta*sintheta);
336 
337  rndm2 = angRndm->GenRandPhi();
338  Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
339  sinphi = std::sin(Phi);
340  cosphi = std::cos(Phi);
341 
342  px = -sintheta * cosphi;
343  py = -sintheta * sinphi;
344  pz = -costheta;
345 
346  // for volume and ponit source use mother or user defined co-ordinates
347  // for plane and surface source user surface-normal or userdefined co-ordinates
348  //
349  G4double finx, finy, finz;
350  if (posDist->SourcePosType == "Point" || posDist->SourcePosType == "Volume") {
351  if (UserAngRef){
352  // Apply Rotation Matrix
353  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
354  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
355  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
356  } else {
357  finx = px;
358  finy = py;
359  finz = pz;
360  }
361  } else { // for plane and surface source
362  if (UserAngRef){
363  // Apply Rotation Matrix
364  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
365  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
366  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
367  } else {
368  finx = (px*posDist->SideRefVec1.x()) + (py*posDist->SideRefVec2.x()) + (pz*posDist->SideRefVec3.x());
369  finy = (px*posDist->SideRefVec1.y()) + (py*posDist->SideRefVec2.y()) + (pz*posDist->SideRefVec3.y());
370  finz = (px*posDist->SideRefVec1.z()) + (py*posDist->SideRefVec2.z()) + (pz*posDist->SideRefVec3.z());
371  }
372  }
373  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
374  finx = finx/ResMag;
375  finy = finy/ResMag;
376  finz = finz/ResMag;
377 
378  particle_momentum_direction.setX(finx);
379  particle_momentum_direction.setY(finy);
380  particle_momentum_direction.setZ(finz);
381 
382  // particle_momentum_direction now contains unit momentum vector.
383  if(verbosityLevel >= 1)
384  {
385  G4cout << "Resultant cosine-law unit momentum vector " << particle_momentum_direction << G4endl;
386  }
387 }
388 
389 void G4SPSAngDistribution::GeneratePlanarFlux()
390 {
391  // particle_momentum_direction now contains unit momentum vector.
392  // nothing need be done here as the m-directions have been set directly
393  // under this option
394  if(verbosityLevel >= 1)
395  {
396  G4cout << "Resultant Planar wave momentum vector " << particle_momentum_direction << G4endl;
397  }
398 }
399 
400 void G4SPSAngDistribution::GenerateUserDefFlux()
401 {
402  G4double rndm, px, py, pz, pmag;
403 
404  if(UserDistType == "NULL")
405  G4cout << "Error: UserDistType undefined" << G4endl;
406  else if(UserDistType == "theta") {
407  Theta = 10.;
408  while(Theta > MaxTheta || Theta < MinTheta)
409  Theta = GenerateUserDefTheta();
410  Phi = 10.;
411  while(Phi > MaxPhi || Phi < MinPhi) {
412  rndm = angRndm->GenRandPhi();
413  Phi = twopi * rndm;
414  }
415  }
416  else if(UserDistType == "phi") {
417  Theta = 10.;
418  while(Theta > MaxTheta || Theta < MinTheta)
419  {
420  rndm = angRndm->GenRandTheta();
421  Theta = std::acos(1. - (2. * rndm));
422  }
423  Phi = 10.;
424  while(Phi > MaxPhi || Phi < MinPhi)
425  Phi = GenerateUserDefPhi();
426  }
427  else if(UserDistType == "both")
428  {
429  Theta = 10.;
430  while(Theta > MaxTheta || Theta < MinTheta)
431  Theta = GenerateUserDefTheta();
432  Phi = 10.;
433  while(Phi > MaxPhi || Phi < MinPhi)
434  Phi = GenerateUserDefPhi();
435  }
436  px = -std::sin(Theta) * std::cos(Phi);
437  py = -std::sin(Theta) * std::sin(Phi);
438  pz = -std::cos(Theta);
439 
440  pmag = std::sqrt((px*px) + (py*py) + (pz*pz));
441 
442  if(!UserWRTSurface) {
443  G4double finx, finy, finz;
444  if (UserAngRef) {
445  // Apply Rotation Matrix
446  // x * AngRef1, y * AngRef2 and z * AngRef3
447  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
448  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
449  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
450  } else { // use mother co-ordinates
451  finx = px;
452  finy = py;
453  finz = pz;
454  }
455  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
456  finx = finx/ResMag;
457  finy = finy/ResMag;
458  finz = finz/ResMag;
459 
460  particle_momentum_direction.setX(finx);
461  particle_momentum_direction.setY(finy);
462  particle_momentum_direction.setZ(finz);
463  }
464  else { // UserWRTSurface = true
465  G4double pxh = px/pmag;
466  G4double pyh = py/pmag;
467  G4double pzh = pz/pmag;
468  if(verbosityLevel > 1) {
469  G4cout <<"SideRefVecs " <<posDist->SideRefVec1<<posDist->SideRefVec2<<posDist->SideRefVec3<<G4endl;
470  G4cout <<"Raw Unit vector "<<pxh<<","<<pyh<<","<<pzh<<G4endl;
471  }
472  G4double resultx = (pxh*posDist->SideRefVec1.x()) + (pyh*posDist->SideRefVec2.x()) +
473  (pzh*posDist->SideRefVec3.x());
474 
475  G4double resulty = (pxh*posDist->SideRefVec1.y()) + (pyh*posDist->SideRefVec2.y()) +
476  (pzh*posDist->SideRefVec3.y());
477 
478  G4double resultz = (pxh*posDist->SideRefVec1.z()) + (pyh*posDist->SideRefVec2.z()) +
479  (pzh*posDist->SideRefVec3.z());
480 
481  G4double ResMag = std::sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz));
482  resultx = resultx/ResMag;
483  resulty = resulty/ResMag;
484  resultz = resultz/ResMag;
485 
486  particle_momentum_direction.setX(resultx);
487  particle_momentum_direction.setY(resulty);
488  particle_momentum_direction.setZ(resultz);
489  }
490 
491  // particle_momentum_direction now contains unit momentum vector.
492  if(verbosityLevel > 0 )
493  {
494  G4cout << "Final User Defined momentum vector " << particle_momentum_direction << G4endl;
495  }
496 }
497 
498 G4double G4SPSAngDistribution::GenerateUserDefTheta()
499 {
500  // Create cumulative histogram if not already done so. Then use RandFlat
501  //::shoot to generate the output Theta value.
502  if(UserDistType == "NULL" || UserDistType == "phi")
503  {
504  // No user defined theta distribution
505  G4cout << "Error ***********************" << G4endl;
506  G4cout << "UserDistType = " << UserDistType << G4endl;
507  return (0.);
508  }
509  else
510  {
511  // UserDistType = theta or both and so a theta distribution
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));
539  }
540 }
541 
542 G4double G4SPSAngDistribution::GenerateUserDefPhi()
543 {
544  // Create cumulative histogram if not already done so. Then use RandFlat
545  //::shoot to generate the output Theta value.
546 
547  if(UserDistType == "NULL" || UserDistType == "theta")
548  {
549  // No user defined phi distribution
550  G4cout << "Error ***********************" << G4endl;
551  G4cout << "UserDistType = " << UserDistType << G4endl;
552  return(0.);
553  }
554  else
555  {
556  // UserDistType = phi or both and so a phi distribution
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));
585  }
586 }
587 //
589 {
590  if (atype == "theta") {
591  UDefThetaH = IPDFThetaH = ZeroPhysVector ;
592  IPDFThetaExist = false ;}
593  else if (atype == "phi"){
594  UDefPhiH = IPDFPhiH = ZeroPhysVector ;
595  IPDFPhiExist = false ;}
596  else {
597  G4cout << "Error, histtype not accepted " << G4endl;
598  }
599 }
600 
601 
603 {
604  // Angular stuff
605  if(AngDistType == "iso")
606  GenerateIsotropicFlux();
607  else if(AngDistType == "cos")
608  GenerateCosineLawFlux();
609  else if(AngDistType == "planar")
610  GeneratePlanarFlux();
611  else if(AngDistType == "beam1d" || AngDistType == "beam2d" )
612  GenerateBeamFlux();
613  else if(AngDistType == "user")
614  GenerateUserDefFlux();
615  else if(AngDistType == "focused")
616  GenerateFocusedFlux();
617  else
618  G4cout << "Error: AngDistType has unusual value" << G4endl;
619  return particle_momentum_direction;
620 }
621 
622 
623 
624 
625 
626 
627 
628 
629