MuonSeedGenerator/src/RPCSeedrecHitFinder.cc

0001 /**

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

0003  *

0004  */
0005
0006 #include "RecoMuon/MuonSeedGenerator/src/RPCSeedrecHitFinder.h"
0007 #include <DataFormats/TrackingRecHit/interface/TrackingRecHit.h>
0008 #include <DataFormats/RPCRecHit/interface/RPCRecHit.h>
0009
0010 using namespace std;
0011 using namespace edm;
0012
0013 // Comparator function must be global function?

0014 // Could not be included in .h file, or there will be 2 lessPhi() functions in both RPCSeedGenerator and RPCSeedFinder module

0015 bool lessPhi(const MuonTransientTrackingRecHit::ConstMuonRecHitPointer& it1,
0016              const MuonTransientTrackingRecHit::ConstMuonRecHitPointer& it2) {
0017   // Don't need to use value() in Geom::Phi to short

0018   return (it1->globalPosition().phi() < it2->globalPosition().phi());
0019 }
0020
0021 RPCSeedrecHitFinder::RPCSeedrecHitFinder() {
0022   // Initiate the member

0023   isLayerset = false;
0024   isConfigured = false;
0025   isInputset = false;
0026   isOutputset = false;
0027   BxRange = 0;
0028   MaxDeltaPhi = 0;
0029   ClusterSet.clear();
0030   LayersinRPC.clear();
0031   therecHits.clear();
0032 }
0033
0034 RPCSeedrecHitFinder::~RPCSeedrecHitFinder() {}
0035
0036 void RPCSeedrecHitFinder::configure(const edm::ParameterSet& iConfig) {
0037   // Set the configuration

0038   BxRange = iConfig.getParameter<unsigned int>("BxRange");
0039   MaxDeltaPhi = iConfig.getParameter<double>("MaxDeltaPhi");
0040   ClusterSet = iConfig.getParameter<std::vector<int> >("ClusterSet");
0041
0042   // Set the signal open

0043   isConfigured = true;
0044 }
0045
0046 void RPCSeedrecHitFinder::setInput(MuonRecHitContainer (&recHits)[RPCLayerNumber]) {
0047   for (unsigned int i = 0; i < RPCLayerNumber; i++)
0048     recHitsRPC[i] = &recHits[i];
0049   isInputset = true;
0050 }
0051
0052 void RPCSeedrecHitFinder::unsetInput() { isInputset = false; }
0053 void RPCSeedrecHitFinder::setOutput(RPCSeedFinder* Seed) {
0054   theSeed = Seed;
0055   isOutputset = true;
0056 }
0057
0058 void RPCSeedrecHitFinder::setLayers(const std::vector<unsigned int>& Layers) {
0059   LayersinRPC = Layers;
0060   isLayerset = true;
0061 }
0062
0063 void RPCSeedrecHitFinder::fillrecHits() {
0064   if (isLayerset == false || isConfigured == false || isOutputset == false || isInputset == false) {
0065     cout << "Not set the IO or not configured yet" << endl;
0066     return;
0067   }
0068   cout << "Now fill recHits from Layers: ";
0069   for (unsigned int k = 0; k < LayersinRPC.size(); k++)
0070     cout << LayersinRPC[k] << " ";
0071   cout << endl;
0072   unsigned int LayerIndex = 0;
0073   therecHits.clear();
0074   complete(LayerIndex);
0075
0076   // Unset the signal

0077   LayersinRPC.clear();
0078   therecHits.clear();
0079   isLayerset = false;
0080 }
0081
0082 void RPCSeedrecHitFinder::complete(unsigned int LayerIndex) {
0083   for (MuonRecHitContainer::const_iterator it = recHitsRPC[LayersinRPC[LayerIndex]]->begin();
0084        it != recHitsRPC[LayersinRPC[LayerIndex]]->end();
0085        it++) {
0086     cout << "Completing layer[" << LayersinRPC[LayerIndex] << "]." << endl;
0087
0088     // Check validation

0089     if (!(*it)->isValid())
0090       continue;
0091
0092     // Check BX range, be sure there is only RPCRecHit in the MuonRecHitContainer when use the dynamic_cast

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

0095     RPCRecHit* thisRPCRecHit = dynamic_cast<RPCRecHit*>(thisTrackingRecHit);
0096     int BX = thisRPCRecHit->BunchX();
0097     int ClusterSize = thisRPCRecHit->clusterSize();
0098     delete thisTrackingRecHit;
0099     // Check BX

0100     if ((unsigned int)abs(BX) > BxRange)
0101       continue;
0102     // Check cluster size

0103     bool Clustercheck = false;
0104     if (ClusterSet.empty())
0105       Clustercheck = true;
0106     for (std::vector<int>::const_iterator CluIter = ClusterSet.begin(); CluIter != ClusterSet.end(); CluIter++)
0107       if (ClusterSize == (*CluIter))
0108         Clustercheck = true;
0109     if (Clustercheck != true)
0110       continue;
0111     // Check the recHits Phi range

0112     GlobalPoint pos = (*it)->globalPosition();
0113     double Phi = pos.phi();
0114     cout << "Phi: " << Phi << endl;
0115     // The recHits should locate in some phi range

0116     therecHits.push_back(*it);
0117     double deltaPhi = getdeltaPhifromrecHits();
0118     cout << "Delta phi: " << deltaPhi << endl;
0119     therecHits.pop_back();
0120     if (deltaPhi > MaxDeltaPhi)
0121       continue;
0122
0123     // If pass all, add to the seed

0124     therecHits.push_back(*it);
0125     cout << "RecHit's global position: " << pos.x() << ", " << pos.y() << ", " << pos.z() << endl;
0126
0127     // Check if this recHit is the last one in the seed

0128     // If it is the last one, calculate the seed

0129     if (LayerIndex == (LayersinRPC.size() - 1)) {
0130       cout << "Check and fill one seed." << endl;
0131       checkandfill();
0132     }
0133     // If it is not the last one, continue to fill the seed from other layers

0134     else
0135       complete(LayerIndex + 1);
0136
0137     // Remember to pop the recHit before add another one from the same layer!

0138     therecHits.pop_back();
0139   }
0140 }
0141
0142 double RPCSeedrecHitFinder::getdeltaPhifromrecHits() {
0143   ConstMuonRecHitContainer sortRecHits = therecHits;
0144   sort(sortRecHits.begin(), sortRecHits.end(), lessPhi);
0145   cout << "Sorted recHit's Phi: ";
0146   for (ConstMuonRecHitContainer::const_iterator iter = sortRecHits.begin(); iter != sortRecHits.end(); iter++)
0147     cout << (*iter)->globalPosition().phi() << ", ";
0148   cout << endl;
0149   // Calculate the deltaPhi, take care Geom::Phi always in range [-pi,pi)

0150   // In case of some deltaPhi larger then Pi, use value() in Geom::Phi to get the true value in radians of Phi, then do the calculation

0151   double deltaPhi = 0;
0152   if (sortRecHits.size() <= 1)
0153     return deltaPhi;
0154   if (sortRecHits.size() == 2) {
0155     ConstMuonRecHitContainer::const_iterator iter1 = sortRecHits.begin();
0156     ConstMuonRecHitContainer::const_iterator iter2 = sortRecHits.begin();
0157     iter2++;
0158     deltaPhi = (((*iter2)->globalPosition().phi().value() - (*iter1)->globalPosition().phi().value()) > M_PI)
0159                    ? (2 * M_PI - ((*iter2)->globalPosition().phi().value() - (*iter1)->globalPosition().phi().value()))
0160                    : ((*iter2)->globalPosition().phi().value() - (*iter1)->globalPosition().phi().value());
0161     return deltaPhi;
0162   } else {
0163     deltaPhi = 2 * M_PI;
0164     int n = 0;
0165     for (ConstMuonRecHitContainer::const_iterator iter = sortRecHits.begin(); iter != sortRecHits.end(); iter++) {
0166       cout << "Before this loop deltaPhi is " << deltaPhi << endl;
0167       n++;
0168       double deltaPhi_more = 0;
0169       double deltaPhi_less = 0;
0170       if (iter == sortRecHits.begin()) {
0171         cout << "Calculateing frist loop..." << endl;
0172         ConstMuonRecHitContainer::const_iterator iter_more = ++iter;
0173         --iter;
0174         ConstMuonRecHitContainer::const_iterator iter_less = sortRecHits.end();
0175         --iter_less;
0176         cout << "more_Phi: " << (*iter_more)->globalPosition().phi()
0177              << ", less_Phi: " << (*iter_less)->globalPosition().phi()
0178              << ", iter_Phi: " << (*iter)->globalPosition().phi() << endl;
0179         deltaPhi_more =
0180             (2 * M_PI) - ((*iter_more)->globalPosition().phi().value() - (*iter)->globalPosition().phi().value());
0181         deltaPhi_less = (*iter_less)->globalPosition().phi().value() - (*iter)->globalPosition().phi().value();
0182       } else if (iter == (--sortRecHits.end())) {
0183         cout << "Calculateing last loop..." << endl;
0184         ConstMuonRecHitContainer::const_iterator iter_less = --iter;
0185         ++iter;
0186         ConstMuonRecHitContainer::const_iterator iter_more = sortRecHits.begin();
0187         cout << "more_Phi: " << (*iter_more)->globalPosition().phi()
0188              << ", less_Phi: " << (*iter_less)->globalPosition().phi()
0189              << ", iter_Phi: " << (*iter)->globalPosition().phi() << endl;
0190         deltaPhi_less =
0191             (2 * M_PI) - ((*iter)->globalPosition().phi().value() - (*iter_less)->globalPosition().phi().value());
0192         deltaPhi_more = (*iter)->globalPosition().phi().value() - (*iter_more)->globalPosition().phi().value();
0193       } else {
0194         cout << "Calculateing " << n << "st loop..." << endl;
0195         ConstMuonRecHitContainer::const_iterator iter_less = --iter;
0196         ++iter;
0197         ConstMuonRecHitContainer::const_iterator iter_more = ++iter;
0198         --iter;
0199         cout << "more_Phi: " << (*iter_more)->globalPosition().phi()
0200              << ", less_Phi: " << (*iter_less)->globalPosition().phi()
0201              << ", iter_Phi: " << (*iter)->globalPosition().phi() << endl;
0202         deltaPhi_less =
0203             (2 * M_PI) - ((*iter)->globalPosition().phi().value() - (*iter_less)->globalPosition().phi().value());
0204         deltaPhi_more =
0205             (2 * M_PI) - ((*iter_more)->globalPosition().phi().value() - (*iter)->globalPosition().phi().value());
0206       }
0207       if (deltaPhi > deltaPhi_more)
0208         deltaPhi = deltaPhi_more;
0209       if (deltaPhi > deltaPhi_less)
0210         deltaPhi = deltaPhi_less;
0211
0212       cout << "For this loop deltaPhi_more is " << deltaPhi_more << endl;
0213       cout << "For this loop deltaPhi_less is " << deltaPhi_less << endl;
0214       cout << "For this loop deltaPhi is " << deltaPhi << endl;
0215     }
0216     return deltaPhi;
0217   }
0218 }
0219
0220 void RPCSeedrecHitFinder::checkandfill() {
0221   if (therecHits.size() >= 3) {
0222     theSeed->setrecHits(therecHits);
0223     theSeed->seed();
0224   } else
0225     cout << "Layer less than 3, could not fill a RPCSeedFinder" << endl;
0226 }