TrackingTools/TrackAssociator/src/TrackDetMatchInfo.cc - master - githubmirror/cmssw

TrackingTools/TrackAssociator/src/TrackDetMatchInfo.cc
588 строк23 KB
Malik Shahzad Muzaffar
added missing TowerHO in the case statement
30 окт 2019, 13:38
30 окт 2019, 13:38f0b6f49
#include <map>
#include "TrackingTools/TrackAssociator/interface/TrackDetMatchInfo.h"
#include "DetIdInfo.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "Math/VectorUtil.h"
#include <algorithm>

///////////////////////////

std::string TrackDetMatchInfo::dumpGeometry(const DetId& id) {
  if (!caloGeometry || !caloGeometry->getSubdetectorGeometry(id) ||
      !caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)) {
    throw cms::Exception("FatalError") << "Failed to access geometry for DetId: " << id.rawId();
  }
  std::ostringstream oss;

  const CaloCellGeometry::CornersVec& points = caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)->getCorners();
  for (CaloCellGeometry::CornersVec::const_iterator point = points.begin(); point != points.end(); ++point)
    oss << "(" << point->z() << ", " << point->perp() << ", " << point->eta() << ", " << point->phi() << "), \t";
  return oss.str();
}

GlobalPoint TrackDetMatchInfo::getPosition(const DetId& id) {
  // this part might be slow
  if (!caloGeometry || !caloGeometry->getSubdetectorGeometry(id) ||
      !caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)) {
    throw cms::Exception("FatalError") << "Failed to access geometry for DetId: " << id.rawId();
    return GlobalPoint(0, 0, 0);
  }
  return caloGeometry->getSubdetectorGeometry(id)->getGeometry(id)->getPosition();
}

double TrackDetMatchInfo::crossedEnergy(EnergyType type) {
  double energy(0);
  switch (type) {
    case EcalRecHits: {
      for (std::vector<const EcalRecHit*>::const_iterator hit = crossedEcalRecHits.begin();
           hit != crossedEcalRecHits.end();
           hit++)
        energy += (*hit)->energy();
    } break;
    case HcalRecHits: {
      for (std::vector<const HBHERecHit*>::const_iterator hit = crossedHcalRecHits.begin();
           hit != crossedHcalRecHits.end();
           hit++)
        energy += (*hit)->energy();
    } break;
    case HORecHits: {
      for (std::vector<const HORecHit*>::const_iterator hit = crossedHORecHits.begin(); hit != crossedHORecHits.end();
           hit++)
        energy += (*hit)->energy();
    } break;
    case TowerTotal: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        energy += (*hit)->energy();
    } break;
    case TowerEcal: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        energy += (*hit)->emEnergy();
    } break;
    case TowerHcal: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        energy += (*hit)->hadEnergy();
    } break;
    case TowerHO: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        energy += (*hit)->outerEnergy();
    } break;
    default:
      throw cms::Exception("FatalError") << "Unknown calo energy type: " << type;
  }
  return energy;
}

bool TrackDetMatchInfo::insideCone(const DetId& id, const double dR) {
  GlobalPoint idPosition = getPosition(id);
  if (idPosition.mag() < 0.01)
    return false;

  math::XYZVector idPositionRoot(idPosition.x(), idPosition.y(), idPosition.z());
  math::XYZVector trackP3(stateAtIP.momentum().x(), stateAtIP.momentum().y(), stateAtIP.momentum().z());
  return ROOT::Math::VectorUtil::DeltaR(trackP3, idPositionRoot) < dR;
}

double TrackDetMatchInfo::coneEnergy(double dR, EnergyType type) {
  double energy(0);
  switch (type) {
    case EcalRecHits: {
      for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end(); hit++)
        if (insideCone((*hit)->detid(), dR))
          energy += (*hit)->energy();
    } break;
    case HcalRecHits: {
      for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++)
        if (insideCone((*hit)->detid(), dR))
          energy += (*hit)->energy();
    } break;
    case HORecHits: {
      for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++)
        if (insideCone((*hit)->detid(), dR))
          energy += (*hit)->energy();
    } break;
    case TowerTotal: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        if (insideCone((*hit)->id(), dR))
          energy += (*hit)->energy();
    } break;
    case TowerEcal: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        if (insideCone((*hit)->id(), dR))
          energy += (*hit)->emEnergy();
    } break;
    case TowerHcal: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        if (insideCone((*hit)->id(), dR))
          energy += (*hit)->hadEnergy();
    } break;
    case TowerHO: {
      for (std::vector<const CaloTower*>::const_iterator hit = crossedTowers.begin(); hit != crossedTowers.end(); hit++)
        if (insideCone((*hit)->id(), dR))
          energy += (*hit)->outerEnergy();
    } break;
    default:
      throw cms::Exception("FatalError") << "Unknown calo energy type: " << type;
  }
  return energy;
}

//////////////////////////////////////////////////

double TrackDetMatchInfo::nXnEnergy(const DetId& id, EnergyType type, int gridSize) {
  double energy(0);
  if (id.rawId() == 0)
    return 0.;
  switch (type) {
    case TowerTotal:
    case TowerHcal:
    case TowerEcal:
    case TowerHO: {
      if (id.det() != DetId::Calo) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected CaloTower, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      CaloTowerDetId centerId(id);
      for (std::vector<const CaloTower*>::const_iterator hit = towers.begin(); hit != towers.end(); hit++) {
        CaloTowerDetId neighborId((*hit)->id());
        int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                       (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
        int dPhi = abs(centerId.iphi() - neighborId.iphi());
        if (abs(72 - dPhi) < dPhi)
          dPhi = 72 - dPhi;
        if (dEta <= gridSize && dPhi <= gridSize) {
          switch (type) {
            case TowerTotal:
              energy += (*hit)->energy();
              break;
            case TowerEcal:
              energy += (*hit)->emEnergy();
              break;
            case TowerHcal:
              energy += (*hit)->hadEnergy();
              break;
            case TowerHO:
              energy += (*hit)->outerEnergy();
              break;
            default:
              throw cms::Exception("FatalError") << "Unknown calo tower energy type: " << type;
          }
        }
      }
    } break;
    case EcalRecHits: {
      if (id.det() != DetId::Ecal || (id.subdetId() != EcalBarrel && id.subdetId() != EcalEndcap)) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected EcalBarrel or EcalEndcap, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      // Since the ECAL granularity is small and the gap between EE and EB is significant,
      // energy is computed only within the system that contains the central element
      if (id.subdetId() == EcalBarrel) {
        EBDetId centerId(id);
        for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
             hit++) {
          if ((*hit)->id().subdetId() != EcalBarrel)
            continue;
          EBDetId neighborId((*hit)->id());
          int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                         (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
          int dPhi = abs(centerId.iphi() - neighborId.iphi());
          if (abs(360 - dPhi) < dPhi)
            dPhi = 360 - dPhi;
          if (dEta <= gridSize && dPhi <= gridSize) {
            energy += (*hit)->energy();
          }
        }
      } else {
        // Endcap
        EEDetId centerId(id);
        for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
             hit++) {
          if ((*hit)->id().subdetId() != EcalEndcap)
            continue;
          EEDetId neighborId((*hit)->id());
          if (centerId.zside() == neighborId.zside() && abs(centerId.ix() - neighborId.ix()) <= gridSize &&
              abs(centerId.iy() - neighborId.iy()) <= gridSize) {
            energy += (*hit)->energy();
          }
        }
      }
    } break;
    case HcalRecHits: {
      if (id.det() != DetId::Hcal || (id.subdetId() != HcalBarrel && id.subdetId() != HcalEndcap)) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected HE or HB, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      HcalDetId centerId(id);
      for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++) {
        HcalDetId neighborId((*hit)->id());
        int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                       (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
        int dPhi = abs(centerId.iphi() - neighborId.iphi());
        if (abs(72 - dPhi) < dPhi)
          dPhi = 72 - dPhi;
        if (dEta <= gridSize && dPhi <= gridSize) {
          energy += (*hit)->energy();
        }
      }
    } break;
    case HORecHits: {
      if (id.det() != DetId::Hcal || (id.subdetId() != HcalOuter)) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected HO, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      HcalDetId centerId(id);
      for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++) {
        HcalDetId neighborId((*hit)->id());
        int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                       (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
        int dPhi = abs(centerId.iphi() - neighborId.iphi());
        if (abs(72 - dPhi) < dPhi)
          dPhi = 72 - dPhi;
        if (dEta <= gridSize && dPhi <= gridSize) {
          energy += (*hit)->energy();
        }
      }
    } break;
    default:
      throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
  }
  return energy;
}

double TrackDetMatchInfo::nXnEnergy(EnergyType type, int gridSize) {
  switch (type) {
    case TowerTotal:
    case TowerHcal:
    case TowerEcal:
    case TowerHO:
      if (crossedTowerIds.empty())
        return 0;
      return nXnEnergy(crossedTowerIds.front(), type, gridSize);
      break;
    case EcalRecHits:
      if (crossedEcalIds.empty())
        return 0;
      return nXnEnergy(crossedEcalIds.front(), type, gridSize);
      break;
    case HcalRecHits:
      if (crossedHcalIds.empty())
        return 0;
      return nXnEnergy(crossedHcalIds.front(), type, gridSize);
      break;
    case HORecHits:
      if (crossedHOIds.empty())
        return 0;
      return nXnEnergy(crossedHOIds.front(), type, gridSize);
      break;
    default:
      throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
  }
  return -999;
}

TrackDetMatchInfo::TrackDetMatchInfo()
    : trkGlobPosAtEcal(0, 0, 0),
      trkGlobPosAtHcal(0, 0, 0),
      trkGlobPosAtHO(0, 0, 0),
      trkMomAtEcal(0, 0, 0),
      trkMomAtHcal(0, 0, 0),
      trkMomAtHO(0, 0, 0),
      isGoodEcal(false),
      isGoodHcal(false),
      isGoodCalo(false),
      isGoodHO(false),
      isGoodMuon(false),
      simTrack(nullptr),
      ecalTrueEnergy(-999),
      hcalTrueEnergy(-999) {}

DetId TrackDetMatchInfo::findMaxDeposition(EnergyType type) {
  DetId id;
  float maxEnergy = -9999;
  switch (type) {
    case EcalRecHits: {
      for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end(); hit++)
        if ((*hit)->energy() > maxEnergy) {
          maxEnergy = (*hit)->energy();
          id = (*hit)->detid();
        }
    } break;
    case HcalRecHits: {
      for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++)
        if ((*hit)->energy() > maxEnergy) {
          maxEnergy = (*hit)->energy();
          id = (*hit)->detid();
        }
    } break;
    case HORecHits: {
      for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++)
        if ((*hit)->energy() > maxEnergy) {
          maxEnergy = (*hit)->energy();
          id = (*hit)->detid();
        }
    } break;
    case TowerTotal:
    case TowerEcal:
    case TowerHcal:
    case TowerHO: {
      for (std::vector<const CaloTower*>::const_iterator hit = towers.begin(); hit != towers.end(); hit++) {
        double energy = 0;
        switch (type) {
          case TowerTotal:
            energy = (*hit)->energy();
            break;
          case TowerEcal:
            energy = (*hit)->emEnergy();
            break;
          case TowerHcal:
            energy = (*hit)->hadEnergy();
            break;
          case TowerHO:
            energy = (*hit)->energy();
            break;
          default:
            throw cms::Exception("FatalError") << "Unknown calo tower energy type: " << type;
        }
        if (energy > maxEnergy) {
          maxEnergy = energy;
          id = (*hit)->id();
        }
      }
    } break;
    default:
      throw cms::Exception("FatalError")
          << "Maximal energy deposition: unkown or not implemented energy type requested, type:" << type;
  }
  return id;
}

DetId TrackDetMatchInfo::findMaxDeposition(const DetId& id, EnergyType type, int gridSize) {
  double energy_max(0);
  DetId id_max;
  if (id.rawId() == 0)
    return id_max;
  switch (type) {
    case TowerTotal:
    case TowerHcal:
    case TowerEcal:
    case TowerHO: {
      if (id.det() != DetId::Calo) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected CaloTower, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      CaloTowerDetId centerId(id);
      for (std::vector<const CaloTower*>::const_iterator hit = towers.begin(); hit != towers.end(); hit++) {
        CaloTowerDetId neighborId((*hit)->id());
        int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                       (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
        int dPhi = abs(centerId.iphi() - neighborId.iphi());
        if (abs(72 - dPhi) < dPhi)
          dPhi = 72 - dPhi;
        if (dEta <= gridSize && dPhi <= gridSize) {
          switch (type) {
            case TowerTotal:
              if (energy_max < (*hit)->energy()) {
                energy_max = (*hit)->energy();
                id_max = (*hit)->id();
              }
              break;
            case TowerEcal:
              if (energy_max < (*hit)->emEnergy()) {
                energy_max = (*hit)->emEnergy();
                id_max = (*hit)->id();
              }
              break;
            case TowerHcal:
              if (energy_max < (*hit)->hadEnergy()) {
                energy_max = (*hit)->hadEnergy();
                id_max = (*hit)->id();
              }
              break;
            case TowerHO:
              if (energy_max < (*hit)->outerEnergy()) {
                energy_max = (*hit)->outerEnergy();
                id_max = (*hit)->id();
              }
              break;
            default:
              throw cms::Exception("FatalError") << "Unknown calo tower energy type: " << type;
          }
        }
      }
    } break;
    case EcalRecHits: {
      if (id.det() != DetId::Ecal || (id.subdetId() != EcalBarrel && id.subdetId() != EcalEndcap)) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected EcalBarrel or EcalEndcap, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      // Since the ECAL granularity is small and the gap between EE and EB is significant,
      // energy is computed only within the system that contains the central element
      if (id.subdetId() == EcalBarrel) {
        EBDetId centerId(id);
        for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
             hit++) {
          if ((*hit)->id().subdetId() != EcalBarrel)
            continue;
          EBDetId neighborId((*hit)->id());
          int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                         (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
          int dPhi = abs(centerId.iphi() - neighborId.iphi());
          if (abs(360 - dPhi) < dPhi)
            dPhi = 360 - dPhi;
          if (dEta <= gridSize && dPhi <= gridSize) {
            if (energy_max < (*hit)->energy()) {
              energy_max = (*hit)->energy();
              id_max = (*hit)->id();
            }
          }
        }
      } else {
        // Endcap
        EEDetId centerId(id);
        for (std::vector<const EcalRecHit*>::const_iterator hit = ecalRecHits.begin(); hit != ecalRecHits.end();
             hit++) {
          if ((*hit)->id().subdetId() != EcalEndcap)
            continue;
          EEDetId neighborId((*hit)->id());
          if (centerId.zside() == neighborId.zside() && abs(centerId.ix() - neighborId.ix()) <= gridSize &&
              abs(centerId.iy() - neighborId.iy()) <= gridSize) {
            if (energy_max < (*hit)->energy()) {
              energy_max = (*hit)->energy();
              id_max = (*hit)->id();
            }
          }
        }
      }
    } break;
    case HcalRecHits: {
      if (id.det() != DetId::Hcal || (id.subdetId() != HcalBarrel && id.subdetId() != HcalEndcap)) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected HE or HB, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      HcalDetId centerId(id);
      for (std::vector<const HBHERecHit*>::const_iterator hit = hcalRecHits.begin(); hit != hcalRecHits.end(); hit++) {
        HcalDetId neighborId((*hit)->id());
        int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                       (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
        int dPhi = abs(centerId.iphi() - neighborId.iphi());
        if (abs(72 - dPhi) < dPhi)
          dPhi = 72 - dPhi;
        if (dEta <= gridSize && dPhi <= gridSize) {
          if (energy_max < (*hit)->energy()) {
            energy_max = (*hit)->energy();
            id_max = (*hit)->id();
          }
        }
      }
    } break;
    case HORecHits: {
      if (id.det() != DetId::Hcal || (id.subdetId() != HcalOuter)) {
        throw cms::Exception("FatalError") << "Wrong DetId. Expected HO, but found:\n"
                                           << DetIdInfo::info(id, nullptr) << "\n";
      }
      HcalDetId centerId(id);
      for (std::vector<const HORecHit*>::const_iterator hit = hoRecHits.begin(); hit != hoRecHits.end(); hit++) {
        HcalDetId neighborId((*hit)->id());
        int dEta = abs((centerId.ieta() < 0 ? centerId.ieta() + 1 : centerId.ieta()) -
                       (neighborId.ieta() < 0 ? neighborId.ieta() + 1 : neighborId.ieta()));
        int dPhi = abs(centerId.iphi() - neighborId.iphi());
        if (abs(72 - dPhi) < dPhi)
          dPhi = 72 - dPhi;
        if (dEta <= gridSize && dPhi <= gridSize) {
          if (energy_max < (*hit)->energy()) {
            energy_max = (*hit)->energy();
            id_max = (*hit)->id();
          }
        }
      }
    } break;
    default:
      throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
  }
  return id_max;
}

DetId TrackDetMatchInfo::findMaxDeposition(EnergyType type, int gridSize) {
  DetId id_max;
  switch (type) {
    case TowerTotal:
    case TowerHcal:
    case TowerEcal:
    case TowerHO:
      if (crossedTowerIds.empty())
        return id_max;
      return findMaxDeposition(crossedTowerIds.front(), type, gridSize);
      break;
    case EcalRecHits:
      if (crossedEcalIds.empty())
        return id_max;
      return findMaxDeposition(crossedEcalIds.front(), type, gridSize);
      break;
    case HcalRecHits:
      if (crossedHcalIds.empty())
        return id_max;
      return findMaxDeposition(crossedHcalIds.front(), type, gridSize);
      break;
    case HORecHits:
      if (crossedHOIds.empty())
        return id_max;
      return findMaxDeposition(crossedHOIds.front(), type, gridSize);
      break;
    default:
      throw cms::Exception("FatalError") << "Unkown or not implemented energy type requested, type:" << type;
  }
  return id_max;
}

////////////////////////////////////////////////////////////////////////
// Obsolete
//

double TrackDetMatchInfo::ecalConeEnergy() { return coneEnergy(999, EcalRecHits); }

double TrackDetMatchInfo::hcalConeEnergy() { return coneEnergy(999, HcalRecHits); }

double TrackDetMatchInfo::hoConeEnergy() { return coneEnergy(999, HcalRecHits); }

double TrackDetMatchInfo::ecalCrossedEnergy() { return crossedEnergy(EcalRecHits); }

double TrackDetMatchInfo::hcalCrossedEnergy() { return crossedEnergy(HcalRecHits); }

double TrackDetMatchInfo::hoCrossedEnergy() { return crossedEnergy(HORecHits); }

int TrackDetMatchInfo::numberOfSegments() const {
  int numSegments = 0;
  for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
    numSegments += chamber->segments.size();
  return numSegments;
}

int TrackDetMatchInfo::numberOfSegmentsInStation(int station) const {
  int numSegments = 0;
  for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
    if (chamber->station() == station)
      numSegments += chamber->segments.size();
  return numSegments;
}

int TrackDetMatchInfo::numberOfSegmentsInStation(int station, int detector) const {
  int numSegments = 0;
  for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
    if (chamber->station() == station && chamber->detector() == detector)
      numSegments += chamber->segments.size();
  return numSegments;
}

int TrackDetMatchInfo::numberOfSegmentsInDetector(int detector) const {
  int numSegments = 0;
  for (std::vector<TAMuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); chamber++)
    if (chamber->detector() == detector)
      numSegments += chamber->segments.size();
  return numSegments;
}