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