GEMCSCSegAlgoRR.cc

Go to the documentation of this file.

#include "GEMCSCSegAlgoRR.h"
#include "GEMCSCSegFit.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Geometry/CSCGeometry/interface/CSCLayer.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"

#include <algorithm>
#include <cmath>
#include <string>
#include <iostream>
#include <sstream>

GEMCSCSegAlgoRR::GEMCSCSegAlgoRR(const edm::ParameterSet& ps)
    : GEMCSCSegmentAlgorithm(ps), myName("GEMCSCSegAlgoRR"), sfit_(nullptr) {
  debug = ps.getUntrackedParameter<bool>("GEMCSCDebug");
  minHitsPerSegment = ps.getParameter<unsigned int>("minHitsPerSegment");
  preClustering = ps.getParameter<bool>("preClustering");
  dXclusBoxMax = ps.getParameter<double>("dXclusBoxMax");
  dYclusBoxMax = ps.getParameter<double>("dYclusBoxMax");
  preClustering_useChaining = ps.getParameter<bool>("preClusteringUseChaining");
  dPhiChainBoxMax = ps.getParameter<double>("dPhiChainBoxMax");
  dThetaChainBoxMax = ps.getParameter<double>("dThetaChainBoxMax");
  dRChainBoxMax = ps.getParameter<double>("dRChainBoxMax");
  maxRecHitsInCluster = ps.getParameter<int>("maxRecHitsInCluster");
}

GEMCSCSegAlgoRR::~GEMCSCSegAlgoRR() {}

std::vector<GEMCSCSegment> GEMCSCSegAlgoRR::run(const std::map<uint32_t, const CSCLayer*>& csclayermap,
                                                const std::map<uint32_t, const GEMEtaPartition*>& gemrollmap,
                                                const std::vector<const CSCSegment*>& cscsegments,
                                                const std::vector<const GEMRecHit*>& gemrechits) {
  // This function is called for each CSC Chamber (ME1/1) associated with a GEM chamber in which GEM Rechits were found
  // This means that the csc segment vector can contain more than one segment and that all combinations with the gemrechits have to be tried

  // assign the maps < detid, geometry object > to the class variables
  theCSCLayers_ = csclayermap;
  theGEMEtaParts_ = gemrollmap;

  // LogDebug info about reading of CSC Chamber map and GEM Eta Partition map
  edm::LogVerbatim("GEMCSCSegAlgoRR") << "[GEMCSCSegAlgoRR::run] cached the csclayermap and the gemrollmap";

  // --- LogDebug for CSC Layer map -----------------------------------------
  std::stringstream csclayermapss;
  csclayermapss << "[GEMCSCSegAlgoRR::run] :: csclayermap :: elements [" << std::endl;
  for (std::map<uint32_t, const CSCLayer*>::const_iterator mapIt = theCSCLayers_.begin(); mapIt != theCSCLayers_.end();
       ++mapIt) {
    csclayermapss << "[CSC DetId " << mapIt->first << " =" << CSCDetId(mapIt->first) << ", CSC Layer " << mapIt->second
                  << " =" << (mapIt->second)->id() << "]," << std::endl;
  }
  csclayermapss << "]" << std::endl;
  std::string csclayermapstr = csclayermapss.str();
  edm::LogVerbatim("GEMCSCSegAlgoRR") << csclayermapstr;
  // --- End LogDebug -------------------------------------------------------

  // --- LogDebug for GEM Eta Partition map ---------------------------------
  std::stringstream gemetapartmapss;
  gemetapartmapss << "[GEMCSCSegAlgoRR::run] :: gemetapartmap :: elements [" << std::endl;
  for (std::map<uint32_t, const GEMEtaPartition*>::const_iterator mapIt = theGEMEtaParts_.begin();
       mapIt != theGEMEtaParts_.end();
       ++mapIt) {
    gemetapartmapss << "[GEM DetId " << mapIt->first << " =" << GEMDetId(mapIt->first) << ", GEM EtaPart "
                    << mapIt->second << "]," << std::endl;
  }
  gemetapartmapss << "]" << std::endl;
  std::string gemetapartmapstr = gemetapartmapss.str();
  edm::LogVerbatim("GEMCSCSegAlgoRR") << gemetapartmapstr;
  // --- End LogDebug -------------------------------------------------------

  std::vector<GEMCSCSegment> segmentvectorfinal;
  std::vector<GEMCSCSegment> segmentvectorchamber;
  std::vector<GEMCSCSegment> segmentvectorfromfit;
  std::vector<const TrackingRecHit*> chainedRecHits;

  // From the GEMCSCSegmentBuilder we get
  // - a collection of CSC Segments belonging all to the same chamber
  // - a collection of GEM RecHits belonging to GEM Eta-Partitions close to this CSC
  //
  // Now we have to loop over all CSC Segments
  // and see to which CSC segments we can match the GEM rechits
  // if matching fails we will still keep those segments in the GEMCSC segment collection
  // but we will assign -1 to the matched parameter --> 2B implemented in the DataFormat
  // 1 for matched, 0 for no GEM chamber available and -1 for not matched

  // empty the temporary gemcsc segments vector
  // segmentvectortmp.clear();

  for (std::vector<const CSCSegment*>::const_iterator cscSegIt = cscsegments.begin(); cscSegIt != cscsegments.end();
       ++cscSegIt) {
    // chain hits :: make a vector of TrackingRecHits
    //               that contains the CSC and GEM rechits
    //               and that can be given to the fitter
    chainedRecHits = this->chainHitsToSegm(*cscSegIt, gemrechits);

    // if gemrechits are associated, run the buildSegments step
    if (chainedRecHits.size() > (*cscSegIt)->specificRecHits().size()) {
      segmentvectorfromfit = this->buildSegments(*cscSegIt, chainedRecHits);
    }

    // if no gemrechits are associated, wrap the existing CSCSegment in a GEMCSCSegment class
    else {
      std::vector<const GEMRecHit*> gemRecHits_noGEMrh;  // empty GEMRecHit vector
      GEMCSCSegment tmp(*cscSegIt,
                        gemRecHits_noGEMrh,
                        (*cscSegIt)->localPosition(),
                        (*cscSegIt)->localDirection(),
                        (*cscSegIt)->parametersError(),
                        (*cscSegIt)->chi2());
      segmentvectorfromfit.push_back(tmp);
    }

    segmentvectorchamber.insert(segmentvectorchamber.end(), segmentvectorfromfit.begin(), segmentvectorfromfit.end());
    segmentvectorfromfit.clear();
  }

  // add the found gemcsc segments to the collection of all gemcsc segments and return
  segmentvectorfinal.insert(segmentvectorfinal.end(), segmentvectorchamber.begin(), segmentvectorchamber.end());
  segmentvectorchamber.clear();
  edm::LogVerbatim("GEMCSCSegAlgoRR") << "[GEMCSCSegAlgoRR::run] GEMCSC Segments fitted or wrapped, size of vector = "
                                      << segmentvectorfinal.size();
  return segmentvectorfinal;
}

std::vector<const TrackingRecHit*> GEMCSCSegAlgoRR::chainHitsToSegm(const CSCSegment* cscsegment,
                                                                    const std::vector<const GEMRecHit*>& gemrechits) {
  std::vector<const TrackingRecHit*> chainedRecHits;

  // working with layers makes it here a bit more difficult:
  // the CSC segment points to the chamber, which we cannot ask from the map
  // the CSC rechits point to the layers, from which we can ask the chamber

  auto segLP = cscsegment->localPosition();
  auto segLD = cscsegment->localDirection();
  auto cscrhs = cscsegment->specificRecHits();

  // Loop over the CSC Segment rechits
  // and save a copy in the chainedRecHits vector
  for (auto crh = cscrhs.begin(); crh != cscrhs.end(); crh++) {
    chainedRecHits.push_back(crh->clone());
  }

  // now ask the layer id of the first CSC rechit
  std::vector<const TrackingRecHit*>::const_iterator trhIt = chainedRecHits.begin();
  // make sure pointer is valid
  if (trhIt == chainedRecHits.end()) {
    edm::LogVerbatim("GEMCSCSegFit")
        << "[GEMCSCSegFit::chainHitsToSegm] CSC segment has zero rechits ... end function here";
    return chainedRecHits;
  }
  const CSCLayer* cscLayer = theCSCLayers_.find((*trhIt)->rawId())->second;
  // now ask the chamber id of the first CSC rechit
  if (!cscLayer) {
    edm::LogVerbatim("GEMCSCSegFit")
        << "[GEMCSCSegFit::chainHitsToSegm] CSC rechit ID was not found back in the CSCLayerMap ... end function here";
    return chainedRecHits;
  }
  const CSCChamber* cscChamber = cscLayer->chamber();

  // For non-empty GEM rechit vector
  if (!gemrechits.empty()) {
    float Dphi_min_l1 = 999;
    float Dphi_min_l2 = 999;
    float Dtheta_min_l1 = 999;
    float Dtheta_min_l2 = 999;

    std::vector<const GEMRecHit*>::const_iterator grhIt = gemrechits.begin();

    const GEMRecHit* gemrh_min_l1 = *grhIt;
    const GEMRecHit* gemrh_min_l2 = *grhIt;

    // Loop over GEM rechits from the EnsembleGEMHitContainer
    for (grhIt = gemrechits.begin(); grhIt != gemrechits.end(); ++grhIt) {
      // get GEM Rechit Local & Global Position
      auto rhLP = (*grhIt)->localPosition();
      const GEMEtaPartition* rhRef = theGEMEtaParts_.find((*grhIt)->gemId())->second;
      auto rhGP = rhRef->toGlobal(rhLP);
      // get GEM Rechit Local position w.r.t. the CSC Chamber
      // --> we are interessed in the z-coordinate
      auto rhLP_inSegmRef = cscChamber->toLocal(rhGP);
      // calculate the extrapolation of the CSC segment to the GEM plane (z-coord)
      // to get x- and y- coordinate
      float xe = 0.0, ye = 0.0, ze = 0.0;
      if (segLD.z() != 0) {
        xe = segLP.x() + segLD.x() * rhLP_inSegmRef.z() / segLD.z();
        ye = segLP.y() + segLD.y() * rhLP_inSegmRef.z() / segLD.z();
        ze = rhLP_inSegmRef.z();
      }
      // 3D extrapolated point in the GEM plane
      LocalPoint extrPoint(xe, ye, ze);

      // get GEM Rechit Global Position, but obtained from CSC Chamber
      // --> this should be the same as rhGP = rhRef->toGlobal(rhLP);
      auto rhGP_fromSegmRef = cscChamber->toGlobal(rhLP_inSegmRef);
      float phi_rh = rhGP_fromSegmRef.phi();
      float theta_rh = rhGP_fromSegmRef.theta();
      // get Extrapolat Global Position, also obtained from CSC Chamber
      auto extrPoinGP_fromSegmRef = cscChamber->toGlobal(extrPoint);
      float phi_ext = extrPoinGP_fromSegmRef.phi();
      float theta_ext = extrPoinGP_fromSegmRef.theta();

      // GEM 1st Layer :: Search for GEM Rechit with smallest Delta Eta and Delta Phi
      //                  and save it inside gemrh_min_l1
      if ((*grhIt)->gemId().layer() == 1) {
        float Dphi_l1 = fabs(phi_ext - phi_rh);
        float Dtheta_l1 = fabs(theta_ext - theta_rh);
        if (Dphi_l1 <= Dphi_min_l1) {
          Dphi_min_l1 = Dphi_l1;
          Dtheta_min_l1 = Dtheta_l1;
          gemrh_min_l1 = *grhIt;
        }
      }
      // GEM 2nd Layer :: Search for GEM Rechit with smallest Delta Eta and Delta Phi
      //                  and save it inside gemrh_min_l2
      if ((*grhIt)->gemId().layer() == 2) {
        float Dphi_l2 = fabs(phi_ext - phi_rh);
        float Dtheta_l2 = fabs(theta_ext - theta_rh);
        if (Dphi_l2 <= Dphi_min_l2) {
          Dphi_min_l2 = Dphi_l2;
          Dtheta_min_l2 = Dtheta_l2;
          gemrh_min_l2 = *grhIt;
        }
      }

    }  // end loop over GEM Rechits

    // Check whether GEM rechit with smallest delta eta and delta phi
    // w.r.t. the extrapolation of the CSC segment is within the
    // maxima given by the configuration of the algorithm
    bool phiRequirementOK_l1 = Dphi_min_l1 < dPhiChainBoxMax;
    bool thetaRequirementOK_l1 = Dtheta_min_l1 < dThetaChainBoxMax;
    if (phiRequirementOK_l1 && thetaRequirementOK_l1 && gemrh_min_l1 != nullptr) {
      chainedRecHits.push_back(gemrh_min_l1->clone());
    }
    bool phiRequirementOK_l2 = Dphi_min_l2 < dPhiChainBoxMax;
    bool thetaRequirementOK_l2 = Dtheta_min_l2 < dThetaChainBoxMax;
    if (phiRequirementOK_l2 && thetaRequirementOK_l2 && gemrh_min_l2 != nullptr) {
      chainedRecHits.push_back(gemrh_min_l2->clone());
    }
  }  // End check > 0 GEM rechits

  return chainedRecHits;
}

std::vector<GEMCSCSegment> GEMCSCSegAlgoRR::buildSegments(const CSCSegment* cscsegment,
                                                          const std::vector<const TrackingRecHit*>& rechits) {
  std::vector<GEMCSCSegment> gemcscsegmentvector;
  std::vector<const GEMRecHit*> gemrechits;

  // Leave the function if the amount of rechits in the EnsembleHitContainer < min required
  if (rechits.size() < minHitsPerSegment) {
    return gemcscsegmentvector;
  }

  // Extract the GEMRecHits from the TrackingRecHit vector
  for (std::vector<const TrackingRecHit*>::const_iterator trhIt = rechits.begin(); trhIt != rechits.end(); ++trhIt) {
    if (DetId((*trhIt)->rawId()).subdetId() == MuonSubdetId::GEM) {
      gemrechits.push_back(((const GEMRecHit*)*trhIt));
    }
  }

  // The actual fit on all hits of the vector of the selected Tracking RecHits:
  delete sfit_;
  sfit_ = new GEMCSCSegFit(theCSCLayers_, theGEMEtaParts_, rechits);
  sfit_->fit();

  // obtain all information necessary to make the segment:
  LocalPoint protoIntercept = sfit_->intercept();
  LocalVector protoDirection = sfit_->localdir();
  AlgebraicSymMatrix protoErrors = sfit_->covarianceMatrix();
  double protoChi2 = sfit_->chi2();

  edm::LogVerbatim("GEMCSCSegAlgoRR")
      << "[GEMCSCSegAlgoRR::buildSegments] fit done, will now try to make GEMCSC Segment from CSC Segment in "
      << cscsegment->cscDetId() << " made of " << cscsegment->specificRecHits().size()
      << " rechits and with chi2 = " << cscsegment->chi2() << " and with " << gemrechits.size() << " GEM RecHits";

  // save all information inside GEMCSCSegment
  GEMCSCSegment tmp(cscsegment, gemrechits, protoIntercept, protoDirection, protoErrors, protoChi2);

  edm::LogVerbatim("GEMCSCSegAlgoRR") << "[GEMCSCSegAlgoRR::buildSegments] GEMCSC Segment made";
  gemcscsegmentvector.push_back(tmp);
  return gemcscsegmentvector;
}