CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
LowPtClusterShapeSeedComparitor.cc
Go to the documentation of this file.
2 
6 
11 
14 
16 
17 inline float sqr(float x) { return x*x; }
18 
19 using namespace std;
20 
21 /*****************************************************************************/
23  (const edm::ParameterSet& ps)
24 {
25 
26 }
27 
28 /*****************************************************************************/
30 {
31 }
32 
33 /*****************************************************************************/
35  (const Global2DVector& a, const Global2DVector& b)
36 {
37  return a.x() * b.y() - a.y() * b.x();
38 }
39 
40 /*****************************************************************************/
42  (const vector<GlobalPoint> & g)
43 {
44  // Get 2d points
45  vector<Global2DVector> p;
46  for(vector<GlobalPoint>::const_iterator ig = g.begin();
47  ig!= g.end(); ig++)
48  p.push_back( Global2DVector(ig->x(), ig->y()) );
49 
50  //
51  vector<GlobalVector> globalDirs;
52 
53  // Determine circle
54  CircleFromThreePoints circle(g[0],g[1],g[2]);
55 
56  if(circle.curvature() != 0.)
57  {
58  Global2DVector c (circle.center().x(), circle.center().y());
59 
60  float rad2 = (p[0] - c).mag2();
61  float area = fabsf(areaParallelogram(p[1] - p[0], p[1] - c));
62 
63  float a12;
64  if(area >= rad2) a12 = M_PI/2;
65  else a12 = asin(area / rad2);
66 
67  float slope = (g[1].z() - g[0].z()) / a12;
68 
69  float cotTheta = slope * circle.curvature(); // == sinhEta
70  float coshEta = sqrt(1 + sqr(cotTheta)); // == 1/sinTheta
71 
72  // Calculate globalDirs
73  float sinTheta = 1. / coshEta;
74  float cosTheta = cotTheta * sinTheta;
75 
76  int dir;
77  if(areaParallelogram(p[0] - c, p[1] - c) > 0) dir = 1; else dir = -1;
78 
79  float curvature = circle.curvature();
80 
81  for(vector<Global2DVector>::const_iterator ip = p.begin();
82  ip!= p.end(); ip++)
83  {
84  Global2DVector v = (*ip - c)*curvature*dir;
85  globalDirs.push_back(GlobalVector(-v.y()*sinTheta,
86  v.x()*sinTheta,
87  cosTheta));
88  }
89  }
90  else{
91  LogDebug("LowPtClusterShapeSeedComparitor")<<"the curvature is null:"
92  <<"\n point1: "<<g[0]
93  <<"\n point2: "<<g[1]
94  <<"\n point3: "<<g[2];
95  }
96  return globalDirs;
97 }
98 
99 /*****************************************************************************/
102 {
103  vector<GlobalPoint> globalPoss;
104 
105  for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator recHit = thits.begin();
106  recHit!= thits.end();
107  recHit++)
108  {
109  DetId detId = (*recHit)->hit()->geographicalId();
110 
111  GlobalPoint gpos = (*recHit)->globalPosition();
112 
113  globalPoss.push_back(gpos);
114  }
115 
116  return globalPoss;
117 }
118 
119 /*****************************************************************************/
121  const edm::EventSetup &es)
122 //(const reco::Track* track, const vector<const TrackingRecHit *> & recHits) const
123 {
124 
125  // Get cluster shape hit filter
127  es.get<CkfComponentsRecord>().get("ClusterShapeHitFilter",shape);
128  theFilter = shape.product();
129 
131 
132  // Get global positions
133  vector<GlobalPoint> globalPoss = getGlobalPoss(thits);
134 
135  // Get global directions
136  vector<GlobalVector> globalDirs = getGlobalDirs(globalPoss);
137 
138  bool ok = true;
139 
140  // Check whether shape of pixel cluster is compatible
141  // with local track direction
142 
143  if (globalDirs.size()!=globalPoss.size() || globalDirs.size()!=thits.size())
144  {
145  LogDebug("LowPtClusterShapeSeedComparitor")<<"not enough global dir calculated:"
146  <<"\nnHits: "<<thits.size()
147  <<"\nnPos: "<<globalPoss.size()
148  <<"\nnDir: "<<globalDirs.size()
149  <<" will say the seed is good anyway.";
150  return true;
151  }
152 
153  for(int i = 0; i < 3; i++)
154  {
155  const SiPixelRecHit* pixelRecHit =
156  dynamic_cast<const SiPixelRecHit *>(thits[i]->hit());
157 
158  if (!pixelRecHit){
159  edm::LogError("LowPtClusterShapeSeedComparitor")<<"this is not a pixel cluster";
160  ok = false; break;
161  }
162 
163  if(!pixelRecHit->isValid())
164  {
165  ok = false; break;
166  }
167 
168  LogDebug("LowPtClusterShapeSeedComparitor")<<"about to compute compatibility."
169  <<"hit ptr: "<<pixelRecHit
170  <<"global direction:"<< globalDirs[i];
171 
172 
173  if(! theFilter->isCompatible(*pixelRecHit, globalDirs[i]) )
174  {
175  LogTrace("LowPtClusterShapeSeedComparitor")
176  << " clusShape is not compatible"
177  << HitInfo::getInfo(*thits[i]->hit());
178 
179  ok = false; break;
180  }
181  }
182 
183  return ok;
184 }
185 
#define LogDebug(id)
static std::string getInfo(const DetId &id)
Definition: HitInfo.cc:24
int i
Definition: DBlmapReader.cc:9
T y() const
Definition: PV2DBase.h:40
const RecHits & container() const
Definition: SeedingHitSet.h:18
static const double slope[3]
float areaParallelogram(const Global2DVector &a, const Global2DVector &b)
virtual bool compatible(const SeedingHitSet &hits, const edm::EventSetup &es)
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
T curvature(T InversePt, const edm::EventSetup &iSetup)
LowPtClusterShapeSeedComparitor(const edm::ParameterSet &ps)
T sqrt(T t)
Definition: SSEVec.h:28
T y() const
Cartesian y coordinate.
Definition: Basic2DVector.h:65
std::vector< GlobalVector > getGlobalDirs(const std::vector< GlobalPoint > &globalPoss)
#define LogTrace(id)
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: DetId.h:20
#define M_PI
Definition: BFit3D.cc:3
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
bool isValid() const
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
Square< F >::type sqr(const F &f)
Definition: Square.h:13
dbl *** dir
Definition: mlp_gen.cc:35
T x() const
Definition: PV2DBase.h:39
mathSSE::Vec4< T > v
std::vector< GlobalPoint > getGlobalPoss(const TransientTrackingRecHit::ConstRecHitContainer &recHits)
T x() const
Cartesian x coordinate.
Definition: Basic2DVector.h:62
Global3DVector GlobalVector
Definition: GlobalVector.h:10
Our base class.
Definition: SiPixelRecHit.h:27