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