SpecialSeedGenerators/src/GenericTripletGenerator.cc

0001 #include "RecoTracker/SpecialSeedGenerators/interface/GenericTripletGenerator.h"
0002 #include "RecoTracker/TkSeedGenerator/interface/FastCircle.h"
0003 typedef SeedingHitSet::ConstRecHitPointer SeedingHit;
0004
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include <map>
0008
0009 GenericTripletGenerator::GenericTripletGenerator(const edm::ParameterSet& conf, edm::ConsumesCollector& iC)
0010     : theSeedingLayerToken(iC.consumes<SeedingLayerSetsHits>(conf.getParameter<edm::InputTag>("LayerSrc"))) {
0011   edm::LogInfo("CtfSpecialSeedGenerator|GenericTripletGenerator") << "Constructing GenericTripletGenerator";
0012 }
0013
0014 const OrderedSeedingHits& GenericTripletGenerator::run(const TrackingRegion& region,
0015                                                        const edm::Event& e,
0016                                                        const edm::EventSetup& es) {
0017   hitTriplets.clear();
0018   hitTriplets.reserve(0);
0019   edm::Handle<SeedingLayerSetsHits> hlayers;
0020   e.getByToken(theSeedingLayerToken, hlayers);
0021   const SeedingLayerSetsHits& layers = *hlayers;
0022   if (layers.numberOfLayersInSet() != 3)
0023     throw cms::Exception("CtfSpecialSeedGenerator")
0024         << "You are using " << layers.numberOfLayersInSet() << " layers in set instead of 3 ";
0025   std::map<float, OrderedHitTriplet> radius_triplet_map;
0026   for (SeedingLayerSetsHits::SeedingLayerSet ls : layers) {
0027     auto innerHits = region.hits(ls[0]);
0028     //std::cout << "innerHits.size()=" << innerHits.size() << std::endl;
0029     auto middleHits = region.hits(ls[1]);
0030     //std::cout << "middleHits.size()=" << middleHits.size() << std::endl;
0031     auto outerHits = region.hits(ls[2]);
0032     //std::cout << "outerHits.size()=" << outerHits.size() << std::endl;
0033     //std::cout << "trying " << innerHits.size()*middleHits.size()*outerHits.size() << " combinations "<<std::endl;
0034     for (auto iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); iOuterHit++) {
0035       for (auto iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); iMiddleHit++) {
0036         for (auto iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); iInnerHit++) {
0037           //GlobalPoint innerpos  = ls[0].hitBuilder()->build(&(**iInnerHit))->globalPosition();
0038           //GlobalPoint middlepos = ls[1].hitBuilder()->build(&(**iMiddleHit))->globalPosition();
0039           //GlobalPoint outerpos  = ls[2].hitBuilder()->build(&(**iOuterHit))->globalPosition();
0040           //FastCircle circle(innerpos,
0041           //          middlepos,
0042           //                  outerpos);
0043           //do a first check on the radius of curvature to reduce the combinatorics
0044           OrderedHitTriplet oht(&(**iInnerHit), &(**iMiddleHit), &(**iOuterHit));
0045           std::pair<bool, float> val_radius = qualityFilter(oht, radius_triplet_map, ls);
0046           if (val_radius.first) {
0047             //if (circle.rho() > 200 || circle.rho() == 0) { //0 radius means straight line
0048             //hitTriplets.push_back(OrderedHitTriplet(*iInnerHit,
0049             //                  *iMiddleHit,
0050             //                  *iOuterHit));
0051             hitTriplets.push_back(oht);
0052             radius_triplet_map.insert(std::make_pair(val_radius.second, oht));
0053           }
0054         }
0055       }
0056     }
0057   }
0058   //std::cout << "ending with " << hitTriplets.size() << " triplets" << std::endl;
0059   return hitTriplets;
0060 }
0061
0062 std::pair<bool, float> GenericTripletGenerator::qualityFilter(const OrderedHitTriplet& oht,
0063                                                               const std::map<float, OrderedHitTriplet>& map,
0064                                                               const SeedingLayerSetsHits::SeedingLayerSet& ls) const {
0065   //first check the radius
0066   GlobalPoint innerpos = oht.inner()->globalPosition();
0067   GlobalPoint middlepos = oht.middle()->globalPosition();
0068   GlobalPoint outerpos = oht.outer()->globalPosition();
0069   std::vector<const TrackingRecHit*> ohttrh;
0070   ohttrh.push_back(&(*(oht.inner())));
0071   ohttrh.push_back(&(*(oht.middle())));
0072   ohttrh.push_back(&(*(oht.outer())));
0073   std::vector<const TrackingRecHit*>::const_iterator ioht;
0074   //now chech that the change in phi is reasonable. the limiting angular distance is the one in case
0075   //one of the two points is a tangent point.
0076   float limit_phi_distance1 = sqrt((middlepos.x() - outerpos.x()) * (middlepos.x() - outerpos.x()) +
0077                                    (middlepos.y() - outerpos.y()) * (middlepos.y() - outerpos.y())) /
0078                               middlepos.mag();  //actually this is the tangent of the limiting angle
0079   float limit_phi_distance2 = sqrt((middlepos.x() - innerpos.x()) * (middlepos.x() - innerpos.x()) +
0080                                    (middlepos.y() - innerpos.y()) * (middlepos.y() - innerpos.y())) /
0081                               innerpos.mag();
0082   //if (fabs(tan(outerpos.phi()-middlepos.phi()))>limit_phi_distance1 ||
0083   //    fabs(tan(innerpos.phi()-middlepos.phi()))>limit_phi_distance2) {
0084   if (fabs(outerpos.phi() - middlepos.phi()) > fabs(atan(limit_phi_distance1)) ||
0085       fabs(innerpos.phi() - middlepos.phi()) > fabs(atan(limit_phi_distance2))) {
0086     //std::cout << "rejected because phi" << std::endl;
0087     return std::make_pair(false, 0.);
0088   }
0089   //should we add a control on the r-z plane too?
0090   /*
0091     //now check that there is no big change in the r-z projection
0092     float dz1 = outerpos.z()-middlepos.z();
0093     float dr1 = sqrt(outerpos.x()*outerpos.x()+outerpos.y()*outerpos.y())-
0094     sqrt(middlepos.x()*middlepos.x()+middlepos.y()*middlepos.y());
0095     float dz2 = middlepos.z()-innerpos.z();
0096     float dr2 = sqrt(middlepos.x()*middlepos.x()+middlepos.y()*middlepos.y())-
0097     sqrt(innerpos.x()*innerpos.x()+innerpos.y()*innerpos.y());
0098     float tan1 = dz1/dr1;
0099     float tan2 = dz2/dr2;
0100     //how much should we allow? should we make it configurable?
0101     if (fabs(tan1-tan2)/tan1>0.5){
0102     //std::cout << "rejected because z" << std::endl;
0103     return std::make_pair(false, 0.);
0104     }
0105     */
0106   //now check the radius is not too small
0107   FastCircle circle(innerpos, middlepos, outerpos);
0108   if (circle.rho() < 200 && circle.rho() != 0)
0109     return std::make_pair(false, circle.rho());  //to small radius
0110   //now check if at least 2 hits are shared with an existing triplet
0111   //look for similar radii in the map
0112   std::map<float, OrderedHitTriplet>::const_iterator lower_bound = map.lower_bound((1 - 0.01) * circle.rho());
0113   std::map<float, OrderedHitTriplet>::const_iterator upper_bound = map.upper_bound((1 + 0.01) * circle.rho());
0114   std::map<float, OrderedHitTriplet>::const_iterator iter;
0115   for (iter = lower_bound; iter != upper_bound && iter->first <= upper_bound->first; iter++) {
0116     int shared = 0;
0117     std::vector<const TrackingRecHit*> curtrh;
0118     curtrh.push_back(&*(iter->second.inner()));
0119     curtrh.push_back(&*(iter->second.middle()));
0120     curtrh.push_back(&*(iter->second.outer()));
0121     std::vector<const TrackingRecHit*>::const_iterator curiter;
0122     for (curiter = curtrh.begin(); curiter != curtrh.end(); curiter++) {
0123       for (ioht = ohttrh.begin(); ioht != ohttrh.end(); ioht++) {
0124         if ((*ioht)->geographicalId() == (*curiter)->geographicalId() &&
0125             ((*ioht)->localPosition() - (*curiter)->localPosition()).mag() < 1e-5)
0126           shared++;
0127       }
0128     }
0129     if (shared > 1)
0130       return std::make_pair(false, circle.rho());
0131   }
0132
0133   return std::make_pair(true, circle.rho());
0134 }