ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SPSAngDistribution.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4SPSAngDistribution.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
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 
54 
56  : Theta(0.), Phi(0.), posDist(0),angRndm(0)
57 {
58  // Angular distribution Variables
61 
62  AngDistType = "planar";
66  MinTheta = 0.;
67  MaxTheta = pi;
68  MinPhi = 0.;
69  MaxPhi = twopi;
70  DR = 0.;
71  DX = 0.;
72  DY = 0.;
73  FocusPoint = G4ThreeVector(0., 0., 0.);
74  UserDistType = "NULL";
75  UserWRTSurface = true;
76  UserAngRef = false;
77  IPDFThetaExist = false;
78  IPDFPhiExist = false;
79  verbosityLevel = 0 ;
80 
82 }
83 
85 {
87 }
88 
89 //
91 {
92  G4AutoLock l(&mutex);
93  if(atype != "iso" && atype != "cos" && atype != "user" && atype != "planar"
94  && atype != "beam1d" && atype != "beam2d" && atype != "focused")
95  G4cout << "Error, distribution must be iso, cos, planar, beam1d, beam2d, focused or user" << G4endl;
96  else
97  AngDistType = atype;
98  if (AngDistType == "cos") MaxTheta = pi/2. ;
99  if (AngDistType == "user") {
103  IPDFPhiExist = false ;
104  }
105 }
106 
108 {
109  G4AutoLock l(&mutex);
110  if(refname == "angref1")
111  AngRef1 = ref.unit(); // x'
112  else if(refname == "angref2")
113  AngRef2 = ref.unit(); // vector in x'y' plane
114 
115  // User defines x' (AngRef1) and a vector in the x'y'
116  // plane (AngRef2). Then, AngRef1 x AngRef2 = AngRef3
117  // the z' vector. Then, AngRef3 x AngRef1 = AngRef2
118  // which will now be y'.
119 
120  AngRef3 = AngRef1.cross(AngRef2); // z'
121  AngRef2 = AngRef3.cross(AngRef1); // y'
122  UserAngRef = true ;
123  if(verbosityLevel == 2)
124  {
125  G4cout << "Angular distribution rotation axes " << AngRef1 << " " << AngRef2 << " " << AngRef3 << G4endl;
126  }
127 }
128 
130 {
131  G4AutoLock l(&mutex);
132  MinTheta = mint;
133 }
134 
136 {
137  G4AutoLock l(&mutex);
138  MinPhi = minp;
139 }
140 
142 {
143  G4AutoLock l(&mutex);
144  MaxTheta = maxt;
145 }
146 
148 {
149  G4AutoLock l(&mutex);
150  MaxPhi = maxp;
151 }
152 
154 {
155  G4AutoLock l(&mutex);
156  DR = r;
157 }
158 
160 {
161  G4AutoLock l(&mutex);
162  DX = r;
163 }
164 
166 {
167  G4AutoLock l(&mutex);
168  DY = r;
169 }
170 
172 {
173  G4AutoLock l(&mutex);
174  particle_momentum_direction = aMomentumDirection.unit();
175 
176 }
177 
179 {
180  G4AutoLock l(&mutex);
181  posDist = a;
182 }
183 
185 {
186  G4AutoLock l(&mutex);
187  angRndm = a;
188 }
189 
191 {
192  G4AutoLock l(&mutex);
193  verbosityLevel = a;
194 }
195 
196 
198 {
199  G4AutoLock l(&mutex);
200  if(UserDistType == "NULL") UserDistType = "theta";
201  if(UserDistType == "phi") UserDistType = "both";
202  G4double thi, val;
203  thi = input.x();
204  val = input.y();
205  if(verbosityLevel >= 1)
206  G4cout << "In UserDefAngTheta" << G4endl;
207  UDefThetaH.InsertValues(thi, val);
208 }
209 
216 
218 {
219  G4AutoLock l(&mutex);
220  if(UserDistType == "NULL") UserDistType = "phi";
221  if(UserDistType == "theta") UserDistType = "both";
222  G4double phhi, val;
223  phhi = input.x();
224  val = input.y();
225  if(verbosityLevel >= 1)
226  G4cout << "In UserDefAngPhi" << G4endl;
227  UDefPhiH.InsertValues(phhi, val);
228 }
229 
231 {
232  G4AutoLock l(&mutex);
233  FocusPoint = input;
234 }
235 
237 {
238  G4AutoLock l(&mutex);
239  // This is only applied in user mode?
240  // if UserWRTSurface = true then the user wants momenta with respect
241  // to the surface normals.
242  // When doing this theta has to be 0-90 only otherwise there will be
243  // errors, which currently are flagged anywhere.
244  UserWRTSurface = wrtSurf;
245 }
246 
248 {
249  G4AutoLock l(&mutex);
250  // if UserAngRef = true the angular distribution is defined wrt
251  // the user defined co-ordinates
252  UserAngRef = userang;
253 }
254 
256 {
257  G4double theta, phi;
258  G4double px, py, pz;
259  if (AngDistType == "beam1d")
260  {
261  theta = G4RandGauss::shoot(0.0,DR);
262  phi = twopi * G4UniformRand();
263  }
264  else
265  {
266  px = G4RandGauss::shoot(0.0,DX);
267  py = G4RandGauss::shoot(0.0,DY);
268  theta = std::sqrt (px*px + py*py);
269  if (theta != 0.) {
270  phi = std::acos(px/theta);
271  if ( py < 0.) phi = -phi;
272  } else {
273  phi = 0.0;
274  }
275  }
276  px = -std::sin(theta) * std::cos(phi);
277  py = -std::sin(theta) * std::sin(phi);
278  pz = -std::cos(theta);
279  G4double finx, finy, finz ;
280  finx = px, finy =py, finz =pz;
281  if (UserAngRef){
282  // Apply Angular Rotation Matrix
283  // x * AngRef1, y * AngRef2 and z * AngRef3
284  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
285  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
286  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
287  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
288  finx = finx/ResMag;
289  finy = finy/ResMag;
290  finz = finz/ResMag;
291  }
292  mom.setX(finx);
293  mom.setY(finy);
294  mom.setZ(finz);
295 
296  // particle_momentum_direction now holds unit momentum vector.
297  if(verbosityLevel >= 1)
298  G4cout << "Generating beam vector: " << mom << G4endl;
299 }
300 
302 {
303  mom = (FocusPoint - posDist->GetParticlePos()).unit();
304  //
305  // particle_momentum_direction now holds unit momentum vector.
306  if(verbosityLevel >= 1)
307  G4cout << "Generating focused vector: " << mom << G4endl;
308 }
309 
311 {
312  // generates isotropic flux.
313  // No vectors are needed.
314  G4double rndm, rndm2;
315  G4double px, py, pz;
316 
317  //
318  G4double sintheta, sinphi,costheta,cosphi;
319  rndm = angRndm->GenRandTheta();
320  costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
321  sintheta = std::sqrt(1. - costheta*costheta);
322 
323  rndm2 = angRndm->GenRandPhi();
324  Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
325  sinphi = std::sin(Phi);
326  cosphi = std::cos(Phi);
327 
328  px = -sintheta * cosphi;
329  py = -sintheta * sinphi;
330  pz = -costheta;
331 
332  // for volume and ponit source use mother or user defined co-ordinates
333  // for plane and surface source user surface-normal or userdefined co-ordinates
334  //
335  G4double finx, finy, finz;
336  if (posDist->GetSourcePosType() == "Point" || posDist->GetSourcePosType() == "Volume") {
337  if (UserAngRef){
338  // Apply Rotation Matrix
339  // x * AngRef1, y * AngRef2 and z * AngRef3
340  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
341  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
342  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
343  } else {
344  finx = px;
345  finy = py;
346  finz = pz;
347  }
348  } else { // for plane and surface source
349  if (UserAngRef){
350  // Apply Rotation Matrix
351  // x * AngRef1, y * AngRef2 and z * AngRef3
352  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
353  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
354  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
355  } else {
356  finx = (px*posDist->GetSideRefVec1().x()) + (py*posDist->GetSideRefVec2().x()) + (pz*posDist->GetSideRefVec3().x());
357  finy = (px*posDist->GetSideRefVec1().y()) + (py*posDist->GetSideRefVec2().y()) + (pz*posDist->GetSideRefVec3().y());
358  finz = (px*posDist->GetSideRefVec1().z()) + (py*posDist->GetSideRefVec2().z()) + (pz*posDist->GetSideRefVec3().z());
359  }
360  }
361  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
362  finx = finx/ResMag;
363  finy = finy/ResMag;
364  finz = finz/ResMag;
365 
366  mom.setX(finx);
367  mom.setY(finy);
368  mom.setZ(finz);
369 
370  // particle_momentum_direction now holds unit momentum vector.
371  if(verbosityLevel >= 1)
372  G4cout << "Generating isotropic vector: " << mom << G4endl;
373 }
374 
376 {
377  // Method to generate flux distributed with a cosine law
378  G4double px, py, pz;
379  G4double rndm, rndm2;
380  //
381  G4double sintheta, sinphi,costheta,cosphi;
382  rndm = angRndm->GenRandTheta();
383  sintheta = std::sqrt( rndm * (std::sin(MaxTheta)*std::sin(MaxTheta) - std::sin(MinTheta)*std::sin(MinTheta) )
384  +std::sin(MinTheta)*std::sin(MinTheta) );
385  costheta = std::sqrt(1. -sintheta*sintheta);
386 
387  rndm2 = angRndm->GenRandPhi();
388  Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
389  sinphi = std::sin(Phi);
390  cosphi = std::cos(Phi);
391 
392  px = -sintheta * cosphi;
393  py = -sintheta * sinphi;
394  pz = -costheta;
395 
396  // for volume and ponit source use mother or user defined co-ordinates
397  // for plane and surface source user surface-normal or userdefined co-ordinates
398  //
399  G4double finx, finy, finz;
400  if (posDist->GetSourcePosType() == "Point" || posDist->GetSourcePosType() == "Volume") {
401  if (UserAngRef){
402  // Apply Rotation Matrix
403  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
404  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
405  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
406  } else {
407  finx = px;
408  finy = py;
409  finz = pz;
410  }
411  } else { // for plane and surface source
412  if (UserAngRef){
413  // Apply Rotation Matrix
414  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
415  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
416  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
417  } else {
418  finx = (px*posDist->GetSideRefVec1().x()) + (py*posDist->GetSideRefVec2().x()) + (pz*posDist->GetSideRefVec3().x());
419  finy = (px*posDist->GetSideRefVec1().y()) + (py*posDist->GetSideRefVec2().y()) + (pz*posDist->GetSideRefVec3().y());
420  finz = (px*posDist->GetSideRefVec1().z()) + (py*posDist->GetSideRefVec2().z()) + (pz*posDist->GetSideRefVec3().z());
421  }
422  }
423  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
424  finx = finx/ResMag;
425  finy = finy/ResMag;
426  finz = finz/ResMag;
427 
428  mom.setX(finx);
429  mom.setY(finy);
430  mom.setZ(finz);
431 
432  // particle_momentum_direction now contains unit momentum vector.
433  if(verbosityLevel >= 1)
434  {
435  G4cout << "Resultant cosine-law unit momentum vector " << mom << G4endl;
436  }
437 }
438 
440 {
441  // particle_momentum_direction now contains unit momentum vector.
442  // nothing need be done here as the m-directions have been set directly
443  // under this option
444  if(verbosityLevel >= 1)
445  {
446  G4cout << "Resultant Planar wave momentum vector " << mom << G4endl;
447  }
448 }
449 
451 {
452  G4double rndm, px, py, pz, pmag;
453 
454  if(UserDistType == "NULL")
455  G4cout << "Error: UserDistType undefined" << G4endl;
456  else if(UserDistType == "theta") {
457  Theta = 10.;
458  while(Theta > MaxTheta || Theta < MinTheta)
460  Phi = 10.;
461  while(Phi > MaxPhi || Phi < MinPhi) {
462  rndm = angRndm->GenRandPhi();
463  Phi = twopi * rndm;
464  }
465  }
466  else if(UserDistType == "phi") {
467  Theta = 10.;
468  while(Theta > MaxTheta || Theta < MinTheta)
469  {
470  rndm = angRndm->GenRandTheta();
471  Theta = std::acos(1. - (2. * rndm));
472  }
473  Phi = 10.;
474  while(Phi > MaxPhi || Phi < MinPhi)
476  }
477  else if(UserDistType == "both")
478  {
479  Theta = 10.;
480  while(Theta > MaxTheta || Theta < MinTheta)
482  Phi = 10.;
483  while(Phi > MaxPhi || Phi < MinPhi)
485  }
486  px = -std::sin(Theta) * std::cos(Phi);
487  py = -std::sin(Theta) * std::sin(Phi);
488  pz = -std::cos(Theta);
489 
490  pmag = std::sqrt((px*px) + (py*py) + (pz*pz));
491 
492  if(!UserWRTSurface) {
493  G4double finx, finy, finz;
494  if (UserAngRef) {
495  // Apply Rotation Matrix
496  // x * AngRef1, y * AngRef2 and z * AngRef3
497  finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x());
498  finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y());
499  finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z());
500  } else { // use mother co-ordinates
501  finx = px;
502  finy = py;
503  finz = pz;
504  }
505  G4double ResMag = std::sqrt((finx*finx) + (finy*finy) + (finz*finz));
506  finx = finx/ResMag;
507  finy = finy/ResMag;
508  finz = finz/ResMag;
509 
510  mom.setX(finx);
511  mom.setY(finy);
512  mom.setZ(finz);
513  }
514  else { // UserWRTSurface = true
515  G4double pxh = px/pmag;
516  G4double pyh = py/pmag;
517  G4double pzh = pz/pmag;
518  if(verbosityLevel > 1) {
520  G4cout <<"Raw Unit vector "<<pxh<<","<<pyh<<","<<pzh<<G4endl;
521  }
522  G4double resultx = (pxh*posDist->GetSideRefVec1().x()) + (pyh*posDist->GetSideRefVec2().x()) +
523  (pzh*posDist->GetSideRefVec3().x());
524 
525  G4double resulty = (pxh*posDist->GetSideRefVec1().y()) + (pyh*posDist->GetSideRefVec2().y()) +
526  (pzh*posDist->GetSideRefVec3().y());
527 
528  G4double resultz = (pxh*posDist->GetSideRefVec1().z()) + (pyh*posDist->GetSideRefVec2().z()) +
529  (pzh*posDist->GetSideRefVec3().z());
530 
531  G4double ResMag = std::sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz));
532  resultx = resultx/ResMag;
533  resulty = resulty/ResMag;
534  resultz = resultz/ResMag;
535 
536  mom.setX(resultx);
537  mom.setY(resulty);
538  mom.setZ(resultz);
539  }
540 
541  // particle_momentum_direction now contains unit momentum vector.
542  if(verbosityLevel > 0 )
543  {
544  G4cout << "Final User Defined momentum vector " << particle_momentum_direction << G4endl;
545  }
546 }
547 
549 {
550  // Create cumulative histogram if not already done so. Then use RandFlat
551  //::shoot to generate the output Theta value.
552  if(UserDistType == "NULL" || UserDistType == "phi")
553  {
554  // No user defined theta distribution
555  G4cout << "Error ***********************" << G4endl;
556  G4cout << "UserDistType = " << UserDistType << G4endl;
557  return (0.);
558  }
559  else
560  {
561  // UserDistType = theta or both and so a theta distribution
562  // is defined. This should be integrated if not already done.
563  G4AutoLock l(&mutex);
564  if(IPDFThetaExist == false)
565  {
566  // IPDF has not been created, so create it
567  G4double bins[1024],vals[1024], sum;
568  G4int ii;
570  bins[0] = UDefThetaH.GetLowEdgeEnergy(size_t(0));
571  vals[0] = UDefThetaH(size_t(0));
572  sum = vals[0];
573  for(ii=1;ii<maxbin;ii++)
574  {
575  bins[ii] = UDefThetaH.GetLowEdgeEnergy(size_t(ii));
576  vals[ii] = UDefThetaH(size_t(ii)) + vals[ii-1];
577  sum = sum + UDefThetaH(size_t(ii));
578  }
579  for(ii=0;ii<maxbin;ii++)
580  {
581  vals[ii] = vals[ii]/sum;
582  IPDFThetaH.InsertValues(bins[ii], vals[ii]);
583  }
584  // Make IPDFThetaExist = true
585  IPDFThetaExist = true;
586  }
587  l.unlock();
588  // IPDF has been create so carry on
589  G4double rndm = G4UniformRand();
590  return(IPDFThetaH.GetEnergy(rndm));
591  }
592 }
593 
595 {
596  // Create cumulative histogram if not already done so. Then use RandFlat
597  //::shoot to generate the output Theta value.
598 
599  if(UserDistType == "NULL" || UserDistType == "theta")
600  {
601  // No user defined phi distribution
602  G4cout << "Error ***********************" << G4endl;
603  G4cout << "UserDistType = " << UserDistType << G4endl;
604  return(0.);
605  }
606  else
607  {
608  // UserDistType = phi or both and so a phi distribution
609  // is defined. This should be integrated if not already done.
610  G4AutoLock l(&mutex);
611  if(IPDFPhiExist == false)
612  {
613  // IPDF has not been created, so create it
614  G4double bins[1024],vals[1024], sum;
615  G4int ii;
616  G4int maxbin = G4int(UDefPhiH.GetVectorLength());
617  bins[0] = UDefPhiH.GetLowEdgeEnergy(size_t(0));
618  vals[0] = UDefPhiH(size_t(0));
619  sum = vals[0];
620  for(ii=1;ii<maxbin;ii++)
621  {
622  bins[ii] = UDefPhiH.GetLowEdgeEnergy(size_t(ii));
623  vals[ii] = UDefPhiH(size_t(ii)) + vals[ii-1];
624  sum = sum + UDefPhiH(size_t(ii));
625  }
626  for(ii=0;ii<maxbin;ii++)
627  {
628  vals[ii] = vals[ii]/sum;
629  IPDFPhiH.InsertValues(bins[ii], vals[ii]);
630  }
631  // Make IPDFPhiExist = true
632  IPDFPhiExist = true;
633  }
634  l.unlock();
635  // IPDF has been create so carry on
636  G4double rndm = G4UniformRand();
637  return(IPDFPhiH.GetEnergy(rndm));
638  }
639 }
640 //
642 {
643  G4AutoLock l(&mutex);
644  if (atype == "theta") {
646  IPDFThetaExist = false ;}
647  else if (atype == "phi"){
649  IPDFPhiExist = false ;}
650  else {
651  G4cout << "Error, histtype not accepted " << G4endl;
652  }
653 }
654 
655 
657 {
658  //Local copy for thread safety
660  // Angular stuff
661  if(AngDistType == "iso")
662  GenerateIsotropicFlux(localM);
663  else if(AngDistType == "cos")
664  GenerateCosineLawFlux(localM);
665  else if(AngDistType == "planar")
666  GeneratePlanarFlux(localM);
667  else if(AngDistType == "beam1d" || AngDistType == "beam2d" )
668  GenerateBeamFlux(localM);
669  else if(AngDistType == "user")
670  GenerateUserDefFlux(localM);
671  else if(AngDistType == "focused")
672  GenerateFocusedFlux(localM);
673  else
674  G4cout << "Error: AngDistType has unusual value" << G4endl;
675  return localM;
676 }
677 
678 
679 
680 
681 
682 
683 
684 
685