5 #include <phparameter/PHParameters.h>
11 #include <Geant4/G4AssemblyVolume.hh>
12 #include <Geant4/G4Box.hh>
13 #include <Geant4/G4Colour.hh>
14 #include <Geant4/G4Cons.hh>
15 #include <Geant4/G4Element.hh>
16 #include <Geant4/G4IntersectionSolid.hh>
17 #include <Geant4/G4LogicalBorderSurface.hh>
18 #include <Geant4/G4LogicalSkinSurface.hh>
19 #include <Geant4/G4LogicalVolume.hh>
20 #include <Geant4/G4Material.hh>
21 #include <Geant4/G4OpticalSurface.hh>
22 #include <Geant4/G4PVPlacement.hh>
23 #include <Geant4/G4ProcessManager.hh>
24 #include <Geant4/G4RotationMatrix.hh>
25 #include <Geant4/G4RunManager.hh>
26 #include <Geant4/G4Sphere.hh>
27 #include <Geant4/G4SubtractionSolid.hh>
28 #include <Geant4/G4SystemOfUnits.hh>
29 #include <Geant4/G4ThreeVector.hh>
30 #include <Geant4/G4Transform3D.hh>
31 #include <Geant4/G4Trap.hh>
32 #include <Geant4/G4Trd.hh>
33 #include <Geant4/G4Tubs.hh>
34 #include <Geant4/G4UImanager.hh>
35 #include <Geant4/G4UnionSolid.hh>
36 #include <Geant4/G4VUserDetectorConstruction.hh>
37 #include <Geant4/G4VisAttributes.hh>
48 const std::string& dnam)
50 , m_Params(parameters)
64 if (!volume)
return 0;
323 if (
Verbosity()) std::cout <<
"barrel radius = " <<
fRadius <<
" mm" << std::endl;
376 double gluethickness = 0.05;
379 double dirclength =
fBarL[2] * (nparts - 1) +
fBarS[2] + gluethickness * nparts;
417 G4Box* gGlue =
new G4Box(
"gGlue", fBar[0] / 2., fBar[1] / 2., 0.5 * gluethickness);
423 for (
int i = 0; i <
fNBar; i++)
426 for (
int j = 0; j < nparts; j++)
428 if (j < (nparts - 1))
430 double z = 0.5 * dirclength - 0.5 *
fBarL[2] - (
fBarL[2] + gluethickness) * j;
442 double z = 0.5 * dirclength - (nparts - 1) *
fBarL[2] - 0.5 *
fBarS[2] - gluethickness * j;
469 double lensMinThikness = 2;
485 double lensMinThikness = 2;
495 double cr2 = sqrt(
fLens[1] *
fLens[1] / 4. + fBar[0] * fBar[0] / 4.);
498 if (
Verbosity()) std::cout <<
"bad lens" << std::endl;
501 fLens[2] = (2 * lensMinThikness + r2 - sqrt(r2 * r2 - cr2 * cr2) + lensMinThikness);
502 if (
Verbosity()) std::cout <<
"lens thickness =" <<
fLens[2] <<
" mm" << std::endl;
504 G4ThreeVector zTrans1(0, 0, -r1 -
fLens[2] / 2. + r1 - sqrt(r1 * r1 - cr2 * cr2) + lensMinThikness);
505 G4ThreeVector zTrans2(0, 0, -r2 -
fLens[2] / 2. + r2 - sqrt(r2 * r2 - cr2 * cr2) + lensMinThikness * 2);
534 double lensMinThikness = 2.0;
540 double layer12 = lensMinThikness * 2;
549 G4ThreeVector zTrans1(0, 0, -r1 -
fLens[2] / 2. + r1 - sqrt(r1 * r1 - shight / 2. * shight / 2.) + lensMinThikness);
550 G4ThreeVector zTrans2(0, 0, -r2 -
fLens[2] / 2. + r2 - sqrt(r2 * r2 - shight / 2. * shight / 2.) + layer12);
590 double shifth = -0.5 * (dirclength +
fLens[2]);
597 else if (fNBar == 1 &&
fLensId == 3)
599 for (
int i = 0; i < 11; i++)
613 for (
int i = 0; i <
fNBar; i++)
672 for (
int i = 0; i <
fNBoxes; i++)
674 double tphi = dphi * i;
730 double WaveLength[
num];
731 double PhotonEnergy[
num];
732 double PMTReflectivity[
num];
733 double EfficiencyMirrors[
num];
734 const double LambdaE = 2.0 * 3.14159265358979323846 * 1.973269602e-16 *
m *
GeV;
735 for (
int i = 0; i <
num; i++)
737 WaveLength[i] = (300 + i * 10) *
nanometer;
738 PhotonEnergy[num - (i + 1)] = LambdaE / WaveLength[i];
739 PMTReflectivity[i] = 0.;
740 EfficiencyMirrors[i] = 0;
746 double QuantumEfficiencyIdial[
num] =
747 {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
748 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
749 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
750 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
753 double QuantumEfficiencyB[
num] =
754 {0., 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, 0.06,
755 0.07, 0.09, 0.1, 0.13, 0.15, 0.17, 0.2, 0.24, 0.26, 0.28, 0.282, 0.284, 0.286,
756 0.288, 0.29, 0.28, 0.26, 0.24, 0.22, 0.20, 0.18, 0.15, 0.13, 0.12, 0.10};
759 double QuantumEfficiencyPMT[
num] =
760 {0.001, 0.002, 0.004, 0.007, 0.011, 0.015, 0.020, 0.026, 0.033, 0.040, 0.045,
761 0.056, 0.067, 0.085, 0.109, 0.129, 0.138, 0.147, 0.158, 0.170,
762 0.181, 0.188, 0.196, 0.203, 0.206, 0.212, 0.218, 0.219, 0.225, 0.230,
763 0.228, 0.222, 0.217, 0.210, 0.199, 0.177};
769 for (
int k = 0;
k < 36;
k++)
771 QuantumEfficiencyB[
k] = QuantumEfficiencyB[
k] * 0.45;
772 QuantumEfficiencyPMT[
k] = QuantumEfficiencyPMT[
k] * .7;
780 PhotocatodBurleMPT->
AddProperty(
"EFFICIENCY", PhotonEnergy, QuantumEfficiencyB, num);
781 PhotocatodBurleMPT->
AddProperty(
"REFLECTIVITY", PhotonEnergy, PMTReflectivity, num);
794 PhotocatodHamamatsuMPT->
AddProperty(
"REFLECTIVITY", PhotonEnergy, PMTReflectivity, num);
808 double ReflectivityMirrorBar[
num] = {
809 0.8755, 0.882, 0.889, 0.895, 0.9, 0.9025, 0.91, 0.913, 0.9165, 0.92, 0.923,
810 0.9245, 0.9285, 0.932, 0.934, 0.935, 0.936, 0.9385, 0.9395, 0.94,
811 0.9405, 0.9405, 0.9405, 0.9405, 0.94, 0.9385, 0.936, 0.934,
812 0.931, 0.9295, 0.928, 0.928, 0.921, 0.92, 0.927, 0.9215};
815 MirrorMPT->
AddProperty(
"REFLECTIVITY", PhotonEnergy, ReflectivityMirrorBar, num);
816 MirrorMPT->
AddProperty(
"EFFICIENCY", PhotonEnergy, EfficiencyMirrors, num);
834 double a,
z, density;
836 int ncomponents, natoms;
904 static const double LambdaE = 2.0 * 3.14159265358979323846 * 1.973269602e-16 *
m *
GeV;
906 double WaveLength[
num];
908 double AirAbsorption[
num];
909 double AirRefractiveIndex[
num];
910 double PhotonEnergy[
num];
913 double PhotonEnergyNlak33a[76] = {1, 1.2511, 1.26386, 1.27687, 1.29016, 1.30372, 1.31758, 1.33173, 1.34619, 1.36097, 1.37607, 1.39152, 1.40731, 1.42347, 1.44, 1.45692, 1.47425, 1.49199, 1.51016, 1.52878, 1.54787, 1.56744, 1.58751, 1.6081, 1.62923, 1.65092, 1.6732, 1.69609, 1.71961, 1.7438, 1.76868, 1.79427, 1.82062, 1.84775, 1.87571, 1.90452, 1.93423, 1.96488, 1.99652, 2.0292, 2.06296, 2.09787, 2.13398, 2.17135, 2.21006, 2.25017, 2.29176, 2.33492, 2.37973, 2.42631, 2.47473, 2.52514, 2.57763, 2.63236, 2.68946, 2.7491, 2.81143, 2.87666, 2.94499, 3.01665, 3.09187, 3.17095, 3.25418, 3.34189, 3.43446, 3.53231, 3.6359, 3.74575, 3.86244, 3.98663, 4.11908, 4.26062, 4.41225, 4.57506, 4.75035, 4.93961};
916 double en_PbF2[] = {1.55, 1.569, 1.59, 1.61, 1.631, 1.653, 1.675, 1.698, 1.722, 1.746, 1.771, 1.797, 1.823, 1.851, 1.879, 1.907, 1.937, 1.968, 2, 2.033, 2.066, 2.101, 2.138, 2.175, 2.214, 2.254, 2.296, 2.339, 2.384, 2.431, 2.48, 2.53, 2.583, 2.638, 2.695, 2.755, 2.818, 2.883, 2.952, 3.024, 3.1, 3.179, 3.263, 3.351, 3.444, 3.542, 3.647, 3.757, 3.875, 3.999, 4.133, 4.275, 4.428, 4.592, 4.769, 4.959};
917 double ab_PbF2[] = {407, 403.3, 379.1, 406.3, 409.7, 408.9, 406.7, 404.7, 391.7, 397.7, 409.6, 403.7, 403.8, 409.7, 404.9, 404.2, 407.1, 411.1, 403.1, 406.1, 415.4, 399.1, 405.8, 408.2, 385.7, 405.6, 405.2, 401.6, 402.6, 407.1, 417.7, 401.1, 389.9, 411.9, 400.9, 398.3, 402.1, 408.7, 384.8, 415.8, 413.1, 385.7, 353.7, 319.1, 293.6, 261.9, 233.6, 204.4, 178.3, 147.6, 118.2, 78.7, 51.6, 41.5, 24.3, 8.8};
918 double ref_PbF2[] = {1.749, 1.749, 1.75, 1.75, 1.751, 1.752, 1.752, 1.753, 1.754, 1.754, 1.755, 1.756, 1.757, 1.757, 1.758, 1.759, 1.76, 1.761, 1.762, 1.764, 1.765, 1.766, 1.768, 1.769, 1.771, 1.772, 1.774, 1.776, 1.778, 1.78, 1.782, 1.785, 1.787, 1.79, 1.793, 1.796, 1.8, 1.804, 1.808, 1.813, 1.818, 1.824, 1.83, 1.837, 1.845, 1.854, 1.865, 1.877, 1.892, 1.91, 1.937, 1.991, 1.38, 1.915, 1.971, 2.019};
920 const int n_Sapphire = 57;
921 double en_Sapphire[] = {1.02212, 1.05518, 1.09045, 1.12610, 1.16307, 1.20023, 1.23984, 1.28043, 1.32221, 1.36561, 1.41019, 1.45641, 1.50393, 1.55310, 1.60393, 1.65643, 1.71059, 1.76665, 1.82437, 1.88397, 1.94576, 2.00946, 2.07504, 2.14283, 2.21321, 2.28542, 2.36025, 2.43727, 2.51693, 2.59924, 2.68480, 2.77245, 2.86337, 2.95693, 3.05379, 3.15320, 3.25674, 3.36273, 3.47294, 3.58646, 3.70433, 3.82549, 3.94979, 4.07976, 4.21285, 4.35032, 4.49380, 4.64012, 4.79258, 4.94946, 5.11064, 5.27816, 5.44985, 5.62797, 5.81266, 6.00407, 6.19920};
922 double ref_Sapphire[] = {1.75188, 1.75253, 1.75319, 1.75382, 1.75444, 1.75505, 1.75567, 1.75629, 1.75691, 1.75754, 1.75818, 1.75883, 1.75949, 1.76017, 1.76088, 1.76160, 1.76235, 1.76314, 1.76395, 1.76480, 1.76570, 1.76664, 1.76763, 1.76867, 1.76978, 1.77095, 1.77219, 1.77350, 1.77490, 1.77639, 1.77799, 1.77968, 1.78150, 1.78343, 1.78551, 1.78772, 1.79011, 1.79265, 1.79540, 1.79835, 1.80154, 1.80497, 1.80864, 1.81266, 1.81696, 1.82163, 1.82674, 1.83223, 1.83825, 1.84480, 1.85191, 1.85975, 1.86829, 1.87774, 1.88822, 1.89988, 1.91270};
924 double ab_Sapphire[n_Sapphire];
926 double transmitance_Sapphire[] = {0.845, 0.844, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.843, 0.842, 0.842, 0.842, 0.842, 0.838, 0.838, 0.838, 0.838, 0.838, 0.838, 0.836, 0.836, 0.836, 0.836, 0.834, 0.834, 0.833, 0.832, 0.832, 0.831, 0.831, 0.83, 0.828, 0.828, 0.828, 0.828, 0.828, 0.828, 0.828, 0.828, 0.828, 0.825, 0.80, 0.67, 0.47, 0.23, 0.22};
927 for (
int i = 0; i < n_Sapphire; i++) ab_Sapphire[i] = (-1) / log(transmitance_Sapphire[i] + 2 * 0.077007) * 2 *
mm;
932 double KamLandOilAbsorption[
num] =
933 {0.97469022, 0.976603956, 0.978511548, 0.980400538, 0.982258449, 0.984072792,
934 0.985831062, 0.987520743, 0.989129303, 0.990644203, 0.992052894,
935 0.993342822, 0.994501428, 0.995516151, 0.996374433, 0.997063719,
936 0.997571464, 0.997885132, 0.997992205, 0.997880183, 0.997536591,
937 0.99, 0.98, 0.97, 0.96, 0.94, 0.93, 0.924507, 0.89982, 0.883299,
938 0.85657, 0.842637, 0.77020213, 0.65727, 0.324022, 0.019192};
941 double QuartzAbsorption[
num] =
942 {0.999572036, 0.999544661, 0.999515062, 0.999483019, 0.999448285,
943 0.999410586, 0.999369611, 0.999325013, 0.999276402, 0.999223336,
944 0.999165317, 0.999101778, 0.999032079, 0.998955488, 0.998871172,
945 0.998778177, 0.99867541, 0.998561611, 0.998435332, 0.998294892, 0.998138345,
946 0.997963425, 0.997767484, 0.997547418,
947 0.99729958, 0.99701966, 0.99670255, 0.996342167, 0.995931242, 0.995461041,
948 0.994921022, 0.994298396, 0.993577567, 0.992739402, 0.991760297, 0.990610945};
951 double EpotekAbsorption[
num] =
952 {0.99999999, 0.99999999, 0.99999999, 0.99999999, 0.99999999,
953 0.99999999, 0.99999999, 0.99999999, 0.99999999, 0.99999999,
954 0.99999999, 0.99999999, 0.99999999, 0.99999999, 0.99999999,
955 0.99999999, 0.99999999, 0.99999999, 0.99999999, 0.99999999,
956 0.99999999, 0.99999999, 0.99999999, 0.99999999, 0.99999999,
957 0.9999, 0.9998, 0.9995, 0.999, 0.998, 0.997, 0.996, 0.9955, 0.993,
961 double Nlak33aAbsorption[76] = {371813, 352095, 331021, 310814, 291458, 272937, 255238, 238342, 222234, 206897, 192313, 178463, 165331, 152896, 141140, 130043, 119585, 109747, 100507, 91846.3, 83743.1, 76176.7, 69126.1, 62570.2, 56488, 50858.3, 45660.1, 40872.4, 36474.6, 32445.8, 28765.9, 25414.6, 22372.2, 19619.3, 17136.9, 14906.5, 12910.2, 11130.3, 9550.13, 8153.3, 6924.25, 5848.04, 4910.46, 4098.04, 3398.06, 2798.54, 2288.32, 1856.99, 1494.92, 1193.28, 943.973, 739.657, 573.715, 440.228, 333.94, 250.229, 185.064, 134.967, 96.9664, 68.5529, 47.6343, 32.4882, 21.7174, 14.2056, 9.07612, 5.65267, 3.4241, 2.01226, 1.14403, 0.62722, 0.330414, 0.166558, 0.0799649, 0.0363677, 0.0155708, 0.00623089};
963 double EpotekThickness = 0.001 * 2.54 *
cm;
964 for (
int i = 0; i <
num; i++)
966 WaveLength[i] = (300 + i * 10) *
nanometer;
968 AirAbsorption[i] = 4. *
cm;
969 AirRefractiveIndex[i] = 1.;
970 PhotonEnergy[num - (i + 1)] = LambdaE / WaveLength[i];
977 EpotekAbsorption[i] = (-1) / log(EpotekAbsorption[i]) * EpotekThickness;
978 QuartzAbsorption[i] = (-1) / log(QuartzAbsorption[i]) * 100 *
cm;
979 KamLandOilAbsorption[i] = (-1) / log(KamLandOilAbsorption[i]) * 50 *
cm;
986 double QuartzRefractiveIndex[
num] = {
987 1.456535, 1.456812, 1.4571, 1.457399, 1.457712, 1.458038, 1.458378,
988 1.458735, 1.459108, 1.4595, 1.459911, 1.460344, 1.460799, 1.46128,
989 1.461789, 1.462326, 1.462897, 1.463502, 1.464146, 1.464833,
990 1.465566, 1.46635, 1.46719, 1.468094, 1.469066, 1.470116, 1.471252, 1.472485,
991 1.473826, 1.475289, 1.476891, 1.478651, 1.480592, 1.482739, 1.485127, 1.487793};
993 double EpotekRefractiveIndex[
num] = {
994 1.554034, 1.555575, 1.55698, 1.558266, 1.559454, 1.56056, 1.561604,
995 1.562604, 1.563579, 1.564547, 1.565526, 1.566536, 1.567595,
996 1.568721, 1.569933, 1.57125, 1.57269, 1.574271, 1.576012,
997 1.577932, 1.580049, 1.582381, 1.584948, 1.587768, 1.590859,
998 1.59424, 1.597929, 1.601946, 1.606307, 1.611033, 1.616141, 1.621651, 1.62758,
999 1.633947, 1.640771, 1.64807};
1001 double KamLandOilRefractiveIndex[
num] = {
1002 1.433055, 1.433369, 1.433698, 1.434045, 1.434409, 1.434793, 1.435198,
1003 1.435626, 1.436077, 1.436555, 1.4371, 1.4376, 1.4382, 1.4388, 1.4395,
1004 1.4402, 1.4409, 1.4415, 1.4425, 1.4434, 1.4444, 1.4455, 1.4464, 1.4479, 1.4501,
1005 1.450428, 1.451976, 1.453666, 1.455513, 1.45754, 1.45977, 1.462231, 1.464958,
1006 1.467991, 1.471377, 1.475174};
1008 double Nlak33aRefractiveIndex[76] = {1.73816, 1.73836, 1.73858, 1.73881, 1.73904, 1.73928, 1.73952, 1.73976, 1.74001, 1.74026, 1.74052, 1.74078, 1.74105, 1.74132, 1.7416, 1.74189, 1.74218, 1.74249, 1.74279, 1.74311, 1.74344, 1.74378, 1.74412, 1.74448, 1.74485, 1.74522, 1.74562, 1.74602, 1.74644, 1.74687, 1.74732, 1.74779, 1.74827, 1.74878, 1.7493, 1.74985, 1.75042, 1.75101, 1.75163, 1.75228, 1.75296, 1.75368, 1.75443, 1.75521, 1.75604, 1.75692, 1.75784, 1.75882, 1.75985, 1.76095, 1.76211, 1.76335, 1.76467, 1.76608, 1.76758, 1.7692, 1.77093, 1.77279, 1.7748, 1.77698, 1.77934, 1.7819, 1.7847, 1.78775, 1.79111, 1.79481, 1.79889, 1.80343, 1.8085, 1.81419, 1.82061, 1.8279, 1.83625, 1.84589, 1.85713, 1.87039};
1015 static bool once =
true;
1019 std::cout << __PRETTY_FUNCTION__ <<
" : warning parameter disable_photon_sim = " <<
m_Params->
get_int_param(
"disable_photon_sim")
1020 <<
" and photon simulation is disabled in DIRC!" << std::endl;
1025 QuartzMPT->
AddProperty(
"RINDEX", PhotonEnergy, QuartzRefractiveIndex, num);
1026 QuartzMPT->
AddProperty(
"ABSLENGTH", PhotonEnergy, QuartzAbsorption, num);
1033 AirMPT->
AddProperty(
"RINDEX", PhotonEnergy, AirRefractiveIndex, num);
1034 AirMPT->
AddProperty(
"ABSLENGTH", PhotonEnergy, AirAbsorption, num);
1040 KamLandOilMPT->
AddProperty(
"RINDEX", PhotonEnergy, KamLandOilRefractiveIndex, num);
1041 KamLandOilMPT->
AddProperty(
"ABSLENGTH", PhotonEnergy, KamLandOilAbsorption, num);
1047 Nlak33aMPT->
AddProperty(
"RINDEX", PhotonEnergyNlak33a, Nlak33aRefractiveIndex, 76);
1048 Nlak33aMPT->
AddProperty(
"ABSLENGTH", PhotonEnergyNlak33a, Nlak33aAbsorption, 76);
1053 PbF2MPT->
AddProperty(
"RINDEX", en_PbF2, ref_PbF2, n_PbF2);
1054 PbF2MPT->
AddProperty(
"ABSLENGTH", en_PbF2, ab_PbF2, n_PbF2);
1059 SapphireMPT->
AddProperty(
"RINDEX", en_Sapphire, ref_Sapphire, n_Sapphire);
1060 SapphireMPT->
AddProperty(
"ABSLENGTH", en_Sapphire, ab_Sapphire, n_Sapphire);
1065 EpotekMPT->
AddProperty(
"RINDEX", PhotonEnergy, EpotekRefractiveIndex, num);
1066 EpotekMPT->
AddProperty(
"ABSLENGTH", PhotonEnergy, EpotekAbsorption, num);
1162 double QuantumEfficiencyIdial[
num] =
1163 {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1164 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1165 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1166 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
1169 double QuantumEfficiencyB[
num] =
1170 {0., 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, 0.06,
1171 0.07, 0.09, 0.1, 0.13, 0.15, 0.17, 0.2, 0.24, 0.26, 0.28, 0.282, 0.284, 0.286,
1172 0.288, 0.29, 0.28, 0.26, 0.24, 0.22, 0.20, 0.18, 0.15, 0.13, 0.12, 0.10};
1175 double QuantumEfficiencyPMT[
num] =
1176 {0.001, 0.002, 0.004, 0.007, 0.011, 0.015, 0.020, 0.026, 0.033, 0.040, 0.045,
1177 0.056, 0.067, 0.085, 0.109, 0.129, 0.138, 0.147, 0.158, 0.170,
1178 0.181, 0.188, 0.196, 0.203, 0.206, 0.212, 0.218, 0.219, 0.225, 0.230,
1179 0.228, 0.222, 0.217, 0.210, 0.199, 0.177};
1190 std::cout <<
"EIC Dirc Detector:" << std::endl;
1191 if (what ==
"ALL" || what ==
"VOLUME")
1193 std::cout <<
"Version 0.1" << std::endl;
1194 std::cout <<
"Parameters:" << std::endl;