MuonDetLayerGeometry.cc

Go to the documentation of this file.

#include <RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h>

#include <FWCore/Utilities/interface/Exception.h>
#include <TrackingTools/DetLayers/interface/DetLayer.h>
#include <DataFormats/MuonDetId/interface/CSCDetId.h>
#include <DataFormats/MuonDetId/interface/DTChamberId.h>
#include <DataFormats/MuonDetId/interface/RPCDetId.h>
#include <DataFormats/MuonDetId/interface/GEMDetId.h>
#include <DataFormats/MuonDetId/interface/ME0DetId.h>

#include <Utilities/General/interface/precomputed_value_sort.h>
#include <DataFormats/GeometrySurface/interface/GeometricSorting.h>

#include <algorithm>

using namespace std;
using namespace geomsort;

MuonDetLayerGeometry::MuonDetLayerGeometry() {}

MuonDetLayerGeometry::~MuonDetLayerGeometry() {
  for (vector<const DetLayer*>::const_iterator it = allDetLayers.begin(); it != allDetLayers.end(); ++it) {
    delete *it;
  }
}

void MuonDetLayerGeometry::addCSCLayers(const pair<vector<DetLayer*>, vector<DetLayer*> >& csclayers) {
  for (auto const it : csclayers.first) {
    cscLayers_fw.push_back(it);
    allForward.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
  }

  for (auto const it : csclayers.second) {
    cscLayers_bk.push_back(it);
    allBackward.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
  }
}

void MuonDetLayerGeometry::addGEMLayers(const pair<vector<DetLayer*>, vector<DetLayer*> >& gemlayers) {
  LogDebug("Muon|RecoMuon|MuonDetLayerGeometry") << "Adding GEMlayers " << std::endl;
  for (auto const it : gemlayers.first) {
    gemLayers_fw.push_back(it);
    allForward.push_back(it);
    detLayersMap[makeDetLayerId(it)] = it;
    LogDebug("Muon|RecoMuon|MuonDetLayerGeometry") << "Adding GEMforward " << std::endl;
  }
  for (auto const it : gemlayers.second) {
    gemLayers_bk.push_back(it);
    allBackward.push_back(it);
    detLayersMap[makeDetLayerId(it)] = it;
    LogDebug("Muon|RecoMuon|MuonDetLayerGeometry") << "Adding GEMbackward " << std::endl;
  }
}

void MuonDetLayerGeometry::addME0Layers(const pair<vector<DetLayer*>, vector<DetLayer*> >& me0layers) {
  LogDebug("Muon|RecoMuon|MuonDetLayerGeometry") << "Adding ME0layers " << std::endl;

  for (auto const it : me0layers.first) {
    me0Layers_fw.push_back(it);
    allForward.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
    LogDebug("Muon|RecoMuon|MuonDetLayerGeometry") << "Adding ME0forward " << std::endl;
  }
  for (auto const it : me0layers.second) {
    me0Layers_bk.push_back(it);
    allBackward.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
    LogDebug("Muon|RecoMuon|MuonDetLayerGeometry") << "Adding ME0backward " << std::endl;
  }
}

void MuonDetLayerGeometry::addRPCLayers(const vector<DetLayer*>& barrelLayers,
                                        const pair<vector<DetLayer*>, vector<DetLayer*> >& endcapLayers) {
  for (auto const it : barrelLayers) {
    rpcLayers_barrel.push_back(it);
    allBarrel.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
  }
  for (auto const it : endcapLayers.first) {
    rpcLayers_fw.push_back(it);
    allForward.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
  }

  for (auto const it : endcapLayers.second) {
    rpcLayers_bk.push_back(it);
    allBackward.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
  }
}

void MuonDetLayerGeometry::addDTLayers(const vector<DetLayer*>& dtlayers) {
  for (auto const it : dtlayers) {
    dtLayers.push_back(it);
    allBarrel.push_back(it);

    detLayersMap[makeDetLayerId(it)] = it;
  }
}

DetId MuonDetLayerGeometry::makeDetLayerId(const DetLayer* detLayer) const {
  if (detLayer->subDetector() == GeomDetEnumerators::CSC) {
    CSCDetId id(detLayer->basicComponents().front()->geographicalId().rawId());

    if (id.station() == 1) {
      if (id.ring() == 1 || id.ring() == 4)
        return CSCDetId(id.endcap(), 1, 1, 0, 0);
      else if (id.ring() == 2 || id.ring() == 3)
        return CSCDetId(id.endcap(), 1, 2, 0, 0);
      else
        throw cms::Exception("InvalidCSCRing") << " Invalid CSC Ring: " << id.ring() << endl;
    } else
      return CSCDetId(id.endcap(), id.station(), 0, 0, 0);

  } else if (detLayer->subDetector() == GeomDetEnumerators::DT) {
    DTChamberId id(detLayer->basicComponents().front()->geographicalId().rawId());
    return DTChamberId(0, id.station(), 0);
  } else if (detLayer->subDetector() == GeomDetEnumerators::RPCBarrel ||
             detLayer->subDetector() == GeomDetEnumerators::RPCEndcap) {
    RPCDetId id(detLayer->basicComponents().front()->geographicalId().rawId());
    return RPCDetId(id.region(), 0, id.station(), 0, id.layer(), 0, 0);
  } else if (detLayer->subDetector() == GeomDetEnumerators::GEM) {
    GEMDetId id(detLayer->basicComponents().front()->geographicalId().rawId());
    return GEMDetId(id.region(), 1, id.station(), id.layer(), 0, 0);
  } else if (detLayer->subDetector() == GeomDetEnumerators::ME0) {
    ME0DetId id(detLayer->basicComponents().front()->geographicalId().rawId());
    return ME0DetId(id.region(), id.layer(), 0, 0);
  } else
    throw cms::Exception("InvalidModuleIdentification");  // << detLayer->module();
}

const vector<const DetLayer*>& MuonDetLayerGeometry::allDTLayers() const { return dtLayers; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allCSCLayers() const { return cscLayers_all; }

const vector<const DetLayer*>& MuonDetLayerGeometry::forwardCSCLayers() const { return cscLayers_fw; }

const vector<const DetLayer*>& MuonDetLayerGeometry::backwardCSCLayers() const { return cscLayers_bk; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allGEMLayers() const { return gemLayers_all; }

const vector<const DetLayer*>& MuonDetLayerGeometry::forwardGEMLayers() const { return gemLayers_fw; }

const vector<const DetLayer*>& MuonDetLayerGeometry::backwardGEMLayers() const { return gemLayers_bk; }



const vector<const DetLayer*>& MuonDetLayerGeometry::allME0Layers() const { return me0Layers_all; }

const vector<const DetLayer*>& MuonDetLayerGeometry::forwardME0Layers() const { return me0Layers_fw; }

const vector<const DetLayer*>& MuonDetLayerGeometry::backwardME0Layers() const { return me0Layers_bk; }


const vector<const DetLayer*>& MuonDetLayerGeometry::allRPCLayers() const { return rpcLayers_all; }

const vector<const DetLayer*>& MuonDetLayerGeometry::barrelRPCLayers() const { return rpcLayers_barrel; }

const vector<const DetLayer*>& MuonDetLayerGeometry::endcapRPCLayers() const { return rpcLayers_endcap; }

const vector<const DetLayer*>& MuonDetLayerGeometry::forwardRPCLayers() const { return rpcLayers_fw; }

const vector<const DetLayer*>& MuonDetLayerGeometry::backwardRPCLayers() const { return rpcLayers_bk; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allLayers() const { return allDetLayers; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allBarrelLayers() const { return allBarrel; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allEndcapLayers() const { return allEndcap; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allForwardLayers() const { return allForward; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allBackwardLayers() const { return allBackward; }

//
const vector<const DetLayer*>& MuonDetLayerGeometry::allEndcapCscGemLayers() const { return allEndcapCscGem; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allCscGemForwardLayers() const { return allCscGemForward; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allCscGemBackwardLayers() const { return allCscGemBackward; }


const vector<const DetLayer*>& MuonDetLayerGeometry::allEndcapCscME0Layers() const { return allEndcapCscME0; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allCscME0ForwardLayers() const { return allCscME0Forward; }

const vector<const DetLayer*>& MuonDetLayerGeometry::allCscME0BackwardLayers() const { return allCscME0Backward; }


const DetLayer* MuonDetLayerGeometry::idToLayer(const DetId& detId) const {
  DetId id;

  if (detId.subdetId() == MuonSubdetId::CSC) {
    CSCDetId cscId(detId.rawId());

    if (cscId.station() == 1) {
      if (cscId.ring() == 1 || cscId.ring() == 4)
        id = CSCDetId(cscId.endcap(), 1, 1, 0, 0);
      else if (cscId.ring() == 2 || cscId.ring() == 3)
        id = CSCDetId(cscId.endcap(), 1, 2, 0, 0);
      else
        throw cms::Exception("InvalidCSCRing") << " Invalid CSC Ring: " << cscId.ring() << endl;
    } else
      id = CSCDetId(cscId.endcap(), cscId.station(), 0, 0, 0);
  }

  else if (detId.subdetId() == MuonSubdetId::DT) {
    DTChamberId dtId(detId.rawId());
    id = DTChamberId(0, dtId.station(), 0);
  } else if (detId.subdetId() == MuonSubdetId::RPC) {
    RPCDetId rpcId(detId.rawId());
    id = RPCDetId(rpcId.region(), 0, rpcId.station(), 0, rpcId.layer(), 0, 0);
  } else if (detId.subdetId() == MuonSubdetId::GEM) {
    GEMDetId gemId(detId.rawId());
    id = GEMDetId(gemId.region(), 1, gemId.station(), gemId.layer(), 0, 0);
  } else if (detId.subdetId() == MuonSubdetId::ME0) {
    ME0DetId me0Id(detId.rawId());
    id = ME0DetId(me0Id.region(), me0Id.layer(), 0, 0);
    LogDebug("Muon|RecoMuon|MuonDetLayerGeometry")
        << " Found an ME0DetId:  " << me0Id.rawId() << ",id: " << id.rawId() << std::endl;

  } else
    throw cms::Exception("InvalidSubdetId") << detId.subdetId();

  std::map<DetId, const DetLayer*>::const_iterator layer = detLayersMap.find(id);
  if (layer == detLayersMap.end())
    return nullptr;
  return layer->second;
}

// Quick way to sort barrel det layers by increasing R,
// do not abuse!
#include <TrackingTools/DetLayers/interface/BarrelDetLayer.h>
struct ExtractBarrelDetLayerR {
  typedef Surface::Scalar result_type;
  result_type operator()(const DetLayer* p) const {
    const BarrelDetLayer* bdl = dynamic_cast<const BarrelDetLayer*>(p);
    if (bdl)
      return bdl->specificSurface().radius();
    else
      return -1.;
  }
};

void MuonDetLayerGeometry::sortLayers() {
  // The following are filled inside-out, no need to re-sort
  // precomputed_value_sort(dtLayers.begin(), dtLayers.end(),ExtractR<DetLayer,float>());
  // precomputed_value_sort(cscLayers_fw.begin(), cscLayers_fw.end(),ExtractAbsZ<DetLayer,float>());
  // precomputed_value_sort(cscLayers_bk.begin(), cscLayers_bk.end(),ExtractAbsZ<DetLayer,float>());
  // precomputed_value_sort(rpcLayers_fw.begin(), rpcLayers_fw.end(),ExtractAbsZ<DetLayer,float>());
  // precomputed_value_sort(rpcLayers_bk.begin(), rpcLayers_bk.end(),ExtractAbsZ<DetLayer,float>());
  // precomputed_value_sort(rpcLayers_barrel.begin(), rpcLayers_barrel.end(), ExtractR<DetLayer,float>());

  // Sort these inside-out
  precomputed_value_sort(allBarrel.begin(), allBarrel.end(), ExtractBarrelDetLayerR());
  precomputed_value_sort(allBackward.begin(), allBackward.end(), ExtractAbsZ<DetLayer, float>());
  precomputed_value_sort(allForward.begin(), allForward.end(), ExtractAbsZ<DetLayer, float>());

  // Build more complicated vectors with correct sorting

  //cscLayers_all: from -Z to +Z
  cscLayers_all.reserve(cscLayers_bk.size() + cscLayers_fw.size());
  std::copy(cscLayers_bk.begin(), cscLayers_bk.end(), back_inserter(cscLayers_all));
  std::reverse(cscLayers_all.begin(), cscLayers_all.end());
  std::copy(cscLayers_fw.begin(), cscLayers_fw.end(), back_inserter(cscLayers_all));

  //gemLayers_all: from -Z to +Z
  gemLayers_all.reserve(gemLayers_bk.size() + gemLayers_fw.size());
  std::copy(gemLayers_bk.begin(), gemLayers_bk.end(), back_inserter(gemLayers_all));
  std::reverse(gemLayers_all.begin(), gemLayers_all.end());
  std::copy(gemLayers_fw.begin(), gemLayers_fw.end(), back_inserter(gemLayers_all));

  //me0Layers_all: from -Z to +Z
  me0Layers_all.reserve(me0Layers_bk.size() + me0Layers_fw.size());
  std::copy(me0Layers_bk.begin(), me0Layers_bk.end(), back_inserter(me0Layers_all));
  std::reverse(me0Layers_all.begin(), me0Layers_all.end());
  std::copy(me0Layers_fw.begin(), me0Layers_fw.end(), back_inserter(me0Layers_all));

  //rpcLayers_endcap: from -Z to +Z
  rpcLayers_endcap.reserve(rpcLayers_bk.size() + rpcLayers_fw.size());
  std::copy(rpcLayers_bk.begin(), rpcLayers_bk.end(), back_inserter(rpcLayers_endcap));
  std::reverse(rpcLayers_endcap.begin(), rpcLayers_endcap.end());
  std::copy(rpcLayers_fw.begin(), rpcLayers_fw.end(), back_inserter(rpcLayers_endcap));

  //rpcLayers_all: order is bw, barrel, fw
  rpcLayers_all.reserve(rpcLayers_bk.size() + rpcLayers_barrel.size() + rpcLayers_fw.size());
  std::copy(rpcLayers_bk.begin(), rpcLayers_bk.end(), back_inserter(rpcLayers_all));
  std::reverse(rpcLayers_all.begin(), rpcLayers_all.end());
  std::copy(rpcLayers_barrel.begin(), rpcLayers_barrel.end(), back_inserter(rpcLayers_all));
  std::copy(rpcLayers_fw.begin(), rpcLayers_fw.end(), back_inserter(rpcLayers_all));

  // allEndcap: order is  all bw, all fw
  allEndcap.reserve(allBackward.size() + allForward.size());
  std::copy(allBackward.begin(), allBackward.end(), back_inserter(allEndcap));
  std::reverse(allEndcap.begin(), allEndcap.end());
  std::copy(allForward.begin(), allForward.end(), back_inserter(allEndcap));

  // allEndcapCSCGEM: order is  all bw, all fw
  allEndcapCscGem.reserve(cscLayers_bk.size() + cscLayers_fw.size() + gemLayers_bk.size() + gemLayers_fw.size());
  std::copy(cscLayers_bk.begin(), cscLayers_bk.end(), back_inserter(allEndcapCscGem));
  std::copy(gemLayers_bk.begin(), gemLayers_bk.end(), back_inserter(allEndcapCscGem));
  std::reverse(allEndcapCscGem.begin(), allEndcapCscGem.end());
  std::copy(cscLayers_fw.begin(), cscLayers_fw.end(), back_inserter(allEndcapCscGem));
  std::copy(gemLayers_fw.begin(), gemLayers_fw.end(), back_inserter(allEndcapCscGem));

  // allCscGemForward
  allCscGemForward.reserve(cscLayers_fw.size() + gemLayers_fw.size());
  std::copy(cscLayers_fw.begin(), cscLayers_fw.end(), back_inserter(allCscGemForward));
  std::copy(gemLayers_fw.begin(), gemLayers_fw.end(), back_inserter(allCscGemForward));

  // allCscGemBackward
  allCscGemBackward.reserve(cscLayers_bk.size() + gemLayers_bk.size());
  std::copy(cscLayers_bk.begin(), cscLayers_bk.end(), back_inserter(allCscGemBackward));
  std::copy(gemLayers_bk.begin(), gemLayers_bk.end(), back_inserter(allCscGemBackward));

  // allCscME0Forward
  allCscME0Forward.reserve(cscLayers_fw.size() + me0Layers_fw.size());
  std::copy(cscLayers_fw.begin(), cscLayers_fw.end(), back_inserter(allCscME0Forward));
  std::copy(me0Layers_fw.begin(), me0Layers_fw.end(), back_inserter(allCscME0Forward));

  // allCscME0Backward
  allCscME0Backward.reserve(cscLayers_bk.size() + me0Layers_bk.size());
  std::copy(cscLayers_bk.begin(), cscLayers_bk.end(), back_inserter(allCscME0Backward));
  std::copy(me0Layers_bk.begin(), me0Layers_bk.end(), back_inserter(allCscME0Backward));

  // allEndcapCSCME0: order is  all bw, all fw
  allEndcapCscME0.reserve(cscLayers_bk.size() + cscLayers_fw.size() + me0Layers_bk.size() + me0Layers_fw.size());
  std::copy(cscLayers_bk.begin(), cscLayers_bk.end(), back_inserter(allEndcapCscME0));
  std::copy(me0Layers_bk.begin(), me0Layers_bk.end(), back_inserter(allEndcapCscME0));
  std::reverse(allEndcapCscME0.begin(), allEndcapCscME0.end());
  std::copy(cscLayers_fw.begin(), cscLayers_fw.end(), back_inserter(allEndcapCscME0));
  std::copy(me0Layers_fw.begin(), me0Layers_fw.end(), back_inserter(allEndcapCscME0));

  // allDetLayers: order is  all bw, all barrel, all fw
  allDetLayers.reserve(allBackward.size() + allBarrel.size() + allForward.size());
  std::copy(allBackward.begin(), allBackward.end(), back_inserter(allDetLayers));
  std::reverse(allDetLayers.begin(), allDetLayers.end());
  std::copy(allBarrel.begin(), allBarrel.end(), back_inserter(allDetLayers));
  std::copy(allForward.begin(), allForward.end(), back_inserter(allDetLayers));

  // number layers
  int sq = 0;
  for (auto l : allDetLayers)
    (*const_cast<DetLayer*>(l)).setSeqNum(sq++);
}