CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/MuonAnalysis/MuonAssociators/src/MatcherUsingTracksAlgorithm.cc

Go to the documentation of this file.
00001 #include "MuonAnalysis/MuonAssociators/interface/MatcherUsingTracksAlgorithm.h"
00002 
00003 #include "DataFormats/Math/interface/deltaR.h"
00004 #include "DataFormats/Math/interface/deltaPhi.h"
00005 
00006 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00007 #include "DataFormats/TrackReco/interface/Track.h"
00008 
00009 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00011 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00013 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
00014 #include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"
00015 
00016 #include "FWCore/Utilities/interface/Exception.h"
00017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00018 
00019 
00020 MatcherUsingTracksAlgorithm::MatcherUsingTracksAlgorithm(const edm::ParameterSet & iConfig) :
00021     whichTrack1_(None),     whichTrack2_(None), 
00022     whichState1_(AtVertex), whichState2_(AtVertex),
00023     srcCut_(iConfig.existsAs<std::string>("srcPreselection") ? iConfig.getParameter<std::string>("srcPreselection") : ""),
00024     matchedCut_(iConfig.existsAs<std::string>("matchedPreselection") ? iConfig.getParameter<std::string>("matchedPreselection") : ""),
00025     requireSameCharge_(iConfig.existsAs<bool>("requireSameCharge") ? iConfig.getParameter<bool>("requireSameCharge") : false)
00026 {
00027     std::string algo = iConfig.getParameter<std::string>("algorithm");
00028     if      (algo == "byTrackRef")           { algo_ = ByTrackRef; }
00029     else if (algo == "byPropagatingSrc")     { algo_ = ByPropagatingSrc; }
00030     else if (algo == "byPropagatingMatched") { algo_ = ByPropagatingMatched; }
00031     else if (algo == "byDirectComparison")   { algo_ = ByDirectComparison; }
00032     else throw cms::Exception("Configuration") << "Value '" << algo << "' for algorithm not yet implemented.\n";
00033 
00034     getConf(iConfig, "src",     whichTrack1_, whichState1_);
00035     getConf(iConfig, "matched", whichTrack2_, whichState2_);
00036 
00037     if (algo_ == ByTrackRef) {
00038         // validate the config
00039         if (whichTrack1_ == None || whichTrack2_ == None) throw cms::Exception("Configuration") << "Algorithm 'byTrackRef' needs tracks not to be 'none'.\n";
00040     } else if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingMatched || algo_ == ByDirectComparison) {
00041         // read matching cuts
00042         maxLocalPosDiff_      = iConfig.getParameter<double>("maxDeltaLocalPos");
00043         maxGlobalMomDeltaR_   = iConfig.getParameter<double>("maxDeltaR");
00044         maxGlobalMomDeltaEta_ = iConfig.existsAs<double>("maxDeltaEta") ? iConfig.getParameter<double>("maxDeltaEta") : maxGlobalMomDeltaR_;
00045         maxGlobalMomDeltaPhi_ = iConfig.existsAs<double>("maxDeltaPhi") ? iConfig.getParameter<double>("maxDeltaPhi") : maxGlobalMomDeltaR_;
00046         maxGlobalDPtRel_      = iConfig.getParameter<double>("maxDeltaPtRel");
00047 
00048         // choice of sorting variable
00049         std::string sortBy = iConfig.getParameter<std::string>("sortBy");
00050         if      (sortBy == "deltaLocalPos") sortBy_ = LocalPosDiff;
00051         else if (sortBy == "deltaPtRel")    sortBy_ = GlobalDPtRel;
00052         else if (sortBy == "deltaR")        sortBy_ = GlobalMomDeltaR;
00053         else if (sortBy == "deltaEta")      sortBy_ = GlobalMomDeltaEta;
00054         else if (sortBy == "deltaPhi")      sortBy_ = GlobalMomDeltaPhi;
00055         else if (sortBy == "chi2")          sortBy_ = Chi2;
00056         else throw cms::Exception("Configuration") << "Parameter 'sortBy' must be one of: deltaLocalPos, deltaPtRel, deltaR, chi2.\n";
00057         // validate the config
00058         if (algo_ == ByPropagatingSrc) {
00059             if (whichTrack2_ == None || whichState2_ == AtVertex) {
00060                 algo_ = ByPropagatingSrcTSCP;
00061                 //throw cms::Exception("Configuration") << "Destination track must be non-null, and state must not be 'AtVertex' (not yet).\n";
00062             }
00063         } else if (algo_ == ByPropagatingMatched) {
00064             if (whichTrack1_ == None || whichState1_ == AtVertex) {
00065                 algo_ = ByPropagatingMatchedTSCP;
00066                 //throw cms::Exception("Configuration") << "Destination track must be non-null, and state must not be 'AtVertex' (not yet).\n";
00067             }
00068         } else if (algo_ == ByDirectComparison) {
00069             bool firstAtVertex = (whichTrack1_ == None || whichState1_ == AtVertex);
00070             bool   secAtVertex = (whichTrack2_ == None || whichState2_ == AtVertex);
00071             if (firstAtVertex) {
00072                 if (!secAtVertex)  throw cms::Exception("Configuration") << "When using 'byDirectComparison' with 'src' at vertex (or None), 'matched' must be at vertex too.\n";
00073             } else {
00074                 if (secAtVertex)  throw cms::Exception("Configuration") << "When using 'byDirectComparison' with 'src' not at vertex, 'matched' can't be at vertex or None.\n";
00075                 if (whichState1_ != whichState2_) throw cms::Exception("Configuration") << "You can't use 'byDirectComparison' with non-matching states.\n";
00076             }
00077         }
00078 
00079         useChi2_ = iConfig.existsAs<bool>("computeChi2") ? iConfig.getParameter<bool>("computeChi2") : false;
00080         if (useChi2_) {
00081            if (whichTrack1_ == None && whichTrack2_ == None) throw cms::Exception("Configuration") << "Can't compute chi2s if both tracks are set to 'none'.\n";
00082            maxChi2_ = iConfig.getParameter<double>("maxChi2");
00083            chi2DiagonalOnly_ = iConfig.getParameter<bool>("chi2DiagonalOnly");
00084            if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingMatched) {
00085                chi2UseVertex_ = iConfig.getParameter<bool>("chi2UsePosition");
00086            } else {
00087                chi2UseVertex_ = iConfig.getParameter<bool>("chi2UseVertex");
00088                if (algo_ == ByDirectComparison) {
00089                     std::string choice = iConfig.getParameter<std::string>("chi2MomentumForDxy");
00090                     if (choice == "src") chi2FirstMomentum_ = true;
00091                     else if (choice != "matched") throw cms::Exception("Configuration") << "chi2MomentumForDxy must be 'src' or 'matched'\n";
00092                }
00093            }
00094         } else maxChi2_ = 1;
00095 
00096         if (sortBy_ == Chi2 && !useChi2_) throw cms::Exception("Configuration") << "Can't sort by chi2s if 'computeChi2s' is not set to true.\n";
00097     }
00098 }
00099 
00100 void 
00101 MatcherUsingTracksAlgorithm::getConf(const edm::ParameterSet & iConfig, const std::string &whatFor, WhichTrack &whichTrack, WhichState &whichState) {
00102     std::string s_whichTrack = iConfig.getParameter<std::string>(whatFor+"Track");
00103     if      (s_whichTrack == "none")    { whichTrack = None; }
00104     else if (s_whichTrack == "tracker") { whichTrack = TrackerTk; }
00105     else if (s_whichTrack == "muon")    { whichTrack = MuonTk; }
00106     else if (s_whichTrack == "global")  { whichTrack = GlobalTk; }
00107     else throw cms::Exception("Configuration") << "Parameter 'useTrack' must be 'none', 'tracker', 'muon', 'global'\n";
00108     if ((whichTrack != None) && (algo_ != ByTrackRef)) {
00109         std::string s_whichState = iConfig.getParameter<std::string>(whatFor+"State");
00110         if      (s_whichState == "atVertex")  { whichState = AtVertex; }
00111         else if (s_whichState == "innermost") { whichState = Innermost; }
00112         else if (s_whichState == "outermost") { whichState = Outermost; }
00113         else throw cms::Exception("Configuration") << "Parameter 'useState' must be 'atVertex', 'innermost', 'outermost'\n";
00114     }
00115 }
00116 
00119 bool 
00120 MatcherUsingTracksAlgorithm::match(const reco::Candidate &c1, const reco::Candidate &c2, float &deltR, float &deltEta, float &deltPhi, float &deltaLocalPos, float &deltaPtRel, float &chi2) const {
00121     if (!(srcCut_(c1) && matchedCut_(c2))) return false;
00122     if (requireSameCharge_ && (c1.charge() != c2.charge())) return false;
00123     switch (algo_) {
00124         case ByTrackRef: { 
00125             reco::TrackRef t1 = getTrack(c1, whichTrack1_); 
00126             reco::TrackRef t2 = getTrack(c2, whichTrack2_); 
00127             if (t1.isNonnull()) {
00128                 if (t1 == t2) return true; 
00129                 if (t1.id() != t2.id()) { edm::LogWarning("MatcherUsingTracksAlgorithm") << "Trying to match by reference tracks coming from different collections.\n"; }
00130             } 
00131             }
00132         case ByPropagatingSrc: {
00133             FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
00134             TrajectoryStateOnSurface target = targetState(c2, whichTrack2_, whichState2_);
00135             return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00136             }
00137         case ByPropagatingMatched: { 
00138             FreeTrajectoryState start = startingState(c2, whichTrack2_, whichState2_);
00139             TrajectoryStateOnSurface target = targetState(c1, whichTrack1_, whichState1_);
00140             return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00141             }
00142         case ByPropagatingSrcTSCP: {
00143             FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
00144             FreeTrajectoryState target = startingState(c2, whichTrack2_, whichState2_);
00145             return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00146             }
00147         case ByPropagatingMatchedTSCP: { 
00148             FreeTrajectoryState start = startingState(c2, whichTrack2_, whichState2_);
00149             FreeTrajectoryState target = startingState(c1, whichTrack1_, whichState1_);
00150             return matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00151             }
00152 
00153         case ByDirectComparison: {
00154             FreeTrajectoryState start = startingState(c1, whichTrack1_, whichState1_);
00155             FreeTrajectoryState otherstart = startingState(c2, whichTrack2_, whichState2_);
00156             return matchByDirectComparison(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2);
00157             }
00158     }
00159     return false;
00160 }
00161 
00165 int 
00166 MatcherUsingTracksAlgorithm::match(const reco::Candidate &c1, const edm::View<reco::Candidate> &c2s, float &deltR, float &deltEta, float &deltPhi, float &deltaLocalPos, float &deltaPtRel, float &chi2) const {
00167     if (!srcCut_(c1)) return -1;
00168     
00169     // working and output variables
00170     FreeTrajectoryState start; TrajectoryStateOnSurface target;
00171     int match = -1;
00172 
00173     // pre-fetch some states if needed
00174     if (algo_ == ByPropagatingSrc || algo_ == ByPropagatingSrcTSCP || algo_ == ByPropagatingMatchedTSCP || algo_ == ByDirectComparison) {
00175         start = startingState(c1, whichTrack1_, whichState1_);
00176     } else if (algo_ == ByPropagatingMatched) target = targetState(c1, whichTrack1_, whichState1_);
00177 
00178     // loop on the  collection
00179     edm::View<reco::Candidate>::const_iterator it, ed; int i;
00180     for (it = c2s.begin(), ed = c2s.end(), i = 0; it != ed; ++it, ++i) {
00181         if (!matchedCut_(*it)) continue;
00182         if (requireSameCharge_ && (c1.charge() != it->charge()) ) continue;
00183         bool exit = false;
00184         switch (algo_) {
00185             case ByTrackRef: {
00186                      reco::TrackRef t1 = getTrack(c1, whichTrack1_); 
00187                      reco::TrackRef t2 = getTrack(*it, whichTrack2_); 
00188                      if (t1.isNonnull()) {
00189                          if (t1 == t2) { match = i; exit = true; }
00190                          if (t1.id() != t2.id()) { edm::LogWarning("MatcherUsingTracksAlgorithm") << "Trying to match by reference tracks coming from different collections.\n"; }
00191                      } 
00192                  } break;
00193             case ByPropagatingSrc: 
00194             case ByPropagatingMatched: {
00195                      if (algo_ == ByPropagatingMatched)  start = startingState(*it, whichTrack2_, whichState2_);
00196                      else if (algo_ == ByPropagatingSrc) target = targetState(*it, whichTrack2_, whichState2_);
00197                      if (matchWithPropagation(start, target, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) { 
00198                          match = i; 
00199                      }
00200                  } break;
00201             case ByDirectComparison: {
00202                      FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
00203                      if (matchByDirectComparison(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
00204                          match = i;
00205                      }
00206                  } break;
00207             case ByPropagatingSrcTSCP: {
00208                     FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
00209                     if (matchWithPropagation(start, otherstart, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
00210                          match = i;
00211                     }
00212                 } break;
00213             case ByPropagatingMatchedTSCP: {
00214                     FreeTrajectoryState otherstart = startingState(*it, whichTrack2_, whichState2_);
00215                     if (matchWithPropagation(otherstart, start, deltR, deltEta, deltPhi, deltaLocalPos, deltaPtRel, chi2)) {
00216                          match = i;
00217                     }
00218                 } break;
00219         }
00220         if (exit) break;
00221     }
00222 
00223     return match;
00224 }
00225 
00226 void 
00227 MatcherUsingTracksAlgorithm::init(const edm::EventSetup & iSetup) {
00228     iSetup.get<IdealMagneticFieldRecord>().get(magfield_);
00229     iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator_);
00230     iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_);
00231 }
00232 
00233 
00234 reco::TrackRef
00235 MatcherUsingTracksAlgorithm::getTrack(const reco::Candidate &reco, WhichTrack whichTrack) const {
00236     reco::TrackRef tk;
00237     const reco::RecoCandidate *rc = dynamic_cast<const reco::RecoCandidate *>(&reco);
00238     if (rc == 0) throw cms::Exception("Invalid Data") << "Input object is not a RecoCandidate.\n";
00239     switch (whichTrack) {
00240         case TrackerTk: tk = rc->track();          break; 
00241         case MuonTk   : tk = rc->standAloneMuon(); break;
00242         case GlobalTk : tk = rc->combinedMuon();   break;
00243         default: break; // just to make gcc happy
00244     }
00245     return tk;
00246 }
00247 
00248 FreeTrajectoryState 
00249 MatcherUsingTracksAlgorithm::startingState(const reco::Candidate &reco, WhichTrack whichTrack, WhichState whichState) const {
00250     FreeTrajectoryState ret;
00251     if (whichTrack != None) {
00252         reco::TrackRef tk = getTrack(reco, whichTrack);
00253         if (tk.isNull()) {
00254             ret = FreeTrajectoryState();
00255         } else {
00256             switch (whichState) {
00257                 case AtVertex:  ret = trajectoryStateTransform::initialFreeState(*tk, magfield_.product()); break;
00258                 case Innermost: ret = trajectoryStateTransform::innerFreeState(  *tk, magfield_.product()); break;
00259                 case Outermost: ret = trajectoryStateTransform::outerFreeState(  *tk, magfield_.product()); break;
00260             }
00261         }
00262     } else {
00263         ret = FreeTrajectoryState(  GlobalPoint( reco.vx(), reco.vy(), reco.vz()),
00264                                     GlobalVector(reco.px(), reco.py(), reco.pz()),
00265                                     reco.charge(),
00266                                     magfield_.product());
00267     }
00268     return ret;
00269 }
00270 
00271 TrajectoryStateOnSurface 
00272 MatcherUsingTracksAlgorithm::targetState(const reco::Candidate &reco, WhichTrack whichTrack, WhichState whichState) const {
00273     TrajectoryStateOnSurface ret;
00274     reco::TrackRef tk = getTrack(reco, whichTrack);
00275     if (tk.isNonnull()) {
00276         switch (whichState) {
00277             case Innermost: ret = trajectoryStateTransform::innerStateOnSurface(  *tk, *geometry_, magfield_.product()); break;
00278             case Outermost: ret = trajectoryStateTransform::outerStateOnSurface(  *tk, *geometry_, magfield_.product()); break;
00279             default: break; // just to make gcc happy
00280         }
00281     }
00282     return ret;
00283 }
00284 
00285 bool
00286 MatcherUsingTracksAlgorithm::matchWithPropagation(const FreeTrajectoryState &start, 
00287                                               const TrajectoryStateOnSurface &target, 
00288                                               float &lastDeltaR,
00289                                               float &lastDeltaEta,
00290                                               float &lastDeltaPhi, 
00291                                               float &lastDeltaLocalPos,
00292                                               float &lastGlobalDPtRel,
00293                                               float &lastChi2) const 
00294 {
00295     if ((start.momentum().mag() == 0) || !target.isValid()) return false;
00296 
00297     TrajectoryStateOnSurface tsos = propagator_->propagate(start, target.surface());
00298 
00299     bool isBest = false;
00300     if (tsos.isValid()) {
00301         float thisLocalPosDiff = (tsos.localPosition()-target.localPosition()).mag();
00302         float thisGlobalMomDeltaR = deltaR(tsos.globalMomentum(), target.globalMomentum());
00303         float thisGlobalMomDeltaPhi = fabs(deltaPhi(tsos.globalMomentum().phi(), target.globalMomentum().phi()));
00304         float thisGlobalMomDeltaEta = fabs(tsos.globalMomentum().eta() - target.globalMomentum().eta());
00305         float thisGlobalDPtRel = (tsos.globalMomentum().perp() - target.globalMomentum().perp())/target.globalMomentum().perp();
00306 
00307         if ((thisLocalPosDiff       < maxLocalPosDiff_) &&
00308             (thisGlobalMomDeltaR    < maxGlobalMomDeltaR_) &&
00309             (thisGlobalMomDeltaEta  < maxGlobalMomDeltaEta_) &&
00310             (thisGlobalMomDeltaPhi  < maxGlobalMomDeltaPhi_) &&
00311             (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
00312             float thisChi2         = useChi2_ ? getChi2(target, tsos, chi2DiagonalOnly_, chi2UseVertex_) : 0;
00313             if (thisChi2 >= maxChi2_) return false;
00314             switch (sortBy_) {
00315                 case LocalPosDiff:      isBest = (thisLocalPosDiff      < lastDeltaLocalPos); break;
00316                 case GlobalMomDeltaR:   isBest = (thisGlobalMomDeltaR   < lastDeltaR);        break;
00317                 case GlobalMomDeltaEta: isBest = (thisGlobalMomDeltaEta < lastDeltaEta);      break;
00318                 case GlobalMomDeltaPhi: isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi);      break;
00319                 case GlobalDPtRel:      isBest = (thisGlobalDPtRel      < lastGlobalDPtRel);  break;
00320                 case Chi2:              isBest = (thisChi2              < lastChi2);          break;
00321             }
00322             if (isBest) {
00323                 lastDeltaLocalPos = thisLocalPosDiff;
00324                 lastDeltaR        = thisGlobalMomDeltaR;
00325                 lastDeltaEta      = thisGlobalMomDeltaEta;
00326                 lastDeltaPhi      = thisGlobalMomDeltaPhi;
00327                 lastGlobalDPtRel  = thisGlobalDPtRel;
00328                 lastChi2          = thisChi2;
00329             } 
00330         }  // if match
00331     }
00332 
00333     return isBest;
00334 }
00335 
00336 bool
00337 MatcherUsingTracksAlgorithm::matchWithPropagation(const FreeTrajectoryState &start, 
00338                                               const FreeTrajectoryState &target, 
00339                                               float &lastDeltaR,
00340                                               float &lastDeltaEta,
00341                                               float &lastDeltaPhi, 
00342                                               float &lastDeltaLocalPos,
00343                                               float &lastGlobalDPtRel,
00344                                               float &lastChi2) const 
00345 {
00346     if ((start.momentum().mag() == 0) || (target.momentum().mag() == 0)) return false;
00347     TSCPBuilderNoMaterial propagator;
00348     /*2.2.X*/ try { 
00349     TrajectoryStateClosestToPoint tscp = propagator(start, target.position());
00350     // if (!tscp.isValid()) return false;  // in 3.1.X
00351     
00352     bool isBest = false;
00353     float thisLocalPosDiff = (tscp.position()-target.position()).mag();
00354     float thisGlobalMomDeltaR   = deltaR(tscp.momentum(), target.momentum());
00355     float thisGlobalMomDeltaPhi = fabs(deltaPhi(tscp.momentum().phi(), target.momentum().phi()));
00356     float thisGlobalMomDeltaEta = fabs(tscp.momentum().eta() - target.momentum().eta());
00357     float thisGlobalDPtRel = (tscp.momentum().perp() - target.momentum().perp())/target.momentum().perp();
00358 
00359     if ((thisLocalPosDiff       < maxLocalPosDiff_) &&
00360             (thisGlobalMomDeltaR    < maxGlobalMomDeltaR_) &&
00361             (thisGlobalMomDeltaEta  < maxGlobalMomDeltaEta_) &&
00362             (thisGlobalMomDeltaPhi  < maxGlobalMomDeltaPhi_) &&
00363             (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
00364         float thisChi2         = useChi2_ ? getChi2(target, tscp, chi2DiagonalOnly_, chi2UseVertex_) : 0;
00365         if (thisChi2 >= maxChi2_) return false;
00366         switch (sortBy_) {
00367             case LocalPosDiff:      isBest = (thisLocalPosDiff      < lastDeltaLocalPos); break;
00368             case GlobalMomDeltaR:   isBest = (thisGlobalMomDeltaR   < lastDeltaR);        break;
00369             case GlobalMomDeltaEta: isBest = (thisGlobalMomDeltaEta < lastDeltaEta);      break;
00370             case GlobalMomDeltaPhi: isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi);      break;
00371             case GlobalDPtRel:      isBest = (thisGlobalDPtRel      < lastGlobalDPtRel);  break;
00372             case Chi2:              isBest = (thisChi2              < lastChi2);          break;
00373         }
00374         if (isBest) {
00375             lastDeltaLocalPos = thisLocalPosDiff;
00376             lastDeltaR        = thisGlobalMomDeltaR;
00377             lastDeltaEta      = thisGlobalMomDeltaEta;
00378             lastDeltaPhi      = thisGlobalMomDeltaPhi;
00379             lastGlobalDPtRel  = thisGlobalDPtRel;
00380             lastChi2          = thisChi2;
00381         } 
00382     }  // if match
00383 
00384     return isBest;
00385     /*2.2.X*/ } catch (const TrajectoryStateException &err) { return false; } 
00386 }
00387 
00388 
00389 bool
00390 MatcherUsingTracksAlgorithm::matchByDirectComparison(const FreeTrajectoryState &start, 
00391                                                  const FreeTrajectoryState &target, 
00392                                                  float &lastDeltaR,
00393                                                  float &lastDeltaEta,
00394                                                  float &lastDeltaPhi, 
00395                                                  float &lastDeltaLocalPos,
00396                                                  float &lastGlobalDPtRel,
00397                                                  float &lastChi2) const 
00398 {
00399     if ((start.momentum().mag() == 0) || target.momentum().mag() == 0) return false;
00400 
00401     bool isBest = false;
00402     float thisLocalPosDiff = (start.position()-target.position()).mag();
00403     float thisGlobalMomDeltaR   = deltaR(start.momentum(), target.momentum());
00404     float thisGlobalMomDeltaPhi = fabs(deltaPhi(start.momentum().phi(), target.momentum().phi()));
00405     float thisGlobalMomDeltaEta = fabs(start.momentum().eta() - target.momentum().eta());
00406     float thisGlobalDPtRel = (start.momentum().perp() - target.momentum().perp())/target.momentum().perp();
00407 
00408     if ((thisLocalPosDiff       < maxLocalPosDiff_) &&
00409             (thisGlobalMomDeltaR    < maxGlobalMomDeltaR_) &&
00410             (thisGlobalMomDeltaEta  < maxGlobalMomDeltaEta_) &&
00411             (thisGlobalMomDeltaPhi  < maxGlobalMomDeltaPhi_) &&
00412             (fabs(thisGlobalDPtRel) < maxGlobalDPtRel_)) {
00413         float thisChi2 = useChi2_ ? getChi2(start, target, chi2DiagonalOnly_, chi2UseVertex_, chi2FirstMomentum_) : 0;
00414         if (thisChi2 >= maxChi2_) return false;
00415         switch (sortBy_) {
00416             case LocalPosDiff:      isBest = (thisLocalPosDiff      < lastDeltaLocalPos); break;
00417             case GlobalMomDeltaR:   isBest = (thisGlobalMomDeltaR   < lastDeltaR);        break;
00418             case GlobalMomDeltaEta: isBest = (thisGlobalMomDeltaEta < lastDeltaEta);      break;
00419             case GlobalMomDeltaPhi: isBest = (thisGlobalMomDeltaPhi < lastDeltaPhi);      break;
00420             case GlobalDPtRel:      isBest = (thisGlobalDPtRel      < lastGlobalDPtRel);  break;
00421             case Chi2:              isBest = (thisChi2              < lastChi2);          break;
00422         }
00423         if (isBest) {
00424             lastDeltaLocalPos = thisLocalPosDiff;
00425             lastDeltaR        = thisGlobalMomDeltaR;
00426             lastDeltaEta      = thisGlobalMomDeltaEta;
00427             lastDeltaPhi      = thisGlobalMomDeltaPhi;
00428             lastGlobalDPtRel  = thisGlobalDPtRel;
00429             lastChi2          = thisChi2;
00430         } 
00431     }  // if match
00432 
00433     return isBest;
00434 }
00435 
00436 double
00437 MatcherUsingTracksAlgorithm::getChi2(const FreeTrajectoryState &start, const FreeTrajectoryState &other, bool diagonalOnly, bool useVertex, bool useFirstMomentum) {
00438     if (!start.hasError() && !other.hasError()) throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
00439     AlgebraicSymMatrix55 cov;
00440     if (start.hasError()) cov += start.curvilinearError().matrix();
00441     if (other.hasError()) cov += other.curvilinearError().matrix();
00442     cropAndInvert(cov, diagonalOnly, !useVertex);
00443     GlobalVector p1 = start.momentum(), p2 = other.momentum();
00444     GlobalPoint  x1 = start.position(), x2 = other.position();
00445     GlobalVector p = useFirstMomentum ? p1 : p2; double pt = p.perp(), pm = p.mag();
00446     double dsz = (x1.z()-x2.z())*pt/pm - ( (x1.x()-x2.x())*p.x() + (x1.y()-x2.y())*p.y() ) / pt * p.z()/pm;
00447     double dxy = ( -(x1.x()-x2.x())*p.y() + (x1.y()-x2.y())*p.x() )/pt;
00448     AlgebraicVector5 diff(start.charge()/p1.mag()-other.charge()/p2.mag(),
00449                           p1.theta()-p2.theta(),
00450                           (p1.phi()-p2.phi()).value(), 
00451                           dxy,
00452                           dsz);
00453     return ROOT::Math::Similarity(diff, cov);
00454 }
00455 
00456 double
00457 MatcherUsingTracksAlgorithm::getChi2(const FreeTrajectoryState &start, const TrajectoryStateClosestToPoint &other, bool diagonalOnly, bool useVertex) {
00458     if (!start.hasError() && !other.hasError()) throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
00459     double pt; // needed by pgconvert
00460     AlgebraicSymMatrix55 cov;
00461     if (start.hasError()) cov += PerigeeConversions::ftsToPerigeeError(start).covarianceMatrix();
00462     if (other.hasError()) cov += other.perigeeError().covarianceMatrix();
00463     cropAndInvert(cov, diagonalOnly, !useVertex);
00464     AlgebraicVector5 pgpar1 = PerigeeConversions::ftsToPerigeeParameters(start,other.referencePoint(),pt).vector();
00465     AlgebraicVector5 pgpar2 = other.perigeeParameters().vector();
00466     AlgebraicVector5 diff(pgpar1-pgpar2);
00467     return ROOT::Math::Similarity(diff, cov);
00468 }
00469 
00470 double
00471 MatcherUsingTracksAlgorithm::getChi2(const TrajectoryStateOnSurface &start, const TrajectoryStateOnSurface &other, bool diagonalOnly, bool usePosition) {
00472     if (!start.hasError() && !other.hasError()) throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
00473     AlgebraicSymMatrix55 cov;
00474     if (start.hasError()) cov += start.localError().matrix();
00475     if (other.hasError()) cov += other.localError().matrix();
00476     cropAndInvert(cov, diagonalOnly, !usePosition);
00477     AlgebraicVector5 diff(start.localParameters().mixedFormatVector() - other.localParameters().mixedFormatVector());
00478     return ROOT::Math::Similarity(diff, cov);
00479 }
00480 
00481 void
00482 MatcherUsingTracksAlgorithm::cropAndInvert(AlgebraicSymMatrix55 &cov, bool diagonalOnly, bool top3by3only) {
00483     if (!top3by3only) {
00484         if (diagonalOnly) {
00485             for (size_t i = 0; i < 5; ++i) { for (size_t j = i+1; j < 5; ++j) {
00486                 cov(i,j) = 0;
00487             } }
00488         }
00489         cov.Invert();
00490     } else {
00491         // get 3x3 covariance
00492         AlgebraicSymMatrix33 momCov = cov.Sub<AlgebraicSymMatrix33>(0,0); // get 3x3 matrix
00493         if (diagonalOnly) { momCov(0,1) = 0; momCov(0,2) = 0; momCov(1,2) = 0; }
00494         // invert        
00495         momCov.Invert();
00496         // place it
00497         cov.Place_at(momCov,0,0);
00498         // zero the rest
00499         for (size_t i = 3; i < 5; ++i) { for (size_t j = i; j < 5; ++j) { 
00500             cov(i,j) = 0;
00501         } }
00502     }
00503 }
00504 
00505