MuonSeedGenerator/src/RPCCosmicSeedrecHitFinder.cc

0001 /**

0002  *  See header file for a description of this class.

0003  *

0004  */
0005
0006 #include "RecoMuon/MuonSeedGenerator/src/RPCCosmicSeedrecHitFinder.h"
0007 #include <DataFormats/TrackingRecHit/interface/TrackingRecHit.h>
0008 #include <DataFormats/RPCRecHit/interface/RPCRecHit.h>
0009 #include <FWCore/Framework/interface/ESHandle.h>
0010 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0011
0012 using namespace std;
0013 using namespace edm;
0014
0015 MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator find(
0016     MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator firstIter,
0017     MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator lastIter,
0018     const MuonTransientTrackingRecHit::MuonRecHitPointer& recHitRef) {
0019   MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator index = lastIter;
0020   for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iter = firstIter; iter != lastIter; iter++)
0021     if ((*iter) == recHitRef)
0022       index = iter;
0023   return index;
0024 }
0025
0026 RPCCosmicSeedrecHitFinder::RPCCosmicSeedrecHitFinder() {
0027   // Initiate the member

0028   isLayerset = false;
0029   isConfigured = false;
0030   isInputset = false;
0031   isOutputset = false;
0032   isEdgeset = false;
0033   BxRange = 0;
0034   MaxDeltaPhi = 0;
0035   ClusterSet.clear();
0036   innerBounds.clear();
0037   isLayersmixed = false;
0038   LayersinRPC.clear();
0039   therecHits.clear();
0040   isOuterLayerfilled = false;
0041 }
0042
0043 RPCCosmicSeedrecHitFinder::~RPCCosmicSeedrecHitFinder() {}
0044
0045 void RPCCosmicSeedrecHitFinder::configure(const edm::ParameterSet& iConfig) {
0046   // Set the configuration

0047   BxRange = iConfig.getParameter<unsigned int>("BxRange");
0048   MaxDeltaPhi = iConfig.getParameter<double>("MaxDeltaPhi");
0049   ClusterSet = iConfig.getParameter<std::vector<int> >("ClusterSet");
0050
0051   // Set the signal open

0052   isConfigured = true;
0053 }
0054
0055 void RPCCosmicSeedrecHitFinder::setInput(MuonRecHitContainer (&recHits)[RPCLayerNumber]) {
0056   for (unsigned int i = 0; i < RPCLayerNumber; i++) {
0057     AllrecHits[i].clear();
0058     AllrecHits[i] = recHits[i];
0059   }
0060
0061   // Set the signal open

0062   isInputset = true;
0063 }
0064
0065 void RPCCosmicSeedrecHitFinder::setEdge(const RPCGeometry& rpcGeometry) {
0066   // Find all chamber in RB1in and collect their surface

0067   const std::vector<DetId>& AllRPCId = rpcGeometry.detIds();
0068   for (std::vector<DetId>::const_iterator it = AllRPCId.begin(); it != AllRPCId.end(); it++) {
0069     RPCDetId RPCId(it->rawId());
0070     int Region = RPCId.region();
0071     int Station = RPCId.station();
0072     int Layer = RPCId.layer();
0073     if (Region == 0 && Station == 1 && Layer == 1) {
0074       const BoundPlane RPCChamberSurface = rpcGeometry.chamber(RPCId)->surface();
0075       innerBounds.push_back(RPCChamberSurface);
0076     }
0077   }
0078
0079   // Set the signal open

0080   isEdgeset = true;
0081 }
0082
0083 void RPCCosmicSeedrecHitFinder::unsetEdge() {
0084   // Clear all surfaces of chambers in RB1in

0085   innerBounds.clear();
0086   isEdgeset = false;
0087 }
0088
0089 void RPCCosmicSeedrecHitFinder::unsetInput() {
0090   for (unsigned int i = 0; i < RPCLayerNumber; i++)
0091     AllrecHits[i].clear();
0092   isInputset = false;
0093 }
0094
0095 void RPCCosmicSeedrecHitFinder::setOutput(RPCSeedFinder* Seed) {
0096   theSeed = Seed;
0097   // Set the signal open

0098   isOutputset = true;
0099 }
0100
0101 void RPCCosmicSeedrecHitFinder::setLayers(const std::vector<unsigned int>& Layers) {
0102   LayersinRPC = Layers;
0103   // Set the signal open

0104   isLayerset = true;
0105 }
0106
0107 void RPCCosmicSeedrecHitFinder::fillrecHits() {
0108   if (isLayerset == false || isConfigured == false || isOutputset == false || isInputset == false ||
0109       isEdgeset == false) {
0110     cout << "Not set the IO or not configured yet" << endl;
0111     return;
0112   }
0113
0114   therecHits.clear();
0115
0116   if (LayersinRPC.empty()) {
0117     cout << "Not set with any layers" << endl;
0118     LayersinRPC.clear();
0119     therecHits.clear();
0120     isLayerset = false;
0121   }
0122
0123   // check the layers, 1=all barrel, 2=all endcap, 3=mix

0124   unsigned int Component = LayerComponent();
0125   if (Component == 3)
0126     isLayersmixed = true;
0127   else
0128     isLayersmixed = false;
0129
0130   GlobalVector initVector(0, 0, 0);
0131   const MuonRecHitPointer recHitRef;
0132   isOuterLayerfilled = false;
0133   complete(initVector, recHitRef);
0134
0135   // Unset the signal

0136   LayersinRPC.clear();
0137   therecHits.clear();
0138   isLayerset = false;
0139 }
0140
0141 int RPCCosmicSeedrecHitFinder::LayerComponent() {
0142   bool isBarrel = false;
0143   bool isEndcap = false;
0144   for (std::vector<unsigned int>::const_iterator it = LayersinRPC.begin(); it != LayersinRPC.end(); it++) {
0145     if ((*it) < BarrelLayerNumber)
0146       isBarrel = true;
0147     if ((*it) >= BarrelLayerNumber && (*it) < (BarrelLayerNumber + EachEndcapLayerNumber * 2))
0148       isEndcap = true;
0149   }
0150   if (isBarrel == true && isEndcap == true)
0151     return 3;
0152   if (isEndcap == true)
0153     return 2;
0154   if (isBarrel == true)
0155     return 1;
0156   return 0;
0157 }
0158
0159 bool RPCCosmicSeedrecHitFinder::complete(const GlobalVector& lastSegment, const MuonRecHitPointer& lastrecHitRef) {
0160   bool isrecHitsfound = false;
0161
0162   for (unsigned int i = 0; i < RPCLayerNumber; i++)
0163     for (MuonRecHitContainer::const_iterator it = AllrecHits[i].begin(); it != AllrecHits[i].end(); it++) {
0164       cout << "Finding recHits from " << i << " th layer" << endl;
0165       // information for recHits

0166       GlobalPoint currentPosition = (*it)->globalPosition();
0167       int currentBX;
0168
0169       // Check validation

0170       if (!(*it)->isValid())
0171         continue;
0172
0173       // Check BX range, be sure there is only RPCRecHit in the MuonRecHitContainer when use the dynamic_cast

0174       TrackingRecHit* thisTrackingRecHit = (*it)->hit()->clone();
0175       // Should also delete the RPCRecHit object cast by dynamic_cast<> ?

0176       RPCRecHit* thisRPCRecHit = dynamic_cast<RPCRecHit*>(thisTrackingRecHit);
0177       currentBX = thisRPCRecHit->BunchX();
0178       int ClusterSize = thisRPCRecHit->clusterSize();
0179       delete thisTrackingRecHit;
0180       // Check BX

0181       if ((unsigned int)abs(currentBX) > BxRange)
0182         continue;
0183
0184       // Check cluster size

0185       bool Clustercheck = false;
0186       if (ClusterSet.empty())
0187         Clustercheck = true;
0188       for (std::vector<int>::const_iterator CluIter = ClusterSet.begin(); CluIter != ClusterSet.end(); CluIter++)
0189         if (ClusterSize == (*CluIter))
0190           Clustercheck = true;
0191       if (Clustercheck != true)
0192         continue;
0193
0194       cout << "Candidate recHit's position: " << currentPosition.x() << ", " << currentPosition.y() << ", "
0195            << currentPosition.z() << ". BX : " << currentBX << endl;
0196       // Fill 1st recHit from outer layers and rest recHits from all layers

0197       if (isOuterLayerfilled == false) {
0198         // Pick out the recHit from outer layers and fill it

0199         if (!isouterLayer(*it))
0200           continue;
0201
0202         // If pass all, add to the seed

0203         GlobalVector currentSegment = GlobalVector(0, 0, 0);
0204         cout << "1st recHit's global position: " << currentPosition.x() << ", " << currentPosition.y() << ", "
0205              << currentPosition.z() << ". BX: " << currentBX << endl;
0206         isrecHitsfound = true;
0207         therecHits.push_back(*it);
0208         isOuterLayerfilled = true;
0209         complete(currentSegment, *it);
0210         // Remember to pop the recHit before add another one from the same layer!

0211         therecHits.pop_back();
0212         isOuterLayerfilled = false;
0213       } else {
0214         GlobalPoint lastPosition = lastrecHitRef->globalPosition();
0215         TrackingRecHit* lastTrackingRecHit = lastrecHitRef->hit()->clone();
0216         // Should also delete the RPCRecHit object cast by dynamic_cast<> ?

0217         RPCRecHit* lastRPCRecHit = dynamic_cast<RPCRecHit*>(lastTrackingRecHit);
0218         int lastBX = lastRPCRecHit->BunchX();
0219         delete lastTrackingRecHit;
0220
0221         // Check the Y coordinate, shoule be lower than current one

0222         if (currentPosition.y() >= lastPosition.y())
0223           continue;
0224
0225         // Check the BX, should be larger than current one

0226         if (currentBX < lastBX)
0227           continue;
0228
0229         // If be the 2nd recHit, just fill it

0230         bool isinsideRegion = isinsideAngleRange(lastSegment, lastPosition, currentPosition);
0231         cout << "Check isinsideRegion: " << isinsideRegion << endl;
0232         if (!isinsideRegion)
0233           continue;
0234
0235         // If cross the edge the recHit should belong to another seed

0236         bool iscrossanyEdge = iscorssEdge(lastrecHitRef, *it);
0237         cout << "Check iscrossanyEdge: " << iscrossanyEdge << endl;
0238         if (iscrossanyEdge)
0239           continue;
0240
0241         // If pass all, add to the seed

0242         unsigned int NumberinSeed = therecHits.size();
0243         GlobalVector currentSegment = (GlobalVector)(currentPosition - lastPosition);
0244         cout << (NumberinSeed + 1) << "th recHit's global position: " << currentPosition.x() << ", "
0245              << currentPosition.y() << ", " << currentPosition.z() << ". BX: " << currentBX << endl;
0246         isrecHitsfound = true;
0247         therecHits.push_back(*it);
0248
0249         // if could not find next recHit in the search path, and have enough recHits already, that is the candidate

0250         bool findNext = complete(currentSegment, *it);
0251         if (findNext == false && therecHits.size() > 3) {
0252           for (ConstMuonRecHitContainer::const_iterator iter = therecHits.begin(); iter != therecHits.end(); iter++)
0253             cout << "Find recHit in seed candidate : " << (*iter)->globalPosition().x() << ", "
0254                  << (*iter)->globalPosition().y() << ", " << (*iter)->globalPosition().z() << endl;
0255           checkandfill();
0256         }
0257
0258         // Remember to pop the recHit before add another one from the same layer!

0259         therecHits.pop_back();
0260       }
0261     }
0262   return isrecHitsfound;
0263 }
0264
0265 bool RPCCosmicSeedrecHitFinder::isouterLayer(const MuonRecHitPointer& recHitRef) {
0266   bool isinsideLayers = false;
0267   for (std::vector<unsigned int>::const_iterator it = LayersinRPC.begin(); it != LayersinRPC.end(); it++) {
0268     MuonRecHitContainer::const_iterator index = find(AllrecHits[*it].begin(), AllrecHits[*it].end(), recHitRef);
0269     if (index != AllrecHits[*it].end())
0270       isinsideLayers = true;
0271   }
0272   return isinsideLayers;
0273 }
0274
0275 bool RPCCosmicSeedrecHitFinder::isinsideAngleRange(const GlobalVector& lastSegment,
0276                                                    const GlobalPoint& lastPosition,
0277                                                    const GlobalPoint& currentPosition) {
0278   bool isinsideAngle = true;
0279   GlobalVector SegVec = currentPosition - lastPosition;
0280   if (lastSegment.mag() != 0)
0281     if (fabs((lastSegment.phi() - SegVec.phi()).value()) > MaxDeltaPhi)
0282       isinsideAngle = false;
0283
0284   return isinsideAngle;
0285 }
0286
0287 bool RPCCosmicSeedrecHitFinder::iscorssEdge(const MuonRecHitPointer& lastrecHitRef,
0288                                             const MuonRecHitPointer& currentrecHitRef) {
0289   bool iscorss = false;
0290
0291   // Check if 2 recHits corss the inner bounds

0292   GlobalPoint lastPosition = lastrecHitRef->globalPosition();
0293   GlobalPoint currentPosition = currentrecHitRef->globalPosition();
0294   GlobalPoint testPosition((lastPosition.x() + currentPosition.x()) / 2,
0295                            (lastPosition.y() + currentPosition.y()) / 2,
0296                            (lastPosition.z() + currentPosition.z()) / 2);
0297   /*

0298     for(std::vector<BoundPlane>::const_iterator it = innerBounds.begin(); it != innerBounds.end(); it++) {

0299         //SurfaceOrientation::Side TestSide0 = it->side(currentPosition, 0);

0300         //SurfaceOrientation::Side TestSide1 = it->side(lastPosition, 0);

0301         SurfaceOrientation::Side TestSide = it->side(testPosition, 0);

0302         //cout << "Side of currentPosition: " << TestSide0 << ", Side of lastPosition: " << TestSide1 << ", Side of middlePosition: " << TestSide << endl;

0303         //if(TestSide != SurfaceOrientation::positiveSide)

0304             //iscorss = true;

0305     }

0306     */
0307
0308   // Check when mixLayer is not set

0309   if (isLayersmixed == false) {
0310     DetId lastId = lastrecHitRef->geographicalId();
0311     RPCDetId lastRPCId(lastId.rawId());
0312     int lastRegion = lastRPCId.region();
0313     DetId currentId = currentrecHitRef->geographicalId();
0314     RPCDetId currentRPCId(currentId.rawId());
0315     int currentRegion = currentRPCId.region();
0316     // Check if 2 recHits from different regions

0317     if (lastRegion != currentRegion)
0318       iscorss = true;
0319   }
0320
0321   return iscorss;
0322 }
0323
0324 void RPCCosmicSeedrecHitFinder::checkandfill() {
0325   if (therecHits.size() >= 3) {
0326     theSeed->setrecHits(therecHits);
0327     theSeed->seed();
0328   } else
0329     cout << "Layer less than 3, could not fill a RPCSeedFinder" << endl;
0330 }