EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JsonGeometryConverter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JsonGeometryConverter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
32 
33 #include <cstdio>
34 #include <fstream>
35 #include <iostream>
36 #include <map>
37 #include <sstream>
38 #include <stdexcept>
39 #include <string>
40 
41 #include <boost/algorithm/string.hpp>
42 #include <boost/algorithm/string/finder.hpp>
43 #include <boost/algorithm/string/iter_find.hpp>
44 
45 namespace {
46 
47 using json = nlohmann::json;
48 
49 // helper functions to encode/decode indefinite material
50 //
51 // encoded either as `null` for vacuum or to an array of material parameters
52 
53 json encodeMaterial(const Acts::Material& material) {
54  if (!material) {
55  return nullptr;
56  }
57  json encoded = json::array();
58  for (unsigned i = 0; i < material.parameters().size(); ++i) {
59  encoded.push_back(material.parameters()[i]);
60  }
61  return encoded;
62 }
63 
64 Acts::Material decodeMaterial(const json& encoded) {
65  if (encoded.is_null()) {
66  return {};
67  }
69  Acts::Material::ParametersVector::Zero();
70  for (auto i = params.size(); 0 < i--;) {
71  // .at(...) ensures bound checks
72  params[i] = encoded.at(i);
73  }
74  return Acts::Material(params);
75 }
76 
77 // helper functions to encode/decode concrete material slabs
78 //
79 // encoded as an object w/ two entries: `material` and `thickness`
80 
81 json encodeMaterialSlab(const Acts::MaterialSlab& slab) {
82  return {
83  {"material", encodeMaterial(slab.material())},
84  {"thickness", slab.thickness()},
85  };
86 }
87 
88 Acts::MaterialSlab decodeMaterialSlab(const json& encoded) {
89  return Acts::MaterialSlab(decodeMaterial(encoded.at("material")),
90  encoded.at("thickness").get<float>());
91 }
92 
93 } // namespace
94 
97  : m_cfg(std::move(cfg)) {
98  // Validate the configuration
99  if (!m_cfg.logger) {
100  throw std::invalid_argument("Missing logger");
101  }
102 }
103 
104 std::pair<std::map<Acts::GeometryIdentifier,
105  std::shared_ptr<const Acts::ISurfaceMaterial>>,
106  std::map<Acts::GeometryIdentifier,
107  std::shared_ptr<const Acts::IVolumeMaterial>>>
109  auto& j = materialmaps;
110  // The return maps
111  std::pair<SurfaceMaterialMap, VolumeMaterialMap> maps;
112  ACTS_VERBOSE("j2a: Reading material maps from json file.");
113  ACTS_VERBOSE("j2a: Found entries for " << j.count(m_cfg.volkey)
114  << " volume(s).");
115  // Structured binding
116  for (auto& [key, value] : j.items()) {
117  // Check if this the volume key
118  if (key == m_cfg.volkey) {
119  // Get the volume json
120  auto volj = value;
121  for (auto& [vkey, vvalue] : volj.items()) {
122  // Create the volume id
123  int vid = std::stoi(vkey);
124  Acts::GeometryIdentifier volumeID;
125  volumeID.setVolume(vid);
126  ACTS_VERBOSE("j2a: -> Found Volume " << vid);
127  // Loop through the information in the volume
128  for (auto& [vckey, vcvalue] : vvalue.items()) {
129  if (vckey == m_cfg.boukey and m_cfg.processBoundaries and
130  not vcvalue.empty()) {
131  ACTS_VERBOSE("j2a: --> BoundarySurface(s) to be parsed");
132  for (auto& [bkey, bvalue] : vcvalue.items()) {
133  // Create the boundary id
134  int bid = std::stoi(bkey);
135  Acts::GeometryIdentifier boundaryID(volumeID);
136  boundaryID.setBoundary(bid);
137  ACTS_VERBOSE("j2a: ---> Found boundary surface " << bid);
138  if (bvalue[m_cfg.mapkey] == true) {
139  auto boumat = jsonToSurfaceMaterial(bvalue);
140  maps.first[boundaryID] =
141  std::shared_ptr<const ISurfaceMaterial>(boumat);
142  }
143  }
144  } else if (vckey == m_cfg.laykey) {
145  ACTS_VERBOSE("j2a: --> Layer(s) to be parsed");
146  // Loop over layers and repeat
147  auto layj = vcvalue;
148  for (auto& [lkey, lvalue] : layj.items()) {
149  // Create the layer id
150  int lid = std::stoi(lkey);
151  Acts::GeometryIdentifier layerID(volumeID);
152  layerID.setLayer(lid);
153  ACTS_VERBOSE("j2a: ---> Found Layer " << lid);
154  // Finally loop over layer components
155  for (auto& [lckey, lcvalue] : lvalue.items()) {
156  if (lckey == m_cfg.repkey and m_cfg.processRepresenting and
157  not lcvalue.empty()) {
158  ACTS_VERBOSE("j2a: ----> Found representing surface");
159  if (lcvalue[m_cfg.mapkey] == true) {
160  auto repmat = jsonToSurfaceMaterial(lcvalue);
161  maps.first[layerID] =
162  std::shared_ptr<const Acts::ISurfaceMaterial>(repmat);
163  }
164  } else if (lckey == m_cfg.appkey and m_cfg.processApproaches and
165  not lcvalue.empty()) {
166  ACTS_VERBOSE("j2a: ----> Found approach surface(s)");
167  // Loop over approach surfaces
168  for (auto& [askey, asvalue] : lcvalue.items()) {
169  // Create the layer id, todo set to max value
170  int aid = (askey == "*") ? 0 : std::stoi(askey);
171  Acts::GeometryIdentifier approachID(layerID);
172  approachID.setApproach(aid);
173  ACTS_VERBOSE("j2a: -----> Approach surface " << askey);
174  if (asvalue[m_cfg.mapkey] == true) {
175  auto appmat = jsonToSurfaceMaterial(asvalue);
176  maps.first[approachID] =
177  std::shared_ptr<const Acts::ISurfaceMaterial>(appmat);
178  }
179  }
180  } else if (lckey == m_cfg.senkey and m_cfg.processSensitives and
181  not lcvalue.empty()) {
182  ACTS_VERBOSE("j2a: ----> Found sensitive surface(s)");
183  // Loop over sensitive surfaces
184  for (auto& [sskey, ssvalue] : lcvalue.items()) {
185  // Create the layer id, todo set to max value
186  int sid = (sskey == "*") ? 0 : std::stoi(sskey);
187  Acts::GeometryIdentifier senisitiveID(layerID);
188  senisitiveID.setSensitive(sid);
189  ACTS_VERBOSE("j2a: -----> Sensitive surface " << sskey);
190  if (ssvalue[m_cfg.mapkey] == true) {
191  auto senmat = jsonToSurfaceMaterial(ssvalue);
192  maps.first[senisitiveID] =
193  std::shared_ptr<const Acts::ISurfaceMaterial>(senmat);
194  }
195  }
196  }
197  }
198  }
199 
200  } else if (m_cfg.processVolumes and vckey == m_cfg.matkey and
201  not vcvalue.empty()) {
202  ACTS_VERBOSE("--> VolumeMaterial to be parsed");
203  if (vcvalue[m_cfg.mapkey] == true) {
204  auto intermat = jsonToVolumeMaterial(vcvalue);
205  maps.second[volumeID] =
206  std::shared_ptr<const Acts::IVolumeMaterial>(intermat);
207  }
208  }
209  }
210  }
211  } else if (key == m_cfg.geoversion) {
212  ACTS_VERBOSE("Detector version: " << m_cfg.geoversion);
213  }
214  }
215 
216  // Return the filled maps
217  return maps;
218 }
219 
223  const DetectorMaterialMaps& maps) {
224  DetectorRep detRep;
225  // Collect all GeometryIdentifiers per VolumeID for the formatted output
226  for (auto& [key, value] : maps.first) {
227  geo_id_value vid = key.volume();
228  auto volRep = detRep.volumes.find(vid);
229  if (volRep == detRep.volumes.end()) {
230  detRep.volumes.insert({vid, VolumeRep()});
231  volRep = detRep.volumes.find(vid);
232  volRep->second.volumeID = key;
233  }
234  geo_id_value lid = key.layer();
235  if (lid != 0) {
236  // we are on a layer, get the layer rep
237  auto layRep = volRep->second.layers.find(lid);
238  if (layRep == volRep->second.layers.end()) {
239  volRep->second.layers.insert({lid, LayerRep()});
240  layRep = volRep->second.layers.find(lid);
241  layRep->second.layerID = key;
242  }
243  // now insert appropriately
244  geo_id_value sid = key.sensitive();
245  geo_id_value aid = key.approach();
246  if (sid != 0) {
247  layRep->second.sensitives.insert({sid, value.get()});
248  } else if (aid != 0) {
249  layRep->second.approaches.insert({aid, value.get()});
250  } else {
251  layRep->second.representing = value.get();
252  }
253 
254  } else {
255  // not on a layer can only be a boundary surface
256  geo_id_value bid = key.boundary();
257  volRep->second.boundaries.insert({bid, value.get()});
258  }
259  }
260  for (auto& [key, value] : maps.second) {
261  // find the volume representation
262  geo_id_value vid = key.volume();
263  auto volRep = detRep.volumes.find(vid);
264  if (volRep == detRep.volumes.end()) {
265  detRep.volumes.insert({vid, VolumeRep()});
266  volRep = detRep.volumes.find(vid);
267  volRep->second.volumeID = key;
268  }
269  volRep->second.material = value.get();
270  }
271  // convert the detector representation to json format
272  return detectorRepToJson(detRep);
273 }
274 
277  json detectorj;
278  ACTS_VERBOSE("a2j: Writing json from detector representation");
279  ACTS_VERBOSE("a2j: Found entries for " << detRep.volumes.size()
280  << " volume(s).");
281 
282  json volumesj;
283  for (auto& [key, value] : detRep.volumes) {
284  json volj;
285  ACTS_VERBOSE("a2j: -> Writing Volume: " << key);
286  volj[m_cfg.namekey] = value.volumeName;
287  std::ostringstream svolumeID;
288  svolumeID << value.volumeID;
289  volj[m_cfg.geometryidkey] = svolumeID.str();
290  if (m_cfg.processVolumes && value.material) {
291  volj[m_cfg.matkey] = volumeMaterialToJson(*value.material);
292  }
293  // Write the layers
294  if (not value.layers.empty()) {
295  ACTS_VERBOSE("a2j: ---> Found " << value.layers.size() << " layer(s) ");
296  json layersj;
297  for (auto& [lkey, lvalue] : value.layers) {
298  ACTS_VERBOSE("a2j: ----> Convert layer " << lkey);
299  json layj;
300  std::ostringstream slayerID;
301  slayerID << lvalue.layerID;
302  layj[m_cfg.geometryidkey] = slayerID.str();
303  // First check for approaches
304  if (not lvalue.approaches.empty() and m_cfg.processApproaches) {
305  ACTS_VERBOSE("a2j: -----> Found " << lvalue.approaches.size()
306  << " approach surface(s)");
307  json approachesj;
308  for (auto& [akey, avalue] : lvalue.approaches) {
309  ACTS_VERBOSE("a2j: ------> Convert approach surface " << akey);
310  approachesj[std::to_string(akey)] = surfaceMaterialToJson(*avalue);
311  if (lvalue.approacheSurfaces.find(akey) !=
312  lvalue.approacheSurfaces.end())
313  addSurfaceToJson(approachesj[std::to_string(akey)],
314  lvalue.approacheSurfaces.at(akey));
315  }
316  // Add to the layer json
317  layj[m_cfg.appkey] = approachesj;
318  }
319  // Then check for sensitive
320  if (not lvalue.sensitives.empty() and m_cfg.processSensitives) {
321  ACTS_VERBOSE("a2j: -----> Found " << lvalue.sensitives.size()
322  << " sensitive surface(s)");
323  json sensitivesj;
324  for (auto& [skey, svalue] : lvalue.sensitives) {
325  ACTS_VERBOSE("a2j: ------> Convert sensitive surface " << skey);
326  sensitivesj[std::to_string(skey)] = surfaceMaterialToJson(*svalue);
327  if (lvalue.sensitiveSurfaces.find(skey) !=
328  lvalue.sensitiveSurfaces.end())
329  addSurfaceToJson(sensitivesj[std::to_string(skey)],
330  lvalue.sensitiveSurfaces.at(skey));
331  }
332  // Add to the layer json
333  layj[m_cfg.senkey] = sensitivesj;
334  }
335  // Finally check for representing
336  if (lvalue.representing != nullptr and m_cfg.processRepresenting) {
337  ACTS_VERBOSE("a2j: ------> Convert representing surface ");
338  layj[m_cfg.repkey] = surfaceMaterialToJson(*lvalue.representing);
339  if (lvalue.representingSurface != nullptr)
340  addSurfaceToJson(layj[m_cfg.repkey], lvalue.representingSurface);
341  }
342  layersj[std::to_string(lkey)] = layj;
343  }
344  volj[m_cfg.laykey] = layersj;
345  }
346  // Write the boundary surfaces
347  if (not value.boundaries.empty()) {
348  ACTS_VERBOSE("a2j: ---> Found " << value.boundaries.size()
349  << " boundary/ies ");
350  json boundariesj;
351  for (auto& [bkey, bvalue] : value.boundaries) {
352  ACTS_VERBOSE("a2j: ----> Convert boundary " << bkey);
353  boundariesj[std::to_string(bkey)] = surfaceMaterialToJson(*bvalue);
354  if (value.boundarySurfaces.find(bkey) != value.boundarySurfaces.end())
355  addSurfaceToJson(boundariesj[std::to_string(bkey)],
356  value.boundarySurfaces.at(bkey));
357  }
358  volj[m_cfg.boukey] = boundariesj;
359  }
360 
361  volumesj[std::to_string(key)] = volj;
362  }
363  // Assign the volume json to the detector json
364  detectorj[m_cfg.volkey] = volumesj;
365 
366  return detectorj;
367 }
368 
372  Acts::ISurfaceMaterial* sMaterial = nullptr;
373  // The bin utility for deescribing the data
374  Acts::BinUtility bUtility;
375  for (auto& [key, value] : material.items()) {
376  if (key == m_cfg.transfokeys and not value.empty()) {
377  bUtility = Acts::BinUtility(jsonToTransform(value));
378  break;
379  }
380  }
381  // Convert the material
382  Acts::MaterialSlabMatrix mpMatrix;
383  // Structured binding
384  for (auto& [key, value] : material.items()) {
385  // Check json keys
386  if (key == m_cfg.bin0key and not value.empty()) {
387  bUtility += jsonToBinUtility(value);
388  } else if (key == m_cfg.bin1key and not value.empty()) {
389  bUtility += jsonToBinUtility(value);
390  }
391  if (key == m_cfg.datakey and not value.empty()) {
392  mpMatrix = jsonToMaterialMatrix(value);
393  }
394  }
395 
396  // We have protoMaterial
397  if (mpMatrix.empty()) {
398  sMaterial = new Acts::ProtoSurfaceMaterial(bUtility);
399  } else if (bUtility.bins() == 1) {
400  sMaterial = new Acts::HomogeneousSurfaceMaterial(mpMatrix[0][0]);
401  } else {
402  sMaterial = new Acts::BinnedSurfaceMaterial(bUtility, mpMatrix);
403  }
404  // return what you have
405  return sMaterial;
406 }
407 
410  const json& material) {
411  Acts::IVolumeMaterial* vMaterial = nullptr;
412  // The bin utility for deescribing the data
413  Acts::BinUtility bUtility;
414  for (auto& [key, value] : material.items()) {
415  if (key == m_cfg.transfokeys and not value.empty()) {
416  bUtility = Acts::BinUtility(jsonToTransform(value));
417  break;
418  }
419  }
420  // Convert the material
421  std::vector<Material> mmat;
422  // Structured binding
423  for (auto& [key, value] : material.items()) {
424  // Check json keys
425  if (key == m_cfg.bin0key and not value.empty()) {
426  bUtility += jsonToBinUtility(value);
427  } else if (key == m_cfg.bin1key and not value.empty()) {
428  bUtility += jsonToBinUtility(value);
429  } else if (key == m_cfg.bin2key and not value.empty()) {
430  bUtility += jsonToBinUtility(value);
431  }
432  if (key == m_cfg.datakey and not value.empty()) {
433  for (const auto& bin : value) {
434  mmat.push_back(decodeMaterial(bin));
435  }
436  }
437  }
438 
439  // We have protoMaterial
440  if (mmat.empty()) {
441  vMaterial = new Acts::ProtoVolumeMaterial(bUtility);
442  } else if (mmat.size() == 1) {
443  vMaterial = new Acts::HomogeneousVolumeMaterial(mmat[0]);
444  } else {
445  if (bUtility.dimensions() == 2) {
446  std::function<Acts::Vector2D(Acts::Vector3D)> transfoGlobalToLocal;
447  Acts::Grid2D grid = createGrid2D(bUtility, transfoGlobalToLocal);
448 
451  Acts::Grid2D::index_t nBins = grid.numLocalBins();
452 
453  Acts::EAxis axis1(min[0], max[0], nBins[0]);
454  Acts::EAxis axis2(min[1], max[1], nBins[1]);
455 
456  // Build the grid and fill it with data
457  MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
458 
459  for (size_t bin = 0; bin < mmat.size(); bin++) {
460  mGrid.at(bin) = mmat[bin].parameters();
461  }
462  MaterialMapper<MaterialGrid2D> matMap(transfoGlobalToLocal, mGrid);
463  vMaterial =
465  std::move(matMap), bUtility);
466  } else if (bUtility.dimensions() == 3) {
467  std::function<Acts::Vector3D(Acts::Vector3D)> transfoGlobalToLocal;
468  Acts::Grid3D grid = createGrid3D(bUtility, transfoGlobalToLocal);
469 
472  Acts::Grid3D::index_t nBins = grid.numLocalBins();
473 
474  Acts::EAxis axis1(min[0], max[0], nBins[0]);
475  Acts::EAxis axis2(min[1], max[1], nBins[1]);
476  Acts::EAxis axis3(min[2], max[2], nBins[2]);
477 
478  // Build the grid and fill it with data
479  MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
480 
481  for (size_t bin = 0; bin < mmat.size(); bin++) {
482  mGrid.at(bin) = mmat[bin].parameters();
483  }
484  MaterialMapper<MaterialGrid3D> matMap(transfoGlobalToLocal, mGrid);
485  vMaterial =
487  std::move(matMap), bUtility);
488  }
489  }
490  // return what you have
491  return vMaterial;
492 }
493 
496  DetectorRep detRep;
497  convertToRep(detRep, *tGeometry.highestTrackingVolume());
498  return detectorRepToJson(detRep);
499 }
500 
502  DetectorRep& detRep, const Acts::TrackingVolume& tVolume) {
503  // The writer reader volume representation
504  VolumeRep volRep;
505  volRep.volumeName = tVolume.volumeName();
506  // there are confined volumes
507  if (tVolume.confinedVolumes() != nullptr) {
508  // get through the volumes
509  auto& volumes = tVolume.confinedVolumes()->arrayObjects();
510  // loop over the volumes
511  for (auto& vol : volumes) {
512  // recursive call
513  convertToRep(detRep, *vol);
514  }
515  }
516  // there are dense volumes
517  if (m_cfg.processDenseVolumes && !tVolume.denseVolumes().empty()) {
518  // loop over the volumes
519  for (auto& vol : tVolume.denseVolumes()) {
520  // recursive call
521  convertToRep(detRep, *vol);
522  }
523  }
524  // Get the volume Id
525  Acts::GeometryIdentifier volumeID = tVolume.geometryId();
526  geo_id_value vid = volumeID.volume();
527 
528  // Write the material if there's one
529  if (tVolume.volumeMaterial() != nullptr) {
530  volRep.material = tVolume.volumeMaterial();
531  } else if (m_cfg.processnonmaterial == true) {
532  Acts::BinUtility bUtility = DefaultBin(tVolume);
533  Acts::IVolumeMaterial* bMaterial = new Acts::ProtoVolumeMaterial(bUtility);
534  volRep.material = bMaterial;
535  }
536  // there are confied layers
537  if (tVolume.confinedLayers() != nullptr) {
538  // get the layers
539  auto& layers = tVolume.confinedLayers()->arrayObjects();
540  // loop of the volumes
541  for (auto& lay : layers) {
542  auto layRep = convertToRep(*lay);
543  if (layRep) {
544  // it's a valid representation so let's go with it
545  Acts::GeometryIdentifier layerID = lay->geometryId();
546  geo_id_value lid = layerID.layer();
547  volRep.layers.insert({lid, std::move(layRep)});
548  }
549  }
550  }
551  // Let's finally check the boundaries
552  for (auto& bsurf : tVolume.boundarySurfaces()) {
553  // the surface representation
554  auto& bssfRep = bsurf->surfaceRepresentation();
555  if (bssfRep.surfaceMaterial() != nullptr) {
556  Acts::GeometryIdentifier boundaryID = bssfRep.geometryId();
557  geo_id_value bid = boundaryID.boundary();
558  // Ignore if the volumeID is not correct (i.e. shared boundary)
559  // if (boundaryID.value(Acts::GeometryIdentifier::volume_mask) == vid){
560  volRep.boundaries[bid] = bssfRep.surfaceMaterial();
561  volRep.boundarySurfaces[bid] = &bssfRep;
562  // }
563  } else if (m_cfg.processnonmaterial == true) {
564  // if no material suface exist add a default one for the mapping
565  // configuration
566  Acts::GeometryIdentifier boundaryID = bssfRep.geometryId();
567  geo_id_value bid = boundaryID.boundary();
568  Acts::BinUtility bUtility = DefaultBin(bssfRep);
569  Acts::ISurfaceMaterial* bMaterial =
570  new Acts::ProtoSurfaceMaterial(bUtility);
571  volRep.boundaries[bid] = bMaterial;
572  volRep.boundarySurfaces[bid] = &bssfRep;
573  }
574  }
575  // Write if it's good
576  if (volRep) {
577  volRep.volumeName = tVolume.volumeName();
578  volRep.volumeID = volumeID;
579  detRep.volumes.insert({vid, std::move(volRep)});
580  }
581  return;
582 }
583 
585  const Acts::Layer& tLayer) {
586  LayerRep layRep;
587  // fill layer ID information
588  layRep.layerID = tLayer.geometryId();
589  if (m_cfg.processSensitives and tLayer.surfaceArray() != nullptr) {
590  for (auto& ssf : tLayer.surfaceArray()->surfaces()) {
591  if (ssf != nullptr && ssf->surfaceMaterial() != nullptr) {
592  Acts::GeometryIdentifier sensitiveID = ssf->geometryId();
593  geo_id_value sid = sensitiveID.sensitive();
594  layRep.sensitives.insert({sid, ssf->surfaceMaterial()});
595  layRep.sensitiveSurfaces.insert({sid, ssf});
596  } else if (m_cfg.processnonmaterial == true) {
597  // if no material suface exist add a default one for the mapping
598  // configuration
599  Acts::GeometryIdentifier sensitiveID = ssf->geometryId();
600  geo_id_value sid = sensitiveID.sensitive();
601  Acts::BinUtility sUtility = DefaultBin(*ssf);
602  Acts::ISurfaceMaterial* sMaterial =
603  new Acts::ProtoSurfaceMaterial(sUtility);
604  layRep.sensitives.insert({sid, sMaterial});
605  layRep.sensitiveSurfaces.insert({sid, ssf});
606  }
607  }
608  }
609  // the representing
610  if (!(tLayer.surfaceRepresentation().geometryId() == GeometryIdentifier())) {
611  if (tLayer.surfaceRepresentation().surfaceMaterial() != nullptr) {
613  layRep.representingSurface = &tLayer.surfaceRepresentation();
614  } else if (m_cfg.processnonmaterial == true) {
615  // if no material suface exist add a default one for the mapping
616  // configuration
617  Acts::BinUtility rUtility = DefaultBin(tLayer.surfaceRepresentation());
618  Acts::ISurfaceMaterial* rMaterial =
619  new Acts::ProtoSurfaceMaterial(rUtility);
620  layRep.representing = rMaterial;
621  layRep.representingSurface = &tLayer.surfaceRepresentation();
622  }
623  }
624  // the approach
625  if (tLayer.approachDescriptor() != nullptr) {
626  for (auto& asf : tLayer.approachDescriptor()->containedSurfaces()) {
627  // get the surface and check for material
628  if (asf->surfaceMaterial() != nullptr) {
629  Acts::GeometryIdentifier approachID = asf->geometryId();
630  geo_id_value aid = approachID.approach();
631  layRep.approaches.insert({aid, asf->surfaceMaterial()});
632  layRep.approacheSurfaces.insert({aid, asf});
633  } else if (m_cfg.processnonmaterial == true) {
634  // if no material suface exist add a default one for the mapping
635  // configuration
636  Acts::GeometryIdentifier approachID = asf->geometryId();
637  geo_id_value aid = approachID.approach();
638  Acts::BinUtility aUtility = DefaultBin(*asf);
639  Acts::ISurfaceMaterial* aMaterial =
640  new Acts::ProtoSurfaceMaterial(aUtility);
641  layRep.approaches.insert({aid, aMaterial});
642  layRep.approacheSurfaces.insert({aid, asf});
643  }
644  }
645  }
646  // return the layer representation
647  return layRep;
648 }
649 
651  const Acts::ISurfaceMaterial& sMaterial) {
652  json smj;
653  // A bin utility needs to be written
654  const Acts::BinUtility* bUtility = nullptr;
655  // Check if we have a proto material
656  auto psMaterial = dynamic_cast<const Acts::ProtoSurfaceMaterial*>(&sMaterial);
657  if (psMaterial != nullptr) {
658  // Type is proto material
659  smj[m_cfg.typekey] = "proto";
660  // by default the protoMaterial is not used for mapping
661  smj[m_cfg.mapkey] = false;
662  bUtility = &(psMaterial->binUtility());
663  } else {
664  // Now check if we have a homogeneous material
665  auto hsMaterial =
666  dynamic_cast<const Acts::HomogeneousSurfaceMaterial*>(&sMaterial);
667  if (hsMaterial != nullptr) {
668  // type is homogeneous
669  smj[m_cfg.typekey] = "homogeneous";
670  smj[m_cfg.mapkey] = true;
671  if (m_cfg.writeData) {
672  smj[m_cfg.datakey] = json::array({
673  json::array({
674  encodeMaterialSlab(hsMaterial->materialSlab(0, 0)),
675  }),
676  });
677  }
678  } else {
679  // Only option remaining: BinnedSurface material
680  auto bsMaterial =
681  dynamic_cast<const Acts::BinnedSurfaceMaterial*>(&sMaterial);
682  if (bsMaterial != nullptr) {
683  // type is binned
684  smj[m_cfg.typekey] = "binned";
685  smj[m_cfg.mapkey] = true;
686  bUtility = &(bsMaterial->binUtility());
687  // convert the data
688  // get the material matrix
689  if (m_cfg.writeData) {
690  json mmat = json::array();
691  for (const auto& mpVector : bsMaterial->fullMaterial()) {
692  json mvec = json::array();
693  for (const auto& mp : mpVector) {
694  mvec.push_back(encodeMaterialSlab(mp));
695  }
696  mmat.push_back(std::move(mvec));
697  }
698  smj[m_cfg.datakey] = std::move(mmat);
699  }
700  }
701  }
702  }
703  // add the bin utility
704  if (bUtility != nullptr && !bUtility->binningData().empty()) {
705  std::vector<std::string> binkeys = {m_cfg.bin0key, m_cfg.bin1key};
706  // loop over dimensions and write
707  auto& binningData = bUtility->binningData();
708  // loop over the dimensions
709  for (size_t ibin = 0; ibin < binningData.size(); ++ibin) {
710  json binj;
711  auto cbData = binningData[ibin];
712  binj.push_back(Acts::binningValueNames[cbData.binvalue]);
713  if (cbData.option == Acts::closed) {
714  binj.push_back("closed");
715  } else {
716  binj.push_back("open");
717  }
718  binj.push_back(cbData.bins());
719  // If protoMaterial has a non uniform binning (non default) then it is
720  // used by default in the mapping
721  if (smj[m_cfg.typekey] == "proto" && cbData.bins() > 1)
722  smj[m_cfg.mapkey] = true;
723  // If it's not a proto map, write min / max
724  if (smj[m_cfg.typekey] != "proto") {
725  std::pair<double, double> minMax = {cbData.min, cbData.max};
726  binj.push_back(minMax);
727  }
728  smj[binkeys[ibin]] = binj;
729  }
730  std::vector<double> transfo;
731  Acts::Transform3D transfo_matrix = bUtility->transform();
732  if (not transfo_matrix.isApprox(Acts::Transform3D::Identity())) {
733  for (int i = 0; i < 4; i++) {
734  for (int j = 0; j < 4; j++) {
735  transfo.push_back(transfo_matrix(j, i));
736  }
737  }
738  smj[m_cfg.transfokeys] = transfo;
739  }
740  }
741  return smj;
742 }
743 
745  const Acts::IVolumeMaterial& vMaterial) {
746  json smj;
747  // A bin utility needs to be written
748  const Acts::BinUtility* bUtility = nullptr;
749  // Check if we have a proto material
750  auto pvMaterial = dynamic_cast<const Acts::ProtoVolumeMaterial*>(&vMaterial);
751  if (pvMaterial != nullptr) {
752  // Type is proto material
753  smj[m_cfg.typekey] = "proto";
754  // by default the protoMaterial is not used for mapping
755  smj[m_cfg.mapkey] = false;
756  bUtility = &(pvMaterial->binUtility());
757  } else {
758  // Now check if we have a homogeneous material
759  auto hvMaterial =
760  dynamic_cast<const Acts::HomogeneousVolumeMaterial*>(&vMaterial);
761  if (hvMaterial != nullptr) {
762  // type is homogeneous
763  smj[m_cfg.typekey] = "homogeneous";
764  smj[m_cfg.mapkey] = true;
765  if (m_cfg.writeData) {
766  // array of encoded materials w/ one entry
767  smj[m_cfg.datakey] = json::array({
768  encodeMaterial(hvMaterial->material({0, 0, 0})),
769  });
770  }
771  } else {
772  // Only option remaining: material map
773  auto bvMaterial2D = dynamic_cast<
775  &vMaterial);
776  // Now check if we have a 2D map
777  if (bvMaterial2D != nullptr) {
778  // type is binned
779  smj[m_cfg.typekey] = "interpolated2D";
780  smj[m_cfg.mapkey] = true;
781  bUtility = &(bvMaterial2D->binUtility());
782  // convert the data
783  if (m_cfg.writeData) {
784  json mmat = json::array();
785  MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
786  for (size_t bin = 0; bin < grid.size(); bin++) {
787  mmat.push_back(encodeMaterial(grid.at(bin)));
788  }
789  smj[m_cfg.datakey] = std::move(mmat);
790  }
791  } else {
792  // Only option remaining: material map
793  auto bvMaterial3D = dynamic_cast<const Acts::InterpolatedMaterialMap<
794  MaterialMapper<MaterialGrid3D>>*>(&vMaterial);
795  // Now check if we have a 3D map
796  if (bvMaterial3D != nullptr) {
797  // type is binned
798  smj[m_cfg.typekey] = "interpolated3D";
799  smj[m_cfg.mapkey] = true;
800  bUtility = &(bvMaterial3D->binUtility());
801  // convert the data
802  if (m_cfg.writeData) {
803  json mmat = json::array();
804  MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
805  for (size_t bin = 0; bin < grid.size(); bin++) {
806  mmat.push_back(encodeMaterial(grid.at(bin)));
807  }
808  smj[m_cfg.datakey] = std::move(mmat);
809  }
810  }
811  }
812  }
813  }
814  // add the bin utility
815  if (bUtility != nullptr && !bUtility->binningData().empty()) {
816  std::vector<std::string> binkeys = {m_cfg.bin0key, m_cfg.bin1key,
817  m_cfg.bin2key};
818  // loop over dimensions and write
819  auto& binningData = bUtility->binningData();
820  // loop over the dimensions
821  for (size_t ibin = 0; ibin < binningData.size(); ++ibin) {
822  json binj;
823  auto cbData = binningData[ibin];
824  binj.push_back(Acts::binningValueNames[cbData.binvalue]);
825  if (cbData.option == Acts::closed) {
826  binj.push_back("closed");
827  } else {
828  binj.push_back("open");
829  }
830  binj.push_back(cbData.bins());
831  // If protoMaterial has a non uniform binning (non default) then it is
832  // used by default in the mapping
833  if (smj[m_cfg.typekey] == "proto" && cbData.bins() > 1)
834  smj[m_cfg.mapkey] = true;
835  // If it's not a proto map, write min / max
836  if (smj[m_cfg.typekey] != "proto") {
837  std::pair<double, double> minMax = {cbData.min, cbData.max};
838  binj.push_back(minMax);
839  }
840  smj[binkeys[ibin]] = binj;
841  }
842  std::vector<double> transfo;
843  Acts::Transform3D transfo_matrix = bUtility->transform();
844  for (int i = 0; i < 4; i++) {
845  for (int j = 0; j < 4; j++) {
846  transfo.push_back(transfo_matrix(j, i));
847  }
848  }
849  smj[m_cfg.transfokeys] = transfo;
850  }
851  return smj;
852 }
853 
855  const Surface* surface) {
856  // Get the ID of the surface (redundant but help readability)
857  std::ostringstream SurfaceID;
858  SurfaceID << surface->geometryId();
859  sjson[m_cfg.surfacegeometryidkey] = SurfaceID.str();
860 
861  // Cast the surface bound to both disk and cylinder
862  const Acts::SurfaceBounds& surfaceBounds = surface->bounds();
863  auto sTransform = surface->transform(GeometryContext());
864 
865  const Acts::RadialBounds* radialBounds =
866  dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds);
867  const Acts::CylinderBounds* cylinderBounds =
868  dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds);
869  const Acts::AnnulusBounds* annulusBounds =
870  dynamic_cast<const Acts::AnnulusBounds*>(&surfaceBounds);
871 
872  if (radialBounds != nullptr) {
873  sjson[m_cfg.surfacetypekey] = "Disk";
874  sjson[m_cfg.surfacepositionkey] = sTransform.translation().z();
875  sjson[m_cfg.surfacerangekey] = {radialBounds->rMin(), radialBounds->rMax()};
876  }
877  if (cylinderBounds != nullptr) {
878  sjson[m_cfg.surfacetypekey] = "Cylinder";
879  sjson[m_cfg.surfacepositionkey] = cylinderBounds->get(CylinderBounds::eR);
880  sjson[m_cfg.surfacerangekey] = {
881  -1 * cylinderBounds->get(CylinderBounds::eHalfLengthZ),
882  cylinderBounds->get(CylinderBounds::eHalfLengthZ)};
883  }
884  if (annulusBounds != nullptr) {
885  sjson[m_cfg.surfacetypekey] = "Annulus";
886  sjson[m_cfg.surfacepositionkey] = sTransform.translation().z();
887  sjson[m_cfg.surfacerangekey] = {
888  {annulusBounds->rMin(), annulusBounds->rMax()},
889  {annulusBounds->phiMin(), annulusBounds->phiMax()}};
890  }
891 }
892 
895  const json& data) {
896  Acts::MaterialSlabMatrix mpMatrix;
897  // the input data must be array[array[object]]
898  for (auto& outer : data) {
899  Acts::MaterialSlabVector mpVector;
900  for (auto& inner : outer) {
901  mpVector.emplace_back(decodeMaterialSlab(inner));
902  }
903  mpMatrix.push_back(std::move(mpVector));
904  }
905  return mpMatrix;
906 }
907 
910  const json& bin) {
911  if (bin.size() >= 3) {
912  // finding the iterator position to determine the binning value
913  auto bit = std::find(Acts::binningValueNames.begin(),
914  Acts::binningValueNames.end(), bin[0]);
915  size_t indx = std::distance(Acts::binningValueNames.begin(), bit);
917  Acts::BinningOption bopt = bin[1] == "open" ? Acts::open : Acts::closed;
918  unsigned int bins = bin[2];
919  float min = 0;
920  float max = 0;
921  if (bin.size() >= 4 && bin[3].size() == 2) {
922  min = bin[3][0];
923  max = bin[3][1];
924  }
925  return Acts::BinUtility(bins, min, max, bopt, bval);
926  }
927  return Acts::BinUtility();
928 }
929 
932  const json& transfo) {
934  int i = 0;
935  int j = 0;
936  for (auto& element : transfo) {
937  transform(j, i) = element;
938  j++;
939  if (j == 4) {
940  i++;
941  j = 0;
942  }
943  }
944  return transform;
945 }
946 
948  const Acts::Surface& surface) {
949  Acts::BinUtility bUtility;
950 
951  const Acts::SurfaceBounds& surfaceBounds = surface.bounds();
952  const Acts::RadialBounds* radialBounds =
953  dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds);
954  const Acts::CylinderBounds* cylinderBounds =
955  dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds);
956  const Acts::AnnulusBounds* annulusBounds =
957  dynamic_cast<const Acts::AnnulusBounds*>(&surfaceBounds);
958  const Acts::RectangleBounds* rectangleBounds =
959  dynamic_cast<const Acts::RectangleBounds*>(&surfaceBounds);
960 
961  if (radialBounds != nullptr) {
962  bUtility += BinUtility(
963  1,
964  radialBounds->get(RadialBounds::eAveragePhi) -
965  radialBounds->get(RadialBounds::eHalfPhiSector),
966  radialBounds->get(RadialBounds::eAveragePhi) +
967  radialBounds->get(RadialBounds::eHalfPhiSector),
968  (radialBounds->get(RadialBounds::eHalfPhiSector) - M_PI) < s_epsilon
969  ? Acts::closed
970  : Acts::open,
971  Acts::binPhi);
972  bUtility += BinUtility(1, radialBounds->rMin(), radialBounds->rMax(),
974  return bUtility;
975  }
976  if (cylinderBounds != nullptr) {
977  bUtility += BinUtility(
978  1,
979  cylinderBounds->get(CylinderBounds::eAveragePhi) -
980  cylinderBounds->get(CylinderBounds::eHalfPhiSector),
981  cylinderBounds->get(CylinderBounds::eAveragePhi) +
982  cylinderBounds->get(CylinderBounds::eHalfPhiSector),
983  (cylinderBounds->get(CylinderBounds::eHalfPhiSector) - M_PI) < s_epsilon
984  ? Acts::closed
985  : Acts::open,
986  Acts::binPhi);
987  bUtility +=
988  BinUtility(1, -1 * cylinderBounds->get(CylinderBounds::eHalfLengthZ),
989  cylinderBounds->get(CylinderBounds::eHalfLengthZ),
991  return bUtility;
992  }
993  if (annulusBounds != nullptr) {
994  bUtility += BinUtility(1, annulusBounds->get(AnnulusBounds::eMinPhiRel),
995  annulusBounds->get(AnnulusBounds::eMaxPhiRel),
997  bUtility += BinUtility(1, annulusBounds->rMin(), annulusBounds->rMax(),
999  return bUtility;
1000  }
1001  if (rectangleBounds != nullptr) {
1002  bUtility += BinUtility(1, rectangleBounds->get(RectangleBounds::eMinX),
1003  rectangleBounds->get(RectangleBounds::eMaxX),
1005  bUtility += BinUtility(1, rectangleBounds->get(RectangleBounds::eMinY),
1006  rectangleBounds->get(RectangleBounds::eMaxY),
1008  return bUtility;
1009  }
1010  ACTS_INFO(
1011  "No corresponding bound found for the surface : " << surface.name());
1012  return bUtility;
1013 }
1014 
1016  const Acts::TrackingVolume& volume) {
1017  Acts::BinUtility bUtility;
1018 
1019  auto cyBounds =
1020  dynamic_cast<const CylinderVolumeBounds*>(&(volume.volumeBounds()));
1021  auto cutcylBounds =
1022  dynamic_cast<const CutoutCylinderVolumeBounds*>(&(volume.volumeBounds()));
1023  auto cuBounds =
1024  dynamic_cast<const CuboidVolumeBounds*>(&(volume.volumeBounds()));
1025 
1026  if (cyBounds != nullptr) {
1027  bUtility += BinUtility(1, cyBounds->get(CylinderVolumeBounds::eMinR),
1028  cyBounds->get(CylinderVolumeBounds::eMaxR),
1030  bUtility += BinUtility(
1031  1, -cyBounds->get(CylinderVolumeBounds::eHalfPhiSector),
1032  cyBounds->get(CylinderVolumeBounds::eHalfPhiSector),
1034  ? Acts::closed
1035  : Acts::open,
1036  Acts::binPhi);
1037  bUtility +=
1039  cyBounds->get(CylinderVolumeBounds::eHalfLengthZ),
1041  return bUtility;
1042  }
1043  if (cutcylBounds != nullptr) {
1044  bUtility +=
1045  BinUtility(1, cutcylBounds->get(CutoutCylinderVolumeBounds::eMinR),
1046  cutcylBounds->get(CutoutCylinderVolumeBounds::eMaxR),
1048  bUtility += BinUtility(1, -M_PI, M_PI, Acts::closed, Acts::binPhi);
1049  bUtility += BinUtility(
1050  1, -cutcylBounds->get(CutoutCylinderVolumeBounds::eHalfLengthZ),
1052  Acts::binZ);
1053  return bUtility;
1054  } else if (cuBounds != nullptr) {
1055  bUtility += BinUtility(1, -cuBounds->get(CuboidVolumeBounds::eHalfLengthX),
1056  cuBounds->get(CuboidVolumeBounds::eHalfLengthX),
1058  bUtility += BinUtility(1, -cuBounds->get(CuboidVolumeBounds::eHalfLengthY),
1059  cuBounds->get(CuboidVolumeBounds::eHalfLengthY),
1061  bUtility += BinUtility(1, -cuBounds->get(CuboidVolumeBounds::eHalfLengthZ),
1062  cuBounds->get(CuboidVolumeBounds::eHalfLengthZ),
1064  return bUtility;
1065  }
1066  ACTS_INFO(
1067  "No corresponding bound found for the volume : " << volume.volumeName());
1068  return bUtility;
1069 }