96 G4cout <<
"Low energy photon Compton model is constructed " <<
G4endl;
124 G4cout <<
"Calling G4LowEPPolarizedComptonModel::Initialise()" <<
G4endl;
133 char* path = std::getenv(
"G4LEDATA");
139 for(
G4int i=0; i<numOfCouples; ++i) {
145 for (
G4int j=0; j<nelm; ++j) {
171 G4cout <<
"G4LowEPPolarizedComptonModel is initialized " <<
G4endl
199 G4cout <<
"G4LowEPPolarizedComptonModel::ReadData()"
202 if(
data[Z]) {
return; }
203 const char* datadir = path;
206 datadir = std::getenv(
"G4LEDATA");
209 G4Exception(
"G4LowEPPolarizedComptonModel::ReadData()",
211 "Environment variable G4LEDATA not defined");
219 data[
Z]->SetSpline(
false);
221 std::ostringstream ost;
222 ost << datadir <<
"/livermore/comp/ce-cs-" << Z <<
".dat";
223 std::ifstream
fin(ost.str().c_str());
228 ed <<
"G4LowEPPolarizedComptonModel data file <" << ost.str().c_str()
229 <<
"> is not opened!" <<
G4endl;
230 G4Exception(
"G4LowEPPolarizedComptonModel::ReadData()",
232 ed,
"G4LEDATA version should be G4EMLOW6.34 or later");
236 G4cout <<
"File " << ost.str()
237 <<
" is opened by G4LowEPPolarizedComptonModel" <<
G4endl;
255 G4cout <<
"G4LowEPPolarizedComptonModel::ComputeCrossSectionPerAtom()"
263 if(intZ < 1 || intZ >
maxZ) {
return cs; }
273 if(!pv) {
return cs; }
280 if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*
e1)*pv->
Value(e1); }
281 else if(GammaEnergy <= e2) { cs = pv->
Value(GammaEnergy)/GammaEnergy; }
282 else if(GammaEnergy > e2) { cs = pv->
Value(e2)/GammaEnergy; }
297 G4double g4d_limit = std::pow(10.,-g4d_order);
317 G4cout <<
"G4LowEPPolarizedComptonModel::SampleSecondaries() E(MeV)= "
337 if (!(photonPolarization0.
isOrthogonal(photonDirection0, 1
e-6))||(photonPolarization0.
mag()==0))
344 if ((photonPolarization0.
howOrthogonal(photonDirection0) !=0) && (photonPolarization0.
howOrthogonal(photonDirection0) > g4d_limit))
356 G4double LowEPPCepsilon0 = 1. / (1. + 2. * e0m);
357 G4double LowEPPCepsilon0Sq = LowEPPCepsilon0 * LowEPPCepsilon0;
358 G4double alpha1 = -std::log(LowEPPCepsilon0);
359 G4double alpha2 = 0.5 * (1. - LowEPPCepsilon0Sq);
371 G4cout <<
"Started loop to sample gamma energy" <<
G4endl;
379 LowEPPCepsilonSq = LowEPPCepsilon * LowEPPCepsilon;
383 LowEPPCepsilonSq = LowEPPCepsilon0Sq + (1. - LowEPPCepsilon0Sq) *
G4UniformRand();
384 LowEPPCepsilon = std::sqrt(LowEPPCepsilonSq);
387 oneCosT = (1. - LowEPPCepsilon) / ( LowEPPCepsilon * e0m);
388 sinT2 = oneCosT * (2. - oneCosT);
389 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/
cm);
391 gReject = (1. - LowEPPCepsilon * sinT2 / (1. + LowEPPCepsilonSq)) * scatteringFunction;
396 G4double sinTheta = std::sqrt(sinT2);
398 G4double dirx = sinTheta * std::cos(phi);
399 G4double diry = sinTheta * std::sin(phi);
420 const G4int maxDopplerIterations = 1000;
422 G4double pEIncident = photonEnergy0 ;
438 G4cout <<
"Started loop to sample photon energy and electron direction" <<
G4endl;
467 G4double ePSI = ePAU * momentum_au_to_nat;
470 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c;
480 G4double systemE = eEIncident + pEIncident;
483 G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
485 G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
486 G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
488 pERecoil = (numerator/denominator);
489 eERecoil = systemE - pERecoil;
490 CE_emission_flag = pEIncident - pERecoil;
491 }
while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
504 G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
508 G4double sinAlpha = std::sin(e_alpha);
509 G4double cosAlpha = std::cos(e_alpha);
510 G4double sinBeta = std::sin(e_beta);
511 G4double cosBeta = std::cos(e_beta);
513 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
514 G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
516 G4double var_A = pERecoil*u_p_temp*sinTheta;
517 G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
521 G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
523 G4double var_D = var_D1*var_D2 + var_D3;
526 G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
529 G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
530 G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
533 G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
538 G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
539 G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
542 G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
544 G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
545 G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
557 if ((diff < 0.0) && (
abs(diff / diff1) < g4d_limit) && (
abs(diff / diff2) < g4d_limit) )
563 G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
564 G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
571 if(X_p >1){X_p=1;}
if(X_p<-1){X_p=-1;}
572 if(X_m >1){X_m=1;}
if(X_m<-1){X_m=-1;}
582 if (sol_select < 0.5)
584 ThetaE = std::acos(X_p);
586 if (sol_select > 0.5)
588 ThetaE = std::acos(X_m);
591 cosThetaE = std::cos(ThetaE);
592 sinThetaE = std::sin(ThetaE);
593 G4double Theta = std::acos(cosTheta);
596 G4double iSinThetaE = std::sqrt(1+std::tan((
pi/2.0)-ThetaE)*std::tan((
pi/2.0)-ThetaE));
597 G4double iSinTheta = std::sqrt(1+std::tan((
pi/2.0)-Theta)*std::tan((
pi/2.0)-Theta));
598 G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
600 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
606 }
while ( (iteration <= maxDopplerIterations) && (
abs(cosPhiE) > 1));
609 if (iteration >= maxDopplerIterations)
611 pERecoil = photonEnergy0 ;
622 photonPolarization0,photonPolarization1);
633 G4double eDirX = sinThetaE * std::cos(phi+PhiE);
634 G4double eDirY = sinThetaE * std::sin(phi+PhiE);
637 G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
641 photonPolarization0);
644 eDirection,eKineticEnergy) ;
645 fvect->push_back(dp);
664 size_t nbefore = fvect->size();
668 size_t nafter = fvect->size();
669 if(nafter > nbefore) {
670 for (
size_t i=nbefore; i<nafter; ++i) {
672 if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
675 bindingE -= ((*fvect)[i])->GetKineticEnergy();
691 G4Exception(
"G4LowEPPolarizedComptonModel::SampleSecondaries()",
729 G4AutoLock l(&LowEPPolarizedComptonModelMutex);
739 { 0, 0., 0., 0., 0., 0., 0., 0., 0.},
740 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143.374, 50.787, -5.951, 0.2304418},
741 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53.649, 13.892, -0.948, 0.006996759},
742 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62.090, 21.462, -2.453, 0.093416},
743 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127.906, 46.491, -5.614, 0.2262103},
744 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131.153, 47.132, -5.619, 0.2233819},
745 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128.143, 45.379, -5.325, 0.2083009},
746 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131.048, 46.314, -5.421, 0.2114925},
747 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128.225, 44.818, -5.183, 0.1997155},
748 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122.112, 42.103, -4.796, 0.1819099},
749 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110.143, 37.225, -4.143, 0.1532094},
750 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41.137, 12.313, -1.152, 0.03384553},
751 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53.549, 17.420, -1.840, 0.06431849},
752 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66.243, 22.297, -2.460, 0.09045854},
753 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78.271, 26.757, -3.008, 0.1128195},
754 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85.069, 29.164, -3.291, 0.1239113},
755 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93.640, 32.274, -3.665, 0.1388633},
756 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98.534, 33.958, -3.857, 0.1461557},
757 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100.077, 34.379, -3.891, 0.1468902},
758 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53.819, 17.528, -1.851, 0.0648722},
759 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52.221, 17.169, -1.832, 0.06502094},
760 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55.365, 18.276, -1.961, 0.07002778},
761 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57.412, 18.957, -2.036, 0.07278856},
762 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58.478, 19.270, -2.065, 0.07362722},
763 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62.192, 20.358, -2.167, 0.07666583},
764 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58.007, 18.924, -2.003, 0.0704305},
765 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61.176, 20.067, -2.141, 0.0760269},
766 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61.909, 20.271, -2.159, 0.07653559},
767 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62.402, 20.391, -2.167, 0.07664243},
768 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67.305, 21.954, -2.331, 0.0823267},
769 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62.064, 20.136, -2.122, 0.07437589},
770 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61.068, 19.808, -2.086, 0.07307488},
771 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63.126, 20.553, -2.175, 0.07660222},
772 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65.674, 21.445, -2.278, 0.08054694},
773 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69.457, 22.811, -2.442, 0.08709536},
774 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72.283, 23.808, -2.558, 0.09156808},
775 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74.696, 24.641, -2.653, 0.09516597},
776 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46.235, 14.519, -1.458, 0.04837659},
777 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41.784, 13.065, -1.300, 0.04267703},
778 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44.609, 14.114, -1.429, 0.0479348},
779 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47.142, 15.019, -1.536, 0.0521347},
780 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53.114, 17.052, -1.766, 0.06079426},
781 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54.590, 17.550, -1.822, 0.06290335},
782 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51.272, 16.423, -1.694, 0.05806108},
783 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58.314, 18.839, -1.969, 0.0684608},
784 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59.674, 19.295, -2.020, 0.07037294},
785 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72.228, 23.693, -2.532, 0.09017969},
786 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60.890, 19.647, -2.053, 0.07138694},
787 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56.152, 18.002, -1.863, 0.06410123},
788 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56.208, 18.052, -1.872, 0.06456884},
789 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58.478, 18.887, -1.973, 0.06860356},
790 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60.708, 19.676, -2.066, 0.07225841},
791 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63.341, 20.632, -2.180, 0.0767412},
792 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66.339, 21.716, -2.310, 0.08191981},
793 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67.649, 22.151, -2.357, 0.08357825},
794 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45.302, 14.219, -1.423, 0.04712317},
795 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39.825, 12.363, -1.214, 0.03931009},
796 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38.952, 11.982, -1.160, 0.03681554},
797 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41.959, 13.118, -1.302, 0.04271291},
798 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40.528, 12.555, -1.230, 0.03971407},
799 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39.986, 12.329, -1.200, 0.03843737},
800 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42.756, 13.362, -1.327, 0.0436124},
801 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42.682, 13.314, -1.319, 0.04322509},
802 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42.399, 13.185, -1.301, 0.04243861},
803 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43.226, 13.475, -1.335, 0.04377341},
804 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43.232, 13.456, -1.330, 0.04347536},
805 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42.047, 12.990, -1.270, 0.04095499},
806 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42.405, 13.106, -1.283, 0.04146024},
807 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41.974, 12.926, -1.259, 0.040435},
808 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42.132, 12.967, -1.262, 0.04048908},
809 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42.582, 13.122, -1.280, 0.04119599},
810 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42.586, 13.115, -1.278, 0.04107587},
811 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43.708, 13.509, -1.324, 0.04286491},
812 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44.838, 13.902, -1.369, 0.04457132},
813 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45.545, 14.137, -1.395, 0.04553459},
814 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46.426, 14.431, -1.428, 0.04678218},
815 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47.513, 14.806, -1.471, 0.04842566},
816 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48.225, 15.042, -1.497, 0.04938364},
817 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53.256, 16.739, -1.687, 0.05645173},
818 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53.900, 16.946, -1.709, 0.05723134},
819 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49.801, 15.536, -1.547, 0.05103522},
820 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49.651, 15.512, -1.548, 0.05123203},
821 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51.021, 16.018, -1.609, 0.05364831},
822 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52.693, 16.612, -1.679, 0.05638698},
823 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54.415, 17.238, -1.754, 0.05935566},
824 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56.083, 17.834, -1.824, 0.06206068},
825 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57.860, 18.463, -1.898, 0.0649633},
826 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39.973, 12.164, -1.162, 0.0364598},
827 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34.591, 10.338, -0.956, 0.0287409},
828 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36.002, 10.867, -1.021, 0.03136835},
829 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37.660, 11.475, -1.095, 0.03435334},
830 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37.268, 11.301, -1.071, 0.0330539},
831 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37.695, 11.438, -1.085, 0.03376669},
832 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39.010, 11.927, -1.146, 0.03630848},
833 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37.192, 11.229, -1.057, 0.0325621},
834 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37.589, 11.363, -1.072, 0.03312393},
835 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39.573, 12.095, -1.162, 0.03680527},
836 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40.007, 12.242, -1.178, 0.03737377},
837 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39.509, 12.041, -1.152, 0.03629023},
838 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39.939, 12.183, -1.168, 0.03690464},
839 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39.973, 12.180, -1.166, 0.036773}
863 b = energyRate + 1/energyRate;
865 phiProbability = 1 - (a/
b)*(std::cos(phi)*std::cos(phi));
870 while ( rand2 > phiProbability );
906 c.
setX(std::cos(angle)*(a0.
x())+std::sin(angle)*b0.
x());
907 c.
setY(std::cos(angle)*(a0.
y())+std::sin(angle)*b0.
y());
908 c.
setZ(std::cos(angle)*(a0.
z())+std::sin(angle)*b0.
z());
932 return photonPolarization - photonPolarization.
dot(photonDirection)/photonDirection.
dot(photonDirection) * photonDirection;
946 G4double sinTheta = std::sqrt(sinT2);
948 G4double normalisation = std::sqrt(1. - cosP2*sinT2);
963 if (rand1<(LowEPPCepsilon+1.0/LowEPPCepsilon-2)/(2.0*(LowEPPCepsilon+1.0/LowEPPCepsilon)-4.0*sinT2*cosP2))
978 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
982 G4double xParallel = normalisation*cosBeta;
983 G4double yParallel = -(sinT2*cosPhi*sinPhi)*cosBeta/normalisation;
984 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
986 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
987 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
989 G4double xTotal = (xParallel + xPerpendicular);
990 G4double yTotal = (yParallel + yPerpendicular);
991 G4double zTotal = (zParallel + zPerpendicular);
993 photonPolarization1.
setX(xTotal);
994 photonPolarization1.
setY(yTotal);
995 photonPolarization1.
setZ(zTotal);
997 return photonPolarization1;
1016 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
1021 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
1040 edirection = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();