SiTrackerMRHTools/src/SiTrackerMultiRecHitUpdator.cc

0001 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h"
0002 #include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h"
0003 #include "DataFormats/Math/interface/invertPosDefMatrix.h"
0004 #include "DataFormats/Math/interface/ProjectMatrix.h"
0005 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0006
0007 #include "RecoTracker/SiTrackerMRHTools/interface/GenericProjectedRecHit2D.h"
0008 #include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h"
0009 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0010
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0012 #include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
0013 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 SiTrackerMultiRecHitUpdator::SiTrackerMultiRecHitUpdator(const TransientTrackingRecHitBuilder* builder,
0018                                                          const TrackingRecHitPropagator* hitpropagator,
0019                                                          const float Chi2Cut1D,
0020                                                          const float Chi2Cut2D,
0021                                                          const std::vector<double>& anAnnealingProgram,
0022                                                          bool debug)
0023     : theBuilder(builder),
0024       theHitPropagator(hitpropagator),
0025       theChi2Cut1D(Chi2Cut1D),
0026       theChi2Cut2D(Chi2Cut2D),
0027       theAnnealingProgram(anAnnealingProgram),
0028       debug_(debug) {
0029   theHitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(builder)->cloner();
0030 }
0031
0032 TransientTrackingRecHit::RecHitPointer SiTrackerMultiRecHitUpdator::buildMultiRecHit(
0033     const std::vector<const TrackingRecHit*>& rhv,
0034     const TrajectoryStateOnSurface& tsos,
0035     MeasurementDetWithData& measDet,
0036     float annealing) const {
0037   LogTrace("SiTrackerMultiRecHitUpdator")
0038       << "Calling SiTrackerMultiRecHitUpdator::buildMultiRecHit with AnnealingFactor: " << annealing;
0039
0040   TransientTrackingRecHit::ConstRecHitContainer tcomponents;
0041   for (std::vector<const TrackingRecHit*>::const_iterator iter = rhv.begin(); iter != rhv.end(); iter++) {
0042     TransientTrackingRecHit::RecHitPointer transient = theBuilder->build(*iter);
0043     if (transient->isValid())
0044       tcomponents.push_back(transient);
0045   }
0046   return update(tcomponents, tsos, measDet, annealing);
0047 }
0048
0049 TransientTrackingRecHit::RecHitPointer SiTrackerMultiRecHitUpdator::update(
0050     TransientTrackingRecHit::ConstRecHitPointer original,
0051     const TrajectoryStateOnSurface& tsos,
0052     MeasurementDetWithData& measDet,
0053     double annealing) const {
0054   LogTrace("SiTrackerMultiRecHitUpdator")
0055       << "Calling SiTrackerMultiRecHitUpdator::update with AnnealingFactor: " << annealing;
0056
0057   if (!tsos.isValid()) {
0058     //return original->clone();
0059     throw cms::Exception("SiTrackerMultiRecHitUpdator") << "!!! MultiRecHitUpdator::update(..): tsos NOT valid!!! ";
0060   }
0061
0062   //check if to clone is the right thing
0063   if (original->isValid()) {
0064     if (original->transientHits().empty()) {
0065       return theHitCloner.makeShared(original, tsos);
0066     }
0067   } else {
0068     return theHitCloner.makeShared(original, tsos);
0069   }
0070
0071   TransientTrackingRecHit::ConstRecHitContainer tcomponents = original->transientHits();
0072   return update(tcomponents, tsos, measDet, annealing);
0073 }
0074
0075 /*------------------------------------------------------------------------------------------------------------------------*/
0076 TransientTrackingRecHit::RecHitPointer SiTrackerMultiRecHitUpdator::update(
0077     TransientTrackingRecHit::ConstRecHitContainer& tcomponents,
0078     const TrajectoryStateOnSurface& tsos,
0079     MeasurementDetWithData& measDet,
0080     double annealing) const {
0081   if (tcomponents.empty()) {
0082     LogTrace("SiTrackerMultiRecHitUpdator") << "Empty components vector passed to SiTrackerMultiRecHitUpdator::update, "
0083                                                "returning an InvalidTransientRecHit ";
0084     return std::make_shared<InvalidTrackingRecHit>(measDet.mdet().geomDet(), TrackingRecHit::missing);
0085   }
0086
0087   if (!tsos.isValid()) {
0088     LogTrace("SiTrackerMultiRecHitUpdator")
0089         << "SiTrackerMultiRecHitUpdator::update: tsos NOT valid!!!, returning an InvalidTransientRecHit";
0090     return std::make_shared<InvalidTrackingRecHit>(measDet.mdet().geomDet(), TrackingRecHit::missing);
0091   }
0092
0093   std::vector<TransientTrackingRecHit::RecHitPointer> updatedcomponents;
0094   const GeomDet* geomdet = nullptr;
0095
0096   //running on all over the MRH components
0097   for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iter = tcomponents.begin();
0098        iter != tcomponents.end();
0099        iter++) {
0100     //the first rechit must belong to the same surface of TSOS
0101     if (iter == tcomponents.begin()) {
0102       if (&((*iter)->det()->surface()) != &(tsos.surface())) {
0103         throw cms::Exception("SiTrackerMultiRecHitUpdator")
0104             << "the Trajectory state and the first rechit "
0105                "passed to the SiTrackerMultiRecHitUpdator lay on different surfaces!: state lays on surface "
0106             << tsos.surface().position() << " hit with detid " << (*iter)->det()->geographicalId().rawId()
0107             << " lays on surface " << (*iter)->det()->surface().position();
0108       }
0109
0110       geomdet = (*iter)->det();
0111     }
0112
0113     //if the rechit does not belong to the surface of the tsos
0114     //GenericProjectedRecHit2D is used to prepagate
0115     if (&((*iter)->det()->surface()) != &(tsos.surface())) {
0116       TransientTrackingRecHit::RecHitPointer cloned =
0117           theHitPropagator->project<GenericProjectedRecHit2D>(*iter, *geomdet, tsos, theBuilder);
0118       //if it is used a sensor by sensor grouping this should not appear
0119       if (cloned->isValid())
0120         updatedcomponents.push_back(cloned);
0121
0122     } else {
0123       TransientTrackingRecHit::RecHitPointer cloned = theHitCloner.makeShared(*iter, tsos);
0124       if (cloned->isValid()) {
0125         updatedcomponents.push_back(cloned);
0126       }
0127     }
0128   }
0129
0130   std::vector<std::pair<const TrackingRecHit*, float> > mymap;
0131   std::vector<std::pair<const TrackingRecHit*, float> > normmap;
0132
0133   double a_sum = 0, c_sum = 0;
0134
0135   for (std::vector<TransientTrackingRecHit::RecHitPointer>::iterator ihit = updatedcomponents.begin();
0136        ihit != updatedcomponents.end();
0137        ihit++) {
0138     double a_i = ComputeWeight(tsos, *(*ihit), false, annealing);  //exp(-0.5*Chi2)
0139     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t a_i:" << a_i;
0140     //double c_i = ComputeWeight(tsos, *(*ihit), true, annealing);  //exp(-0.5*theChi2Cut/annealing)/(2.*M_PI*sqrt(det));
0141     //LogTrace("SiTrackerMultiRecHitUpdator")<< "\t\t c_i:" << c_i ;
0142     mymap.push_back(std::pair<const TrackingRecHit*, float>((*ihit)->hit(), a_i));
0143
0144     a_sum += a_i;
0145     //with the new definition, the cut weight is computed only once
0146     if (ihit == updatedcomponents.begin())
0147       c_sum = ComputeWeight(tsos, *(*ihit), true, annealing);  //exp(-0.5*theChi2Cut/annealing)
0148   }
0149   double total_sum = a_sum + c_sum;
0150   LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t c_sum:" << c_sum;
0151   LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t total sum:" << total_sum;
0152
0153   unsigned int counter = 0;
0154   bool invalid = true;
0155   for (std::vector<TransientTrackingRecHit::RecHitPointer>::iterator ihit = updatedcomponents.begin();
0156        ihit != updatedcomponents.end();
0157        ihit++) {
0158     double p = ((mymap[counter].second) / total_sum > 1.e-12 ? (mymap[counter].second) / total_sum : 1.e-12);
0159     //ORCA: float p = ((mymap[counter].second)/total_sum > 0.01 ? (mymap[counter].second)/total_sum : 1.e-6);
0160
0161 #ifdef EDM_ML_DEBUG
0162     auto const& tmp = *mymap[counter].first;
0163     LogTrace("SiTrackerMultiRecHitUpdator")
0164         << "  Component hit type " << typeid(tmp).name() << " and dim:" << tmp.dimension() << " position (PRECISE!!!)"
0165         << tmp.localPosition() << " error " << tmp.localPositionError() << " with weight " << p;
0166 #endif
0167
0168     if (p > 10e-6) {
0169       invalid = false;
0170       normmap.push_back(std::pair<const TrackingRecHit*, float>(mymap[counter].first, p));
0171     }
0172
0173     counter++;
0174   }
0175
0176   if (!invalid) {
0177     SiTrackerMultiRecHitUpdator::LocalParameters param = calcParameters(tsos, normmap);
0178     SiTrackerMultiRecHit updated(param.first, param.second, *normmap.front().first->det(), normmap, annealing);
0179     LogTrace("SiTrackerMultiRecHitUpdator") << " Updated Hit position " << updated.localPosition() << " updated error "
0180                                             << updated.localPositionError() << std::endl;
0181
0182     return std::make_shared<SiTrackerMultiRecHit>(
0183         param.first, param.second, *normmap.front().first->det(), normmap, annealing);
0184
0185   } else {
0186     LogTrace("SiTrackerMultiRecHitUpdator")
0187         << " No hits with weight (> 10e-6) have been found for this MRH." << std::endl;
0188     return std::make_shared<InvalidTrackingRecHit>(measDet.mdet().geomDet(), TrackingRecHit::missing);
0189   }
0190 }
0191
0192 //---------------------------------------------------------------------------------------------------------------
0193 double SiTrackerMultiRecHitUpdator::ComputeWeight(const TrajectoryStateOnSurface& tsos,
0194                                                   const TransientTrackingRecHit& aRecHit,
0195                                                   bool CutWeight,
0196                                                   double annealing) const {
0197   switch (aRecHit.dimension()) {
0198     case 1:
0199       return ComputeWeight<1>(tsos, aRecHit, CutWeight, annealing);
0200     case 2:
0201       return ComputeWeight<2>(tsos, aRecHit, CutWeight, annealing);
0202     case 3:
0203       return ComputeWeight<3>(tsos, aRecHit, CutWeight, annealing);
0204     case 4:
0205       return ComputeWeight<4>(tsos, aRecHit, CutWeight, annealing);
0206     case 5:
0207       return ComputeWeight<5>(tsos, aRecHit, CutWeight, annealing);
0208   }
0209   throw cms::Exception("Rec hit of invalid dimension (not 1,2,3,4,5)")
0210       << "The value was " << aRecHit.dimension() << ", type is " << typeid(aRecHit).name() << "\n";
0211 }
0212
0213 //---------------------------------------------------------------------------------------------------------------
0214 template <unsigned int N>
0215 double SiTrackerMultiRecHitUpdator::ComputeWeight(const TrajectoryStateOnSurface& tsos,
0216                                                   const TransientTrackingRecHit& aRecHit,
0217                                                   bool CutWeight,
0218                                                   double annealing) const {
0219   //  typedef typename AlgebraicROOTObject<N,5>::Matrix MatN5;
0220   //  typedef typename AlgebraicROOTObject<5,N>::Matrix Mat5N;
0221   //  typedef typename AlgebraicROOTObject<N,N>::SymMatrix SMatNN;
0222   //  typedef typename AlgebraicROOTObject<N>::Vector VecN;
0223
0224   // define variables that will be used to setup the KfComponentsHolder
0225   ProjectMatrix<double, 5, N> pf;
0226   typename AlgebraicROOTObject<N>::Vector r, rMeas;
0227   typename AlgebraicROOTObject<N, N>::SymMatrix R, RMeas, W;
0228   AlgebraicVector5 x = tsos.localParameters().vector();
0229   const AlgebraicSymMatrix55& C = (tsos.localError().matrix());
0230
0231   // setup the holder with the correct dimensions and get the values
0232   KfComponentsHolder holder;
0233   holder.template setup<N>(&r, &R, &pf, &rMeas, &RMeas, x, C);
0234   aRecHit.getKfComponents(holder);
0235
0236   typename AlgebraicROOTObject<N>::Vector diff = r - rMeas;
0237
0238   if (!CutWeight) {
0239     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t r:" << r;
0240     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t tsospos:" << rMeas;
0241     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t diff:" << diff;
0242     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t R:" << R;
0243     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t RMeas:" << RMeas;
0244   }
0245
0246   R += RMeas;  //assume that TSOS is predicted || comb one
0247   if (!CutWeight)
0248     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t R+RMeas:" << R;
0249
0250   //computing chi2 with the smoothTsos
0251   // SMatNN R_smooth = R - RMeas;
0252   // if(!CutWeight)  LogTrace("SiTrackerMultiRecHitUpdator")<< "\t\t R-=Rmeas:" << R_smooth ;
0253   // bool ierr2_bis = invertPosDefMatrix(R_smooth);
0254   // double Chi2_smooth = ROOT::Math::Similarity(diff, R_smooth);
0255
0256   //computing chi2 with the smoothTsos
0257   //SMatNN R_pred   = R + RMeas;
0258   //if(!CutWeight)  LogTrace("SiTrackerMultiRecHitUpdator")<< "\t\t R+=Rmeas:" << R_pred   ;
0259   //bool ierr2_bis = invertPosDefMatrix(R_pred  );
0260   //double Chi2_pred   = ROOT::Math::Similarity(diff, R_pred  );
0261
0262   //Det2 method will preserve the content of the Matrix
0263   //and return true when the calculation is successfull
0264   double det;
0265   bool ierr = R.Det2(det);
0266
0267   bool ierr2 = invertPosDefMatrix(R);  //ierr will be set to true when inversion is successfull
0268   double Chi2 = ROOT::Math::Similarity(diff, R);
0269
0270   if (!ierr || !ierr2) {
0271     LogTrace("SiTrackerMultiRecHitUpdator") << "SiTrackerMultiRecHitUpdator::ComputeWeight: W not valid!" << std::endl;
0272     LogTrace("SiTrackerMultiRecHitUpdator") << "V: " << R << " AnnealingFactor: " << annealing << std::endl;
0273   }
0274
0275   if (!CutWeight) {
0276     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t det:" << det;
0277     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t Chi2:" << Chi2;
0278     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t Chi2/ann:" << Chi2 / annealing;
0279   }
0280
0281   double temp_weight = 0.0;
0282   if (CutWeight && N == 1) {
0283     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t Chi2Cut1D:" << theChi2Cut1D;
0284     temp_weight = exp(-0.5 * theChi2Cut1D / annealing);
0285   } else if (CutWeight && N == 2) {
0286     LogTrace("SiTrackerMultiRecHitUpdator") << "\t\t Chi2Cut2D:" << theChi2Cut2D;
0287     temp_weight = exp(-0.5 * theChi2Cut2D / annealing);
0288   }
0289
0290   if (!CutWeight) {
0291     temp_weight = exp(-0.5 * Chi2 / annealing);
0292   }
0293
0294   return temp_weight;
0295 }
0296
0297 //-----------------------------------------------------------------------------------------------------------
0298 SiTrackerMultiRecHitUpdator::LocalParameters SiTrackerMultiRecHitUpdator::calcParameters(
0299     const TrajectoryStateOnSurface& tsos, std::vector<std::pair<const TrackingRecHit*, float> >& aHitMap) const {
0300   //supposing all the hits inside of a MRH have the same dimension
0301   LogTrace("SiTrackerMultiRecHitUpdator")
0302       << "SiTrackerMultiRecHitUpdator::LocalParameters: dim first recHit: " << aHitMap[0].first->dimension()
0303       << std::endl;
0304   switch (aHitMap[0].first->dimension()) {
0305     case 1:
0306       return calcParameters<1>(tsos, aHitMap);
0307     case 2:
0308       return calcParameters<2>(tsos, aHitMap);
0309   }
0310   throw cms::Exception("Rec hit of invalid dimension for computing MRH (not 1,2)")
0311       << "The value was " << aHitMap[0].first->dimension() << ", type is " << typeid(aHitMap[0].first).name() << "\n";
0312 }
0313
0314 //-----------------------------------------------------------------------------------------------------------
0315 template <unsigned int N>
0316 SiTrackerMultiRecHitUpdator::LocalParameters SiTrackerMultiRecHitUpdator::calcParameters(
0317     const TrajectoryStateOnSurface& tsos, std::vector<std::pair<const TrackingRecHit*, float> >& aHitMap) const {
0318   typedef typename AlgebraicROOTObject<N, N>::SymMatrix SMatNN;
0319   typedef typename AlgebraicROOTObject<N>::Vector VecN;
0320
0321   VecN m_sum;
0322   SMatNN W_sum;
0323   LocalPoint position;
0324   LocalError error;
0325
0326   //for TID and TEC the correlation is really high -> need to be scorrelated and then correlated again
0327   float s = 0.1;
0328
0329   for (std::vector<std::pair<const TrackingRecHit*, float> >::const_iterator ihit = aHitMap.begin();
0330        ihit != aHitMap.end();
0331        ihit++) {
0332     // define variables that will be used to setup the KfComponentsHolder
0333     ProjectMatrix<double, 5, N> pf;
0334     typename AlgebraicROOTObject<N>::Vector r, rMeas;
0335     typename AlgebraicROOTObject<N, N>::SymMatrix V, VMeas, Wtemp;
0336     AlgebraicVector5 x = tsos.localParameters().vector();
0337     const AlgebraicSymMatrix55& C = (tsos.localError().matrix());
0338
0339     // setup the holder with the correct dimensions and get the values
0340     KfComponentsHolder holder;
0341     holder.template setup<N>(&r, &V, &pf, &rMeas, &VMeas, x, C);
0342     (ihit->first)->getKfComponents(holder);
0343
0344     LogTrace("SiTrackerMultiRecHitUpdator") << "\t position: " << r;
0345     LogTrace("SiTrackerMultiRecHitUpdator") << "\t error: " << V;
0346
0347     //scorrelation  in TID and TEC
0348     if (N == 2 && TIDorTEChit(ihit->first)) {
0349       V(0, 1) = V(1, 0) = V(0, 1) * s;
0350       //      V(1,0) = V(1,0)*s;
0351       LogTrace("SiTrackerMultiRecHitUpdator") << "\t error scorr: " << V;
0352     }
0353
0354     bool ierr = invertPosDefMatrix(V);
0355     if (!ierr) {
0356       edm::LogError("SiTrackerMultiRecHitUpdator") << "SiTrackerMultiRecHitUpdator::calcParameters: V not valid!";
0357     }
0358     LogTrace("SiTrackerMultiRecHitUpdator") << "\t inverse error: " << V;
0359
0360     //compute m_sum and W_sum
0361     m_sum += (ihit->second * V * r);
0362     W_sum += (ihit->second * V);
0363   }
0364
0365   bool ierr_sum = invertPosDefMatrix(W_sum);
0366   if (!ierr_sum) {
0367     edm::LogError("SiTrackerMultiRecHitUpdator") << "SiTrackerMultiRecHitUpdator::calcParameters: W_sum not valid!";
0368   }
0369
0370   LogTrace("SiTrackerMultiRecHitUpdator") << "\t inverse total error: " << W_sum;
0371   typename AlgebraicROOTObject<N>::Vector parameters = W_sum * m_sum;
0372   if (N == 1) {
0373     position = LocalPoint(parameters(0), 0.f);
0374     error = LocalError(W_sum(0, 0), 0.f, std::numeric_limits<float>::max());
0375   } else if (N == 2) {
0376     position = LocalPoint(parameters(0), parameters(1));
0377     //ri-correlation  in TID and TEC
0378     if (TIDorTEChit(aHitMap.at(0).first))
0379       error = LocalError(W_sum(0, 0), W_sum(0, 1) / s, W_sum(1, 1));
0380     else
0381       error = LocalError(W_sum(0, 0), W_sum(0, 1), W_sum(1, 1));
0382   }
0383
0384   return std::make_pair(position, error);
0385 }
0386
0387 bool SiTrackerMultiRecHitUpdator::TIDorTEChit(const TrackingRecHit* const& hit) const {
0388   DetId hitId = hit->geographicalId();
0389
0390   if (hitId.det() == DetId::Tracker &&
0391       (hitId.subdetId() == StripSubdetector::TEC || hitId.subdetId() == StripSubdetector::TID)) {
0392     return true;
0393   }
0394
0395   return false;
0396 }