SiStripFineDelayTLA.cc

Go to the documentation of this file.

#include <memory>
#include <string>
#include <iostream>
#include <TMath.h>
#include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTLA.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "DataFormats/GeometryVector/interface/LocalVector.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
#include <Geometry/CommonTopologies/interface/Topology.h>
#include <Geometry/CommonTopologies/interface/StripTopology.h>
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include <TrackingTools/PatternTools/interface/Trajectory.h>
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "DQM/SiStripCommissioningSources/plugins/tracking/SimpleTrackRefitter.h"

using namespace std;
SiStripFineDelayTLA::SiStripFineDelayTLA(edm::ParameterSet const& conf) : 
  conf_(conf)
{
  refitter_ = new SimpleTrackRefitter(conf);
}

void SiStripFineDelayTLA::init(const edm::Event& e, const edm::EventSetup& es)
{
  // get geometry
  edm::ESHandle<TrackerGeometry> estracker;
  es.get<TrackerDigiGeometryRecord>().get(estracker);
  tracker=&(*estracker);
  // the refitter
  refitter_->setServices(es);
}

// Virtual destructor needed.
SiStripFineDelayTLA::~SiStripFineDelayTLA() 
{  
  delete refitter_;
}  

std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const reco::Track& theT)
{
  std::vector<Trajectory> traj = refitter_->refitTrack(theT);
  return findtrackangle(traj);
}

std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const TrajectorySeed& seed, const reco::Track& theT){
  vector<Trajectory> traj = refitter_->refitTrack(seed,theT);
  return findtrackangle(traj);
}       

std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const std::vector<Trajectory>& trajVec)
{
  if (trajVec.size()) {
  return findtrackangle(trajVec.front()); }
  std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > hitangleassociation;
  return hitangleassociation;
}

std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> > SiStripFineDelayTLA::findtrackangle(const Trajectory& traj)
{
  std::vector<std::pair< std::pair<DetId, LocalPoint> ,float> >hitangleassociation;
  std::vector<TrajectoryMeasurement> TMeas=traj.measurements();
  std::vector<TrajectoryMeasurement>::iterator itm;
  for (itm=TMeas.begin();itm!=TMeas.end();itm++){
    TrajectoryStateOnSurface tsos=itm->updatedState();
    const TransientTrackingRecHit::ConstRecHitPointer thit=itm->recHit();
    const SiStripMatchedRecHit2D* matchedhit=dynamic_cast<const SiStripMatchedRecHit2D*>((*thit).hit());
    const SiStripRecHit2D* hit=dynamic_cast<const SiStripRecHit2D*>((*thit).hit());
    LocalVector trackdirection=tsos.localDirection();
    if(matchedhit){//if matched hit...
    GluedGeomDet * gdet=(GluedGeomDet *)tracker->idToDet(matchedhit->geographicalId());
    GlobalVector gtrkdir=gdet->toGlobal(trackdirection);
    // trackdirection on mono det
    // this the pointer to the mono hit of a matched hit 
    const SiStripRecHit2D *monohit=matchedhit->monoHit();
    const GeomDetUnit * monodet=gdet->monoDet();
    LocalVector monotkdir=monodet->toLocal(gtrkdir);
    if(monotkdir.z()!=0){
      // the local angle (mono)
          float localpitch = ((StripTopology*)(&monodet->topology()))->localPitch(tsos.localPosition());
          float thickness = ((((((monohit->geographicalId())>>25)&0x7f)==0xd)||
                         ((((monohit->geographicalId())>>25)&0x7f)==0xe))&&
                       ((((monohit->geographicalId())>>5)&0x7)>4)) ? 0.0500 : 0.0320;
          float angle = computeAngleCorr(monotkdir, localpitch, thickness);
      hitangleassociation.push_back(make_pair(make_pair(monohit->geographicalId(),monohit->localPosition()), angle)); 
      // trackdirection on stereo det
      // this the pointer to the stereo hit of a matched hit 
      const SiStripRecHit2D *stereohit=matchedhit->stereoHit();
      const GeomDetUnit * stereodet=gdet->stereoDet(); 
      LocalVector stereotkdir=stereodet->toLocal(gtrkdir);
      if(stereotkdir.z()!=0){
        // the local angle (stereo)
            float localpitch = ((StripTopology*)(&stereodet->topology()))->localPitch(tsos.localPosition());
            float thickness = ((((((stereohit->geographicalId())>>25)&0x7f)==0xd)||
                           ((((stereohit->geographicalId())>>25)&0x7f)==0xe))&&
                         ((((stereohit->geographicalId())>>5)&0x7)>4)) ? 0.0500 : 0.0320;
            float angle = computeAngleCorr(stereotkdir, localpitch, thickness);
        hitangleassociation.push_back(make_pair(make_pair(stereohit->geographicalId(),stereohit->localPosition()), angle)); 
      }
    }
    }
    else if(hit){
        GeomDetUnit * det=(GeomDetUnit *)tracker->idToDet(hit->geographicalId());
    //  hit= pointer to the rechit
    if(trackdirection.z()!=0){
          // the local angle (single hit)
          float localpitch = ((StripTopology*)(&det->topology()))->localPitch(tsos.localPosition());
          float thickness = ((((((hit->geographicalId())>>25)&0x7f)==0xd)||
                         ((((hit->geographicalId())>>25)&0x7f)==0xe))&&
                       ((((hit->geographicalId())>>5)&0x7)>4)) ? 0.0500 : 0.0320;
          float angle = computeAngleCorr(trackdirection, localpitch, thickness);
      hitangleassociation.push_back(make_pair(make_pair(hit->geographicalId(),hit->localPosition()), angle)); 
    }
    }
  }
  return hitangleassociation;
}

double SiStripFineDelayTLA::computeAngleCorr(const LocalVector& v, double pitch, double thickness)
{
  double v_xy = sqrt(v.x()*v.x()+v.y()*v.y());
  double L = fabs(thickness*v_xy/v.z());
  double Lmax = fabs(pitch/v.x()*v_xy);
  if(L<Lmax) {
    LogDebug("SiStripFineDelayTLA ") << L << " vs " << Lmax 
        << " Signal contained in strip. Correction is " << v.z()/v.mag();
    return v.z()/v.mag();
  } else {
    LogDebug("SiStripFineDelayTLA ") << L << " vs " << Lmax
       << " Signal not contained in strip. Correction is " << thickness/pitch*v.x()/v_xy*v.z()/v.mag()
       << " instead of " << v.z()/v.mag();
    return thickness/pitch*v.x()/v_xy*v.z()/v.mag();
  }
}