GEMCSCSegmentBuilder.cc

Go to the documentation of this file.

 
#include <RecoLocalMuon/GEMCSCSegment/plugins/GEMCSCSegmentBuilder.h>
#include <RecoLocalMuon/GEMCSCSegment/plugins/GEMCSCSegmentAlgorithm.h>
#include <RecoLocalMuon/GEMCSCSegment/plugins/GEMCSCSegmentBuilderPluginFactory.h>
 
#include "FWCore/Framework/interface/ESHandle.h"
#include <FWCore/Utilities/interface/Exception.h>
#include <FWCore/MessageLogger/interface/MessageLogger.h>
 
#include <Geometry/Records/interface/MuonGeometryRecord.h>
#include <Geometry/CSCGeometry/interface/CSCChamberSpecs.h>
#include <Geometry/CSCGeometry/interface/CSCLayer.h>
#include <Geometry/CSCGeometry/interface/CSCGeometry.h>
#include <Geometry/GEMGeometry/interface/GEMGeometry.h>
#include <Geometry/GEMGeometry/interface/GEMEtaPartition.h>
 
#include <DataFormats/MuonDetId/interface/CSCDetId.h>
#include <DataFormats/MuonDetId/interface/GEMDetId.h>
#include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
#include <DataFormats/CSCRecHit/interface/CSCRangeMapAccessor.h>
#include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
#include <DataFormats/CSCRecHit/interface/CSCSegmentCollection.h>
#include <DataFormats/CSCRecHit/interface/CSCSegment.h>
#include <DataFormats/GEMRecHit/interface/GEMRecHit.h>
#include <DataFormats/GEMRecHit/interface/GEMRecHitCollection.h>
 
GEMCSCSegmentBuilder::GEMCSCSegmentBuilder(const edm::ParameterSet& ps)
    :  // Ask factory to build this algorithm, giving it appropriate ParameterSet
      algo{GEMCSCSegmentBuilderPluginFactory::get()->create(ps.getParameter<std::string>("algo_name"),
                                                            ps.getParameter<edm::ParameterSet>("algo_psets"))},
      gemgeom_{nullptr},
      cscgeom_{nullptr} {
  edm::LogVerbatim("GEMCSCSegmentBuilder")
      << "[GEMCSCSegmentBuilder :: ctor] algorithm name: " << ps.getParameter<std::string>("algo_name");
}
 
GEMCSCSegmentBuilder::~GEMCSCSegmentBuilder() {
  edm::LogVerbatim("GEMCSCSegmentBuilder") << "[GEMCSCSegmentBuilder :: dstor] deleted the algorithm";
}
 
void GEMCSCSegmentBuilder::LinkGEMRollsToCSCChamberIndex(const GEMGeometry* gemGeo, const CSCGeometry* cscGeo) {
  for (TrackingGeometry::DetContainer::const_iterator it = gemGeo->dets().begin(); it < gemGeo->dets().end(); it++) {
    const GEMChamber* ch = dynamic_cast<const GEMChamber*>(*it);
    if (ch != nullptr) {
      std::vector<const GEMEtaPartition*> rolls = (ch->etaPartitions());
      for (std::vector<const GEMEtaPartition*>::const_iterator r = rolls.begin(); r != rolls.end(); ++r) {
        GEMDetId gemId = (*r)->id();
        int region = gemId.region();
        if (region != 0) {
          int station = gemId.station();
          int ring = gemId.ring();
          int gemchamber = gemId.chamber();
          int layer = gemId.layer();
          int cscring = ring;
          int cscstation = station;
          int cscchamber = gemchamber;
          int csclayer = layer;
          CSCStationIndex ind(region, cscstation, cscring, cscchamber, csclayer);
          std::set<GEMDetId> myrolls;
          if (rollstoreCSC.find(ind) != rollstoreCSC.end())
            myrolls = rollstoreCSC[ind];
          myrolls.insert(gemId);
          rollstoreCSC[ind] = myrolls;
        }
      }
    }
  }
 
  // edm::LogVerbatim to print details of std::map< CSCIndex, std::set<GEMRolls> >
  for (std::map<CSCStationIndex, std::set<GEMDetId> >::iterator mapit = rollstoreCSC.begin();
       mapit != rollstoreCSC.end();
       ++mapit) {
    CSCStationIndex map_first = mapit->first;
    std::set<GEMDetId> map_secnd = mapit->second;
    std::stringstream GEMRollsstream;
    for (std::set<GEMDetId>::iterator setit = map_secnd.begin(); setit != map_secnd.end(); ++setit) {
      GEMRollsstream << "[ GEM Id: " << setit->rawId() << " (" << *setit << ")"
                     << "]," << std::endl;
    }
    std::string GEMRollsstr = GEMRollsstream.str();
    edm::LogVerbatim("GEMCSCSegmentBuilder")
        << "[GEMCSCSegmentBuilder :: LinkGEMRollsToCSCChamberIndex] CSC Station Index :: [" << map_first.region() << ","
        << map_first.station() << "," << map_first.ring() << "," << map_first.chamber() << "," << map_first.layer()
        << "] has following GEM rolls: [" << GEMRollsstr << "]" << std::endl;
  }
}
 
void GEMCSCSegmentBuilder::build(const GEMRecHitCollection* recHits,
                                 const CSCSegmentCollection* cscsegments,
                                 GEMCSCSegmentCollection& oc) {
  edm::LogVerbatim("GEMCSCSegmentBuilder")
      << "[GEMCSCSegmentBuilder :: build] Total number of GEM rechits in this event: " << recHits->size()
      << " Total number of CSC segments in this event: " << cscsegments->size();
 
  // Structure of the build function:
  // --------------------------------
  // The build function is called once per Event
  // It loops over CSC Segment Collection and
  //   - identifies the CSC Chambers with segments
  //   - searches for rechits in GEM rolls associared
  //   - creates CSC segment collection with GEM RecHits associated (1)
  //   - creates CSC segment collection without GEM RecHits associated (2)
  // then
  //   - passes CSC segment collection (1) and GEM rechits to segment builder
  //   - makes a copy of CSC segment collection (2) and push it in the GEMCSC Segment collection
 
  // GEM Roll can contain more than one rechit     --> therefore create map <detId, vector<GEMRecHit> >
  // CSC Chamber can contain more than one segment --> therefore create map <detId, vector<CSCSegment> >
  std::map<uint32_t, std::vector<GEMRecHit*> > gemRecHitCollection;
  std::map<uint32_t, std::vector<const CSCSegment*> > cscSegColl_GEM11;
  std::map<uint32_t, std::vector<const CSCSegment*> > cscSegColl_noGEM;
 
  // === Loop over CSC Segment Collection ===
  // ========================================
  for (CSCSegmentCollection::const_iterator segmIt = cscsegments->begin(); segmIt != cscsegments->end(); ++segmIt) {
    CSCDetId CSCId = segmIt->cscDetId();
 
    // Search for Matches between GEM Roll and CSC Chamber
    //   - only matches between ME1/1(a&b) and GE1/1
    //   - notation: CSC ring 1 = ME1/1b; CSC ring 4 = ME1/1a
 
    // Case A :: ME1/1 Segments can have GEM Rechits associated to them
    // ================================================================
    if (CSCId.station() == 1 && (CSCId.ring() == 1 || CSCId.ring() == 4)) {
      edm::LogVerbatim("GEMCSCSegmentBuilder")
          << "[GEMCSCSegmentBuilder :: build] Found ME1/1 Segment in " << CSCId.rawId() << " = " << CSCId << std::endl;
 
      // 1) Save the CSC Segment in CSC segment collection
      // -------------------------------------------------
      // get CSC segment vector associated to this CSCDetId
      // if no vector is associated yet to this CSCDetId, create empty vector
      // then add the segment to the vector
      // and assign the vector again to the CSCDetId
      std::vector<const CSCSegment*> cscsegmentvector = cscSegColl_GEM11[CSCId.rawId()];
      cscsegmentvector.push_back(segmIt->clone());
      cscSegColl_GEM11[CSCId.rawId()] = cscsegmentvector;
 
      // 2) Make a vector of GEM DetIds associated to this CSC chamber of the segment
      // ----------------------------------------------------------------------------
      std::vector<GEMDetId> rollsForThisCSCvector;
 
      int cscRegion = CSCId.endcap();
      int cscStation = CSCId.station();
      int cscChamber = CSCId.chamber();
      int gemRegion = 1;
      if (cscRegion == 2)
        gemRegion = -1;  // CSC Endcaps are numbered 1,2; GEM Endcaps are numbered +1, -1
      int gemRing = 1;   // this we can hardcode, only GEMs in Ring 1
      int gemStation = cscStation;
 
      int gem1stChamber = cscChamber;
      // Just adding also neighbouring chambers here is not enough to get the GEM-CSC segment looking at overlapping chambers
      // Need to disentangle the use of the CSC chamber and GEM roll in the GEMCSCSegFit class
      // For now just ignore the neighbouring chambers and keep code commented out ... at later point we will include it again
      // int gem2ndChamber = gem1stChamber+1; if(gem2ndChamber>36) gem2ndChamber-=36; // neighbouring GEM chamber X+1
      // int gem3rdChamber = gem1stChamber-1; if(gem2ndChamber<1)  gem2ndChamber+=36; // neighbouring GEM chamber X-1
 
      std::vector<CSCStationIndex> indexvector;
      CSCStationIndex index11(gemRegion, gemStation, gemRing, gem1stChamber, 1);  // GEM Chamber Layer 1
      CSCStationIndex index12(gemRegion, gemStation, gemRing, gem1stChamber, 2);  // GEM Chamber Layer 2
      indexvector.push_back(index11);
      indexvector.push_back(index12);
      // for now not inserting neighbouring chambers and keep code commented out ... at later point we will include it again
      // CSCStationIndex index21(gemRegion,gemStation,gemRing,gem2ndChamber,1);         CSCStationIndex index22(gemRegion,gemStation,gemRing,gem2ndChamber,2);
      // CSCStationIndex index31(gemRegion,gemStation,gemRing,gem3rdChamber,1);         CSCStationIndex index32(gemRegion,gemStation,gemRing,gem3rdChamber,2);
      // indexvector.push_back(index21); indexvector.push_back(index22);
      // indexvector.push_back(index31); indexvector.push_back(index32);
 
      for (std::vector<CSCStationIndex>::iterator cscIndexIt = indexvector.begin(); cscIndexIt != indexvector.end();
           ++cscIndexIt) {
        std::set<GEMDetId> rollsForThisCSC = rollstoreCSC[*cscIndexIt];
        for (std::set<GEMDetId>::iterator gemRollIt = rollsForThisCSC.begin(); gemRollIt != rollsForThisCSC.end();
             ++gemRollIt) {
          rollsForThisCSCvector.push_back(*gemRollIt);
        }
      }
 
      // 3) Loop over GEM Rechits
      // ------------------------
      edm::LogVerbatim("GEMCSCSegmentBuilder")
          << "[GEMCSCSegmentBuilder :: build] Start Loop over GEM Rechits :: size = " << recHits->size()
          << " and number of Rolls for this CSC :: " << rollsForThisCSCvector.size() << std::endl;
      for (GEMRecHitCollection::const_iterator hitIt = recHits->begin(); hitIt != recHits->end(); ++hitIt) {
        GEMDetId gemIdfromHit = hitIt->gemId();
 
        // Loop over GEM rolls being pointed by a CSC segment and look for a match
        for (std::vector<GEMDetId>::iterator gemRollIt = rollsForThisCSCvector.begin();
             gemRollIt != rollsForThisCSCvector.end();
             ++gemRollIt) {
          GEMDetId gemIdfromCSC(*gemRollIt);
          if (gemIdfromHit == GEMDetId(gemIdfromCSC.rawId())) {
            // get GEM RecHit vector associated to this CSC DetId
            // if no vector is associated yet to this CSC DetId, create empty vector
            // then add the rechits to the vector
            // and assign the vector again to the CSC DetId
            std::vector<GEMRecHit*> gemrechitvector = gemRecHitCollection[CSCId.rawId()];
            // check whether this hit was already filled in the gemrechit vector
            bool hitfound = false;
            for (std::vector<GEMRecHit*>::const_iterator it = gemrechitvector.begin(); it != gemrechitvector.end();
                 ++it) {
              if (hitIt->gemId() == (*it)->gemId() && hitIt->localPosition() == (*it)->localPosition())
                hitfound = true;
            }
            if (!hitfound)
              gemrechitvector.push_back(hitIt->clone());
            gemRecHitCollection[CSCId.rawId()] = gemrechitvector;
            edm::LogVerbatim("GEMCSCSegmentBuilder")
                << "[GEMCSCSegmentBuilder :: build] GEM Rechit in " << hitIt->gemId() << "[" << hitIt->gemId().rawId()
                << "] added to CSC segment found in " << CSCId << " [" << CSCId.rawId() << "]" << std::endl;
          }
        }  // end Loop over GEM Rolls
      }    // end Loop over GEM RecHits
    }      // end requirement of CSC segment in ME1/1
 
    // Case B :: all other CSC Chambers have no GEM Chamber associated
    // ===============================================================
    else if (!(CSCId.station() == 1 && (CSCId.ring() == 1 || CSCId.ring() == 4))) {
      edm::LogVerbatim("GEMCSCSegmentBuilder") << "[GEMCSCSegmentBuilder :: build] Found non-ME1/1 Segment in "
                                               << CSCId.rawId() << " = " << CSCId << std::endl;
 
      // get CSC segment vector associated to this CSCDetId
      // if no vector is associated yet to this CSCDetId, create empty vector
      // then add the segment to the vector
      // and assign the vector again to the CSCDetId
      std::vector<const CSCSegment*> cscsegmentvector_noGEM = cscSegColl_noGEM[CSCId.rawId()];
      cscsegmentvector_noGEM.push_back(segmIt->clone());
      cscSegColl_noGEM[CSCId.rawId()] = cscsegmentvector_noGEM;
    }
 
    else {
    }  // no other option
 
  }  // end Loop over csc segments
 
  // === Now pass CSC Segments and GEM RecHits to the Segment Algorithm ===
  // ======================================================================
 
  // Here we do the first modification to the code of Raffaella
  //
  // Raffaella's code loops over the CSC Segments and assign GEM Rechits to them
  // Then it loops over the GemRechits and does the segmentbuilding using the
  // GEM Rechits and the assigned CSC segment. Now this relationship is not a one-to-one relationship (bijection)
  // In case there are more than 1 segment in the CSC, we have no control that the right supersegment is built
  //
  // Now we have a particular case where two muons passed through ME1/1,
  // resulting in 4 segments (2 real ones and 2 ghosts)
  // In this case the gemrechits are assigned to all 4 CSC segments
  // but the supersegment is only constructed once, instead of 4 trials,
  // resulting in 2 supersegments
  //
  // This is something I will change in the GEMCSCSegAlgoRR
  // In GEMCSCSegmentBuilder we just give all CSC segments and all GEM Rechits
  // belonging to one CSC chamber to the GEMCSC Segment Algo
 
  // case A :: Look for GEM Rechits associated to CSC Segments
  // =========================================================
  edm::LogVerbatim("GEMCSCSegmentBuilder")
      << "[GEMCSCSegmentBuilder :: build] Run over the gemRecHitCollection (size = " << gemRecHitCollection.size()
      << ") and build GEMCSC Segments";
  for (auto gemHitIt = gemRecHitCollection.begin(); gemHitIt != gemRecHitCollection.end(); ++gemHitIt) {
    CSCDetId cscId(gemHitIt->first);
 
    // hits
    std::vector<const CSCSegment*> cscSegments = cscSegColl_GEM11[cscId.rawId()];
    std::vector<const GEMRecHit*> gemRecHits;
    // dets
    std::map<uint32_t, const CSCLayer*> cscLayers;
    std::map<uint32_t, const GEMEtaPartition*> gemEtaPartitions;
 
    // fill dets for CSC
    std::vector<const CSCLayer*> cscLayerVector = cscgeom_->chamber(cscId.rawId())->layers();
    for (std::vector<const CSCLayer*>::const_iterator layIt = cscLayerVector.begin(); layIt != cscLayerVector.end();
         ++layIt) {
      cscLayers[(*layIt)->id()] = cscgeom_->layer((*layIt)->id());
    }
 
    // fill dets & hits for GEM
    for (auto rechit = gemHitIt->second.begin(); rechit != gemHitIt->second.end(); ++rechit) {
      GEMDetId gemid = (*rechit)->gemId();
      gemRecHits.push_back(*rechit);
      gemEtaPartitions[gemid.rawId()] = gemgeom_->etaPartition(gemid.rawId());
    }
 
    // LogDebug
    edm::LogVerbatim("GEMCSCSegmentBuilder")
        << "[GEMCSCSegmentBuilder :: build] ask SegmentAlgo to be run :: CSC DetId " << cscId.rawId() << " = " << cscId
        << " with " << cscSegments.size() << " CSC segments and " << gemRecHits.size() << " GEM rechits";
 
    // Ask the Segment Algorithm to build a GEMCSC Segment
    std::vector<GEMCSCSegment> segmentvector = algo->run(cscLayers, gemEtaPartitions, cscSegments, gemRecHits);
 
    // --- LogDebug --------------------------------------------------------
    edm::LogVerbatim("GEMCSCSegmentBuilder") << "[GEMCSCSegmentBuilder :: build] SegmentAlgo ran, now trying to add "
                                                "the segments from the returned GEMCSCSegment vector (with size = "
                                             << segmentvector.size() << ") to the master collection";
    std::stringstream segmentvectorss;
    segmentvectorss << "[GEMCSCSegmentBuilder :: build] :: GEMCSC segmentvector :: elements [" << std::endl;
    for (std::vector<GEMCSCSegment>::const_iterator segIt = segmentvector.begin(); segIt != segmentvector.end();
         ++segIt) {
      segmentvectorss << "[GEMCSC Segment details: \n" << *segIt << "]," << std::endl;
    }
    segmentvectorss << "]";
    std::string segmentvectorstr = segmentvectorss.str();
    edm::LogVerbatim("GEMCSCSegmentBuilder") << segmentvectorstr;
    // --- End LogDebug ----------------------------------------------------
 
    // Add the segments to master collection
    edm::LogVerbatim("GEMCSCSegmentBuilder")
        << "[GEMCSCSegmentBuilder :: build] |--> GEMCSC Segments created, now try to add to the collection :: CSC "
           "DetId: "
        << cscId.rawId() << " = " << cscId << " added " << segmentvector.size() << " GEMCSC Segments";
    oc.put(cscId, segmentvector.begin(), segmentvector.end());
    edm::LogVerbatim("GEMCSCSegmentBuilder")
        << "[GEMCSCSegmentBuilder :: build] |--> GEMCSC Segments added to the collection";
  }
 
  // Case B :: push CSC segments without GEM rechits
  // associated to them into the GEMCSC Segment Collection
  // =====================================================
  // Also for the case where no GEM chamber is associated to the CSC chambers
  // fill also GEMCSC segment, but this is basically a copy of the CSC segment
  // this will allow us to use only the GEMCSC segment collection in STA Mu Reco
  edm::LogVerbatim("GEMCSCSegmentBuilder")
      << "[GEMCSCSegmentBuilder :: build] Run over the CSC Segment collection without GEM (size = "
      << cscSegColl_noGEM.size() << " CSC chambers with segments) and wrap GEMCSC Segments";
  for (auto cscSegIt = cscSegColl_noGEM.begin(); cscSegIt != cscSegColl_noGEM.end(); ++cscSegIt) {
    CSCDetId cscId(cscSegIt->first);
 
    std::vector<const GEMRecHit*> gemRecHits_noGEM;  // make empty gem rechits vector
    std::vector<const CSCSegment*> cscSegments_noGEM = cscSegColl_noGEM[cscSegIt->first];
    std::vector<GEMCSCSegment> segmentvector;
 
    edm::LogVerbatim("GEMCSCSegmentBuilder")
        << "[GEMCSCSegmentBuilder :: build] |--> Run over the CSC Segment vector without GEM (size = "
        << cscSegments_noGEM.size() << " CSC segments) and wrap GEMCSC Segments";
    for (std::vector<const CSCSegment*>::iterator it = cscSegments_noGEM.begin(); it != cscSegments_noGEM.end(); ++it) {
      // wrap the CSC segment in the GEMCSC segment class
      GEMCSCSegment tmp(*it,
                        gemRecHits_noGEM,
                        (*it)->localPosition(),
                        (*it)->localDirection(),
                        (*it)->parametersError(),
                        (*it)->chi2());
      segmentvector.push_back(tmp);
    }
 
    // add the segments to the master collection
    edm::LogVerbatim("GEMCSCSegmentBuilder")
        << "[GEMCSCSegmentBuilder :: build] |--> CSC Segments wrapped, now try to add to the collection :: CSC DetId: "
        << cscId.rawId() << " = " << cscId << " added " << segmentvector.size() << " GEMCSC Segments";
    oc.put(cscId, segmentvector.begin(), segmentvector.end());
    edm::LogVerbatim("GEMCSCSegmentBuilder")
        << "[GEMCSCSegmentBuilder :: build] |--> CSC Segments added to the collection";
  }
 
  edm::LogVerbatim("GEMCSCSegmentBuilder") << "[GEMCSCSegmentBuilder :: build] job done !!!";
}
 
void GEMCSCSegmentBuilder::setGeometry(const GEMGeometry* gemgeom, const CSCGeometry* cscgeom) {
  gemgeom_ = gemgeom;
  cscgeom_ = cscgeom;
}