File indexing completed on 2023-03-17 11:19:23
0001
0002
0003
0004
0005
0006
0007
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
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
0046
0047 GEMCSCSegAlgoRR::~GEMCSCSegAlgoRR() {}
0048
0049
0050
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
0057
0058
0059
0060 theCSCLayers_ = csclayermap;
0061 theGEMEtaParts_ = gemrollmap;
0062
0063
0064 edm::LogVerbatim("GEMCSCSegAlgoRR") << "[GEMCSCSegAlgoRR::run] cached the csclayermap and the gemrollmap";
0065
0066
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
0078
0079
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
0092
0093 std::vector<GEMCSCSegment> segmentvectorfinal;
0094 std::vector<GEMCSCSegment> segmentvectorchamber;
0095 std::vector<GEMCSCSegment> segmentvectorfromfit;
0096 std::vector<const TrackingRecHit*> chainedRecHits;
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111 for (std::vector<const CSCSegment*>::const_iterator cscSegIt = cscsegments.begin(); cscSegIt != cscsegments.end();
0112 ++cscSegIt) {
0113
0114
0115
0116 chainedRecHits = this->chainHitsToSegm(*cscSegIt, gemrechits);
0117
0118
0119 if (chainedRecHits.size() > (*cscSegIt)->specificRecHits().size()) {
0120 segmentvectorfromfit = this->buildSegments(*cscSegIt, chainedRecHits);
0121 }
0122
0123
0124 else {
0125 std::vector<const GEMRecHit*> gemRecHits_noGEMrh;
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
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
0149
0150
0151
0152
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
0159
0160
0161
0162 auto segLP = cscsegment->localPosition();
0163 auto segLD = cscsegment->localDirection();
0164 auto cscrhs = cscsegment->specificRecHits();
0165
0166
0167
0168 for (auto crh = cscrhs.begin(); crh != cscrhs.end(); crh++) {
0169 chainedRecHits.push_back(crh->clone());
0170 }
0171
0172
0173 std::vector<const TrackingRecHit*>::const_iterator trhIt = chainedRecHits.begin();
0174
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
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
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
0202 for (grhIt = gemrechits.begin(); grhIt != gemrechits.end(); ++grhIt) {
0203
0204 auto rhLP = (*grhIt)->localPosition();
0205 const GEMEtaPartition* rhRef = theGEMEtaParts_.find((*grhIt)->gemId())->second;
0206 auto rhGP = rhRef->toGlobal(rhLP);
0207
0208
0209 auto rhLP_inSegmRef = cscChamber->toLocal(rhGP);
0210
0211
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
0219 LocalPoint extrPoint(xe, ye, ze);
0220
0221
0222
0223 auto rhGP_fromSegmRef = cscChamber->toGlobal(rhLP_inSegmRef);
0224 float phi_rh = rhGP_fromSegmRef.phi();
0225 float theta_rh = rhGP_fromSegmRef.theta();
0226
0227 auto extrPoinGP_fromSegmRef = cscChamber->toGlobal(extrPoint);
0228 float phi_ext = extrPoinGP_fromSegmRef.phi();
0229 float theta_ext = extrPoinGP_fromSegmRef.theta();
0230
0231
0232
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
0243
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 }
0255
0256
0257
0258
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 }
0270
0271 return chainedRecHits;
0272 }
0273
0274
0275
0276
0277
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
0285 if (rechits.size() < minHitsPerSegment) {
0286 return gemcscsegmentvector;
0287 }
0288
0289
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
0297 delete sfit_;
0298 sfit_ = new GEMCSCSegFit(theCSCLayers_, theGEMEtaParts_, rechits);
0299 sfit_->fit();
0300
0301
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
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 }