Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:10 

0001 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0002 #include "Calibration/IsolatedParticles/interface/FindEtaPhi.h"
0003 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
0004 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0005 
0006 #include <algorithm>
0007 #include <iostream>
0008 
0009 namespace spr {
0010 
0011   template <typename T>
0012   std::vector<std::pair<DetId, double> > eHCALmatrixCell(const HcalTopology* topology,
0013                                                          const DetId& det,
0014                                                          edm::Handle<T>& hits,
0015                                                          int ieta,
0016                                                          int iphi,
0017                                                          bool includeHO,
0018                                                          double hbThr,
0019                                                          double heThr,
0020                                                          double hfThr,
0021                                                          double hoThr,
0022                                                          bool debug) {
0023     std::vector<DetId> dets(1, det);
0024     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, debug);
0025     return spr::energyDetIdHCAL(vdets, hits, hbThr, heThr, hfThr, hoThr, debug);
0026   }
0027 
0028   //*** Returns energy in NxN Hcal towers
0029   template <typename T>
0030   std::pair<double, int> eHCALmatrixTotal(const HcalTopology* topology,
0031                                           const DetId& det0,
0032                                           edm::Handle<T>& hits,
0033                                           int ieta,
0034                                           int iphi,
0035                                           bool includeHO,
0036                                           double hbThr,
0037                                           double heThr,
0038                                           double hfThr,
0039                                           double hoThr,
0040                                           bool debug) {
0041     HcalDetId hcid0(det0.rawId());
0042     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0043     DetId det(hcid.rawId());
0044     if (debug)
0045       edm::LogVerbatim("IsoTrack") << "Inside eHCALmatrixTotal " << ieta << "X" << iphi << " Inclusion of HO Flag "
0046                                    << includeHO;
0047 
0048     spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug);
0049 
0050     // Get maximum of all the trials
0051     double energySum = 0;
0052     int itrym = 0;
0053     for (int itry = 0; itry < etaphi.ntrys; itry++) {
0054       double energy = energyHCALmatrixTotal(topology,
0055                                             det,
0056                                             hits,
0057                                             etaphi.ietaE[itry],
0058                                             etaphi.ietaW[itry],
0059                                             etaphi.iphiN[itry],
0060                                             etaphi.iphiS[itry],
0061                                             includeHO,
0062                                             hbThr,
0063                                             heThr,
0064                                             hfThr,
0065                                             hoThr,
0066                                             debug);
0067       if (energy > energySum) {
0068         energySum = energy;
0069         itrym = itry;
0070       }
0071       if (debug)
0072         edm::LogVerbatim("IsoTrack") << "eHCALMatrix::Total energy " << energy << " Max " << energySum << "in trial # "
0073                                      << itrym;
0074     }
0075 
0076     return std::pair<double, int>(energySum, itrym);
0077   }
0078 
0079   //*** Returns vector of hits in 3x3 or 5x5 Hcal towers
0080   template <typename T>
0081   double energyHCALmatrix(const HcalTopology* topology,
0082                           const DetId& det0,
0083                           edm::Handle<T>& hits,
0084                           int ieta,
0085                           int iphi,
0086                           bool includeHO,
0087                           double hbThr,
0088                           double heThr,
0089                           double hfThr,
0090                           double hoThr,
0091                           bool debug) {
0092     HcalDetId hcid0(det0.rawId());
0093     HcalDetId hcid(hcid0.subdet(), hcid0.ieta(), hcid0.iphi(), 1);
0094     DetId det(hcid.rawId());
0095     if (debug)
0096       edm::LogVerbatim("IsoTrack") << "Inside energyHCALmatrix " << 2 * ieta + 1 << "X" << 2 * iphi + 1
0097                                    << " Inclusion of HO Flag " << includeHO;
0098 
0099     double energy = 0;
0100     if (ieta > 2) {
0101       edm::LogVerbatim("IsoTrack") << " no matrix > 5x5 is coded ! ";
0102       return energy;
0103     }
0104     // used dets during building the matrix
0105     std::vector<DetId> dets;
0106     dets.clear();
0107 
0108     // central tower
0109     std::vector<DetId> vNeighboursDetId(1, det);
0110     energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0111 
0112     if (ieta > 0) {
0113       // go to east from central tower
0114       vNeighboursDetId.clear();
0115       if (debug)
0116         edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: det " << (HcalDetId)det;
0117       vNeighboursDetId = topology->east(det);
0118       if (debug) {
0119         std::ostringstream st1;
0120         st1 << "Neighbour:: East";
0121         for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
0122           st1 << " " << (HcalDetId)vNeighboursDetId[i];
0123         edm::LogVerbatim("IsoTrack") << st1.str();
0124       }
0125       energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0126       if (ieta == 2) {
0127         for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
0128           std::vector<DetId> vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]);
0129           if (debug) {
0130             std::ostringstream st1;
0131             st1 << "Neighbour:: East";
0132             for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
0133               st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
0134             edm::LogVerbatim("IsoTrack") << st1.str();
0135           }
0136           energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0137         }
0138       }
0139 
0140       // go to west from central tower
0141       vNeighboursDetId.clear();
0142       vNeighboursDetId = topology->west(det);
0143       if (debug) {
0144         std::ostringstream st1;
0145         st1 << "Neighbour:: West";
0146         for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
0147           st1 << " " << (HcalDetId)vNeighboursDetId[i];
0148         edm::LogVerbatim("IsoTrack") << st1.str();
0149       }
0150       energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0151       if (ieta == 2) {
0152         for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
0153           std::vector<DetId> vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]);
0154           if (debug) {
0155             std::ostringstream st1;
0156             st1 << "Neighbour:: West";
0157             for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
0158               st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
0159             edm::LogVerbatim("IsoTrack") << st1.str();
0160           }
0161           energy += energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0162         }
0163       }
0164     }
0165 
0166     // do steps in phi to north
0167     DetId detnorth = det;
0168     for (int inorth = 0; inorth < iphi; inorth++) {
0169       if (debug)
0170         edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: detNorth " << (HcalDetId)detnorth;
0171       std::vector<DetId> NorthDetId = topology->north(detnorth);
0172       if (debug) {
0173         std::ostringstream st1;
0174         st1 << "Neighbour:: North";
0175         for (unsigned int i = 0; i < NorthDetId.size(); i++)
0176           st1 << " " << (HcalDetId)NorthDetId[i];
0177         edm::LogVerbatim("IsoTrack") << st1.str();
0178       }
0179       if (!(NorthDetId.empty())) {
0180         energy += energyHCAL(NorthDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0181         if (ieta > 0) {
0182           // go to east
0183           vNeighboursDetId.clear();
0184           vNeighboursDetId = topology->east(NorthDetId[0]);
0185           if (debug) {
0186             std::ostringstream st1;
0187             st1 << "Neighbour:: East";
0188             for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
0189               st1 << " " << (HcalDetId)vNeighboursDetId[i];
0190             edm::LogVerbatim("IsoTrack") << st1.str();
0191           }
0192           energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0193           if (ieta == 2) {
0194             for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
0195               std::vector<DetId> vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]);
0196               if (debug) {
0197                 std::ostringstream st1;
0198                 st1 << "Neighbour:: East";
0199                 for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
0200                   st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
0201                 edm::LogVerbatim("IsoTrack") << st1.str();
0202               }
0203               energy +=
0204                   energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0205             }
0206           }
0207 
0208           // go to west
0209           vNeighboursDetId.clear();
0210           vNeighboursDetId = topology->west(NorthDetId[0]);
0211           if (debug) {
0212             std::ostringstream st1;
0213             st1 << "Neighbour:: West";
0214             for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
0215               st1 << " " << (HcalDetId)vNeighboursDetId[i];
0216             edm::LogVerbatim("IsoTrack") << st1.str();
0217           }
0218           energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0219           if (ieta == 2) {
0220             for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
0221               std::vector<DetId> vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]);
0222               if (debug) {
0223                 std::ostringstream st1;
0224                 st1 << "Neighbour:: West";
0225                 for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
0226                   st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
0227                 edm::LogVerbatim("IsoTrack") << st1.str();
0228               }
0229               energy +=
0230                   energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0231             }
0232           }
0233         }
0234         detnorth = NorthDetId[0];
0235       } else {
0236         break;
0237       }
0238     }
0239 
0240     // do steps in phi to south
0241     DetId detsouth = det;
0242     for (int isouth = 0; isouth < iphi; isouth++) {
0243       if (debug)
0244         edm::LogVerbatim("IsoTrack") << "hitHCALMatrix:: detSouth " << (HcalDetId)detsouth;
0245       vNeighboursDetId.clear();
0246       std::vector<DetId> SouthDetId = topology->south(detsouth);
0247       if (!(SouthDetId.empty())) {
0248         energy += energyHCAL(SouthDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0249         if (debug) {
0250           std::ostringstream st1;
0251           st1 << "Neighbour:: South";
0252           for (unsigned int i = 0; i < SouthDetId.size(); i++)
0253             st1 << " " << (HcalDetId)SouthDetId[i];
0254           edm::LogVerbatim("IsoTrack") << st1.str();
0255         }
0256 
0257         // go to east
0258         if (ieta > 0) {
0259           // go to east from central tower
0260           vNeighboursDetId.clear();
0261           vNeighboursDetId = topology->east(SouthDetId[0]);
0262           if (debug) {
0263             std::ostringstream st1;
0264             st1 << "Neighbour:: East";
0265             for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
0266               st1 << " " << (HcalDetId)vNeighboursDetId[i];
0267             edm::LogVerbatim("IsoTrack") << st1.str();
0268           }
0269           energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0270           if (ieta == 2) {
0271             for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
0272               std::vector<DetId> vNeighboursDetIdc = topology->east(vNeighboursDetId[ii]);
0273               if (debug) {
0274                 std::ostringstream st1;
0275                 st1 << "Neighbour:: East";
0276                 for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
0277                   st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
0278                 edm::LogVerbatim("IsoTrack") << st1.str();
0279               }
0280               energy +=
0281                   energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0282             }
0283           }
0284 
0285           // go to west
0286           vNeighboursDetId.clear();
0287           vNeighboursDetId = topology->west(SouthDetId[0]);
0288           if (debug) {
0289             std::ostringstream st1;
0290             st1 << "Neighbour:: West";
0291             for (unsigned int i = 0; i < vNeighboursDetId.size(); i++)
0292               st1 << " " << (HcalDetId)vNeighboursDetId[i];
0293             edm::LogVerbatim("IsoTrack") << st1.str();
0294           }
0295           energy += energyHCAL(vNeighboursDetId, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0296           if (ieta == 2) {
0297             for (int ii = 0; ii < (int)vNeighboursDetId.size(); ii++) {
0298               std::vector<DetId> vNeighboursDetIdc = topology->west(vNeighboursDetId[ii]);
0299               if (debug) {
0300                 std::ostringstream st1;
0301                 st1 << "Neighbour:: West";
0302                 for (unsigned int i = 0; i < vNeighboursDetIdc.size(); i++)
0303                   st1 << " " << (HcalDetId)vNeighboursDetIdc[i];
0304                 edm::LogVerbatim("IsoTrack") << st1.str();
0305               }
0306               energy +=
0307                   energyHCAL(vNeighboursDetIdc, dets, topology, hits, includeHO, hbThr, heThr, hfThr, hoThr, debug);
0308             }
0309           }
0310         }
0311         detsouth = SouthDetId[0];
0312       } else {
0313         break;
0314       }
0315     }
0316 
0317     return energy;
0318   }
0319 
0320   template <typename T>
0321   double energyHCAL(std::vector<DetId>& vNeighboursDetId,
0322                     std::vector<DetId>& dets,
0323                     const HcalTopology* topology,
0324                     edm::Handle<T>& hits,
0325                     bool includeHO,
0326                     double hbThr,
0327                     double heThr,
0328                     double hfThr,
0329                     double hoThr,
0330                     bool debug) {
0331     double energySum = 0;
0332     int kHit = 0;
0333     for (int i = 0; i < (int)vNeighboursDetId.size(); i++) {
0334       int n = std::count(dets.begin(), dets.end(), vNeighboursDetId[i]);
0335       if (n != 0)
0336         continue;
0337       dets.push_back(vNeighboursDetId[i]);
0338       std::vector<typename T::const_iterator> hit = spr::findHit(hits, vNeighboursDetId[i]);
0339 
0340       double energy = 0;
0341       int subdet = ((HcalDetId)(vNeighboursDetId[i].rawId())).subdet();
0342       double eThr = hbThr;
0343       if (subdet == (int)(HcalEndcap))
0344         eThr = heThr;
0345       else if (subdet == (int)(HcalForward))
0346         eThr = hfThr;
0347       else if (subdet == (int)(HcalOuter))
0348         eThr = hoThr;
0349       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0350         if (hit[ihit] != hits->end()) {
0351           energy += hit[ihit]->energy();
0352           kHit++;
0353         }
0354       }
0355       if (energy > eThr)
0356         energySum += energy;
0357 
0358       HcalDetId depth = vNeighboursDetId[i];
0359       // max. three depths can be in endcap and we go to 2nd and 3rd from 1st where we are now
0360       for (int idepth = 0; idepth < 2; idepth++) {
0361         std::vector<DetId> vUpDetId = topology->up(depth);
0362         if (!(vUpDetId.empty())) {
0363           if (includeHO || vUpDetId[0].subdetId() != (int)(HcalOuter)) {
0364             int n = std::count(dets.begin(), dets.end(), vUpDetId[0]);
0365             if (n == 0) {
0366               dets.push_back(vUpDetId[0]);
0367               hit = spr::findHit(hits, vUpDetId[0]);
0368               energy = 0;
0369               subdet = ((HcalDetId)(vUpDetId[0].rawId())).subdet();
0370               eThr = hbThr;
0371               if (subdet == (int)(HcalEndcap))
0372                 eThr = heThr;
0373               else if (subdet == (int)(HcalForward))
0374                 eThr = hfThr;
0375               else if (subdet == (int)(HcalOuter))
0376                 eThr = hoThr;
0377               for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0378                 if (hit[ihit] != hits->end()) {
0379                   energy += hit[ihit]->energy();
0380                   kHit++;
0381                 }
0382               }
0383               if (energy > eThr)
0384                 energySum += energy;
0385             }
0386           }
0387           depth = vUpDetId[0];
0388         }
0389       }
0390     }
0391     if (debug)
0392       edm::LogVerbatim("IsoTrack") << "energyHCAL:: Energy " << energySum << " from " << kHit << " hits";
0393     return energySum;
0394   }
0395 
0396   template <typename T>
0397   std::vector<std::pair<DetId, double> > energyDetIdHCAL(std::vector<DetId>& vdets,
0398                                                          edm::Handle<T>& hits,
0399                                                          double hbThr,
0400                                                          double heThr,
0401                                                          double hfThr,
0402                                                          double hoThr,
0403                                                          bool debug) {
0404     std::vector<std::pair<DetId, double> > energyDetId;
0405     int khit = 0;
0406     unsigned int maxsize = vdets.size();
0407     for (unsigned int i = 0; i < maxsize; i++) {
0408       std::vector<typename T::const_iterator> hit = spr::findHit(hits, vdets[i]);
0409       double energy = 0;
0410       bool ok = false;
0411       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0412         if (hit[ihit] != hits->end()) {
0413           energy += hit[ihit]->energy();
0414           khit++;
0415           ok = true;
0416           if (debug)
0417             edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours:: Hit " << khit << " " << (HcalDetId)vdets[i]
0418                                          << " energy " << hit[ihit]->energy();
0419         }
0420       }
0421       if (debug)
0422         edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours::Detector " << i << " " << (HcalDetId)vdets[i]
0423                                      << " energy " << energy;
0424       if (ok) {
0425         int subdet = ((HcalDetId)(vdets[i].rawId())).subdet();
0426         double eThr = hbThr;
0427         if (subdet == (int)(HcalEndcap))
0428           eThr = heThr;
0429         else if (subdet == (int)(HcalForward))
0430           eThr = hfThr;
0431         else if (subdet == (int)(HcalOuter))
0432           eThr = hoThr;
0433         if (energy > eThr)
0434           energyDetId.push_back(std::pair<DetId, double>(vdets[i], energy));
0435       }
0436     }
0437 
0438     if (debug)
0439       edm::LogVerbatim("IsoTrack") << "energyDetIdHCALneighbours::EnergyDetId buffer size " << energyDetId.size();
0440     return energyDetId;
0441   }
0442 
0443 }  // namespace spr
[ 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!