GEMCSCSegment/plugins/GEMCSCSegAlgoRR.cc

0001 /**
0002  * \file GEMCSCSegAlgoRR.cc
0003  * originally based on CSCSegAlgoDF.cc
0004  * modified by Piet Verwilligen
0005  * to use GEMCSCSegFit class
0006  *
0007  *  \authors: Raffaella Radogna
0008  */
0009
0010 #include "GEMCSCSegAlgoRR.h"
0011 #include "GEMCSCSegFit.h"
0012
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0017 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0018 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0019 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0020
0021 #include <algorithm>
0022 #include <cmath>
0023 #include <string>
0024 #include <iostream>
0025 #include <sstream>
0026
0027 /**
0028  *  Constructor
0029  */
0030 GEMCSCSegAlgoRR::GEMCSCSegAlgoRR(const edm::ParameterSet& ps)
0031     : GEMCSCSegmentAlgorithm(ps), myName("GEMCSCSegAlgoRR"), sfit_(nullptr) {
0032   debug = ps.getUntrackedParameter<bool>("GEMCSCDebug");
0033   minHitsPerSegment = ps.getParameter<unsigned int>("minHitsPerSegment");
0034   preClustering = ps.getParameter<bool>("preClustering");
0035   dXclusBoxMax = ps.getParameter<double>("dXclusBoxMax");
0036   dYclusBoxMax = ps.getParameter<double>("dYclusBoxMax");
0037   preClustering_useChaining = ps.getParameter<bool>("preClusteringUseChaining");
0038   dPhiChainBoxMax = ps.getParameter<double>("dPhiChainBoxMax");
0039   dThetaChainBoxMax = ps.getParameter<double>("dThetaChainBoxMax");
0040   dRChainBoxMax = ps.getParameter<double>("dRChainBoxMax");
0041   maxRecHitsInCluster = ps.getParameter<int>("maxRecHitsInCluster");
0042 }
0043
0044 /**
0045  * Destructor
0046  */
0047 GEMCSCSegAlgoRR::~GEMCSCSegAlgoRR() {}
0048
0049 /**
0050  * Run the algorithm
0051  */
0052 std::vector<GEMCSCSegment> GEMCSCSegAlgoRR::run(const std::map<uint32_t, const CSCLayer*>& csclayermap,
0053                                                 const std::map<uint32_t, const GEMEtaPartition*>& gemrollmap,
0054                                                 const std::vector<const CSCSegment*>& cscsegments,
0055                                                 const std::vector<const GEMRecHit*>& gemrechits) {
0056   // This function is called for each CSC Chamber (ME1/1) associated with a GEM chamber in which GEM Rechits were found
0057   // This means that the csc segment vector can contain more than one segment and that all combinations with the gemrechits have to be tried
0058
0059   // assign the maps < detid, geometry object > to the class variables
0060   theCSCLayers_ = csclayermap;
0061   theGEMEtaParts_ = gemrollmap;
0062
0063   // LogDebug info about reading of CSC Chamber map and GEM Eta Partition map
0064   edm::LogVerbatim("GEMCSCSegAlgoRR") << "[GEMCSCSegAlgoRR::run] cached the csclayermap and the gemrollmap";
0065
0066   // --- LogDebug for CSC Layer map -----------------------------------------
0067   std::stringstream csclayermapss;
0068   csclayermapss << "[GEMCSCSegAlgoRR::run] :: csclayermap :: elements [" << std::endl;
0069   for (std::map<uint32_t, const CSCLayer*>::const_iterator mapIt = theCSCLayers_.begin(); mapIt != theCSCLayers_.end();
0070        ++mapIt) {
0071     csclayermapss << "[CSC DetId " << mapIt->first << " =" << CSCDetId(mapIt->first) << ", CSC Layer " << mapIt->second
0072                   << " =" << (mapIt->second)->id() << "]," << std::endl;
0073   }
0074   csclayermapss << "]" << std::endl;
0075   std::string csclayermapstr = csclayermapss.str();
0076   edm::LogVerbatim("GEMCSCSegAlgoRR") << csclayermapstr;
0077   // --- End LogDebug -------------------------------------------------------
0078
0079   // --- LogDebug for GEM Eta Partition map ---------------------------------
0080   std::stringstream gemetapartmapss;
0081   gemetapartmapss << "[GEMCSCSegAlgoRR::run] :: gemetapartmap :: elements [" << std::endl;
0082   for (std::map<uint32_t, const GEMEtaPartition*>::const_iterator mapIt = theGEMEtaParts_.begin();
0083        mapIt != theGEMEtaParts_.end();
0084        ++mapIt) {
0085     gemetapartmapss << "[GEM DetId " << mapIt->first << " =" << GEMDetId(mapIt->first) << ", GEM EtaPart "
0086                     << mapIt->second << "]," << std::endl;
0087   }
0088   gemetapartmapss << "]" << std::endl;
0089   std::string gemetapartmapstr = gemetapartmapss.str();
0090   edm::LogVerbatim("GEMCSCSegAlgoRR") << gemetapartmapstr;
0091   // --- End LogDebug -------------------------------------------------------
0092
0093   std::vector<GEMCSCSegment> segmentvectorfinal;
0094   std::vector<GEMCSCSegment> segmentvectorchamber;
0095   std::vector<GEMCSCSegment> segmentvectorfromfit;
0096   std::vector<const TrackingRecHit*> chainedRecHits;
0097
0098   // From the GEMCSCSegmentBuilder we get
0099   // - a collection of CSC Segments belonging all to the same chamber
0100   // - a collection of GEM RecHits belonging to GEM Eta-Partitions close to this CSC
0101   //
0102   // Now we have to loop over all CSC Segments
0103   // and see to which CSC segments we can match the GEM rechits
0104   // if matching fails we will still keep those segments in the GEMCSC segment collection
0105   // but we will assign -1 to the matched parameter --> 2B implemented in the DataFormat
0106   // 1 for matched, 0 for no GEM chamber available and -1 for not matched
0107
0108   // empty the temporary gemcsc segments vector
0109   // segmentvectortmp.clear();
0110
0111   for (std::vector<const CSCSegment*>::const_iterator cscSegIt = cscsegments.begin(); cscSegIt != cscsegments.end();
0112        ++cscSegIt) {
0113     // chain hits :: make a vector of TrackingRecHits
0114     //               that contains the CSC and GEM rechits
0115     //               and that can be given to the fitter
0116     chainedRecHits = this->chainHitsToSegm(*cscSegIt, gemrechits);
0117
0118     // if gemrechits are associated, run the buildSegments step
0119     if (chainedRecHits.size() > (*cscSegIt)->specificRecHits().size()) {
0120       segmentvectorfromfit = this->buildSegments(*cscSegIt, chainedRecHits);
0121     }
0122
0123     // if no gemrechits are associated, wrap the existing CSCSegment in a GEMCSCSegment class
0124     else {
0125       std::vector<const GEMRecHit*> gemRecHits_noGEMrh;  // empty GEMRecHit vector
0126       GEMCSCSegment tmp(*cscSegIt,
0127                         gemRecHits_noGEMrh,
0128                         (*cscSegIt)->localPosition(),
0129                         (*cscSegIt)->localDirection(),
0130                         (*cscSegIt)->parametersError(),
0131                         (*cscSegIt)->chi2());
0132       segmentvectorfromfit.push_back(tmp);
0133     }
0134
0135     segmentvectorchamber.insert(segmentvectorchamber.end(), segmentvectorfromfit.begin(), segmentvectorfromfit.end());
0136     segmentvectorfromfit.clear();
0137   }
0138
0139   // add the found gemcsc segments to the collection of all gemcsc segments and return
0140   segmentvectorfinal.insert(segmentvectorfinal.end(), segmentvectorchamber.begin(), segmentvectorchamber.end());
0141   segmentvectorchamber.clear();
0142   edm::LogVerbatim("GEMCSCSegAlgoRR") << "[GEMCSCSegAlgoRR::run] GEMCSC Segments fitted or wrapped, size of vector = "
0143                                       << segmentvectorfinal.size();
0144   return segmentvectorfinal;
0145 }
0146
0147 /**
0148  * Chain hits :: make a TrackingRecHit vector containing both CSC and GEM rechits
0149  *               take the hits from the CSCSegment and add a clone of them
0150  *               take the hits from the GEM RecHit vector and
0151  *               check whether they are compatible with the CSC segment extrapolation
0152  *               to the GEM plane. If they are compatible, add these GEM rechits
0153  */
0154 std::vector<const TrackingRecHit*> GEMCSCSegAlgoRR::chainHitsToSegm(const CSCSegment* cscsegment,
0155                                                                     const std::vector<const GEMRecHit*>& gemrechits) {
0156   std::vector<const TrackingRecHit*> chainedRecHits;
0157
0158   // working with layers makes it here a bit more difficult:
0159   // the CSC segment points to the chamber, which we cannot ask from the map
0160   // the CSC rechits point to the layers, from which we can ask the chamber
0161
0162   auto segLP = cscsegment->localPosition();
0163   auto segLD = cscsegment->localDirection();
0164   auto cscrhs = cscsegment->specificRecHits();
0165
0166   // Loop over the CSC Segment rechits
0167   // and save a copy in the chainedRecHits vector
0168   for (auto crh = cscrhs.begin(); crh != cscrhs.end(); crh++) {
0169     chainedRecHits.push_back(crh->clone());
0170   }
0171
0172   // now ask the layer id of the first CSC rechit
0173   std::vector<const TrackingRecHit*>::const_iterator trhIt = chainedRecHits.begin();
0174   // make sure pointer is valid
0175   if (trhIt == chainedRecHits.end()) {
0176     edm::LogVerbatim("GEMCSCSegFit")
0177         << "[GEMCSCSegFit::chainHitsToSegm] CSC segment has zero rechits ... end function here";
0178     return chainedRecHits;
0179   }
0180   const CSCLayer* cscLayer = theCSCLayers_.find((*trhIt)->rawId())->second;
0181   // now ask the chamber id of the first CSC rechit
0182   if (!cscLayer) {
0183     edm::LogVerbatim("GEMCSCSegFit")
0184         << "[GEMCSCSegFit::chainHitsToSegm] CSC rechit ID was not found back in the CSCLayerMap ... end function here";
0185     return chainedRecHits;
0186   }
0187   const CSCChamber* cscChamber = cscLayer->chamber();
0188
0189   // For non-empty GEM rechit vector
0190   if (!gemrechits.empty()) {
0191     float Dphi_min_l1 = 999;
0192     float Dphi_min_l2 = 999;
0193     float Dtheta_min_l1 = 999;
0194     float Dtheta_min_l2 = 999;
0195
0196     std::vector<const GEMRecHit*>::const_iterator grhIt = gemrechits.begin();
0197
0198     const GEMRecHit* gemrh_min_l1 = *grhIt;
0199     const GEMRecHit* gemrh_min_l2 = *grhIt;
0200
0201     // Loop over GEM rechits from the EnsembleGEMHitContainer
0202     for (grhIt = gemrechits.begin(); grhIt != gemrechits.end(); ++grhIt) {
0203       // get GEM Rechit Local & Global Position
0204       auto rhLP = (*grhIt)->localPosition();
0205       const GEMEtaPartition* rhRef = theGEMEtaParts_.find((*grhIt)->gemId())->second;
0206       auto rhGP = rhRef->toGlobal(rhLP);
0207       // get GEM Rechit Local position w.r.t. the CSC Chamber
0208       // --> we are interessed in the z-coordinate
0209       auto rhLP_inSegmRef = cscChamber->toLocal(rhGP);
0210       // calculate the extrapolation of the CSC segment to the GEM plane (z-coord)
0211       // to get x- and y- coordinate
0212       float xe = 0.0, ye = 0.0, ze = 0.0;
0213       if (segLD.z() != 0) {
0214         xe = segLP.x() + segLD.x() * rhLP_inSegmRef.z() / segLD.z();
0215         ye = segLP.y() + segLD.y() * rhLP_inSegmRef.z() / segLD.z();
0216         ze = rhLP_inSegmRef.z();
0217       }
0218       // 3D extrapolated point in the GEM plane
0219       LocalPoint extrPoint(xe, ye, ze);
0220
0221       // get GEM Rechit Global Position, but obtained from CSC Chamber
0222       // --> this should be the same as rhGP = rhRef->toGlobal(rhLP);
0223       auto rhGP_fromSegmRef = cscChamber->toGlobal(rhLP_inSegmRef);
0224       float phi_rh = rhGP_fromSegmRef.phi();
0225       float theta_rh = rhGP_fromSegmRef.theta();
0226       // get Extrapolat Global Position, also obtained from CSC Chamber
0227       auto extrPoinGP_fromSegmRef = cscChamber->toGlobal(extrPoint);
0228       float phi_ext = extrPoinGP_fromSegmRef.phi();
0229       float theta_ext = extrPoinGP_fromSegmRef.theta();
0230
0231       // GEM 1st Layer :: Search for GEM Rechit with smallest Delta Eta and Delta Phi
0232       //                  and save it inside gemrh_min_l1
0233       if ((*grhIt)->gemId().layer() == 1) {
0234         float Dphi_l1 = fabs(phi_ext - phi_rh);
0235         float Dtheta_l1 = fabs(theta_ext - theta_rh);
0236         if (Dphi_l1 <= Dphi_min_l1) {
0237           Dphi_min_l1 = Dphi_l1;
0238           Dtheta_min_l1 = Dtheta_l1;
0239           gemrh_min_l1 = *grhIt;
0240         }
0241       }
0242       // GEM 2nd Layer :: Search for GEM Rechit with smallest Delta Eta and Delta Phi
0243       //                  and save it inside gemrh_min_l2
0244       if ((*grhIt)->gemId().layer() == 2) {
0245         float Dphi_l2 = fabs(phi_ext - phi_rh);
0246         float Dtheta_l2 = fabs(theta_ext - theta_rh);
0247         if (Dphi_l2 <= Dphi_min_l2) {
0248           Dphi_min_l2 = Dphi_l2;
0249           Dtheta_min_l2 = Dtheta_l2;
0250           gemrh_min_l2 = *grhIt;
0251         }
0252       }
0253
0254     }  // end loop over GEM Rechits
0255
0256     // Check whether GEM rechit with smallest delta eta and delta phi
0257     // w.r.t. the extrapolation of the CSC segment is within the
0258     // maxima given by the configuration of the algorithm
0259     bool phiRequirementOK_l1 = Dphi_min_l1 < dPhiChainBoxMax;
0260     bool thetaRequirementOK_l1 = Dtheta_min_l1 < dThetaChainBoxMax;
0261     if (phiRequirementOK_l1 && thetaRequirementOK_l1 && gemrh_min_l1 != nullptr) {
0262       chainedRecHits.push_back(gemrh_min_l1->clone());
0263     }
0264     bool phiRequirementOK_l2 = Dphi_min_l2 < dPhiChainBoxMax;
0265     bool thetaRequirementOK_l2 = Dtheta_min_l2 < dThetaChainBoxMax;
0266     if (phiRequirementOK_l2 && thetaRequirementOK_l2 && gemrh_min_l2 != nullptr) {
0267       chainedRecHits.push_back(gemrh_min_l2->clone());
0268     }
0269   }  // End check > 0 GEM rechits
0270
0271   return chainedRecHits;
0272 }
0273
0274 /**
0275  * This algorithm uses a Minimum Spanning Tree (ST) approach to build
0276  * endcap muon track segments from the rechits in the 6 layers of a CSC
0277  * and the 2 layers inside a GE1/1 GEM. Fit is implemented in GEMCSCSegFit.cc
0278  */
0279 std::vector<GEMCSCSegment> GEMCSCSegAlgoRR::buildSegments(const CSCSegment* cscsegment,
0280                                                           const std::vector<const TrackingRecHit*>& rechits) {
0281   std::vector<GEMCSCSegment> gemcscsegmentvector;
0282   std::vector<const GEMRecHit*> gemrechits;
0283
0284   // Leave the function if the amount of rechits in the EnsembleHitContainer < min required
0285   if (rechits.size() < minHitsPerSegment) {
0286     return gemcscsegmentvector;
0287   }
0288
0289   // Extract the GEMRecHits from the TrackingRecHit vector
0290   for (std::vector<const TrackingRecHit*>::const_iterator trhIt = rechits.begin(); trhIt != rechits.end(); ++trhIt) {
0291     if (DetId((*trhIt)->rawId()).subdetId() == MuonSubdetId::GEM) {
0292       gemrechits.push_back(((const GEMRecHit*)*trhIt));
0293     }
0294   }
0295
0296   // The actual fit on all hits of the vector of the selected Tracking RecHits:
0297   delete sfit_;
0298   sfit_ = new GEMCSCSegFit(theCSCLayers_, theGEMEtaParts_, rechits);
0299   sfit_->fit();
0300
0301   // obtain all information necessary to make the segment:
0302   LocalPoint protoIntercept = sfit_->intercept();
0303   LocalVector protoDirection = sfit_->localdir();
0304   AlgebraicSymMatrix protoErrors = sfit_->covarianceMatrix();
0305   double protoChi2 = sfit_->chi2();
0306
0307   edm::LogVerbatim("GEMCSCSegAlgoRR")
0308       << "[GEMCSCSegAlgoRR::buildSegments] fit done, will now try to make GEMCSC Segment from CSC Segment in "
0309       << cscsegment->cscDetId() << " made of " << cscsegment->specificRecHits().size()
0310       << " rechits and with chi2 = " << cscsegment->chi2() << " and with " << gemrechits.size() << " GEM RecHits";
0311
0312   // save all information inside GEMCSCSegment
0313   GEMCSCSegment tmp(cscsegment, gemrechits, protoIntercept, protoDirection, protoErrors, protoChi2);
0314
0315   edm::LogVerbatim("GEMCSCSegAlgoRR") << "[GEMCSCSegAlgoRR::buildSegments] GEMCSC Segment made";
0316   gemcscsegmentvector.push_back(tmp);
0317   return gemcscsegmentvector;
0318 }