Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:46:21

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