CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiTrackerMultiRecHitUpdator.cc
Go to the documentation of this file.
4 //#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
10 
11 
13  const TrackingRecHitPropagator* hitpropagator,
14  const float Chi2Cut,
15  const std::vector<double>& anAnnealingProgram):
16  theBuilder(builder),
17  theHitPropagator(hitpropagator),
18  theChi2Cut(Chi2Cut),
19  theAnnealingProgram(anAnnealingProgram){}
20 //theAnnealingStep(0),
21 //theIsUpdating(true){}
22 
23 
26  float annealing) const{
28  for (std::vector<const TrackingRecHit*>::const_iterator iter = rhv.begin(); iter != rhv.end(); iter++){
30  if (transient->isValid()) tcomponents.push_back(transient);
31  }
32 
33  return update(tcomponents, tsos, annealing);
34 
35 }
36 
39  double annealing) const{
40  LogTrace("SiTrackerMultiRecHitUpdator") << "Calling SiTrackerMultiRecHitUpdator::update with AnnealingFactor: " << annealing;
41  if (original->isValid())
42  LogTrace("SiTrackerMultiRecHitUpdator") << "Original Hit position " << original->localPosition() << " original error "
43  << original->parametersError();
44  else LogTrace("SiTrackerMultiRecHitUpdator") << "Invalid hit";
45 
46  if(!tsos.isValid()) {
47  //return original->clone();
48  throw cms::Exception("SiTrackerMultiRecHitUpdator") << "!!! MultiRecHitUpdator::update(..): tsos NOT valid!!! ";
49  }
50 
51  //check if to clone is the right thing
52  if (original->transientHits().empty()) return original->clone(tsos);
53 
54  TransientTrackingRecHit::ConstRecHitContainer tcomponents = original->transientHits();
55  return update(tcomponents, tsos, annealing);
56 }
57 
60  double annealing) const{
61 
62  if (tcomponents.empty()){
63  LogTrace("SiTrackerMultiRecHitUpdator") << "Empty components vector passed to SiTrackerMultiRecHitUpdator::update, returning an InvalidTransientRecHit ";
65  }
66 
67  if(!tsos.isValid()) {
68  LogTrace("SiTrackerMultiRecHitUpdator")<<"SiTrackerMultiRecHitUpdator::update: tsos NOT valid!!!, returning an InvalidTransientRecHit";
70  }
71 
72  std::vector<TransientTrackingRecHit::RecHitPointer> updatedcomponents;
73  const GeomDet* geomdet = 0;
74  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iter = tcomponents.begin(); iter != tcomponents.end(); iter++){
75  if (iter == tcomponents.begin()) {
76  if (&((*iter)->det()->surface())!=&(tsos.surface())){
77  throw cms::Exception("SiTrackerMultiRecHitUpdator") << "the Trajectory state and the first rechit passed to the SiTrackerMultiRecHitUpdator lay on different surfaces!: state lays on surface " << tsos.surface().position() << " hit with detid " << (*iter)->det()->geographicalId().rawId() << " lays on surface " << (*iter)->det()->surface().position();
78  }
79  geomdet = (*iter)->det();
80  LogTrace("SiTrackerMultiRecHitUpdator") << "Current reference surface located at " << geomdet->surface().position();
81  // LogTrace("SiTrackerMultiRecHitUpdator")<< "TSOS position " << tsos.localPosition();
82  }
83  if (&((*iter)->det()->surface())!=&(tsos.surface())){
85  // LogTrace("SiTrackerMultiRecHitUpdator") << "hit propagated";
86 
87  if (cloned->isValid()) updatedcomponents.push_back(cloned);
88  } else {
89  TransientTrackingRecHit::RecHitPointer cloned = (*iter)->clone(tsos);
90  if (cloned->isValid()) updatedcomponents.push_back(cloned);
91  }
92  }
93  // LogTrace("SiTrackerMultiRecHitUpdator") << "hit cloned";
94  int ierr;
95  std::vector<std::pair<const TrackingRecHit*, float> > mymap;
96  std::vector<std::pair<const TrackingRecHit*, float> > normmap;
97 
98  double a_sum=0, c_sum=0;
99 
100  AlgebraicVector2 tsospos;
101  tsospos[0]=tsos.localPosition().x();
102  tsospos[1]=tsos.localPosition().y();
103  LogTrace("SiTrackerMultiRecHitUpdator")<< "TSOS position " << tsos.localPosition();
104  for(std::vector<TransientTrackingRecHit::RecHitPointer>::iterator ihit = updatedcomponents.begin(); ihit != updatedcomponents.end(); ihit++) {
105  AlgebraicVector2 r(asSVector<2>((*ihit)->parameters()) - tsospos);
106  AlgebraicSymMatrix22 V(asSMatrix<2>((*ihit)->parametersError()));
107  V *= annealing;//assume that TSOS is smoothed one
108  //V += me.measuredError(*ihit);// result = b*V + H*C*H.T()
109  AlgebraicSymMatrix22 W(V.Inverse(ierr));
110  double det;
111  bool ierr2=!(V.Det2(det));
112 
113  if(ierr != 0|| ierr2) {
114  LogTrace("SiTrackerMultiRecHitUpdator")<<"MultiRecHitUpdator::update: W not valid!"<<std::endl;
115  LogTrace("SiTrackerMultiRecHitUpdator")<<"V: "<<V<<" AnnealingFactor: "<<annealing<<std::endl;
116  }
117  double Chi2 = ROOT::Math::Similarity(r,W);// Chi2 = r.T()*W*r
118  double a_i = exp(-0.5*Chi2)/(2.*M_PI*sqrt(det));
119  mymap.push_back(std::pair<const TrackingRecHit*, float>((*ihit)->hit(), a_i));
120  double c_i = exp(-0.5*theChi2Cut/annealing)/(2.*M_PI*sqrt(det));
121  a_sum += a_i;
122  c_sum += c_i;
123  }
124  double total_sum = a_sum + c_sum;
125 
126  unsigned int counter = 0;
128  for(std::vector<TransientTrackingRecHit::RecHitPointer>::iterator ihit = updatedcomponents.begin(); ihit != updatedcomponents.end(); ihit++) {
129  //uncomment lines below to have like ORCA
130  double p = ((mymap[counter].second)/total_sum > 1.e-6 ? (mymap[counter].second)/total_sum : 1.e-6);
131  //float p = ((mymap[counter].second)/total_sum > 0.01 ? (mymap[counter].second)/total_sum : 1.e-6);
132  normmap.push_back(std::pair<const TrackingRecHit*, float>(mymap[counter].first, p));
133  //let's store the weight in the component TransientTrackingRecHit too
134  (*ihit)->setWeight(p);
135  (*ihit)->setAnnealingFactor(annealing);
136  finalcomponents.push_back(*ihit);
137  LogTrace("SiTrackerMultiRecHitUpdator")<< "Component hit type " << typeid(*mymap[counter].first).name()
138  << " position " << mymap[counter].first->localPosition()
139  << " error " << mymap[counter].first->localPositionError()
140  << " with weight " << p;
141  counter++;
142  }
143 
144  mymap = normmap;
145  // LocalError er = calcParametersError(finalcomponents);
146  // LocalPoint p = calcParameters(finalcomponents, er);
148  SiTrackerMultiRecHit updated(param.first, param.second, normmap.front().first->geographicalId(), normmap);
149  LogTrace("SiTrackerMultiRecHitUpdator") << "Updated Hit position " << updated.localPosition() << " updated error " << updated.parametersError() << std::endl;
150  //return new SiTrackerMultiRecHit(normmap);
151  return TSiTrackerMultiRecHit::build(geomdet, &updated, finalcomponents, annealing);
152 }
153 
154 
156  AlgebraicSymMatrix22 W_sum;
157  AlgebraicVector2 m_sum;
158  int ierr;
159  for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = map.begin(); ihit != map.end(); ihit ++) {
160  AlgebraicVector2 m(asSVector<2>((*ihit)->parameters()));
161  AlgebraicSymMatrix22 V(asSMatrix<2>((*ihit)->parametersError()));
162  AlgebraicSymMatrix22 W(V.Inverse(ierr));
163 
164  if(ierr != 0) {
165  edm::LogError("SiTrackerMultiRecHitUpdator")<<"MultiRecHit::checkParameters: W not valid!"<<std::endl;
166  }
167 
168  else {
169  W_sum += ((*ihit)->weight()*W);
170  m_sum += ((*ihit)->weight()*(W*m));
171  }
172  }
173  AlgebraicSymMatrix22 V_sum= W_sum.Inverse(ierr);
174  AlgebraicVector2 parameters = V_sum*m_sum;
175  LocalError error=LocalError(V_sum(0,0), V_sum(0,1), V_sum(1,1));
177  return std::make_pair(position,error);
178 }
179 
181  AlgebraicSymMatrix22 W_sum;
182  int ierr;
183  for(TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = map.begin(); ihit != map.end(); ihit ++) {
184  AlgebraicSymMatrix22 V(asSMatrix<2>((*ihit)->parametersError()));
185  AlgebraicSymMatrix22 W(V.Inverse(ierr));
186 
187  if(ierr != 0) {
188  edm::LogError("SiTrackerMultiRecHitUpdator")<<"MultiRecHit::checkParametersError: W not valid!"<<std::endl;
189  }
190 
191  else W_sum += ((*ihit)->weight()*W);
192  }
193  AlgebraicSymMatrix22 parametersError = W_sum.Inverse(ierr);
194  return LocalError(parametersError(0,0), parametersError(0,1), parametersError(1,1));
195 }
196 
198  AlgebraicVector2 m_sum;
199  int ierr;
200  for( TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = map.begin(); ihit != map.end(); ihit ++) {
201  AlgebraicVector2 m(asSVector<2>((*ihit)->parameters()));
202  AlgebraicSymMatrix22 V(asSMatrix<2>((*ihit)->parametersError()));
203  AlgebraicSymMatrix22 W(V.Inverse(ierr));
204 
205  if(ierr != 0) {
206  edm::LogError("SiTrackerMultiRecHitUpdator")<<"MultiRecHit::checkParameters: W not valid!"<<std::endl;
207  }
208 
209  //m_sum += ihit->weight()*(W*m);
210  else m_sum += ((*ihit)->weight()*(W*m));
211  }
212  AlgebraicSymMatrix22 V_sum;
213 
214  V_sum(0,0) = er.xx();
215  V_sum(0,1) = er.xy();
216  V_sum(1,1) = er.yy();
217  //AlgebraicSymMatrix V_sum(parametersError());
218  AlgebraicVector2 parameters = V_sum*m_sum;
219  return LocalPoint(parameters(0), parameters(1));
220 }
221 
float xx() const
Definition: LocalError.h:19
virtual TransientTrackingRecHit::RecHitPointer buildMultiRecHit(const std::vector< const TrackingRecHit * > &rhv, TrajectoryStateOnSurface tsos, float annealing=1.) const
ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2 > > AlgebraicSymMatrix22
static RecHitPointer build(const GeomDet *geom, Type type=TrackingRecHit::missing, const DetLayer *layer=0)
T y() const
Definition: PV3DBase.h:57
LocalPoint calcParameters(TransientTrackingRecHit::ConstRecHitContainer &map, const LocalError &er) const
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
std::pair< LocalPoint, LocalError > LocalParameters
const TrackingRecHitPropagator * theHitPropagator
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
dictionary map
Definition: Association.py:160
virtual TransientTrackingRecHit::RecHitPointer update(TransientTrackingRecHit::ConstRecHitPointer original, TrajectoryStateOnSurface tsos, double annealing=1.) const
U second(std::pair< T, U > const &p)
float xy() const
Definition: LocalError.h:20
float yy() const
Definition: LocalError.h:21
static RecHitPointer build(const GeomDet *geom, const SiTrackerMultiRecHit *rh, const ConstRecHitContainer &components, float annealing=1.)
T sqrt(T t)
Definition: SSEVec.h:28
virtual RecHitPointer build(const TrackingRecHit *p) const =0
build a tracking rechit from an existing rechit
const TransientTrackingRecHitBuilder * theBuilder
bool first
Definition: L1TdeRCT.cc:79
#define LogTrace(id)
std::vector< ConstRecHitPointer > ConstRecHitContainer
#define M_PI
Definition: BFit3D.cc:3
const Surface & surface() const
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
LocalError calcParametersError(TransientTrackingRecHit::ConstRecHitContainer &map) const
Definition: Chi2.h:17
T x() const
Definition: PV3DBase.h:56
const PositionType & position() const
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
ROOT::Math::SVector< double, 2 > AlgebraicVector2
SiTrackerMultiRecHitUpdator(const TransientTrackingRecHitBuilder *builder, const TrackingRecHitPropagator *hitpropagator, const float Chi2Cut, const std::vector< double > &anAnnealingProgram)
TransientTrackingRecHit::RecHitPointer project(const TransientTrackingRecHit::ConstRecHitPointer hit, const GeomDet &det, const TrajectoryStateOnSurface ts) const