ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MakeActsGeometry.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MakeActsGeometry.cc
1 
8 #include "MakeActsGeometry.h"
9 
10 #include <trackbase/TrkrDefs.h>
11 
12 #include <intt/CylinderGeomIntt.h>
13 #include <intt/InttDefs.h>
14 
15 #include <mvtx/CylinderGeom_Mvtx.h>
16 #include <mvtx/MvtxDefs.h>
17 
19 
20 #include <tpc/TpcDefs.h>
21 
23 
26 #include <g4detectors/PHG4CylinderGeom.h> // for PHG4CylinderGeom
28 
29 #include <phgeom/PHGeomUtility.h>
30 #include <phgeom/PHGeomIOTGeo.h>
31 #include <phgeom/PHGeomTGeo.h>
32 
34 #include <phool/PHCompositeNode.h>
35 #include <phool/PHDataNode.h>
36 #include <phool/PHNode.h>
37 #include <phool/PHNodeIterator.h>
38 #include <phool/PHObject.h>
39 #include <phool/getClass.h>
40 #include <phool/phool.h>
41 
50 
51 #include <ActsExamples/Detector/IBaseDetector.hpp>
52 #include <ActsExamples/EventData/Track.hpp>
53 #include <ActsExamples/Framework/AlgorithmContext.hpp>
54 #include <ActsExamples/Framework/IContextDecorator.hpp>
55 #include <ActsExamples/Framework/WhiteBoard.hpp>
56 #include <ActsExamples/Geometry/CommonGeometry.hpp>
57 #include <ActsExamples/Options/CommonOptions.hpp>
58 #include <ActsExamples/Plugins/Obj/ObjWriterOptions.hpp>
59 #include <ActsExamples/Utilities/Options.hpp>
60 
61 #include <TGeoManager.h>
62 #include <TMatrixT.h>
63 #include <TObject.h>
64 #include <TSystem.h>
65 #include <TVector3.h>
66 
67 #include <cmath>
68 #include <cstddef>
69 #include <cstdlib>
70 #include <iostream>
71 #include <map>
72 #include <memory>
73 #include <utility>
74 #include <vector>
75 
76 namespace
77 {
79  TrackingVolumePtr find_volume_by_name( const Acts::TrackingVolume* master, const std::string& name )
80  {
81  // skip if name is not composite
82  if( master->volumeName().empty() || master->volumeName()[0] != '{' ) return nullptr;
83 
84  // loop over children
85  for( const auto& child:master->confinedVolumes()->arrayObjects() )
86  {
87  if( child->volumeName() == name ) return child;
88  else if( auto found = find_volume_by_name( child.get(), name ) ) return found;
89  }
90 
91  // not found
92  return nullptr;
93  }
94 }
95 
96 MakeActsGeometry::MakeActsGeometry(const std::string &name)
97 : SubsysReco(name)
99 
101 {
103 }
104 
106 {
107 
110 
121 
122  // fill ActsSurfaceMap content
127 
128  // fill TPC volume ids
129  for( const auto& [hitsetid, surfaceVector]:m_clusterSurfaceMapTpcEdit )
130  for( const auto& surface:surfaceVector )
131  { m_surfMaps->tpcVolumeIds.insert( surface->geometryId().volume() ); }
132 
133  // fill Micromegas volume ids
134  for( const auto& [hitsetid, surface]:m_clusterSurfaceMapMmEdit )
135  { m_surfMaps->micromegasVolumeIds.insert( surface->geometryId().volume() ); }
136 
137  // print
138  if( Verbosity() )
139  {
140  for( const auto& id:m_surfMaps->tpcVolumeIds )
141  { std::cout << "MakeActsGeometry::InitRun - TPC volume id: " << id << std::endl; }
142 
143  for( const auto& id:m_surfMaps->micromegasVolumeIds )
144  { std::cout << "MakeActsGeometry::InitRun - Micromegas volume id: " << id << std::endl; }
145  }
146 
148 }
149 
151 {
153 }
155 {
157 }
158 
160 {
161 
163  // this also adds the micromegas surfaces
164  // Do this before anything else, so that the geometry is finalized
165 
166  // This should be done only on the first tracking pass, to avoid adding surfaces twice
167  if(fake_surfaces)
168  editTPCGeometry(topNode);
169  else
170  std::cout << " NOT adding TPC and MMs surfaces" << std::endl;
171 
172  // need to get nodes first, in order to be able to build the proper micromegas geometry
173  if(getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
175 
176  // In case we did not call EditTpcGeometry, we still want to make the MMs surface map
178  {
179  m_buildMMs = true;
180  std::cout << " WILL add MMs surfaces to ActsSurfaceMap " << std::endl;
181  }
182 
185 
188 
190  //makeTGeoNodeMap(topNode);
191 
193  if(Verbosity() > 3)
194  {
195  PHGeomUtility::ExportGeomtry(topNode, "sPHENIXActsGeom.root");
196  PHGeomUtility::ExportGeomtry(topNode, "sPHENIXActsGeom.gdml");
197  }
198 
200 }
201 
203 {
204  // Because we reset and rebuild the geomNode, we do edits of the TPC and Micromegas geometry in the same module
205 
206  PHGeomTGeo *geomNode = PHGeomUtility::GetGeomTGeoNode(topNode, true);
207  assert(geomNode);
208 
210  if(geomNode->isValid())
211  {
212  geomNode->Reset();
213  }
214  PHGeomIOTGeo *dstGeomIO = PHGeomUtility::GetGeomIOTGeoNode(topNode, false);
215  assert(dstGeomIO);
216  assert(dstGeomIO->isValid());
217 
218  TGeoManager *geoManager = dstGeomIO->ConstructTGeoManager();
219  geomNode->SetGeometry(geoManager);
220  assert(geoManager);
221 
222  TGeoVolume *World_vol = geoManager->GetTopVolume();
223 
224  // TPC geometry edits
225  //===============
226 
227  TGeoNode *tpc_envelope_node = nullptr;
228  TGeoNode *tpc_gas_north_node = nullptr;
229 
230  // find tpc north gas volume at path of World*/tpc_envelope*
231  if (Verbosity())
232  {
233  std::cout << "EditTPCGeometry - searching under volume: ";
234  World_vol->Print();
235  }
236  for (int i = 0; i < World_vol->GetNdaughters(); i++)
237  {
238  TString node_name = World_vol->GetNode(i)->GetName();
239 
240  if (node_name.BeginsWith("tpc_envelope"))
241  {
242  if (Verbosity())
243  std::cout << "EditTPCGeometry - found " << node_name << std::endl;
244 
245  tpc_envelope_node = World_vol->GetNode(i);
246  break;
247  }
248  }
249  assert(tpc_envelope_node);
250 
251  // find tpc north gas volume at path of World*/tpc_envelope*/tpc_gas_north*
252  TGeoVolume *tpc_envelope_vol = tpc_envelope_node->GetVolume();
253  assert(tpc_envelope_vol);
254  if (Verbosity())
255  {
256  std::cout << "EditTPCGeometry - searching under volume: ";
257  tpc_envelope_vol->Print();
258  }
259 
260  for (int i = 0; i < tpc_envelope_vol->GetNdaughters(); i++)
261  {
262  TString node_name = tpc_envelope_vol->GetNode(i)->GetName();
263 
264  if (node_name.BeginsWith("tpc_gas_north"))
265  {
266  if (Verbosity())
267  std::cout << "EditTPCGeometry - found " << node_name << std::endl;
268 
269  tpc_gas_north_node = tpc_envelope_vol->GetNode(i);
270  break;
271  }
272  }
273 
274  assert(tpc_gas_north_node);
275  TGeoVolume *tpc_gas_north_vol = tpc_gas_north_node->GetVolume();
276  assert(tpc_gas_north_vol);
277 
278  if (Verbosity())
279  {
280  std::cout << "EditTPCGeometry - gas volume: ";
281  tpc_gas_north_vol->Print();
282  }
283 
284  // adds surfaces to the underlying volume, so both north and south placements get them
285  addActsTpcSurfaces(tpc_gas_north_vol, geoManager);
286 
287  // Micromegas geometry edits
288  // These will only be made if the Micromegas nodes are found in the node tree
289  //====================
290  // The micromegas detectors have both layers in the same tile. The inner and outer sides are mirrored
291  // The detector details are (printed from Init method), where the thickness corresponds to the drift volume:
292  // layer 55: Phi segmented, radius: 82.2565 cm, thickness: 0.3 cm, zmin: -110cm, zmax: 110cm, pitch: 0.0976562cm
293  // layer: 55 volume: MICROMEGAS_55_Gas2_inner_phys
294  // layer 56: radius: z segmented, 82.6998 cm, thickness: 0.3 cm, zmin: -110cm, zmax: 110cm, pitch: 0.195312cm
295  // layer: 56 volume: MICROMEGAS_55_Gas2_outer_phys
296 
297  TGeoNode *micromegas_envelope_node = nullptr;
298  for (int i = 0; i < World_vol->GetNdaughters(); i++)
299  {
300  TString node_name = World_vol->GetNode(i)->GetName();
301 
302  if (node_name.BeginsWith("MICROMEGAS"))
303  {
304  if (Verbosity())
305  std::cout << "EditTPCGeometry - found micromegas node " << node_name << std::endl;
306 
307  micromegas_envelope_node = World_vol->GetNode(i);
308  break;
309  }
310  }
311 
312  /*
313  need to load micromegas geometry already now because it is needed for
314  defining the volumes relevant for acts
315  */
316  m_geomContainerMicromegas = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL" );
317 
318  if(micromegas_envelope_node && m_geomContainerMicromegas)
319  {
320 
322  m_buildMMs = true;
323 
324  TGeoVolume *micromegas_envelope_vol = micromegas_envelope_node->GetVolume();
325  assert(micromegas_envelope_vol);
326 
327  // Get inner and outer volume and edit them
328  for (int i = 0; i < micromegas_envelope_vol->GetNdaughters(); i++)
329  {
330  TString node_name = micromegas_envelope_vol->GetNode(i)->GetName();
331 
332  // this gets both inner and outer
333  if (node_name.BeginsWith("MICROMEGAS_55_Gas2"))
334  {
335  if (Verbosity())
336  { std::cout << "EditTPCGeometry - found Micromegas node " << node_name << std::endl; }
337 
338  auto micromegas_node = micromegas_envelope_vol->GetNode(i);
339  const int mm_layer = node_name.BeginsWith("MICROMEGAS_55_Gas2_inner") ? 0:1;
340 
341  auto micromegas_vol = micromegas_node->GetVolume();
342  assert(micromegas_vol);
343 
344  addActsMicromegasSurfaces(mm_layer, micromegas_vol, geoManager);
345  }
346  }
347  }
348 
349  // done
350  geoManager->CloseGeometry();
351 
352  // save the edited geometry to DST persistent IO node for downstream DST files
354 
355 }
356 
357 //________________________________________________________________________________
358 void MakeActsGeometry::addActsMicromegasSurfaces( int mm_layer, TGeoVolume *micromegas_vol, TGeoManager *geoManager)
359 {
360  // The input micromegas_vol is either the inner (mm_layer 0) or outer (mm_layer 1) drift volume
361 
362  // load medium
363  TGeoMedium *micromegas_medium = micromegas_vol->GetMedium();
364  assert(micromegas_medium);
365 
366  // get first layer
367  int first_layer_mm = static_cast<CylinderGeomMicromegas*>(m_geomContainerMicromegas->GetFirstLayerGeom())->get_layer();
368 
369  // load relevant geometry
370  const auto layergeom = static_cast<CylinderGeomMicromegas*>(m_geomContainerMicromegas->GetLayerGeom(first_layer_mm + mm_layer));
371 
372  // loop over micromegas tiles
373  // there will be one volume defined per tile
374  for( size_t tileid = 0; tileid < layergeom->get_tiles_count(); ++tileid )
375  {
376 
377  // get relevant tile
378  const auto& tile = layergeom->get_tile(tileid);
379 
380  // volume name
381  const auto volume_name = Form( "micromegas_measurement_%i_%zu", mm_layer, tileid );
382 
383  // create volume
384  /*
385  * in acts local coordinates, x axis is the normal to the surface
386  * y and z are the measurement directions
387  */
388  auto micromegas_measurement_vol = geoManager->MakeBox(
389  volume_name,
390  micromegas_medium,
391  (layergeom->get_thickness() - 0.1)/2,
393  tile.m_sizeZ/2 );
394 
395  micromegas_measurement_vol->SetLineColor(kBlack);
396  micromegas_measurement_vol->SetFillColor(kYellow);
397  micromegas_measurement_vol->SetVisibility(kTRUE);
398 
399  // get position of the center in global frame
400  const TVector3 global_center = layergeom->get_world_from_local_coords( tileid, { 0, 0, 0 });
401 
402  // create relevant rotation
403  const auto rotation_name = Form( "micromegas_rotation_%i_%zu", mm_layer, tileid );
404  auto rotation = new TGeoRotation(rotation_name);
405  rotation->RotateZ( tile.m_centerPhi*180./M_PI );
406  auto micromegas_measurement_location = new TGeoCombiTrans( global_center.x(), global_center.y(), global_center.z(), rotation );
407  micromegas_vol->AddNode(micromegas_measurement_vol, 1, micromegas_measurement_location);
408  }
409 }
410 
411 void MakeActsGeometry::addActsTpcSurfaces(TGeoVolume *tpc_gas_vol,
412  TGeoManager *geoManager)
413 {
414  TGeoMedium *tpc_gas_medium = tpc_gas_vol->GetMedium();
415  assert(tpc_gas_medium);
416 
417  TGeoVolume *tpc_gas_measurement_vol[m_nTpcLayers];
418  double tan_half_phi = tan(m_surfStepPhi / 2.0);
419 
420  for(unsigned int ilayer = 0; ilayer < m_nTpcLayers; ++ilayer)
421  {
422  // make a box for this layer
423  char bname[500];
424  sprintf(bname,"tpc_gas_measurement_%u",ilayer);
425 
426  // Because we use a box, not a section of a cylinder, we need this to prevent overlaps
427  // set the nominal r*phi dimension of the box so they just touch at the inner edge when placed
428  double box_r_phi = 2.0 * tan_half_phi *
429  (m_layerRadius[ilayer] - m_layerThickness[ilayer] / 2.0);
430 
431  tpc_gas_measurement_vol[ilayer] = geoManager->MakeBox(bname, tpc_gas_medium,
433  box_r_phi*half_width_clearance_phi,
435 
436  tpc_gas_measurement_vol[ilayer]->SetLineColor(kBlack);
437  tpc_gas_measurement_vol[ilayer]->SetFillColor(kYellow);
438  tpc_gas_measurement_vol[ilayer]->SetVisibility(kTRUE);
439 
440  if(Verbosity() > 30)
441  {
442  std::cout << Verbosity() << " Made box for layer " << ilayer
443  << " with dx " << m_layerThickness[ilayer] << " dy "
444  << box_r_phi << " ref arc "
445  << m_surfStepPhi * m_layerRadius[ilayer] << " dz "
446  << m_surfStepZ << std::endl;
447  tpc_gas_measurement_vol[ilayer]->Print();
448  }
449  }
450 
451  int copy = 0;
452  for (unsigned int iz = 0; iz < m_nSurfZ; ++iz)
453  {
454  // The (half) tpc gas volume is 105.5 cm long and is symmetric around (x,y,z) = (0,0,0) in its frame
455  double z_center = -105.5/2.0 + m_surfStepZ / 2.0 + (double) iz * m_surfStepZ;
456 
457  for (unsigned int imod = 0; imod < m_nTpcModulesPerLayer; ++imod)
458  {
459  for (unsigned int iphi = 0; iphi < m_nSurfPhi; ++iphi)
460  {
461 
462  double min_phi = m_modulePhiStart +
463  (double) imod * m_moduleStepPhi +
464  (double) iphi * m_surfStepPhi;
465  double phi_center = min_phi + m_surfStepPhi / 2.0;
466  double phi_center_degrees = phi_center * 180.0 / M_PI;
467 
468  for (unsigned int ilayer = 0; ilayer < m_nTpcLayers; ++ilayer)
469  {
470  copy++;
471 
472  // place copies of the gas volume to fill up the layer
473 
474  double x_center = m_layerRadius[ilayer] * cos(phi_center);
475  double y_center = m_layerRadius[ilayer] * sin(phi_center);
476 
477  char rot_name[500];
478  sprintf(rot_name,"tpc_gas_rotation_%i", copy);
479  TGeoCombiTrans *tpc_gas_measurement_location
480  = new TGeoCombiTrans(x_center, y_center, z_center,
481  new TGeoRotation(rot_name,
482  phi_center_degrees,
483  0, 0));
484 
485  tpc_gas_vol->AddNode(tpc_gas_measurement_vol[ilayer],
486  copy, tpc_gas_measurement_location);
487 
488  if(Verbosity() > 30 && ilayer == 30)
489  {
490  std::cout << " Made copy " << copy << " iz " << iz
491  << " imod " << imod << " ilayer " << ilayer
492  << " iphi " << iphi << std::endl;
493  std::cout << " x_center " << x_center
494  << " y_center " << y_center
495  << " z_center " << z_center
496  << " phi_center_degrees " << phi_center_degrees
497  << std::endl;
498  }
499  }
500  }
501  }
502  }
503 }
504 
509 {
510 
511  // define int argc and char* argv to provide options to processGeometry
512  const int argc = 20;
513  char* arg[argc];
514 
515  if(Verbosity() > 0)
516  std::cout << PHWHERE << "Magnetic field " << m_magField
517  << " with rescale " << m_magFieldRescale << std::endl;
518 
519  std::string responseFile, materialFile;
520  setMaterialResponseFile(responseFile, materialFile);
521 
522  // Response file contains arguments necessary for geometry building
523  std::string argstr[argc]{
524  "-n1", "-l0",
525  "--response-file", responseFile,
526  "--mat-input-type","file",
527  "--mat-input-file", materialFile,
528  "--bf-values","0","0", m_magField,
529  "--bf-bscalor"};
530 
531  argstr[13] = std::to_string(m_magFieldRescale);
532 
533 
535  if(m_magField.find(".root") != std::string::npos)
536  {
537  if(m_magField.find("2d") != std::string::npos)
538  {
539  m_magFieldRescale = 1;
540  }
541 
542  m_magField = std::string(getenv("CALIBRATIONROOT")) +
543  std::string("/Field/Map/sphenix3dbigmapxyz.root");
544 
545  argstr[8] = "--bf-map";
546  argstr[9] = m_magField;
547  argstr[10]= "--bf-name";
548  argstr[11] = "fieldmap";
549  argstr[12] = "--bf-lscalor";
550  argstr[13] = "10";
551  argstr[14] = "--bf-bscalor";
552  argstr[15] = std::to_string(m_magFieldRescale);
553 
554  }
555 
556  if(Verbosity() > 0)
557  std::cout << "Mag field now " << m_magField << " with rescale "
558  << m_magFieldRescale << std::endl;
559 
560  // Set vector of chars to arguments needed
561  for (int i = 0; i < argc; ++i)
562  {
563  if(Verbosity() > 0)
564  std::cout << argstr[i] << ", ";
565  // need a copy, since .c_str() returns a const char * and process geometry will not take a const
566  arg[i] = strdup(argstr[i].c_str());
567  }
568 
569  // We replicate the relevant functionality of
570  //acts/Examples/Run/Common/src/GeometryExampleBase::ProcessGeometry() in MakeActsGeometry()
571  // so we get access to the results. The layer builder magically gets the TGeoManager
572 
573  makeGeometry(argc, arg, m_detector);
574 
575  for(int i=0; i<argc; i++)
576  free(arg[i]);
577 
578 }
579 void MakeActsGeometry::setMaterialResponseFile(std::string& responseFile,
580  std::string& materialFile)
581 {
582 
583  responseFile = "tgeo-sphenix.response";
584  materialFile = "sphenix-material.json";
585  if(m_buildMMs)
586  materialFile = "sphenix-mm-material.json";
587 
589  std::ifstream file;
590 
591  file.open(responseFile);
592  if(!file)
593  {
594  if(m_buildMMs)
595  responseFile = std::string(getenv("OFFLINE_MAIN")) +
596  std::string("/share/tgeo-sphenix-mms.response");
597  else
598  responseFile = std::string(getenv("OFFLINE_MAIN")) +
599  std::string("/share/tgeo-sphenix.response");
600  }
601 
602  file.open(materialFile);
603  if(!file)
604  {
605  std::cout << materialFile
606  << " not found locally, use repo version"
607  << std::endl;
608 
609  if(m_buildMMs)
610  materialFile = std::string(getenv("CALIBRATIONROOT")) +
611  std::string("/ACTS/sphenix-mm-material.json");
612  else
613  materialFile = std::string(getenv("CALIBRATIONROOT")) +
614  std::string("/ACTS/sphenix-material.json");
615  }
616 
617  if(Verbosity() > -1)
618  {
619  std::cout << "using Acts material file : " << materialFile
620  << std::endl;
621  std::cout << "Using Acts TGeoResponse file : " << responseFile
622  << std::endl;
623  }
624 
625  return;
626 
627 }
628 void MakeActsGeometry::makeGeometry(int argc, char* argv[], ActsExamples::IBaseDetector &detector)
629 {
630 
638 
640  detector.addOptions(desc);
641  auto vm = ActsExamples::Options::parse(desc, argc, argv);
642 
644  auto geometry = ActsExamples::Geometry::build(vm, detector);
646  m_tGeometry = geometry.first;
647  m_contextDecorators = geometry.second;
648 
650 
651  size_t ievt = 0;
652  size_t ialg = 0;
653 
655  ActsExamples::WhiteBoard eventStore(
656  Acts::getDefaultLogger("EventStore#" + std::to_string(ievt),
658 
660  ActsExamples::AlgorithmContext context(ialg, ievt, eventStore);
661 
662  m_calibContext = context.calibContext;
663  m_magFieldContext = context.magFieldContext;
664  m_geoCtxt = context.geoContext;
665 
666  unpackVolumes();
667 
668  return;
669 }
670 
672 {
675  auto vol = m_tGeometry->highestTrackingVolume();
676 
677  if(Verbosity() > 10 )
678  { std::cout << "MakeActsGeometry::unpackVolumes - top volume: " << vol->volumeName() << std::endl; }
679 
680  if(m_buildMMs)
681  {
682  // micromegas
683  //std::cout << "Before: m_clusterSurfaceMapMm size " << m_clusterSurfaceMapMmEdit.size() << std::endl;
684  auto mmBarrel = find_volume_by_name( vol, "MICROMEGAS::Barrel" );
685  assert( mmBarrel );
686  makeMmMapPairs(mmBarrel);
687  //std::cout << "After: m_clusterSurfaceMapMm size " << m_clusterSurfaceMapMmEdit.size() << std::endl;
688  }
689 
690  {
691  // MVTX
692  //std::cout << "Before: Mvtx: m_clusterSurfaceMapSilicon size " << m_clusterSurfaceMapSilicon.size() << std::endl;
693  auto mvtxBarrel = find_volume_by_name( vol, "MVTX::Barrel" );
694  assert( mvtxBarrel );
695  makeMvtxMapPairs(mvtxBarrel);
696  //std::cout << "After: Mvtx: m_clusterSurfaceMapSilicon size " << m_clusterSurfaceMapSilicon.size() << std::endl;
697  }
698 
699  {
700  // INTT
701  //std::cout << "Before: INTT: m_clusterSurfaceMapSilicon size " << m_clusterSurfaceMapSilicon.size() << std::endl;
702  auto inttBarrel = find_volume_by_name( vol, "Silicon::Barrel" );
703  assert( inttBarrel );
704  makeInttMapPairs(inttBarrel);
705  //std::cout << "After: INTT: m_clusterSurfaceMapSilicon size " << m_clusterSurfaceMapSilicon.size() << std::endl;
706  }
707 
708  {
709  // TPC
710  //std::cout << "Before: m_clusterSurfaceMapTpc size " << m_clusterSurfaceMapTpcEdit.size() << std::endl;
711  auto tpcBarrel = find_volume_by_name( vol, "TPC::Barrel" );
712  assert( tpcBarrel );
713  makeTpcMapPairs(tpcBarrel);
714  //std::cout << "After: m_clusterSurfaceMapTpc size " << m_clusterSurfaceMapTpcEdit.size() << std::endl;
715  }
716 
717  return;
718 }
719 
721 {
722  if(Verbosity() > 10)
723  { std::cout << "MakeActsGeometry::makeTpcMapPairs - tpcVolume: " << tpcVolume->volumeName() << std::endl; }
724 
725  auto tpcLayerArray = tpcVolume->confinedLayers();
726  auto tpcLayerVector = tpcLayerArray->arrayObjects();
727 
729  for(unsigned int i = 0; i < tpcLayerVector.size(); i++)
730  {
731  auto surfaceArray = tpcLayerVector.at(i)->surfaceArray();
732  if(surfaceArray == NULL){
733  continue;
734  }
737  auto surfaceVector = surfaceArray->surfaces();
738  for( unsigned int j = 0; j < surfaceVector.size(); j++)
739  {
740  auto surf = surfaceVector.at(j)->getSharedPtr();
741  auto vec3d = surf->center(m_geoCtxt);
742 
744  std::vector<double> world_center = {vec3d(0) / 10.0,
745  vec3d(1) / 10.0,
746  vec3d(2) / 10.0};
747 
749  unsigned int layer = TrkrDefs::getLayer(hitsetkey);
750 
753  //std::map<TrkrDefs::hitsetkey, std::vector<Surface>>::iterator mapIter;
754  std::map<unsigned int, std::vector<Surface>>::iterator mapIter;
755  //mapIter = m_clusterSurfaceMapTpcEdit.find(hitsetkey);
756  mapIter = m_clusterSurfaceMapTpcEdit.find(layer);
757 
758  if(mapIter != m_clusterSurfaceMapTpcEdit.end())
759  {
760  mapIter->second.push_back(surf);
761  }
762  else
763  {
765  std::vector<Surface> dumvec;
766  dumvec.push_back(surf);
767  std::pair<unsigned int, std::vector<Surface>> tmp =
768  std::make_pair(layer, dumvec);
769  m_clusterSurfaceMapTpcEdit.insert(tmp);
770  }
771 
772  }
773  }
774 
775 }
776 
777 //____________________________________________________________________________________________
779 {
780  if(Verbosity() > 10)
781  { std::cout << "MakeActsGeometry::makeMmMapPairs - mmVolume: " << mmVolume->volumeName() << std::endl; }
782  const auto mmLayerArray = mmVolume->confinedLayers();
783  const auto mmLayerVector = mmLayerArray->arrayObjects();
784 
786  for(unsigned int i = 0; i < mmLayerVector.size(); i++)
787  {
788  auto surfaceArray = mmLayerVector.at(i)->surfaceArray();
789  if(!surfaceArray) continue;
790 
793  const auto surfaceVector = surfaceArray->surfaces();
794  for( unsigned int j = 0; j < surfaceVector.size(); j++)
795  {
796  auto surface = surfaceVector.at(j)->getSharedPtr();
797  auto vec3d = surface->center(m_geoCtxt);
798 
800  TVector3 world_center(
801  vec3d(0)/Acts::UnitConstants::cm,
802  vec3d(1)/Acts::UnitConstants::cm,
803  vec3d(2)/Acts::UnitConstants::cm
804  );
805 
806  // get relevant layer
807  int layer = -1;
808  CylinderGeomMicromegas* layergeom = nullptr;
809  const auto range = m_geomContainerMicromegas->get_begin_end();
810  for( auto iter = range.first; iter != range.second; ++iter )
811  {
812  auto this_layergeom = static_cast<CylinderGeomMicromegas*>( iter->second );
813  if(this_layergeom->check_radius(world_center))
814  {
815  layer = iter->first;
816  layergeom = this_layergeom;
817  break;
818  }
819  }
820 
821  if( !layergeom )
822  {
823  std::cout << "MakeActsGeometry::makeMmMapPairs - could not file CylinderGeomMicromegas matching ACTS surface" << std::endl;
824  continue;
825  }
826 
827  // get matching tile
828  int tileid = layergeom->find_tile_planar( world_center );
829  if( tileid < 0 )
830  {
831  std::cout << "MakeActsGeometry::makeMmMapPairs - could not file Micromegas tile matching ACTS surface" << std::endl;
832  continue;
833  }
834 
835  // get segmentation type
836  const auto segmentation_type = layergeom->get_segmentation_type();
837 
838  // create hitset key and insert surface in map
839  const auto hitsetkey = MicromegasDefs::genHitSetKey(layer, segmentation_type, tileid);
840  const auto [iter, inserted] = m_clusterSurfaceMapMmEdit.insert( std::make_pair( hitsetkey, surface ) );
841  assert( inserted );
842  }
843  }
844 }
845 
847 {
848 
849  if(Verbosity() > 10)
850  { std::cout << "MakeActsGeometry::makeInttMapPairs - inttVolume: " << inttVolume->volumeName() << std::endl; }
851 
852  auto inttLayerArray = inttVolume->confinedLayers();
853 
854  auto inttLayerVector = inttLayerArray->arrayObjects();
855  // inttLayerVector is a std::vector<LayerPtr>
856  for (unsigned int i = 0; i < inttLayerVector.size(); i++)
857  {
858  // Get the ACTS::SurfaceArray from each INTT LayerPtr being iterated over
859  auto surfaceArray = inttLayerVector.at(i)->surfaceArray();
860  if (surfaceArray == NULL)
861  continue;
862 
863  // surfaceVector is an Acts::SurfaceVector, vector of surfaces
864  auto surfaceVector = surfaceArray->surfaces();
865 
866  for (unsigned int j = 0; j < surfaceVector.size(); j++)
867  {
868  auto surf = surfaceVector.at(j)->getSharedPtr();
869  auto vec3d = surf->center(m_geoCtxt);
870 
871  double ref_rad[4] = {7.188, 7.732, 9.680, 10.262};
872 
873  std::vector<double> world_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0}; // convert from mm to cm
874 
878  double layer_rad = sqrt(pow(world_center[0], 2) + pow(world_center[1], 2));
879 
880  unsigned int layer = 0;
881  for (unsigned int i = 0; i < 4; ++i)
882  {
883  if (fabs(layer_rad - ref_rad[i]) < 0.1)
884  layer = i + 3;
885  }
886 
888 
889  // Add this surface to the map
890  std::pair<TrkrDefs::hitsetkey, Surface> tmp = make_pair(hitsetkey, surf);
891  m_clusterSurfaceMapSilicon.insert(tmp);
892 
893  if (Verbosity() > 10)
894  {
895  // check it is in there
896  unsigned int ladderPhi = InttDefs::getLadderPhiId(hitsetkey);
897  unsigned int ladderZ = InttDefs::getLadderZId(hitsetkey);
898 
899  std::cout << "Layer radius " << layer_rad << " layer " << layer << " ladderPhi " << ladderPhi << " ladderZ " << ladderZ
900  << " recover surface from m_clusterSurfaceMapSilicon " << std::endl;
901  std::cout << " surface type " << surf->type() << std::endl;
902  auto assoc_layer = surf->associatedLayer();
903  std::cout << std::endl
904  << " Layer type " << assoc_layer->layerType() << std::endl;
905 
906  auto assoc_det_element = surf->associatedDetectorElement();
907  if (assoc_det_element != nullptr)
908  {
909  std::cout << " Associated detElement has non-null pointer "
910  << assoc_det_element << std::endl;
911  std::cout << std::endl
912  << " Associated detElement found, thickness = "
913  << assoc_det_element->thickness() << std::endl;
914  }
915  else
916  std::cout << std::endl
917  << " Associated detElement is nullptr " << std::endl;
918  }
919  }
920  }
921 
922 }
923 
925 {
926 
927  if(Verbosity() > 10)
928  { std::cout << "MakeActsGeometry::makeMvtxMapPairs - mvtxVolume: " << mvtxVolume->volumeName() << std::endl; }
929 
930  // Now get the LayerArrays corresponding to each volume
931  auto mvtxBarrelLayerArray = mvtxVolume->confinedLayers(); // the barrel is all we care about
932 
933  // This is a BinnedArray<LayerPtr>, but we care only about the sensitive surfaces
934  auto mvtxLayerVector1 = mvtxBarrelLayerArray->arrayObjects(); // the barrel is all we care about
935 
936  // mvtxLayerVector has size 3, but only index 1 has any surfaces since the clayersplits option was removed
937  // the layer number has to be deduced from the radius
938  for (unsigned int i = 0; i < mvtxLayerVector1.size(); ++i)
939  {
940  // Get the Acts::SurfaceArray from each MVTX LayerPtr being iterated over
941  auto surfaceArray = mvtxLayerVector1.at(i)->surfaceArray();
942  if (surfaceArray == NULL)
943  continue;
944 
945  double ref_rad[3] = {2.556, 3.359, 4.134};
946 
947  // surfaceVector is an Acts::SurfaceVector, vector of surfaces
948  //std::vector<const Surface*>
949  auto surfaceVector = surfaceArray->surfaces();
950  for (unsigned int j = 0; j < surfaceVector.size(); j++)
951  {
952  auto surf = surfaceVector.at(j)->getSharedPtr();
953  auto vec3d = surf->center(m_geoCtxt);
954  std::vector<double> world_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0}; // convert from mm to cm
955  double layer_rad = sqrt(pow(world_center[0], 2) + pow(world_center[1], 2));
956  unsigned int layer = 0;
957  for (unsigned int i = 0; i < 3; ++i)
958  {
959  if (fabs(layer_rad - ref_rad[i]) < 0.1)
960  layer = i;
961  }
962 
964 
965  // Add this surface to the map
966  std::pair<TrkrDefs::hitsetkey, Surface> tmp = make_pair(hitsetkey, surf);
967 
968  m_clusterSurfaceMapSilicon.insert(tmp);
969 
970  if (Verbosity() > 10)
971  {
972  unsigned int stave = MvtxDefs::getStaveId(hitsetkey);
973  unsigned int chip = MvtxDefs::getChipId(hitsetkey);
974 
975  // check it is in there
976  std::cout << "Layer radius " << layer_rad << " Layer "
977  << layer << " stave " << stave << " chip " << chip
978  << " recover surface from m_clusterSurfaceMapSilicon "
979  << std::endl;
980 
981  std::map<TrkrDefs::hitsetkey, Surface>::iterator surf_iter = m_clusterSurfaceMapSilicon.find(hitsetkey);
982  std::cout << " surface type " << surf_iter->second->type()
983  << std::endl;
984  auto assoc_layer = surf->associatedLayer();
985  std::cout << std::endl << " Layer type "
986  << assoc_layer->layerType() << std::endl;
987 
988  auto assoc_det_element = surf->associatedDetectorElement();
989  if (assoc_det_element != nullptr)
990  {
991  std::cout << " Associated detElement has non-null pointer "
992  << assoc_det_element << std::endl;
993  std::cout << std::endl
994  << " Associated detElement found, thickness = "
995  << assoc_det_element->thickness() << std::endl;
996  }
997  else
998  std::cout << std::endl
999  << " Associated detElement is nullptr " << std::endl;
1000  }
1001  }
1002  }
1003 }
1004 
1006 {
1007  // Look up TPC surface index values from world position of surface center
1008  // layer
1009  unsigned int layer = 999;
1010  double layer_rad = sqrt(pow(world[0],2) + pow(world[1],2));
1011  for(unsigned int ilayer=0;ilayer<m_nTpcLayers;++ilayer)
1012  {
1013  double tpc_ref_radius_low =
1014  m_layerRadius[ilayer] - m_layerThickness[ilayer] / 2.0;
1015  double tpc_ref_radius_high =
1016  m_layerRadius[ilayer] + m_layerThickness[ilayer] / 2.0;
1017  if(layer_rad >= tpc_ref_radius_low && layer_rad < tpc_ref_radius_high)
1018  {
1019  layer = ilayer;
1020  break;
1021  }
1022  }
1023 
1024  if(layer >= m_nTpcLayers)
1025  {
1026  std::cout << PHWHERE
1027  << "Error: undefined layer, do nothing world = "
1028  << world[0] << " " << world[1] << " " << world[2]
1029  << " layer " << layer << std::endl;
1031  }
1032 
1033  layer += 7;
1034 
1035  // side - from world z sign
1036  unsigned int side;
1037  if(world[2] < 0)
1038  side = 0;
1039  else
1040  side = 1;
1041 
1042  // readout module
1043  unsigned int readout_mod = 999;
1044  double phi_world = atan2(world[1], world[0]);
1045  for(unsigned int imod=0; imod<m_nTpcModulesPerLayer; ++imod)
1046  {
1047  double min_phi = m_modulePhiStart + (double) imod * m_moduleStepPhi;
1048  double max_phi = m_modulePhiStart + (double) (imod+1) * m_moduleStepPhi;
1049  if(phi_world >=min_phi && phi_world < max_phi)
1050  {
1051  readout_mod = imod;
1052  break;
1053  }
1054  }
1055  if(readout_mod >= m_nTpcModulesPerLayer)
1056  {
1057  std::cout << PHWHERE
1058  << "Error: readout_mod is undefined, do nothing phi_world = "
1059  << phi_world << std::endl;
1061  }
1062 
1063  TrkrDefs::hitsetkey hitset_key = TpcDefs::genHitSetKey(layer, readout_mod, side);
1064  if(Verbosity() > 3)
1065  if(layer == 30)
1066  std::cout << " world = " << world[0] << " " << world[1]
1067  << " " << world[2] << " phi_world "
1068  << phi_world*180/3.14159 << " layer " << layer
1069  << " readout_mod " << readout_mod << " side " << side
1070  << " hitsetkey " << hitset_key<< std::endl;
1071 
1072  return hitset_key;
1073 }
1074 
1076 {
1077  // Look up the MVTX sensor index values from the world position of the surface center
1078 
1079  CylinderGeom_Mvtx *layergeom = dynamic_cast<CylinderGeom_Mvtx *>(m_geomContainerMvtx->GetLayerGeom(layer));
1080  if (!layergeom)
1081  {
1082  std::cout << PHWHERE << "Did not get layergeom for layer "
1083  << layer << std::endl;
1084  return 0;
1085  }
1086 
1087  unsigned int stave = 0;
1088  unsigned int chip = 0;
1089  layergeom->get_sensor_indices_from_world_coords(world, stave, chip);
1090 
1091  double check_pos[3] = {0, 0, 0};
1092  layergeom->find_sensor_center(stave, 0, 0, chip, check_pos);
1093 
1094  TrkrDefs::hitsetkey mvtx_hitsetkey = MvtxDefs::genHitSetKey(layer, stave, chip);
1095 
1096  return mvtx_hitsetkey;
1097 }
1098 
1100 {
1101  // Look up the INTT sensor index values from the world position of the surface center
1102 
1103  CylinderGeomIntt *layergeom = dynamic_cast<CylinderGeomIntt *>(m_geomContainerIntt->GetLayerGeom(layer));
1104  if (!layergeom)
1105  {
1106  std::cout << PHWHERE << "Did not get layergeom for layer "
1107  << layer << std::endl;
1108  return 0;
1109  }
1110 
1111  double location[3] = {world[0], world[1], world[2]};
1112  int segment_z_bin = 0;
1113  int segment_phi_bin = 0;
1114  layergeom->find_indices_from_segment_center(segment_z_bin,
1115  segment_phi_bin, location);
1116 
1117  double check_pos[3] = {0, 0, 0};
1118  layergeom->find_segment_center(segment_z_bin, segment_phi_bin, check_pos);
1119 
1120  TrkrDefs::hitsetkey intt_hitsetkey = InttDefs::genHitSetKey(layer, segment_z_bin, segment_phi_bin);
1121 
1122  return intt_hitsetkey;
1123 }
1124 
1126 {
1127  // This is just a diagnostic method
1128  // it lets you list all of the nodes by setting print_sensors = true
1129 
1130  if (!m_geoManager)
1131  {
1132  std::cout << PHWHERE << " Did not find TGeoManager, quit! " << std::endl;
1133  return;
1134  }
1135  TGeoVolume *topVol = m_geoManager->GetTopVolume();
1136  TObjArray *nodeArray = topVol->GetNodes();
1137 
1138  TIter iObj(nodeArray);
1139  while (TObject *obj = iObj())
1140  {
1141  TGeoNode *node = dynamic_cast<TGeoNode *>(obj);
1142  std::string node_str = node->GetName();
1143 
1144  std::string mvtx("MVTX_Wrapper");
1145  std::string intt("ladder");
1146  std::string intt_ext("ladderext");
1147  std::string tpc("tpc_envelope");
1148  std::string micromegas("MICROMEGAS_55");
1149 
1150  if (node_str.compare(0, mvtx.length(), mvtx) == 0) // is it in the MVTX?
1151  {
1152  if (Verbosity() > 2)
1153  std::cout << " node " << node->GetName() << " is the MVTX wrapper"
1154  << std::endl;
1155 
1157  TObjArray *mvtxArray = node->GetNodes();
1158  TIter mvtxObj(mvtxArray);
1159  while(TObject *mvtx = mvtxObj())
1160  {
1161  TGeoNode *mvtxNode = dynamic_cast<TGeoNode *>(mvtx);
1162  if(Verbosity() > 2)
1163  std::cout << "mvtx node name is " << mvtxNode->GetName()
1164  << std::endl;
1165  std::string mvtxav1("av_1");
1166  std::string mvtxNodeName = mvtxNode->GetName();
1167 
1169  if(mvtxNodeName.compare(0, mvtxav1.length(), mvtxav1) == 0)
1170  getMvtxKeyFromNode(mvtxNode);
1171  }
1172  }
1173  else if (node_str.compare(0, intt.length(), intt) == 0) // is it in the INTT?
1174  {
1175  // We do not want the "ladderext" nodes
1176  if (node_str.compare(0, intt_ext.length(), intt_ext) == 0)
1177  continue;
1178 
1179  if (Verbosity() > 2)
1180  std::cout << " node " << node->GetName() << " is in the INTT"
1181  << std::endl;
1182  getInttKeyFromNode(node);
1183  }
1186  else if (node_str.compare(0, tpc.length(), tpc) == 0) // is it in the TPC?
1187  {
1188  if(Verbosity() > 2)
1189  std::cout << " node " << node->GetName()
1190  << " is in the TPC " << std::endl;
1191  }
1192  else if (node_str.compare(0, micromegas.length(), micromegas) == 0) // is it in the Micromegas?
1193  {
1194  if(Verbosity() > 2)
1195  std::cout << " node " << node->GetName()
1196  << " is in the MMs " << std::endl;
1197  }
1198  else
1199  continue;
1200 
1201  bool print_sensor_paths = false; // normally false
1202  if (print_sensor_paths)
1203  {
1204  // Descends the node tree to find the active silicon nodes - used for information only
1205  std::cout << " Top Node is " << node->GetName() << " volume name is " << node->GetVolume()->GetName() << std::endl;
1206 
1207  int nmax_print = 20;
1208  isActive(node, nmax_print);
1209  }
1210  }
1211 }
1212 
1214 {
1215  int layer = -1; // sPHENIX layer number
1216  int itype = -1; // specifies inner (0) or outer (1) sensor
1217  int ladder_phi = -1; // copy number of ladder in phi
1218  int zposneg = -1; // specifies positive (1) or negative (0) z
1219  int ladder_z = -1; // 0-3, from most negative z to most positive
1220 
1221  std::string s = gnode->GetName();
1222  std::string delimiter = "_";
1223  std::string posz("posz");
1224  std::string negz("negz");
1225 
1226  size_t pos = 0;
1227  std::string token;
1228 
1229  int counter = 0;
1230  while ((pos = s.find(delimiter)) != std::string::npos)
1231  {
1232  token = s.substr(0, pos);
1233 
1234  s.erase(0, pos + delimiter.length());
1235  if (counter == 1)
1236  layer = std::atoi(token.c_str()) + 3;
1237  if (counter == 2)
1238  itype = std::atoi(token.c_str());
1239  if (counter == 3)
1240  {
1241  ladder_phi = std::atoi(token.c_str());
1242  if (s.compare(0, negz.length(), negz) == 0) zposneg = 0;
1243  if (s.compare(0, posz.length(), posz) == 0) zposneg = 1;
1244  }
1245  counter++;
1246  }
1247 
1248  ladder_z = itype + zposneg * 2;
1249 
1250  // The active sensor is a daughter of gnode
1251  int ndaught = gnode->GetNdaughters();
1252  if (ndaught == 0)
1253  {
1254  std::cout << PHWHERE << "OOPS: Did not find INTT sensor! Quit."
1255  << std::endl;
1256  exit(1);
1257  }
1258 
1259  std::string intt_refactive("siactive");
1260  TGeoNode *sensor_node = 0;
1261  for (int i = 0; i < ndaught; ++i)
1262  {
1263  std::string node_str = gnode->GetDaughter(i)->GetName();
1264 
1265  if (node_str.compare(0, intt_refactive.length(), intt_refactive) == 0)
1266  {
1267  sensor_node = gnode->GetDaughter(i);
1268  break;
1269  }
1270  }
1271 
1272  // unique key identifying this sensor
1273  TrkrDefs::hitsetkey node_key = InttDefs::genHitSetKey(layer, ladder_z,
1274  ladder_phi);
1275 
1276  std::pair<TrkrDefs::hitsetkey, TGeoNode *> tmp = std::make_pair(
1277  node_key, sensor_node);
1278  m_clusterNodeMap.insert(tmp);
1279 
1280  if (Verbosity() > 1)
1281  std::cout << " INTT layer " << layer << " ladder_phi " << ladder_phi
1282  << " itype " << itype << " zposneg " << zposneg
1283  << " ladder_z " << ladder_z << " name "
1284  << sensor_node->GetName() << std::endl;
1285 
1286  return;
1287 }
1288 void MakeActsGeometry::getTpcKeyFromNode(TGeoNode * /*gnode*/)
1289 {
1290 
1291 }
1293 {
1294  int counter = 0;
1295  int impr = -1; // stave number, 1-48 in TGeo
1296  int layer = -1;
1297  int stave = -1; // derived from impr
1298  int chip = -1; // 9 chips per stave
1299 
1300  std::string s = gnode->GetName();
1301  std::string delimiter = "_";
1302 
1303  size_t pos = 0;
1304  std::string token;
1305 
1306  while ((pos = s.find(delimiter)) != std::string::npos)
1307  {
1308  token = s.substr(0, pos);
1309  //std::cout << token << std::endl;
1310  s.erase(0, pos + delimiter.length());
1311  if (counter == 3)
1312  impr = std::atoi(token.c_str());
1313 
1314  counter++;
1315  }
1316 
1317  // extract layer and stave info from impr
1318  // int staves_in_layer[3] = {12, 16, 20};
1319  // note - impr stave count starts from 1, not 0, but TrkrCluster counting starts from 0, so we reduce it by 1 here
1320  impr -= 1;
1321 
1322  if (impr < 12)
1323  {
1324  layer = 0;
1325  stave = impr;
1326  }
1327  else if (impr > 11 && impr < 28)
1328  {
1329  layer = 1;
1330  stave = impr - 12;
1331  }
1332  else
1333  {
1334  layer = 2;
1335  stave = impr - 28;
1336  }
1337 
1338  // Now descend node tree to find chip ID's - there are multiple chips per stave
1339  TGeoNode *module_node = gnode->GetDaughter(0);
1340  int mnd = module_node->GetNdaughters();
1341  std::string mvtx_chip("MVTXChip");
1342  for (int i = 0; i < mnd; ++i)
1343  {
1344  std::string dstr = module_node->GetDaughter(i)->GetName();
1345  if (dstr.compare(0, mvtx_chip.length(), mvtx_chip) == 0)
1346  {
1347  if (Verbosity() > 3)
1348  std::cout << "Found MVTX layer " << layer << " stave " << stave
1349  << " chip " << i << " with node name "
1350  << module_node->GetDaughter(i)->GetName() << std::endl;
1351 
1352  // Make key for this chip
1353  TrkrDefs::hitsetkey node_key = MvtxDefs::genHitSetKey(layer,
1354  stave, i);
1355 
1356  // add sensor node to map
1357  TGeoNode *sensor_node = module_node->GetDaughter(i)->GetDaughter(0);
1358  std::pair<TrkrDefs::hitsetkey, TGeoNode *> tmp = std::make_pair(
1359  node_key, sensor_node);
1360  m_clusterNodeMap.insert(tmp);
1361 
1362  if (Verbosity() > 3)
1363  std::cout << " MVTX layer " << layer << " stave " << stave
1364  << " chip " << chip << " name "
1365  << sensor_node->GetName() << std::endl;
1366  }
1367  }
1368 
1369  return;
1370 }
1371 
1372 void MakeActsGeometry::isActive(TGeoNode *gnode, int nmax_print)
1373 {
1374  // Not used in analysis: diagnostic, for looking at the node tree only.
1375  // Recursively searches gnode for silicon sensors, prints out heirarchy
1376 
1377  std::string node_str = gnode->GetName();
1378 
1379  std::string intt_refactive("siactive");
1380  std::string mvtx_refactive("MVTXSensor");
1381  std::string tpc_refactive("tpc_gas_measurement");
1382  std::string micromegas_refactive("MICROMEGAS_55");
1383 
1384  if (node_str.compare(0, intt_refactive.length(), intt_refactive) == 0)
1385  {
1386  std::cout << " ******* Found INTT active volume, node is "
1387  << gnode->GetName()
1388  << " volume name is " << gnode->GetVolume()->GetName()
1389  << std::endl;
1390 
1391  return;
1392  }
1393  else if (node_str.compare(0, mvtx_refactive.length(), mvtx_refactive) == 0)
1394  {
1395  std::cout << " ******* Found MVTX active volume, node is "
1396  << gnode->GetName()
1397  << " volume name is " << gnode->GetVolume()->GetName()
1398  << std::endl;
1399 
1400  return;
1401  }
1402  else if (node_str.compare(0, tpc_refactive.length(), tpc_refactive) == 0)
1403  {
1404  if(nprint_tpc < nmax_print)
1405  {
1406  std::cout << " ******* Found TPC active volume, node is "
1407  << gnode->GetName()
1408  << " volume name is " << gnode->GetVolume()->GetName()
1409  << std::endl;
1410  }
1411  nprint_tpc++;
1412 
1413  return;
1414  }
1415  else if (node_str.compare(0, micromegas_refactive.length(), micromegas_refactive) == 0)
1416  {
1417  std::cout << " ******* Found Micromegas active volume, node is "
1418  << gnode->GetName()
1419  << " volume name is " << gnode->GetVolume()->GetName()
1420  << std::endl;
1421 
1422  return;
1423  }
1424  else
1425  {
1426  if(nprint_tpc < nmax_print)
1427  {
1428  std::cout << " ******* Found node "
1429  << gnode->GetName()
1430  << " volume name is " << gnode->GetVolume()->GetName()
1431  << std::endl;
1432  }
1433  nprint_tpc++;
1434 
1435  return;
1436  }
1437 
1438 
1439  int ndaught = gnode->GetNdaughters();
1440 
1441  for (int i = 0; i < ndaught; ++i)
1442  {
1443  std::cout << " " << gnode->GetVolume()->GetName()
1444  << " daughter " << i << " has name "
1445  << gnode->GetDaughter(i)->GetVolume()->GetName() << std::endl;
1446 
1447  isActive(gnode->GetDaughter(i), nmax_print);
1448  }
1449 }
1450 
1451 
1453 {
1456  m_surfStepZ = (m_maxSurfZ - m_minSurfZ) / (double) m_nSurfZ;
1457  m_moduleStepPhi = 2.0 * M_PI / 12.0;
1459  m_surfStepPhi = 2.0 * M_PI / (double) (m_nSurfPhi * m_nTpcModulesPerLayer);
1460  for(unsigned int isector = 0; isector < 3; ++isector)
1461  {
1462  layer_thickness_sector[isector] =
1463  (m_maxRadius[isector] - m_minRadius[isector]) / 16.0;
1464 
1465  for(unsigned int ilayer =0; ilayer < 16; ++ilayer)
1466  {
1467  m_layerRadius[isector*16 + ilayer] =
1468  m_minRadius[isector] + layer_thickness_sector[isector] *
1469  (double) ilayer + layer_thickness_sector[isector] / 2.0;
1470 
1471  m_layerThickness[isector*16 + ilayer] =
1472  layer_thickness_sector[isector];
1473  }
1474  }
1475 }
1476 
1478 {
1479  PHNodeIterator iter(topNode);
1480 
1482  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1483 
1485  if (!dstNode)
1486  {
1487  std::cerr << "DST Node missing, quitting" << std::endl;
1488  throw std::runtime_error("failed to find DST node in PHActsSourceLinks::createNodes");
1489  }
1490 
1492  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
1493 
1495  if (!svtxNode)
1496  {
1497  svtxNode = new PHCompositeNode("SVTX");
1498  dstNode->addNode(svtxNode);
1499  }
1500 
1501  m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
1502  "ActsSurfaceMaps");
1503  if(!m_surfMaps)
1504  {
1505  m_surfMaps = new ActsSurfaceMaps();
1507  = new PHDataNode<ActsSurfaceMaps>(m_surfMaps, "ActsSurfaceMaps");
1508  svtxNode->addNode(node);
1509 
1510  }
1511 
1512  m_actsGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
1513  "ActsTrackingGeometry");
1514  if(!m_actsGeometry)
1515  {
1518  = new PHDataNode<ActsTrackingGeometry>(m_actsGeometry, "ActsTrackingGeometry");
1519  svtxNode->addNode(tGeoNode);
1520  }
1521 
1523 }
1524 
1525 /*
1526  * GetNodes():
1527  * Get all the all the required nodes off the node tree
1528  */
1530 {
1532  if (!m_geoManager)
1533  {
1534  std::cout << PHWHERE << " Did not find TGeoManager, quit! "
1535  << std::endl;
1537  }
1538 
1540  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
1541  if (!m_geomContainerMvtx)
1542  {
1543  std::cout << PHWHERE
1544  << " CYLINDERGEOM_MVTX node not found on node tree"
1545  << std::endl;
1547  }
1548 
1550  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1551  if (!m_geomContainerTpc)
1552  {
1553  std::cout << PHWHERE
1554  << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX"
1555  << std::endl;
1557  }
1558 
1560  PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
1561  if (!m_geomContainerIntt)
1562  {
1563  std::cout << PHWHERE
1564  << " CYLINDERGEOM_INTT node not found on node tree"
1565  << std::endl;
1567  }
1568 
1569  // load micromegas geometry
1570  // do not abort if not found
1571  m_geomContainerMicromegas = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
1573  {
1574  std::cout << PHWHERE
1575  << " CYLINDERGEOM_MICROMEGAS_FULL node not found on node tree"
1576  << std::endl;
1577  }
1578 
1580 }