Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-08 03:34:57 

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: HGCalNumberingScheme.cc
0003 // Description: Numbering scheme for High Granularity Calorimeter
0004 ///////////////////////////////////////////////////////////////////////////////
0005 #include "SimG4CMS/Calo/interface/HGCalNumberingScheme.h"
0006 #include "SimG4CMS/Calo/interface/CaloSimUtils.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/FileInPath.h"
0009 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0010 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0011 #include "Geometry/HGCalCommonData/interface/HGCalCassette.h"
0012 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0013 #include "Geometry/HGCalCommonData/interface/HGCalTileIndex.h"
0014 #include "Geometry/HGCalCommonData/interface/HGCalWaferIndex.h"
0015 #include <array>
0016 #include <fstream>
0017 #include <iostream>
0018 
0019 //#define EDM_ML_DEBUG
0020 
0021 HGCalNumberingScheme::HGCalNumberingScheme(const HGCalDDDConstants& hgc,
0022                                            const DetId::Detector& det,
0023                                            const std::string& name,
0024                                            const std::string& fileName)
0025     : hgcons_(hgc), mode_(hgc.geomMode()), det_(det), name_(name) {
0026 #ifdef EDM_ML_DEBUG
0027   edm::LogVerbatim("HGCSim") << "Creating HGCalNumberingScheme for " << name_ << " Det " << det_ << " Mode " << mode_
0028                              << ":" << HGCalGeometryMode::Hexagon8Full << ":" << HGCalGeometryMode::Hexagon8 << ":"
0029                              << HGCalGeometryMode::Hexagon8File << ":" << HGCalGeometryMode::Hexagon8Module << ":"
0030                              << ":" << HGCalGeometryMode::Hexagon8Cassette << ":" << HGCalGeometryMode::Trapezoid << ":"
0031                              << HGCalGeometryMode::TrapezoidFile << ":" << HGCalGeometryMode::TrapezoidModule << ":"
0032                              << HGCalGeometryMode::TrapezoidCassette << ":" << HGCalGeometryMode::TrapezoidFineCell;
0033 #endif
0034   firstLayer_ = hgcons_.getLayerOffset();
0035   if (!fileName.empty()) {
0036     edm::FileInPath filetmp1("SimG4CMS/Calo/data/" + fileName);
0037     const std::string& filetmp2 = filetmp1.fullPath();
0038     std::ifstream fInput(filetmp2.c_str());
0039     if (!fInput.good()) {
0040       edm::LogVerbatim("HGCalSim") << "Cannot open file " << filetmp2;
0041     } else {
0042       char buffer[80];
0043       while (fInput.getline(buffer, 80)) {
0044         std::vector<std::string> items = CaloSimUtils::splitString(std::string(buffer));
0045         if (items.size() == 3) {
0046           if (hgcons_.waferHexagon8File()) {
0047             int layer = std::atoi(items[0].c_str());
0048             int waferU = std::atoi(items[1].c_str());
0049             int waferV = std::atoi(items[2].c_str());
0050             indices_.emplace_back(HGCalWaferIndex::waferIndex(layer, waferU, waferV, false));
0051           } else if (hgcons_.tileTrapezoid()) {
0052             int layer = std::atoi(items[0].c_str());
0053             int ring = std::atoi(items[1].c_str());
0054             int iphi = std::atoi(items[2].c_str());
0055             indices_.emplace_back(HGCalTileIndex::tileIndex(layer, ring, iphi));
0056           }
0057         } else if (items.size() == 1) {
0058           int dumpdet = std::atoi(items[0].c_str());
0059           dumpDets_.emplace_back(dumpdet);
0060         } else if (items.size() == 4) {
0061           int idet = std::atoi(items[0].c_str());
0062           int layer = std::atoi(items[1].c_str());
0063           int zside = std::atoi(items[2].c_str());
0064           int cassette = std::atoi(items[3].c_str());
0065           dumpCassette_.emplace_back(HGCalCassette::cassetteIndex(idet, layer, zside, cassette));
0066         }
0067       }
0068 #ifdef EDM_ML_DEBUG
0069       edm::LogVerbatim("HGCalSim") << "Reads in " << indices_.size() << ":" << dumpDets_.size() << ":"
0070                                    << dumpCassette_.size() << " component information from " << filetmp2
0071                                    << " Layer Offset " << firstLayer_;
0072 #endif
0073       fInput.close();
0074     }
0075   }
0076 }
0077 
0078 HGCalNumberingScheme::~HGCalNumberingScheme() {
0079 #ifdef EDM_ML_DEBUG
0080   edm::LogVerbatim("HGCSim") << "Deleting HGCalNumberingScheme";
0081 #endif
0082 }
0083 
0084 uint32_t HGCalNumberingScheme::getUnitID(int layer, int module, int cell, int iz, const G4ThreeVector& pos, double& wt) {
0085   // module is the copy number of the wafer as placed in the layer
0086   uint32_t index(0);
0087   wt = 1.0;
0088 #ifdef EDM_ML_DEBUG
0089   edm::LogVerbatim("HGCSim") << "HGCalNumberingScheme:: input Layer " << layer << " Module " << module << " Cell "
0090                              << cell << " iz " << iz << " Position " << pos;
0091 #endif
0092   if (hgcons_.waferHexagon8()) {
0093     int cellU(0), cellV(0), waferType(-1), waferU(0), waferV(0);
0094     bool debug(false);
0095     if (cell >= 0) {
0096       waferType = HGCalTypes::getUnpackedType(module);
0097       waferU = HGCalTypes::getUnpackedU(module);
0098       waferV = HGCalTypes::getUnpackedV(module);
0099       cellU = HGCalTypes::getUnpackedCellU(cell);
0100       cellV = HGCalTypes::getUnpackedCellV(cell);
0101     } else if (mode_ != HGCalGeometryMode::Hexagon8) {
0102       int zside = (pos.z() > 0) ? 1 : -1;
0103       double xx = zside * pos.x();
0104       int wU = HGCalTypes::getUnpackedU(module);
0105       int wV = HGCalTypes::getUnpackedV(module);
0106       if (!indices_.empty()) {
0107         int indx = HGCalWaferIndex::waferIndex(firstLayer_ + layer, wU, wV, false);
0108         if (std::find(indices_.begin(), indices_.end(), indx) != indices_.end())
0109           debug = true;
0110       }
0111       if (!dumpDets_.empty()) {
0112         if ((std::find(dumpDets_.begin(), dumpDets_.end(), det_) != dumpDets_.end()) &&
0113             (hgcons_.waferInfo(layer, wU, wV).part != HGCalTypes::WaferFull))
0114           debug = true;
0115       }
0116       if (!dumpCassette_.empty()) {
0117         int indx = HGCalWaferIndex::waferIndex(firstLayer_ + layer, wU, wV, false);
0118         auto ktr = hgcons_.getParameter()->waferInfoMap_.find(indx);
0119         if (ktr != hgcons_.getParameter()->waferInfoMap_.end()) {
0120           if (std::find(dumpCassette_.begin(),
0121                         dumpCassette_.end(),
0122                         HGCalCassette::cassetteIndex(det_, firstLayer_ + layer, zside, (ktr->second).cassette)) !=
0123               dumpCassette_.end())
0124             debug = true;
0125         }
0126       }
0127       hgcons_.waferFromPosition(xx, pos.y(), zside, layer, waferU, waferV, cellU, cellV, waferType, wt, false, debug);
0128     }
0129     if (waferType >= 0) {
0130       if (hgcons_.waferHexagon8File()) {
0131         int type = hgcons_.waferType(layer, waferU, waferV, true);
0132         if (type != waferType) {
0133 #ifdef EDM_ML_DEBUG
0134           edm::LogVerbatim("HGCSim") << "HGCalNumberingScheme:: " << name_ << " Layer|u|v|Index|module|cell " << layer
0135                                      << ":" << waferU << ":" << waferV << ":"
0136                                      << HGCalWaferIndex::waferIndex(layer, waferU, waferV, false) << ":" << module
0137                                      << ":" << cell << " has a type mismatch " << waferType << ":" << type;
0138 #endif
0139           if (type != HGCSiliconDetId::HGCalLD300)
0140             waferType = type;
0141         }
0142       }
0143       index = HGCSiliconDetId(det_, iz, waferType, layer, waferU, waferV, cellU, cellV).rawId();
0144       if (debug) {
0145         int zside = (pos.z() > 0) ? 1 : -1;
0146         double xx = zside * pos.x();
0147         edm::LogVerbatim("HGCSim") << "OK Input " << det_ << ":" << zside << ":" << layer << ":" << xx << ":" << pos.y()
0148                                    << ":" << module << ":" << cell << " WaferType " << waferType << " Wafer " << waferU
0149                                    << ":" << waferV << " Cell " << cellU << ":" << cellV << " wt " << wt << " Mode "
0150                                    << mode_;
0151       }
0152 #ifdef EDM_ML_DEBUG
0153     } else {
0154       edm::LogVerbatim("HGCSim") << "Bad WaferType " << waferType << " for Layer:u:v " << layer << ":" << waferU << ":"
0155                                  << waferV;
0156 #endif
0157     }
0158   } else if (hgcons_.tileTrapezoid()) {
0159     std::array<int, 3> id = hgcons_.assignCellTrap(pos.x(), pos.y(), pos.z(), layer, false);
0160 #ifdef EDM_ML_DEBUG
0161     edm::LogVerbatim("HGCSim") << "Trapezoid Position Layer " << layer << " Position " << pos.x() << ":" << pos.y()
0162                                << ":" << pos.z() << " ID " << id[0] << ":" << id[1] << ":" << id[2];
0163 #endif
0164     if (id[2] >= 0) {
0165       std::pair<int, int> typm = hgcons_.tileType(layer, id[0], 0);
0166       HGCScintillatorDetId detId(id[2], layer, iz * id[0], id[1], false, 0);
0167       if (typm.first >= 0) {
0168         detId.setType(typm.first);
0169         detId.setSiPM(typm.second);
0170       }
0171       index = detId.rawId();
0172       bool debug(false);
0173       if (!indices_.empty()) {
0174         int indx = HGCalTileIndex::tileIndex(layer, id[0], id[1]);
0175         if (std::find(indices_.begin(), indices_.end(), indx) != indices_.end())
0176           debug = true;
0177       }
0178       if (debug)
0179         edm::LogVerbatim("HGCSim") << "Radius/Phi " << id[0] << ":" << id[1] << " Type " << id[2] << ":" << typm.first
0180                                    << " SiPM " << typm.second << ":" << hgcons_.tileSiPM(typm.second) << " Layer "
0181                                    << layer << " z " << iz << " ID " << detId << " wt " << wt << " position " << pos
0182                                    << " R " << pos.perp();
0183 #ifdef EDM_ML_DEBUG
0184     } else {
0185       edm::LogVerbatim("HGCSim") << "Radius/Phi " << id[0] << ":" << id[1] << " Type " << id[2] << " Layer|iz " << layer
0186                                  << ":" << iz << " for i/p Layer " << layer << " module " << module << " cell " << cell
0187                                  << " iz " << iz << " pos " << pos << " wt " << wt << " ERROR";
0188 #endif
0189     }
0190   }
0191 #ifdef EDM_ML_DEBUG
0192   bool matchOnly = hgcons_.waferHexagon8Module();
0193   bool debug = hgcons_.waferHexagon8File();
0194   if (debug)
0195     edm::LogVerbatim("HGCSim") << "HGCalNumberingScheme::i/p " << det_ << ":" << layer << ":" << module << ":" << cell
0196                                << ":" << iz << ":" << pos.x() << ":" << pos.y() << ":" << pos.z() << " ID " << std::hex
0197                                << index << std::dec << " wt " << wt;
0198   bool ok = checkPosition(index, pos, matchOnly, debug);
0199   if (matchOnly && (!ok))
0200     edm::LogVerbatim("HGCSim") << "HGCalNumberingScheme::i/p " << det_ << ":" << layer << ":" << module << ":" << cell
0201                                << ":" << iz << ":" << pos.x() << ":" << pos.y() << ":" << pos.z() << " ID " << std::hex
0202                                << index << std::dec << " wt " << wt << " flag " << ok << " ERROR";
0203 #endif
0204   return index;
0205 }
0206 
0207 bool HGCalNumberingScheme::checkPosition(uint32_t index, const G4ThreeVector& pos, bool matchOnly, bool debug) const {
0208   std::pair<float, float> xy;
0209   bool ok(false), iok(true);
0210   double z1(0), tolR(14.0), tolZ(1.0);
0211   int lay(-1);
0212   if (index == 0) {
0213   } else if (DetId(index).det() == DetId::HGCalHSi) {
0214     HGCSiliconDetId id = HGCSiliconDetId(index);
0215     lay = id.layer();
0216     xy = hgcons_.locateCell(
0217         id.zside(), lay, id.waferU(), id.waferV(), id.cellU(), id.cellV(), false, true, false, false, false);
0218     z1 = hgcons_.waferZ(lay, false);
0219     ok = true;
0220     tolR = 14.0;
0221     tolZ = 1.0;
0222   } else if (DetId(index).det() == DetId::HGCalHSc) {
0223     HGCScintillatorDetId id = HGCScintillatorDetId(index);
0224     lay = id.layer();
0225     xy = hgcons_.locateCellTrap(id.zside(), lay, id.ietaAbs(), id.iphi(), false, false);
0226     z1 = hgcons_.waferZ(lay, false);
0227     ok = true;
0228     tolR = 50.0;
0229     tolZ = 5.0;
0230   }
0231   if (ok) {
0232     double r1 = std::sqrt(xy.first * xy.first + xy.second * xy.second);
0233     double r2 = pos.perp();
0234     double z2 = std::abs(pos.z());
0235     std::pair<double, double> zrange = hgcons_.rangeZ(false);
0236     std::pair<double, double> rrange = hgcons_.rangeR(z2, false);
0237     bool match = (std::abs(r1 - r2) < tolR) && (std::abs(z1 - z2) < tolZ);
0238     bool inok = ((r2 >= rrange.first) && (r2 <= rrange.second) && (z2 >= zrange.first) && (z2 <= zrange.second));
0239     bool outok = ((r1 >= rrange.first) && (r1 <= rrange.second) && (z1 >= zrange.first) && (z1 <= zrange.second));
0240     std::string ck = (((r1 < rrange.first - tolR) || (r1 > rrange.second + tolR) || (z1 < zrange.first - tolZ) ||
0241                        (z1 > zrange.second + tolZ))
0242                           ? "***** ERROR *****"
0243                           : "");
0244     if (matchOnly && match)
0245       ck = "";
0246     if (!ck.empty())
0247       iok = false;
0248     if (!(match && inok && outok) || debug) {
0249       edm::LogVerbatim("HGCSim") << "HGCalNumberingScheme::Detector " << det_ << " Layer " << lay << " R " << r2 << ":"
0250                                  << r1 << ":" << rrange.first << ":" << rrange.second << " Z " << z2 << ":" << z1 << ":"
0251                                  << zrange.first << ":" << zrange.second << " Match " << match << ":" << inok << ":"
0252                                  << outok << " " << ck;
0253       edm::LogVerbatim("HGCSim") << "Original " << pos.x() << ":" << pos.y() << " return " << xy.first << ":"
0254                                  << xy.second;
0255       if ((DetId(index).det() == DetId::HGCalEE) || (DetId(index).det() == DetId::HGCalHSi)) {
0256         int zside = (pos.z() > 0) ? 1 : -1;
0257         double wt(0), xx(zside * pos.x());
0258         int waferU, waferV, cellU, cellV, waferType;
0259         hgcons_.waferFromPosition(xx, pos.y(), zside, lay, waferU, waferV, cellU, cellV, waferType, wt, false, true);
0260         xy = hgcons_.locateCell(zside, lay, waferU, waferV, cellU, cellV, false, true, false, false, true);
0261         double dx = (xx - xy.first);
0262         double dy = (pos.y() - xy.second);
0263         double dR = std::sqrt(dx * dx + dy * dy);
0264         if (dR > tolR) {
0265           ck = " ***** ERROR *****";
0266           iok = false;
0267         } else {
0268           ck = "";
0269         }
0270         edm::LogVerbatim("HGCSim") << "HGCalNumberingScheme " << HGCSiliconDetId(index) << " original position " << xx
0271                                    << ":" << pos.y() << " derived " << xy.first << ":" << xy.second << " Difference "
0272                                    << dR << ck;
0273       }
0274     }
0275   }
0276   return iok;
0277 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!