ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ForwardDualReadoutDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ForwardDualReadoutDetector.cc
1 #include "PHG4ForwardDualReadoutDetector.h"
2 #include "PHG4ForwardDualReadoutDisplayAction.h"
3 #include "PHG4ForwardDualReadoutSteppingAction.h"
4 
5 #include <g4main/PHG4Detector.h> // for PHG4Detector
6 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
7 #include <g4main/PHG4Subsystem.h>
8 
9 #include <Geant4/G4Box.hh>
10 #include <Geant4/G4Tubs.hh>
11 #include <Geant4/G4SubtractionSolid.hh>
12 #include <Geant4/G4Cons.hh>
13 #include <Geant4/G4LogicalVolume.hh>
14 #include <Geant4/G4PVPlacement.hh>
15 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
16 #include <Geant4/G4String.hh> // for G4String
17 #include <Geant4/G4SystemOfUnits.hh>
18 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
19 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
20 #include <Geant4/G4Types.hh> // for G4double, G4int
21 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
22 #include <Geant4/G4PVParameterised.hh>
23 #include <Geant4/G4PVReplica.hh>
24 #include <Geant4/G4NistManager.hh>
25 #include <Geant4/G4OpticalSurface.hh>
26 #include <Geant4/G4LogicalSkinSurface.hh>
27 #include <Geant4/G4LogicalBorderSurface.hh>
28 
29 #include <TSystem.h>
30 
31 #include <cmath>
32 #include <cstdlib>
33 #include <fstream>
34 #include <iostream>
35 #include <sstream>
36 #include <utility> // for pair, make_pair
37 
38 class G4VSolid;
39 class PHCompositeNode;
40 
41 using namespace std;
42 
43 //_______________________________________________________________________
44 PHG4ForwardDualReadoutDetector::PHG4ForwardDualReadoutDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const std::string& dnam)
45  : PHG4Detector(subsys, Node, dnam)
46  , m_DisplayAction(dynamic_cast<PHG4ForwardDualReadoutDisplayAction*>(subsys->GetDisplayAction()))
47  , m_SteppingAction(0)
48  , _place_in_x(0.0 * mm)
49  , _place_in_y(0.0 * mm)
50  , _place_in_z(4000.0 * mm)
51  , _center_offset_x(0.0 * mm)
52  , _center_offset_y(0.0 * mm)
53  , _quadratic_detector(0)
54  , _rot_in_x(0.0)
55  , _rot_in_y(0.0)
56  , _rot_in_z(0.0)
57  , _rMin1(50 * mm)
58  , _rMax1(2620 * mm)
59  , _rMin2(50 * mm)
60  , _rMax2(3369 * mm)
61  , _dZ(1000 * mm)
62  , _sPhi(0)
63  , _dPhi(2 * M_PI)
64  , _tower_type(0)
65  , _tower_readout(0.5 * mm)
66  , _tower_dx(100 * mm)
67  , _tower_dy(100 * mm)
68  , _tower_dz(1000.0 * mm)
69  , _scintFiber_diam(1.0 * mm)
70  , _cerenkovFiber_diam(1.0 * mm)
71  , _cerenkovFiber_material(0)
72  , _tower_makeNotched(0)
73  , _absorber_Material(0)
74  , _materialScintillator("G4_POLYSTYRENE")
75  , _materialAbsorber("G4_Fe")
76  , _active(1)
77  , _absorberactive(0)
78  , _layer(0)
79  , _blackhole(0)
80  , _towerlogicnameprefix("hdrcaloTower")
81  , _superdetector("NONE")
82  , _mapping_tower_file("")
83 {
84 }
85 //_______________________________________________________________________
86 int PHG4ForwardDualReadoutDetector::IsInForwardDualReadout(G4VPhysicalVolume* volume) const
87 {
88  if (volume->GetName().find(_towerlogicnameprefix) != string::npos)
89  {
90  if (volume->GetName().find("scintillator") != string::npos)
91  {
92  if (_active)
93  return 1;
94  else
95  return 0;
96  }
97  else if (volume->GetName().find("cherenkov") != string::npos)
98  {
99  if (_active)
100  return 1;
101  else
102  return 0;
103  }
104  //only record energy in actual absorber- drop energy lost in air gaps inside drcalo envelope
105  else if (volume->GetName().find("absorber") != string::npos)
106  {
107  if (_absorberactive)
108  return -1;
109  else
110  return 0;
111  }
112  else if (volume->GetName().find("envelope") != string::npos)
113  {
114  return 0;
115  }
116  }
117 
118  return 0;
119 }
120 
121 //_______________________________________________________________________
122 void PHG4ForwardDualReadoutDetector::ConstructMe(G4LogicalVolume* logicWorld)
123 {
124  if (Verbosity() > 0)
125  {
126  cout << "PHG4ForwardDualReadoutDetector: Begin Construction" << endl;
127  }
128 
129  //Read parameters for detector construction from file
130  ParseParametersFromTable();
131 
132  //Create the cone envelope = 'world volume' for the calorimeter
133  G4Material* Air = GetDetectorMaterial("G4_AIR");
134 
135  G4VSolid* drcalo_envelope_solid;
136  if(_quadratic_detector){
137  // box with round cutout in the middle
138  G4VSolid* beampipe_cutout = new G4Cons("hdrcalo_beampipe_cutout",
139  0, _rMin1,
140  0, _rMin1,
141  _dZ / 2.0,
142  _sPhi, _dPhi);
143  drcalo_envelope_solid = new G4Box("hdrcalo_envelope_solid_precut",
144  _rMax1,
145  _rMax1,
146  _tower_dz / 2.0);
147  drcalo_envelope_solid = new G4SubtractionSolid(G4String("hdrcalo_envelope_solid"), drcalo_envelope_solid, beampipe_cutout
148  , 0 ,G4ThreeVector( 0 , 0 ,0.));
149  } else {
150  drcalo_envelope_solid = new G4Cons("hdrcalo_envelope_solid",
151  _rMin1, _rMax1,
152  _rMin2, _rMax2,
153  _dZ / 2.0,
154  _sPhi, _dPhi);
155  }
156 
157  G4LogicalVolume* drcalo_envelope_log = new G4LogicalVolume(drcalo_envelope_solid, Air, G4String("hdrcalo_envelope"), 0, 0, 0);
158 
159  m_DisplayAction->AddVolume(drcalo_envelope_log, "FdrcaloEnvelope");
160 
161  //Define rotation attributes for envelope cone
162  G4RotationMatrix drcalo_rotm;
163  drcalo_rotm.rotateX(_rot_in_x);
164  drcalo_rotm.rotateY(_rot_in_y);
165  drcalo_rotm.rotateZ(_rot_in_z);
166 
167  //Place envelope cone in simulation
168  ostringstream name_envelope;
169  name_envelope.str("");
170  name_envelope << _towerlogicnameprefix << "_envelope" << endl;
171 
172  new G4PVPlacement(G4Transform3D(drcalo_rotm, G4ThreeVector(_place_in_x, _place_in_y, _place_in_z)),
173  drcalo_envelope_log, name_envelope.str().c_str(), logicWorld, 0, false, OverlapCheck());
174 
175 
176  G4LogicalVolume* singletower;
177  if(_tower_type==5) singletower = ConstructTowerFCStyle(0);
178  else singletower = ConstructTower(0); //4x4 fibre tower with 2 scint and 2 Cherenkov fibres
179 
180 
181  G4Material* material_air = GetDetectorMaterial("G4_AIR");
182  // number of towers in radial direction (on y axis)
183  int rowNtow = (int) ( (_rMax1-(_tower_dy/2)) / _tower_dy);
184  for(int row=rowNtow;row>=-rowNtow;row--){
185  // pythagoras -> get available length in circular mother volume for towers
186  // divide given length by tower width -> get number of towers that can be placed
187  int currRowNtow = (int) ( ( 2* sqrt(pow(_rMax1,2)-pow( (abs(row)*_tower_dy) ,2)) ) / _tower_dy );
188  if(currRowNtow==0) continue;
189  // we want an odd number of towers to be symmetrically centered around 0
190  if ( currRowNtow % 2 == 0) currRowNtow-=1;
191 
192  if( ( (abs(row)*_tower_dy) ) < _rMin1 ){ // _rMin1
193  if(_center_offset_x!=0){
194  // pythagoras -> get available length in circular mother volume for towers
195  // divide given length by tower width -> get number of towers that can be placed
196  int currRowNtowInner = (int) ( ( 2* sqrt(pow(_rMin1,2)-pow( (abs(row)*_tower_dy) ,2)) ) / _tower_dy );
197  // we want an odd number of towers to be symmetrically centered around 0
198  if(_quadratic_detector){
199  if ( currRowNtowInner % 2 == 0) currRowNtowInner+=1;
200  currRowNtowInner+=1;
201  } else {
202  if ( currRowNtowInner % 2 == 0) currRowNtowInner+=1;
203  }
204  int offsetrows = (int) ( _center_offset_x / _tower_dy );
205  if ( offsetrows % 2 != 1) offsetrows-=1;
206 
207  int currRowNtowMod = currRowNtow;
208  if(_quadratic_detector){
209  currRowNtowMod = 2*rowNtow;
210  }
211 
212  // create mother volume with space for currRowNtow towers along x-axis
213  auto DRCalRowLeftSolid = new G4Box("DRCalRowLeftBox" + std::to_string(row), ((currRowNtowMod - currRowNtowInner) / 2 + offsetrows) * _tower_dx / 2.0,_tower_dy / 2.0,_tower_dz / 2.0);
214  auto DRCalRowLeftLogical = new G4LogicalVolume(DRCalRowLeftSolid,material_air,"DRCalRowLeftLogical" + std::to_string(row));
215  // replicate singletower tower design currRowNtow times along x-axis
216  new G4PVReplica("DRCalRowLeftPhysical" + std::to_string(row),singletower,DRCalRowLeftLogical,
217  kXAxis,((currRowNtowMod - currRowNtowInner) / 2 + offsetrows),_tower_dx);
218 
219  ostringstream name_row_twr_left;
220  name_row_twr_left.str("");
221  name_row_twr_left << _towerlogicnameprefix << "_row_" << row << "_left" << endl;
222  new G4PVPlacement(0, G4ThreeVector( - ( ( currRowNtowInner / 2.0 ) * _tower_dx ) - ( ((currRowNtowMod - currRowNtowInner) / 2 - offsetrows) * _tower_dx / 2.0 ), (row*_tower_dy), 0),
223  DRCalRowLeftLogical, name_row_twr_left.str().c_str(), drcalo_envelope_log, 0, false, OverlapCheck());
224 
225  // create mother volume with space for currRowNtow towers along x-axis
226  auto DRCalRowRightSolid = new G4Box("DRCalRowRightBox" + std::to_string(row), ((currRowNtowMod - currRowNtowInner) / 2 - offsetrows ) * _tower_dx / 2.0,_tower_dy / 2.0,_tower_dz / 2.0);
227  auto DRCalRowRightLogical = new G4LogicalVolume(DRCalRowRightSolid,material_air,"DRCalRowRightLogical" + std::to_string(row));
228  // replicate singletower tower design currRowNtow times along x-axis
229  new G4PVReplica("DRCalRowRightPhysical" + std::to_string(row),singletower,DRCalRowRightLogical,
230  kXAxis,((currRowNtowMod - currRowNtowInner) / 2 - offsetrows ),_tower_dx);
231 
232  ostringstream name_row_twr_right;
233  name_row_twr_right.str("");
234  name_row_twr_right << _towerlogicnameprefix << "_row_" << row << "_right" << endl;
235  new G4PVPlacement(0, G4ThreeVector( ( ( currRowNtowInner / 2.0 ) * _tower_dx ) + ( ((currRowNtowMod - currRowNtowInner) / 2 + offsetrows) * _tower_dx / 2.0 ), (row*_tower_dy), 0),
236  DRCalRowRightLogical, name_row_twr_right.str().c_str(), drcalo_envelope_log, 0, false, OverlapCheck());
237  } else {
238  // pythagoras -> get available length in circular mother volume for towers
239  // divide given length by tower width -> get number of towers that can be placed
240  int currRowNtowInner = (int) ( ( 2* sqrt(pow(_rMin1,2)-pow( (abs(row)*_tower_dy) - (_tower_dy/2.0) ,2)) ) / _tower_dy );
241  // we want an odd number of towers to be symmetrically centered around 0
242  if ( currRowNtowInner % 2 != 0) currRowNtowInner-=1;
243  // currRowNtowInner+=2;
244  int currRowNtowMod = currRowNtow;
245  if(_quadratic_detector){
246  currRowNtowMod = rowNtow;
247  }
248  // create mother volume with space for currRowNtow towers along x-axis
249  auto DRCalRowSolid = new G4Box("DRCalRowBox" + std::to_string(row), (currRowNtowMod - currRowNtowInner) / 2 * _tower_dx / 2.0,_tower_dy / 2.0,_tower_dz / 2.0);
250  auto DRCalRowLogical = new G4LogicalVolume(DRCalRowSolid,material_air,"DRCalRowLogical" + std::to_string(row));
251  // replicate singletower tower design currRowNtow times along x-axis
252  new G4PVReplica("DRCalRowPhysical" + std::to_string(row),singletower,DRCalRowLogical,
253  kXAxis,(currRowNtowMod - currRowNtowInner) / 2,_tower_dx);
254 
255  ostringstream name_row_twr;
256  name_row_twr.str("");
257  name_row_twr << _towerlogicnameprefix << "_row_" << row << "_left" << endl;
258  new G4PVPlacement(0, G4ThreeVector( - ( ( currRowNtowInner / 2.0 ) * _tower_dx ) - ( (currRowNtowMod - currRowNtowInner) / 2 * _tower_dx / 2.0 ), (row*_tower_dy), 0),
259  DRCalRowLogical, name_row_twr.str().c_str(), drcalo_envelope_log, 0, false, OverlapCheck());
260 
261  ostringstream name_row_twr2;
262  name_row_twr2.str("");
263  name_row_twr2 << _towerlogicnameprefix << "_row_" << row << "_left" << endl;
264  new G4PVPlacement(0, G4ThreeVector( ( ( currRowNtowInner / 2.0 ) * _tower_dx ) + ( (currRowNtowMod - currRowNtowInner) / 2 * _tower_dx / 2.0 ), (row*_tower_dy), 0),
265  DRCalRowLogical, name_row_twr2.str().c_str(), drcalo_envelope_log, 0, false, OverlapCheck());
266  }
267 
268  } else {
269  if(_quadratic_detector){
270  // cout << currRowNtow << endl;
271  // create mother volume with space for currRowNtow towers along x-axis
272  auto DRCalRowSolid = new G4Box("DRCalRowBox" + std::to_string(row), 2 * rowNtow * _tower_dx / 2.0,_tower_dy / 2.0,_tower_dz / 2.0);
273  auto DRCalRowLogical = new G4LogicalVolume(DRCalRowSolid,material_air,"DRCalRowLogical" + std::to_string(row));
274  // replicate singletower tower design currRowNtow times along x-axis
275  new G4PVReplica("DRCalRowPhysical" + std::to_string(row),singletower,DRCalRowLogical,
276  kXAxis,2 * rowNtow,_tower_dx);
277 
278  ostringstream name_row_twr;
279  name_row_twr.str("");
280  name_row_twr << _towerlogicnameprefix << "_row_" << row << endl;
281  new G4PVPlacement(0, G4ThreeVector(0, (row*_tower_dy), 0),
282  DRCalRowLogical, name_row_twr.str().c_str(), drcalo_envelope_log, 0, false, OverlapCheck());
283  } else {
284  // create mother volume with space for currRowNtow towers along x-axis
285  // cout << currRowNtow << endl;
286  auto DRCalRowSolid = new G4Box("DRCalRowBox" + std::to_string(row), currRowNtow * _tower_dx / 2.0,_tower_dy / 2.0,_tower_dz / 2.0);
287  auto DRCalRowLogical = new G4LogicalVolume(DRCalRowSolid,material_air,"DRCalRowLogical" + std::to_string(row));
288  // replicate singletower tower design currRowNtow times along x-axis
289  new G4PVReplica("DRCalRowPhysical" + std::to_string(row),singletower,DRCalRowLogical,
290  kXAxis,currRowNtow,_tower_dx);
291 
292  ostringstream name_row_twr;
293  name_row_twr.str("");
294  name_row_twr << _towerlogicnameprefix << "_row_" << row << endl;
295  new G4PVPlacement(0, G4ThreeVector(0, (row*_tower_dy), 0),
296  DRCalRowLogical, name_row_twr.str().c_str(), drcalo_envelope_log, 0, false, OverlapCheck());
297  }
298  }
299  }
300 
301  return;
302 }
303 
304 //_______________________________________________________________________
305 G4LogicalVolume*
306 PHG4ForwardDualReadoutDetector::ConstructTower(int type)
307 {
308  if (Verbosity() > 0)
309  {
310  cout << "PHG4ForwardDualReadoutDetector: Build logical volume for single tower..." << endl;
311  }
312 
313  //create logical volume for single tower
314  G4Material* material_air = GetDetectorMaterial("G4_AIR");
315  // constructed tower base element
316  G4VSolid* base_tower_solid = new G4Box(G4String("base_tower_solid"),
317  _tower_dx / 2.0,
318  _tower_dy / 2.0,
319  _tower_dz / 2.0);
320 
321  G4LogicalVolume* base_tower_logic = new G4LogicalVolume(base_tower_solid,
322  material_air,
323  "base_tower_logic",
324  0, 0, 0);
325  int maxsubtow = (int) ( (_tower_dx) / (_tower_readout));
326  G4double addtowsize = (_tower_dx - (maxsubtow * _tower_readout))/maxsubtow;
327  // 2x2 fiber tower base element
328  G4VSolid* single_tower_solid = new G4Box(G4String("single_tower_solid"),
329  (_tower_readout + addtowsize) / 2.0,
330  (_tower_readout + addtowsize) / 2.0,
331  _tower_dz / 2.0);
332 
333  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid,
334  material_air,
335  "single_tower_logic",
336  0, 0, 0);
337 
338  m_DisplayAction->AddVolume(single_tower_logic, "FdrcaloEnvelope");
339  //create geometry volumes to place inside single_tower
340 
341  G4double diameter_fiber = _scintFiber_diam;
342  G4double diameter_fiber_cherenkov = _cerenkovFiber_diam;
343  G4double airgap = 0.1 * mm;
344 
345  // fibre cutout
346  G4VSolid* single_cutout_tube = new G4Tubs(G4String("single_cutout_tube"),
347  0,
348  ( diameter_fiber + airgap ) / 2.0,
349  1.03 * _tower_dz / 1.0, //make it 1.03 times longer to ensure full cutout
350  0.,2*M_PI*rad);
351  G4VSolid* single_cutout_tube_cherenkov = new G4Tubs(G4String("single_cutout_tube_cherenkov"),
352  0,
353  ( diameter_fiber_cherenkov + airgap ) / 2.0,
354  1.03 * _tower_dz / 1.0, //make it 1.03 times longer to ensure full cutout
355  0.,2*M_PI*rad);
356  // notch cutout
357  G4VSolid* single_cutout_box = new G4Box(G4String("single_cutout_box"),
358  ( diameter_fiber + airgap ) / 2.0,
359  1.03 * ( diameter_fiber + airgap ) / 4.0, //make it 1.03 times longer to ensure full cutout
360  1.03 * _tower_dz / 1.0);
361  G4VSolid* single_cutout_box_cherenkov = new G4Box(G4String("single_cutout_box_cherenkov"),
362  ( diameter_fiber_cherenkov + airgap ) / 2.0,
363  1.03 * ( diameter_fiber_cherenkov + airgap ) / 4.0, //make it 1.03 times longer to ensure full cutout
364  1.03 * _tower_dz / 1.0);
365  // absorber base object
366  G4VSolid* solid_absorber_cher = new G4Box(G4String("solid_absorber_temp_cher"),
367  (_tower_readout + addtowsize) / 4.0,
368  (_tower_readout + addtowsize) / 4.0,
369  _tower_dz / 2.0);
370  G4VSolid* solid_absorber_scin = new G4Box(G4String("solid_absorber_temp_scin"),
371  (_tower_readout + addtowsize) / 4.0,
372  (_tower_readout + addtowsize) / 4.0,
373  _tower_dz / 2.0);
374 
375  if(_tower_makeNotched){
376  // cut out fiber hole
377  solid_absorber_cher = new G4SubtractionSolid(G4String("solid_absorber_cher_f1"), solid_absorber_cher, single_cutout_tube_cherenkov
378  , 0 ,G4ThreeVector( 0 , ( (_tower_readout + addtowsize) / 4.0 ) - ( ( diameter_fiber_cherenkov + airgap ) / 2.0 ) ,0.)); // top right
379  solid_absorber_scin = new G4SubtractionSolid(G4String("solid_absorber_scin_f1"), solid_absorber_scin, single_cutout_tube
380  , 0 ,G4ThreeVector( 0 , ( (_tower_readout + addtowsize) / 4.0 ) - ( ( diameter_fiber + airgap ) / 2.0 ) ,0.)); // top right
381  // cut out notch
382  solid_absorber_cher = new G4SubtractionSolid(G4String("solid_absorber_cher_box1"), solid_absorber_cher, single_cutout_box_cherenkov
383  , 0 ,G4ThreeVector( 0 , ( (_tower_readout + addtowsize) / 4.0 ) - ( ( diameter_fiber_cherenkov + airgap ) / 4.0 ) ,0.));
384  solid_absorber_scin = new G4SubtractionSolid(G4String("solid_absorber_scin_box1"), solid_absorber_scin, single_cutout_box
385  , 0 ,G4ThreeVector( 0 , ( (_tower_readout + addtowsize) / 4.0 ) - ( ( diameter_fiber + airgap ) / 4.0 ) ,0.));
386  } else {
387  solid_absorber_cher = new G4SubtractionSolid(G4String("solid_absorber_temp_cher_f1"), solid_absorber_cher, single_cutout_tube_cherenkov
388  , 0 ,G4ThreeVector( 0 , 0 ,0.)); // top right
389  solid_absorber_scin = new G4SubtractionSolid(G4String("solid_absorber_temp_scin_f1"), solid_absorber_scin, single_cutout_tube
390  , 0 ,G4ThreeVector( 0 , 0 ,0.)); // top right
391  }
392  G4VSolid* solid_scintillator = new G4Tubs(G4String("single_scintillator_fiber"),
393  0,
394  diameter_fiber / 2.0,
395  _tower_dz / 2.0,
396  0.,2*M_PI*rad);
397 
398 
399  G4VSolid* solid_cherenkov = new G4Tubs(G4String("single_cherenkov_fiber"),
400  0,
401  diameter_fiber_cherenkov / 2.0,
402  _tower_dz / 2.0,
403  0.,2*M_PI*rad);
404 
405 
406  G4Material* material_scintillator = GetScintillatorMaterial();
407 
408 
409  G4NistManager* man = G4NistManager::Instance();
410  G4Material* material_absorber;
411  if(_absorber_Material==0)material_absorber = man->FindOrBuildMaterial(_materialAbsorber.c_str());
412  else if(_absorber_Material==1)material_absorber = man->FindOrBuildMaterial("G4_W");
413  else if(_absorber_Material==2)material_absorber = man->FindOrBuildMaterial("G4_Cu");
414  else if(_absorber_Material==3)material_absorber = man->FindOrBuildMaterial("G4_Pb");
415  else material_absorber = man->FindOrBuildMaterial(_materialAbsorber.c_str());
416 
417 
418  G4LogicalVolume* logic_absorber_cher = new G4LogicalVolume(solid_absorber_cher,
419  material_absorber,
420  "absorber_solid_logic_cher",
421  0, 0, 0);
422 
423  G4LogicalVolume* logic_absorber_scin = new G4LogicalVolume(solid_absorber_scin,
424  material_absorber,
425  "absorber_solid_logic_scin",
426  0, 0, 0);
427 
428  G4LogicalVolume* logic_scint = new G4LogicalVolume(solid_scintillator,
429  material_scintillator,
430  "hdrcalo_single_scintillator_fiber_logic",
431  0, 0, 0);
432 
433  G4Material *material_cherenkov;
434  if(_cerenkovFiber_material==0) material_cherenkov = GetPMMAMaterial();
435  else if(_cerenkovFiber_material==1) material_cherenkov = GetQuartzMaterial();
436  else material_cherenkov = GetPMMAMaterial();
437 
438  G4LogicalVolume* logic_cherenk = new G4LogicalVolume(solid_cherenkov,
439  material_cherenkov,
440  "hdrcalo_single_cherenkov_fiber_logic",
441  0, 0, 0);
442 
443  m_DisplayAction->AddVolume(logic_absorber_cher, "Absorber");
444  m_DisplayAction->AddVolume(logic_absorber_scin, "Absorber");
445 
446  m_DisplayAction->AddVolume(logic_scint, "Scintillator");
447  m_DisplayAction->AddVolume(logic_cherenk, "Cherenkov");
448 
449  //place physical volumes for absorber and scintillator fiber
450 
451  ostringstream name_absorber;
452  name_absorber.str("");
453  name_absorber << _towerlogicnameprefix << "absorbersolid" << endl;
454 
455  ostringstream name_scintillator;
456  name_scintillator.str("");
457  name_scintillator << _towerlogicnameprefix << "singlescintillatorfiber" << endl;
458 
459  ostringstream name_cherenkov;
460  name_cherenkov.str("");
461  name_cherenkov << _towerlogicnameprefix << "_single_cherenkov_fiber" << endl;
462 
463 
464  new G4PVPlacement(0, G4ThreeVector( (_tower_readout + addtowsize) / 4.0, (_tower_readout + addtowsize) / 4.0 , 0),
465  logic_absorber_cher,
466  name_absorber.str().c_str()+std::to_string(1),
467  single_tower_logic,
468  0, 0, OverlapCheck());
469 
470  new G4PVPlacement(0, G4ThreeVector( (_tower_readout + addtowsize) / 4.0, -(_tower_readout + addtowsize) / 4.0 , 0),
471  logic_absorber_scin,
472  name_absorber.str().c_str()+std::to_string(2),
473  single_tower_logic,
474  0, 0, OverlapCheck());
475 
476  new G4PVPlacement(0, G4ThreeVector( -(_tower_readout + addtowsize) / 4.0, (_tower_readout + addtowsize) / 4.0 , 0),
477  logic_absorber_scin,
478  name_absorber.str().c_str()+std::to_string(3),
479  single_tower_logic,
480  0, 0, OverlapCheck());
481 
482  new G4PVPlacement(0, G4ThreeVector( -(_tower_readout + addtowsize) / 4.0, -(_tower_readout + addtowsize) / 4.0 , 0),
483  logic_absorber_cher,
484  name_absorber.str().c_str()+std::to_string(4),
485  single_tower_logic,
486  0, 0, OverlapCheck());
487 
488  if(_tower_makeNotched){
489  // place scintillator fibers (top left, bottom right)
490  new G4PVPlacement(0, G4ThreeVector( -(_tower_readout + addtowsize) / 4.0, ( (_tower_readout + addtowsize) / 2.0 ) - ( ( diameter_fiber + airgap ) / 2.0 ) , 0),
491  logic_scint,
492  name_scintillator.str().c_str()+std::to_string(1),
493  single_tower_logic,
494  0, 0, OverlapCheck());
495  new G4PVPlacement(0, G4ThreeVector( (_tower_readout + addtowsize) / 4.0, - ( ( diameter_fiber + airgap ) / 2.0 ) , 0),
496  logic_scint,
497  name_scintillator.str().c_str()+std::to_string(2),
498  single_tower_logic,
499  0, 0, OverlapCheck());
500 
501  // place cherenkov fibers (top right, bottom left)
502  new G4PVPlacement(0, G4ThreeVector( (_tower_readout + addtowsize) / 4.0 , ( (_tower_readout + addtowsize) / 2.0 ) - ( ( diameter_fiber_cherenkov + airgap ) / 2.0 ) , 0),
503  logic_cherenk,
504  name_cherenkov.str().c_str()+std::to_string(1),
505  single_tower_logic,
506  0, 0, OverlapCheck());
507  new G4PVPlacement(0, G4ThreeVector( -(_tower_readout + addtowsize) / 4.0 , - ( ( diameter_fiber_cherenkov + airgap ) / 2.0 ) , 0),
508  logic_cherenk,
509  name_cherenkov.str().c_str()+std::to_string(2),
510  single_tower_logic,
511  0, 0, OverlapCheck());
512  } else {
513  // place scintillator fibers (top left, bottom right)
514  new G4PVPlacement(0, G4ThreeVector( -(_tower_readout + addtowsize) / 4.0, (_tower_readout + addtowsize) / 4.0 , 0),
515  logic_scint,
516  name_scintillator.str().c_str()+std::to_string(1),
517  single_tower_logic,
518  0, 0, OverlapCheck());
519  new G4PVPlacement(0, G4ThreeVector( (_tower_readout + addtowsize) / 4.0, - (_tower_readout + addtowsize) / 4.0 , 0),
520  logic_scint,
521  name_scintillator.str().c_str()+std::to_string(2),
522  single_tower_logic,
523  0, 0, OverlapCheck());
524 
525  // place cherenkov fibers (top right, bottom left)
526  new G4PVPlacement(0, G4ThreeVector( (_tower_readout + addtowsize) / 4.0 , (_tower_readout + addtowsize) / 4.0 , 0),
527  logic_cherenk,
528  name_cherenkov.str().c_str()+std::to_string(3),
529  single_tower_logic,
530  0, 0, OverlapCheck());
531  new G4PVPlacement(0, G4ThreeVector( -(_tower_readout + addtowsize) / 4.0 , - (_tower_readout + addtowsize) / 4.0 , 0),
532  logic_cherenk,
533  name_cherenkov.str().c_str()+std::to_string(4),
534  single_tower_logic,
535  0, 0, OverlapCheck());
536  }
537 
538 
539  int rowNtow = (int) ( (_tower_dx) / (_tower_readout + addtowsize));
540  for(int row=(rowNtow / 2);row>=-rowNtow / 2;row--){
541  // create mother volume with space for currRowNtow towers along x-axis
542  auto DRCalRowSolid = new G4Box("DRCalRowBoxBase" + std::to_string(row), _tower_dx / 2.0, (_tower_readout + addtowsize) / 2.0, _tower_dz / 2.0);
543  auto DRCalRowLogical = new G4LogicalVolume(DRCalRowSolid,material_air,"DRCalRowLogicalBase" + std::to_string(row));
544  // replicate singletower tower design currRowNtow times along x-axis
545  new G4PVReplica("DRCalRowPhysicalBase" + std::to_string(row),single_tower_logic,DRCalRowLogical,
546  kXAxis, rowNtow, _tower_readout + addtowsize);
547 
548  m_DisplayAction->AddVolume(DRCalRowLogical, "FdrcaloEnvelope");
549  ostringstream name_row_twr;
550  name_row_twr.str("");
551  name_row_twr << _towerlogicnameprefix << "_row_" << row << endl;
552  new G4PVPlacement(0, G4ThreeVector(0, (row * (_tower_readout + addtowsize)), 0),
553  DRCalRowLogical, name_row_twr.str().c_str(), base_tower_logic, 0, false, OverlapCheck());
554  }
555 
556  if (Verbosity() > 0)
557  {
558  cout << "PHG4ForwardDualReadoutDetector: Building logical volume for single tower done." << endl;
559  }
560 
561  return base_tower_logic;
562 }
563 
564 //_______________________________________________________________________
565 G4LogicalVolume*
566 PHG4ForwardDualReadoutDetector::ConstructTowerFCStyle(int type)
567 {
568  if (Verbosity() > 0)
569  {
570  cout << "PHG4ForwardDualReadoutDetector: Build logical volume for single tower..." << endl;
571  }
572 
573  //create logical volume for single tower
574  G4Material* material_air = GetDetectorMaterial("G4_AIR");
575  // 2x2 tower base element
576  G4VSolid* base_tower_solid = new G4Box(G4String("base_tower_solid"),
577  _tower_dx / 2.0,
578  _tower_dy / 2.0,
579  _tower_dz / 2.0);
580 
581  G4LogicalVolume* base_tower_logic = new G4LogicalVolume(base_tower_solid,
582  material_air,
583  "base_tower_logic",
584  0, 0, 0);
585  G4double copperTubeDiam = _tower_readout / 2;
586  int maxsubtow = (int) ( (_tower_dx) / (2 * copperTubeDiam));
587  G4double addtowsize = (_tower_dx - (maxsubtow * 2 * copperTubeDiam))/maxsubtow;
588  // 2x2 tower base element
589  G4VSolid* single_tower_solid = new G4Box(G4String("single_tower_solid"),
590  (2 * copperTubeDiam + addtowsize) / 2.0,
591  (2 * copperTubeDiam + addtowsize) / 2.0,
592  _tower_dz / 2.0);
593 
594  G4LogicalVolume* single_tower_logic = new G4LogicalVolume(single_tower_solid,
595  material_air,
596  "single_tower_logic",
597  0, 0, 0);
598 
599  m_DisplayAction->AddVolume(single_tower_logic, "FdrcaloEnvelope");
600  //create geometry volumes to place inside single_tower
601 
602  G4double diameter_fiber = _scintFiber_diam;
603  G4double diameter_fiber_cherenkov = _cerenkovFiber_diam;
604 
605  // fibre cutout
606  G4VSolid* solid_absorber = new G4Tubs(G4String("ttl_copper_tube_solid"),
607  0,
608  copperTubeDiam / 2.0,
609  _tower_dz / 2.0,
610  0.,2*M_PI*rad);
611 
612  G4VSolid* solid_scintillator = new G4Tubs(G4String("single_scintillator_fiber"),
613  0,
614  diameter_fiber / 2.0,
615  _tower_dz / 2.0,
616  0.,2*M_PI*rad);
617 
618 
619  G4VSolid* solid_cherenkov = new G4Tubs(G4String("single_cherenkov_fiber"),
620  0,
621  diameter_fiber_cherenkov / 2.0,
622  _tower_dz / 2.0,
623  0.,2*M_PI*rad);
624 
625 
626  G4Material* material_scintillator = GetScintillatorMaterial();
627 
628 
629  G4NistManager* man = G4NistManager::Instance();
630  G4Material* material_absorber;
631  if(_absorber_Material==0)material_absorber = man->FindOrBuildMaterial(_materialAbsorber.c_str());
632  else if(_absorber_Material==1)material_absorber = man->FindOrBuildMaterial("G4_W");
633  else if(_absorber_Material==2)material_absorber = man->FindOrBuildMaterial("G4_Cu");
634  else if(_absorber_Material==3)material_absorber = man->FindOrBuildMaterial("G4_Pb");
635  else material_absorber = man->FindOrBuildMaterial(_materialAbsorber.c_str());
636 
637 
638  G4LogicalVolume* logic_absorber = new G4LogicalVolume(solid_absorber,
639  material_absorber,
640  "absorber_solid_logic",
641  0, 0, 0);
642 
643  G4LogicalVolume* logic_scint = new G4LogicalVolume(solid_scintillator,
644  material_scintillator,
645  "hdrcalo_single_scintillator_fiber_logic",
646  0, 0, 0);
647 
648  G4Material *material_cherenkov;
649  if(_cerenkovFiber_material==0) material_cherenkov = GetPMMAMaterial();
650  else if(_cerenkovFiber_material==1) material_cherenkov = GetQuartzMaterial();
651  else material_cherenkov = GetPMMAMaterial();
652 
653  G4LogicalVolume* logic_cherenk = new G4LogicalVolume(solid_cherenkov,
654  material_cherenkov,
655  "hdrcalo_single_cherenkov_fiber_logic",
656  0, 0, 0);
657 
658  m_DisplayAction->AddVolume(logic_absorber, "Absorber");
659  m_DisplayAction->AddVolume(logic_scint, "Scintillator");
660  m_DisplayAction->AddVolume(logic_cherenk, "Cherenkov");
661 
662  //place physical volumes for absorber and scintillator fiber
663 
664  ostringstream name_absorber;
665  name_absorber.str("");
666  name_absorber << _towerlogicnameprefix << "absorbersolid" << endl;
667 
668  ostringstream name_scintillator;
669  name_scintillator.str("");
670  name_scintillator << _towerlogicnameprefix << "singlescintillatorfiber" << endl;
671 
672  ostringstream name_cherenkov;
673  name_cherenkov.str("");
674  name_cherenkov << _towerlogicnameprefix << "_single_cherenkov_fiber" << endl;
675 
676  // place copper rods
677  new G4PVPlacement(0, G4ThreeVector( -copperTubeDiam/2, copperTubeDiam/2 , 0),
678  logic_absorber,
679  name_absorber.str().c_str(),
680  single_tower_logic,
681  0, 0, OverlapCheck());
682  new G4PVPlacement(0, G4ThreeVector( -copperTubeDiam/2, -copperTubeDiam/2 , 0),
683  logic_absorber,
684  name_absorber.str().c_str(),
685  single_tower_logic,
686  0, 0, OverlapCheck());
687  new G4PVPlacement(0, G4ThreeVector( copperTubeDiam/2, copperTubeDiam/2 , 0),
688  logic_absorber,
689  name_absorber.str().c_str(),
690  single_tower_logic,
691  0, 0, OverlapCheck());
692  new G4PVPlacement(0, G4ThreeVector( copperTubeDiam/2, -copperTubeDiam/2 , 0),
693  logic_absorber,
694  name_absorber.str().c_str(),
695  single_tower_logic,
696  0, 0, OverlapCheck());
697 
698  // place scintillator fibers (top left, bottom right)
699  new G4PVPlacement(0, G4ThreeVector( 0, 0 , 0),
700  logic_scint,
701  name_scintillator.str().c_str(),
702  single_tower_logic,
703  0, 0, OverlapCheck());
704  new G4PVPlacement(0, G4ThreeVector( copperTubeDiam, copperTubeDiam , 0),
705  logic_scint,
706  name_scintillator.str().c_str(),
707  single_tower_logic,
708  0, 0, OverlapCheck());
709 
710  // place cherenkov fibers (top right, bottom left)
711  new G4PVPlacement(0, G4ThreeVector( copperTubeDiam, 0 , 0),
712  logic_cherenk,
713  name_cherenkov.str().c_str(),
714  single_tower_logic,
715  0, 0, OverlapCheck());
716  new G4PVPlacement(0, G4ThreeVector( 0, copperTubeDiam , 0),
717  logic_cherenk,
718  name_cherenkov.str().c_str(),
719  single_tower_logic,
720  0, 0, OverlapCheck());
721 
722 
723  int rowNtow = (int) ( (_tower_dx) / (2 * copperTubeDiam + addtowsize));
724  for(int row=(rowNtow / 2);row>=-rowNtow / 2;row--){
725  // create mother volume with space for currRowNtow towers along x-axis
726  auto DRCalRowSolid = new G4Box("DRCalRowBox", _tower_dx / 2.0, (2 * copperTubeDiam + addtowsize) / 2.0, _tower_dz / 2.0);
727  auto DRCalRowLogical = new G4LogicalVolume(DRCalRowSolid,material_air,"DRCalRowLogical");
728  // replicate singletower tower design currRowNtow times along x-axis
729  new G4PVReplica("DRCalRowPhysical",single_tower_logic,DRCalRowLogical,
730  kXAxis, rowNtow, 2 * copperTubeDiam + addtowsize);
731 
732  m_DisplayAction->AddVolume(DRCalRowLogical, "FdrcaloEnvelope");
733  ostringstream name_row_twr;
734  name_row_twr.str("");
735  name_row_twr << _towerlogicnameprefix << "_row_" << row << endl;
736  new G4PVPlacement(0, G4ThreeVector(0, (row * (2 * copperTubeDiam + addtowsize)), 0),
737  DRCalRowLogical, name_row_twr.str().c_str(), base_tower_logic, 0, false, OverlapCheck());
738  }
739 
740  // new G4PVPlacement(0, G4ThreeVector( 0, 0 , 0),
741  // single_tower_logic,
742  // name_cherenkov.str().c_str(),
743  // base_tower_logic,
744  // 0, 0, OverlapCheck());
745 
746  if (Verbosity() > 0)
747  {
748  cout << "PHG4ForwardDualReadoutDetector: Building logical volume for single tower done." << endl;
749  }
750 
751  return base_tower_logic;
752 }
753 
754 //_______________________________________________________________________
755 G4Material*
756 PHG4ForwardDualReadoutDetector::GetScintillatorMaterial()
757 {
758  if (Verbosity() > 0)
759  {
760  cout << "PHG4ForwardDualReadoutDetector: Making Scintillator material..." << endl;
761  }
762 
763 
764  G4MaterialPropertiesTable *tab = new G4MaterialPropertiesTable();
765 
766  const G4int ntab = 31;
767  tab->AddConstProperty("FASTTIMECONSTANT", 2.8*ns); // was 6
768  // tab->AddConstProperty("SCINTILLATIONYIELD", 13.9/keV); // was 200/MEV nominal 10
769  tab->AddConstProperty("SCINTILLATIONYIELD", 200/MeV); // was 200/MEV nominal, should maybe be 13.9/keV
770  tab->AddConstProperty("RESOLUTIONSCALE", 1.0);
771 
772  G4double opt_en[] =
773  { 1.37760*eV, 1.45864*eV, 1.54980*eV, 1.65312*eV, 1.71013*eV, 1.77120*eV, 1.83680*eV, 1.90745*eV, 1.98375*eV, 2.06640*eV,
774  2.10143*eV, 2.13766*eV, 2.17516*eV, 2.21400*eV, 2.25426*eV, 2.29600*eV, 2.33932*eV, 2.38431*eV, 2.43106*eV, 2.47968*eV,
775  2.53029*eV, 2.58300*eV, 2.63796*eV, 2.69531*eV, 2.75520*eV, 2.81782*eV, 2.88335*eV, 2.95200*eV, 3.09960*eV, 3.54241*eV,
776  4.13281*eV }; // 350 - 800 nm
777  G4double scin_fast[] =
778  { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0003, 0.0008, 0.0032,
779  0.0057, 0.0084, 0.0153, 0.0234, 0.0343, 0.0604, 0.0927, 0.1398, 0.2105, 0.2903,
780  0.4122, 0.5518, 0.7086, 0.8678, 1.0000, 0.8676, 0.2311, 0.0033, 0.0012, 0.0000,
781  0 };
782  tab->AddProperty("FASTCOMPONENT", opt_en, scin_fast, ntab);
783 
784  G4double opt_r[] =
785  { 1.5749, 1.5764, 1.5782, 1.5803, 1.5815, 1.5829, 1.5845, 1.5862, 1.5882, 1.5904,
786  1.5914, 1.5924, 1.5935, 1.5947, 1.5959, 1.5972, 1.5986, 1.6000, 1.6016, 1.6033,
787  1.6051, 1.6070, 1.6090, 1.6112, 1.6136, 1.6161, 1.6170, 1.6230, 1.62858, 1.65191,
788  1.69165 };
789  tab->AddProperty("RINDEX", opt_en, opt_r, ntab);
790 
791  G4double opt_abs[] =
792  { 2.714*m, 3.619*m, 5.791*m, 4.343*m, 7.896*m, 5.429*m, 36.19*m, 17.37*m, 36.19*m, 5.429*m,
793  13.00*m, 14.50*m, 16.00*m, 18.00*m, 16.50*m, 17.00*m, 14.00*m, 16.00*m, 15.00*m, 14.50*m,
794  13.00*m, 12.00*m, 10.00*m, 8.000*m, 7.238*m, 4.000*m, 1.200*m, 0.500*m, 0.200*m, 0.200*m,
795  0.100*m };
796  tab->AddProperty("ABSLENGTH", opt_en, opt_abs, ntab);
797 
798  G4double density;
799  G4int ncomponents;
800  // G4Material* material_G4_POLYSTYRENE = GetDetectorMaterial(_materialScintillator.c_str());
801  G4Material* material_G4_POLYSTYRENE = new G4Material("G4_POLYSTYRENE", density = 1.05 * g / cm3, ncomponents = 2);
802  material_G4_POLYSTYRENE->AddElement(G4Element::GetElement("C"), 8);
803  material_G4_POLYSTYRENE->AddElement(G4Element::GetElement("H"), 8);
804  material_G4_POLYSTYRENE->GetIonisation()->SetBirksConstant(0.126*mm/MeV);
805  material_G4_POLYSTYRENE->SetMaterialPropertiesTable(tab);
806 
807  if (Verbosity() > 0)
808  {
809  cout << "PHG4ForwardDualReadoutDetector: Making Scintillator material done." << endl;
810  }
811 
812  return material_G4_POLYSTYRENE;
813 }
814 
815 //_______________________________________________________________________
816 G4Material*
817 PHG4ForwardDualReadoutDetector::GetPMMAMaterial()
818 {
819  if (Verbosity() > 0)
820  {
821  cout << "PHG4ForwardDualReadoutDetector: Making PMMA material..." << endl;
822  }
823 
824  G4double density;
825  G4int ncomponents;
826 
827  G4Material* material_PMMA = new G4Material("PMMA", density = 1.18 * g / cm3, ncomponents = 3);
828  material_PMMA->AddElement(G4Element::GetElement("C"), 5);
829  material_PMMA->AddElement(G4Element::GetElement("H"), 8);
830  material_PMMA->AddElement(G4Element::GetElement("O"), 2);
831 
832  const G4int nEntries = 31;
833 
834  G4double photonEnergy[nEntries] =
835  { 1.37760*eV, 1.45864*eV, 1.54980*eV, 1.65312*eV, 1.71013*eV, 1.77120*eV, 1.83680*eV, 1.90745*eV, 1.98375*eV, 2.06640*eV,
836  2.10143*eV, 2.13766*eV, 2.17516*eV, 2.21400*eV, 2.25426*eV, 2.29600*eV, 2.33932*eV, 2.38431*eV, 2.43106*eV, 2.47968*eV,
837  2.53029*eV, 2.58300*eV, 2.63796*eV, 2.69531*eV, 2.75520*eV, 2.81782*eV, 2.88335*eV, 2.95200*eV, 3.09960*eV, 3.54241*eV,
838  4.13281*eV };
839  G4double refractiveIndexWLSfiber[nEntries] =
840  { 1.4852, 1.4859, 1.4867, 1.4877, 1.4882, 1.4888, 1.4895, 1.4903, 1.4911, 1.4920,
841  1.4924, 1.4929, 1.4933, 1.4938, 1.4943, 1.4948, 1.4954, 1.4960, 1.4966, 1.4973,
842  1.4981, 1.4989, 1.4997, 1.5006, 1.5016, 1.5026, 1.5038, 1.5050, 1.5052, 1.5152,
843  1.5306 };
844 
845  G4double absWLSfiber[nEntries] =
846  { 0.414*m, 0.965*m, 2.171*m, 4.343*m, 1.448*m, 4.343*m, 14.48*m, 21.71*m, 8.686*m, 39.48*m,
847  48.25*m, 54.29*m, 57.91*m, 54.29*m, 33.40*m, 31.02*m, 43.43*m, 43.43*m, 41.36*m, 39.48*m,
848  37.76*m, 36.19*m, 36.19*m, 33.40*m, 31.02*m, 28.95*m, 25.55*m, 24.13*m, 21.71*m, 2.171*m,
849  0.434*m };
850 
851 
852  // Add entries into properties table
853  G4MaterialPropertiesTable* mptWLSfiber = new G4MaterialPropertiesTable();
854  mptWLSfiber->AddProperty("RINDEX",photonEnergy,refractiveIndexWLSfiber,nEntries);
855  mptWLSfiber->AddProperty("ABSLENGTH",photonEnergy,absWLSfiber,nEntries);
856  material_PMMA->SetMaterialPropertiesTable(mptWLSfiber);
857  if (Verbosity() > 0)
858  {
859  cout << "PHG4ForwardDualReadoutDetector: Making PMMA material done." << endl;
860  }
861 
862  return material_PMMA;
863 }
864 
865 //_______________________________________________________________________
866 G4Material*
867 PHG4ForwardDualReadoutDetector::GetQuartzMaterial()
868 {
869  if (Verbosity() > 0)
870  {
871  cout << "PHG4ForwardDualReadoutDetector: Making Scintillator material..." << endl;
872  }
873 
874  G4MaterialPropertiesTable* mptWLSfiber = new G4MaterialPropertiesTable();
875 
876  G4double density;
877  G4int ncomponents;
878 
879  // G4Material* material_Quartz = GetDetectorMaterial("Quartz");
880  G4Material *material_Quartz = new G4Material("Quartz", density = 2.200 * g / cm3, ncomponents = 2);
881  material_Quartz->AddElement(G4Element::GetElement("Si"), 1);
882  material_Quartz->AddElement(G4Element::GetElement("O"), 2);
883 
884  const G4int nEntriesQuartz = 279;
885 
886  G4double photonEnergyQuartz[nEntriesQuartz] =
887  { 5.9040*eV, 5.7667*eV, 5.6356*eV, 5.5104*eV, 5.3906*eV, 5.2759*eV, 5.1660*eV, 5.0606*eV, 4.9594*eV, 4.8621*eV,
888  4.7686*eV, 4.6786*eV, 4.5920*eV, 4.5085*eV, 4.4280*eV, 4.3503*eV, 4.2753*eV, 4.2029*eV, 4.1328*eV, 4.0651*eV,
889  3.9995*eV, 3.9360*eV, 3.8745*eV, 3.8149*eV, 3.7571*eV, 3.7010*eV, 3.6466*eV, 3.5937*eV, 3.5424*eV, 3.4925*eV,
890  3.4440*eV, 3.3968*eV, 3.3509*eV, 3.3062*eV, 3.2627*eV, 3.2204*eV, 3.1791*eV, 3.1388*eV, 3.0996*eV, 3.0613*eV,
891  3.0240*eV, 2.9876*eV, 2.9520*eV, 2.9173*eV, 2.8834*eV, 2.8502*eV, 2.8178*eV, 2.7862*eV, 2.7552*eV, 2.7249*eV,
892  2.6953*eV, 2.6663*eV, 2.6380*eV, 2.6102*eV, 2.5830*eV, 2.5564*eV, 2.5303*eV, 2.5047*eV, 2.4797*eV, 2.4551*eV,
893  2.4311*eV, 2.4075*eV, 2.3843*eV, 2.3616*eV, 2.3393*eV, 2.3175*eV, 2.2960*eV, 2.2749*eV, 2.2543*eV, 2.2339*eV,
894  2.2140*eV, 2.1944*eV, 2.1752*eV, 2.1562*eV, 2.1377*eV, 2.1194*eV, 2.1014*eV, 2.0838*eV, 2.0664*eV, 2.0493*eV,
895  2.0325*eV, 2.0160*eV, 1.9997*eV, 1.9837*eV, 1.9680*eV, 1.9525*eV, 1.9373*eV, 1.9222*eV, 1.9074*eV, 1.8929*eV,
896  1.8785*eV, 1.8644*eV, 1.8505*eV, 1.8368*eV, 1.8233*eV, 1.8100*eV, 1.7969*eV, 1.7839*eV, 1.7712*eV, 1.7463*eV,
897  1.7220*eV, 1.6984*eV, 1.6755*eV, 1.6531*eV, 1.6314*eV, 1.6102*eV, 1.5895*eV, 1.5694*eV, 1.5498*eV, 1.5307*eV,
898  1.5120*eV, 1.4938*eV, 1.4760*eV, 1.4586*eV, 1.4417*eV, 1.4251*eV, 1.4089*eV, 1.3931*eV, 1.3776*eV, 1.3625*eV,
899  1.3477*eV, 1.3332*eV, 1.3190*eV, 1.3051*eV, 1.2915*eV, 1.2782*eV, 1.2651*eV, 1.2524*eV, 1.2398*eV, 1.2276*eV,
900  1.2155*eV, 1.2037*eV, 1.1922*eV, 1.1808*eV, 1.1697*eV, 1.1587*eV, 1.1480*eV, 1.1375*eV, 1.1271*eV, 1.1170*eV,
901  1.1070*eV, 1.0972*eV, 1.0876*eV, 1.0781*eV, 1.0688*eV, 1.0597*eV, 1.0507*eV, 1.0419*eV, 1.0332*eV, 1.0247*eV,
902  1.0163*eV, 1.0080*eV, 0.9999*eV, 0.9919*eV, 0.9840*eV, 0.9763*eV, 0.9686*eV, 0.9611*eV, 0.9537*eV, 0.9464*eV,
903  0.9393*eV, 0.9322*eV, 0.9253*eV, 0.9184*eV, 0.9116*eV, 0.9050*eV, 0.8984*eV, 0.8920*eV, 0.8856*eV, 0.8793*eV,
904  0.8731*eV, 0.8670*eV, 0.8610*eV, 0.8551*eV, 0.8492*eV, 0.8434*eV, 0.8377*eV, 0.8321*eV, 0.8266*eV, 0.8211*eV,
905  0.8157*eV, 0.8104*eV, 0.8051*eV, 0.7999*eV, 0.7948*eV, 0.7897*eV, 0.7847*eV, 0.7798*eV, 0.7749*eV, 0.7701*eV,
906  0.7653*eV, 0.7606*eV, 0.7560*eV, 0.7514*eV, 0.7469*eV, 0.7424*eV, 0.7380*eV, 0.7336*eV, 0.7293*eV, 0.7251*eV,
907  0.7208*eV, 0.7167*eV, 0.7126*eV, 0.7085*eV, 0.7045*eV, 0.7005*eV, 0.6965*eV, 0.6926*eV, 0.6888*eV, 0.6850*eV,
908  0.6812*eV, 0.6775*eV, 0.6738*eV, 0.6702*eV, 0.6666*eV, 0.6630*eV, 0.6595*eV, 0.6560*eV, 0.6525*eV, 0.6491*eV,
909  0.6458*eV, 0.6424*eV, 0.6391*eV, 0.6358*eV, 0.6326*eV, 0.6294*eV, 0.6262*eV, 0.6230*eV, 0.6199*eV, 0.6168*eV,
910  0.6138*eV, 0.6108*eV, 0.6078*eV, 0.6048*eV, 0.6019*eV, 0.5990*eV, 0.5961*eV, 0.5932*eV, 0.5904*eV, 0.5876*eV,
911  0.5848*eV, 0.5821*eV, 0.5794*eV, 0.5767*eV, 0.5740*eV, 0.5714*eV, 0.5687*eV, 0.5661*eV, 0.5636*eV, 0.5610*eV,
912  0.5585*eV, 0.5560*eV, 0.5535*eV, 0.5510*eV, 0.5486*eV, 0.5462*eV, 0.5438*eV, 0.5414*eV, 0.5391*eV, 0.5367*eV,
913  0.5344*eV, 0.5321*eV, 0.5298*eV, 0.5276*eV, 0.5254*eV, 0.5231*eV, 0.5209*eV, 0.5188*eV, 0.5166*eV, 0.5145*eV,
914  0.5123*eV, 0.5102*eV, 0.5081*eV, 0.5061*eV, 0.5040*eV, 0.5020*eV, 0.4999*eV, 0.4979*eV, 0.4959*eV};
915  G4double refractiveIndexQuartz[nEntriesQuartz] =
916  { 1.5384, 1.5332, 1.5285, 1.5242, 1.5202, 1.5166, 1.5133, 1.5103, 1.5074, 1.5048,
917  1.5024, 1.5001, 1.4980, 1.4960, 1.4942, 1.4924, 1.4908, 1.4892, 1.4878, 1.4864,
918  1.4851, 1.4839, 1.4827, 1.4816, 1.4806, 1.4796, 1.4787, 1.4778, 1.4769, 1.4761,
919  1.4753, 1.4745, 1.4738, 1.4731, 1.4725, 1.4719, 1.4713, 1.4707, 1.4701, 1.4696,
920  1.4691, 1.4686, 1.4681, 1.4676, 1.4672, 1.4668, 1.4663, 1.4660, 1.4656, 1.4652,
921  1.4648, 1.4645, 1.4641, 1.4638, 1.4635, 1.4632, 1.4629, 1.4626, 1.4623, 1.4621,
922  1.4618, 1.4615, 1.4613, 1.4610, 1.4608, 1.4606, 1.4603, 1.4601, 1.4599, 1.4597,
923  1.4595, 1.4593, 1.4591, 1.4589, 1.4587, 1.4586, 1.4584, 1.4582, 1.4580, 1.4579,
924  1.4577, 1.4576, 1.4574, 1.4572, 1.4571, 1.4570, 1.4568, 1.4567, 1.4565, 1.4564,
925  1.4563, 1.4561, 1.4560, 1.4559, 1.4558, 1.4556, 1.4555, 1.4554, 1.4553, 1.4551,
926  1.4549, 1.4546, 1.4544, 1.4542, 1.4540, 1.4539, 1.4537, 1.4535, 1.4533, 1.4531,
927  1.4530, 1.4528, 1.4527, 1.4525, 1.4523, 1.4522, 1.4520, 1.4519, 1.4518, 1.4516,
928  1.4515, 1.4513, 1.4512, 1.4511, 1.4509, 1.4508, 1.4507, 1.4505, 1.4504, 1.4503,
929  1.4502, 1.4500, 1.4499, 1.4498, 1.4497, 1.4496, 1.4494, 1.4493, 1.4492, 1.4491,
930  1.4490, 1.4489, 1.4487, 1.4486, 1.4485, 1.4484, 1.4483, 1.4482, 1.4481, 1.4479,
931  1.4478, 1.4477, 1.4476, 1.4475, 1.4474, 1.4473, 1.4471, 1.4470, 1.4469, 1.4468,
932  1.4467, 1.4466, 1.4465, 1.4464, 1.4462, 1.4461, 1.4460, 1.4459, 1.4458, 1.4457,
933  1.4455, 1.4454, 1.4453, 1.4452, 1.4451, 1.4450, 1.4449, 1.4447, 1.4446, 1.4445,
934  1.4444, 1.4443, 1.4441, 1.4440, 1.4439, 1.4438, 1.4437, 1.4435, 1.4434, 1.4433,
935  1.4432, 1.4431, 1.4429, 1.4428, 1.4427, 1.4426, 1.4424, 1.4423, 1.4422, 1.4420,
936  1.4419, 1.4418, 1.4417, 1.4415, 1.4414, 1.4413, 1.4411, 1.4410, 1.4409, 1.4407,
937  1.4406, 1.4405, 1.4403, 1.4402, 1.4401, 1.4399, 1.4398, 1.4397, 1.4395, 1.4394,
938  1.4392, 1.4391, 1.4389, 1.4388, 1.4387, 1.4385, 1.4384, 1.4382, 1.4381, 1.4379,
939  1.4378, 1.4376, 1.4375, 1.4373, 1.4372, 1.4370, 1.4369, 1.4367, 1.4366, 1.4364,
940  1.4363, 1.4361, 1.4360, 1.4358, 1.4357, 1.4355, 1.4353, 1.4352, 1.4350, 1.4349,
941  1.4347, 1.4345, 1.4344, 1.4342, 1.4340, 1.4339, 1.4337, 1.4335, 1.4334, 1.4332,
942  1.4330, 1.4328, 1.4327, 1.4325, 1.4323, 1.4322, 1.4320, 1.4318, 1.4316, 1.4314,
943  1.4313, 1.4311, 1.4309, 1.4307, 1.4305, 1.4304, 1.4302, 1.4300, 1.4298 };
944  mptWLSfiber->AddProperty("RINDEX",photonEnergyQuartz,refractiveIndexQuartz,nEntriesQuartz);
945 
946 
947  const G4int nEntries_Quartz = 14;
948  G4double PhotonEnergy_Quartz[nEntries_Quartz] =
949  { 2.21*eV, 2.30*eV, 2.38*eV, 2.48*eV, 2.58*eV, 2.70*eV, 2.82*eV, 2.95*eV, 3.10*eV, 3.26*eV,
950  3.44*eV, 3.65*eV, 3.88*eV, 4.13*eV };
951  G4double Quartz_Abs[nEntries_Quartz] =
952  { 550.7*mm, 530.7*mm, 590.1*mm, 490.7*mm, 470.7*mm, 520.3*mm, 500.0*mm, 470.7*mm, 450.5*mm, 270.5*mm,
953  190.1*mm, 60.9*mm, 10.6*mm, 4.0*mm};
954  mptWLSfiber->AddProperty("ABSLENGTH", PhotonEnergy_Quartz, Quartz_Abs, nEntries_Quartz);
955 
956  material_Quartz->SetMaterialPropertiesTable(mptWLSfiber);
957 
958  if (Verbosity() > 0)
959  {
960  cout << "PHG4ForwardDualReadoutDetector: Making Scintillator material done." << endl;
961  }
962 
963  return material_Quartz;
964 }
965 
966 
967 int PHG4ForwardDualReadoutDetector::ParseParametersFromTable()
968 {
969  //Open the datafile, if it won't open return an error
970  ifstream istream_mapping;
971  istream_mapping.open(_mapping_tower_file);
972  if (!istream_mapping.is_open())
973  {
974  std::cout << "ERROR in PHG4ForwardHcalDetector: Failed to open mapping file " << _mapping_tower_file << std::endl;
975  gSystem->Exit(1);
976  }
977 
978  //loop over lines in file
979  string line_mapping;
980  while (getline(istream_mapping, line_mapping))
981  {
982  //Skip lines starting with / including a '#'
983  if (line_mapping.find("#") != string::npos)
984  {
985  if (Verbosity() > 0)
986  {
987  std::cout << "PHG4ForwardHcalDetector: SKIPPING line in mapping file: " << line_mapping << std::endl;
988  }
989  continue;
990  }
991 
992  istringstream iss(line_mapping);
993  //If this line is not a comment and not a tower, save parameter as string / value.
994  string parname;
995  G4double parval;
996 
997  //read string- break if error
998  if (!(iss >> parname >> parval))
999  {
1000  cout << "ERROR in PHG4ForwardHcalDetector: Failed to read line in mapping file " << _mapping_tower_file << std::endl;
1001  gSystem->Exit(1);
1002  }
1003 
1004  m_GlobalParameterMap.insert(make_pair(parname, parval));
1005  }
1006 
1007  //Update member variables for global parameters based on parsed parameter file
1008  std::map<string, G4double>::iterator parit;
1009 
1010 
1011  parit = m_GlobalParameterMap.find("Gtype");
1012  if (parit != m_GlobalParameterMap.end())
1013  {
1014  _tower_type = parit->second;
1015  }
1016 
1017  parit = m_GlobalParameterMap.find("Gtower_readout");
1018  if (parit != m_GlobalParameterMap.end())
1019  {
1020  _tower_readout = parit->second * cm;
1021  m_SteppingAction->SetTowerReadout(_tower_readout);
1022  }
1023 
1024  parit = m_GlobalParameterMap.find("Gtower_dx");
1025  if (parit != m_GlobalParameterMap.end())
1026  {
1027  _tower_dx = parit->second * cm;
1028  m_SteppingAction->SetTowerSize(_tower_dx);
1029  }
1030 
1031  parit = m_GlobalParameterMap.find("Gtower_dy");
1032  if (parit != m_GlobalParameterMap.end())
1033  {
1034  _tower_dy = parit->second * cm;
1035  }
1036 
1037  parit = m_GlobalParameterMap.find("Gtower_dz");
1038  if (parit != m_GlobalParameterMap.end())
1039  {
1040  _tower_dz = parit->second * cm;
1041  }
1042 
1043  // new start
1044  parit = m_GlobalParameterMap.find("Scint_Diam");
1045  if (parit != m_GlobalParameterMap.end())
1046  {
1047  _scintFiber_diam = parit->second * cm;
1048  }
1049 
1050  parit = m_GlobalParameterMap.find("Cerenkov_Diam");
1051  if (parit != m_GlobalParameterMap.end())
1052  {
1053  _cerenkovFiber_diam = parit->second * cm;
1054  }
1055 
1056  parit = m_GlobalParameterMap.find("Cerenkov_Material");
1057  if (parit != m_GlobalParameterMap.end())
1058  {
1059  _cerenkovFiber_material = parit->second;
1060  }
1061 
1062  parit = m_GlobalParameterMap.find("NotchCutout");
1063  if (parit != m_GlobalParameterMap.end())
1064  {
1065  _tower_makeNotched = parit->second;
1066  }
1067 
1068  parit = m_GlobalParameterMap.find("Absorber_Material");
1069  if (parit != m_GlobalParameterMap.end())
1070  {
1071  _absorber_Material = parit->second;
1072  }
1073  // new end
1074 
1075  parit = m_GlobalParameterMap.find("Gr1_inner");
1076  if (parit != m_GlobalParameterMap.end())
1077  {
1078  _rMin1 = parit->second * cm;
1079  }
1080 
1081  parit = m_GlobalParameterMap.find("Gr1_outer");
1082  if (parit != m_GlobalParameterMap.end())
1083  {
1084  _rMax1 = parit->second * cm;
1085  m_SteppingAction->SetDetectorSize(_rMax1);
1086  }
1087 
1088  parit = m_GlobalParameterMap.find("Gr2_inner");
1089  if (parit != m_GlobalParameterMap.end())
1090  {
1091  _rMin2 = parit->second * cm;
1092  }
1093 
1094  parit = m_GlobalParameterMap.find("Gr2_outer");
1095  if (parit != m_GlobalParameterMap.end())
1096  {
1097  _rMax2 = parit->second * cm;
1098  }
1099 
1100  parit = m_GlobalParameterMap.find("Gdz");
1101  if (parit != m_GlobalParameterMap.end())
1102  {
1103  _dZ = parit->second * cm;
1104  }
1105 
1106  parit = m_GlobalParameterMap.find("Gx0");
1107  if (parit != m_GlobalParameterMap.end())
1108  {
1109  _place_in_x = parit->second * cm;
1110  }
1111 
1112  parit = m_GlobalParameterMap.find("Gy0");
1113  if (parit != m_GlobalParameterMap.end())
1114  {
1115  _place_in_y = parit->second * cm;
1116  }
1117 
1118  parit = m_GlobalParameterMap.find("Gz0");
1119  if (parit != m_GlobalParameterMap.end())
1120  {
1121  _place_in_z = parit->second * cm;
1122  }
1123 
1124  parit = m_GlobalParameterMap.find("Center_Offset_x");
1125  if (parit != m_GlobalParameterMap.end())
1126  {
1127  _center_offset_x = parit->second * cm;
1128  }
1129 
1130  parit = m_GlobalParameterMap.find("Center_Offset_y");
1131  if (parit != m_GlobalParameterMap.end())
1132  {
1133  _center_offset_y = parit->second * cm;
1134  }
1135  parit = m_GlobalParameterMap.find("Quadratic_Detector");
1136  if (parit != m_GlobalParameterMap.end())
1137  {
1138  _quadratic_detector = parit->second;
1139  }
1140 
1141  parit = m_GlobalParameterMap.find("Grot_x");
1142  if (parit != m_GlobalParameterMap.end())
1143  {
1144  _rot_in_x = parit->second * cm;
1145  }
1146 
1147  parit = m_GlobalParameterMap.find("Grot_y");
1148  if (parit != m_GlobalParameterMap.end())
1149  {
1150  _rot_in_y = parit->second * cm;
1151  }
1152 
1153  parit = m_GlobalParameterMap.find("Grot_z");
1154  if (parit != m_GlobalParameterMap.end())
1155  {
1156  _rot_in_z = parit->second * cm;
1157  }
1158 
1159  return 0;
1160 }