ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4EICMvtxDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4EICMvtxDetector.cc
1 #include "PHG4EICMvtxDetector.h"
2 
3 #include "PHG4MvtxDefs.h"
5 
7 
9 
10 #include <phparameter/PHParameters.h>
11 #include <phparameter/PHParametersContainer.h>
12 
13 #include <g4main/PHG4Detector.h> // for PHG4Detector
14 #include <g4main/PHG4DisplayAction.h> // for PHG4DisplayAction
15 #include <g4main/PHG4Subsystem.h> // for PHG4Subsystem
16 
17 #include <phool/PHCompositeNode.h>
18 #include <phool/PHIODataNode.h>
19 #include <phool/PHNode.h> // for PHNode
20 #include <phool/PHNodeIterator.h> // for PHNodeIterator
21 #include <phool/PHObject.h> // for PHObject
22 #include <phool/getClass.h>
23 
24 #include <Geant4/G4AssemblyVolume.hh>
25 #include <Geant4/G4GDMLParser.hh>
26 #include <Geant4/G4GDMLReadStructure.hh> // for G4GDMLReadStructure
27 #include <Geant4/G4LogicalVolume.hh>
28 #include <Geant4/G4Material.hh>
29 #include <Geant4/G4RotationMatrix.hh> // for G4RotationMatrix
30 #include <Geant4/G4String.hh> // for G4String
31 #include <Geant4/G4SystemOfUnits.hh>
32 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
33 #include <Geant4/G4Transform3D.hh> // for G4Transform3D
34 #include <Geant4/G4Types.hh> // for G4double
35 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
36 
37 #include <cmath>
38 #include <cstdio> // for sprintf
39 #include <iostream> // for operator<<, basic...
40 #include <memory>
41 #include <sstream>
42 #include <utility> // for pair, make_pair
43 #include <vector> // for vector, vector<>:...
44 
45 using namespace std;
46 
47 PHG4EICMvtxDetector::PHG4EICMvtxDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const PHParametersContainer* _paramsContainer, const std::string& dnam)
48  : PHG4Detector(subsys, Node, dnam)
49  , m_DisplayAction(dynamic_cast<PHG4MvtxDisplayAction*>(subsys->GetDisplayAction()))
50  , m_ParamsContainer(_paramsContainer)
51  , m_StaveGeometryFile(_paramsContainer->GetParameters(PHG4MvtxDefs::GLOBAL)->get_string_param("stave_geometry_file"))
52 {
53  // Verbosity(2);
54 
55  if (Verbosity() > 0)
56  cout << "PHG4EICMvtxDetector constructor called" << endl;
57 
58  if (Verbosity() > 10)
59  cout << " cm " << cm << " mm " << mm << endl;
60  for (int ilayer = 0; ilayer < n_Layers; ++ilayer)
61  {
62  const PHParameters* params = m_ParamsContainer->GetParameters(ilayer);
63  m_IsLayerActive[ilayer] = params->get_int_param("active");
64  m_IsLayerAbsorberActive[ilayer] = params->get_int_param("absorberactive");
65  m_IsBlackHole[ilayer] = params->get_int_param("blackhole");
66  m_N_staves[ilayer] = params->get_int_param("N_staves");
67  m_nominal_radius[ilayer] = params->get_double_param("layer_nominal_radius");
68  m_nominal_phitilt[ilayer] = params->get_double_param("phitilt");
69  m_nominal_phi0[ilayer] = params->get_double_param("phi0");
70  }
71  /*
72  const PHParameters* alpide_params = m_ParamsContainer->GetParameters(PHG4MvtxDefs::ALPIDE_SEGMENTATION);
73  m_PixelX = alpide_params->get_double_param("pixel_x");
74  m_PixelZ = alpide_params->get_double_param("pixel_z");
75  m_PixelThickness = alpide_params->get_double_param("pixel_thickness");
76  */
77  if (Verbosity() > 0)
78  {
79  cout << "PHG4EICMvtxDetector constructor: making Mvtx detector. " << endl;
80  }
81 }
82 
83 //_______________________________________________________________
84 //_______________________________________________________________
86 {
87  // Is this volume one of the sensors?
88  // Checks if pointer matches one of our stored sensors for this layer
89  if (m_SensorPV.find(volume) != m_SensorPV.end())
90  {
91  if (Verbosity() > 0)
92  {
93  cout << " -- PHG4MvtxTDetector::IsSensor --" << endl;
94  cout << " volume Name : " << volume->GetName() << endl;
95  cout << " -----------------------------------------" << endl;
96  }
97  return 1;
98  }
99 
100  return 0;
101 }
102 
104 {
105  // Does this stave belong to this layer?
106  // Since the Assembly volume read from GDML does not give unique pointers
107  // to sensors, we need to check the stave, which is unique
108  auto iter = m_StavePV.find(volume);
109  if (iter != m_StavePV.end())
110  {
111  tie(layer, stave) = iter->second;
112  if (Verbosity() > 0)
113  {
114  cout << " -- PHG4EICMvtxDetector::IsInMvtx --" << endl;
115  cout << " layer: " << layer << endl;
116  cout << " stave: " << stave << endl;
117  cout << " volume Name : " << volume->GetName() << endl;
118  cout << " stave Name : " << iter->first->GetName() << endl;
119  cout << " -----------------------------------------" << endl;
120  }
121  return 1;
122  }
123 
124  return 0;
125 }
126 
127 int PHG4EICMvtxDetector::get_layer(int index) const
128 {
129  // Get Mvtx layer from stave index in the Mvtx
130  // Mvtx stave index start from 0, i.e index = 0 for stave 0 in layer 0
131  int lay = 0;
132  while (!(index < m_N_staves[lay]))
133  {
134  index -= m_N_staves[lay];
135  lay++;
136  }
137  return lay;
138 }
139 
140 int PHG4EICMvtxDetector::get_stave(int index) const
141 {
142  // Get stave index in the layer from stave index in the Mvtx
143  int lay = 0;
144  while (!(index < m_N_staves[lay]))
145  {
146  index -= m_N_staves[lay];
147  lay++;
148  }
149  return index;
150 }
151 
153 {
154  // This is called from PHG4PhenixDetector::Construct()
155 
156  if (Verbosity() > 0)
157  {
158  cout << endl
159  << "PHG4EICMvtxDetector::Construct called for Mvtx " << endl;
160  }
161  // the tracking layers are placed directly in the world volume, since some layers are (touching) double layers
162  // this reads in the ITS stave geometry from a file and constructs the layer from it
163  ConstructMvtx(logicWorld);
164 
165  // This object provides the strip center locations when given the ladder segment and strip internal cordinates in the sensor
166  AddGeometryNode();
167  return;
168 }
169 
171 {
172  if (Verbosity() > 0)
173  {
174  cout << " PHG4EICMvtxDetector::ConstructMvtx:" << endl;
175  cout << endl;
176  }
177  //===================================
178  // Import the stave physical volume here
179  //===================================
180 
181  // import the staves from the gemetry file
182  std::unique_ptr<G4GDMLReadStructure> reader(new G4GDMLReadStructure());
183  G4GDMLParser gdmlParser(reader.get());
184  gdmlParser.Read(m_StaveGeometryFile, false);
185 
186  // figure out which assembly we want
187  char assemblyname[500];
188  sprintf(assemblyname, "MVTXStave");
189 
190  if (Verbosity() > 0)
191  {
192  cout << "Geting the stave assembly named " << assemblyname << endl;
193  }
194  G4AssemblyVolume* av_ITSUStave = reader->GetAssembly(assemblyname);
195 
196  for (unsigned short ilayer = 0; ilayer < n_Layers; ++ilayer)
197  {
198  if (m_IsLayerActive[ilayer])
199  {
200  if (Verbosity() > 0)
201  {
202  cout << endl;
203  cout << " Constructing Layer " << ilayer << endl;
204  }
205  ConstructMvtx_Layer(ilayer, av_ITSUStave, trackerenvelope);
206  }
207  }
208  FillPVArray(av_ITSUStave);
209  SetDisplayProperty(av_ITSUStave);
210 
211  return 0;
212 }
213 
215 {
216  //=========================================
217  // Now we populate the whole layer with the staves
218  //==========================================
219 
220  int N_staves = m_N_staves[layer];
221  G4double layer_nominal_radius = m_nominal_radius[layer];
222  G4double phitilt = m_nominal_phitilt[layer];
223  G4double phi0 = m_nominal_phi0[layer]; //YCM: azimuthal offset for the first stave
224 
225  if (N_staves < 0)
226  {
227  // The user did not specify how many staves to use for this layer, so we calculate the best value
228  // We choose a phistep that fits N_staves into this radius, but with an arclength separation of AT LEAST arcstep
229  // ideally, the radius would be chosen so that numstaves = N_staves exactly, to get the closest spacing of staves possible without overlaps
230  double arcstep = 12.25;
231  double numstaves = 2.0 * M_PI * layer_nominal_radius / arcstep; // this is just to print out
232  N_staves = int(2.0 * M_PI * layer_nominal_radius / arcstep); // this is the number of staves used
233 
234  // Also suggest the ideal radius for this layer
235  if (Verbosity() > 0)
236  {
237  cout << " Calculated N_staves for layer " /*<< layer*/
238  << " layer_nominal_radius " << layer_nominal_radius
239  << " ITS arcstep " << arcstep
240  << " circumference divided by arcstep " << numstaves
241  << " N_staves " << N_staves
242  << endl;
243  cout << "A radius for this layer of " << (double) N_staves * arcstep / (2.0 * M_PI) + 0.01 << " or "
244  << (double) (N_staves + 1) * arcstep / (2.0 * M_PI) + 0.01 << " would produce perfect stave spacing" << endl;
245  }
246  }
247 
248  G4double phistep = get_phistep(layer); // this produces even stave spacing
249  double z_location = 0.0;
250 
251  if (Verbosity() > 0)
252  {
253  cout << " layer " /*<< layer*/
254  << " layer_nominal_radius " << layer_nominal_radius
255  << " N_staves " << N_staves
256  << " phistep " << phistep
257  << " phitilt " << phitilt
258  << " phi0 " << phi0
259  << endl;
260  }
261 
262  // The stave starts out at (0,0,0) and oriented so that the sensors face upward in y
263  // So we need to rotate the sensor 90 degrees before placing it using phi_offset
264  // note that the gdml file uses a negative phi_offset - different coord system, apparently - the following works
265  double phi_offset = M_PI / 2.0;
266 
267  for (int iphi = 0; iphi < N_staves; iphi++)
268  {
269  // Place the ladder segment envelopes at the correct z and phi
270  // This is the azimuthal angle at which we place the stave
271  // phi0 is the azimuthal offset for the first stave
272  G4double phi_rotation = phi0 + (double) iphi * phistep;
273 
274  G4RotationMatrix Ra;
275  G4ThreeVector Ta;
276 
277  if (Verbosity() > 0)
278  {
279  cout << "phi_offset = " << phi_offset << " iphi " << iphi << " phi_rotation = " << phi_rotation << " phitilt " << phitilt << endl;
280  }
281  // It is first rotated in phi by the azimuthal angle phi_rotation, plus the 90 degrees needed to point the face of the sensor at the origin, plus the tilt (if a tilt is appropriate)
282 
283  // note - if this is layer 0-2, phitilt is the additional tilt for clearance. Otherwise it is zero
284  Ra.rotateZ(phi_rotation + phi_offset + phitilt);
285  // Then translated as follows
286 
287  Ta.setX(layer_nominal_radius * cos(phi_rotation));
288  Ta.setY(layer_nominal_radius * sin(phi_rotation));
289  Ta.setZ(z_location);
290 
291  if (Verbosity() > 0)
292  {
293  cout << " iphi " << iphi << " phi_rotation " << phi_rotation
294  << " x " << layer_nominal_radius * cos(phi_rotation)
295  << " y " << layer_nominal_radius * sin(phi_rotation)
296  << " z " << z_location
297  << endl;
298  }
299  G4Transform3D Tr(Ra, Ta);
300 
301  av_ITSUStave->MakeImprint(trackerenvelope, Tr, 0, OverlapCheck());
302  }
303 
304  if (Verbosity() > 0)
305  {
306  cout << "This layer has a total of " << N_staves << " staves" << endl;
307  }
308  return 0;
309 }
310 
312 {
313  // cout <<"SetDisplayProperty - G4AssemblyVolume w/ TotalImprintedVolumes "<<av->TotalImprintedVolumes()
314  // <<"/"<<av->GetImprintsCount()<<endl;
315 
316  std::vector<G4VPhysicalVolume*>::iterator it = av->GetVolumesIterator();
317 
318  int nDaughters = av->TotalImprintedVolumes();
319  for (int i = 0; i < nDaughters; ++i, ++it)
320  {
321  // cout <<"SetDisplayProperty - AV["<<i<<"] = "<<(*it)->GetName()<<endl;
322  G4VPhysicalVolume* pv = (*it);
323 
324  G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
325  SetDisplayProperty(worldLogical);
326  }
327 }
328 
330 {
331  string material_name(lv->GetMaterial()->GetName());
332 
333  if (Verbosity() >= 50)
334  {
335  cout << "SetDisplayProperty - LV " << lv->GetName() << " built with "
336  << material_name << endl;
337  }
338  vector<string> matname = {"SI", "KAPTON", "ALUMINUM", "Carbon", "M60J3K", "WATER"};
339  bool found = false;
340  for (string nam : matname)
341  {
342  if (material_name.find(nam) != std::string::npos)
343  {
344  m_DisplayAction->AddVolume(lv, nam);
345  if (Verbosity() >= 50)
346  {
347  cout << "SetDisplayProperty - LV " << lv->GetName() << " display with " << nam << endl;
348  }
349  found = true;
350  break;
351  }
352  }
353  if (!found)
354  {
355  m_DisplayAction->AddVolume(lv, "ANYTHING_ELSE");
356  }
357  int nDaughters = lv->GetNoDaughters();
358  for (int i = 0; i < nDaughters; ++i)
359  {
360  G4VPhysicalVolume* pv = lv->GetDaughter(i);
361 
362  // cout <<"SetDisplayProperty - PV["<<i<<"] = "<<pv->GetName()<<endl;
363 
364  G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
365  SetDisplayProperty(worldLogical);
366  }
367 }
368 
370 {
371  int active = 0;
372  for (auto& isAct : m_IsLayerActive)
373  {
374  active |= isAct;
375  }
376  if (active) // At least one layer is active
377  {
378  ostringstream geonode;
379  geonode << "CYLINDERGEOM_" << ((m_SuperDetector != "NONE") ? m_SuperDetector : m_Detector);
380  PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str().c_str());
381  if (!geo)
382  {
383  geo = new PHG4CylinderGeomContainer();
384  PHNodeIterator iter(topNode());
385  PHCompositeNode* runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
386  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str().c_str(), "PHObject");
387  runNode->addNode(newNode);
388  }
389  // here in the detector class we have internal units(mm), convert to cm
390  // before putting into the geom object
391  for (unsigned short ilayer = 0; ilayer < n_Layers; ++ilayer)
392  {
393  CylinderGeom_Mvtx* mygeom = new CylinderGeom_Mvtx(ilayer,
394  m_N_staves[ilayer],
395  m_nominal_radius[ilayer] / cm,
396  get_phistep(ilayer) / rad,
397  m_nominal_phitilt[ilayer] / rad,
398  m_nominal_phi0[ilayer] / rad
399  );
400  geo->AddLayerGeom(ilayer, mygeom);
401  } // loop per layers
402  if (Verbosity())
403  {
404  geo->identify();
405  }
406  } //is active
407 } // AddGeometryNode
408 
410 {
411  if (Verbosity() > 0)
412  {
413  cout << "-- FillPVArray --" << endl;
414  }
415  std::vector<G4VPhysicalVolume*>::iterator it = av->GetVolumesIterator();
416 
417  int nDaughters = av->TotalImprintedVolumes();
418  for (int i = 0; i < nDaughters; ++i, ++it)
419  {
420  G4VPhysicalVolume* pv = (*it);
421 
422  G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
423  // we only care about the staves, which contain the sensors, not the structures
424  if (pv->GetName().find("MVTXHalfStave_pv") != string::npos)
425  {
426  int layer = get_layer(m_StavePV.size());
427  int stave = get_stave(m_StavePV.size());
428 
429  m_StavePV.insert(make_pair(pv, make_tuple(layer, stave)));
430 
431  if (Verbosity() > 0)
432  {
433  cout << "Mvtx layer id " << layer << endl;
434  cout << "Stave in layer id " << stave << endl;
435  cout << "Mvtx stave count " << m_StavePV.size() << endl;
436  cout << "FillPVArray - AV[" << i << "] = " << (*it)->GetName() << endl;
437  cout << " LV[" << i << "] = " << worldLogical->GetName() << endl;
438  }
439 
440  FindSensor(worldLogical);
441  }
442  else // in case of stave structure
443  {
444  if (Verbosity() > 0)
445  {
446  cout << "FillPVArray - AV[" << i << "] = " << (*it)->GetName() << endl;
447  cout << " LV[" << i << "] = " << worldLogical->GetName() << endl;
448  }
449  }
450  }
451 }
452 
454 {
455  int nDaughters = lv->GetNoDaughters();
456  for (int i = 0; i < nDaughters; ++i)
457  {
458  G4VPhysicalVolume* pv = lv->GetDaughter(i);
459  if (Verbosity() > 0)
460  {
461  cout << " PV[" << i << "]: " << pv->GetName() << endl;
462  }
463  if (pv->GetName().find("MVTXSensor_") != string::npos)
464  {
465  m_SensorPV.insert(pv);
466 
467  if (Verbosity() > 0)
468  {
469  cout << " Adding Sensor Vol <" << pv->GetName() << " (" << m_SensorPV.size() << ")>" << endl;
470  }
471  }
472 
473  G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
474 
475  if (Verbosity() > 0)
476  {
477  cout << " LV[" << i << "]: " << worldLogical->GetName() << endl;
478  }
479  FindSensor(worldLogical);
480  }
481 }