FWGeometry.cc

Go to the documentation of this file.

#include "TFile.h"
#include "TTree.h"
#include "TEveGeoNode.h"
#include "TPRegexp.h"
#include "TSystem.h"
#include "TGeoArb8.h"
#include "TObjArray.h"

#include "Fireworks/Core/interface/FWGeometry.h"
#include "Fireworks/Core/interface/fwLog.h"
#include "DataFormats/DetId/interface/DetId.h"

// AMT deprication of tracker specific DetIds
#include "CalibTracker/StandaloneTrackerTopology/interface/StandaloneTrackerTopology.h"

#include <iostream>
#include <cassert>
#include <sstream>
#include <stdexcept>
#include <algorithm>

FWGeometry::FWGeometry(void) : m_producerVersion(0) {}

FWGeometry::~FWGeometry(void) {}

TFile* FWGeometry::findFile(const char* fileName) {
  std::string searchPath = ".";

  if (gSystem->Getenv("CMSSW_SEARCH_PATH")) {
    TString paths = gSystem->Getenv("CMSSW_SEARCH_PATH");

    TObjArray* tokens = paths.Tokenize(":");
    for (int i = 0; i < tokens->GetEntries(); ++i) {
      TObjString* path = (TObjString*)tokens->At(i);
      searchPath += ":";
      searchPath += static_cast<const char*>(path->GetString());
      if (gSystem->Getenv("CMSSW_VERSION"))
        searchPath += "/Fireworks/Geometry/data/";
    }
  }

  TString fn = fileName;
  const char* fp = gSystem->FindFile(searchPath.c_str(), fn, kFileExists);
  return fp ? TFile::Open(fp) : nullptr;
}

void FWGeometry::loadMap(const char* fileName) {
  TFile* file = findFile(fileName);
  if (!file) {
    throw std::runtime_error("ERROR: failed to find geometry file. Initialization failed.");
    return;
  }

  TTree* tree = static_cast<TTree*>(file->Get("idToGeo"));
  if (!tree) {
    throw std::runtime_error("ERROR: cannot find detector id map in the file. Initialization failed.");
    return;
  }

  unsigned int id;
  Float_t points[24];
  Float_t topology[9];
  Float_t shape[5];
  Float_t translation[3];
  Float_t matrix[9];
  bool loadPoints = tree->GetBranch("points") != nullptr;
  bool loadParameters = tree->GetBranch("topology") != nullptr;
  bool loadShape = tree->GetBranch("shape") != nullptr;
  bool loadTranslation = tree->GetBranch("translation") != nullptr;
  bool loadMatrix = tree->GetBranch("matrix") != nullptr;
  tree->SetBranchAddress("id", &id);
  if (loadPoints)
    tree->SetBranchAddress("points", &points);
  if (loadParameters)
    tree->SetBranchAddress("topology", &topology);
  if (loadShape)
    tree->SetBranchAddress("shape", &shape);
  if (loadTranslation)
    tree->SetBranchAddress("translation", &translation);
  if (loadMatrix)
    tree->SetBranchAddress("matrix", &matrix);

  unsigned int treeSize = tree->GetEntries();
  if (m_idToInfo.size() != treeSize)
    m_idToInfo.resize(treeSize);
  for (unsigned int i = 0; i < treeSize; ++i) {
    tree->GetEntry(i);

    m_idToInfo[i].id = id;
    if (loadPoints) {
      for (unsigned int j = 0; j < 24; ++j)
        m_idToInfo[i].points[j] = points[j];
    }
    if (loadParameters) {
      for (unsigned int j = 0; j < 9; ++j)
        m_idToInfo[i].parameters[j] = topology[j];
    }
    if (loadShape) {
      for (unsigned int j = 0; j < 5; ++j)
        m_idToInfo[i].shape[j] = shape[j];
    }
    if (loadTranslation) {
      for (unsigned int j = 0; j < 3; ++j)
        m_idToInfo[i].translation[j] = translation[j];
    }
    if (loadMatrix) {
      for (unsigned int j = 0; j < 9; ++j)
        m_idToInfo[i].matrix[j] = matrix[j];
    }
  }

  m_versionInfo.productionTag = static_cast<TNamed*>(file->Get("tag"));
  m_versionInfo.cmsswVersion = static_cast<TNamed*>(file->Get("CMSSW_VERSION"));
  m_versionInfo.extraDetectors = static_cast<TObjArray*>(file->Get("ExtraDetectors"));

  TString path = file->GetPath();
  if (path.EndsWith(":/"))
    path.Resize(path.Length() - 2);

  if (m_versionInfo.productionTag)
    fwLog(fwlog::kInfo) << Form(
        "Load %s %s from %s\n", tree->GetName(), m_versionInfo.productionTag->GetTitle(), path.Data());
  else
    fwLog(fwlog::kInfo) << Form("Load %s from %s\n", tree->GetName(), path.Data());

  TNamed* producerInfo = static_cast<TNamed*>(file->Get("PRODUCER_VERSION"));
  if (producerInfo) {
    m_producerVersion = atoi(producerInfo->GetTitle());
  }

  TNamed* ttopology = static_cast<TNamed*>(file->Get("TrackerTopology"));
  if (ttopology) {
    std::string xml = ttopology->GetTitle();
    m_trackerTopology = std::unique_ptr<TrackerTopology>(
        new TrackerTopology(StandaloneTrackerTopology::fromTrackerParametersXMLString(xml)));
  }

  file->Close();
}

void FWGeometry::initMap(const FWRecoGeom::InfoMap& map) {
  FWRecoGeom::InfoMapItr begin = map.begin();
  FWRecoGeom::InfoMapItr end = map.end();
  unsigned int mapSize = map.size();
  if (m_idToInfo.size() != mapSize)
    m_idToInfo.resize(mapSize);
  unsigned int i = 0;
  for (FWRecoGeom::InfoMapItr it = begin; it != end; ++it, ++i) {
    m_idToInfo[i].id = it->id;
    for (unsigned int j = 0; j < 24; ++j)
      m_idToInfo[i].points[j] = it->points[j];
    for (unsigned int j = 0; j < 9; ++j)
      m_idToInfo[i].parameters[j] = it->topology[j];
    for (unsigned int j = 0; j < 5; ++j)
      m_idToInfo[i].shape[j] = it->shape[j];
    for (unsigned int j = 0; j < 3; ++j)
      m_idToInfo[i].translation[j] = it->translation[j];
    for (unsigned int j = 0; j < 9; ++j)
      m_idToInfo[i].matrix[j] = it->matrix[j];
  }
}

const TGeoMatrix* FWGeometry::getMatrix(unsigned int id) const {
  std::map<unsigned int, TGeoMatrix*>::iterator mit = m_idToMatrix.find(id);
  if (mit != m_idToMatrix.end())
    return mit->second;

  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
    return nullptr;
  } else {
    const GeomDetInfo& info = *it;
    TGeoTranslation trans(info.translation[0], info.translation[1], info.translation[2]);
    TGeoRotation rotation;
    const Double_t matrix[9] = {info.matrix[0],
                                info.matrix[1],
                                info.matrix[2],
                                info.matrix[3],
                                info.matrix[4],
                                info.matrix[5],
                                info.matrix[6],
                                info.matrix[7],
                                info.matrix[8]};
    rotation.SetMatrix(matrix);

    m_idToMatrix[id] = new TGeoCombiTrans(trans, rotation);
    return m_idToMatrix[id];
  }
}

std::vector<unsigned int> FWGeometry::getMatchedIds(Detector det, SubDetector subdet) const {
  std::vector<unsigned int> ids;
  unsigned int mask = (det << 4) | (subdet);
  for (IdToInfoItr it = m_idToInfo.begin(), itEnd = m_idToInfo.end(); it != itEnd; ++it) {
    if (FWGeometry::match_id(*it, mask))
      ids.push_back((*it).id);
  }

  return ids;
}

std::vector<unsigned int> FWGeometry::getMatchedIds(Detector det) const {
  std::vector<unsigned int> ids;

  for (const auto& it : m_idToInfo) {
    if (((it.id >> kDetOffset) & 0xF) != det)
      continue;
    // select only the 1st layer
    // if(((it->id>>17)&0x1F) != 12) continue;

    // select only the first cell of each wafer
    if (det != HGCalHSc) {
      bool flag(false);

      for (const auto& id_it : ids) {
        flag = (~0x3FF & it.id) == (~0x3FF & id_it);
        if (flag)
          break;
      }

      if (flag)
        continue;
    }

    ids.push_back(it.id);
  }
  return ids;
}

TGeoShape* FWGeometry::getShape(unsigned int id) const {
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geoemtry found for id " << id << std::endl;
    return nullptr;
  } else {
    return getShape(*it);
  }
}

TGeoShape* FWGeometry::getShape(const GeomDetInfo& info) const {
  TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
  TGeoShape* geoShape = nullptr;
  if (info.shape[0] == 1) {
    geoShape = new TGeoTrap(info.shape[3],  //dz
                            0,              //theta
                            0,              //phi
                            info.shape[4],  //dy1
                            info.shape[1],  //dx1
                            info.shape[2],  //dx2
                            0,              //alpha1
                            info.shape[4],  //dy2
                            info.shape[1],  //dx3
                            info.shape[2],  //dx4
                            0);             //alpha2
  } else
    geoShape = new TGeoBBox(info.shape[1], info.shape[2], info.shape[3]);

  return geoShape;
}

TEveGeoShape* FWGeometry::getEveShape(unsigned int id) const {
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geoemtry found for id " << id << std::endl;
    return nullptr;
  } else {
    const GeomDetInfo& info = *it;
    double array[16] = {info.matrix[0],
                        info.matrix[3],
                        info.matrix[6],
                        0.,
                        info.matrix[1],
                        info.matrix[4],
                        info.matrix[7],
                        0.,
                        info.matrix[2],
                        info.matrix[5],
                        info.matrix[8],
                        0.,
                        info.translation[0],
                        info.translation[1],
                        info.translation[2],
                        1.};
    TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
    TEveGeoShape* shape = new TEveGeoShape(TString::Format("RecoGeom Id=%u", id));
    TGeoShape* geoShape = getShape(info);
    shape->SetShape(geoShape);
    // Set transformation matrix from a column-major array
    shape->SetTransMatrix(array);
    return shape;
  }
}

TEveGeoShape* FWGeometry::getHGCSiliconEveShape(unsigned int id) const {
#if 0 
   const unsigned int type = (id>>26)&0x3;
   // select the middle cell of each waifer
   id &= ~0x3FF;
   id |= (type == 0) ? 0x16B : 0xE7;
#else
  float sideToSideWaferSize = 16.7441f;
  float dx = sideToSideWaferSize / 2;
  float sidey = dx / sqrt(3);
  float dy = 2 * sidey;

  int waferUint = (id >> 10) & 0xF;
  int waferVint = (id >> 15) & 0xF;
  float waferU = ((id >> 14) & 0x1) ? -sideToSideWaferSize * waferUint : sideToSideWaferSize * waferUint;
  float waferV = ((id >> 19) & 0x1) ? -sideToSideWaferSize * waferVint : sideToSideWaferSize * waferVint;

  float waferX = (-2 * waferU + waferV) / 2;
  float waferY = waferV * sqrt(3) / 2;
#endif
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
    return nullptr;
  }

  GeomDetInfo info = *it;

  TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
  TEveGeoShape* shape = new TEveGeoShape(TString::Format("RecoGeom Id=%u", id));

  float dz = fabs(info.points[14] - info.points[2]) * 0.5;

  info.translation[2] = (info.points[14] + info.points[2]) / 2.0f;
  info.translation[0] = waferX * ((0 < info.translation[2]) - (info.translation[2] < 0));
  info.translation[1] = waferY;

  TGeoXtru* geoShape = new TGeoXtru(2);
  Double_t x[6] = {-dx, -dx, 0.0, dx, dx, 0.0};
  Double_t y[6] = {-sidey, sidey, dy, sidey, -sidey, -dy};
  geoShape->DefinePolygon(6, x, y);
  geoShape->DefineSection(0, -dz);
  geoShape->DefineSection(1, dz);

  shape->SetShape(geoShape);
  double array[16] = {info.matrix[0],
                      info.matrix[3],
                      info.matrix[6],
                      0.,
                      info.matrix[1],
                      info.matrix[4],
                      info.matrix[7],
                      0.,
                      info.matrix[2],
                      info.matrix[5],
                      info.matrix[8],
                      0.,
                      info.translation[0],
                      info.translation[1],
                      info.translation[2],
                      1.};
  // Set transformation matrix from a column-major array
  shape->SetTransMatrix(array);
  return shape;
}

TEveGeoShape* FWGeometry::getHGCScintillatorEveShape(unsigned int id) const {
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
    return nullptr;
  }

  GeomDetInfo info = *it;

  TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
  TEveGeoShape* shape = new TEveGeoShape(TString::Format("RecoGeom Id=%u", id));

  TGeoXtru* geoShape = new TGeoXtru(2);
  Double_t x[4] = {info.points[0], info.points[3], info.points[6], info.points[9]};
  Double_t y[4] = {info.points[1], info.points[4], info.points[7], info.points[10]};

  bool isNeg = info.shape[3] < 0;
  geoShape->DefinePolygon(4, x, y);
  geoShape->DefineSection(0, isNeg * info.shape[3]);
  geoShape->DefineSection(1, !isNeg * info.shape[3]);
  info.translation[2] = info.points[2];

  shape->SetShape(geoShape);
  double array[16] = {info.matrix[0],
                      info.matrix[3],
                      info.matrix[6],
                      0.,
                      info.matrix[1],
                      info.matrix[4],
                      info.matrix[7],
                      0.,
                      info.matrix[2],
                      info.matrix[5],
                      info.matrix[8],
                      0.,
                      info.translation[0],
                      info.translation[1],
                      info.translation[2],
                      1.};
  // Set transformation matrix from a column-major array
  shape->SetTransMatrix(array);
  return shape;
}

const float* FWGeometry::getCorners(unsigned int id) const {
  // reco geometry points
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
    return nullptr;
  } else {
    return (*it).points;
  }
}

const float* FWGeometry::getParameters(unsigned int id) const {
  // reco geometry parameters
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
    return nullptr;
  } else {
    return (*it).parameters;
  }
}

const float* FWGeometry::getShapePars(unsigned int id) const {
  // reco geometry parameters
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
    return nullptr;
  } else {
    return (*it).shape;
  }
}

void FWGeometry::localToGlobal(unsigned int id, const float* local, float* global, bool translatep) const {
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
  } else {
    localToGlobal(*it, local, global, translatep);
  }
}

void FWGeometry::localToGlobal(
    unsigned int id, const float* local1, float* global1, const float* local2, float* global2, bool translatep) const {
  IdToInfoItr it = FWGeometry::find(id);
  if (it == m_idToInfo.end()) {
    fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
  } else {
    localToGlobal(*it, local1, global1, translatep);
    localToGlobal(*it, local2, global2, translatep);
  }
}

FWGeometry::IdToInfoItr FWGeometry::find(unsigned int id) const {
  FWGeometry::IdToInfoItr begin = m_idToInfo.begin();
  FWGeometry::IdToInfoItr end = m_idToInfo.end();
  return std::lower_bound(begin, end, id);
}

void FWGeometry::localToGlobal(const GeomDetInfo& info, const float* local, float* global, bool translatep) const {
  for (int i = 0; i < 3; ++i) {
    global[i] = translatep ? info.translation[i] : 0;
    global[i] += local[0] * info.matrix[3 * i] + local[1] * info.matrix[3 * i + 1] + local[2] * info.matrix[3 * i + 2];
  }
}

//______________________________________________________________________________

bool FWGeometry::VersionInfo::haveExtraDet(const char* det) const {
  return (extraDetectors && extraDetectors->FindObject(det)) ? true : false;
}