TrackInfoProducer/src/TrackInfoProducerAlgorithm.cc

0001 #include "AnalysisAlgos/TrackInfoProducer/interface/TrackInfoProducerAlgorithm.h"
0002 #include "AnalysisDataFormats/TrackInfo/interface/TrackingRecHitInfo.h"
0003 #include "DataFormats/TrackReco/interface/Track.h"
0004 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0005 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0006 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0007
0008 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0010 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0011 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0012 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0013
0014 using namespace reco;
0015
0016 void TrackInfoProducerAlgorithm::run(const edm::Ref<std::vector<Trajectory> > traj_iterator,
0017                                      TrackRef track,
0018                                      TrackInfo& output,
0019                                      const TrackerGeometry* tracker) {
0020   std::vector<TrajectoryMeasurement> measurements = traj_iterator->measurements();
0021
0022   std::vector<TrajectoryMeasurement>::iterator traj_mes_iterator;
0023   //edm::LogInfo("TrackInfoProducer") << "Number of Measurements: "<<measurements.size();
0024   TrackInfo::TrajectoryInfo trajinfo;
0025   for (traj_mes_iterator = measurements.begin(); traj_mes_iterator != measurements.end();
0026        ++traj_mes_iterator) {  //loop on measurements
0027
0028     TrajectoryStateOnSurface fwdtsos = traj_mes_iterator->forwardPredictedState();
0029     TrajectoryStateOnSurface bwdtsos = traj_mes_iterator->backwardPredictedState();
0030     TrajectoryStateOnSurface updatedtsos = traj_mes_iterator->updatedState();
0031     TrajectoryStateCombiner statecombiner;
0032     TrajectoryStateOnSurface combinedtsos = statecombiner.combine(fwdtsos, bwdtsos);
0033
0034     ConstRecHitPointer ttrh = traj_mes_iterator->recHit();
0035     LocalPoint pos;
0036     if (ttrh->isValid())
0037       pos = ttrh->hit()->localPosition();
0038     unsigned int detid = ttrh->hit()->geographicalId().rawId();
0039
0040     TrackingRecHitRef thehitref;
0041     TrackingRecHit const* thehitptr = nullptr;
0042     int i = 0, j = 0;
0043
0044     for (auto const& thehit : track->recHits()) {
0045       i++;
0046       LocalPoint hitpos;
0047       if (thehit->isValid())
0048         hitpos = thehit->localPosition();
0049       if (thehit->geographicalId().rawId() == detid && (hitpos - pos).mag() < 1e-4) {
0050         thehitptr = thehit;
0051         thehitref = track->extra()->recHitRef(i - 1);
0052         j++;
0053         break;
0054       }
0055     }
0056
0057     PTrajectoryStateOnDet const& fwdptsod = trajectoryStateTransform::persistentState(fwdtsos, detid);
0058     PTrajectoryStateOnDet const& bwdptsod = trajectoryStateTransform::persistentState(bwdtsos, detid);
0059     PTrajectoryStateOnDet const& updatedptsod = trajectoryStateTransform::persistentState(updatedtsos, detid);
0060     PTrajectoryStateOnDet const& combinedptsod = trajectoryStateTransform::persistentState(combinedtsos, detid);
0061
0062     const ProjectedSiStripRecHit2D* phit = dynamic_cast<const ProjectedSiStripRecHit2D*>(thehitptr);
0063     const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(thehitptr);
0064
0065     RecHitType type = Single;
0066     LocalVector monofwd, stereofwd;
0067     LocalVector monobwd, stereobwd;
0068     LocalVector monoco, stereoco;
0069     LocalVector monoup, stereoup;
0070
0071     LocalPoint pmonofwd, pstereofwd;
0072     LocalPoint pmonobwd, pstereobwd;
0073     LocalPoint pmonoco, pstereoco;
0074     LocalPoint pmonoup, pstereoup;
0075     if (matchedhit) {
0076       type = Matched;
0077       const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(tracker->idToDet(matchedhit->geographicalId()));
0078
0079       GlobalVector gtrkdirfwd = gdet->toGlobal(fwdptsod.parameters().momentum());
0080       GlobalVector gtrkdirbwd = gdet->toGlobal(bwdptsod.parameters().momentum());
0081       GlobalVector gtrkdirup = gdet->toGlobal(updatedptsod.parameters().momentum());
0082       GlobalVector gtrkdirco = gdet->toGlobal(combinedptsod.parameters().momentum());
0083
0084       const GeomDetUnit* monodet = gdet->monoDet();
0085
0086       monofwd = monodet->toLocal(gtrkdirfwd);
0087       monobwd = monodet->toLocal(gtrkdirbwd);
0088       monoup = monodet->toLocal(gtrkdirup);
0089       monoco = monodet->toLocal(gtrkdirco);
0090
0091       pmonofwd = project(gdet, monodet, fwdptsod.parameters().position(), monofwd);
0092       pmonobwd = project(gdet, monodet, bwdptsod.parameters().position(), monobwd);
0093       pmonoup = project(gdet, monodet, updatedptsod.parameters().position(), monoup);
0094       pmonoco = project(gdet, monodet, combinedptsod.parameters().position(), monoco);
0095
0096       const GeomDetUnit* stereodet = gdet->stereoDet();
0097
0098       stereofwd = stereodet->toLocal(gtrkdirfwd);
0099       stereobwd = stereodet->toLocal(gtrkdirbwd);
0100       stereoup = stereodet->toLocal(gtrkdirup);
0101       stereoco = stereodet->toLocal(gtrkdirco);
0102
0103       pstereofwd = project(gdet, stereodet, fwdptsod.parameters().position(), stereofwd);
0104       pstereobwd = project(gdet, stereodet, bwdptsod.parameters().position(), stereobwd);
0105       pstereoup = project(gdet, stereodet, updatedptsod.parameters().position(), stereoup);
0106       pstereoco = project(gdet, stereodet, combinedptsod.parameters().position(), stereoco);
0107
0108     } else if (phit) {
0109       type = Projected;
0110       const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(tracker->idToDet(phit->geographicalId()));
0111
0112       GlobalVector gtrkdirfwd = gdet->toGlobal(fwdptsod.parameters().momentum());
0113       GlobalVector gtrkdirbwd = gdet->toGlobal(bwdptsod.parameters().momentum());
0114       GlobalVector gtrkdirup = gdet->toGlobal(updatedptsod.parameters().momentum());
0115       GlobalVector gtrkdirco = gdet->toGlobal(combinedptsod.parameters().momentum());
0116       const SiStripRecHit2D& originalhit = phit->originalHit();
0117       const GeomDetUnit* det;
0118       if (!StripSubdetector(originalhit.geographicalId().rawId()).stereo()) {
0119         det = gdet->monoDet();
0120         monofwd = det->toLocal(gtrkdirfwd);
0121         monobwd = det->toLocal(gtrkdirbwd);
0122         monoup = det->toLocal(gtrkdirup);
0123         monoco = det->toLocal(gtrkdirco);
0124         pmonofwd = project(gdet, det, fwdptsod.parameters().position(), monofwd);
0125         pmonobwd = project(gdet, det, bwdptsod.parameters().position(), monobwd);
0126         pmonoup = project(gdet, det, updatedptsod.parameters().position(), monoup);
0127         pmonoco = project(gdet, det, combinedptsod.parameters().position(), monoco);
0128       } else {
0129         det = gdet->stereoDet();
0130         stereofwd = det->toLocal(gtrkdirfwd);
0131         stereobwd = det->toLocal(gtrkdirbwd);
0132         stereoup = det->toLocal(gtrkdirup);
0133         stereoco = det->toLocal(gtrkdirco);
0134         pstereofwd = project(gdet, det, fwdptsod.parameters().position(), stereofwd);
0135         pstereobwd = project(gdet, det, bwdptsod.parameters().position(), stereobwd);
0136         pstereoup = project(gdet, det, updatedptsod.parameters().position(), stereoup);
0137         pstereoco = project(gdet, det, combinedptsod.parameters().position(), stereoco);
0138       }
0139     }
0140     TrackingRecHitInfo::TrackingStates states;
0141     if (!forwardPredictedStateTag_.empty())
0142       states.insert(std::make_pair(
0143           FwPredicted,
0144           TrackingStateInfo(std::make_pair(monofwd, stereofwd), std::make_pair(pmonofwd, pstereofwd), fwdptsod)));
0145     if (!backwardPredictedStateTag_.empty())
0146       states.insert(std::make_pair(
0147           BwPredicted,
0148           TrackingStateInfo(std::make_pair(monobwd, stereobwd), std::make_pair(pmonobwd, pstereobwd), bwdptsod)));
0149     if (!updatedStateTag_.empty())
0150       states.insert(std::make_pair(
0151           Updated,
0152           TrackingStateInfo(std::make_pair(monoup, stereoup), std::make_pair(pmonoup, pstereoup), updatedptsod)));
0153     if (!combinedStateTag_.empty())
0154       states.insert(std::make_pair(
0155           Combined,
0156           TrackingStateInfo(std::make_pair(monoco, stereoco), std::make_pair(pmonoco, pstereoco), combinedptsod)));
0157
0158     TrackingRecHitInfo tkRecHitInfo(type, states);
0159
0160     if (j != 0) {
0161       trajinfo.insert(std::make_pair(thehitref, tkRecHitInfo));
0162     }
0163     //      else  edm::LogInfo("TrackInfoProducer") << "RecHit not associated ";
0164   }
0165   output = TrackInfo((traj_iterator->seed()), trajinfo);
0166 }
0167
0168 LocalPoint TrackInfoProducerAlgorithm::project(const GeomDet* det,
0169                                                const GeomDet* projdet,
0170                                                LocalPoint position,
0171                                                LocalVector trackdirection) const {
0172   GlobalPoint globalpoint = (det->surface()).toGlobal(position);
0173
0174   // position of the initial and final point of the strip in glued local coordinates
0175   LocalPoint projposition = (projdet->surface()).toLocal(globalpoint);
0176
0177   //correct the position with the track direction
0178
0179   float scale = -projposition.z() / trackdirection.z();
0180
0181   projposition += scale * trackdirection;
0182
0183   return projposition;
0184 }