Merge pull request #12882 from VinInn/PhiWindows

Fix phi windows computation

CalibTracker/SiStripHitEfficiency/src/HitEff.cc

@@ -408,7 +408,7 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es){
 	  const TrajectoryStateOnSurface tsosTEC9 = itm->updatedState();
 
 	  // check if track on positive or negative z
-	  if (!iidd ==  StripSubdetector::TEC) cout << "there is a problem with TEC 9 extrapolation" << endl;
+	  if (DEBUG)  if (!(iidd ==  StripSubdetector::TEC)) cout << "there is a problem with TEC 9 extrapolation" << endl;
 
 	  //cout << " tec9 id = " << iidd << " and side = " << tTopo->tecSide(iidd) << endl;
 	  vector<TrajectoryMeasurement> tmp;

DataFormats/Math/interface/normalizedPhi.h

@@ -1,11 +1,30 @@
 #ifndef Math_notmalizedPhi_h
 #define Math_notmalizedPhi_h
 #include "DataFormats/Math/interface/deltaPhi.h"
-/* return a value of phi into interval [-pi,+pi]
- *
- */
+
+// return a value of phi into interval [-pi,+pi]
 template<typename T>
 inline
 T normalizedPhi(T phi) { return reco::reduceRange(phi);}
 
+// cernlib V306
+template<typename T>
+inline 
+T proxim(T b, T a) {
+        constexpr T c1 = 2.*M_PI;
+        constexpr T c2 = 1/c1;
+        return b+c1*std::round(c2*(a-b));
+}
+
+template<typename T>
+inline
+bool checkPhiInRange(T phi, T phi1, T phi2) {
+    constexpr T c1 = 2.*M_PI;
+    phi1 = normalizedPhi(phi1);
+    phi2 = proxim(phi2,phi1);
+    // phi & phi1 are in [-pi,pi] range...
+    return ( (phi1 <= phi) && (phi <= phi2) ) ||
+           ( (phi1 <= phi+c1) && (phi+c1 <= phi2) );
+}
+
 #endif

RecoEgamma/EgammaElectronProducers/python/gsfElectronCkfTrackCandidateMaker_cff.py

@@ -43,6 +43,7 @@
 gsfElectronChi2.ComponentName = 'gsfElectronChi2'
 gsfElectronChi2.MaxChi2 = 100000.
 gsfElectronChi2.nSigma = 3.
+gsfElectronChi2.MaxDispacement = 100
 gsfElectronChi2.MaxSagitta = -1
 
 TrajectoryFilterForPixelMatchGsfElectrons = cms.PSet(

RecoEgamma/EgammaPhotonProducers/python/chi2EstimatorForInOutFit_cfi.py

@@ -5,4 +5,5 @@
 Chi2MeasurementEstimatorForInOut.ComponentName = 'Chi2ForInOut'
 Chi2MeasurementEstimatorForInOut.MaxChi2 = 100.
 Chi2MeasurementEstimatorForInOut.nSigma = 3.
+Chi2MeasurementEstimatorForInOut.MaxDispacement = 100
 Chi2MeasurementEstimatorForInOut.MaxSagitta=-1

RecoEgamma/EgammaPhotonProducers/python/chi2EstimatorForOutInFit_cfi.py

@@ -5,5 +5,6 @@
 Chi2MeasurementEstimatorForOutIn.ComponentName = 'Chi2ForOutIn'
 Chi2MeasurementEstimatorForOutIn.MaxChi2 = 500.
 Chi2MeasurementEstimatorForOutIn.nSigma = 3.
+Chi2MeasurementEstimatorForOutIn.MaxDispacement = 100
 Chi2MeasurementEstimatorForOutIn.MaxSagitta = -1. 
 

RecoEgamma/EgammaPhotonProducers/python/looseChi2Estimator_cfi.py

@@ -5,4 +5,5 @@
 chi2CutForConversionTrajectoryBuilder.ComponentName = 'eleLooseChi2'
 chi2CutForConversionTrajectoryBuilder.MaxChi2 = 100000.
 chi2CutForConversionTrajectoryBuilder.nSigma = 3.
+chi2CutForConversionTrajectoryBuilder.MaxDispacement = 100.
 chi2CutForConversionTrajectoryBuilder.MaxSagitta = -1

RecoPixelVertexing/PixelTriplets/plugins/PixelTripletHLTGenerator.cc

@@ -20,28 +20,29 @@
 #include "RecoPixelVertexing/PixelTriplets/plugins/KDTreeLinkerAlgo.h" //amend to point at your copy...
 #include "RecoPixelVertexing/PixelTriplets/plugins/KDTreeLinkerTools.h"
 
-#include "FWCore/Utilities/interface/isFinite.h"
 #include "CommonTools/Utils/interface/DynArray.h"
 
+#include "DataFormats/Math/interface/normalizedPhi.h"
 
 #include<cstdio>
+#include<iostream>
 
 using pixelrecoutilities::LongitudinalBendingCorrection;
-typedef PixelRecoRange<float> Range;
+using Range=PixelRecoRange<float>;
 
 using namespace std;
 using namespace ctfseeding;
 
+
 PixelTripletHLTGenerator:: PixelTripletHLTGenerator(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
   : HitTripletGeneratorFromPairAndLayers(cfg),
     useFixedPreFiltering(cfg.getParameter<bool>("useFixedPreFiltering")),
     extraHitRZtolerance(cfg.getParameter<double>("extraHitRZtolerance")),
     extraHitRPhitolerance(cfg.getParameter<double>("extraHitRPhitolerance")),
     useMScat(cfg.getParameter<bool>("useMultScattering")),
-    useBend(cfg.getParameter<bool>("useBending"))
-{
-  dphi =  (useFixedPreFiltering) ?  cfg.getParameter<double>("phiPreFiltering") : 0;
-
+    useBend(cfg.getParameter<bool>("useBending")),
+    dphi(useFixedPreFiltering ?  cfg.getParameter<double>("phiPreFiltering") : 0) 
+{  
   edm::ParameterSet comparitorPSet =
     cfg.getParameter<edm::ParameterSet>("SeedComparitorPSet");
   std::string comparitorName = comparitorPSet.getParameter<std::string>("ComponentName");
@@ -87,8 +88,13 @@ void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
 
   declareDynArray(KDTreeLinkerAlgo<unsigned int>,size, hitTree);
   float rzError[size]; //save maximum errors
-  constexpr float maxphi = 2.*M_PI, minphi = -maxphi; // increase to cater for any range
-
+
+
+  const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI)/4.f : float(M_PI)/8.f; // FIXME move to config?? 
+  const float maxphi = M_PI+maxDelphi, minphi = -maxphi; // increase to cater for any range
+  const float safePhi = M_PI-maxDelphi; // sideband
+
+
   // fill the prediction vector
   for (int il=0; il<size; ++il) {
     thirdHitMap[il] = &(*theLayerCache)(thirdLayers[il], region, ev, es);
@@ -101,7 +107,7 @@ void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
                          *thirdLayers[il].detLayer(), useMScat, useBend);
 
     layerTree.clear();
-    float minv=999999.0, maxv= -999999.0; // Initialise to extreme values in case no hits
+    float minv=999999.0f, maxv= -minv; // Initialise to extreme values in case no hits
     float maxErr=0.0f;
     for (unsigned int i=0; i!=hits.size(); ++i) {
       auto angle = hits.phi(i);
@@ -111,11 +117,9 @@ void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
       float myerr = hits.dv[i];
       maxErr = std::max(maxErr,myerr);
       layerTree.emplace_back(i, angle, v); // save it
-      //FIXME the side bands do not need to be so large 
-      if (angle < 0)  // wrap all points in phi
-	{ layerTree.emplace_back(i, angle+Geom::ftwoPi(), v);}
-      else
-	{ layerTree.emplace_back(i, angle-Geom::ftwoPi(), v);}
+      // populate side-bands
+      if (angle>safePhi) layerTree.emplace_back(i, angle-Geom::ftwoPi(), v);
+      else if (angle<-safePhi) layerTree.emplace_back(i, angle+Geom::ftwoPi(), v);
     }
     KDTreeBox phiZ(minphi, maxphi, minv-0.01f, maxv+0.01f);  // declare our bounds
     //add fudge factors in case only one hit and also for floating-point inaccuracy
@@ -191,54 +195,50 @@ void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
 
 	// std::cout << "++R " << radius.min() << " " << radius.max()  << std::endl;
 
-	/*
-	Range rPhi1m = predictionRPhitmp(radius.max(), -1);
-	Range rPhi1p = predictionRPhitmp(radius.max(),  1);
-	Range rPhi2m = predictionRPhitmp(radius.min(), -1);
-	Range rPhi2p = predictionRPhitmp(radius.min(),  1);
-	Range rPhi1 = rPhi1m.sum(rPhi1p);
-	Range rPhi2 = rPhi2m.sum(rPhi2p);
-	
-	
-	auto rPhi1N = predictionRPhitmp(radius.max());
-	auto rPhi2N = predictionRPhitmp(radius.min());
-
-	std::cout << "VI " 
-		  << rPhi1N.first <<'/'<< rPhi1.first << ' '
-		  << rPhi1N.second <<'/'<< rPhi1.second << ' '
-		  << rPhi2N.first <<'/'<< rPhi2.first << ' '
-		  << rPhi2N.second <<'/'<< rPhi2.second
-		  << std::endl;
-	
-	*/
-
 	auto rPhi1 = predictionRPhitmp(radius.max());
+        bool ok1 = !rPhi1.empty();
+        if (ok1) {
+          correction.correctRPhiRange(rPhi1);
+          rPhi1.first  /= radius.max();
+          rPhi1.second /= radius.max();
+        }
 	auto rPhi2 = predictionRPhitmp(radius.min());
-
-
-	correction.correctRPhiRange(rPhi1);
-	correction.correctRPhiRange(rPhi2);
-	rPhi1.first  /= radius.max();
-	rPhi1.second /= radius.max();
-	rPhi2.first  /= radius.min();
-	rPhi2.second /= radius.min();
-	phiRange = mergePhiRanges(rPhi1,rPhi2);
+        bool ok2 = !rPhi2.empty();
+       	if (ok2) {
+	  correction.correctRPhiRange(rPhi2);
+	  rPhi2.first  /= radius.min();
+	  rPhi2.second /= radius.min();
+        }
+
+        if (ok1) { 
+          rPhi1.first = normalizedPhi(rPhi1.first);
+          rPhi1.second = proxim(rPhi1.second,rPhi1.first);
+          if(ok2) {
+            rPhi2.first = proxim(rPhi2.first,rPhi1.first);
+            rPhi2.second = proxim(rPhi2.second,rPhi1.first);
+            phiRange = rPhi1.sum(rPhi2);
+          } else phiRange=rPhi1;
+        } else if(ok2) {
+          rPhi2.first = normalizedPhi(rPhi2.first);
+          rPhi2.second = proxim(rPhi2.second,rPhi2.first);
+          phiRange=rPhi2;
+        } else continue;
+
       }
-      
+
       constexpr float nSigmaRZ = 3.46410161514f; // std::sqrt(12.f); // ...and continue as before
       constexpr float nSigmaPhi = 3.f;
 
       foundNodes.clear(); // Now recover hits in bounding box...
-      //FIXME this needs a cleanup 
       float prmin=phiRange.min(), prmax=phiRange.max();
-      if( edm::isNotFinite(prmax) | edm::isNotFinite(prmin) ) continue;
-      if ((prmax-prmin) > Geom::ftwoPi())
-	{ prmax=Geom::fpi(); prmin = -Geom::fpi();}
-      else
-	{ while (prmax>maxphi) { prmin -= Geom::ftwoPi(); prmax -= Geom::ftwoPi();}
-	  while (prmin<minphi) { prmin += Geom::ftwoPi(); prmax += Geom::ftwoPi();}
-	  // This needs range -twoPi to +twoPi to work
-	}
+
+
+      if (prmax-prmin>maxDelphi) {
+        auto prm = phiRange.mean();
+        prmin = prm - 0.5f*maxDelphi;
+        prmax = prm + 0.5f*maxDelphi;
+      }
+ 
       if (barrelLayer)
 	{
 	  Range regMax = predictionRZ.detRange();
@@ -287,6 +287,7 @@ void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
 
 	float ir = 1.f/hits.rv(KDdata);
 	float phiErr = nSigmaPhi * hits.drphi[KDdata]*ir;
+        bool nook=true;
 	for (int icharge=-1; icharge <=1; icharge+=2) {
 	  Range rangeRPhi = predictionRPhi(hits.rv(KDdata), icharge);
 	  correction.correctRPhiRange(rangeRPhi);
@@ -298,28 +299,13 @@ void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
 	    } else {
 	      LogDebug("RejectedTriplet") << "rejected triplet from comparitor ";
 	    }
-	    break;
+	    nook=false; break;
 	  } 
 	}
+        if (nook) LogDebug("RejectedTriplet") << "rejected triplet from second phicheck " << p3_phi;
       }
     }
   }
   // std::cout << "triplets " << result.size() << std::endl;
 }
 
-bool PixelTripletHLTGenerator::checkPhiInRange(float phi, float phi1, float phi2) const
-{
-  while (phi > phi2) phi -=  Geom::ftwoPi();
-  while (phi < phi1) phi +=  Geom::ftwoPi();
-  return (  (phi1 <= phi) && (phi <= phi2) );
-}  
-
-std::pair<float,float> PixelTripletHLTGenerator::mergePhiRanges(const std::pair<float,float>& r1,
-								const std::pair<float,float>& r2) const 
-{ float r2_min=r2.first;
-  float r2_max=r2.second;
-  while (r1.first-r2_min > Geom::fpi()) { r2_min += Geom::ftwoPi(); r2_max += Geom::ftwoPi();}
-  while (r1.first-r2_min < -Geom::fpi()) { r2_min -= Geom::ftwoPi();  r2_max -= Geom::ftwoPi(); }
-
-  return std::make_pair(min(r1.first,r2_min),max(r1.second,r2_max));
-}

RecoPixelVertexing/PixelTriplets/plugins/PixelTripletHLTGenerator.h

@@ -32,17 +32,12 @@ typedef CombinedHitTripletGenerator::LayerCacheType       LayerCacheType;
                             const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers) override;
 
 private:
-  bool checkPhiInRange(float phi, float phi1, float phi2) const;
-  std::pair<float,float> mergePhiRanges(
-      const std::pair<float,float> &r1, const std::pair<float,float> &r2) const; 
-
-private:
-  bool useFixedPreFiltering;
-  float extraHitRZtolerance;
-  float extraHitRPhitolerance;
-  bool useMScat;
-  bool useBend;
-  float dphi;
+  const bool useFixedPreFiltering;
+  const float extraHitRZtolerance;
+  const float extraHitRPhitolerance;
+  const bool useMScat;
+  const bool useBend;
+  const float dphi;
   std::unique_ptr<SeedComparitor> theComparitor;
 
 };