CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
StripCPEgeometric.cc
Go to the documentation of this file.
4 #include <numeric>
5 
7 localParameters( const SiStripCluster& cluster, const GeomDetUnit& det, const LocalTrajectoryParameters& ltp) const {
8  StripCPE::Param const& p = param(det);
9 
10  const LocalPoint& pos = ltp.position();
11  LocalVector track = ltp.momentum();
12  track *= (track.z()<0) ? fabs(p.thickness/track.z()) :
13  (track.z()>0) ? -fabs(p.thickness/track.z()) :
14  p.maxLength/track.mag() ;
15 
16  const float fullProjection = p.coveredStrips( track+p.drift, pos );
17  stats_t<float> projection;
18  {
19  const float absProj = fabs( fullProjection );
20  const float minProj = 2*p.thickness*tan_diffusion_angle/p.topology->localPitch(pos);
21  const float projection_rel_err2 = thickness_rel_err2 + p.pitch_rel_err2;
22  projection = stats_t<float>::from_relative_uncertainty2( std::max(absProj, minProj), projection_rel_err2);
23  }
24 
25  const std::vector<stats_t<float> > Q = reco::InverseCrosstalkMatrix::unfold( cluster.amplitudes(), xtalk1[p.moduleGeom] );
26  const stats_t<float> strip = cluster.firstStrip() + offset_from_firstStrip( Q, projection );
27 
28  const float corrected = strip() - 0.5*(1-shift[p.moduleGeom]) * fullProjection
29  + 0.5*p.coveredStrips(track, ltp.position());
30  const float error2 = std::max( strip.error2(), minimum_uncertainty_squared );
31 
32  return std::make_pair( p.topology->localPosition( corrected, ltp.vector()),
33  p.topology->localError( corrected, error2, ltp.vector()));
34 }
35 
36 
38 offset_from_firstStrip( const std::vector<stats_t<float> >& Q, const stats_t<float>& proj) const {
39  WrappedCluster wc(Q);
40  if( useNPlusOne( wc, proj)) wc.addSuppressedEdgeStrip(); else
41  while( useNMinusOne( wc, proj)) wc.dropSmallerEdgeStrip();
42 
43  if( proj() < wc.N-2 ) return stats_t<float>( wc.middle(), pow(wc.N-proj(),2) / 12.);
44  if( wc.deformed() ) return stats_t<float>( wc.centroid()(), 1 / 12.);
45  if( proj > wc.maxProjection() ) return stats_t<float>( wc.centroid()(), 1 / 12.);
46 
47  if( ambiguousSize( wc, proj) ) {
48  const stats_t<float> probably = geometric_position( wc, proj);
50  const stats_t<float> maybe = geometric_position( wc, proj);
51  return stats_t<float>( probably(), std::max( probably.error2(), float(maybe.error2() + pow( probably()-maybe() ,2)/12 )) );
52  }
53  return geometric_position( wc, proj);
54 }
55 
58  const stats_t<float> x = wc.middle() + 0.5 * proj * wc.eta();
59  return wc.N==1
60  ? stats_t<float>( x(), pow( 1-0.82*proj(), 2 ) / 12 )
61  : stats_t<float>( x(), scaling_squared * x.error2() ) ;
62 }
63 
64 inline
66 useNPlusOne(const WrappedCluster& wc, const stats_t<float>& proj) const
67 { return wc.maxProjection() < proj && proj() < wc.N+1 ; }
68 
69 inline
71 useNMinusOne(const WrappedCluster& wc, const stats_t<float>& proj) const {
72  if( proj() > wc.N-1) return false;
73  if( wc.smallerEdgeStrip() < 0 ) return true;
74  if( proj() < wc.N-3) return false;
75  if( proj() < wc.N-2) return true;
76  if( wc.eta().sigmaFrom(0) < 3) return false;
77 
78  WrappedCluster wcTest(wc); wcTest.dropSmallerEdgeStrip();
79  if( proj >= wcTest.maxProjection() ) return false;
80  if( wc.sign()*wc.eta()() > 1./(wc.N-1) ) return true;
81 
82  return wc.smallerEdgeStrip().sigmaFrom(0) < noise_threshold;
83 }
84 
85 inline
87 ambiguousSize( const WrappedCluster& wc, const stats_t<float>& proj) const
88 { return
89  proj() < wc.N-1 &&
90  wc.smallerEdgeStrip()()>0 &&
91  wc.smallerEdgeStrip().sigmaFrom(0) < maybe_noise_threshold; }
92 
94 WrappedCluster(const std::vector<stats_t<float> >& Q) :
95  N(Q.size()-2),
96  clusterFirst(Q.begin()+1),
97  first(clusterFirst)
98  {}
99 
100 inline
103  if( eta().sigmaFrom(0) < 1 ) { first--; N+=2; }
104  else if( *first > last() ) { first--; N+=1; }
105  else { N+=1; }
106 }
107 
108 inline
111  if( *first < last() ) { first++; N-=1; } else
112  if( last() < *first ) { N-=1; } else
113  { first++; N-=2; }
114 }
115 
116 inline
118 middle() const
119 { return (first-clusterFirst) + N/2.;}
120 
121 inline
123 centroid() const {
124  stats_t<float> sumXQ(0);
125  for(std::vector<stats_t<float> >::const_iterator i = first; i<first+N; i++) sumXQ += (i-clusterFirst)*(*i);
126  return sumXQ/sumQ() + 0.5;
127 }
128 
129 inline
131 sumQ() const
132 { return accumulate(first, first+N, stats_t<float>(0));}
133 
134 inline
136 eta() const
137 { return ( last() - *first ) / sumQ() ;}
138 
139 inline
141 deformed() const
142 { return N>2 && std::max((*first)(),last()()) > accumulate(first+1,first+N-1,stats_t<float>(0))() / (N-2);}
143 
144 inline
147 { return N * (1 + sign()*eta() ).inverse(); }
148 
149 inline
152 { return std::min(*first, last()); }
153 
154 inline
156 sign() const
157 { return ( *first < last() ) ? 1 : -1; }
int i
Definition: DBlmapReader.cc:9
void strip(std::string &input, const std::string &blanks=" \n\t")
Definition: stringTools.cc:16
WrappedCluster(const std::vector< stats_t< float > > &)
float pitch_rel_err2
Definition: StripCPE.h:43
LocalPoint position() const
Local x and y position coordinates.
float thickness
Definition: StripCPE.h:43
StripTopology const * topology
Definition: StripCPE.h:41
stats_t< float > smallerEdgeStrip() const
const float scaling_squared
#define min(a, b)
Definition: mlp_lapack.h:161
T error2() const
uint16_t firstStrip() const
T eta() const
std::pair< LocalPoint, LocalError > LocalValues
stats_t< float > centroid() const
StripClusterParameterEstimator::LocalValues localParameters(const SiStripCluster &, const GeomDetUnit &, const LocalTrajectoryParameters &) const
std::vector< double > xtalk1
Definition: StripCPE.h:36
AlgebraicVector5 vector() const
virtual float localPitch(const LocalPoint &) const =0
stats_t< float > geometric_position(const WrappedCluster &, const stats_t< float > &) const
T mag() const
Definition: PV3DBase.h:66
SiStripDetId::ModuleGeometry moduleGeom
Definition: StripCPE.h:45
const T & max(const T &a, const T &b)
static stats_t from_relative_uncertainty2(T q, T re2)
T z() const
Definition: PV3DBase.h:63
stats_t< float > sumQ() const
stats_t< float > offset_from_firstStrip(const std::vector< stats_t< float > > &, const stats_t< float > &) const
std::vector< double > shift
Definition: StripCPE.h:35
LocalVector momentum() const
Momentum vector in the local frame.
float coveredStrips(const LocalVector &, const LocalPoint &) const
Definition: StripCPE.cc:78
bool first
Definition: L1TdeRCT.cc:94
const float maybe_noise_threshold
#define N
Definition: blowfish.cc:9
bool useNPlusOne(const WrappedCluster &, const stats_t< float > &) const
const float minimum_uncertainty_squared
bool ambiguousSize(const WrappedCluster &, const stats_t< float > &) const
const float tan_diffusion_angle
virtual LocalError localError(float strip, float stripErr2) const =0
#define begin
Definition: vmac.h:31
float maxLength
Definition: StripCPE.h:43
Param const & param(const GeomDetUnit &det) const
Definition: StripCPE.h:48
LocalVector drift
Definition: StripCPE.h:42
virtual LocalPoint localPosition(float strip) const =0
const float thickness_rel_err2
Definition: DDAxes.h:10
const float noise_threshold
stats_t< float > maxProjection() const
tuple size
Write out results.
const std::vector< uint8_t > & amplitudes() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
static std::vector< stats_t< float > > unfold(const std::vector< uint8_t > &q, const float x)
bool useNMinusOne(const WrappedCluster &, const stats_t< float > &) const