ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EicDircDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EicDircDetector.cc
1 #include "G4EicDircDetector.h"
2 
4 
5 #include <phparameter/PHParameters.h>
6 
7 #include <g4main/PHG4Detector.h> // for PHG4Detector
8 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
9 #include <g4main/PHG4Subsystem.h>
10 
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>
38 
39 #include <cmath>
40 #include <iostream> // for operator<<, endl, bas...
41 
42 class G4VSolid;
43 class PHCompositeNode;
44 
46  PHCompositeNode* Node,
48  const std::string& dnam)
49  : PHG4Detector(subsys, Node, dnam)
50  , m_Params(parameters)
51  , m_DisplayAction(dynamic_cast<G4EicDircDisplayAction*>(subsys->GetDisplayAction()))
52 {
53 }
54 
55 //_______________________________________________________________
56 
58 {
59  /*std::map<G4VPhysicalVolume *, int>::const_iterator iter = m_PhysicalVolumes_active.find(volume);
60  if(iter != m_PhysicalVolumes_active.end())
61  {
62  return iter->second;
63  }*/
64  if (!volume) return 0;
65  std::map<G4LogicalVolume*, int>::const_iterator iter = m_LogicalVolumes_active.find(volume->GetLogicalVolume());
66 
67  if (iter != m_LogicalVolumes_active.end())
68  {
69  return iter->second;
70  }
71 
72  return 0;
73 }
74 
76 {
77  // ---------- DIRC supoort stucture ---------------
78 
79  //G4Material *Air = GetDetectorMaterial("G4_AIR");
80 
81  // positions
82  /*G4double rMin = m_Params->get_double_param("rMin"); // center location of Al support plate
83  G4double det_height = 2.1 * cm;
84  G4double place_z = m_Params->get_double_param("place_z");
85  G4double detlength = m_Params->get_double_param("length");
86 
87  G4VSolid *dirc_envelope_solid = new G4Cons("dirc_envelope_solid",
88  rMin - det_height / 2 - 2 * cm, rMin + det_height / 2 + 2 * cm,
89  rMin - det_height / 2 - 2 * cm, rMin + det_height / 2 + 2 * cm,
90  detlength / 2.0,
91  0, 2 * M_PI);
92 
93  DetectorLog_Det = new G4LogicalVolume(dirc_envelope_solid, Air, name_base + "_Log");
94 
95  G4VPhysicalVolume* wDetectorLog_Det = new G4PVPlacement(0, G4ThreeVector(0, 0, place_z), DetectorLog_Det, name_base + "_Physical", logicWorld, false, 0, overlapcheck_sector); // FullEnvelope
96  m_PhysicalVolumes_active[wDetectorLog_Det] = 20;
97 
98  // Single module with length based on readout (contains 14 LGADs [counting across both sides] in x-direction and 6 in z-direction)
99  G4double baseplate_length = 43.1 * mm;
100  //G4double baseplate_width = 56.5 * mm / 2;
101  G4double segmentlength = 6 * baseplate_length; //(detlength - 10 * cm) / 6;//m_Params->get_double_param("length");
102 
103  G4VSolid *sol_module_envelope = new G4Trd("sol_module_envelope",
104  sin(M_PI / 12.) * rMin, sin(M_PI / 12.) * (rMin + det_height),
105  segmentlength / 2, segmentlength / 2,
106  det_height / 2);
107 
108  log_module_envelope = new G4LogicalVolume(sol_module_envelope, Air, "log_module_envelope");
109 
110  G4double cooling_plate_height = 6.35 * mm;
111  G4double support_height = 7 * cm;
112 
113  // G4Material* mat_carbonfiber = new G4Material("CarbonFiberSupport", 1.44 * g / cm3, 1);
114  // mat_carbonfiber->AddElement(G4Element::GetElement("C"), 1);
115  // G4double density; //z=mean number of protons;
116  // G4int ncomponents;
117  // carbon+epoxy material
118  // G4Material *cfrp_intt = new G4Material("CFRP_INTT", density = 1.69 * g / cm3, ncomponents = 3);
119  // cfrp_intt->AddElement(G4Element::GetElement("C"), 10);
120  // cfrp_intt->AddElement(G4Element::GetElement("H"), 6);
121  // cfrp_intt->AddElement(G4Element::GetElement("O"), 1);
122 
123  // SUPPORT STRUCTURES
124  G4double support_width = 1 * mm;
125  // build components of single segment here
126  G4VSolid *Sol_End_Support = new G4Trd("Sol_End_Support",
127  sin(M_PI / 12.) * (rMin - support_height * 0.9) - 2 * mm, sin(M_PI / 12.) * (rMin) -4 * mm, // x1, x2
128  support_width / 2, support_width / 2, // length
129  support_height * 0.73 / 2); // height
130 
131  Log_End_Support = new G4LogicalVolume(Sol_End_Support, GetDetectorMaterial("G4_Fe"), "Log_End_Support_Raw");
132 
133  // place End side and back side support structure for the segment
134  // RegisterPhysicalVolume( new G4PVPlacement(0, G4ThreeVector(0, segmentlength/2-support_width/2, -support_height/2), Log_End_Support,
135  // "Front_Support_Physical", log_module_envelope, false, 0, overlapcheck_sector), false);
136  // RegisterPhysicalVolume( new G4PVPlacement(0, G4ThreeVector(0, -segmentlength/2+support_width/2, -support_height/2), Log_End_Support,
137  // "Back_Support_Physical", log_module_envelope, false, 0, overlapcheck_sector), false);
138 
139 
140  // place longitudinal supports left, middle and right side of sector
141  G4VSolid *Sol_Longitudinal_Support = new G4Trd("Sol_Longitudinal_Support",
142  support_width / 2, support_width / 2, // x1, x2
143  segmentlength / 2 - 1 * mm, segmentlength / 2 - 1 * mm, // length
144  support_height * 0.73 / 2); // height
145 
146  Log_Longitudinal_Support = new G4LogicalVolume(Sol_Longitudinal_Support, GetDetectorMaterial("G4_Fe"), "Log_Longitudinal_Support_Raw");
147 
148  G4RotationMatrix *supportrot = new G4RotationMatrix();
149  supportrot->rotateY(-M_PI / 12.);
150  if (rMin < 85 * cm)
151  {
152  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Left = new G4PVPlacement(supportrot, G4ThreeVector(sin(M_PI / 12.) * (rMin - support_height / 2) - support_width / 2, 0, -support_height / 2), Log_Longitudinal_Support, "Mother_Segment_Raw_Physical_Left", log_module_envelope, false, 0, overlapcheck_sector); // Support
153  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Left] = 21;
154 
155  G4RotationMatrix *supportrot2 = new G4RotationMatrix();
156  supportrot2->rotateY(M_PI / 12.);
157 
158  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Right = new G4PVPlacement(supportrot2, G4ThreeVector(-sin(M_PI / 12.) * (rMin - support_height / 2) + support_width / 2, 0, -support_height / 2), Log_Longitudinal_Support, "Mother_Segment_Raw_Physical_Right", log_module_envelope, false, 0, overlapcheck_sector);
159  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Right] = 22;
160  }
161 
162  G4double modulesep = 1 * mm;
163  G4double moduleShift = -8 * mm;
164  if (rMin < 85 * cm) moduleShift = -3 * mm;
165  if (rMin < 66 * cm) moduleShift = -1 * mm;
166  if (rMin < 55 * cm) moduleShift = 4 * mm;
167 
168  for (int isec = 0; isec < 12; isec++)
169  {
170  // if(isec!=3 && isec!=4)continue; // NOTE REMOVE
171  // if(isec!=3)continue; // NOTE REMOVE
172  G4RotationMatrix *motherrot = new G4RotationMatrix();
173  motherrot->rotateX(M_PI / 2);
174  motherrot->rotateY((isec - 3) * 2 * M_PI / 12.);
175  // // central segments
176  G4VPhysicalVolume* wlog_module_envelope = new G4PVPlacement(motherrot, G4ThreeVector((rMin - det_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin - det_height / 2 + moduleShift) * sin(isec * 2 * M_PI / 12.), 0 * modulesep), log_module_envelope, "Mother_Segment_Raw_Physical_Center_" + std::to_string(isec), DetectorLog_Det, false, 0, overlapcheck_sector); // ModuleEnvelope
177 
178  m_PhysicalVolumes_active[wlog_module_envelope] = 23;
179 
180  for (int ilen = 1; ilen < ((detlength / 2 - segmentlength / 2) / segmentlength); ilen++)
181  {
182  G4RotationMatrix *supfinalrot = new G4RotationMatrix();
183  // supfinalrot->rotateX(M_PI/2);
184  supfinalrot->rotateX(M_PI / 2);
185  supfinalrot->rotateY((isec - 3) * 2 * M_PI / 12.);
186  if (ilen == 2 || (ilen == 7))
187  {
188  if (rMin < 85 * cm)
189  {
190  G4VPhysicalVolume* wFront_Support_Physical_1 = new G4PVPlacement(supfinalrot, G4ThreeVector((rMin - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * sin(isec* 2 * M_PI / 12.), ilen * segmentlength + segmentlength / 2 + ilen * modulesep), Log_End_Support, "Front_Support_Physical_1_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
191  m_PhysicalVolumes_active[wFront_Support_Physical_1] = 24;
192 
193  G4VPhysicalVolume* wFront_Support_Physical_2 = new G4PVPlacement(supfinalrot, G4ThreeVector((rMin - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * sin(isec* 2 * M_PI / 12.), -(ilen * segmentlength + segmentlength / 2 + ilen * modulesep)), Log_End_Support, "Front_Support_Physical_2_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
194 
195  m_PhysicalVolumes_active[wFront_Support_Physical_2] = 25;
196 
197  }
198  }
199 
200  // forward segments
201  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Fwd = new G4PVPlacement(motherrot, G4ThreeVector((rMin - det_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin - det_height / 2 + moduleShift) * sin(isec * 2 * M_PI / 12.), ilen * segmentlength + ilen * modulesep), log_module_envelope, "Mother_Segment_Raw_Physical_Fwd_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
202 
203  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Fwd] = 26;
204 
205 
206  // backward segments
207  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Bwd = new G4PVPlacement(motherrot, G4ThreeVector((rMin - det_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin - det_height / 2 + moduleShift) * sin(isec * 2 * M_PI / 12.), -ilen * segmentlength - ilen * modulesep), log_module_envelope, "Mother_Segment_Raw_Physical_Bwd_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
208 
209  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Bwd] = 27;
210 
211  }
212  }
213 
214  // -------- INNER FRAME -----------------
215 
216  G4double rMin_inner = m_Params->get_double_param("rMin_inner"); // center location of Al support plate
217 
218  G4VSolid *sol_module_envelope_inner = new G4Trd("sol_module_envelope_inner", sin(M_PI / 12.) * rMin_inner, sin(M_PI / 12.) * (rMin_inner + det_height), segmentlength / 2, segmentlength / 2, det_height / 2);
219 
220  log_module_envelope_inner = new G4LogicalVolume(sol_module_envelope_inner, Air, "log_module_envelope_inner");
221 
222  // SUPPORT STRUCTURES
223 
224  // build components of single segment here
225  G4VSolid *Sol_End_Support_inner = new G4Trd("Sol_End_Support_inner", sin(M_PI / 12.) * (rMin_inner - support_height * 0.9) - 2 * mm, sin(M_PI / 12.) * (rMin_inner) -4 * mm, support_width / 2, support_width / 2, support_height * 0.73 / 2);
226 
227  Log_End_Support_inner = new G4LogicalVolume(Sol_End_Support_inner, GetDetectorMaterial("G4_Fe"), "Log_End_Support_Raw_inner");
228 
229  if (rMin_inner < 85 * cm)
230  {
231  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Left_inner = new G4PVPlacement(supportrot, G4ThreeVector(sin(M_PI / 12.) * (rMin_inner - support_height / 2) - support_width / 2, 0, -support_height / 2), Log_Longitudinal_Support, "Mother_Segment_Raw_Physical_Left_inner", log_module_envelope_inner, false, 0, overlapcheck_sector); // Support
232  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Left_inner] = 28;
233 
234  G4RotationMatrix *supportrot2 = new G4RotationMatrix();
235  supportrot2->rotateY(M_PI / 12.);
236 
237  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Right_inner = new G4PVPlacement(supportrot2, G4ThreeVector(-sin(M_PI / 12.) * (rMin_inner - support_height / 2) + support_width / 2, 0, -support_height / 2), Log_Longitudinal_Support, "Mother_Segment_Raw_Physical_Right_inner", log_module_envelope_inner, false, 0, overlapcheck_sector);
238  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Right_inner] = 29;
239  }
240 
241  //G4double modulesep = 1 * mm;
242  //G4double moduleShift = -8 * mm;
243  if (rMin_inner < 85 * cm) moduleShift = -3 * mm;
244  //if (rMin < 66 * cm) moduleShift = -1 * mm;
245  //if (rMin < 55 * cm) moduleShift = 4 * mm;
246 
247  for (int isec = 0; isec < 12; isec++)
248  {
249  // if(isec!=3 && isec!=4)continue; // NOTE REMOVE
250  // if(isec!=3)continue; // NOTE REMOVE
251  G4RotationMatrix *motherrot = new G4RotationMatrix();
252  motherrot->rotateX(M_PI / 2);
253  motherrot->rotateY((isec - 3) * 2 * M_PI / 12.);
254  // // central segments
255  G4VPhysicalVolume* wlog_module_envelope_inner = new G4PVPlacement(motherrot, G4ThreeVector((rMin_inner - det_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin_inner - det_height / 2 + moduleShift) * sin(isec * 2 * M_PI / 12.), 0 * modulesep), log_module_envelope_inner, "Mother_Segment_Raw_Physical_Center_inner_" + std::to_string(isec), DetectorLog_Det, false, 0, overlapcheck_sector); // ModuleEnvelope
256 
257  m_PhysicalVolumes_active[wlog_module_envelope_inner] = 30;
258 
259  for (int ilen = 1; ilen < ((detlength / 2 - segmentlength / 2) / segmentlength); ilen++)
260  {
261  G4RotationMatrix *supfinalrot = new G4RotationMatrix();
262  // supfinalrot->rotateX(M_PI/2);
263  supfinalrot->rotateX(M_PI / 2);
264  supfinalrot->rotateY((isec - 3) * 2 * M_PI / 12.);
265  if (ilen == 2 || (ilen == 7))
266  {
267  if (rMin_inner < 85 * cm)
268  {
269  G4VPhysicalVolume* wFront_Support_Physical_inner_1 = new G4PVPlacement(supfinalrot, G4ThreeVector((rMin_inner - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin_inner - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * sin(isec* 2 * M_PI / 12.), ilen * segmentlength + segmentlength / 2 + ilen * modulesep), Log_End_Support_inner, "Front_Support_Physical_inner_1_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
270  m_PhysicalVolumes_active[wFront_Support_Physical_inner_1] = 31;
271 
272  G4VPhysicalVolume* wFront_Support_Physical_inner_2 = new G4PVPlacement(supfinalrot, G4ThreeVector((rMin_inner - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin_inner - support_height / 2 - det_height / 2 - cooling_plate_height / 2 + moduleShift) * sin(isec* 2 * M_PI / 12.), -(ilen * segmentlength + segmentlength / 2 + ilen * modulesep)), Log_End_Support_inner, "Front_Support_Physical_inner_2_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
273 
274  m_PhysicalVolumes_active[wFront_Support_Physical_inner_2] = 32;
275 
276  }
277  }
278 
279  // forward segments
280  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Fwd_inner = new G4PVPlacement(motherrot, G4ThreeVector((rMin_inner - det_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin_inner - det_height / 2 + moduleShift) * sin(isec * 2 * M_PI / 12.), ilen * segmentlength + ilen * modulesep), log_module_envelope_inner, "Mother_Segment_Raw_Physical_Fwd_inner_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
281 
282  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Fwd_inner] = 33;
283 
284 
285  // backward segments
286  G4VPhysicalVolume* wMother_Segment_Raw_Physical_Bwd_inner = new G4PVPlacement(motherrot, G4ThreeVector((rMin_inner - det_height / 2 + moduleShift) * cos(isec * 2 * M_PI / 12.), (rMin_inner - det_height / 2 + moduleShift) * sin(isec * 2 * M_PI / 12.), -ilen * segmentlength - ilen * modulesep), log_module_envelope_inner, "Mother_Segment_Raw_Physical_Bwd_inner_" + std::to_string(isec) + "_" + std::to_string(ilen), DetectorLog_Det, false, 0, overlapcheck_sector);
287 
288  m_PhysicalVolumes_active[wMother_Segment_Raw_Physical_Bwd_inner] = 34;
289 
290  }
291  }*/
292 
293  // -------------- DIRC ----------------------
294 
295  fGeomType = m_Params->get_int_param("Geom_type");
296  fLensId = m_Params->get_int_param("Lens_id");
297  fNBar = m_Params->get_double_param("NBars");
298 
299  if (Verbosity()) std::cout << "Nbars = " << fNBar << std::endl;
300 
301  fNRow = m_Params->get_int_param("MCP_rows");
302  fNCol = m_Params->get_int_param("MCP_columns");
303 
304  fPrizm[0] = m_Params->get_double_param("Prizm_width");
305  fPrizm[1] = m_Params->get_double_param("Prizm_length");
306  fPrizm[3] = m_Params->get_double_param("Prizm_height_at_lens");
307  fPrizm[2] = fPrizm[3] + (fPrizm[1] * tan(32 * deg));
308 
309  fBar[0] = fBarL[0] = fBarS[0] = m_Params->get_double_param("Bar_thickness");
310  fBar[1] = fBarL[1] = fBarS[1] = m_Params->get_double_param("Bar_width");
311  fBarL[2] = m_Params->get_double_param("BarL_length");
312  fBarS[2] = m_Params->get_double_param("BarS_length");
313 
314  fNBoxes = m_Params->get_int_param("NBoxes");
315  fRadius = m_Params->get_double_param("Radius");
316  zshift = m_Params->get_double_param("z_shift");
317 
318  //-------------------
319 
320  fBarsGap = 0.15;
321 
322  if (Verbosity()) std::cout << "Nbarboxes = " << fNBoxes << std::endl;
323  if (Verbosity()) std::cout << "barrel radius = " << fRadius << " mm" << std::endl;
324 
325  fMirror[0] = m_Params->get_double_param("Mirror_height");
326  fMirror[1] = fPrizm[0];
327  fMirror[2] = 1;
328 
329  fMcpTotal[0] = fMcpTotal[1] = 53 + 4;
330  fMcpTotal[2] = 1;
331  fMcpActive[0] = fMcpActive[1] = 53;
332  fMcpActive[2] = 1;
333  fLens[0] = fLens[1] = 40;
334  fLens[2] = 10;
335 
336  fBoxWidth = fPrizm[0];
337 
338  if (fGeomType == 1) fNBoxes = 1;
339 
340  fFd[0] = fBoxWidth;
341  fFd[1] = fPrizm[2];
342  fFd[2] = 1;
343 
344  if (fLensId == 0 || fLensId == 10)
345  {
346  fLens[2] = 0;
347  }
348  if (fLensId == 2)
349  {
350  fLens[0] = fPrizm[3];
351  fLens[1] = 175;
352  fLens[2] = 14.4;
353  }
354  if (fLensId == 3)
355  {
356  fLens[0] = fPrizm[3];
357  fLens[1] = fPrizm[0] / fNBar;
358  fLens[2] = 15;
359  if (fNBar == 1) fLens[1] = fPrizm[0] / 11.;
360  }
361 
362  if (fLensId == 6)
363  {
364  fLens[0] = fPrizm[3];
365  fLens[1] = fPrizm[0];
366  fLens[2] = 12;
367  }
368 
369  fPrtRot = new G4RotationMatrix();
370 
371  //-------------------------
372  DefineMaterials();
373 
374  // ------------- Volumes --------------
375 
376  double gluethickness = 0.05;
377  int nparts = m_Params->get_int_param("Bar_pieces");
378 
379  double dirclength = fBarL[2] * (nparts - 1) + fBarS[2] + gluethickness * nparts;
380 
381  // The DIRC
382  //G4Box* gDirc = new G4Box("gDirc",205.962, 193.251, 0.5*dirclength+314.565);
383  //lDirc = new G4LogicalVolume(gDirc,defaultMaterial,"lDirc",0,0,0);
384 
385  G4Box* gFd = new G4Box("gFd", 0.5 * fFd[1], 0.5 * fFd[0], 0.5 * fFd[2]);
386  lFd = new G4LogicalVolume(gFd, defaultMaterial, "lFd", 0, 0, 0);
388 
389  double dphi = 360 * deg / (double) fNBoxes;
390  //G4VPhysicalVolume* phy;
391 
392  G4AssemblyVolume* assemblySector = new G4AssemblyVolume();
393 
394  /*for(int i=0; i<fNBoxes; i++){
395  double tphi = dphi*i;
396  double dx = fRadius * cos(tphi);
397  double dy = fRadius * sin(tphi);
398 
399  G4RotationMatrix *tRot = new G4RotationMatrix();
400  tRot->rotateZ(-tphi);
401  G4ThreeVector g4vec_sector_pos(dx, dy, zshift);
402  //phy = new G4PVPlacement(tRot,G4ThreeVector(dx,dy,zshift),lDirc,"wDirc",logicWorld,false,i,OverlapCheck());
403  assemblySector->MakeImprint(logicWorld, g4vec_sector_pos, tRot, i, OverlapCheck());
404  //m_PhysicalVolumes_active[phy] = 1;
405  }*/
406 
407  // The Bar
408  G4Box* gBarL = new G4Box("gBarL", fBarL[0] / 2., fBarL[1] / 2., fBarL[2] / 2.);
409  lBarL = new G4LogicalVolume(gBarL, BarMaterial, "lBarL", 0, 0, 0);
411 
412  G4Box* gBarS = new G4Box("gBarS", fBarS[0] / 2., fBarS[1] / 2., fBarS[2] / 2.);
413  lBarS = new G4LogicalVolume(gBarS, BarMaterial, "lBarS", 0, 0, 0);
415 
416  // Glue
417  G4Box* gGlue = new G4Box("gGlue", fBar[0] / 2., fBar[1] / 2., 0.5 * gluethickness);
418  lGlue = new G4LogicalVolume(gGlue, epotekMaterial, "lGlue", 0, 0, 0);
420 
421  int id = 0;
422 
423  for (int i = 0; i < fNBar; i++)
424  {
425  double shifty = i * (fBar[1] + fBarsGap) - 0.5 * fBoxWidth + fBar[1] / 2.;
426  for (int j = 0; j < nparts; j++)
427  {
428  if (j < (nparts - 1))
429  {
430  double z = 0.5 * dirclength - 0.5 * fBarL[2] - (fBarL[2] + gluethickness) * j;
431  //pDirc[i] = new G4PVPlacement(0,G4ThreeVector(0,shifty,z),lBarL,"wBar",lDirc,false,id,OverlapCheck());
432  G4ThreeVector g4vec_lBarL_pos(0, shifty, z);
433  assemblySector->AddPlacedVolume(lBarL, g4vec_lBarL_pos, 0);
434  //m_PhysicalVolumes_active[pDirc[i]] = 3;
435  //wGlue = new G4PVPlacement(0,G4ThreeVector(0,shifty,z-0.5*(fBarL[2]+gluethickness)),lGlue,"wGlue", lDirc,false,id,OverlapCheck());
436  G4ThreeVector g4vec_lGlue_pos(0, shifty, z - 0.5 * (fBarL[2] + gluethickness));
437  assemblySector->AddPlacedVolume(lGlue, g4vec_lGlue_pos, 0);
438  //m_PhysicalVolumes_active[wGlue] = 4;
439  }
440  else
441  {
442  double z = 0.5 * dirclength - (nparts - 1) * fBarL[2] - 0.5 * fBarS[2] - gluethickness * j;
443  //pDirc[i] = new G4PVPlacement(0,G4ThreeVector(0,shifty,z),lBarS,"wBar",lDirc,false,id,OverlapCheck());
444  G4ThreeVector g4vec_lBarS_pos(0, shifty, z);
445  assemblySector->AddPlacedVolume(lBarS, g4vec_lBarS_pos, 0);
446  //m_PhysicalVolumes_active[pDirc[i]] = 3;
447  //wGlue = new G4PVPlacement(0,G4ThreeVector(0,shifty,z-0.5*(fBarS[2]+gluethickness)),lGlue,"wGlue", lDirc,false,id,OverlapCheck());
448  G4ThreeVector g4vec_lGlue_pos(0, shifty, z - 0.5 * (fBarS[2] + gluethickness));
449  assemblySector->AddPlacedVolume(lGlue, g4vec_lGlue_pos, 0);
450  //m_PhysicalVolumes_active[wGlue] = 4;
451  }
452  id++;
453  }
454  }
455 
456  // The Mirror
457  G4Box* gMirror = new G4Box("gMirror", fMirror[0] / 2., fMirror[1] / 2., fMirror[2] / 2.);
458  lMirror = new G4LogicalVolume(gMirror, MirrorMaterial, "lMirror", 0, 0, 0);
460  //wMirror =new G4PVPlacement(0,G4ThreeVector(0,0,0.5*dirclength+fMirror[2]/2.),lMirror,"wMirror", lDirc,false,0,OverlapCheck());
461  G4ThreeVector g4vec_lMirror_pos(0, 0, 0.5 * dirclength + fMirror[2] / 2.);
462  assemblySector->AddPlacedVolume(lMirror, g4vec_lMirror_pos, 0);
463  //m_PhysicalVolumes_active[wMirror] = 5;
464 
465  // The Lenses
466  if (fLensId == 2)
467  { // 2-layer lens
468  double lensrad = 70;
469  double lensMinThikness = 2;
470  G4Box* gfbox = new G4Box("Fbox", fLens[0] / 2., fLens[1] / 2., fLens[2] / 2.);
471  G4ThreeVector zTrans(0, 0, -lensrad + fLens[2] / 2. - lensMinThikness);
472 
473  G4Sphere* gsphere = new G4Sphere("Sphere", 0, 70, 0. * deg, 360. * deg, 0. * deg, 380. * deg);
474  G4IntersectionSolid* gLens1 = new G4IntersectionSolid("Fbox*Sphere", gfbox, gsphere, new G4RotationMatrix(), zTrans);
475  G4SubtractionSolid* gLens2 = new G4SubtractionSolid("Fbox-Sphere", gfbox, gsphere, new G4RotationMatrix(), zTrans);
476 
477  lLens1 = new G4LogicalVolume(gLens1, Nlak33aMaterial, "lLens1", 0, 0, 0); //Nlak33aMaterial
479  lLens2 = new G4LogicalVolume(gLens2, BarMaterial, "lLens2", 0, 0, 0);
481  }
482 
483  if (fLensId == 3)
484  { // 3-component spherical lens
485  double lensMinThikness = 2;
486 
487  double r1 = 47.8; //0;
488  double r2 = 29.1; //0;
489 
490  //r1 = (r1==0)? 47.8: r1;
491  //r2 = (r2==0)? 29.1: r2;
492  // r1=80;
493  // r2=35;
494 
495  double cr2 = sqrt(fLens[1] * fLens[1] / 4. + fBar[0] * fBar[0] / 4.);
496  if (cr2 > r2)
497  {
498  if (Verbosity()) std::cout << "bad lens" << std::endl;
499  cr2 = r2;
500  }
501  fLens[2] = (2 * lensMinThikness + r2 - sqrt(r2 * r2 - cr2 * cr2) + lensMinThikness);
502  if (Verbosity()) std::cout << "lens thickness =" << fLens[2] << " mm" << std::endl;
503 
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);
506 
507  G4Box* gfbox = new G4Box("Fbox", fLens[0] / 2., fLens[1] / 2., fLens[2] / 2.);
508  G4Tubs* gfstube = new G4Tubs("ftube", 0, cr2, fLens[2] / 2., 0, 360 * deg);
509 
510  G4Sphere* gsphere1 = new G4Sphere("Sphere1", 0, r1, 0, 360 * deg, 0, 360 * deg);
511  G4Sphere* gsphere2 = new G4Sphere("Sphere2", 0, r2, 0, 360 * deg, 0, 360 * deg);
512 
513  G4IntersectionSolid* gbbox = new G4IntersectionSolid("bbox", gfbox, gfbox, new G4RotationMatrix(), G4ThreeVector(0, 0, -lensMinThikness * 2));
514  G4IntersectionSolid* gsbox = new G4IntersectionSolid("sbox", gfstube, gfbox, new G4RotationMatrix(), G4ThreeVector(0, 0, lensMinThikness * 2));
515 
516  G4UnionSolid* gubox = new G4UnionSolid("unionbox", gbbox, gsbox, new G4RotationMatrix(), G4ThreeVector(0, 0, 0));
517 
518  G4IntersectionSolid* gLens1 = new G4IntersectionSolid("Lens1", gubox, gsphere1, new G4RotationMatrix(), -zTrans1);
519  G4SubtractionSolid* gLenst = new G4SubtractionSolid("temp", gubox, gsphere1, new G4RotationMatrix(), -zTrans1);
520 
521  G4IntersectionSolid* gLens2 = new G4IntersectionSolid("Lens2", gLenst, gsphere2, new G4RotationMatrix(), -zTrans2);
522  G4SubtractionSolid* gLens3 = new G4SubtractionSolid("Lens3", gLenst, gsphere2, new G4RotationMatrix(), -zTrans2);
523 
524  lLens1 = new G4LogicalVolume(gLens1, BarMaterial, "lLens1", 0, 0, 0);
526  lLens2 = new G4LogicalVolume(gLens2, Nlak33aMaterial, "lLens2", 0, 0, 0); //Nlak33aMaterial //PbF2Material //SapphireMaterial
528  lLens3 = new G4LogicalVolume(gLens3, BarMaterial, "lLens3", 0, 0, 0);
530  }
531 
532  if (fLensId == 6)
533  { // 3-component cylindrical lens
534  double lensMinThikness = 2.0;
535 
536  double r1 = 33; //0;
537  double r2 = 24; //0;
538 
539  lensMinThikness = 2;
540  double layer12 = lensMinThikness * 2;
541 
542  // r1 = (r1==0)? 27.45: r1;
543  // r2 = (r2==0)? 20.02: r2;
544 
545  //r1 = (r1==0)? 33: r1;
546  //r2 = (r2==0)? 24: r2;
547  double shight = 25;
548 
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);
551 
552  G4Box* gfbox = new G4Box("fbox", 0.5 * fLens[0], 0.5 * fLens[1], 0.5 * fLens[2]);
553  G4Box* gcbox = new G4Box("cbox", 0.5 * fLens[0], 0.5 * fLens[1] + 1, 0.5 * fLens[2]);
554  G4ThreeVector tTrans1(0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
555  G4ThreeVector tTrans0(-0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
556  G4SubtractionSolid* tubox = new G4SubtractionSolid("tubox", gfbox, gcbox, new G4RotationMatrix(), tTrans1);
557  G4SubtractionSolid* gubox = new G4SubtractionSolid("gubox", tubox, gcbox, new G4RotationMatrix(), tTrans0);
558 
559  G4Tubs* gcylinder1 = new G4Tubs("Cylinder1", 0, r1, 0.5 * fLens[1], 0 * deg, 360 * deg);
560  G4Tubs* gcylinder2 = new G4Tubs("Cylinder2", 0, r2, 0.5 * fLens[1] - 0.5, 0 * deg, 360 * deg);
561  G4Tubs* gcylinder1c = new G4Tubs("Cylinder1c", 0, r1, 0.5 * fLens[1] + 0.5, 0 * deg, 360 * deg);
562  G4Tubs* gcylinder2c = new G4Tubs("Cylinder2c", 0, r2, 0.5 * fLens[1] + 0.5, 0 * deg, 360 * deg);
563  G4RotationMatrix* xRot = new G4RotationMatrix();
564  xRot->rotateX(M_PI / 2. * rad);
565 
566  G4IntersectionSolid* gLens1 = new G4IntersectionSolid("Lens1", gubox, gcylinder1, xRot, zTrans1);
567  G4SubtractionSolid* gLenst = new G4SubtractionSolid("temp", gubox, gcylinder1c, xRot, zTrans1);
568 
569  G4IntersectionSolid* gLens2 = new G4IntersectionSolid("Lens2", gLenst, gcylinder2, xRot, zTrans2);
570  G4SubtractionSolid* gLens3 = new G4SubtractionSolid("Lens3", gLenst, gcylinder2c, xRot, zTrans2);
571 
572  lLens1 = new G4LogicalVolume(gLens1, BarMaterial, "lLens1", 0, 0, 0);
574  lLens2 = new G4LogicalVolume(gLens2, Nlak33aMaterial, "lLens2", 0, 0, 0);
576  lLens3 = new G4LogicalVolume(gLens3, BarMaterial, "lLens3", 0, 0, 0);
578  }
579 
580  if (fLensId == 100)
581  {
582  fLens[2] = 200;
583  G4Box* gLens3 = new G4Box("gLens1", fBar[0] / 2., 0.5 * fBoxWidth, 100);
584  lLens3 = new G4LogicalVolume(gLens3, BarMaterial, "lLens3", 0, 0, 0);
586  }
587 
588  if (fLensId != 0 && fLensId != 10)
589  {
590  double shifth = -0.5 * (dirclength + fLens[2]);
591  if (fLensId == 100)
592  {
593  //phy = new G4PVPlacement(0,G4ThreeVector(0,0,shifth),lLens3,"wLens3", lDirc,false,0,OverlapCheck());
594  G4ThreeVector g4vec_lLens3_pos(0, 0, shifth);
595  assemblySector->AddPlacedVolume(lLens3, g4vec_lLens3_pos, 0);
596  }
597  else if (fNBar == 1 && fLensId == 3)
598  {
599  for (int i = 0; i < 11; i++)
600  {
601  double shifty = i * fLens[1] - 0.5 * (fBoxWidth - fLens[1]);
602  //phy = new G4PVPlacement(0,G4ThreeVector(0,shifty,shifth),lLens1,"wLens1", lDirc,false,i,OverlapCheck());
603  G4ThreeVector g4vec_lLens_pos(0, shifty, shifth);
604  assemblySector->AddPlacedVolume(lLens1, g4vec_lLens_pos, 0);
605  //phy = new G4PVPlacement(0,G4ThreeVector(0,shifty,shifth),lLens2,"wLens2", lDirc,false,i,OverlapCheck());
606  assemblySector->AddPlacedVolume(lLens2, g4vec_lLens_pos, 0);
607  //phy = new G4PVPlacement(0,G4ThreeVector(0,shifty,shifth),lLens3,"wLens3", lDirc,false,i,OverlapCheck());
608  assemblySector->AddPlacedVolume(lLens3, g4vec_lLens_pos, 0);
609  }
610  }
611  else
612  {
613  for (int i = 0; i < fNBar; i++)
614  {
615  double shifty = i * fLens[1] - 0.5 * (fBoxWidth - fLens[1]);
616  if (fLensId != 6)
617  {
618  //phy = new G4PVPlacement(0,G4ThreeVector(0,shifty,shifth),lLens1,"wLens1", lDirc,false,i,OverlapCheck());
619  G4ThreeVector g4vec_lLens_pos(0, shifty, shifth);
620  assemblySector->AddPlacedVolume(lLens1, g4vec_lLens_pos, 0);
621  //m_PhysicalVolumes_active[phy] = 6;
622 
623  //phy = new G4PVPlacement(0,G4ThreeVector(0,shifty,shifth),lLens2,"wLens2", lDirc,false,i,OverlapCheck());
624  assemblySector->AddPlacedVolume(lLens2, g4vec_lLens_pos, 0);
625  //m_PhysicalVolumes_active[phy] = 7;
626 
627  if (fLensId == 3)
628  {
629  //phy = new G4PVPlacement(0,G4ThreeVector(0,shifty,shifth),lLens3,"wLens3", lDirc,false,i,OverlapCheck());
630  assemblySector->AddPlacedVolume(lLens3, g4vec_lLens_pos, 0);
631  }
632  //m_PhysicalVolumes_active[phy] = 8;
633  }
634  }
635  if (fLensId == 6)
636  {
637  double sh = 0;
638  //phy = new G4PVPlacement(0,G4ThreeVector(0,0,shifth-sh),lLens1,"wLens1", lDirc,false,0,OverlapCheck());
639  G4ThreeVector g4vec_lLens_pos(0, 0, shifth - sh);
640  assemblySector->AddPlacedVolume(lLens1, g4vec_lLens_pos, 0);
641  //phy = new G4PVPlacement(0,G4ThreeVector(0,0,shifth-sh),lLens2,"wLens2", lDirc,false,0,OverlapCheck());
642  assemblySector->AddPlacedVolume(lLens2, g4vec_lLens_pos, 0);
643  //phy = new G4PVPlacement(0,G4ThreeVector(0,0,shifth-sh),lLens3,"wLens3", lDirc,false,0,OverlapCheck());
644  assemblySector->AddPlacedVolume(lLens3, g4vec_lLens_pos, 0);
645  }
646  }
647  }
648 
649  // The Prizm
650  G4Trap* gPrizm = new G4Trap("gPrizm", fPrizm[0], fPrizm[1], fPrizm[2], fPrizm[3]);
651  lPrizm = new G4LogicalVolume(gPrizm, BarMaterial, "lPrizm", 0, 0, 0);
653 
654  G4RotationMatrix* xRot = new G4RotationMatrix();
655  //xRot->rotateX(-M_PI/2.*rad); // for lDirc
656  xRot->rotateX(M_PI / 2. * rad);
657 
658  G4RotationMatrix* fdrot = new G4RotationMatrix();
659  double evshiftz = 0.5 * dirclength + fPrizm[1] + fMcpActive[2] / 2. + fLens[2];
660  double evshiftx = 0;
661 
662  fPrismShift = G4ThreeVector((fPrizm[2] + fPrizm[3]) / 4. - 0.5 * fPrizm[3], 0, -(0.5 * (dirclength + fPrizm[1]) + fLens[2]));
663  //phy = new G4PVPlacement(xRot,fPrismShift,lPrizm,"wPrizm", lDirc,false,0,OverlapCheck());
664  assemblySector->AddPlacedVolume(lPrizm, fPrismShift, xRot);
665  //m_PhysicalVolumes_active[phy] = 9;
666 
667  //phy = new G4PVPlacement(fdrot,G4ThreeVector(0.5*fFd[1]-0.5*fPrizm[3]-evshiftx,0,-evshiftz),lFd,"wFd", lDirc,false,0,OverlapCheck());
668  G4ThreeVector g4vec_lFd_pos(0.5 * fFd[1] - 0.5 * fPrizm[3] - evshiftx, 0, -evshiftz);
669  assemblySector->AddPlacedVolume(lFd, g4vec_lFd_pos, fdrot);
670  //m_PhysicalVolumes_active[phy] = 2;
671 
672  for (int i = 0; i < fNBoxes; i++)
673  {
674  double tphi = dphi * i;
675  double dx = fRadius * cos(tphi);
676  double dy = fRadius * sin(tphi);
677 
678  G4RotationMatrix* tRot = new G4RotationMatrix();
679  tRot->rotateZ(tphi);
680  G4ThreeVector g4vec_sector_pos(dx, dy, zshift);
681  //phy = new G4PVPlacement(tRot,G4ThreeVector(dx,dy,zshift),lDirc,"wDirc",logicWorld,false,i,OverlapCheck());
682  assemblySector->MakeImprint(logicWorld, g4vec_sector_pos, tRot, i, OverlapCheck());
683  //m_PhysicalVolumes_active[phy] = 1;
684  }
685 
686  // MCP --
687 
688  G4Box* gMcp = new G4Box("gMcp", fMcpTotal[0] / 2., fMcpTotal[1] / 2., fMcpTotal[2] / 2.);
689  lMcp = new G4LogicalVolume(gMcp, BarMaterial, "lMcp", 0, 0, 0);
691 
692  fNpix1 = 16;
693  fNpix2 = 16;
694 
695  if (Verbosity()) std::cout << "fNpix1=" << fNpix1 << " fNpix2=" << fNpix2 << std::endl;
696 
697  // The MCP Pixel
698  G4Box* gPixel = new G4Box("gPixel", 0.5 * fMcpActive[0] / fNpix1, 0.5 * fMcpActive[1] / fNpix2, fMcpActive[2] / 16.);
699  lPixel = new G4LogicalVolume(gPixel, BarMaterial, "lPixel", 0, 0, 0);
701 /*
702  for (int i = 0; i < fNpix2; i++)
703  {
704  for (int j = 0; j < fNpix1; j++)
705  {
706  double shiftx = i * (fMcpActive[0] / fNpix1) - fMcpActive[0] / 2. + 0.5 * fMcpActive[0] / fNpix1;
707  double shifty = j * (fMcpActive[1] / fNpix2) - fMcpActive[1] / 2. + 0.5 * fMcpActive[1] / fNpix2;
708  G4VPhysicalVolume* phy1 = new G4PVPlacement(0, G4ThreeVector(shiftx, shifty, 0), lPixel, "wPixel", lMcp, false, fNpix2 * i + j, OverlapCheck());
709  }
710  }
711 
712  int mcpId = 0;
713  double gapx = (fPrizm[2] - fNCol * fMcpTotal[0]) / (double) (fNCol + 1);
714  double gapy = (fBoxWidth - fNRow * fMcpTotal[1]) / (double) (fNRow + 1);
715  for (int i = 0; i < fNCol; i++)
716  {
717  for (int j = 0; j < fNRow; j++)
718  {
719  double shiftx = i * (fMcpTotal[0] + gapx) - 0.5 * (fFd[1] - fMcpTotal[0]) + gapx;
720  double shifty = j * (fMcpTotal[1] + gapy) - 0.5 * (fBoxWidth - fMcpTotal[1]) + gapy;
721  //phy = new G4PVPlacement(0,G4ThreeVector(shiftx,shifty,0),lMcp,"wMcp", lFd,false,mcpId++);
722  G4VPhysicalVolume* phy2 = new G4PVPlacement(0, G4ThreeVector(shiftx, shifty, 0), lMcp, "wMcp", lFd, false, mcpId, OverlapCheck());
723  //m_PhysicalVolumes_active[phy] = 10;
724  mcpId++;
725  }
726  }
727 */
728  {
729  const int num = 36;
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++)
736  {
737  WaveLength[i] = (300 + i * 10) * nanometer;
738  PhotonEnergy[num - (i + 1)] = LambdaE / WaveLength[i];
739  PMTReflectivity[i] = 0.;
740  EfficiencyMirrors[i] = 0;
741  }
742 
743  /***************** QUANTUM EFFICIENCY OF BURLE AND HAMAMTSU PMT'S *****/
744 
745  //ideal pmt quantum efficiency
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};
751 
752  // Burle PMT's
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};
757 
758  //hamamatsu pmt quantum efficiency
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};
764 
765  // these quantum efficiencies have to be multiplied by geometry
766  // efficiency of given PMT's
767  // for Hamamatsu by factor 0.7
768  // for Burle by factor 0.45
769  for (int k = 0; k < 36; k++)
770  {
771  QuantumEfficiencyB[k] = QuantumEfficiencyB[k] * 0.45;
772  QuantumEfficiencyPMT[k] = QuantumEfficiencyPMT[k] * .7;
773  }
774 
775  /* define quantum efficiency for burle PMT's - the same efficiency is
776  assigned to pads and also to slots !!!! */
777 
778  //burle pmt - bigger slots => logicPad
779  G4MaterialPropertiesTable* PhotocatodBurleMPT = new G4MaterialPropertiesTable();
780  PhotocatodBurleMPT->AddProperty("EFFICIENCY", PhotonEnergy, QuantumEfficiencyB, num);
781  PhotocatodBurleMPT->AddProperty("REFLECTIVITY", PhotonEnergy, PMTReflectivity, num);
782 
783  G4OpticalSurface* BurlePMTOpSurface =
784  new G4OpticalSurface("BurlePMTOpSurface", glisur, polished,
786  BurlePMTOpSurface->SetMaterialPropertiesTable(PhotocatodBurleMPT);
787 
788  /* hamamatsu pmt's - smaller slots => quantum efficiency again
789  assign to slot and pad */
790 
791  fQuantumEfficiency = QuantumEfficiencyIdial;
792  G4MaterialPropertiesTable* PhotocatodHamamatsuMPT = new G4MaterialPropertiesTable();
793  PhotocatodHamamatsuMPT->AddProperty("EFFICIENCY", PhotonEnergy, fQuantumEfficiency, num);
794  PhotocatodHamamatsuMPT->AddProperty("REFLECTIVITY", PhotonEnergy, PMTReflectivity, num);
795 
796  G4OpticalSurface* HamamatsuPMTOpSurface =
797  new G4OpticalSurface("HamamatsuPMTOpSurface", glisur, polished, dielectric_metal);
798  HamamatsuPMTOpSurface->SetMaterialPropertiesTable(PhotocatodHamamatsuMPT);
799 
800  // // assignment to pad
801  // if(hamamatsu8500)
802  new G4LogicalSkinSurface("HamamatsuPMTSurface", lPixel, HamamatsuPMTOpSurface);
803 
804  // Mirror
805  G4OpticalSurface* MirrorOpSurface =
806  new G4OpticalSurface("MirrorOpSurface", glisur, polished, dielectric_metal);
807 
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};
813 
815  MirrorMPT->AddProperty("REFLECTIVITY", PhotonEnergy, ReflectivityMirrorBar, num);
816  MirrorMPT->AddProperty("EFFICIENCY", PhotonEnergy, EfficiencyMirrors, num);
817 
818  MirrorOpSurface->SetMaterialPropertiesTable(MirrorMPT);
819  new G4LogicalSkinSurface("MirrorSurface", lMirror, MirrorOpSurface);
820  }
821 
823 
824  /*PrtOpBoundaryProcess *fBoundaryProcess = new PrtOpBoundaryProcess();
825  G4ProcessManager *pmanager = G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
826  pmanager->AddDiscreteProcess(fBoundaryProcess);
827  pmanager->SetProcessOrderingToFirst(fBoundaryProcess, idxPostStep);
828  */
829 }
830 
832 {
833  G4String symbol; //a=mass of a mole;
834  double a, z, density; //z=mean number of protons;
835 
836  int ncomponents, natoms;
837  double fractionmass;
838 
839  // define Elements
840  G4Element* H = new G4Element("Hydrogen", symbol = "H", z = 1., a = 1.01 * g / mole);
841  G4Element* C = new G4Element("Carbon", symbol = "C", z = 6., a = 12.01 * g / mole);
842  G4Element* N = new G4Element("Nitrogen", symbol = "N", z = 7., a = 14.01 * g / mole);
843  G4Element* O = new G4Element("Oxygen", symbol = "O", z = 8., a = 16.00 * g / mole);
844  G4Element* Si = new G4Element("Silicon", symbol = "Si", z = 14., a = 28.09 * g / mole);
845  G4Element* Al = new G4Element("Aluminum", symbol = "Al", z = 13., a = 26.98 * g / mole);
846 
847  // quartz material = SiO2
848  G4Material* SiO2 = new G4Material("quartz", density = 2.200 * g / cm3, ncomponents = 2);
849  SiO2->AddElement(Si, natoms = 1);
850  SiO2->AddElement(O, natoms = 2);
851 
852  Nlak33aMaterial = new G4Material("Nlak33a", density = 4.220 * g / cm3, ncomponents = 2);
853  Nlak33aMaterial->AddElement(Si, natoms = 1);
854  Nlak33aMaterial->AddElement(O, natoms = 2);
855 
856  PbF2Material = new G4Material("PbF2", density = 3.97 * g / cm3, ncomponents = 2);
857  PbF2Material->AddElement(Si, natoms = 1);
858  PbF2Material->AddElement(O, natoms = 2);
859 
860  SapphireMaterial = new G4Material("Sapphire", density = 4.220 * g / cm3, ncomponents = 2);
861  SapphireMaterial->AddElement(Al, natoms = 2);
862  SapphireMaterial->AddElement(O, natoms = 3);
863 
864  G4Material* Vacuum = new G4Material("interGalactic", 1., 1.008 * g / mole,
865  1.e-25 * g / cm3, kStateGas,
866  2.73 * kelvin, 3.e-18 * pascal);
867  G4Material* Air = new G4Material("Air", density = 1.290 * mg / cm3, ncomponents = 2);
868  Air->AddElement(N, fractionmass = 0.7);
869  Air->AddElement(O, fractionmass = 0.3);
870 
871  G4Material* Aluminum = new G4Material("Aluminum", density = 2.7 * g / cm3, ncomponents = 1);
872  Aluminum->AddElement(Al, fractionmass = 1.0);
873 
874  G4Material* KamLandOil = new G4Material("KamLandOil", density = 0.914 * g / cm3, ncomponents = 2);
875  KamLandOil->AddElement(C, natoms = 12);
876  KamLandOil->AddElement(H, natoms = 26);
877 
878  G4Material* CarbonFiber = new G4Material("CarbonFiber", density = 0.145 * g / cm3, ncomponents = 1);
879  CarbonFiber->AddElement(C, natoms = 1);
880 
881  /* as I don't know the exact material composition,
882  I will use Epoxyd material composition and add
883  the optical property of Epotek to this material */
884 
885  G4Material* Epotek = new G4Material("Epotek", density = 1.2 * g / cm3, ncomponents = 3);
886 
887  Epotek->AddElement(C, natoms = 3);
888  Epotek->AddElement(H, natoms = 5);
889  Epotek->AddElement(O, natoms = 2);
890 
891  // assign main materials
892  if (fGeomType < 10)
893  defaultMaterial = Vacuum;
894  else
895  defaultMaterial = Air; //Vacuum // material of world
896  frontMaterial = CarbonFiber;
897  BarMaterial = SiO2; // material of all Bars, Quartz and Window
898  OilMaterial = KamLandOil; // material of volume 1,2,3,4
899  MirrorMaterial = Aluminum; // mirror material
900  epotekMaterial = Epotek; // Epotek material - glue between bars
901 
902  // ------------ Generate & Add Material Properties Table ------------
903 
904  static const double LambdaE = 2.0 * 3.14159265358979323846 * 1.973269602e-16 * m * GeV;
905  const int num = 36;
906  double WaveLength[num];
907  //double Absorption[num]; // default value for absorption
908  double AirAbsorption[num]; // absorption value for air
909  double AirRefractiveIndex[num]; // air refractive index
910  double PhotonEnergy[num]; // energy of photons which correspond to the given
911  // refractive or absoprtion values
912 
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};
914 
915  int n_PbF2 = 56;
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};
919 
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};
923 
924  double ab_Sapphire[n_Sapphire];
925  //crystan standard grade 2mm
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;
928 
929  /*************************** ABSORPTION COEFFICIENTS *****************************/
930 
931  // absorption of KamLandOil per 50 cm - from jjv
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};
939 
940  // absorption of quartz per 1 m - from jjv
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};
949 
950  // absorption of epotek per one layer - thicknes 0.001'' - from jjv
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,
958  0.9871, 0.9745};
959 
960  //N-Lak 33a
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};
962 
963  double EpotekThickness = 0.001 * 2.54 * cm;
964  for (int i = 0; i < num; i++)
965  {
966  WaveLength[i] = (300 + i * 10) * nanometer;
967  //Absorption[i]= 100*m; // not true, just due to definiton -> not absorb any
968  AirAbsorption[i] = 4. * cm; // if photon in the air -> kill it immediately
969  AirRefractiveIndex[i] = 1.;
970  PhotonEnergy[num - (i + 1)] = LambdaE / WaveLength[i];
971 
972  /* as the absorption is given per length and G4 needs
973  mean free path length, calculate it here
974  mean free path length - taken as probability equal 1/e
975  that the photon will be absorbed */
976 
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;
980  }
981 
982  /**************************** REFRACTIVE INDEXES ****************************/
983 
984  // only phase refractive indexes are necessary -> g4 calculates group itself !!
985 
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};
992 
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};
1000 
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};
1007 
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};
1009 
1010  /* ASSIGNING REFRACTIVE AND ABSORPTION PROPERTIES TO THE GIVEN MATERIALS */
1011 
1012  if (m_Params->get_int_param("disable_photon_sim") == 0)
1013  {
1014  // option disable photon simulation if explicitly set via macro
1015  static bool once = true;
1016  if (once)
1017  {
1018  once = false;
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;
1021  }
1022 
1023  // Quartz material => Si02
1025  QuartzMPT->AddProperty("RINDEX", PhotonEnergy, QuartzRefractiveIndex, num);
1026  QuartzMPT->AddProperty("ABSLENGTH", PhotonEnergy, QuartzAbsorption, num);
1027 
1028  // assign this parameter table to BAR material
1030 
1031  // Air
1033  AirMPT->AddProperty("RINDEX", PhotonEnergy, AirRefractiveIndex, num);
1034  AirMPT->AddProperty("ABSLENGTH", PhotonEnergy, AirAbsorption, num);
1035  // assign this parameter table to the air
1037 
1038  // KamLandOil
1040  KamLandOilMPT->AddProperty("RINDEX", PhotonEnergy, KamLandOilRefractiveIndex, num);
1041  KamLandOilMPT->AddProperty("ABSLENGTH", PhotonEnergy, KamLandOilAbsorption, num);
1042  // assing this parameter table to the KamLandOil
1043  OilMaterial->SetMaterialPropertiesTable(KamLandOilMPT);
1044 
1045  // N-Lak 33a
1047  Nlak33aMPT->AddProperty("RINDEX", PhotonEnergyNlak33a, Nlak33aRefractiveIndex, 76);
1048  Nlak33aMPT->AddProperty("ABSLENGTH", PhotonEnergyNlak33a, Nlak33aAbsorption, 76);
1050 
1051  // PbF2
1053  PbF2MPT->AddProperty("RINDEX", en_PbF2, ref_PbF2, n_PbF2);
1054  PbF2MPT->AddProperty("ABSLENGTH", en_PbF2, ab_PbF2, n_PbF2);
1056 
1057  // Sapphire
1059  SapphireMPT->AddProperty("RINDEX", en_Sapphire, ref_Sapphire, n_Sapphire);
1060  SapphireMPT->AddProperty("ABSLENGTH", en_Sapphire, ab_Sapphire, n_Sapphire);
1062 
1063  // Epotek Glue
1065  EpotekMPT->AddProperty("RINDEX", PhotonEnergy, EpotekRefractiveIndex, num);
1066  EpotekMPT->AddProperty("ABSLENGTH", PhotonEnergy, EpotekAbsorption, num);
1067  // assign this parameter table to the epotek
1069  }
1070 }
1071 
1073 {
1074  /*G4VisAttributes *waModuleEnvelope = new G4VisAttributes();
1075  waModuleEnvelope->SetVisibility(false);
1076  waModuleEnvelope->SetColour(G4Colour::Red());
1077  waModuleEnvelope->SetForceWireframe(true);
1078  log_module_envelope->SetVisAttributes(waModuleEnvelope);
1079  log_module_envelope_inner->SetVisAttributes(waModuleEnvelope);
1080 
1081  G4VisAttributes *waFullEnvelope = new G4VisAttributes();
1082  waFullEnvelope->SetVisibility(false);
1083  waFullEnvelope->SetColour(G4Colour::Red());
1084  waFullEnvelope->SetForceWireframe(true);
1085  DetectorLog_Det->SetVisAttributes(waFullEnvelope);
1086 
1087  G4VisAttributes *waSupport = new G4VisAttributes();
1088  waSupport->SetColour(G4Colour::Gray());
1089  waSupport->SetVisibility(true);
1090  waSupport->SetForceSolid(true);
1091  Log_End_Support->SetVisAttributes(waSupport);
1092  Log_Longitudinal_Support->SetVisAttributes(waSupport);
1093  Log_End_Support_inner->SetVisAttributes(waSupport);
1094  */
1095  G4Colour DircColour = G4Colour(1., 1.0, 0., 1.);
1096 
1097  /*G4VisAttributes *waDirc = new G4VisAttributes(DircColour);
1098  waDirc->SetForceWireframe(true);
1099  waDirc->SetVisibility(false);
1100  lDirc->SetVisAttributes(waDirc);
1101  */
1102  G4VisAttributes* waFd = new G4VisAttributes(DircColour);
1103  waFd->SetForceWireframe(true);
1104  lFd->SetVisAttributes(waFd);
1105 
1106  G4VisAttributes* waBar = new G4VisAttributes(G4Colour(0., 153./255, 0., 1)); //0.05
1107  waBar->SetVisibility(true);
1108  lBarL->SetVisAttributes(waBar);
1109  lBarS->SetVisAttributes(waBar);
1110 
1111  G4VisAttributes* waGlue = new G4VisAttributes(G4Colour(0., 0.4, 0.9, 1));
1112  waGlue->SetVisibility(true);
1113  lGlue->SetVisAttributes(waGlue);
1114 
1115  G4VisAttributes* waMirror = new G4VisAttributes(G4Colour(1., 1., 0.9, 1));
1116  waMirror->SetVisibility(true);
1117  lMirror->SetVisAttributes(waMirror);
1118 
1119  double transp = 1;
1120  G4VisAttributes* vaLens = new G4VisAttributes(G4Colour(0., 1., 1., transp));
1121  //vaLens->SetForceWireframe(true);
1122  // vaLens->SetForceAuxEdgeVisible(true);
1123  // vaLens->SetForceLineSegmentsPerCircle(10);
1124  // vaLens->SetLineWidth(4);
1125 
1126  if (fLensId == 100) lLens3->SetVisAttributes(vaLens);
1127 
1128  if (fLensId == 2 || fLensId == 3 || fLensId == 6)
1129  {
1130  lLens1->SetVisAttributes(vaLens);
1131  G4VisAttributes* vaLens2 = new G4VisAttributes(G4Colour(0., 1., 1., transp));
1132  vaLens2->SetColour(G4Colour(0., 0.5, 1., transp));
1133  vaLens2->SetForceWireframe(true);
1134  lLens2->SetVisAttributes(vaLens2);
1135  if (fLensId == 3 || fLensId == 6) lLens3->SetVisAttributes(vaLens);
1136  }
1137 
1138  G4VisAttributes* waPrizm = new G4VisAttributes(G4Colour(0., 0.9, 0.9, 1)); //0.4
1139  // waPrizm->SetForceAuxEdgeVisible(true);
1140  // waPrizm->SetForceSolid(true);
1141  lPrizm->SetVisAttributes(waPrizm);
1142 
1143  // G4VisAttributes *waMcp = new G4VisAttributes(G4Colour(0.1,0.1,0.9,0.3));
1144  G4VisAttributes* waMcp = new G4VisAttributes(G4Colour(1.0, 0., 0.1, 1));
1145  // waMcp->SetForceWireframe(true);
1146  waMcp->SetForceSolid(true);
1147  lMcp->SetVisAttributes(waMcp);
1148 
1149  G4VisAttributes* waPixel = new G4VisAttributes(G4Colour(0.7, 0.0, 0.1, 1));
1150  waPixel->SetForceWireframe(true);
1151  if (fMcpLayout == 3)
1152  waPixel->SetVisibility(true);
1153  else
1154  waPixel->SetVisibility(true);
1155  lPixel->SetVisAttributes(waPixel);
1156 }
1157 
1159 {
1160  const int num = 36;
1161  //ideal pmt quantum efficiency
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};
1167 
1168  // Burle PMT's
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};
1173 
1174  //hamamatsu pmt quantum efficiency
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};
1180 
1181  if (id == 0) fQuantumEfficiency = QuantumEfficiencyIdial;
1182  if (id == 1) fQuantumEfficiency = QuantumEfficiencyPMT;
1183  if (id == 2) fQuantumEfficiency = QuantumEfficiencyB;
1184 
1186 }
1187 
1188 void G4EicDircDetector::Print(const std::string& what) const
1189 {
1190  std::cout << "EIC Dirc Detector:" << std::endl;
1191  if (what == "ALL" || what == "VOLUME")
1192  {
1193  std::cout << "Version 0.1" << std::endl;
1194  std::cout << "Parameters:" << std::endl;
1195  m_Params->Print();
1196  }
1197  return;
1198 }