ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4B0ECALDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4B0ECALDetector.cc
1 //____________________________________________________________________________..
2 //
3 // This is a working template for the G4 Construct() method which needs to be implemented
4 // We wedge a method between the G4 Construct() to enable volume hierarchies on the macro
5 // so here it is called ConstructMe() but there is no functional difference
6 // Currently this installs a simple G4Box solid, creates a logical volume from it
7 // and places it. Put your own detector in place (just make sure all active volumes
8 // get inserted into the m_PhysicalVolumesSet)
9 //
10 // Rather than using hardcoded values you should consider using the parameter class
11 // Parameter names and defaults are set in EICG4B0Subsystem::SetDefaultParameters()
12 // Only parameters defined there can be used (also to override in the macro)
13 // to avoids typos.
14 // IMPORTANT: parameters have no inherent units, there is a convention (cm/deg)
15 // but in any case you need to multiply them here with the correct CLHEP/G4 unit
16 //
17 // The place where you put your own detector is marked with
18 // //begin implement your own here://
19 // //end implement your own here://
20 // Do not forget to include the G4 includes for your volumes
21 //____________________________________________________________________________..
22 
23 #include "EICG4B0ECALDetector.h"
24 //#include "EICG4B0SteppingAction.h"
25 //#include "EICG4B0Subsystem.h"
26 
27 #include <phparameter/PHParameters.h>
28 
29 #include <g4main/PHG4Detector.h>
30 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
31 //#include <g4main/PHG4Subsystem.h>
32 
33 #include <TSystem.h>
34 #include <Geant4/G4Box.hh>
35 #include <Geant4/G4Color.hh>
36 #include <Geant4/G4LogicalVolume.hh>
37 #include <Geant4/G4Material.hh>
38 #include <Geant4/G4PVPlacement.hh>
39 #include <Geant4/G4RotationMatrix.hh>
40 #include <Geant4/G4SubtractionSolid.hh>
41 #include <Geant4/G4SystemOfUnits.hh>
42 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
43 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
44 #include <Geant4/G4Tubs.hh>
45 #include <Geant4/G4Types.hh> // for G4double, G4int
46 #include <Geant4/G4UnionSolid.hh>
47 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
48 #include <Geant4/G4VisAttributes.hh>
49 
50 #include <phool/recoConsts.h> //For rc WorldMaterial
51 
52 #include <cmath>
53 #include <cstdlib>
54 #include <fstream>
55 #include <iostream>
56 #include <sstream>
57 #include <utility>
58 
59 class G4VSolid;
60 class PHCompositeNode;
61 
62 //____________________________________________________________________________..
64  PHCompositeNode *Node,
66  const std::string &dnam, const int lyr)
67  : PHG4Detector(subsys, Node, dnam)
68  , m_Params(parameters)
69  , m_Layer(lyr)
70  , _mapping_tower_file("")
71  , m_TowerLogicNamePrefix("B0ECALTower")
72 {
73 }
74 
75 //_______________________________________________________________
77 {
78  G4LogicalVolume *mylogvol = volume->GetLogicalVolume();
79  if (m_LogicalVolSet.find(mylogvol) != m_LogicalVolSet.end())
80  {
81  return 1;
82  }
83  return 0;
84 }
85 
87 {
88  if (IsInDetector(volume))
89  {
90  return 1;
91  }
92  return -1;
93 }
94 
95 //_______________________________________________________________
97 {
98  //begin implement your own here://
99  // Do not forget to multiply the parameters with their respective CLHEP/G4 unit !
100  if (Verbosity() > 0)
101  {
102  std::cout << "EICG4B0ECALDetector: Begin Construction" << std::endl;
103  }
104 
105  std::cout << " !!! length = " << m_Params->get_double_param("length");
106  if (m_Params->get_double_param("spanningAngle") >= 360)
107  {
108  std::cout << " !!! No PACKMAN" << std::endl;
109  return;
110  }
111  //Print("ALL");
112 
113  G4VSolid *solid0 = new G4Tubs("EICG4B0ECALSolid0",
114  0.,
115  m_Params->get_double_param("outer_radius") * cm,
116  m_Params->get_double_param("length") / 2. * cm,
117  m_Params->get_double_param("startAngle") * degree,
118  m_Params->get_double_param("spanningAngle") * degree);
119  G4VSolid *solidPipeHole = new G4Tubs("EICG4B0IonPipeSolid",
120  0.,
121  m_Params->get_double_param("pipe_hole_r") * cm,
122  m_Params->get_double_param("length") * cm,
123  0., 360. * degree);
124  G4VSolid *solidCableHole = new G4Tubs("EICG4B0CableSolid",
125  0.,
126  m_Params->get_double_param("cable_hole") * cm,
127  m_Params->get_double_param("length") * cm,
128  0., 360. * degree);
129  G4VSolid *solidPipeHole1 = new G4Box("EICG4B0PipeSolid1",
130  m_Params->get_double_param("pipe_hole") / 2. * cm,
131  m_Params->get_double_param("pipe_hole_r") * cm,
132  m_Params->get_double_param("length") * cm);
133  G4VSolid *solid1 = new G4Tubs("EICG4B0ECALSolid1",
134  0.,
135  (m_Params->get_double_param("outer_radius") - m_Params->get_double_param("d_radius")) * cm,
136  m_Params->get_double_param("length") / 2. * cm,
137  (m_Params->get_double_param("startAngle") + m_Params->get_double_param("spanningAngle")) * degree,
138  (360 - m_Params->get_double_param("spanningAngle")) * degree);
139  G4UnionSolid *solid10 = new G4UnionSolid("EICG4B0ECALSolid10", solid0, solid1);
140  G4SubtractionSolid *solids = new G4SubtractionSolid("EICG4B0Solid", solid10, solidPipeHole, 0, G4ThreeVector((m_Params->get_double_param("pipe_x") + m_Params->get_double_param("pipe_hole") / 2) * cm, m_Params->get_double_param("pipe_y") * cm, m_Params->get_double_param("pipe_z") * cm));
141  G4SubtractionSolid *solids1 = new G4SubtractionSolid("EICG4B0Solid", solids, solidPipeHole, 0, G4ThreeVector((m_Params->get_double_param("pipe_x") - m_Params->get_double_param("pipe_hole") / 2) * cm, m_Params->get_double_param("pipe_y") * cm, m_Params->get_double_param("pipe_z") * cm));
142  G4SubtractionSolid *solids2 = new G4SubtractionSolid("EICG4B0Solid", solids1, solidCableHole, 0, G4ThreeVector(m_Params->get_double_param("cable_x") * cm, m_Params->get_double_param("cable_y") * cm, m_Params->get_double_param("cable_z") * cm));
143  G4SubtractionSolid *solidB0 = new G4SubtractionSolid("EICG4B0Solid", solids2, solidPipeHole1, 0, G4ThreeVector(m_Params->get_double_param("pipe_x") * cm, m_Params->get_double_param("pipe_y") * cm, m_Params->get_double_param("pipe_z") * cm));
144  G4RotationMatrix *rotm = new G4RotationMatrix();
145  if (_mapping_tower_file.empty())
146  {
147  std::cout << "ERROR in EICG4B0EcalDetector: No mapping file specified. Abort detector construction." << std::endl;
148  // gSystem->Exit(1);
149  }
150  /* Read parameters for detector construction and mappign from file */
154  G4LogicalVolume *b0_ecal_log = new G4LogicalVolume(solidB0,
155  WorldMaterial,
156  "B0ECAL_envelope",
157  0, 0, 0);
158 
159  G4VisAttributes *vis = new G4VisAttributes(G4Color(0.8, 0.4, 0.2, 1.0));
160  vis->SetColor(0.8, 0.4, 0.2, 1.0);
161  b0_ecal_log->SetVisAttributes(vis);
162  /* Place envelope cone in simulation */
163  std::string name_envelope = m_TowerLogicNamePrefix + "_envelope";
164 
165  G4VPhysicalVolume *phy = new G4PVPlacement(rotm, G4ThreeVector(m_Params->get_double_param("place_x") * cm, m_Params->get_double_param("place_y") * cm, m_Params->get_double_param("place_z") * cm),
166  b0_ecal_log, name_envelope, logicWorld, 0, false, OverlapCheck());
167 
168  //Create towers for the B0 Ecal:
169  /* Construct single calorimeter tower */
170  if (Verbosity() > 0)
171  {
172  std::cout << "B0 ECal Envelope Location: x, y, z:" << std::endl;
173  std::cout << "Building Calorimeter from " << m_Params->get_string_param("material") << std::endl;
174  std::cout << m_Params->get_double_param("place_x") << "\t";
175  std::cout << m_Params->get_double_param("place_y") << "\t";
176  std::cout << m_Params->get_double_param("place_z") << std::endl;
177  }
178  G4LogicalVolume *singletower = ConstructTower();
179  /* Place calorimeter tower within envelope */
180  PlaceTower(b0_ecal_log, singletower);
181  m_PhysicalVolumesSet.insert(phy);
182  // hard code detector id to detid
183  m_PhysicalVolumesDet.insert({phy, m_Params->get_double_param("detid") + 1});
184  // m_LogicalVolumesSet.insert(logical);
185  //end implement your own here://
186 
187  return;
188 }
189 //_______________________________________________________________
191 {
192  if (Verbosity() > 0)
193  {
194  std::cout << "EICG4B0ECALDetector: Build logical volume for single tower..." << std::endl;
195  }
196 
197  /* create logical volume for single tower */
198  G4Material *EcalMaterial = GetDetectorMaterial(m_Params->get_string_param("material"));
199  double TowerDx = m_Params->get_double_param("tower_size") * cm;
200  double TowerDy = m_Params->get_double_param("tower_size") * cm;
201  double TowerDz = m_Params->get_double_param("length") * cm;
202  if (Verbosity() > 0)
203  {
204  std::cout << "EICG4B0ECALDetector: Construct tower " << m_Params->get_double_param("tower_size") << "\t"
205  << m_Params->get_double_param("length") << std::endl;
206  }
207  G4VSolid *single_tower_solid = new G4Box("single_tower_solid",
208  TowerDx / 2.0,
209  TowerDy / 2.0,
210  TowerDz / 2.0);
211 
212  G4LogicalVolume *single_tower_logic = new G4LogicalVolume(single_tower_solid,
213  EcalMaterial,
214  "single_tower_logic",
215  0, 0, 0);
216 
217  G4VisAttributes *vis = new G4VisAttributes(G4Color(0.8, 0.4, 0.2, .1));
218  vis->SetForceSolid(false);
219  single_tower_logic->SetVisAttributes(vis);
220  m_LogicalVolSet.insert(single_tower_logic);
221  return single_tower_logic;
222 }
223 
225 {
226  /* Loop over all tower positions in std::vector and place tower */
227  for (std::map<std::string, towerposition>::iterator iterator = m_TowerPositionMap.begin(); iterator != m_TowerPositionMap.end(); ++iterator)
228  {
229  if (Verbosity() > 0)
230  {
231  std::cout << "EICG4B0ECALDetector: Place tower " << iterator->first
232  << " idx_j = " << iterator->second.idx_j << ", idx_k = " << iterator->second.idx_k
233  << " at x = " << iterator->second.x / cm << " , y = " << iterator->second.y / cm << " , z = " << iterator->second.z / cm << std::endl;
234  }
235 
236  int copyno = (iterator->second.idx_j << 16) + iterator->second.idx_k;
237  new G4PVPlacement(0, G4ThreeVector(iterator->second.x, iterator->second.y, iterator->second.z),
238  singletower,
239  iterator->first,
240  b0ecalenvelope,
241  0, copyno, OverlapCheck());
242  // 0, 0, OverlapCheck());
243  }
244 
245  return 0;
246 }
247 
249 {
250  /* Open the datafile, if it won't open return an error */
251  std::ifstream istream_mapping;
252  istream_mapping.open(_mapping_tower_file);
253  //istream_mapping.open("B0ECAL_mapping_v0.txt");
254  if (!istream_mapping.is_open())
255  {
256  std::cout << "ERROR in EICG4B0ECALDetector: Failed to open mapping file " << _mapping_tower_file << std::endl;
257  gSystem->Exit(1);
258  }
259 
260  /* loop over lines in file */
261  std::string line_mapping;
262  while (getline(istream_mapping, line_mapping))
263  {
264  /* Skip lines starting with / including a '#' */
265  if (line_mapping.find("#") != std::string::npos)
266  {
267  if (Verbosity() > 0)
268  {
269  std::cout << "EICG4B0ECALDetector: SKIPPING line in mapping file: " << line_mapping << std::endl;
270  }
271  continue;
272  }
273 
274  std::istringstream iss(line_mapping);
275  unsigned idx_j, idx_k, idx_l;
276  G4double pos_x, pos_y, pos_z;
277  double Gpos_x, Gpos_y, Gpos_z, Gpos_z1;
278  G4double size_x, size_y, size_z;
279  G4double rot_x, rot_y, rot_z;
280  G4double dummy;
281  std::string dummys;
282  //G4double GlobalPlaceInX;
283  //G4double GlobalPlaceInY;
284  //G4double GlobalPlaceInZ;
285 
286  if (line_mapping.find("B0 ") != std::string::npos)
287  {
288  if (!(iss >> dummys >> Gpos_x >> Gpos_y >> Gpos_z >> Gpos_z1))
289  {
290  std::cout << "ERROR in EICG4B0Detector: Failed to read global position from mapping file " << _mapping_tower_file << std::endl;
291  gSystem->Exit(1);
292  }
293  if (m_Params->get_double_param("global_x") != Gpos_x || m_Params->get_double_param("global_y") != Gpos_y || m_Params->get_double_param("global_z") != Gpos_z)
294  {
295  std::cout << std::endl;
296  std::cout << "B0 position changed since mapping file produced or wrong mapping is used " << _mapping_tower_file << std::endl;
297  std::cout << m_Params->get_double_param("global_x") << " " << m_Params->get_double_param("global_y") << " " << m_Params->get_double_param("global_z") << std::endl;
298  std::cout << Gpos_x << " " << Gpos_y << " " << Gpos_z << std::endl;
299  std::cout << m_Params->get_double_param("global_x") - Gpos_x << " " << m_Params->get_double_param("global_y") - Gpos_y << " " << m_Params->get_double_param("global_z") - Gpos_z << std::endl;
300  // gSystem->Exit(1);
301  }
302  // GlobalPlaceInX=pos_x*cm;
303  // GlobalPlaceInY=pos_y*cm;
304  // GlobalPlaceInZ=pos_z*cm;
305  }
306  // /* If line starts with keyword Tower, add to tower positions */
307  else if (line_mapping.find("Tower ") != std::string::npos)
308  {
309  /* read std::string- break if error */
310  if (!(iss >> dummys >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> size_x >> size_y >> size_z >> rot_x >> rot_y >> rot_z >> dummy))
311  //if (!(iss >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> rot_x >> rot_y >> rot_z >> dummy))
312  {
313  std::cout << "ERROR in EICG4B0ECALDetector: Failed to read line in mapping file " << _mapping_tower_file << std::endl;
314  gSystem->Exit(1);
315  }
316 
317  /* Construct unique name for tower */
318  /* Mapping file uses cm, this class uses mm for length */
319  std::ostringstream towername;
320  towername.str("");
321  towername << m_TowerLogicNamePrefix << "_j_" << idx_j << "_k_" << idx_k;
322 
323  /* Add Geant4 units */
324  pos_x = pos_x * cm;
325  pos_y = pos_y * cm;
326  pos_z = pos_z * cm;
327 
328  /* insert tower into tower map */
329  towerposition tower_new;
330  tower_new.x = pos_x;
331  tower_new.y = pos_y;
332  tower_new.z = pos_z;
333  tower_new.idx_j = idx_j;
334  tower_new.idx_k = idx_k;
335  m_TowerPositionMap.insert(make_pair(towername.str(), tower_new));
336  }
337  else
338  std::cout << "ERROR in EICG4B0ECALDetector: Unknown line in Mapping File " << _mapping_tower_file << std::endl;
339  }
340 
341  return 0;
342 }
343 
344 //_______________________________________________________________
345 void EICG4B0ECALDetector::Print(const std::string &what) const
346 {
347  std::cout << "EICG4B0ECAL Detector:" << std::endl;
348  if (what == "ALL" || what == "VOLUME")
349  {
350  std::cout << "Version 0.1" << std::endl;
351  std::cout << "Parameters:" << std::endl;
352  m_Params->Print();
353  }
354  return;
355 }
356 
358 {
359  return m_Params;
360 }