CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HICMeasurementEstimator.cc
Go to the documentation of this file.
1 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
2 #include "RecoHI/HiMuonAlgos/interface/HICMeasurementEstimator.h"
3 #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h"
4 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
5 using namespace edm;
6 using namespace std;
7 using namespace cms;
8 
9 //#define DEBUG
10 
11 std::pair<bool,double>
12 HICMeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos,
13  const TransientTrackingRecHit& aRecHit) const {
14  std::pair<bool,double> flag(false,0.);
15  if(!tsos.isValid()) {
16 #ifdef DEBUG
17  std::cout<<" HICMeasurementEstimator::estimate::trajectory is not valid "<<std::endl;
18 #endif
19  return flag;
20  }
21 
22  switch (aRecHit.dimension()) {
23  case 1: return estimate<1>(tsos,aRecHit);
24  case 2: return estimate<2>(tsos,aRecHit);
25  case 3: return estimate<3>(tsos,aRecHit);
26  case 4: return estimate<4>(tsos,aRecHit);
27  case 5: return estimate<5>(tsos,aRecHit);
28  }
29  throw cms::Exception("RecHit of invalid size (not 1,2,3,4,5)");
30 }
31 
32 template <unsigned int D> std::pair<bool,double>
33 HICMeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos,
34  const TransientTrackingRecHit& aRecHit) const {
35  typedef typename AlgebraicROOTObject<D>::Vector Vec;
36  typedef typename AlgebraicROOTObject<D>::SymMatrix Mat;
37  double est = 0.;
38 // If RecHit is not valid
39  if(!(aRecHit.isValid())) {
40 #ifdef DEBUG
41  std::cout<<" Measurement estimator::RecHit is invalid "<<std::endl;
42 #endif
43  return HitReturnType(false,est);
44  }
45 // Check windows
46 
47  double dphi = fabs(tsos.freeTrajectoryState()->parameters().position().phi() - aRecHit.globalPosition().phi() - thePhiBoundMean);
48  double dz = fabs( tsos.freeTrajectoryState()->parameters().position().z() - aRecHit.globalPosition().z() - theZBoundMean );
49  double dr = fabs( tsos.freeTrajectoryState()->parameters().position().perp() - aRecHit.globalPosition().perp() - theZBoundMean );
50 #ifdef DEBUG
51  std::cout<<" Momentum "<<tsos.freeTrajectoryState()->parameters().momentum().perp()<<" "<<tsos.freeTrajectoryState()->parameters().momentum().z()<<std::endl;
52  std::cout<<" RecHit position r "<<aRecHit.globalPosition().perp()<<" phi "<<aRecHit.globalPosition().phi()<<" "<<aRecHit.globalPosition().z()<<std::endl;
53  std::cout<<" Predicted position "<<tsos.freeTrajectoryState()->parameters().position().perp()<<" "<<tsos.freeTrajectoryState()->parameters().position().phi()<<
54  " "<<tsos.freeTrajectoryState()->parameters().position().z()<<std::endl;
55  std::cout<<" HICMeasurementEstimator::phi "<<dphi<<" "<<thePhiBound<<std::endl;
56  std::cout<<" HICMeasurementEstimator::z "<<dz<<" "<<theZBound<<std::endl;
57  std::cout<<" HICMeasurementEstimator::z "<<dr<<" "<<theZBound<<std::endl;
58 #endif
59  if( dphi > thePhiBound ) {
60 #ifdef DEBUG
61  std::cout<<" HICMeasurementEstimator::phi::failed "<<std::endl;
62 #endif
63  return HitReturnType(false,est);
64  }
65  if( dz > theZBound ) {
66 #ifdef DEBUG
67  std::cout<<" HICMeasurementEstimator::z::failed "<<std::endl;
68 #endif
69  return HitReturnType(false,est);
70  }
71  if( dr > theZBound ) {
72 #ifdef DEBUG
73  std::cout<<" HICMeasurementEstimator::r::failed "<<std::endl;
74 #endif
75  return HitReturnType(false,est);
76  }
77 
78 
79  MeasurementExtractor me(tsos);
80  Vec r = asSVector<D>(aRecHit.parameters()) - me.measuredParameters<D>(aRecHit);
81  Mat R = asSMatrix<D>(aRecHit.parametersError()) + me.measuredError<D>(aRecHit);
82  //int ierr = ! R.Invert(); // if (ierr != 0) throw exception; //
83  R.Invert();
84  est = ROOT::Math::Similarity(r, R);
85 
86  if( est > theChi2Cut )
87  {
88 #ifdef DEBUG
89  std::cout<<" HICMeasurementEstimator::chi2::failed "<<est<<" "<<theChi2Cut<<std::endl;
90 #endif
91 
92  return HitReturnType(false,est);
93  }
94 
95  return HitReturnType(true,est);
96 }
97 
98 bool HICMeasurementEstimator::estimate( const TrajectoryStateOnSurface& ts,
99  const BoundPlane& plane) const
100 {
101 
102 // cout<<" start estimate plane "<<endl;
103  double pi = 4.*atan(1.);
104  double twopi = 2.*pi;
105  float theZError = plane.bounds().length() + 4.;
106  float thePhiError = 2.*plane.bounds().width()/plane.position().perp();
107 // Change 02.07.08
108 // float thePhiError = 4.*plane.bounds().width()/plane.position().perp();
109 
110 #ifdef DEBUG
111  cout<<" ======================================================================================== ";
112  cout<<" Estimate detector:: tsos : detector : Error "<<endl;
113  cout<<" R "<<ts.globalPosition().perp()<<" "<<plane.position().perp()<<" "<<theZError<<endl;
114  cout<<" Phi "<<ts.globalPosition().phi()<<" "<<plane.position().phi()<<" "<<thePhiError<<endl;
115  cout<<" Z "<<ts.globalPosition().z()<<" "<<plane.position().z()<<" "<<theZError<<endl;
116 #endif
117 
118  bool flag = false;
119  if(fabs(ts.globalPosition().perp()-plane.position().perp())<theZError){
120  if(fabs(ts.globalPosition().z()-plane.position().z())<theZError){
121  float phi1 = ts.globalPosition().phi();
122  float phi2 = plane.position().phi();
123  if(phi1<0.) phi1 = twopi+phi1;
124  if(phi2<0.) phi2 = twopi+phi2;
125  float dfi = fabs(phi1-phi2);
126  if(dfi>pi) dfi = twopi-dfi;
127  if(dfi<thePhiError) flag = true;
128  }
129  }
130 #ifdef DEBUG
131  cout<<" Estimate = "<<flag<<endl;
132 #endif
133 
134  return flag;
135 
136 }
137 
138 vector<double> HICMeasurementEstimator::setCuts(Trajectory& traj, const DetLayer* b)
139 {
140  vector<double> theCuts;
141  const DetLayer* a = traj.data().back().layer();
142  const DetLayer* first = traj.data().front().layer();
143 // const DetLayer* last = traj.data().front().layer();
144  thePhiWinMean = 0.;
145  theZWinMean = 0.;
146  thePhiWin = 0.;
147  theZWin = 0.;
148  theNewCut = 11.; // change 5->10 03.07.2008 //10-11 23.06.09
149  theNewCutB = 5.;
150 
151  thePhiWinMeanB = 0.002;
152  theZWinMeanB = 0.;
153  thePhiWinB = 0.008;
154  theZWinB = 17.;
155 
156  theZCutMean = 0.;
157  thePhiCutMean = 0.;
158  thePhiCut = 0.;
159  theZCut = 0.;
160 
161  theLayer = b;
162  theLastLayer = a;
163 
164 
165  theTrajectorySize = traj.data().size();
166 
167 
168  if( theBarrel.size() == 0 || theForward.size() == 0 )
169  {
170 #ifdef DEBUG
171  cout<<" HICMeasurementEstimator::setCuts:: no datector map "<<endl;
172 #endif
173  return theCuts;
174  }
175 
176  if( a->location() == GeomDetEnumerators::barrel )
177  {
178  if( first->location() == GeomDetEnumerators::barrel )
179  {
180  thePhiWin = (*theHICConst).phiwinbar[(*theBarrel.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
181  theZWin = (*theHICConst).zwinbar[(*theBarrel.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
182  thePhiCut = (*theHICConst).phicutbar[(*theBarrel.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
183  theZCut = (*theHICConst).zcutbar[(*theBarrel.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
184 // cout<<" Barrel first -Barrel cuts::layers "<<(*theBarrel.find(first)).second<<" "<<(*theBarrel.find(a)).second<<" "<<(*theBarrel.find(b)).second<<endl;
185 // cout<<" Barrel first -Barrel cuts "<< thePhiWin<<" "<<theZWin <<" "<<thePhiCut <<" "<<theZCut<<endl;
186  }
187  else
188  {
189  if(first->surface().position().z() > 0. )
190  {
191  thePhiWin = (*theHICConst).phiwinfbb[(*theForward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
192  theZWin = (*theHICConst).zwinfbb[(*theForward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
193  thePhiCut = (*theHICConst).phicutfbb[(*theForward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
194  theZCut = (*theHICConst).zcutfbb[(*theForward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
195 // cout<<" Endcap first positive -Barrel cuts::layers "<<(*theForward.find(first)).second<<" "<<(*theBarrel.find(a)).second<<" "<<(*theBarrel.find(b)).second<<endl;
196 // cout<<" Endcap first positive -Barrel cuts "<< thePhiWin<<" "<<theZWin <<" "<<thePhiCut <<" "<<theZCut<<endl;
197  } else {
198  thePhiWin = (*theHICConst).phiwinfbb[(*theBackward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
199  theZWin = (*theHICConst).zwinfbb[(*theBackward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
200  thePhiCut = (*theHICConst).phicutfbb[(*theBackward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
201  theZCut = (*theHICConst).zcutfbb[(*theBackward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
202 // cout<<" Endcap first negative -Barrel cuts::layers "<<(*theBackward.find(first)).second<<" "<<(*theBarrel.find(a)).second<<" "<<(*theBarrel.find(b)).second<<endl;
203 // cout<<" Endcap first negative -Barrel cuts "<< thePhiWin<<" "<<theZWin <<" "<<thePhiCut <<" "<<theZCut<<endl;
204  }
205  } //
206 
207 
208  theCuts.push_back(thePhiWin); theCuts.push_back(theZWin);
209  theCuts.push_back(thePhiCut); theCuts.push_back(theZCut);
210 
211  return theCuts;
212  }
213  if( a->location() == GeomDetEnumerators::endcap && b->location() == GeomDetEnumerators::endcap)
214  {
215  if( a->surface().position().z() > 0. )
216  {
217  thePhiWin = (*theHICConst).phiwinfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theForward.find(b)).second];
218  theZWin = (*theHICConst).zwinfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theForward.find(b)).second];
219  thePhiCut = (*theHICConst).phicutfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theForward.find(b)).second];
220  theZCut = (*theHICConst).zcutfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theForward.find(b)).second];
221  }
222  else
223  {
224  thePhiWin = (*theHICConst).phiwinfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBackward.find(b)).second];
225  theZWin = (*theHICConst).zwinfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBackward.find(b)).second];
226  thePhiCut = (*theHICConst).phicutfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBackward.find(b)).second];
227  theZCut = (*theHICConst).zcutfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBackward.find(b)).second];
228  }
229 
230  if( theLowMult == 1 )
231  {
232  if( b->subDetector() == GeomDetEnumerators::PixelEndcap ) theNewCut = 20.;
233  if( traj.measurements().size() == 1 ) theNewCut = 20.;
234  theNewCutB = 30.;
235  }
236 
237  thePhiWinMeanB = 0.004;
238  thePhiWinB = 0.05;
239 
240 
241  theCuts.push_back(thePhiWin); theCuts.push_back(theZWin);
242  theCuts.push_back(thePhiCut); theCuts.push_back(theZCut);
243 
244  return theCuts;
245  }
246  if( a->location() == GeomDetEnumerators::endcap && b->location() == GeomDetEnumerators::barrel )
247  {
248 
249  if( a->surface().position().z() > 0. )
250  {
251  thePhiWin = (*theHICConst).phiwinbfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theBarrel.find(b)).second];
252  theZWin = (*theHICConst).zwinbfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theBarrel.find(b)).second];
253  thePhiCut = (*theHICConst).phicutbfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
254  theZCut = (*theHICConst).zcutbfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
255  }
256  else
257  {
258  thePhiWin = (*theHICConst).phiwinbfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBarrel.find(b)).second];
259  theZWin = (*theHICConst).zwinbfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBarrel.find(b)).second];
260  thePhiCut = (*theHICConst).phicutbfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
261  theZCut = (*theHICConst).zcutbfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
262  }
263 
264  if( b->subDetector() == GeomDetEnumerators::PixelBarrel) theNewCut = 20.;
265 
266  thePhiWinMeanB = 0.004;
267  thePhiWinB = 0.016;
268 
269  theCuts.push_back(thePhiWin); theCuts.push_back(theZWin);
270  theCuts.push_back(thePhiCut); theCuts.push_back(theZCut);
271 
272  return theCuts;
273  }
274 // cout<<" HICMeasurementEstimator::setCuts::Error: unknown detector layer "<<endl;
275  return theCuts;
276 }
277 
278 void HICMeasurementEstimator::chooseCuts(int& i)
279 {
280 
281  theChi2Cut = theNewCut;
282 // cout<<" Choose Chi2Cut "<<theChi2Cut<<endl;
283 
284  if( i == 1 )
285  {
286  thePhiBound = thePhiWin;
287  theZBound = theZWin;
288  thePhiBoundMean = thePhiWinMean;
289  theZBoundMean = theZWinMean;
290 // cout<<" HICMeasurementEstimator::chooseCuts "<<i<<" "<<thePhiBound<<" "<<theZBound<<endl;
291  }
292  if( i == 2 )
293  {
294  thePhiBound = thePhiWinB;
295  theZBound = theZWinB;
296  thePhiBoundMean = thePhiWinMean;
297  theZBoundMean = theZWinMean;
298 // cout<<" HICMeasurementEstimator::chooseCuts "<<i<<" "<<thePhiBound<<" "<<theZBound<<endl;
299  }
300  if( i == 3 )
301  {
302  thePhiBound = thePhiCut;
303  theZBound = theZCut;
304  thePhiBoundMean = thePhiCutMean;
305  theZBoundMean = theZCutMean;
306 // cout<<" HICMeasurementEstimator::chooseCuts "<<i<<" "<<thePhiBound<<" "<<theZBound<<endl;
307  }
308 
309  theCutType = i;
310 }
311 
312 int HICMeasurementEstimator::getDetectorCode(const DetLayer* a)
313 {
314  int layer = 0;
315  if( a->location() == GeomDetEnumerators::barrel )
316  {
317  layer = (*theBarrel.find(a)).second;
318  }
319  if( a->location() == GeomDetEnumerators::endcap )
320  {
321  if( a->surface().position().z() > 0. ) { layer = 100+(*theForward.find(a)).second;}
322  if( a->surface().position().z() < 0. ) { layer = -100-(*theBackward.find(a)).second;}
323  }
324  return layer;
325 }
326 
327 void HICMeasurementEstimator::setHICDetMap()
328 {
329 
330 #ifdef DEBUG
331  std::cout<<" Set Detector Map... "<<std::endl;
332 #endif
333 
334  int ila=0;
335  for ( std::vector<BarrelDetLayer*>::const_iterator ilayer = bl.begin(); ilayer != bl.end(); ilayer++)
336  {
337  theBarrel[(*ilayer)]=ila;
338  ila++;
339  }
340 //
341 // The same for forward part.
342 //
343 
344  int ilf1 = 0;
345  int ilf2 = 0;
346  for ( vector<ForwardDetLayer*>::const_iterator ilayer = fpos.begin();
347  ilayer != fpos.end(); ilayer++)
348  {
349  theForward[(*ilayer)] = ilf1;
350  ilf1++;
351  }
352  for ( vector<ForwardDetLayer*>::const_iterator ilayer = fneg.begin();
353  ilayer != fneg.end(); ilayer++)
354  {
355 #ifdef DEBUG
356  cout<<" HICDetectorMap::negative layers "<<(**ilayer).position().z()<<" "<<ilf2<<endl;
357 #endif
358  theBackward[(*ilayer)] = ilf2;
359  ilf2++;
360  }
361 
362 }
GeometricSearchDet::surface
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
i
int i
Definition: DBlmapReader.cc:9
TrackingRecHit::dimension
virtual int dimension() const =0
flag
long int flag
Definition: mlp_lapack.h:47
Bounds::length
virtual float length() const =0
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:66
AlgebraicROOTObject::SymMatrix
ROOT::Math::SMatrix< double, D1, D1, ROOT::Math::MatRepSym< double, D1 > > SymMatrix
Definition: AlgebraicROOTObjects.h:66
DetLayer::location
virtual Location location() const =0
Which part of the detector (barrel, endcap)
FreeTrajectoryState::parameters
const GlobalTrajectoryParameters & parameters() const
Definition: FreeTrajectoryState.h:101
TrajectoryStateOnSurface::freeTrajectoryState
FreeTrajectoryState * freeTrajectoryState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:119
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
DetLayer::subDetector
virtual SubDetector subDetector() const =0
The type of detector (PixelBarrel, PixelEndcap, TIB, TOB, TID, TEC, CSC, DT, RPCBarrel, RPCEndcap)
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:128
HICMeasurementEstimator::setCuts
virtual std::vector< double > setCuts(Trajectory &traj, const DetLayer *theCurrentLayer)
Definition: HICMeasurementEstimator.cc:138
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:203
TrackingRecHit::parameters
virtual AlgebraicVector parameters() const =0
Trajectory::measurements
DataContainer const & measurements() const
Definition: Trajectory.h:169
PV3DBase::z
T z() const
Definition: PV3DBase.h:58
BoundSurface::bounds
const Bounds & bounds() const
Definition: BoundSurface.h:89
AlgebraicROOTObject::Vector
ROOT::Math::SVector< double, D1 > Vector
Definition: AlgebraicROOTObjects.h:65
TrackingRecHit::isValid
bool isValid() const
Definition: TrackingRecHit.h:63
b
double b
Definition: hdecay.h:120
HICMeasurementEstimator::getDetectorCode
virtual int getDetectorCode(const DetLayer *a)
Definition: HICMeasurementEstimator.cc:312
Trajectory::data
DataContainer const & data() const
obsolete name, use measurements() instead.
Definition: Trajectory.h:171
TransientTrackingRecHit::globalPosition
virtual GlobalPoint globalPosition() const
Definition: TransientTrackingRecHit.cc:12
a
double a
Definition: hdecay.h:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
HICMeasurementEstimator::estimate
virtual std::pair< bool, double > estimate(const TrajectoryStateOnSurface &, const TransientTrackingRecHit &) const
Definition: HICMeasurementEstimator.cc:12
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
edm::first
T first(std::pair< T, U > const &p)
Definition: ParameterSet.cc:199
Bounds::width
virtual float width() const =0
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:41
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:150
TrackingRecHit::parametersError
virtual AlgebraicSymMatrix parametersError() const =0
Generated for CMSSW Reference Manual by doxygen 1.8.5