ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4MicromegasDetector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4MicromegasDetector.cc
1 
8 
9 #include <phparameter/PHParameters.h>
10 
12 
13 #include <g4main/PHG4Detector.h>
14 
16 
17 #include <phool/getClass.h>
18 #include <phool/PHCompositeNode.h>
19 #include <phool/PHIODataNode.h>
20 #include <phool/PHNode.h> // for PHNode
21 #include <phool/PHNodeIterator.h> // for PHNodeIterator
22 #include <phool/PHObject.h> // for PHObject
23 #include <phool/recoConsts.h>
24 
25 #include <Geant4/G4Tubs.hh>
26 #include <Geant4/G4Color.hh>
27 #include <Geant4/G4LogicalVolume.hh>
28 #include <Geant4/G4Material.hh>
29 #include <Geant4/G4PVPlacement.hh>
30 #include <Geant4/G4SystemOfUnits.hh>
31 #include <Geant4/G4VisAttributes.hh>
32 #include <Geant4/G4String.hh> // for G4String
33 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
34 #include <Geant4/G4Types.hh> // for G4double
35 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
36 #include <Geant4/G4VSolid.hh> // for G4VSolid
37 
38 #include <cmath>
39 #include <iostream>
40 #include <numeric>
41 #include <tuple> // for make_tuple, tuple
42 #include <utility> // for pair, make_pair
43 #include <vector> // for vector
44 
45 //____________________________________________________________________________..
47  : PHG4Detector(subsys, Node, dnam)
48  , m_Params(parameters)
49 {}
50 
51 //_______________________________________________________________
53 { return m_activeVolumes.find( volume ) != m_activeVolumes.end(); }
54 
55 //_______________________________________________________________
57 {
58  const auto iter = m_activeVolumes.find( volume );
59  return iter == m_activeVolumes.end() ? -1:iter->second;
60 }
61 
62 //_______________________________________________________________
64 {
66  construct_micromegas(logicWorld);
68 }
69 
70 //_______________________________________________________________
71 void PHG4MicromegasDetector::Print(const std::string &what) const
72 {
73  std::cout << "PHG4Micromegas Detector:" << std::endl;
74  if (what == "ALL" || what == "VOLUME")
75  {
76  std::cout << "Version 0.1" << std::endl;
77  std::cout << "Parameters:" << std::endl;
78  m_Params->Print();
79  }
80  return;
81 }
82 
83 //_______________________________________________________________
85 {
86  // get the list of NIST materials
87  // ---------------------------------
88  auto G4_N = GetDetectorMaterial("G4_N");
89  auto G4_O = GetDetectorMaterial("G4_O");
90  auto G4_C = GetDetectorMaterial("G4_C");
91  auto G4_H = GetDetectorMaterial("G4_H");
92  auto G4_Si = GetDetectorMaterial("G4_Si");
93  auto G4_Ar = GetDetectorMaterial("G4_Ar");
94  auto G4_Cr = GetDetectorMaterial("G4_Cr");
95  auto G4_Fe = GetDetectorMaterial("G4_Fe");
96  auto G4_Mn = GetDetectorMaterial("G4_Mn");
97  auto G4_Ni = GetDetectorMaterial("G4_Ni");
98  auto G4_Cu = GetDetectorMaterial("G4_Cu");
99 
100  // combine elements
101  // ----------------
102  static const G4double temperature = 298.15*kelvin;
103  static const G4double pressure = 1.*atmosphere;
104 
105  // FR4
106  if (!GetDetectorMaterial("mmg_FR4", false))
107  {
108  auto mmg_FR4 = new G4Material( "mmg_FR4", 1.860*g/cm3, 4, kStateSolid);
109  mmg_FR4->AddMaterial( G4_C, 0.43550 );
110  mmg_FR4->AddMaterial( G4_H, 0.03650 );
111  mmg_FR4->AddMaterial( G4_O, 0.28120 );
112  mmg_FR4->AddMaterial( G4_Si, 0.24680 );
113  }
114 
115  // Kapton
116  if (!GetDetectorMaterial("mmg_Kapton", false))
117  {
118  auto mmg_Kapton = new G4Material( "mmg_Kapton", 1.420*g/cm3, 4, kStateSolid);
119  mmg_Kapton->AddMaterial( G4_C, 0.6911330 );
120  mmg_Kapton->AddMaterial( G4_H, 0.0263620 );
121  mmg_Kapton->AddMaterial( G4_N, 0.0732700 );
122  mmg_Kapton->AddMaterial( G4_O, 0.2092350);
123  }
124 
125  // MMgas
126  if (!GetDetectorMaterial("mmg_Gas", false))
127  {
128  auto mmg_Gas = new G4Material( "mmg_Gas", 0.00170335*g/cm3, 3, kStateGas, temperature, pressure);
129  mmg_Gas->AddMaterial( G4_Ar, 0.900 );
130  mmg_Gas->AddMaterial( G4_C, 0.0826586 );
131  mmg_Gas->AddMaterial( G4_H, 0.0173414 );
132  }
133 
134  // MMMesh
135  if (!GetDetectorMaterial("mmg_Mesh", false))
136  {
137  auto mmg_Mesh = new G4Material( "mmg_Mesh", 2.8548*g/cm3, 5, kStateSolid);
138  mmg_Mesh->AddMaterial( G4_Cr, 0.1900 );
139  mmg_Mesh->AddMaterial( G4_Fe, 0.6800 );
140  mmg_Mesh->AddMaterial( G4_Mn, 0.0200 );
141  mmg_Mesh->AddMaterial( G4_Ni, 0.1000 );
142  mmg_Mesh->AddMaterial( G4_Si, 0.0100 );
143  }
144 
145  // MMStrips
146  if (!GetDetectorMaterial("mmg_Strips", false))
147  { new G4Material( "mmg_Strips", 5.248414*g/cm3, G4_Cu, kStateSolid); }
148 
149  // MMResistivePaste
150  if (!GetDetectorMaterial("mmg_ResistPaste", false))
151  { new G4Material( "mmg_ResistPaste", 0.77906*g/cm3, G4_C, kStateSolid); }
152 
153 }
154 
155 //_______________________________________________________________
157 {
158  // components enumeration
159  /*
160  this describes all the detector onion layers for a single side
161  note that the detector is two sided
162  */
163  enum class Component
164  {
165  PCB,
166  CuStrips,
167  KaptonStrips,
168  ResistiveStrips,
169  Gas1,
170  Mesh,
171  Gas2,
172  DriftCuElectrode,
173  DriftKapton,
174  DriftCarbon
175  };
176 
177  // layer thickness
178  // numbers from M. Vandenbroucke <maxence.vandenbroucke@cea.fr>
179  const std::map<Component,float> layer_thickness =
180  {
181  { Component::PCB, 1.*mm },
182  { Component::CuStrips, 12.*micrometer },
183  { Component::KaptonStrips, 50.*micrometer },
184  { Component::ResistiveStrips, 20.*micrometer },
185  { Component::Gas1, 120.*micrometer },
186  { Component::Mesh, 18.*0.8*micrometer }, // 0.8 correction factor is to account for the mesh denstity@18/45
187  { Component::Gas2, 3.*mm },
188  { Component::DriftCuElectrode, 15.*micrometer },
189  { Component::DriftKapton, 50.*micrometer },
190  { Component::DriftCarbon, 1.*mm }
191  };
192 
193  // materials
194  const std::map<Component,G4Material*> layer_material =
195  {
196  { Component::PCB, GetDetectorMaterial("mmg_FR4") },
197  { Component::CuStrips, GetDetectorMaterial("mmg_Strips") },
198  { Component::KaptonStrips, GetDetectorMaterial("mmg_Kapton") },
199  { Component::ResistiveStrips, GetDetectorMaterial("mmg_ResistPaste" ) },
200  { Component::Gas1, GetDetectorMaterial( "mmg_Gas" ) },
201  { Component::Mesh, GetDetectorMaterial("mmg_Mesh") },
202  { Component::Gas2, GetDetectorMaterial( "mmg_Gas" ) },
203  { Component::DriftCuElectrode, GetDetectorMaterial("G4_Cu") },
204  { Component::DriftKapton, GetDetectorMaterial("mmg_Kapton") },
205  { Component::DriftCarbon, GetDetectorMaterial("G4_C") }
206  };
207 
208  // color
209  const std::map<Component, G4Colour> layer_color =
210  {
211  { Component::PCB, G4Colour::Green()},
212  { Component::CuStrips, G4Colour::Brown()},
213  { Component::KaptonStrips, G4Colour::Brown()},
214  { Component::ResistiveStrips, G4Colour::Black()},
215  { Component::Gas1, G4Colour::Grey()},
216  { Component::Mesh, G4Colour::White()},
217  { Component::Gas2, G4Colour::Grey()},
218  { Component::DriftCuElectrode, G4Colour::Brown()},
219  { Component::DriftKapton, G4Colour::Brown()},
220  { Component::DriftCarbon, G4Colour(150/255., 75/255., 0)}
221  };
222 
223  // setup layers in the correct order, going outwards from beam axis
224  /* same compoment can appear multiple times. Layer names must be unique */
225  using LayerDefinition = std::tuple<Component,std::string>;
226  const std::vector<LayerDefinition> layer_definitions =
227  {
228  // inner side
229  std::make_tuple( Component::DriftCarbon, "DriftCarbon_inner" ),
230  std::make_tuple( Component::DriftKapton, "DriftKapton_inner" ),
231  std::make_tuple( Component::DriftCuElectrode, "DriftCuElectrode_inner" ),
232  std::make_tuple( Component::Gas2, "Gas2_inner" ),
233  std::make_tuple( Component::Mesh, "Mesh_inner" ),
234  std::make_tuple( Component::Gas1, "Gas1_inner" ),
235  std::make_tuple( Component::ResistiveStrips, "ResistiveStrips_inner" ),
236  std::make_tuple( Component::KaptonStrips, "KaptonStrips_inner" ),
237  std::make_tuple( Component::CuStrips, "CuStrips_inner" ),
238 
239  // PCB
240  std::make_tuple( Component::PCB, "PCB" ),
241 
242  // outer side (= inner side, mirrored)
243  std::make_tuple( Component::CuStrips, "CuStrips_outer" ),
244  std::make_tuple( Component::KaptonStrips, "KaptonStrips_outer" ),
245  std::make_tuple( Component::ResistiveStrips, "ResistiveStrips_outer" ),
246  std::make_tuple( Component::Gas1, "Gas1_outer" ),
247  std::make_tuple( Component::Mesh, "Mesh_outer" ),
248  std::make_tuple( Component::Gas2, "Gas2_outer" ),
249  std::make_tuple( Component::DriftCuElectrode, "DriftCuElectrode_outer" ),
250  std::make_tuple( Component::DriftKapton, "DriftKapton_outer" ),
251  std::make_tuple( Component::DriftCarbon, "DriftCarbon_outer" )
252  };
253 
254  // start seting up volumes
255  // get initial radius
256  const double radius = m_Params->get_double_param("mm_radius")*cm;
257  const double length = m_Params->get_double_param("mm_length")*cm;
258 
259  // get total thickness
260  const double thickness = std::accumulate(
261  layer_definitions.begin(), layer_definitions.end(), 0.,
262  [layer_thickness](double value, LayerDefinition layer )
263  { return value + layer_thickness.at(std::get<0>(layer)); } );
264 
265  std::cout << "PHG4MicromegasDetector::ConstructMe - detector thickness is " << thickness/cm << " cm" << std::endl;
266 
267  // create mother volume
268  auto cylinder_solid = new G4Tubs( G4String(GetName()), radius - 0.001*mm, radius + thickness + 0.001*mm, length / 2., 0, M_PI*2);
269 
271  auto cylinder_logic = new G4LogicalVolume( cylinder_solid, GetDetectorMaterial(rc->get_StringFlag("WorldMaterial")), G4String(GetName()) );
272  auto vis = new G4VisAttributes(G4Color(G4Colour::Grey()));
273  vis->SetForceSolid(true);
274  vis->SetVisibility(false);
275  cylinder_logic->SetVisAttributes(vis);
276 
277  // add placement
278  new G4PVPlacement( nullptr, G4ThreeVector(0,0,0), cylinder_logic, G4String(GetName()), logicWorld, false, 0, OverlapCheck() );
279 
280  // keep track of current layer
281  int layer_index = m_FirstLayer;
282 
283  // create detector
284  /* we loop over registered layers and create volumes for each */
285  auto current_radius = radius;
286  for( const auto& layer:layer_definitions )
287  {
288  const Component& type = std::get<0>(layer);
289  const std::string& name = std::get<1>(layer);
290 
291  // layer name
292  G4String cname = G4String(GetName()) + "_" + name;
293 
294  // get thickness, material and name
295  const auto thickness = layer_thickness.at(type);
296  const auto material = layer_material.at(type);
297  const auto color = layer_color.at(type);
298 
299  auto component_solid = new G4Tubs(cname+"_solid", current_radius, current_radius+thickness, length/2, 0, M_PI*2);
300  auto component_logic = new G4LogicalVolume( component_solid, material, cname+"_logic");
301  auto vis = new G4VisAttributes( color );
302  vis->SetForceSolid(true);
303  vis->SetVisibility(true);
304  component_logic->SetVisAttributes(vis);
305 
306  auto component_phys = new G4PVPlacement( nullptr, G4ThreeVector(0,0,0), component_logic, cname+"_phys", cylinder_logic, false, 0, OverlapCheck() );
307 
308  // store active volume
309  if( type == Component::Gas2 ) m_activeVolumes.insert( std::make_pair( component_phys, layer_index++ ) );
310  else m_passiveVolumes.insert( component_phys );
311 
312  // update radius
313  current_radius += thickness;
314  }
315 
316  // print physical layers
317  std::cout << "PHG4MicromegasDetector::ConstructMe - first layer: " << m_FirstLayer << std::endl;
318  for( const auto& pair:m_activeVolumes )
319  { std::cout << "PHG4MicromegasDetector::ConstructMe - layer: " << pair.second << " volume: " << pair.first->GetName() << std::endl; }
320 
321  return;
322 }
323 
324 //_______________________________________________________________
326 {
327  // do nothing if detector is inactive
328  if( !m_Params->get_int_param("active")) return;
329 
330  // find or create geometry node
331  std::string geonode_name = std::string( "CYLINDERGEOM_" ) + m_SuperDetector;
332  auto geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode_name);
333  if (!geonode)
334  {
335  geonode = new PHG4CylinderGeomContainer();
336  PHNodeIterator iter(topNode());
337  auto runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
338  auto newNode = new PHIODataNode<PHObject>(geonode, geonode_name, "PHObject");
339  runNode->addNode(newNode);
340  }
341 
342  // add cylinder objects
343  /* one cylinder is added per physical volume. The dimention correspond to the drift volume */
344  for( const auto& pair:m_activeVolumes )
345  {
346  // store layer and volume
347  const int layer = pair.second;
348  const G4VPhysicalVolume* volume_phys = pair.first;
349 
350  // get solid volume, cast to a tube
351  const auto tub = dynamic_cast<const G4Tubs*>( volume_phys->GetLogicalVolume()->GetSolid() );
352 
353  // create cylinder and match geometry
354  /* note: cylinder segmentation type and pitch is set in PHG4MicromegasHitReco */
355  auto cylinder = new CylinderGeomMicromegas(layer);
356  cylinder->set_radius( (tub->GetInnerRadius()/cm + tub->GetOuterRadius()/cm)/2 );
357  cylinder->set_thickness( tub->GetOuterRadius()/cm - tub->GetInnerRadius()/cm );
358  cylinder->set_zmin( -tub->GetZHalfLength()/cm );
359  cylinder->set_zmax( tub->GetZHalfLength()/cm );
360  geonode->AddLayerGeom(layer, cylinder);
361  }
362 
363 }