File indexing completed on 2023-03-17 11:19:19
0001
0002
0003
0004
0005
0006 #include "RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010
0011 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
0012
0013 #include "RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.h"
0014 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
0015
0016 #include "DataFormats/Common/interface/OwnVector.h"
0017 #include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
0018 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0019 #include "DataFormats/DTRecHit/interface/DTRecHit1DPair.h"
0020 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0021 #include "DataFormats/DTRecHit/interface/DTChamberRecSegment2D.h"
0022
0023 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0024
0025 using namespace std;
0026 using namespace edm;
0027
0028
0029
0030
0031
0032 DTCombinatorialPatternReco4D::DTCombinatorialPatternReco4D(const ParameterSet& pset, ConsumesCollector cc)
0033 : DTRecSegment4DBaseAlgo(pset), theAlgoName("DTCombinatorialPatternReco4D"), theDTGeometryToken(cc.esConsumes()) {
0034
0035 debug = pset.getUntrackedParameter<bool>("debug");
0036
0037
0038 applyT0corr = pset.getParameter<bool>("performT0SegCorrection");
0039
0040 computeT0corr = pset.existsAs<bool>("computeT0Seg") ? pset.getParameter<bool>("computeT0Seg") : true;
0041
0042
0043 theUpdator = new DTSegmentUpdator(pset, cc);
0044
0045
0046
0047
0048
0049
0050 allDTRecHits = pset.getParameter<bool>("AllDTRecHits");
0051
0052
0053
0054 the2DAlgo = new DTCombinatorialPatternReco(pset.getParameter<ParameterSet>("Reco2DAlgoConfig"), cc);
0055 }
0056
0057 DTCombinatorialPatternReco4D::~DTCombinatorialPatternReco4D() {
0058 delete the2DAlgo;
0059 delete theUpdator;
0060 }
0061
0062 void DTCombinatorialPatternReco4D::setES(const EventSetup& setup) {
0063 theDTGeometry = setup.getHandle(theDTGeometryToken);
0064 the2DAlgo->setES(setup);
0065 theUpdator->setES(setup);
0066 }
0067
0068 void DTCombinatorialPatternReco4D::setChamber(const DTChamberId& chId) {
0069
0070 theChamber = theDTGeometry->chamber(chId);
0071 }
0072
0073 void DTCombinatorialPatternReco4D::setDTRecHit1DContainer(Handle<DTRecHitCollection> all1DHits) {
0074 theHitsFromPhi1.clear();
0075 theHitsFromPhi2.clear();
0076 theHitsFromTheta.clear();
0077
0078 DTRecHitCollection::range rangeHitsFromPhi1 =
0079 all1DHits->get(DTRangeMapAccessor::layersBySuperLayer(DTSuperLayerId(theChamber->id(), 1)));
0080 DTRecHitCollection::range rangeHitsFromPhi2 =
0081 all1DHits->get(DTRangeMapAccessor::layersBySuperLayer(DTSuperLayerId(theChamber->id(), 3)));
0082
0083 vector<DTRecHit1DPair> hitsFromPhi1(rangeHitsFromPhi1.first, rangeHitsFromPhi1.second);
0084 vector<DTRecHit1DPair> hitsFromPhi2(rangeHitsFromPhi2.first, rangeHitsFromPhi2.second);
0085 if (debug)
0086 cout << "Number of DTRecHit1DPair in the SL 1 (Phi 1): " << hitsFromPhi1.size() << endl
0087 << "Number of DTRecHit1DPair in the SL 3 (Phi 2): " << hitsFromPhi2.size() << endl;
0088
0089 theHitsFromPhi1 = hitsFromPhi1;
0090 theHitsFromPhi2 = hitsFromPhi2;
0091
0092 if (allDTRecHits) {
0093 DTRecHitCollection::range rangeHitsFromTheta =
0094 all1DHits->get(DTRangeMapAccessor::layersBySuperLayer(DTSuperLayerId(theChamber->id(), 2)));
0095
0096 vector<DTRecHit1DPair> hitsFromTheta(rangeHitsFromTheta.first, rangeHitsFromTheta.second);
0097 if (debug)
0098 cout << "Number of DTRecHit1DPair in the SL 2 (Theta): " << hitsFromTheta.size() << endl;
0099 theHitsFromTheta = hitsFromTheta;
0100 }
0101 }
0102
0103 void DTCombinatorialPatternReco4D::setDTRecSegment2DContainer(Handle<DTRecSegment2DCollection> all2DSegments) {
0104 theSegments2DTheta.clear();
0105
0106 if (!allDTRecHits) {
0107
0108 DTRecSegment2DCollection::range rangeTheta = all2DSegments->get(DTSuperLayerId(theChamber->id(), 2));
0109
0110
0111 vector<DTSLRecSegment2D> segments2DTheta(rangeTheta.first, rangeTheta.second);
0112
0113 if (debug)
0114 cout << "Number of 2D-segments in the second SL (Theta): " << segments2DTheta.size() << endl;
0115 theSegments2DTheta = segments2DTheta;
0116 }
0117 }
0118
0119 OwnVector<DTRecSegment4D> DTCombinatorialPatternReco4D::reconstruct() {
0120 OwnVector<DTRecSegment4D> result;
0121
0122 if (debug) {
0123 cout << "Segments in " << theChamber->id() << endl;
0124 cout << "Reconstructing of the Phi segments" << endl;
0125 }
0126
0127 vector<std::shared_ptr<DTHitPairForFit>> pairPhiOwned;
0128 vector<DTSegmentCand*> resultPhi = buildPhiSuperSegmentsCandidates(pairPhiOwned);
0129
0130 if (debug)
0131 cout << "There are " << resultPhi.size() << " Phi cand" << endl;
0132
0133 if (allDTRecHits) {
0134
0135 const DTSuperLayer* sl = theChamber->superLayer(2);
0136
0137 if (sl) {
0138
0139 if (debug)
0140 cout << "Reconstructing of the Theta segments" << endl;
0141 OwnVector<DTSLRecSegment2D> thetaSegs = the2DAlgo->reconstruct(sl, theHitsFromTheta);
0142 vector<DTSLRecSegment2D> segments2DTheta(thetaSegs.begin(), thetaSegs.end());
0143 theSegments2DTheta = segments2DTheta;
0144 }
0145 }
0146
0147 bool hasZed = false;
0148
0149
0150 if (!theSegments2DTheta.empty()) {
0151 hasZed = !theSegments2DTheta.empty();
0152 if (debug)
0153 cout << "There are " << theSegments2DTheta.size() << " Theta cand" << endl;
0154 } else {
0155 if (debug)
0156 cout << "No Theta SL" << endl;
0157 }
0158
0159
0160 if (debug)
0161 cout << "Building of the concrete DTRecSegment4D" << endl;
0162 if (!resultPhi.empty()) {
0163 for (vector<DTSegmentCand*>::const_iterator phi = resultPhi.begin(); phi != resultPhi.end(); ++phi) {
0164 std::unique_ptr<DTChamberRecSegment2D> superPhi(**phi);
0165
0166 theUpdator->update(superPhi.get(), false);
0167 if (debug)
0168 cout << "superPhi: " << *superPhi << endl;
0169
0170 if (hasZed) {
0171
0172
0173 for (vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin(); zed != theSegments2DTheta.end();
0174 ++zed) {
0175 if (debug)
0176 cout << "Theta: " << *zed << endl;
0177
0178 DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
0179
0180
0181
0182 const DTSuperLayer* zSL = theChamber->superLayer(ZedSegSLId);
0183
0184
0185 LocalPoint zPos(
0186 zed->localPosition().x(), (zSL->toLocal(theChamber->toGlobal(superPhi->localPosition()))).y(), 0.);
0187
0188 const LocalPoint posZInCh = theChamber->toLocal(zSL->toGlobal(zPos));
0189
0190 const LocalVector dirZInCh = theChamber->toLocal(zSL->toGlobal(zed->localDirection()));
0191
0192 DTRecSegment4D* newSeg = new DTRecSegment4D(*superPhi, *zed, posZInCh, dirZInCh);
0193
0194 if (debug)
0195 cout << "Created a 4D seg " << *newSeg << endl;
0196
0197
0198 theUpdator->update(newSeg, false, false);
0199 if (debug)
0200 cout << " seg updated " << *newSeg << endl;
0201
0202 if (!applyT0corr && computeT0corr)
0203 theUpdator->calculateT0corr(newSeg);
0204 if (applyT0corr)
0205 theUpdator->update(newSeg, true, false);
0206
0207 result.push_back(newSeg);
0208 }
0209 } else {
0210
0211
0212 DTRecSegment4D* newSeg = new DTRecSegment4D(*superPhi);
0213
0214 if (debug)
0215 cout << "Created a 4D segment using only the 2D Phi segment " << *newSeg << endl;
0216
0217
0218 if (!applyT0corr && computeT0corr)
0219 theUpdator->calculateT0corr(newSeg);
0220 if (applyT0corr)
0221 theUpdator->update(newSeg, true, false);
0222
0223 result.push_back(newSeg);
0224 }
0225 }
0226 } else {
0227
0228 if (hasZed) {
0229 for (vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin(); zed != theSegments2DTheta.end();
0230 ++zed) {
0231 if (debug)
0232 cout << "Theta: " << *zed << endl;
0233
0234
0235 DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
0236
0237 const LocalPoint posZInCh =
0238 theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localPosition()));
0239 const LocalVector dirZInCh =
0240 theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localDirection()));
0241
0242 DTRecSegment4D* newSeg = new DTRecSegment4D(*zed, posZInCh, dirZInCh);
0243
0244 if (debug)
0245 cout << "Created a 4D segment using only the 2D Theta segment " << *newSeg << endl;
0246
0247 if (!applyT0corr && computeT0corr)
0248 theUpdator->calculateT0corr(newSeg);
0249 if (applyT0corr)
0250 theUpdator->update(newSeg, true, false);
0251
0252 result.push_back(newSeg);
0253 }
0254 }
0255 }
0256
0257 for (vector<DTSegmentCand*>::iterator phi = resultPhi.begin(); phi != resultPhi.end(); ++phi)
0258 delete *phi;
0259
0260 return result;
0261 }
0262
0263 vector<DTSegmentCand*> DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates(
0264 vector<std::shared_ptr<DTHitPairForFit>>& pairPhiOwned) {
0265 DTSuperLayerId slId;
0266
0267 if (!theHitsFromPhi1.empty())
0268 slId = theHitsFromPhi1.front().wireId().superlayerId();
0269 else if (!theHitsFromPhi2.empty())
0270 slId = theHitsFromPhi2.front().wireId().superlayerId();
0271 else {
0272 if (debug)
0273 cout << "DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates: "
0274 << "No Hits in the two Phi SL" << endl;
0275 return vector<DTSegmentCand*>();
0276 }
0277
0278 const DTSuperLayer* sl = theDTGeometry->superLayer(slId);
0279
0280 vector<std::shared_ptr<DTHitPairForFit>> pairPhi1 = the2DAlgo->initHits(sl, theHitsFromPhi1);
0281
0282 vector<std::shared_ptr<DTHitPairForFit>> pairPhi2 = the2DAlgo->initHits(sl, theHitsFromPhi2);
0283
0284 copy(pairPhi2.begin(), pairPhi2.end(), back_inserter(pairPhi1));
0285
0286 pairPhiOwned.swap(pairPhi1);
0287
0288 return the2DAlgo->buildSegments(sl, pairPhiOwned);
0289 }
0290
0291
0292 DTRecSegment4D* DTCombinatorialPatternReco4D::segmentSpecialZed(const DTRecSegment4D* seg) {
0293
0294
0295 const DTSLRecSegment2D* zedSeg = seg->zSegment();
0296 std::vector<DTRecHit1D> hits = zedSeg->specificRecHits();
0297
0298
0299 int nHits = hits.size();
0300 DTRecHit1D middle = hits[static_cast<int>(nHits / 2.)];
0301
0302
0303 LocalPoint posInSL = zedSeg->localPosition();
0304 LocalVector dirInSL = zedSeg->localDirection();
0305 LocalPoint posInMiddleLayer = posInSL + dirInSL * (-posInSL.z()) / cos(dirInSL.theta());
0306
0307
0308 auto hit = std::make_unique<DTRecHit1D>(
0309 middle.wireId(), middle.lrSide(), middle.digiTime(), posInMiddleLayer, zedSeg->localPositionError());
0310
0311 std::vector<DTRecHit1D> newHits(1, *hit);
0312
0313
0314 LocalPoint pos(zedSeg->localPosition());
0315 LocalVector dir(zedSeg->localDirection());
0316 AlgebraicSymMatrix cov(zedSeg->covMatrix());
0317 double chi2(zedSeg->chi2());
0318
0319 auto newZed = std::make_unique<DTSLRecSegment2D>(zedSeg->superLayerId(), pos, dir, cov, chi2, newHits);
0320
0321
0322
0323 DTRecSegment4D* result = new DTRecSegment4D(*seg->phiSegment(), *newZed, seg->localPosition(), seg->localDirection());
0324
0325 delete seg;
0326
0327
0328 return result;
0329 }