RPCRecHit/plugins/CSCSegtoRPC.cc

0001 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0002 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0003 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0004 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0005 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0006 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
0009 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
0010 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
0011 #include "CSCSegtoRPC.h"
0012 #include "CSCStationIndex.h"
0013 #include "CSCObjectMap.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 CSCSegtoRPC::CSCSegtoRPC(edm::ConsumesCollector iC, const edm::ParameterSet& iConfig)
0017     : rpcGeoToken_(iC.esConsumes()), cscGeoToken_(iC.esConsumes()), cscMapToken_(iC.esConsumes()) {
0018   minBX = iConfig.getParameter<int>("minBX");
0019   maxBX = iConfig.getParameter<int>("maxBX");
0020 }
0021
0022 std::unique_ptr<RPCRecHitCollection> CSCSegtoRPC::thePoints(const CSCSegmentCollection* allCSCSegments,
0023                                                             const edm::EventSetup& iSetup,
0024                                                             bool debug,
0025                                                             double eyr) {
0026   edm::ESHandle<RPCGeometry> rpcGeo = iSetup.getHandle(rpcGeoToken_);
0027   edm::ESHandle<CSCGeometry> cscGeo = iSetup.getHandle(cscGeoToken_);
0028   edm::ESHandle<CSCObjectMap> cscMap = iSetup.getHandle(cscMapToken_);
0029
0030   double MaxD = 80.;
0031
0032   if (debug)
0033     LogDebug("CSCSegtoRPC") << "CSC \t Number of CSC Segments in this event = " << allCSCSegments->size() << std::endl;
0034
0035   auto _ThePoints = std::make_unique<RPCRecHitCollection>();
0036   edm::OwnVector<RPCRecHit> RPCPointVector;
0037
0038   if (allCSCSegments->size() == 0) {
0039     if (debug)
0040       LogDebug("CSCSegtoRPC") << "CSC 0 segments skiping event" << std::endl;
0041   } else {
0042     std::map<CSCDetId, int> CSCSegmentsCounter;
0043     CSCSegmentCollection::const_iterator segment;
0044
0045     int segmentsInThisEventInTheEndcap = 0;
0046
0047     for (segment = allCSCSegments->begin(); segment != allCSCSegments->end(); ++segment) {
0048       CSCSegmentsCounter[segment->cscDetId()]++;
0049       segmentsInThisEventInTheEndcap++;
0050     }
0051
0052     float myTime = -9999.;
0053     float myTimeErr = -9999.;
0054     int myBx = -99;
0055     const float bunchCrossTimeDiff = 25.;  // time between bunch crossings
0056
0057     if (debug)
0058       LogDebug("CSCSegtoRPC") << "CSC \t loop over all the CSCSegments " << std::endl;
0059     for (segment = allCSCSegments->begin(); segment != allCSCSegments->end(); ++segment) {
0060       CSCDetId CSCId = segment->cscDetId();
0061
0062       myTime = segment->time();
0063       myBx = round(myTime / bunchCrossTimeDiff);
0064       if (!(myBx <= maxBX && myBx >= minBX))  // rpc read in bx in [-2, 2]
0065         continue;
0066
0067       if (debug)
0068         LogDebug("CSCSegtoRPC") << "CSC \t \t This Segment is in Chamber id: " << CSCId << std::endl;
0069       if (debug)
0070         LogDebug("CSCSegtoRPC") << "CSC \t \t Number of segments in this CSC = " << CSCSegmentsCounter[CSCId]
0071                                 << std::endl;
0072       if (debug)
0073         LogDebug("CSCSegtoRPC")
0074             << "CSC \t \t Is the only one in this CSC? is not ind the ring 1? Are there more than 2 "
0075                "segments in the event?"
0076             << std::endl;
0077
0078       if (CSCSegmentsCounter[CSCId] == 1 && CSCId.ring() != 1 && allCSCSegments->size() >= 2) {
0079         if (debug)
0080           LogDebug("CSCSegtoRPC") << "CSC \t \t yes" << std::endl;
0081         int cscEndCap = CSCId.endcap();
0082         int cscStation = CSCId.station();
0083         int cscRing = CSCId.ring();
0084         int rpcRegion = 1;
0085         if (cscEndCap == 2)
0086           rpcRegion = -1;  //Relacion entre las endcaps
0087         int rpcRing = cscRing;
0088         if (cscRing == 4)
0089           rpcRing = 1;
0090         int rpcStation = cscStation;
0091         int rpcSegment = CSCId.chamber();
0092
0093         LocalPoint segmentPosition = segment->localPosition();
0094         LocalVector segmentDirection = segment->localDirection();
0095         float dz = segmentDirection.z();
0096
0097         if (debug)
0098           LogDebug("CSCSegtoRPC") << "CSC \t \t \t Information about the segment"
0099                                   << "RecHits =" << segment->nRecHits() << "Angle =" << acos(dz) * 180 / 3.1415926
0100                                   << std::endl;
0101
0102         if (debug)
0103           LogDebug("CSCSegtoRPC")
0104               << "CSC \t \t Is a good Segment? dim = 4, 4 <= nRecHits <= 10 Incident angle int range 45 < "
0105               << acos(dz) * 180 / 3.1415926 << " < 135? " << std::endl;
0106
0107         if ((segment->dimension() == 4) && (segment->nRecHits() <= 10 && segment->nRecHits() >= 4)) {
0108           if (debug)
0109             LogDebug("CSCSegtoRPC") << "CSC \t \t yes" << std::endl;
0110           if (debug)
0111             LogDebug("CSCSegtoRPC") << "CSC \t \t CSC Segment Dimension " << segment->dimension() << std::endl;
0112
0113           float Xo = segmentPosition.x();
0114           float Yo = segmentPosition.y();
0115           float Zo = segmentPosition.z();
0116           float dx = segmentDirection.x();
0117           float dy = segmentDirection.y();
0118           float dz = segmentDirection.z();
0119
0120           if (debug)
0121             LogDebug("CSCSegtoRPC") << "Creating the CSCIndex" << std::endl;
0122           CSCStationIndex theindex(rpcRegion, rpcStation, rpcRing, rpcSegment);
0123           if (debug)
0124             LogDebug("CSCSegtoRPC") << "Getting the Rolls for the given index" << std::endl;
0125           std::set<RPCDetId> rollsForThisCSC = cscMap->getRolls(theindex);
0126           if (debug)
0127             LogDebug("CSCSegtoRPC") << "CSC \t \t Getting chamber from Geometry" << std::endl;
0128           const CSCChamber* TheChamber = cscGeo->chamber(CSCId);
0129           if (debug)
0130             LogDebug("CSCSegtoRPC") << "CSC \t \t Getting ID from Chamber" << std::endl;
0131           const CSCDetId TheId = TheChamber->id();
0132
0133           if (debug)
0134             LogDebug("CSCSegtoRPC") << "CSC \t \t Number of rolls for this CSC = " << rollsForThisCSC.size()
0135                                     << std::endl;
0136
0137           if (debug)
0138             LogDebug("CSCSegtoRPC") << "CSC \t \t Printing The Id" << TheId << std::endl;
0139
0140           if (rpcRing != 1) {  //They don't exist in Run3!
0141
0142             assert(!rollsForThisCSC.empty());
0143
0144             if (debug)
0145               LogDebug("CSCSegtoRPC") << "CSC \t \t Loop over all the rolls asociated to this CSC" << std::endl;
0146             for (std::set<RPCDetId>::iterator iteraRoll = rollsForThisCSC.begin(); iteraRoll != rollsForThisCSC.end();
0147                  iteraRoll++) {
0148               const RPCRoll* rollasociated = rpcGeo->roll(*iteraRoll);
0149               RPCDetId rpcId = rollasociated->id();
0150
0151               if (debug)
0152                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t We are in the roll getting the surface" << rpcId << std::endl;
0153               const BoundPlane& RPCSurface = rollasociated->surface();
0154
0155               if (debug)
0156                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t RollID: " << rpcId << std::endl;
0157
0158               if (debug)
0159                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Doing the extrapolation to this roll" << std::endl;
0160               if (debug)
0161                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t CSC Segment Direction in CSCLocal " << segmentDirection
0162                                         << std::endl;
0163               if (debug)
0164                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t CSC Segment Point in CSCLocal " << segmentPosition
0165                                         << std::endl;
0166
0167               GlobalPoint CenterPointRollGlobal = RPCSurface.toGlobal(LocalPoint(0, 0, 0));
0168               if (debug)
0169                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Center (0,0,0) of the Roll in Global" << CenterPointRollGlobal
0170                                         << std::endl;
0171               GlobalPoint CenterPointCSCGlobal = TheChamber->toGlobal(LocalPoint(0, 0, 0));
0172               if (debug)
0173                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Center (0,0,0) of the CSC in Global" << CenterPointCSCGlobal
0174                                         << std::endl;
0175               GlobalPoint segmentPositionInGlobal =
0176                   TheChamber->toGlobal(segmentPosition);  //new way to convert to global
0177               if (debug)
0178                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Segment Position in Global" << segmentPositionInGlobal
0179                                         << std::endl;
0180               LocalPoint CenterRollinCSCFrame = TheChamber->toLocal(CenterPointRollGlobal);
0181
0182               if (debug) {  //to check CSC RPC phi relation!
0183                 float rpcphi = 0;
0184                 float cscphi = 0;
0185
0186                 (CenterPointRollGlobal.barePhi() < 0) ? rpcphi = 2 * 3.141592 + CenterPointRollGlobal.barePhi()
0187                                                       : rpcphi = CenterPointRollGlobal.barePhi();
0188
0189                 (CenterPointCSCGlobal.barePhi() < 0) ? cscphi = 2 * 3.1415926536 + CenterPointCSCGlobal.barePhi()
0190                                                      : cscphi = CenterPointCSCGlobal.barePhi();
0191
0192                 float df = fabs(cscphi - rpcphi);
0193                 float dr = fabs(CenterPointRollGlobal.perp() - CenterPointCSCGlobal.perp());
0194                 float diffz = CenterPointRollGlobal.z() - CenterPointCSCGlobal.z();
0195                 float dfg = df * 180. / 3.14159265;
0196
0197                 if (debug)
0198                   LogDebug("CSCSegtoRPC") << "CSC \t \t \t z of RPC=" << CenterPointRollGlobal.z() << "z of CSC"
0199                                           << CenterPointCSCGlobal.z() << " dfg=" << dfg << std::endl;
0200
0201                 RPCGeomServ rpcsrv(rpcId);
0202
0203                 if (dr > 200. || fabs(dz) > 55. || dfg > 1.) {
0204                   if (debug)
0205                     LogDebug("CSCSegtoRPC")
0206                         << "\t \t \t CSC Station= " << CSCId.station() << " Ring= " << CSCId.ring()
0207                         << " Chamber= " << CSCId.chamber() << " cscphi=" << cscphi * 180 / 3.14159265
0208                         << "\t RPC Station= " << rpcId.station() << " ring= " << rpcId.ring() << " segment =-> "
0209                         << rpcsrv.segment() << " rollphi=" << rpcphi * 180 / 3.14159265 << "\t dfg=" << dfg
0210                         << " dz=" << diffz << " dr=" << dr << std::endl;
0211                 }
0212               }
0213
0214               float D = CenterRollinCSCFrame.z();
0215
0216               float X = Xo + dx * D / dz;
0217               float Y = Yo + dy * D / dz;
0218               float Z = D;
0219
0220               const TrapezoidalStripTopology* top_ =
0221                   dynamic_cast<const TrapezoidalStripTopology*>(&(rollasociated->topology()));
0222               LocalPoint xmin = top_->localPosition(0.);
0223               if (debug)
0224                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t xmin of this  Roll " << xmin << "cm" << std::endl;
0225               LocalPoint xmax = top_->localPosition((float)rollasociated->nstrips());
0226               if (debug)
0227                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t xmax of this  Roll " << xmax << "cm" << std::endl;
0228               float rsize = fabs(xmax.x() - xmin.x());
0229               if (debug)
0230                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Roll Size " << rsize << "cm" << std::endl;
0231               float stripl = top_->stripLength();
0232               float stripw = top_->pitch();
0233
0234               if (debug)
0235                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Strip Lenght " << stripl << "cm" << std::endl;
0236               if (debug)
0237                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Strip Width " << stripw << "cm" << std::endl;
0238
0239               if (debug)
0240                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t X Predicted in CSCLocal= " << X << "cm" << std::endl;
0241               if (debug)
0242                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Y Predicted in CSCLocal= " << Y << "cm" << std::endl;
0243               if (debug)
0244                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Z Predicted in CSCLocal= " << Z << "cm" << std::endl;
0245
0246               float extrapolatedDistance = sqrt((X - Xo) * (X - Xo) + (Y - Yo) * (Y - Yo) + (Z - Zo) * (Z - Zo));
0247
0248               if (debug)
0249                 LogDebug("CSCSegtoRPC") << "CSC \t \t \t Is the distance of extrapolation less than MaxD? ="
0250                                         << extrapolatedDistance << "cm"
0251                                         << " MaxD=" << MaxD << "cm" << std::endl;
0252
0253               if (extrapolatedDistance <= MaxD) {
0254                 if (debug)
0255                   LogDebug("CSCSegtoRPC") << "CSC \t \t \t yes" << std::endl;
0256
0257                 GlobalPoint GlobalPointExtrapolated = TheChamber->toGlobal(LocalPoint(X, Y, Z));
0258                 if (debug)
0259                   LogDebug("CSCSegtoRPC")
0260                       << "CSC \t \t \t Point ExtraPolated in Global" << GlobalPointExtrapolated << std::endl;
0261
0262                 LocalPoint PointExtrapolatedRPCFrame = RPCSurface.toLocal(GlobalPointExtrapolated);
0263
0264                 if (debug)
0265                   LogDebug("CSCSegtoRPC")
0266                       << "CSC \t \t \t Point Extrapolated in RPCLocal" << PointExtrapolatedRPCFrame << std::endl;
0267                 if (debug)
0268                   LogDebug("CSCSegtoRPC") << "CSC \t \t \t Corner of the Roll = (" << rsize * eyr << "," << stripl * eyr
0269                                           << ")" << std::endl;
0270                 if (debug)
0271                   LogDebug("CSCSegtoRPC")
0272                       << "CSC \t \t \t Info About the Point Extrapolated in X Abs ("
0273                       << fabs(PointExtrapolatedRPCFrame.x()) << "," << fabs(PointExtrapolatedRPCFrame.y()) << ","
0274                       << fabs(PointExtrapolatedRPCFrame.z()) << ")" << std::endl;
0275                 if (debug)
0276                   LogDebug("CSCSegtoRPC") << "CSC \t \t \t dz=" << fabs(PointExtrapolatedRPCFrame.z())
0277                                           << " dx=" << fabs(PointExtrapolatedRPCFrame.x())
0278                                           << " dy=" << fabs(PointExtrapolatedRPCFrame.y()) << std::endl;
0279
0280                 if (debug)
0281                   LogDebug("CSCSegtoRPC") << "CSC \t \t \t Does the extrapolation go inside this roll????" << std::endl;
0282
0283                 if (fabs(PointExtrapolatedRPCFrame.z()) < 1. && fabs(PointExtrapolatedRPCFrame.x()) < rsize * eyr &&
0284                     fabs(PointExtrapolatedRPCFrame.y()) < stripl * eyr) {
0285                   if (debug)
0286                     LogDebug("CSCSegtoRPC") << "CSC \t \t \t \t yes" << std::endl;
0287                   if (debug)
0288                     LogDebug("CSCSegtoRPC") << "CSC \t \t \t \t Creating the RecHit" << std::endl;
0289
0290                   RPCRecHit RPCPoint(rpcId, myBx, PointExtrapolatedRPCFrame);
0291                   RPCPoint.setTimeAndError(myTime, myTimeErr);
0292
0293                   if (debug)
0294                     LogDebug("CSCSegtoRPC") << "CSC \t \t \t \t Clearing the vector" << std::endl;
0295                   RPCPointVector.clear();
0296                   if (debug)
0297                     LogDebug("CSCSegtoRPC") << "CSC \t \t \t \t Pushing back" << std::endl;
0298                   RPCPointVector.push_back(RPCPoint);
0299                   if (debug)
0300                     LogDebug("CSCSegtoRPC") << "CSC \t \t \t \t Putting the vector" << std::endl;
0301                   _ThePoints->put(rpcId, RPCPointVector.begin(), RPCPointVector.end());
0302                 }
0303               }
0304             }
0305           }
0306         }
0307       }
0308     }
0309   }
0310   return _ThePoints;
0311 }