Core/src/FWGeometry.cc

0001 #include "TFile.h"
0002 #include "TTree.h"
0003 #include "TEveGeoNode.h"
0004 #include "TEveGeoShape.h"
0005 #include "TPRegexp.h"
0006 #include "TSystem.h"
0007 #include "TGeoArb8.h"
0008 #include "TObjArray.h"
0009 #include "TObjString.h"
0010 #include "TPRegexp.h"
0011
0012 #include "Fireworks/Core/interface/FWGeometry.h"
0013 #include "Fireworks/Core/interface/fwLog.h"
0014 #include "DataFormats/DetId/interface/DetId.h"
0015
0016 // AMT deprication of tracker specific DetIds
0017 #include "CalibTracker/StandaloneTrackerTopology/interface/StandaloneTrackerTopology.h"
0018
0019 #include <cassert>
0020 #include <iostream>
0021 #include <memory>
0022
0023 #include <sstream>
0024 #include <stdexcept>
0025 #include <algorithm>
0026
0027 FWGeometry::FWGeometry(void) : m_producerVersion(0) {}
0028
0029 FWGeometry::~FWGeometry(void) {}
0030
0031 bool FWGeometry::isEmpty() const {
0032   // AMT this is a check if geomtery is not loaded
0033   // e.g. cmsShow starts with no data file and without given explicit argument ( --geometry-file option )
0034
0035   return m_idToInfo.empty();
0036 }
0037
0038 TFile* FWGeometry::findFile(const char* fileName) {
0039   std::string searchPath = ".";
0040
0041   if (gSystem->Getenv("CMSSW_SEARCH_PATH")) {
0042     TString paths = gSystem->Getenv("CMSSW_SEARCH_PATH");
0043
0044     TObjArray* tokens = paths.Tokenize(":");
0045     for (int i = 0; i < tokens->GetEntries(); ++i) {
0046       TObjString* path = (TObjString*)tokens->At(i);
0047       searchPath += ":";
0048       searchPath += static_cast<const char*>(path->GetString());
0049       if (gSystem->Getenv("CMSSW_VERSION"))
0050         searchPath += "/Fireworks/Geometry/data/";
0051     }
0052   }
0053
0054   TString fn = fileName;
0055   const char* fp = gSystem->FindFile(searchPath.c_str(), fn, kFileExists);
0056   return fp ? TFile::Open(fp) : nullptr;
0057 }
0058
0059 void FWGeometry::applyGlobalTag(const std::string& globalTag) {
0060   const std::string fnRun2 = "cmsGeomRun2.root";
0061   const std::string fnRun3 = "cmsGeom2021.root";
0062   const std::string fnSLHC = "cmsGeom2026.root";
0063
0064   TPMERegexp year_re("^[^_]+_[a-zA-Z]*20(\\d\\d)_");
0065   TPMERegexp run_re("^[^_]+_[a-zA-Z]*Run(\\d)_");
0066
0067   TString test = globalTag.c_str();
0068   std::string cfn;
0069   if (year_re.Match(test)) {
0070     TString r = year_re[1];
0071     int year = atoi(r.Data());
0072     if (year < 18) {
0073       cfn = fnRun2;
0074     } else if (year < 21) {
0075       cfn = fnRun3;
0076     } else {
0077       cfn = fnSLHC;
0078     }
0079   } else if (run_re.Match(test)) {
0080     TString rn = run_re[1];
0081     if (rn == "1") {
0082       fwLog(fwlog::kWarning) << "Run1 geometry not included. Using Run2 geometry." << std::endl;
0083       cfn = fnRun2;
0084     } else if (rn == "2") {
0085       cfn = fnRun2;
0086     } else if (rn == "4") {
0087       cfn = fnSLHC;
0088     } else {
0089       fwLog(fwlog::kWarning) << "Detected Run" << rn << ". Using geometry scenario 2021.\n";
0090       cfn = fnRun3;
0091     }
0092   } else {
0093     fwLog(fwlog::kWarning) << "Could not guess geometry from global tag.  Using geometry scenario 2021.\n";
0094     cfn = fnRun3;
0095   }
0096
0097   fwLog(fwlog::kInfo) << "Guessed geometry " << cfn << " from global tag " << globalTag << std::endl;
0098   if (cfn.compare(m_fileName)) {
0099     loadMap(cfn.c_str());
0100   }
0101 }
0102
0103 void FWGeometry::loadMap(const char* iFileName) {
0104   TFile* file = findFile(iFileName);
0105   if (!file) {
0106     throw std::runtime_error("ERROR: failed to find geometry file. Initialization failed.");
0107     return;
0108   }
0109   m_fileName = iFileName;
0110   TTree* tree = static_cast<TTree*>(file->Get("idToGeo"));
0111   if (!tree) {
0112     throw std::runtime_error("ERROR: cannot find detector id map in the file. Initialization failed.");
0113     return;
0114   }
0115
0116   unsigned int id;
0117   Float_t points[24];
0118   Float_t topology[9];
0119   Float_t shape[5];
0120   Float_t translation[3];
0121   Float_t matrix[9];
0122   bool loadPoints = tree->GetBranch("points") != nullptr;
0123   bool loadParameters = tree->GetBranch("topology") != nullptr;
0124   bool loadShape = tree->GetBranch("shape") != nullptr;
0125   bool loadTranslation = tree->GetBranch("translation") != nullptr;
0126   bool loadMatrix = tree->GetBranch("matrix") != nullptr;
0127   tree->SetBranchAddress("id", &id);
0128   if (loadPoints)
0129     tree->SetBranchAddress("points", &points);
0130   if (loadParameters)
0131     tree->SetBranchAddress("topology", &topology);
0132   if (loadShape)
0133     tree->SetBranchAddress("shape", &shape);
0134   if (loadTranslation)
0135     tree->SetBranchAddress("translation", &translation);
0136   if (loadMatrix)
0137     tree->SetBranchAddress("matrix", &matrix);
0138
0139   // reset previous values
0140   m_idToInfo.clear();
0141   for (const auto& p : m_idToMatrix)
0142     delete p.second;
0143   m_trackerTopology.reset();
0144
0145   unsigned int treeSize = tree->GetEntries();
0146   if (m_idToInfo.size() != treeSize)
0147     m_idToInfo.resize(treeSize);
0148   for (unsigned int i = 0; i < treeSize; ++i) {
0149     tree->GetEntry(i);
0150
0151     m_idToInfo[i].id = id;
0152     if (loadPoints) {
0153       for (unsigned int j = 0; j < 24; ++j)
0154         m_idToInfo[i].points[j] = points[j];
0155     }
0156     if (loadParameters) {
0157       for (unsigned int j = 0; j < 9; ++j)
0158         m_idToInfo[i].parameters[j] = topology[j];
0159     }
0160     if (loadShape) {
0161       for (unsigned int j = 0; j < 5; ++j)
0162         m_idToInfo[i].shape[j] = shape[j];
0163     }
0164     if (loadTranslation) {
0165       for (unsigned int j = 0; j < 3; ++j)
0166         m_idToInfo[i].translation[j] = translation[j];
0167     }
0168     if (loadMatrix) {
0169       for (unsigned int j = 0; j < 9; ++j)
0170         m_idToInfo[i].matrix[j] = matrix[j];
0171     }
0172   }
0173
0174   m_versionInfo.productionTag = static_cast<TNamed*>(file->Get("tag"));
0175   m_versionInfo.cmsswVersion = static_cast<TNamed*>(file->Get("CMSSW_VERSION"));
0176   m_versionInfo.extraDetectors = static_cast<TObjArray*>(file->Get("ExtraDetectors"));
0177
0178   TString path = file->GetPath();
0179   if (path.EndsWith(":/"))
0180     path.Resize(path.Length() - 2);
0181
0182   if (m_versionInfo.productionTag)
0183     fwLog(fwlog::kInfo) << Form(
0184         "Load %s %s from %s\n", tree->GetName(), m_versionInfo.productionTag->GetTitle(), path.Data());
0185   else
0186     fwLog(fwlog::kInfo) << Form("Load %s from %s\n", tree->GetName(), path.Data());
0187
0188   TNamed* producerInfo = static_cast<TNamed*>(file->Get("PRODUCER_VERSION"));
0189   if (producerInfo) {
0190     m_producerVersion = atoi(producerInfo->GetTitle());
0191   }
0192
0193   TNamed* ttopology = static_cast<TNamed*>(file->Get("TrackerTopology"));
0194   if (ttopology) {
0195     std::string xml = ttopology->GetTitle();
0196     m_trackerTopology =
0197         std::make_unique<TrackerTopology>(StandaloneTrackerTopology::fromTrackerParametersXMLString(xml));
0198   }
0199
0200   file->Close();
0201 }
0202
0203 void FWGeometry::initMap(const FWRecoGeom::InfoMap& map) {
0204   FWRecoGeom::InfoMapItr begin = map.begin();
0205   FWRecoGeom::InfoMapItr end = map.end();
0206   unsigned int mapSize = map.size();
0207   if (m_idToInfo.size() != mapSize)
0208     m_idToInfo.resize(mapSize);
0209   unsigned int i = 0;
0210   for (FWRecoGeom::InfoMapItr it = begin; it != end; ++it, ++i) {
0211     m_idToInfo[i].id = it->id;
0212     for (unsigned int j = 0; j < 24; ++j)
0213       m_idToInfo[i].points[j] = it->points[j];
0214     for (unsigned int j = 0; j < 9; ++j)
0215       m_idToInfo[i].parameters[j] = it->topology[j];
0216     for (unsigned int j = 0; j < 5; ++j)
0217       m_idToInfo[i].shape[j] = it->shape[j];
0218     for (unsigned int j = 0; j < 3; ++j)
0219       m_idToInfo[i].translation[j] = it->translation[j];
0220     for (unsigned int j = 0; j < 9; ++j)
0221       m_idToInfo[i].matrix[j] = it->matrix[j];
0222   }
0223 }
0224
0225 const TGeoMatrix* FWGeometry::getMatrix(unsigned int id) const {
0226   std::map<unsigned int, TGeoMatrix*>::iterator mit = m_idToMatrix.find(id);
0227   if (mit != m_idToMatrix.end())
0228     return mit->second;
0229
0230   IdToInfoItr it = FWGeometry::find(id);
0231   if (it == m_idToInfo.end()) {
0232     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0233     return nullptr;
0234   } else {
0235     const GeomDetInfo& info = *it;
0236     TGeoTranslation trans(info.translation[0], info.translation[1], info.translation[2]);
0237     TGeoRotation rotation;
0238     const Double_t matrix[9] = {info.matrix[0],
0239                                 info.matrix[1],
0240                                 info.matrix[2],
0241                                 info.matrix[3],
0242                                 info.matrix[4],
0243                                 info.matrix[5],
0244                                 info.matrix[6],
0245                                 info.matrix[7],
0246                                 info.matrix[8]};
0247     rotation.SetMatrix(matrix);
0248
0249     m_idToMatrix[id] = new TGeoCombiTrans(trans, rotation);
0250     return m_idToMatrix[id];
0251   }
0252 }
0253
0254 std::vector<unsigned int> FWGeometry::getMatchedIds(Detector det, SubDetector subdet) const {
0255   std::vector<unsigned int> ids;
0256   unsigned int mask = (det << 4) | (subdet);
0257   for (IdToInfoItr it = m_idToInfo.begin(), itEnd = m_idToInfo.end(); it != itEnd; ++it) {
0258     if (FWGeometry::match_id(*it, mask))
0259       ids.push_back((*it).id);
0260   }
0261
0262   return ids;
0263 }
0264
0265 std::vector<unsigned int> FWGeometry::getMatchedIds(Detector det) const {
0266   std::vector<unsigned int> ids;
0267   for (const auto& it : m_idToInfo) {
0268     if (((it.id >> kDetOffset) & 0xF) != det)
0269       continue;
0270
0271     // select only the fake DetIds that have all the (u,v) bits set at 1.  This
0272     // is used to draw the HGCal Geometry that is wafer-based for the silicon
0273     // part. The Scintillators are treated on a tile-basis.
0274     if (det == HGCalHSc) {
0275       ids.push_back(it.id);
0276     } else {
0277       auto key = 0x3FF;  // 10 bits mask of 1s.
0278       if ((it.id | key) == it.id) {
0279         ids.push_back(it.id);
0280       }
0281     }
0282   }
0283   return ids;
0284 }
0285
0286 TGeoShape* FWGeometry::getShape(unsigned int id) const {
0287   IdToInfoItr it = FWGeometry::find(id);
0288   if (it == m_idToInfo.end()) {
0289     fwLog(fwlog::kWarning) << "no reco geoemtry found for id " << id << std::endl;
0290     return nullptr;
0291   } else {
0292     return getShape(*it);
0293   }
0294 }
0295
0296 TGeoShape* FWGeometry::getShape(const GeomDetInfo& info) const {
0297   TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
0298   TGeoShape* geoShape = nullptr;
0299   if (info.shape[0] == 1) {
0300     geoShape = new TGeoTrap(info.shape[3],  //dz
0301                             0,              //theta
0302                             0,              //phi
0303                             info.shape[4],  //dy1
0304                             info.shape[1],  //dx1
0305                             info.shape[2],  //dx2
0306                             0,              //alpha1
0307                             info.shape[4],  //dy2
0308                             info.shape[1],  //dx3
0309                             info.shape[2],  //dx4
0310                             0);             //alpha2
0311   } else
0312     geoShape = new TGeoBBox(info.shape[1], info.shape[2], info.shape[3]);
0313
0314   return geoShape;
0315 }
0316
0317 TEveGeoShape* FWGeometry::getEveShape(unsigned int id) const {
0318   IdToInfoItr it = FWGeometry::find(id);
0319   if (it == m_idToInfo.end()) {
0320     fwLog(fwlog::kWarning) << "no reco geoemtry found for id " << id << std::endl;
0321     return nullptr;
0322   } else {
0323     const GeomDetInfo& info = *it;
0324     double array[16] = {info.matrix[0],
0325                         info.matrix[3],
0326                         info.matrix[6],
0327                         0.,
0328                         info.matrix[1],
0329                         info.matrix[4],
0330                         info.matrix[7],
0331                         0.,
0332                         info.matrix[2],
0333                         info.matrix[5],
0334                         info.matrix[8],
0335                         0.,
0336                         info.translation[0],
0337                         info.translation[1],
0338                         info.translation[2],
0339                         1.};
0340     TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
0341     TEveGeoShape* shape = new TEveGeoShape(TString::Format("RecoGeom Id=%u", id));
0342     TGeoShape* geoShape = getShape(info);
0343     shape->SetShape(geoShape);
0344     // Set transformation matrix from a column-major array
0345     shape->SetTransMatrix(array);
0346     return shape;
0347   }
0348 }
0349
0350 TEveGeoShape* FWGeometry::getHGCSiliconEveShape(unsigned int id) const {
0351   IdToInfoItr it = FWGeometry::find(id);
0352   if (it == m_idToInfo.end()) {
0353     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0354     return nullptr;
0355   }
0356
0357   GeomDetInfo info = *it;
0358
0359   TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
0360   TEveGeoShape* shape = new TEveGeoShape(TString::Format("RecoGeom Id=%u", id));
0361
0362   TGeoXtru* geoShape = new TGeoXtru(2);
0363   Double_t x[6];
0364   Double_t y[6];
0365   for (unsigned int i = 0; i < 6; ++i) {
0366     x[i] = info.points[i * 3];
0367     y[i] = info.points[3 * i + 1];
0368   }
0369   geoShape->DefinePolygon(6, x, y);
0370   geoShape->DefineSection(0, info.points[2] - 0.0150);  // First plane at the Z position of the wafer, minus 150um
0371   geoShape->DefineSection(1, info.points[2] + 0.0150);  // Second plane at the Z position of the wafer, minus 150um
0372
0373   shape->SetShape(geoShape);
0374   double array[16] = {info.matrix[0],
0375                       info.matrix[3],
0376                       info.matrix[6],
0377                       0.,
0378                       info.matrix[1],
0379                       info.matrix[4],
0380                       info.matrix[7],
0381                       0.,
0382                       info.matrix[2],
0383                       info.matrix[5],
0384                       info.matrix[8],
0385                       0.,
0386                       0.,  // translation x
0387                       0.,  // translation y
0388                       0.,  // translation z
0389                       1.};
0390   // Set transformation matrix from a column-major array
0391   shape->SetTransMatrix(array);
0392   return shape;
0393 }
0394
0395 TEveGeoShape* FWGeometry::getHGCScintillatorEveShape(unsigned int id) const {
0396   IdToInfoItr it = FWGeometry::find(id);
0397   if (it == m_idToInfo.end()) {
0398     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0399     return nullptr;
0400   }
0401
0402   GeomDetInfo info = *it;
0403
0404   TEveGeoManagerHolder gmgr(TEveGeoShape::GetGeoMangeur());
0405   TEveGeoShape* shape = new TEveGeoShape(TString::Format("RecoGeom Id=%u", id));
0406
0407   TGeoXtru* geoShape = new TGeoXtru(2);
0408   Double_t x[4] = {info.points[0], info.points[3], info.points[6], info.points[9]};
0409   Double_t y[4] = {info.points[1], info.points[4], info.points[7], info.points[10]};
0410
0411   bool isNeg = info.shape[3] < 0;
0412   geoShape->DefinePolygon(4, x, y);
0413   geoShape->DefineSection(0, isNeg * info.shape[3]);
0414   geoShape->DefineSection(1, !isNeg * info.shape[3]);
0415   info.translation[2] = info.points[2];
0416
0417   shape->SetShape(geoShape);
0418   double array[16] = {info.matrix[0],
0419                       info.matrix[3],
0420                       info.matrix[6],
0421                       0.,
0422                       info.matrix[1],
0423                       info.matrix[4],
0424                       info.matrix[7],
0425                       0.,
0426                       info.matrix[2],
0427                       info.matrix[5],
0428                       info.matrix[8],
0429                       0.,
0430                       info.translation[0],
0431                       info.translation[1],
0432                       info.translation[2],
0433                       1.};
0434   // Set transformation matrix from a column-major array
0435   shape->SetTransMatrix(array);
0436   return shape;
0437 }
0438
0439 const float* FWGeometry::getCorners(unsigned int id) const {
0440   // reco geometry points
0441   IdToInfoItr it = FWGeometry::find(id);
0442   if (it == m_idToInfo.end()) {
0443     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0444     return nullptr;
0445   } else {
0446     return (*it).points;
0447   }
0448 }
0449
0450 const float* FWGeometry::getParameters(unsigned int id) const {
0451   // reco geometry parameters
0452   IdToInfoItr it = FWGeometry::find(id);
0453   if (it == m_idToInfo.end()) {
0454     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0455     return nullptr;
0456   } else {
0457     return (*it).parameters;
0458   }
0459 }
0460
0461 const float* FWGeometry::getShapePars(unsigned int id) const {
0462   // reco geometry parameters
0463   IdToInfoItr it = FWGeometry::find(id);
0464   if (it == m_idToInfo.end()) {
0465     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0466     return nullptr;
0467   } else {
0468     return (*it).shape;
0469   }
0470 }
0471
0472 void FWGeometry::localToGlobal(unsigned int id, const float* local, float* global, bool translatep) const {
0473   IdToInfoItr it = FWGeometry::find(id);
0474   if (it == m_idToInfo.end()) {
0475     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0476   } else {
0477     localToGlobal(*it, local, global, translatep);
0478   }
0479 }
0480
0481 void FWGeometry::localToGlobal(
0482     unsigned int id, const float* local1, float* global1, const float* local2, float* global2, bool translatep) const {
0483   IdToInfoItr it = FWGeometry::find(id);
0484   if (it == m_idToInfo.end()) {
0485     fwLog(fwlog::kWarning) << "no reco geometry found for id " << id << std::endl;
0486   } else {
0487     localToGlobal(*it, local1, global1, translatep);
0488     localToGlobal(*it, local2, global2, translatep);
0489   }
0490 }
0491
0492 FWGeometry::IdToInfoItr FWGeometry::find(unsigned int id) const {
0493   FWGeometry::IdToInfoItr begin = m_idToInfo.begin();
0494   FWGeometry::IdToInfoItr end = m_idToInfo.end();
0495   return std::lower_bound(begin, end, id);
0496 }
0497
0498 void FWGeometry::localToGlobal(const GeomDetInfo& info, const float* local, float* global, bool translatep) const {
0499   for (int i = 0; i < 3; ++i) {
0500     global[i] = translatep ? info.translation[i] : 0;
0501     global[i] += local[0] * info.matrix[3 * i] + local[1] * info.matrix[3 * i + 1] + local[2] * info.matrix[3 * i + 2];
0502   }
0503 }
0504
0505 //______________________________________________________________________________
0506
0507 bool FWGeometry::VersionInfo::haveExtraDet(const char* det) const {
0508   return (extraDetectors && extraDetectors->FindObject(det)) ? true : false;
0509 }