CaloGeometryTools/src/CaloGeometryHelper.cc

0001 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0002
0003 //#include "FWCore/ParameterSet/interface/ParameterSet.h"
0004
0005 // needed for the debugging
0006 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0007 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0008 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0009 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0010 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0011
0012 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0013 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0014 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0015 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0016 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0017 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0018
0019 #include "FastSimulation/CaloGeometryTools/interface/DistanceToCell.h"
0020 #include "FastSimulation/CaloGeometryTools/interface/Crystal.h"
0021
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023
0024 #include <algorithm>
0025
0026 CaloGeometryHelper::CaloGeometryHelper() : Calorimeter() {
0027   neighbourmapcalculated_ = false;
0028   psLayer1Z_ = 303;
0029   psLayer2Z_ = 307;
0030 }
0031
0032 CaloGeometryHelper::CaloGeometryHelper(const edm::ParameterSet& fastCalo) : Calorimeter(fastCalo) {
0033   //  std::cout << " In the constructor with ParameterSet " << std::endl;
0034   psLayer1Z_ = 303;
0035   psLayer2Z_ = 307;
0036 }
0037
0038 void CaloGeometryHelper::initialize(double bField) {
0039   buildCrystalArray();
0040   buildNeighbourArray();
0041   bfield_ = bField;
0042   preshowerPresent_ = (getEcalPreshowerGeometry() != nullptr);
0043
0044   if (preshowerPresent_) {
0045     ESDetId cps1(getEcalPreshowerGeometry()->getClosestCellInPlane(GlobalPoint(80., 80., 303.), 1));
0046     psLayer1Z_ = getEcalPreshowerGeometry()->getGeometry(cps1)->getPosition().z();
0047     ESDetId cps2(getEcalPreshowerGeometry()->getClosestCellInPlane(GlobalPoint(80., 80., 307.), 2));
0048     psLayer2Z_ = getEcalPreshowerGeometry()->getGeometry(cps2)->getPosition().z();
0049     LogDebug("CaloGeometryTools") << " Preshower layer positions " << psLayer1Z_ << " " << psLayer2Z_ << std::endl;
0050   } else
0051     LogDebug("CaloGeometryTools") << " No preshower present" << std::endl;
0052
0053   //  std::cout << " Preshower layer positions " << psLayer1Z_ << " " << psLayer2Z_ << std::endl;
0054 }
0055
0056 CaloGeometryHelper::~CaloGeometryHelper() { ; }
0057
0058 DetId CaloGeometryHelper::getClosestCell(const XYZPoint& point, bool ecal, bool central) const {
0059   DetId result;
0060   if (ecal) {
0061     if (central) {
0062       //      std::cout << "EcalBarrelGeometry_" << " " << EcalBarrelGeometry_ << std::endl;
0063       result = EcalBarrelGeometry_->getClosestCell(GlobalPoint(point.X(), point.Y(), point.Z()));
0064 #ifdef DEBUGGCC
0065       if (result.null())
0066         return result;
0067       GlobalPoint ip = GlobalPoint(point.X(), point.Y(), point.Z());
0068       GlobalPoint cc = EcalBarrelGeometry_->getGeometry(result)->getPosition();
0069       float deltaeta2 = ip.eta() - cc.eta();
0070       deltaeta2 *= deltaeta2;
0071       float deltaphi2 = acos(cos(ip.phi() - cc.phi()));
0072       deltaphi2 *= deltaphi2;
0073       Histos::instance()->fill("h100", point.eta(), sqrt(deltaeta2 + deltaphi2));
0074 #endif
0075     } else {
0076       result = EcalEndcapGeometry_->getClosestCell(GlobalPoint(point.X(), point.Y(), point.Z()));
0077 #ifdef DEBUGGCC
0078       if (result.null()) {
0079         return result;
0080       }
0081       GlobalPoint ip = GlobalPoint(point.X(), point.Y(), point.Z());
0082       GlobalPoint cc = EcalEndcapGeometry_->getGeometry(result)->getPosition();
0083       Histos::instance()->fill("h110", point.eta(), (ip - cc).perp());
0084 #endif
0085     }
0086   } else {
0087     result = static_cast<const HcalGeometry*>(HcalGeometry_)
0088                  ->getClosestCell(GlobalPoint(point.X(), point.Y(), point.Z()), true);
0089     HcalDetId myDetId(result);
0090
0091     // special patch for HF (this is already a part of HcalGeometry)
0092     if (myDetId.subdetId() == HcalForward) {
0093       int mylayer;
0094       if (fabs(point.Z()) > 1132.) {
0095         mylayer = 2;
0096       } else {
0097         mylayer = 1;
0098       }
0099       HcalDetId myDetId2((HcalSubdetector)myDetId.subdetId(), myDetId.ieta(), myDetId.iphi(), mylayer);
0100       result = myDetId2;
0101       //      return result;
0102     }
0103
0104     // Special patch to correct the HCAL geometry (does not work)
0105     /*
0106       if(result.subdetId()!=HcalEndcap) return result;
0107       if(myDetId.depth()==3) return result;
0108
0109       int ieta=myDetId.ietaAbs();
0110       float azmin=400.458;         /// in sync with BaseParticlePropagator
0111
0112       if(ieta<=17)
0113         return result;
0114       else if(ieta>=18 && ieta<=26)
0115         azmin += 35.0;    // don't consider ieta=18 nose separately
0116       else if(ieta>=27)
0117         azmin += 21.0;
0118
0119       HcalDetId first(HcalEndcap,myDetId.ieta(),myDetId.iphi(),1);
0120       bool layer2=(fabs(point.Z())>azmin);
0121       if(!layer2)
0122         {
0123           return first;
0124         }
0125       else
0126         {
0127           HcalDetId second(HcalEndcap,myDetId.ieta(),myDetId.iphi(),2);
0128       if(second!=HcalDetId()) result=second;
0129     }
0130       */
0131 #ifdef DEBUGGCC
0132     if (result.null()) {
0133       return result;
0134     }
0135     GlobalPoint ip = GlobalPoint(point.x(), point.y(), point.z());
0136     GlobalPoint cc = HcalGeometry_->getPosition(result);
0137     float deltaeta2 = ip.eta() - cc.eta();
0138     deltaeta2 *= deltaeta2;
0139     float deltaphi2 = acos(cos(ip.phi() - cc.phi()));
0140     deltaphi2 *= deltaphi2;
0141
0142     Histos::instance()->fill("h120", point.eta(), sqrt(deltaeta2 + deltaphi2));
0143 #endif
0144   }
0145   return result;
0146 }
0147
0148 void CaloGeometryHelper::getWindow(const DetId& pivot, int s1, int s2, std::vector<DetId>& vec) const {
0149   // currently the getWindow method is the same for EcalBarrelTopology and EndcapTopology
0150   // (implemented in CaloSubDetectorTopology)
0151   // optimized versions are foreseen
0152   vec = getEcalTopology(pivot.subdetId())->getWindow(pivot, s1, s2);
0153   DistanceToCell distance(getEcalGeometry(pivot.subdetId()), pivot);
0154   sort(vec.begin(), vec.end(), distance);
0155 }
0156
0157 void CaloGeometryHelper::buildCrystal(const DetId& cell, Crystal& xtal) const {
0158   if (cell.subdetId() == EcalBarrel) {
0159     xtal = Crystal(cell, &barrelCrystals_[EBDetId(cell).hashedIndex()]);
0160     return;
0161   }
0162   if (cell.subdetId() == EcalEndcap) {
0163     xtal = Crystal(cell, &endcapCrystals_[EEDetId(cell).hashedIndex()]);
0164     return;
0165   }
0166 }
0167
0168 // Build the array of (max)8 neighbors
0169 void CaloGeometryHelper::buildNeighbourArray() {
0170   static const CaloDirection orderedDir[8] = {SOUTHWEST, SOUTH, SOUTHEAST, WEST, EAST, NORTHWEST, NORTH, NORTHEAST};
0171
0172   const unsigned nbarrel = EBDetId::kSizeForDenseIndexing;
0173   // Barrel first. The hashed index runs from 0 to 61199
0174   barrelNeighbours_.resize(nbarrel);
0175
0176   //std::cout << " Building the array of neighbours (barrel) " ;
0177
0178   const std::vector<DetId>& vec(EcalBarrelGeometry_->getValidDetIds(DetId::Ecal, EcalBarrel));
0179   unsigned size = vec.size();
0180   for (unsigned ic = 0; ic < size; ++ic) {
0181     // We get the 9 cells in a square.
0182     std::vector<DetId> neighbours(EcalBarrelTopology_->getWindow(vec[ic], 3, 3));
0183     //      std::cout << " Cell " << EBDetId(vec[ic]) << std::endl;
0184     unsigned nneighbours = neighbours.size();
0185
0186     unsigned hashedindex = EBDetId(vec[ic]).hashedIndex();
0187     if (hashedindex >= nbarrel) {
0188       LogDebug("CaloGeometryTools") << " Array overflow " << std::endl;
0189     }
0190
0191     // If there are 9 cells, it is easy, and this order is know:
0192     //      6  7  8
0193     //      3  4  5
0194     //      0  1  2   (0 = SOUTHWEST)
0195
0196     if (nneighbours == 9) {
0197       //barrelNeighbours_[hashedindex].reserve(8);
0198       unsigned int nn = 0;
0199       for (unsigned in = 0; in < nneighbours; ++in) {
0200         // remove the centre
0201         if (neighbours[in] != vec[ic]) {
0202           barrelNeighbours_[hashedindex][nn] = (neighbours[in]);
0203           nn++;
0204           //          std::cout << " Neighbour " << in << " " << EBDetId(neighbours[in]) << std::endl;
0205         }
0206       }
0207     } else {
0208       DetId central(vec[ic]);
0209       //barrelNeighbours_[hashedindex].resize(8,DetId(0));
0210       for (unsigned idir = 0; idir < 8; ++idir) {
0211         DetId testid = central;
0212         bool status = move(testid, orderedDir[idir], false);
0213         if (status)
0214           barrelNeighbours_[hashedindex][idir] = testid;
0215       }
0216     }
0217   }
0218
0219   // Moved to the endcap
0220
0221   //  std::cout << " done " << size << std::endl;
0222   //  std::cout << " Building the array of neighbours (endcap) " ;
0223
0224   const std::vector<DetId>& vece(EcalEndcapGeometry_->getValidDetIds(DetId::Ecal, EcalEndcap));
0225   size = vece.size();
0226   // There are some holes in the hashedIndex for the EE. Hence the array is bigger than the number
0227   // of crystals
0228   const unsigned nendcap = EEDetId::kSizeForDenseIndexing;
0229
0230   endcapNeighbours_.resize(nendcap);
0231   for (unsigned ic = 0; ic < size; ++ic) {
0232     // We get the 9 cells in a square.
0233     std::vector<DetId> neighbours(EcalEndcapTopology_->getWindow(vece[ic], 3, 3));
0234     unsigned nneighbours = neighbours.size();
0235     // remove the centre
0236     unsigned hashedindex = EEDetId(vece[ic]).hashedIndex();
0237
0238     if (hashedindex >= nendcap) {
0239       LogDebug("CaloGeometryTools") << " Array overflow " << std::endl;
0240     }
0241
0242     if (nneighbours == 9) {
0243       //endcapNeighbours_[hashedindex].reserve(8);
0244       unsigned int nn = 0;
0245       for (unsigned in = 0; in < nneighbours; ++in) {
0246         // remove the centre
0247         if (neighbours[in] != vece[ic]) {
0248           endcapNeighbours_[hashedindex][nn] = (neighbours[in]);
0249           nn++;
0250         }
0251       }
0252     } else {
0253       DetId central(vece[ic]);
0254       //endcapNeighbours_[hashedindex].resize(8,DetId(0));
0255       for (unsigned idir = 0; idir < 8; ++idir) {
0256         DetId testid = central;
0257         bool status = move(testid, orderedDir[idir], false);
0258         if (status)
0259           endcapNeighbours_[hashedindex][idir] = testid;
0260       }
0261     }
0262   }
0263   //  std::cout << " done " << size <<std::endl;
0264   neighbourmapcalculated_ = true;
0265 }
0266
0267 const CaloGeometryHelper::NeiVect& CaloGeometryHelper::getNeighbours(const DetId& detid) const {
0268   return (detid.subdetId() == EcalBarrel) ? barrelNeighbours_[EBDetId(detid).hashedIndex()]
0269                                           : endcapNeighbours_[EEDetId(detid).hashedIndex()];
0270 }
0271
0272 bool CaloGeometryHelper::move(DetId& cell, const CaloDirection& dir, bool fast) const {
0273   DetId originalcell = cell;
0274   if (dir == NONE || cell == DetId(0))
0275     return false;
0276
0277   // Conversion CaloDirection and index in the table
0278   // CaloDirection :NONE,SOUTH,SOUTHEAST,SOUTHWEST,EAST,WEST, NORTHEAST,NORTHWEST,NORTH
0279   // Table : SOUTHWEST,SOUTH,SOUTHEAST,WEST,EAST,NORTHWEST,NORTH, NORTHEAST
0280   static const int calodirections[9] = {-1, 1, 2, 0, 4, 3, 7, 5, 6};
0281
0282   if (fast && neighbourmapcalculated_) {
0283     DetId result = (originalcell.subdetId() == EcalBarrel)
0284                        ? barrelNeighbours_[EBDetId(originalcell).hashedIndex()][calodirections[dir]]
0285                        : endcapNeighbours_[EEDetId(originalcell).hashedIndex()][calodirections[dir]];
0286     bool status = !result.null();
0287     cell = result;
0288     return status;
0289   }
0290
0291   if (dir == NORTH || dir == SOUTH || dir == EAST || dir == WEST) {
0292     return simplemove(cell, dir);
0293   } else {
0294     if (dir == NORTHEAST || dir == NORTHWEST || dir == SOUTHEAST || dir == SOUTHWEST)
0295       return diagonalmove(cell, dir);
0296   }
0297
0298   cell = DetId(0);
0299   return false;
0300 }
0301
0302 bool CaloGeometryHelper::simplemove(DetId& cell, const CaloDirection& dir) const {
0303   std::vector<DetId> neighbours;
0304   if (cell.subdetId() == EcalBarrel)
0305     neighbours = EcalBarrelTopology_->getNeighbours(cell, dir);
0306   else if (cell.subdetId() == EcalEndcap)
0307     neighbours = EcalEndcapTopology_->getNeighbours(cell, dir);
0308
0309   if ((!neighbours.empty()) && (!neighbours[0].null())) {
0310     cell = neighbours[0];
0311     return true;
0312   } else {
0313     cell = DetId(0);
0314     return false;
0315   }
0316 }
0317
0318 bool CaloGeometryHelper::diagonalmove(DetId& cell, const CaloDirection& dir) const {
0319   bool result;
0320   // One has to try both paths
0321   if (dir == NORTHEAST) {
0322     result = simplemove(cell, NORTH);
0323     if (result)
0324       return simplemove(cell, EAST);
0325     else {
0326       result = simplemove(cell, EAST);
0327       if (result)
0328         return simplemove(cell, NORTH);
0329       else
0330         return false;
0331     }
0332   } else if (dir == NORTHWEST) {
0333     result = simplemove(cell, NORTH);
0334     if (result)
0335       return simplemove(cell, WEST);
0336     else {
0337       result = simplemove(cell, WEST);
0338       if (result)
0339         return simplemove(cell, NORTH);
0340       else
0341         return false;
0342     }
0343   } else if (dir == SOUTHEAST) {
0344     result = simplemove(cell, SOUTH);
0345     if (result)
0346       return simplemove(cell, EAST);
0347     else {
0348       result = simplemove(cell, EAST);
0349       if (result)
0350         return simplemove(cell, SOUTH);
0351       else
0352         return false;
0353     }
0354   } else if (dir == SOUTHWEST) {
0355     result = simplemove(cell, SOUTH);
0356     if (result)
0357       return simplemove(cell, WEST);
0358     else {
0359       result = simplemove(cell, SOUTH);
0360       if (result)
0361         return simplemove(cell, WEST);
0362       else
0363         return false;
0364     }
0365   }
0366   cell = DetId(0);
0367   return false;
0368 }
0369
0370 bool CaloGeometryHelper::borderCrossing(const DetId& c1, const DetId& c2) const {
0371   if (c1.subdetId() != c2.subdetId())
0372     return false;
0373
0374   if (c1.subdetId() == EcalBarrel) {
0375     // there is a crack if the two cells don't belong to the same
0376     // module
0377     EBDetId cc1(c1);
0378     EBDetId cc2(c2);
0379     return (cc1.im() != cc2.im() || cc1.ism() != cc2.ism());
0380   }
0381
0382   if (c1.subdetId() == EcalEndcap) {
0383     // there is a crack if the two cells don't belong to the same
0384     // module
0385     return (EEDetId(c1).isc() != EEDetId(c2).isc());
0386   }
0387   return false;
0388 }
0389
0390 void CaloGeometryHelper::buildCrystalArray() {
0391   const unsigned nbarrel = EBDetId::kSizeForDenseIndexing;
0392   // Barrel first. The hashed index runs from 0 to 61199
0393   barrelCrystals_.resize(nbarrel, BaseCrystal());
0394
0395   //std::cout << " Building the array of crystals (barrel) " ;
0396   const std::vector<DetId>& vec(EcalBarrelGeometry_->getValidDetIds(DetId::Ecal, EcalBarrel));
0397   unsigned size = vec.size();
0398   for (unsigned ic = 0; ic < size; ++ic) {
0399     unsigned hashedindex = EBDetId(vec[ic]).hashedIndex();
0400     auto geom = EcalBarrelGeometry_->getGeometry(vec[ic]);
0401     BaseCrystal xtal(vec[ic]);
0402     xtal.setCorners(geom->getCorners(), geom->getPosition());
0403     barrelCrystals_[hashedindex] = xtal;
0404   }
0405
0406   //  std::cout << " done " << size << std::endl;
0407   //  std::cout << " Building the array of crystals (endcap) " ;
0408
0409   const std::vector<DetId>& vece(EcalEndcapGeometry_->getValidDetIds(DetId::Ecal, EcalEndcap));
0410   size = vece.size();
0411   // There are some holes in the hashedIndex for the EE. Hence the array is bigger than the number
0412   // of crystals
0413   const unsigned nendcap = EEDetId::kSizeForDenseIndexing;
0414
0415   endcapCrystals_.resize(nendcap, BaseCrystal());
0416   for (unsigned ic = 0; ic < size; ++ic) {
0417     unsigned hashedindex = EEDetId(vece[ic]).hashedIndex();
0418     auto geom = EcalEndcapGeometry_->getGeometry(vece[ic]);
0419     BaseCrystal xtal(vece[ic]);
0420     xtal.setCorners(geom->getCorners(), geom->getPosition());
0421     endcapCrystals_[hashedindex] = xtal;
0422   }
0423   //  std::cout << " done " << size << std::endl;
0424 }