00001 #include <memory>
00002 #include <string>
00003 #include <iostream>
00004 #include <TMath.h>
00005 #include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTLA.h"
00006
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 #include "FWCore/Framework/interface/ESHandle.h"
00009 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00010 #include "DataFormats/GeometryVector/interface/LocalVector.h"
00011 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00012 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00013 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00014 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
00015 #include <Geometry/CommonTopologies/interface/Topology.h>
00016 #include <Geometry/CommonTopologies/interface/StripTopology.h>
00017 #include "DataFormats/TrackReco/interface/Track.h"
00018 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00019 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00020 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
00021 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00022 #include <TrackingTools/PatternTools/interface/Trajectory.h>
00023 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00024 #include "DQM/SiStripCommissioningSources/plugins/tracking/SimpleTrackRefitter.h"
00025
00026 using namespace std;
00027 SiStripFineDelayTLA::SiStripFineDelayTLA(edm::ParameterSet const& conf) :
00028 conf_(conf)
00029 {
00030 refitter_ = new SimpleTrackRefitter(conf);
00031 }
00032
00033 void SiStripFineDelayTLA::init(const edm::Event& e, const edm::EventSetup& es)
00034 {
00035
00036 edm::ESHandle<TrackerGeometry> estracker;
00037 es.get<TrackerDigiGeometryRecord>().get(estracker);
00038 tracker=&(*estracker);
00039
00040 refitter_->setServices(es);
00041 }
00042
00043
00044 SiStripFineDelayTLA::~SiStripFineDelayTLA()
00045 {
00046 delete refitter_;
00047 }
00048
00049 std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const reco::Track& theT)
00050 {
00051 std::vector<Trajectory> traj = refitter_->refitTrack(theT);
00052 return findtrackangle(traj);
00053 }
00054
00055 std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const TrajectorySeed& seed, const reco::Track& theT){
00056 vector<Trajectory> traj = refitter_->refitTrack(seed,theT);
00057 return findtrackangle(traj);
00058 }
00059
00060 std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const std::vector<Trajectory>& trajVec)
00061 {
00062 if (trajVec.size()) {
00063 return findtrackangle(trajVec.front()); }
00064 std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > hitangleassociation;
00065 return hitangleassociation;
00066 }
00067
00068 std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const Trajectory& traj)
00069 {
00070 std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> >hitangleassociation;
00071 std::vector<TrajectoryMeasurement> TMeas=traj.measurements();
00072 std::vector<TrajectoryMeasurement>::iterator itm;
00073 for (itm=TMeas.begin();itm!=TMeas.end();itm++){
00074 TrajectoryStateOnSurface tsos=itm->updatedState();
00075 const TransientTrackingRecHit::ConstRecHitPointer thit=itm->recHit();
00076 const SiStripMatchedRecHit2D* matchedhit=dynamic_cast<const SiStripMatchedRecHit2D*>((*thit).hit());
00077 const SiStripRecHit2D* hit=dynamic_cast<const SiStripRecHit2D*>((*thit).hit());
00078 LocalVector trackdirection=tsos.localDirection();
00079 if(matchedhit){
00080 GluedGeomDet * gdet=(GluedGeomDet *)tracker->idToDet(matchedhit->geographicalId());
00081 GlobalVector gtrkdir=gdet->toGlobal(trackdirection);
00082
00083
00084 const SiStripRecHit2D *monohit=matchedhit->monoHit();
00085 const GeomDetUnit * monodet=gdet->monoDet();
00086 LocalVector monotkdir=monodet->toLocal(gtrkdir);
00087 if(monotkdir.z()!=0){
00088
00089 float localpitch = ((StripTopology*)(&monodet->topology()))->localPitch(tsos.localPosition());
00090 float thickness = ((((((monohit->geographicalId())>>25)&0x7f)==0xd)||
00091 ((((monohit->geographicalId())>>25)&0x7f)==0xe))&&
00092 ((((monohit->geographicalId())>>5)&0x7)>4)) ? 0.0500 : 0.0320;
00093 float angle = computeAngleCorr(monotkdir, localpitch, thickness);
00094 hitangleassociation.push_back(make_pair(make_pair(monohit->geographicalId(),monohit->localPosition()), angle));
00095
00096
00097 const SiStripRecHit2D *stereohit=matchedhit->stereoHit();
00098 const GeomDetUnit * stereodet=gdet->stereoDet();
00099 LocalVector stereotkdir=stereodet->toLocal(gtrkdir);
00100 if(stereotkdir.z()!=0){
00101
00102 float localpitch = ((StripTopology*)(&stereodet->topology()))->localPitch(tsos.localPosition());
00103 float thickness = ((((((stereohit->geographicalId())>>25)&0x7f)==0xd)||
00104 ((((stereohit->geographicalId())>>25)&0x7f)==0xe))&&
00105 ((((stereohit->geographicalId())>>5)&0x7)>4)) ? 0.0500 : 0.0320;
00106 float angle = computeAngleCorr(stereotkdir, localpitch, thickness);
00107 hitangleassociation.push_back(make_pair(make_pair(stereohit->geographicalId(),stereohit->localPosition()), angle));
00108 }
00109 }
00110 }
00111 else if(hit){
00112 GeomDetUnit * det=(GeomDetUnit *)tracker->idToDet(hit->geographicalId());
00113
00114 if(trackdirection.z()!=0){
00115
00116 float localpitch = ((StripTopology*)(&det->topology()))->localPitch(tsos.localPosition());
00117 float thickness = ((((((hit->geographicalId())>>25)&0x7f)==0xd)||
00118 ((((hit->geographicalId())>>25)&0x7f)==0xe))&&
00119 ((((hit->geographicalId())>>5)&0x7)>4)) ? 0.0500 : 0.0320;
00120 float angle = computeAngleCorr(trackdirection, localpitch, thickness);
00121 hitangleassociation.push_back(make_pair(make_pair(hit->geographicalId(),hit->localPosition()), angle));
00122 }
00123 }
00124 }
00125 return hitangleassociation;
00126 }
00127
00128 double SiStripFineDelayTLA::computeAngleCorr(const LocalVector& v, double pitch, double thickness)
00129 {
00130 double v_xy = sqrt(v.x()*v.x()+v.y()*v.y());
00131 double L = fabs(thickness*v_xy/v.z());
00132 double Lmax = fabs(pitch/v.x()*v_xy);
00133 if(L<Lmax) {
00134 LogDebug("SiStripFineDelayTLA ") << L << " vs " << Lmax
00135 << " Signal contained in strip. Correction is " << v.z()/v.mag();
00136 return v.z()/v.mag();
00137 } else {
00138 LogDebug("SiStripFineDelayTLA ") << L << " vs " << Lmax
00139 << " Signal not contained in strip. Correction is " << thickness/pitch*v.x()/v_xy*v.z()/v.mag()
00140 << " instead of " << v.z()/v.mag();
00141 return thickness/pitch*v.x()/v_xy*v.z()/v.mag();
00142 }
00143 }
00144