File indexing completed on 2024-04-06 12:26:15
0001 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0002 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0003 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0004 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0005 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0006 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0007 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
0008 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
0009 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
0010 #include "DataFormats/DetId/interface/DetId.h"
0011 #include "FWCore/Framework/interface/ConsumesCollector.h"
0012 #include "TracktoRPC.h"
0013 #include "DTStationIndex.h"
0014 #include "DTObjectMap.h"
0015 #include "CSCStationIndex.h"
0016 #include "CSCObjectMap.h"
0017
0018 #include <ctime>
0019 #include <TMath.h>
0020
0021 bool TracktoRPC::ValidRPCSurface(RPCDetId rpcid, LocalPoint LocalP, const RPCGeometry *rpcGeo) {
0022 const GeomDet *whichdet3 = rpcGeo->idToDet(rpcid.rawId());
0023 const RPCRoll *aroll = dynamic_cast<const RPCRoll *>(whichdet3);
0024 float locx = LocalP.x(), locy = LocalP.y();
0025 if (aroll->isBarrel()) {
0026 const Bounds &rollbound = rpcGeo->idToDet((rpcid))->surface().bounds();
0027 float boundlength = rollbound.length();
0028 float boundwidth = rollbound.width();
0029
0030 if (fabs(locx) < boundwidth / 2 && fabs(locy) < boundlength / 2 && locy > -boundlength / 2)
0031 return true;
0032 else
0033 return false;
0034
0035 } else if (aroll->isForward()) {
0036 const Bounds &rollbound = rpcGeo->idToDet((rpcid))->surface().bounds();
0037 float boundlength = rollbound.length();
0038 float boundwidth = rollbound.width();
0039
0040 float nminx = TMath::Pi() * (18 * boundwidth / TMath::Pi() - boundlength) / 18;
0041 float ylimit = ((boundlength) / (boundwidth / 2 - nminx / 2)) * fabs(locx) + boundlength / 2 -
0042 ((boundlength) / (boundwidth / 2 - nminx / 2)) * (boundwidth / 2);
0043 if (ylimit < -boundlength / 2)
0044 ylimit = -boundlength / 2;
0045
0046 if (fabs(locx) < boundwidth / 2 && fabs(locy) < boundlength / 2 && locy > ylimit)
0047 return true;
0048 else
0049 return false;
0050 } else
0051 return false;
0052 }
0053
0054 TracktoRPC::TracktoRPC(const edm::ParameterSet &iConfig, const edm::InputTag &tracklabel, edm::ConsumesCollector iC)
0055 : rpcGeoToken_(iC.esConsumes()),
0056 dtGeoToken_(iC.esConsumes()),
0057 dtMapToken_(iC.esConsumes()),
0058 cscGeoToken_(iC.esConsumes()),
0059 cscMapToken_(iC.esConsumes()),
0060 propagatorToken_(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))) {
0061 if (tracklabel.label().find("cosmic") == 0)
0062 theTrackTransformer = std::make_unique<TrackTransformerForCosmicMuons>(iConfig, iC);
0063 else if (tracklabel.label().find("globalCosmic") == 0)
0064 theTrackTransformer = std::make_unique<TrackTransformerForCosmicMuons>(iConfig, iC);
0065 else
0066 theTrackTransformer = std::make_unique<TrackTransformer>(iConfig, iC);
0067 }
0068
0069 std::unique_ptr<RPCRecHitCollection> TracktoRPC::thePoints(reco::TrackCollection const *alltracks,
0070 edm::EventSetup const &iSetup,
0071 bool debug) {
0072 auto _ThePoints = std::make_unique<RPCRecHitCollection>();
0073
0074
0075 theTrackTransformer->setServices(iSetup);
0076
0077 const RPCGeometry *rpcGeo = &iSetup.getData(rpcGeoToken_);
0078 const DTGeometry *dtGeo = &iSetup.getData(dtGeoToken_);
0079 const DTObjectMap *dtMap = &iSetup.getData(dtMapToken_);
0080 const CSCGeometry *cscGeo = &iSetup.getData(cscGeoToken_);
0081 const CSCObjectMap *cscMap = &iSetup.getData(cscMapToken_);
0082 const Propagator *propagator = &iSetup.getData(propagatorToken_);
0083
0084 edm::OwnVector<RPCRecHit> RPCPointVector;
0085 std::vector<uint32_t> rpcput;
0086 double MaxD = 999.;
0087
0088 for (TrackCollection::const_iterator track = alltracks->begin(); track != alltracks->end(); track++) {
0089 Trajectories trajectories = theTrackTransformer->transform(*track);
0090 if (debug)
0091 std::cout << "Building Trajectory from Track. " << std::endl;
0092
0093 std::vector<uint32_t> rpcrolls;
0094 std::vector<uint32_t> rpcrolls2;
0095 std::map<uint32_t, int> rpcNdtcsc;
0096 std::map<uint32_t, int> rpcrollCounter;
0097
0098 float tInX = track->innerPosition().X(), tInY = track->innerPosition().Y(), tInZ = track->innerPosition().Z();
0099 float tOuX = track->outerPosition().X(), tOuY = track->outerPosition().Y(), tOuZ = track->outerPosition().Z();
0100 if (tInX > tOuX) {
0101 float temp = tOuX;
0102 tOuX = tInX;
0103 tInX = temp;
0104 }
0105 if (tInY > tOuY) {
0106 float temp = tOuY;
0107 tOuY = tInY;
0108 tInY = temp;
0109 }
0110 if (tInZ > tOuZ) {
0111 float temp = tOuZ;
0112 tOuZ = tInZ;
0113 tInZ = temp;
0114 }
0115
0116 if (debug)
0117 std::cout << "in (x,y,z): (" << tInX << ", " << tInY << ", " << tInZ << ")" << std::endl;
0118 if (debug)
0119 std::cout << "out (x,y,z): (" << tOuX << ", " << tOuY << ", " << tOuZ << ")" << std::endl;
0120
0121 if (debug)
0122 std::cout << "1. Search expeted RPC roll detid !!" << std::endl;
0123 for (trackingRecHit_iterator hit = track->recHitsBegin(); hit != track->recHitsEnd(); hit++) {
0124 if ((*hit)->isValid()) {
0125 DetId id = (*hit)->geographicalId();
0126
0127 if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::DT) {
0128 const GeomDet *geomDet = dtGeo->idToDet((*hit)->geographicalId());
0129 const DTLayer *dtlayer = dynamic_cast<const DTLayer *>(geomDet);
0130 if (dtlayer)
0131 for (Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end();
0132 ++trajectory) {
0133 const BoundPlane &DTSurface = dtlayer->surface();
0134 const GlobalPoint dcPoint = DTSurface.toGlobal(LocalPoint(0., 0., 0.));
0135
0136 TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
0137 const TrajectoryStateOnSurface &upd2 = (tMt).updatedState();
0138 if (upd2.isValid()) {
0139 LocalPoint trajLP = upd2.localPosition();
0140 LocalPoint trackLP = (*hit)->localPosition();
0141 float dx = trajLP.x() - trackLP.x(), dy = trajLP.y() - trackLP.y();
0142 if (dx > 10. && dy > 10.)
0143 continue;
0144
0145 DTChamberId dtid(geomDet->geographicalId().rawId());
0146 int dtW = dtid.wheel(), dtS = dtid.sector(), dtT = dtid.station();
0147 if (dtS == 13)
0148 dtS = 4;
0149 if (dtS == 14)
0150 dtS = 10;
0151 DTStationIndex theindex(0, dtW, dtS, dtT);
0152 const std::set<RPCDetId> &rollsForThisDT = dtMap->getRolls(theindex);
0153 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisDT.begin(); iteraRoll != rollsForThisDT.end();
0154 iteraRoll++) {
0155 const RPCRoll *rollasociated = rpcGeo->roll(*iteraRoll);
0156
0157 TrajectoryStateOnSurface ptss =
0158 propagator->propagate(upd2, rpcGeo->idToDet(rollasociated->id())->surface());
0159 if (ptss.isValid())
0160 if (ValidRPCSurface(rollasociated->id().rawId(), ptss.localPosition(), rpcGeo)) {
0161 rpcrollCounter[rollasociated->id().rawId()]++;
0162 bool check = true;
0163 std::vector<uint32_t>::iterator rpcroll;
0164 for (rpcroll = rpcrolls.begin(); rpcroll < rpcrolls.end(); rpcroll++)
0165 if (rollasociated->id().rawId() == *rpcroll)
0166 check = false;
0167 if (check == true) {
0168 rpcrolls.push_back(rollasociated->id().rawId());
0169 RPCGeomServ servId(rollasociated->id().rawId());
0170 if (debug)
0171 std::cout << "1\t Barrel RPC roll" << rollasociated->id().rawId() << " "
0172 << servId.name().c_str() << std::endl;
0173 }
0174 }
0175 }
0176 }
0177 }
0178 } else if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
0179 const GeomDet *geomDet = cscGeo->idToDet((*hit)->geographicalId());
0180 const CSCLayer *csclayer = dynamic_cast<const CSCLayer *>(geomDet);
0181
0182 CSCDetId cscid(geomDet->geographicalId().rawId());
0183 if (csclayer)
0184 for (Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end();
0185 ++trajectory) {
0186 const BoundPlane &CSCSurface = csclayer->surface();
0187 const GlobalPoint dcPoint = CSCSurface.toGlobal(LocalPoint(0., 0., 0.));
0188
0189 TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
0190 const TrajectoryStateOnSurface &upd2 = (tMt).updatedState();
0191
0192 if (upd2.isValid() && cscid.station() != 4 && cscid.ring() != 1) {
0193 LocalPoint trajLP = upd2.localPosition();
0194 LocalPoint trackLP = (*hit)->localPosition();
0195 float dx = trajLP.x() - trackLP.x(), dy = trajLP.y() - trackLP.y();
0196 if (dx > 10. && dy > 10.)
0197 continue;
0198
0199 int En = cscid.endcap(), St = cscid.station(), Ri = cscid.ring();
0200 int rpcSegment = cscid.chamber();
0201 if (En == 2)
0202 En = -1;
0203 if (Ri == 4)
0204 Ri = 1;
0205
0206 CSCStationIndex theindex(En, St, Ri, rpcSegment);
0207 const std::set<RPCDetId> &rollsForThisCSC = cscMap->getRolls(theindex);
0208 for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin();
0209 iteraRoll != rollsForThisCSC.end();
0210 iteraRoll++) {
0211 const RPCRoll *rollasociated = rpcGeo->roll(*iteraRoll);
0212
0213 TrajectoryStateOnSurface ptss =
0214 propagator->propagate(upd2, rpcGeo->idToDet(rollasociated->id())->surface());
0215 if (ptss.isValid())
0216 if (ValidRPCSurface(rollasociated->id().rawId(), ptss.localPosition(), rpcGeo)) {
0217 rpcrollCounter[rollasociated->id().rawId()]++;
0218 bool check = true;
0219 std::vector<uint32_t>::iterator rpcroll;
0220 for (rpcroll = rpcrolls.begin(); rpcroll < rpcrolls.end(); rpcroll++)
0221 if (rollasociated->id().rawId() == *rpcroll)
0222 check = false;
0223 if (check == true) {
0224 rpcrolls.push_back(rollasociated->id().rawId());
0225 RPCGeomServ servId(rollasociated->id().rawId());
0226 if (debug)
0227 std::cout << "1\t Forward RPC roll" << rollasociated->id().rawId() << " "
0228 << servId.name().c_str() << std::endl;
0229 }
0230 }
0231 }
0232 }
0233 }
0234 } else {
0235 if (debug)
0236 std::cout << "1\t The hit is not DT/CSC's. " << std::endl;
0237 }
0238 }
0239 }
0240 if (debug)
0241 std::cout << "First step OK!!\n2. Search nearest DT/CSC sufrace!!" << std::endl;
0242 std::vector<uint32_t>::iterator rpcroll;
0243 for (rpcroll = rpcrolls.begin(); rpcroll < rpcrolls.end(); rpcroll++) {
0244 RPCDetId rpcid(*rpcroll);
0245 const GlobalPoint &rGP = rpcGeo->idToDet(*rpcroll)->surface().toGlobal(LocalPoint(0, 0, 0));
0246 RPCGeomServ servId(rpcid);
0247 int rEn = rpcid.region(), rSe = rpcid.sector(), rWr = rpcid.ring(), rSt = rpcid.station(), rCh = servId.segment();
0248
0249 if (rpcrollCounter[*rpcroll] < 3)
0250 continue;
0251
0252 uint32_t dtcscid = 0;
0253 double distance = MaxD;
0254
0255
0256
0257
0258 for (trackingRecHit_iterator hit = track->recHitsBegin(); hit != track->recHitsEnd(); hit++) {
0259 if ((*hit)->isValid()) {
0260 DetId id = (*hit)->geographicalId();
0261 if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::DT) {
0262 const GeomDet *geomDet = dtGeo->idToDet((*hit)->geographicalId());
0263
0264 const GlobalPoint &dtGP = dtGeo->idToDet((*hit)->geographicalId())->surface().toGlobal(LocalPoint(0, 0, 0));
0265 double dx = rGP.x() - dtGP.x(), dy = rGP.y() - dtGP.y(), dz = rGP.z() - dtGP.z();
0266 double distanceN = sqrt(dx * dx + dy * dy + dz * dz);
0267
0268 DTChamberId dtid(geomDet->geographicalId().rawId());
0269 int Se = dtid.sector(), Wh = dtid.wheel(), St = dtid.station();
0270 if (Se == 13)
0271 Se = 4;
0272 if (Se == 14)
0273 Se = 10;
0274
0275 if (rEn == 0 && (rSe - Se) == 0 && (rWr - Wh) == 0 && (rSt - St) == 0 && distanceN < distance) {
0276 dtcscid = geomDet->geographicalId().rawId();
0277 distance = distanceN;
0278 if (debug)
0279 std::cout << "2\t DT " << dtcscid << " Wheel : " << Wh << " station : " << St << " sector : " << Se
0280 << std::endl;
0281 }
0282 } else if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
0283 const GeomDet *geomDet = cscGeo->idToDet((*hit)->geographicalId());
0284
0285 const GlobalPoint &cscGP =
0286 cscGeo->idToDet((*hit)->geographicalId())->surface().toGlobal(LocalPoint(0, 0, 0));
0287 double dx = rGP.x() - cscGP.x(), dy = rGP.y() - cscGP.y(), dz = rGP.z() - cscGP.z();
0288 double distanceN = sqrt(dx * dx + dy * dy + dz * dz);
0289
0290 CSCDetId cscid(geomDet->geographicalId().rawId());
0291 int En = cscid.endcap(), Ri = cscid.ring(), St = cscid.station(), Ch = cscid.chamber();
0292 if (En == 2)
0293 En = -1;
0294 if (Ri == 4)
0295 Ri = 1;
0296
0297 if ((rEn - En) == 0 && (rSt - St) == 0 && (Ch - rCh) == 0 && rWr != 1 && rSt != 4 && distanceN < distance) {
0298 dtcscid = geomDet->geographicalId().rawId();
0299 distance = distanceN;
0300 if (debug)
0301 std::cout << "2\t CSC " << dtcscid << " region : " << En << " station : " << St << " Ring : " << Ri
0302 << " chamber : " << Ch << std::endl;
0303 }
0304 }
0305 }
0306 }
0307 if (dtcscid != 0 && distance < MaxD) {
0308 rpcrolls2.push_back(*rpcroll);
0309 rpcNdtcsc[*rpcroll] = dtcscid;
0310 }
0311 }
0312 if (debug)
0313 std::cout << "Second step OK!! \n3. Propagate to RPC from DT/CSC!!" << std::endl;
0314
0315 std::vector<uint32_t>::iterator rpcroll2;
0316 for (rpcroll2 = rpcrolls2.begin(); rpcroll2 < rpcrolls2.end(); rpcroll2++) {
0317 bool check = true;
0318 std::vector<uint32_t>::iterator rpcput_;
0319 for (rpcput_ = rpcput.begin(); rpcput_ < rpcput.end(); rpcput_++)
0320 if (*rpcroll2 == *rpcput_)
0321 check = false;
0322
0323 if (check == true) {
0324 uint32_t dtcscid = rpcNdtcsc[*rpcroll2];
0325 DetId id(dtcscid);
0326 if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::DT) {
0327 const GeomDet *geomDet = dtGeo->idToDet(dtcscid);
0328 const DTLayer *dtlayer = dynamic_cast<const DTLayer *>(geomDet);
0329
0330 if (dtlayer)
0331 for (Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end();
0332 ++trajectory) {
0333 const BoundPlane &DTSurface = dtlayer->surface();
0334 const GlobalPoint dcPoint = DTSurface.toGlobal(LocalPoint(0., 0., 0.));
0335
0336 TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
0337 const TrajectoryStateOnSurface &upd2 = (tMt).updatedState();
0338 if (upd2.isValid()) {
0339 TrajectoryStateOnSurface ptss = propagator->propagate(upd2, rpcGeo->idToDet(*rpcroll2)->surface());
0340 if (ptss.isValid())
0341 if (ValidRPCSurface(*rpcroll2, ptss.localPosition(), rpcGeo)) {
0342 float rpcGPX = ptss.globalPosition().x();
0343 float rpcGPY = ptss.globalPosition().y();
0344 float rpcGPZ = ptss.globalPosition().z();
0345
0346 if (tInX > rpcGPX || tOuX < rpcGPX)
0347 continue;
0348 if (tInY > rpcGPY || tOuY < rpcGPY)
0349 continue;
0350 if (tInZ > rpcGPZ || tOuZ < rpcGPZ)
0351 continue;
0352
0353 const GeomDet *geomDet2 = rpcGeo->idToDet(*rpcroll2);
0354 const RPCRoll *aroll = dynamic_cast<const RPCRoll *>(geomDet2);
0355 const RectangularStripTopology *top_ =
0356 dynamic_cast<const RectangularStripTopology *>(&(aroll->topology()));
0357 LocalPoint xmin = top_->localPosition(0.);
0358 LocalPoint xmax = top_->localPosition((float)aroll->nstrips());
0359 float rsize = fabs(xmax.x() - xmin.x());
0360 float stripl = top_->stripLength();
0361
0362 float eyr = 1;
0363
0364 float locx = ptss.localPosition().x(), locy = ptss.localPosition().y(),
0365 locz = ptss.localPosition().z();
0366 if (locx < rsize * eyr && locy < stripl * eyr && locz < 1.) {
0367 RPCRecHit RPCPoint(*rpcroll2, 0, LocalPoint(locx, locy, locz));
0368
0369 RPCGeomServ servId(*rpcroll2);
0370 if (debug)
0371 std::cout << "3\t Barrel Expected RPC " << servId.name().c_str()
0372 << " \tLocalposition X: " << locx << ", Y: " << locy << " GlobalPosition(x,y,z) ("
0373 << rpcGPX << ", " << rpcGPY << ", " << rpcGPZ << ")" << std::endl;
0374 RPCPointVector.clear();
0375 RPCPointVector.push_back(RPCPoint);
0376 _ThePoints->put(*rpcroll2, RPCPointVector.begin(), RPCPointVector.end());
0377 rpcput.push_back(*rpcroll2);
0378 }
0379 }
0380 }
0381 }
0382 } else if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
0383 const GeomDet *geomDet4 = cscGeo->idToDet(dtcscid);
0384 const CSCLayer *csclayer = dynamic_cast<const CSCLayer *>(geomDet4);
0385
0386 if (csclayer)
0387 for (Trajectories::const_iterator trajectory = trajectories.begin(); trajectory != trajectories.end();
0388 ++trajectory) {
0389 const BoundPlane &CSCSurface = csclayer->surface();
0390 const GlobalPoint dcPoint = CSCSurface.toGlobal(LocalPoint(0., 0., 0.));
0391
0392 TrajectoryMeasurement tMt = trajectory->closestMeasurement(dcPoint);
0393 const TrajectoryStateOnSurface &upd2 = (tMt).updatedState();
0394 if (upd2.isValid()) {
0395 TrajectoryStateOnSurface ptss = propagator->propagate(upd2, rpcGeo->idToDet(*rpcroll2)->surface());
0396 if (ptss.isValid())
0397 if (ValidRPCSurface(*rpcroll2, ptss.localPosition(), rpcGeo)) {
0398 float rpcGPX = ptss.globalPosition().x();
0399 float rpcGPY = ptss.globalPosition().y();
0400 float rpcGPZ = ptss.globalPosition().z();
0401
0402 if (tInX > rpcGPX || tOuX < rpcGPX)
0403 continue;
0404 if (tInY > rpcGPY || tOuY < rpcGPY)
0405 continue;
0406 if (tInZ > rpcGPZ || tOuZ < rpcGPZ)
0407 continue;
0408
0409 RPCDetId rpcid(*rpcroll2);
0410 const GeomDet *geomDet3 = rpcGeo->idToDet(*rpcroll2);
0411 const RPCRoll *aroll = dynamic_cast<const RPCRoll *>(geomDet3);
0412 const TrapezoidalStripTopology *top_ =
0413 dynamic_cast<const TrapezoidalStripTopology *>(&(aroll->topology()));
0414 LocalPoint xmin = top_->localPosition(0.);
0415 LocalPoint xmax = top_->localPosition((float)aroll->nstrips());
0416 float rsize = fabs(xmax.x() - xmin.x());
0417 float stripl = top_->stripLength();
0418
0419
0420 float eyr = 1;
0421
0422 float locx = ptss.localPosition().x(), locy = ptss.localPosition().y(),
0423 locz = ptss.localPosition().z();
0424 if (locx < rsize * eyr && locy < stripl * eyr && locz < 1.) {
0425 RPCRecHit RPCPoint(*rpcroll2, 0, LocalPoint(locx, locy, locz));
0426 RPCGeomServ servId(*rpcroll2);
0427 if (debug)
0428 std::cout << "3\t Forward Expected RPC " << servId.name().c_str()
0429 << " \tLocalposition X: " << locx << ", Y: " << locy << " GlobalPosition(x,y,z) ("
0430 << rpcGPX << ", " << rpcGPY << ", " << rpcGPZ << ")" << std::endl;
0431 RPCPointVector.clear();
0432 RPCPointVector.push_back(RPCPoint);
0433 _ThePoints->put(*rpcroll2, RPCPointVector.begin(), RPCPointVector.end());
0434 rpcput.push_back(*rpcroll2);
0435 }
0436 }
0437 }
0438 }
0439 }
0440 }
0441 }
0442 if (debug)
0443 std::cout << "last steps OK!! " << std::endl;
0444 }
0445
0446 return _ThePoints;
0447 }