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