Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:31:37

0001 #include <map>
0002 #include "TrackingTools/TrackAssociator/interface/TrackDetMatchInfo.h"
0003 #include "DetIdInfo.h"
0004 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0005 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0006 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0007 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0008 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "DataFormats/Math/interface/LorentzVector.h"
0011 #include "Math/VectorUtil.h"
0012 #include <algorithm>
0013 
0014 ///////////////////////////
0015 
0016 std::string TrackDetMatchInfo::dumpGeometry(const DetId& id) {
0017   if (!caloGeometry || !caloGeometry->getSubdetectorGeometry(id) ||
0018       !caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)) {
0019     throw cms::Exception("FatalError") << "Failed to access geometry for DetId: " << id.rawId();
0020   }
0021   std::ostringstream oss;
0022 
0023   const CaloCellGeometry::CornersVec& points = caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)->getCorners();
0024   for (CaloCellGeometry::CornersVec::const_iterator point = points.begin(); point != points.end(); ++point)
0025     oss << "(" << point->z() << ", " << point->perp() << ", " << point->eta() << ", " << point->phi() << "), \t";
0026   return oss.str();
0027 }
0028 
0029 GlobalPoint TrackDetMatchInfo::getPosition(const DetId& id) {
0030   // this part might be slow
0031   if (!caloGeometry || !caloGeometry->getSubdetectorGeometry(id) ||
0032       !caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)) {
0033     throw cms::Exception("FatalError") << "Failed to access geometry for DetId: " << id.rawId();
0034     return GlobalPoint(0, 0, 0);
0035   }
0036   return caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)->getPosition();
0037 }
0038 
0039 double TrackDetMatchInfo::crossedEnergy(EnergyType type) {
0040   double energy(0);
0041   switch (type) {
0042     case EcalRecHits: {
0043       for (std::vector<const EcalRecHit*>::const_iterator hit = crossedEcalRecHits.begin();
0044            hit != crossedEcalRecHits.end();
0045            hit++)
0046         energy += (*hit)->energy();
0047     } break;
0048     case HcalRecHits: {
0049       for (std::vector<const HBHERecHit*>::const_iterator hit = crossedHcalRecHits.begin();
0050            hit != crossedHcalRecHits.end();
0051            hit++)
0052         energy += (*hit)->energy();
0053     } break;
0054     case HORecHits: {
0055       for (std::vector<const HORecHit*>::const_iterator hit = crossedHORecHits.begin(); hit != crossedHORecHits.end();
0056            hit++)
0057         energy += (*hit)->energy();
0058     } break;
0059     case TowerTotal: {
0060       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0061         energy += (*hit)->energy();
0062     } break;
0063     case TowerEcal: {
0064       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0065         energy += (*hit)->emEnergy();
0066     } break;
0067     case TowerHcal: {
0068       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0069         energy += (*hit)->hadEnergy();
0070     } break;
0071     case TowerHO: {
0072       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0073         energy += (*hit)->outerEnergy();
0074     } break;
0075     default:
0076       throw cms::Exception("FatalError") << "Unknown calo energy type: " << type;
0077   }
0078   return energy;
0079 }
0080 
0081 bool TrackDetMatchInfo::insideCone(const DetId& id, const double dR) {
0082   GlobalPoint idPosition = getPosition(id);
0083   if (idPosition.mag() < 0.01)
0084     return false;
0085 
0086   math::XYZVector idPositionRoot(idPosition.x(), idPosition.y(), idPosition.z());
0087   math::XYZVector trackP3(stateAtIP.momentum().x(), stateAtIP.momentum().y(), stateAtIP.momentum().z());
0088   return ROOT::Math::VectorUtil::DeltaR(trackP3, idPositionRoot) < dR;
0089 }
0090 
0091 double TrackDetMatchInfo::coneEnergy(double dR, EnergyType type) {
0092   double energy(0);
0093   switch (type) {
0094     case EcalRecHits: {
0095       for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end(); hit++)
0096         if (insideCone((*hit)->detid(), dR))
0097           energy += (*hit)->energy();
0098     } break;
0099     case HcalRecHits: {
0100       for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++)
0101         if (insideCone((*hit)->detid(), dR))
0102           energy += (*hit)->energy();
0103     } break;
0104     case HORecHits: {
0105       for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++)
0106         if (insideCone((*hit)->detid(), dR))
0107           energy += (*hit)->energy();
0108     } break;
0109     case TowerTotal: {
0110       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0111         if (insideCone((*hit)->id(), dR))
0112           energy += (*hit)->energy();
0113     } break;
0114     case TowerEcal: {
0115       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0116         if (insideCone((*hit)->id(), dR))
0117           energy += (*hit)->emEnergy();
0118     } break;
0119     case TowerHcal: {
0120       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0121         if (insideCone((*hit)->id(), dR))
0122           energy += (*hit)->hadEnergy();
0123     } break;
0124     case TowerHO: {
0125       for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
0126         if (insideCone((*hit)->id(), dR))
0127           energy += (*hit)->outerEnergy();
0128     } break;
0129     default:
0130       throw cms::Exception("FatalError") << "Unknown calo energy type: " << type;
0131   }
0132   return energy;
0133 }
0134 
0135 //////////////////////////////////////////////////
0136 
0137 double TrackDetMatchInfo::nXnEnergy(const DetId& id, EnergyType type, int gridSize) {
0138   double energy(0);
0139   if (id.rawId() == 0)
0140     return 0.;
0141   switch (type) {
0142     case TowerTotal:
0143     case TowerHcal:
0144     case TowerEcal:
0145     case TowerHO: {
0146       if (id.det() != DetId::Calo) {
0147         throw cms::Exception("FatalError") << "Wrong DetId. Expected CaloTower, but found:\n"
0148                                            << DetIdInfo::info(id, nullptr) << "\n";
0149       }
0150       CaloTowerDetId centerId(id);
0151       for (std::vector<const CaloTower*>::const_iterator hit = towers.begin(); hit != towers.end(); hit++) {
0152         CaloTowerDetId neighborId((*hit)->id());
0153         int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0154                        (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0155         int dPhi = abs(centerId.iphi() - neighborId.iphi());
0156         if (abs(72 - dPhi) < dPhi)
0157           dPhi = 72 - dPhi;
0158         if (dEta <= gridSize && dPhi <= gridSize) {
0159           switch (type) {
0160             case TowerTotal:
0161               energy += (*hit)->energy();
0162               break;
0163             case TowerEcal:
0164               energy += (*hit)->emEnergy();
0165               break;
0166             case TowerHcal:
0167               energy += (*hit)->hadEnergy();
0168               break;
0169             case TowerHO:
0170               energy += (*hit)->outerEnergy();
0171               break;
0172             default:
0173               throw cms::Exception("FatalError") << "Unknown calo tower energy type: " << type;
0174           }
0175         }
0176       }
0177     } break;
0178     case EcalRecHits: {
0179       if (id.det() != DetId::Ecal || (id.subdetId() != EcalBarrel && id.subdetId() != EcalEndcap)) {
0180         throw cms::Exception("FatalError") << "Wrong DetId. Expected EcalBarrel or EcalEndcap, but found:\n"
0181                                            << DetIdInfo::info(id, nullptr) << "\n";
0182       }
0183       // Since the ECAL granularity is small and the gap between EE and EB is significant,
0184       // energy is computed only within the system that contains the central element
0185       if (id.subdetId() == EcalBarrel) {
0186         EBDetId centerId(id);
0187         for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
0188              hit++) {
0189           if ((*hit)->id().subdetId() != EcalBarrel)
0190             continue;
0191           EBDetId neighborId((*hit)->id());
0192           int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0193                          (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0194           int dPhi = abs(centerId.iphi() - neighborId.iphi());
0195           if (abs(360 - dPhi) < dPhi)
0196             dPhi = 360 - dPhi;
0197           if (dEta <= gridSize && dPhi <= gridSize) {
0198             energy += (*hit)->energy();
0199           }
0200         }
0201       } else {
0202         // Endcap
0203         EEDetId centerId(id);
0204         for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
0205              hit++) {
0206           if ((*hit)->id().subdetId() != EcalEndcap)
0207             continue;
0208           EEDetId neighborId((*hit)->id());
0209           if (centerId.zside() == neighborId.zside() && abs(centerId.ix() - neighborId.ix()) <= gridSize &&
0210               abs(centerId.iy() - neighborId.iy()) <= gridSize) {
0211             energy += (*hit)->energy();
0212           }
0213         }
0214       }
0215     } break;
0216     case HcalRecHits: {
0217       if (id.det() != DetId::Hcal || (id.subdetId() != HcalBarrel && id.subdetId() != HcalEndcap)) {
0218         throw cms::Exception("FatalError") << "Wrong DetId. Expected HE or HB, but found:\n"
0219                                            << DetIdInfo::info(id, nullptr) << "\n";
0220       }
0221       HcalDetId centerId(id);
0222       for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++) {
0223         HcalDetId neighborId((*hit)->id());
0224         int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0225                        (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0226         int dPhi = abs(centerId.iphi() - neighborId.iphi());
0227         if (abs(72 - dPhi) < dPhi)
0228           dPhi = 72 - dPhi;
0229         if (dEta <= gridSize && dPhi <= gridSize) {
0230           energy += (*hit)->energy();
0231         }
0232       }
0233     } break;
0234     case HORecHits: {
0235       if (id.det() != DetId::Hcal || (id.subdetId() != HcalOuter)) {
0236         throw cms::Exception("FatalError") << "Wrong DetId. Expected HO, but found:\n"
0237                                            << DetIdInfo::info(id, nullptr) << "\n";
0238       }
0239       HcalDetId centerId(id);
0240       for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++) {
0241         HcalDetId neighborId((*hit)->id());
0242         int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0243                        (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0244         int dPhi = abs(centerId.iphi() - neighborId.iphi());
0245         if (abs(72 - dPhi) < dPhi)
0246           dPhi = 72 - dPhi;
0247         if (dEta <= gridSize && dPhi <= gridSize) {
0248           energy += (*hit)->energy();
0249         }
0250       }
0251     } break;
0252     default:
0253       throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
0254   }
0255   return energy;
0256 }
0257 
0258 double TrackDetMatchInfo::nXnEnergy(EnergyType type, int gridSize) {
0259   switch (type) {
0260     case TowerTotal:
0261     case TowerHcal:
0262     case TowerEcal:
0263     case TowerHO:
0264       if (crossedTowerIds.empty())
0265         return 0;
0266       return nXnEnergy(crossedTowerIds.front(), type, gridSize);
0267       break;
0268     case EcalRecHits:
0269       if (crossedEcalIds.empty())
0270         return 0;
0271       return nXnEnergy(crossedEcalIds.front(), type, gridSize);
0272       break;
0273     case HcalRecHits:
0274       if (crossedHcalIds.empty())
0275         return 0;
0276       return nXnEnergy(crossedHcalIds.front(), type, gridSize);
0277       break;
0278     case HORecHits:
0279       if (crossedHOIds.empty())
0280         return 0;
0281       return nXnEnergy(crossedHOIds.front(), type, gridSize);
0282       break;
0283     default:
0284       throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
0285   }
0286   return -999;
0287 }
0288 
0289 TrackDetMatchInfo::TrackDetMatchInfo()
0290     : trkGlobPosAtEcal(0, 0, 0),
0291       trkGlobPosAtHcal(0, 0, 0),
0292       trkGlobPosAtHO(0, 0, 0),
0293       trkMomAtEcal(0, 0, 0),
0294       trkMomAtHcal(0, 0, 0),
0295       trkMomAtHO(0, 0, 0),
0296       isGoodEcal(false),
0297       isGoodHcal(false),
0298       isGoodCalo(false),
0299       isGoodHO(false),
0300       isGoodMuon(false),
0301       simTrack(nullptr),
0302       ecalTrueEnergy(-999),
0303       hcalTrueEnergy(-999) {}
0304 
0305 DetId TrackDetMatchInfo::findMaxDeposition(EnergyType type) {
0306   DetId id;
0307   float maxEnergy = -9999;
0308   switch (type) {
0309     case EcalRecHits: {
0310       for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end(); hit++)
0311         if ((*hit)->energy() > maxEnergy) {
0312           maxEnergy = (*hit)->energy();
0313           id = (*hit)->detid();
0314         }
0315     } break;
0316     case HcalRecHits: {
0317       for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++)
0318         if ((*hit)->energy() > maxEnergy) {
0319           maxEnergy = (*hit)->energy();
0320           id = (*hit)->detid();
0321         }
0322     } break;
0323     case HORecHits: {
0324       for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++)
0325         if ((*hit)->energy() > maxEnergy) {
0326           maxEnergy = (*hit)->energy();
0327           id = (*hit)->detid();
0328         }
0329     } break;
0330     case TowerTotal:
0331     case TowerEcal:
0332     case TowerHcal:
0333     case TowerHO: {
0334       for (std::vector<const CaloTower*>::const_iterator hit = towers.begin(); hit != towers.end(); hit++) {
0335         double energy = 0;
0336         switch (type) {
0337           case TowerTotal:
0338             energy = (*hit)->energy();
0339             break;
0340           case TowerEcal:
0341             energy = (*hit)->emEnergy();
0342             break;
0343           case TowerHcal:
0344             energy = (*hit)->hadEnergy();
0345             break;
0346           case TowerHO:
0347             energy = (*hit)->energy();
0348             break;
0349           default:
0350             throw cms::Exception("FatalError") << "Unknown calo tower energy type: " << type;
0351         }
0352         if (energy > maxEnergy) {
0353           maxEnergy = energy;
0354           id = (*hit)->id();
0355         }
0356       }
0357     } break;
0358     default:
0359       throw cms::Exception("FatalError")
0360           << "Maximal energy deposition: unkown or not implemented energy type requested, type:" << type;
0361   }
0362   return id;
0363 }
0364 
0365 DetId TrackDetMatchInfo::findMaxDeposition(const DetId& id, EnergyType type, int gridSize) {
0366   double energy_max(0);
0367   DetId id_max;
0368   if (id.rawId() == 0)
0369     return id_max;
0370   switch (type) {
0371     case TowerTotal:
0372     case TowerHcal:
0373     case TowerEcal:
0374     case TowerHO: {
0375       if (id.det() != DetId::Calo) {
0376         throw cms::Exception("FatalError") << "Wrong DetId. Expected CaloTower, but found:\n"
0377                                            << DetIdInfo::info(id, nullptr) << "\n";
0378       }
0379       CaloTowerDetId centerId(id);
0380       for (std::vector<const CaloTower*>::const_iterator hit = towers.begin(); hit != towers.end(); hit++) {
0381         CaloTowerDetId neighborId((*hit)->id());
0382         int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0383                        (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0384         int dPhi = abs(centerId.iphi() - neighborId.iphi());
0385         if (abs(72 - dPhi) < dPhi)
0386           dPhi = 72 - dPhi;
0387         if (dEta <= gridSize && dPhi <= gridSize) {
0388           switch (type) {
0389             case TowerTotal:
0390               if (energy_max < (*hit)->energy()) {
0391                 energy_max = (*hit)->energy();
0392                 id_max = (*hit)->id();
0393               }
0394               break;
0395             case TowerEcal:
0396               if (energy_max < (*hit)->emEnergy()) {
0397                 energy_max = (*hit)->emEnergy();
0398                 id_max = (*hit)->id();
0399               }
0400               break;
0401             case TowerHcal:
0402               if (energy_max < (*hit)->hadEnergy()) {
0403                 energy_max = (*hit)->hadEnergy();
0404                 id_max = (*hit)->id();
0405               }
0406               break;
0407             case TowerHO:
0408               if (energy_max < (*hit)->outerEnergy()) {
0409                 energy_max = (*hit)->outerEnergy();
0410                 id_max = (*hit)->id();
0411               }
0412               break;
0413             default:
0414               throw cms::Exception("FatalError") << "Unknown calo tower energy type: " << type;
0415           }
0416         }
0417       }
0418     } break;
0419     case EcalRecHits: {
0420       if (id.det() != DetId::Ecal || (id.subdetId() != EcalBarrel && id.subdetId() != EcalEndcap)) {
0421         throw cms::Exception("FatalError") << "Wrong DetId. Expected EcalBarrel or EcalEndcap, but found:\n"
0422                                            << DetIdInfo::info(id, nullptr) << "\n";
0423       }
0424       // Since the ECAL granularity is small and the gap between EE and EB is significant,
0425       // energy is computed only within the system that contains the central element
0426       if (id.subdetId() == EcalBarrel) {
0427         EBDetId centerId(id);
0428         for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
0429              hit++) {
0430           if ((*hit)->id().subdetId() != EcalBarrel)
0431             continue;
0432           EBDetId neighborId((*hit)->id());
0433           int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0434                          (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0435           int dPhi = abs(centerId.iphi() - neighborId.iphi());
0436           if (abs(360 - dPhi) < dPhi)
0437             dPhi = 360 - dPhi;
0438           if (dEta <= gridSize && dPhi <= gridSize) {
0439             if (energy_max < (*hit)->energy()) {
0440               energy_max = (*hit)->energy();
0441               id_max = (*hit)->id();
0442             }
0443           }
0444         }
0445       } else {
0446         // Endcap
0447         EEDetId centerId(id);
0448         for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
0449              hit++) {
0450           if ((*hit)->id().subdetId() != EcalEndcap)
0451             continue;
0452           EEDetId neighborId((*hit)->id());
0453           if (centerId.zside() == neighborId.zside() && abs(centerId.ix() - neighborId.ix()) <= gridSize &&
0454               abs(centerId.iy() - neighborId.iy()) <= gridSize) {
0455             if (energy_max < (*hit)->energy()) {
0456               energy_max = (*hit)->energy();
0457               id_max = (*hit)->id();
0458             }
0459           }
0460         }
0461       }
0462     } break;
0463     case HcalRecHits: {
0464       if (id.det() != DetId::Hcal || (id.subdetId() != HcalBarrel && id.subdetId() != HcalEndcap)) {
0465         throw cms::Exception("FatalError") << "Wrong DetId. Expected HE or HB, but found:\n"
0466                                            << DetIdInfo::info(id, nullptr) << "\n";
0467       }
0468       HcalDetId centerId(id);
0469       for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++) {
0470         HcalDetId neighborId((*hit)->id());
0471         int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0472                        (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0473         int dPhi = abs(centerId.iphi() - neighborId.iphi());
0474         if (abs(72 - dPhi) < dPhi)
0475           dPhi = 72 - dPhi;
0476         if (dEta <= gridSize && dPhi <= gridSize) {
0477           if (energy_max < (*hit)->energy()) {
0478             energy_max = (*hit)->energy();
0479             id_max = (*hit)->id();
0480           }
0481         }
0482       }
0483     } break;
0484     case HORecHits: {
0485       if (id.det() != DetId::Hcal || (id.subdetId() != HcalOuter)) {
0486         throw cms::Exception("FatalError") << "Wrong DetId. Expected HO, but found:\n"
0487                                            << DetIdInfo::info(id, nullptr) << "\n";
0488       }
0489       HcalDetId centerId(id);
0490       for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++) {
0491         HcalDetId neighborId((*hit)->id());
0492         int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
0493                        (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
0494         int dPhi = abs(centerId.iphi() - neighborId.iphi());
0495         if (abs(72 - dPhi) < dPhi)
0496           dPhi = 72 - dPhi;
0497         if (dEta <= gridSize && dPhi <= gridSize) {
0498           if (energy_max < (*hit)->energy()) {
0499             energy_max = (*hit)->energy();
0500             id_max = (*hit)->id();
0501           }
0502         }
0503       }
0504     } break;
0505     default:
0506       throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
0507   }
0508   return id_max;
0509 }
0510 
0511 DetId TrackDetMatchInfo::findMaxDeposition(EnergyType type, int gridSize) {
0512   DetId id_max;
0513   switch (type) {
0514     case TowerTotal:
0515     case TowerHcal:
0516     case TowerEcal:
0517     case TowerHO:
0518       if (crossedTowerIds.empty())
0519         return id_max;
0520       return findMaxDeposition(crossedTowerIds.front(), type, gridSize);
0521       break;
0522     case EcalRecHits:
0523       if (crossedEcalIds.empty())
0524         return id_max;
0525       return findMaxDeposition(crossedEcalIds.front(), type, gridSize);
0526       break;
0527     case HcalRecHits:
0528       if (crossedHcalIds.empty())
0529         return id_max;
0530       return findMaxDeposition(crossedHcalIds.front(), type, gridSize);
0531       break;
0532     case HORecHits:
0533       if (crossedHOIds.empty())
0534         return id_max;
0535       return findMaxDeposition(crossedHOIds.front(), type, gridSize);
0536       break;
0537     default:
0538       throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
0539   }
0540   return id_max;
0541 }
0542 
0543 ////////////////////////////////////////////////////////////////////////
0544 // Obsolete
0545 //
0546 
0547 double TrackDetMatchInfo::ecalConeEnergy() { return coneEnergy(999, EcalRecHits); }
0548 
0549 double TrackDetMatchInfo::hcalConeEnergy() { return coneEnergy(999, HcalRecHits); }
0550 
0551 double TrackDetMatchInfo::hoConeEnergy() { return coneEnergy(999, HcalRecHits); }
0552 
0553 double TrackDetMatchInfo::ecalCrossedEnergy() { return crossedEnergy(EcalRecHits); }
0554 
0555 double TrackDetMatchInfo::hcalCrossedEnergy() { return crossedEnergy(HcalRecHits); }
0556 
0557 double TrackDetMatchInfo::hoCrossedEnergy() { return crossedEnergy(HORecHits); }
0558 
0559 int TrackDetMatchInfo::numberOfSegments() const {
0560   int numSegments = 0;
0561   for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
0562     numSegments += chamber->segments.size();
0563   return numSegments;
0564 }
0565 
0566 int TrackDetMatchInfo::numberOfSegmentsInStation(int station) const {
0567   int numSegments = 0;
0568   for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
0569     if (chamber->station() == station)
0570       numSegments += chamber->segments.size();
0571   return numSegments;
0572 }
0573 
0574 int TrackDetMatchInfo::numberOfSegmentsInStation(int station, int detector) const {
0575   int numSegments = 0;
0576   for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
0577     if (chamber->station() == station && chamber->detector() == detector)
0578       numSegments += chamber->segments.size();
0579   return numSegments;
0580 }
0581 
0582 int TrackDetMatchInfo::numberOfSegmentsInDetector(int detector) const {
0583   int numSegments = 0;
0584   for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
0585     if (chamber->detector() == detector)
0586       numSegments += chamber->segments.size();
0587   return numSegments;
0588 }