SiStripFineDelayHit.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    SiStripFineDelayHit
// Class:      SiStripFineDelayHit
//
//
// Original Author:  Christophe DELAERE
//         Created:  Fri Nov 17 10:52:42 CET 2006
//
//
 
// system include files
#include <memory>
#include <utility>
#include <vector>
#include <algorithm>
 
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
 
#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include <DataFormats/SiStripCommon/interface/ConstantsForRunType.h>
#include <DataFormats/SiStripCommon/interface/SiStripFedKey.h>
#include <CondFormats/SiStripObjects/interface/FedChannelConnection.h>
#include <CondFormats/SiStripObjects/interface/SiStripFedCabling.h>
#include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
#include <CondFormats/SiStripObjects/interface/SiStripNoises.h>
#include <CondFormats/DataRecord/interface/SiStripFedCablingRcd.h>
#include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
#include <CondFormats/SiStripObjects/interface/SiStripNoises.h>
 
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "DataFormats/GeometryVector/interface/LocalVector.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include <Geometry/CommonTopologies/interface/Topology.h>
#include <Geometry/CommonTopologies/interface/StripTopology.h>
 
#include <TrackingTools/PatternTools/interface/Trajectory.h>
 
#include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayHit.h"
#include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTLA.h"
#include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTOF.h"
 
#include "TMath.h"
 
//
// constructors and destructor
//
SiStripFineDelayHit::SiStripFineDelayHit(const edm::ParameterSet& iConfig) : event_(nullptr) {
  //register your products
  produces<edm::DetSetVector<SiStripRawDigi> >("FineDelaySelection");
  //now do what ever other initialization is needed
  anglefinder_ = new SiStripFineDelayTLA(iConfig);
  cosmic_ = iConfig.getParameter<bool>("cosmic");
  field_ = iConfig.getParameter<bool>("MagneticField");
  maxAngle_ = iConfig.getParameter<double>("MaxTrackAngle");
  minTrackP2_ = iConfig.getParameter<double>("MinTrackMomentum") * iConfig.getParameter<double>("MinTrackMomentum");
  maxClusterDistance_ = iConfig.getParameter<double>("MaxClusterDistance");
  /*
   clusterLabel_ = iConfig.getParameter<edm::InputTag>("ClustersLabel");
   trackLabel_ = iConfig.getParameter<edm::InputTag>("TracksLabel");
   seedLabel_  = iConfig.getParameter<edm::InputTag>("SeedsLabel");
   inputModuleLabel_ = iConfig.getParameter<edm::InputTag>( "InputModuleLabel" ) ;
   digiLabel_ = iConfig.getParameter<edm::InputTag>("DigiLabel");
   */
  clustersToken_ =
      consumes<edmNew::DetSetVector<SiStripCluster> >(iConfig.getParameter<edm::InputTag>("ClustersLabel"));
  trackToken_ = consumes<std::vector<Trajectory> >(iConfig.getParameter<edm::InputTag>("TracksLabel"));
  trackCollectionToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("TracksLabel"));
  seedcollToken_ = consumes<TrajectorySeedCollection>(iConfig.getParameter<edm::InputTag>("SeedsLabel"));
  inputModuleToken_ = consumes<SiStripEventSummary>(iConfig.getParameter<edm::InputTag>("InputModuleLabel"));
  digiToken_ = consumes<edm::DetSetVector<SiStripDigi> >(iConfig.getParameter<edm::InputTag>("DigiLabel"));
 
  homeMadeClusters_ = iConfig.getParameter<bool>("NoClustering");
  explorationWindow_ = iConfig.getParameter<uint32_t>("ExplorationWindow");
  noTracking_ = iConfig.getParameter<bool>("NoTracking");
  mode_ = 0;
}
 
SiStripFineDelayHit::~SiStripFineDelayHit() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
  delete anglefinder_;
}
 
//
// member functions
//
SiStripFineDelayHit::DeviceMask SiStripFineDelayHit::deviceMask(const StripSubdetector::SubDetector subdet,
                                                                const int substructure,
                                                                const TrackerTopology* tkrTopo) {
  uint32_t rootDetId = 0;
  uint32_t maskDetId = 0;
 
  switch (subdet) {
    case StripSubdetector::TIB: {
      rootDetId = tkrTopo->tibDetId(substructure, 0, 0, 0, 0, 0).rawId();
      maskDetId = tkrTopo->tibDetId(15, 0, 0, 0, 0, 0).rawId();
      break;
    }
    case StripSubdetector::TID: {
      rootDetId = tkrTopo->tidDetId(substructure > 0 ? 2 : 1, abs(substructure), 0, 0, 0, 0).rawId();
      maskDetId = tkrTopo->tidDetId(3, 15, 0, 0, 0, 0).rawId();
      break;
    }
    case StripSubdetector::TOB: {
      rootDetId = tkrTopo->tobDetId(substructure, 0, 0, 0, 0).rawId();
      maskDetId = tkrTopo->tobDetId(15, 0, 0, 0, 0).rawId();
      break;
    }
    case StripSubdetector::TEC: {
      rootDetId = tkrTopo->tecDetId(substructure > 0 ? 2 : 1, abs(substructure), 0, 0, 0, 0, 0).rawId();
      maskDetId = tkrTopo->tecDetId(3, 15, 0, 0, 0, 0, 0).rawId();
      break;
    }
    default:
      break;
  }
  return std::make_pair(maskDetId, rootDetId);
}
 
std::vector<std::pair<uint32_t, std::pair<double, double> > > SiStripFineDelayHit::detId(
    const TrackerGeometry& tracker,
    const TrackerTopology* tkrTopo,
    const reco::Track* tk,
    const std::vector<Trajectory>& trajVec,
    const StripSubdetector::SubDetector subdet,
    const int substructure) {
  if (substructure == 0xff)
    return detId(tracker, tkrTopo, tk, trajVec, 0, 0);
  // first determine the root detId we are looking for
  DeviceMask mask = deviceMask(subdet, substructure, tkrTopo);
  // then call the method that loops on recHits
  return detId(tracker, tkrTopo, tk, trajVec, mask.first, mask.second);
}
 
std::vector<std::pair<uint32_t, std::pair<double, double> > > SiStripFineDelayHit::detId(
    const TrackerGeometry& tracker,
    const TrackerTopology* tkrTopo,
    const reco::Track* tk,
    const std::vector<Trajectory>& trajVec,
    const uint32_t& maskDetId,
    const uint32_t& rootDetId) {
  bool onDisk = ((maskDetId == tkrTopo->tidDetId(3, 15, 0, 0, 0, 0).rawId()) ||
                 (maskDetId == tkrTopo->tecDetId(3, 15, 0, 0, 0, 0, 0).rawId()));
  std::vector<std::pair<uint32_t, std::pair<double, double> > > result;
  std::vector<uint32_t> usedDetids;
  // now loop on recHits to find the right detId plus the track local angle
  std::vector<std::pair<std::pair<DetId, LocalPoint>, float> > hitangle;
  if (!cosmic_) {
    // use trajectories in event.
    // we have first to find the right trajectory for the considered track.
    for (std::vector<Trajectory>::const_iterator traj = trajVec.begin(); traj < trajVec.end(); ++traj) {
      if (((traj->lastMeasurement().recHit()->geographicalId().rawId() ==
            (*(tk->recHitsEnd() - 1))->geographicalId().rawId()) &&
           (traj->lastMeasurement().recHit()->localPosition().x() == (*(tk->recHitsEnd() - 1))->localPosition().x())) ||
          ((traj->firstMeasurement().recHit()->geographicalId().rawId() ==
            (*(tk->recHitsEnd() - 1))->geographicalId().rawId()) &&
           (traj->firstMeasurement().recHit()->localPosition().x() == (*(tk->recHitsEnd() - 1))->localPosition().x()))) {
        hitangle = anglefinder_->findtrackangle(*traj);
        break;
      }
    }
  } else {
    edm::Handle<TrajectorySeedCollection> seedcoll;
    //    event_->getByLabel(seedLabel_,seedcoll);
    event_->getByToken(seedcollToken_, seedcoll);
    // use trajectories in event.
    hitangle = anglefinder_->findtrackangle(trajVec);
  }
  LogDebug("DetId") << "number of hits for the track: " << hitangle.size();
  std::vector<std::pair<std::pair<DetId, LocalPoint>, float> >::iterator iter;
  // select the interesting DetIds, based on the ID and TLA
  for (iter = hitangle.begin(); iter != hitangle.end(); iter++) {
    // check the detId.
    // if substructure was 0xff, then maskDetId and rootDetId == 0
    // this implies all detids are accepted. (also if maskDetId=rootDetId=0 explicitely).
    // That "unusual" mode of operation allows to analyze also Latency scans
    LogDebug("DetId") << "check the detid: " << std::hex << (iter->first.first.rawId()) << " vs " << rootDetId
                      << " with a mask of " << maskDetId << std::dec << std::endl;
 
    if (((iter->first.first.rawId() & maskDetId) != rootDetId))
      continue;
    if (std::find(usedDetids.begin(), usedDetids.end(), iter->first.first.rawId()) != usedDetids.end())
      continue;
    // check the local angle (extended to the equivalent angle correction)
    LogDebug("DetId") << "check the angle: " << fabs((iter->second));
    if (1 - fabs(fabs(iter->second) - 1) < cos(maxAngle_ / 180. * TMath::Pi()))
      continue;
    // returns the detid + the time of flight to there
    std::pair<uint32_t, std::pair<double, double> > el;
    std::pair<double, double> subel;
    el.first = iter->first.first.rawId();
    // here, we compute the TOF.
    // For cosmics, some track parameters are missing. Parameters are recomputed.
    // for our calculation, the track momemtum at any point is enough:
    // only used without B field or for the sign of Pz.
    double trackParameters[5];
    for (int i = 0; i < 5; i++)
      trackParameters[i] = tk->parameters()[i];
    if (cosmic_)
      SiStripFineDelayTOF::trackParameters(*tk, trackParameters);
    double hit[3];
    const GeomDetUnit* det(tracker.idToDetUnit(iter->first.first));
    Surface::GlobalPoint gp = det->surface().toGlobal(iter->first.second);
    hit[0] = gp.x();
    hit[1] = gp.y();
    hit[2] = gp.z();
    double phit[3];
    phit[0] = tk->momentum().x();
    phit[1] = tk->momentum().y();
    phit[2] = tk->momentum().z();
    subel.first = SiStripFineDelayTOF::timeOfFlight(cosmic_, field_, trackParameters, hit, phit, onDisk);
    subel.second = iter->second;
    el.second = subel;
    // returns the detid + TOF
    result.push_back(el);
    usedDetids.push_back(el.first);
  }
  return result;
}
 
bool SiStripFineDelayHit::rechit(reco::Track* tk, uint32_t det_id) {
  for (trackingRecHit_iterator it = tk->recHitsBegin(); it != tk->recHitsEnd(); it++)
    if ((*it)->geographicalId().rawId() == det_id) {
      return (*it)->isValid();
      break;
    }
  return false;
}
 
// VI January 2012: FIXME
// do not understand what is going on here: each hit has a cluster: by definition will be the closest!
std::pair<const SiStripCluster*, double> SiStripFineDelayHit::closestCluster(
    const TrackerGeometry& tracker,
    const reco::Track* tk,
    const uint32_t& det_id,
    const edmNew::DetSetVector<SiStripCluster>& clusters,
    const edm::DetSetVector<SiStripDigi>& hits) {
  std::pair<const SiStripCluster*, double> result(nullptr, 0.);
  double hitStrip = -1;
  int nstrips = -1;
  // localize the crossing point of the track on the module
  for (trackingRecHit_iterator it = tk->recHitsBegin(); it != tk->recHitsEnd(); it++) {
    LogDebug("closestCluster") << "(*it)->geographicalId().rawId() vs det_id" << (*it)->geographicalId().rawId() << " "
                               << det_id;
    //handle the mono rechits
    if ((*it)->geographicalId().rawId() == det_id) {
      if (!(*it)->isValid())
        continue;
      LogDebug("closestCluster") << " using the single mono hit";
      LocalPoint lp = (*it)->localPosition();
      const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker.idToDet((*it)->geographicalId()));
      MeasurementPoint p = gdu->topology().measurementPosition(lp);
      hitStrip = p.x();
      nstrips = (dynamic_cast<const StripTopology*>(&(gdu->topology())))->nstrips();
      break;
    }
    /* FIXME: local position is not there anymore...
    //handle stereo part of matched hits
    //one could try to cast to SiStripMatchedRecHit2D but it is faster to look at the detid
    else if((det_id - (*it)->geographicalId().rawId())==1) {
      const SiStripMatchedRecHit2D* hit2D = dynamic_cast<const SiStripMatchedRecHit2D*>(&(**it));
      if(!hit2D) continue; // this is a security that should never trigger
      const SiStripRecHit2D* stereo = hit2D->stereoHit();
      if(!stereo) continue; // this is a security that should never trigger
      if(!stereo->isValid()) continue;
      LogDebug("closestCluster") << " using the stereo hit";
      LocalPoint lp = stereo->localPosition();
      const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker.idToDet(stereo->geographicalId()));
      MeasurementPoint p = gdu->topology().measurementPosition(lp);
      hitStrip = p.x();
      nstrips = (dynamic_cast<const StripTopology*>(&(gdu->topology())))->nstrips();
      break;
    }
    //handle mono part of matched hits
    //one could try to cast to SiStripMatchedRecHit2D but it is faster to look at the detid
    else if((det_id - (*it)->geographicalId().rawId())==2) {
      const SiStripMatchedRecHit2D* hit2D = dynamic_cast<const SiStripMatchedRecHit2D*>(&(**it));
      if(!hit2D) continue; // this is a security that should never trigger
      const SiStripRecHit2D* mono = hit2D->monoHit();
      if(!mono) continue; // this is a security that should never trigger
      if(!mono->isValid()) continue;
      LogDebug("closestCluster") << " using the mono hit";
      LocalPoint lp = mono->localPosition();
      const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker.idToDet(mono->geographicalId()));
      MeasurementPoint p = gdu->topology().measurementPosition(lp);
      hitStrip = p.x();
      nstrips = (dynamic_cast<const StripTopology*>(&(gdu->topology())))->nstrips();
      break;
    }
    */
  }
  LogDebug("closestCluster") << " hit strip = " << hitStrip;
  if (hitStrip < 0)
    return result;
  if (homeMadeClusters_) {
    // take the list of digis on the module
    for (edm::DetSetVector<SiStripDigi>::const_iterator DSViter = hits.begin(); DSViter != hits.end(); DSViter++) {
      if (DSViter->id == det_id) {
        // loop from hitstrip-n to hitstrip+n (explorationWindow_) and select the highest strip
        int minStrip = int(round(hitStrip)) - explorationWindow_;
        minStrip = minStrip < 0 ? 0 : minStrip;
        int maxStrip = int(round(hitStrip)) + explorationWindow_ + 1;
        maxStrip = maxStrip >= nstrips ? nstrips - 1 : maxStrip;
        edm::DetSet<SiStripDigi>::const_iterator rangeStart = DSViter->end();
        edm::DetSet<SiStripDigi>::const_iterator rangeStop = DSViter->end();
        for (edm::DetSet<SiStripDigi>::const_iterator digiIt = DSViter->begin(); digiIt != DSViter->end(); ++digiIt) {
          if (digiIt->strip() >= minStrip && rangeStart == DSViter->end())
            rangeStart = digiIt;
          if (digiIt->strip() <= maxStrip)
            rangeStop = digiIt;
        }
        if (rangeStart != DSViter->end()) {
          if (rangeStop != DSViter->end())
            ++rangeStop;
          // build a fake cluster
          LogDebug("closestCluster") << "build a fake cluster ";
          SiStripCluster* newCluster =
              new SiStripCluster(SiStripCluster::SiStripDigiRange(rangeStart, rangeStop));  // /!\ ownership transfered
          result.first = newCluster;
          result.second = fabs(newCluster->barycenter() - hitStrip);
        }
        break;
      }
    }
  } else {
    // loop on the detsetvector<cluster> to find the right one
    for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = clusters.begin(); DSViter != clusters.end();
         DSViter++)
      if (DSViter->id() == det_id) {
        LogDebug("closestCluster") << " detset with the right detid. ";
        edmNew::DetSet<SiStripCluster>::const_iterator begin = DSViter->begin();
        edmNew::DetSet<SiStripCluster>::const_iterator end = DSViter->end();
        //find the cluster close to the hitStrip
        result.second = 1000.;
        for (edmNew::DetSet<SiStripCluster>::const_iterator iter = begin; iter != end; ++iter) {
          double dist = fabs(iter->barycenter() - hitStrip);
          if (dist < result.second) {
            result.second = dist;
            result.first = &(*iter);
          }
        }
        break;
      }
  }
  return result;
}
 
// ------------ method called to produce the data  ------------
void SiStripFineDelayHit::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  using namespace edm;
  // Retrieve commissioning information from "event summary"
  edm::Handle<SiStripEventSummary> runsummary;
  //   iEvent.getByLabel( inputModuleLabel_, runsummary );
  iEvent.getByToken(inputModuleToken_, runsummary);
  if (runsummary->runType() == sistrip::APV_LATENCY)
    mode_ = 2;  // LatencyScan
  else if (runsummary->runType() == sistrip::FINE_DELAY)
    mode_ = 1;  // DelayScan
  else {
    mode_ = 0;  //unknown
    return;
  }
 
  if (noTracking_) {
    produceNoTracking(iEvent, iSetup);
    return;
  }
  event_ = &iEvent;
  // container for the selected hits
  std::vector<edm::DetSet<SiStripRawDigi> > output;
  output.reserve(100);
  // access the tracks
  edm::Handle<reco::TrackCollection> trackCollection;
  //   iEvent.getByLabel(trackLabel_,trackCollection);
  iEvent.getByToken(trackCollectionToken_, trackCollection);
  const reco::TrackCollection* tracks = trackCollection.product();
  edm::ESHandle<TrackerGeometry> tracker;
  iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
  if (!tracks->empty()) {
    anglefinder_->init(iEvent, iSetup);
    LogDebug("produce") << "Found " << tracks->size() << " tracks.";
    // look at the hits if one needs them
    edm::Handle<edm::DetSetVector<SiStripDigi> > hits;
    const edm::DetSetVector<SiStripDigi>* hitSet = nullptr;
    if (homeMadeClusters_) {
      //       iEvent.getByLabel(digiLabel_,hits);
      iEvent.getByToken(digiToken_, hits);
      hitSet = hits.product();
    }
    // look at the clusters
    edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters;
    //     iEvent.getByLabel(clusterLabel_, clusters);
    iEvent.getByToken(clustersToken_, clusters);
    const edmNew::DetSetVector<SiStripCluster>* clusterSet = clusters.product();
    // look at the trajectories if they are in the event
    std::vector<Trajectory> trajVec;
    edm::Handle<std::vector<Trajectory> > TrajectoryCollection;
    //     iEvent.getByLabel(trackLabel_,TrajectoryCollection);
    iEvent.getByToken(trackToken_, TrajectoryCollection);
    trajVec = *(TrajectoryCollection.product());
    // Get TrackerTopology
    edm::ESHandle<TrackerTopology> tTopo;
    iSetup.get<TrackerTopologyRcd>().get(tTopo);
    // loop on tracks
    for (reco::TrackCollection::const_iterator itrack = tracks->begin(); itrack < tracks->end(); itrack++) {
      // first check the track Pt
      if ((itrack->px() * itrack->px() + itrack->py() * itrack->py() + itrack->pz() * itrack->pz()) < minTrackP2_)
        continue;
      // check that we have something in the layer we are interested in
      std::vector<std::pair<uint32_t, std::pair<double, double> > > intersections;
      if (mode_ == 1) {
        // Retrieve and decode commissioning information from "event summary"
        edm::Handle<SiStripEventSummary> summary;
        //         iEvent.getByLabel( inputModuleLabel_, summary );
        iEvent.getByToken(inputModuleToken_, summary);
        uint32_t layerCode = (const_cast<SiStripEventSummary*>(summary.product())->layerScanned()) >> 16;
        StripSubdetector::SubDetector subdet = StripSubdetector::TIB;
        if (((layerCode >> 6) & 0x3) == 0)
          subdet = StripSubdetector::TIB;
        else if (((layerCode >> 6) & 0x3) == 1)
          subdet = StripSubdetector::TOB;
        else if (((layerCode >> 6) & 0x3) == 2)
          subdet = StripSubdetector::TID;
        else if (((layerCode >> 6) & 0x3) == 3)
          subdet = StripSubdetector::TEC;
        int32_t layerIdx = (layerCode & 0xF) * (((layerCode >> 4) & 0x3) ? -1 : 1);
        intersections = detId(*tracker, tTopo.product(), &(*itrack), trajVec, subdet, layerIdx);
      } else {
        // for latency scans, no layer is specified -> no cut on detid
        intersections = detId(*tracker, tTopo.product(), &(*itrack), trajVec);
      }
      LogDebug("produce") << "  Found " << intersections.size() << " interesting intersections." << std::endl;
      for (std::vector<std::pair<uint32_t, std::pair<double, double> > >::iterator it = intersections.begin();
           it < intersections.end();
           it++) {
        std::pair<const SiStripCluster*, double> candidateCluster =
            closestCluster(*tracker, &(*itrack), it->first, *clusterSet, *hitSet);
        if (candidateCluster.first) {
          LogDebug("produce") << "    Found a cluster." << std::endl;
          // cut on the distance
          if (candidateCluster.second > maxClusterDistance_)
            continue;
          LogDebug("produce") << "    The cluster is close enough." << std::endl;
          // build the rawdigi corresponding to the leading strip and save it
          // here, only the leading strip is retained. All other rawdigis in the module are set to 0.
          const auto& amplitudes = candidateCluster.first->amplitudes();
          uint8_t leadingCharge = 0;
          uint8_t leadingStrip = candidateCluster.first->firstStrip();
          uint8_t leadingPosition = 0;
          for (auto amplit = amplitudes.begin(); amplit < amplitudes.end(); amplit++, leadingStrip++) {
            if (leadingCharge < *amplit) {
              leadingCharge = *amplit;
              leadingPosition = leadingStrip;
            }
          }
 
          // look for an existing detset
          std::vector<edm::DetSet<SiStripRawDigi> >::iterator newdsit = output.begin();
          for (; newdsit != output.end() && newdsit->detId() != connectionMap_[it->first]; ++newdsit) {
          }
          // if there is no detset yet, create it.
          if (newdsit == output.end()) {
            edm::DetSet<SiStripRawDigi> newds(connectionMap_[it->first]);
            output.push_back(newds);
            newdsit = output.end() - 1;
          }
 
          LogDebug("produce") << " New Hit...   TOF:" << it->second.first << ", charge: " << int(leadingCharge)
                              << " at " << int(leadingPosition) << "." << std::endl
                              << "Angular correction: " << it->second.second << " giving a final value of "
                              << int(leadingCharge * fabs(it->second.second))
                              << " for fed key = " << connectionMap_[it->first] << " (detid=" << it->first << ")";
          // apply corrections to the leading charge, but only if it has not saturated.
          if (leadingCharge < 255) {
            // correct the leading charge for the crossing angle
            leadingCharge = uint8_t(leadingCharge * fabs(it->second.second));
            // correct for module thickness for TEC and TOB
            if ((((it->first >> 25) & 0x7f) == 0xd) ||
                ((((it->first >> 25) & 0x7f) == 0xe) && (((it->first >> 5) & 0x7) > 4)))
              leadingCharge = uint8_t((leadingCharge * 0.64));
          }
          //code the time of flight in the digi
          unsigned int tof = abs(int(round(it->second.first * 10)));
          tof = tof > 255 ? 255 : tof;
          SiStripRawDigi newSiStrip(leadingCharge + (tof << 8));
          newdsit->push_back(newSiStrip);
          LogDebug("produce") << "New edm::DetSet<SiStripRawDigi> added.";
        }
        if (homeMadeClusters_)
          delete candidateCluster.first;  // we are owner of home-made clusters
      }
    }
  }
  // add the selected hits to the event.
  LogDebug("produce") << "Putting " << output.size() << " new hits in the event.";
  std::unique_ptr<edm::DetSetVector<SiStripRawDigi> > formatedOutput(new edm::DetSetVector<SiStripRawDigi>(output));
  iEvent.put(std::move(formatedOutput), "FineDelaySelection");
}
 
// Simple solution when tracking is not available/ not working
void SiStripFineDelayHit::produceNoTracking(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  event_ = &iEvent;
  // Get TrackerTopology
  edm::ESHandle<TrackerTopology> tTopo;
  iSetup.get<TrackerTopologyRcd>().get(tTopo);
  // container for the selected hits
  std::vector<edm::DetSet<SiStripRawDigi> > output;
  output.reserve(100);
  // Retrieve and decode commissioning information from "event summary"
  edm::Handle<SiStripEventSummary> summary;
  //   iEvent.getByLabel( inputModuleLabel_, summary );
  iEvent.getByToken(inputModuleToken_, summary);
  uint32_t layerCode = (const_cast<SiStripEventSummary*>(summary.product())->layerScanned()) >> 16;
  StripSubdetector::SubDetector subdet = StripSubdetector::TIB;
  if (((layerCode >> 6) & 0x3) == 0)
    subdet = StripSubdetector::TIB;
  else if (((layerCode >> 6) & 0x3) == 1)
    subdet = StripSubdetector::TOB;
  else if (((layerCode >> 6) & 0x3) == 2)
    subdet = StripSubdetector::TID;
  else if (((layerCode >> 6) & 0x3) == 3)
    subdet = StripSubdetector::TEC;
  int32_t layerIdx = (layerCode & 0xF) * (((layerCode >> 4) & 0x3) ? -1 : 1);
  DeviceMask mask = deviceMask(subdet, layerIdx, tTopo.product());
  // look at the clusters
  edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters;
  //   iEvent.getByLabel(clusterLabel_,clusters);
  iEvent.getByToken(clustersToken_, clusters);
  for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
       DSViter++) {
    // check that we are in the layer of interest
    if (mode_ == 1 && ((DSViter->id() & mask.first) != mask.second))
      continue;
    // iterate over clusters
    edmNew::DetSet<SiStripCluster>::const_iterator begin = DSViter->begin();
    edmNew::DetSet<SiStripCluster>::const_iterator end = DSViter->end();
    edm::DetSet<SiStripRawDigi> newds(connectionMap_[DSViter->id()]);
    for (edmNew::DetSet<SiStripCluster>::const_iterator iter = begin; iter != end; ++iter) {
      // build the rawdigi corresponding to the leading strip and save it
      // here, only the leading strip is retained. All other rawdigis in the module are set to 0.
      auto const& amplitudes = iter->amplitudes();
      uint8_t leadingCharge = 0;
      uint8_t leadingStrip = iter->firstStrip();
      uint8_t leadingPosition = 0;
      for (auto amplit = amplitudes.begin(); amplit < amplitudes.end(); amplit++, leadingStrip++) {
        if (leadingCharge < *amplit) {
          leadingCharge = *amplit;
          leadingPosition = leadingStrip;
        }
      }
      // apply some sanity cuts. This is needed since we don't use tracking to clean clusters
      // 1.5< noise <8
      // charge<250
      // 50 > s/n > 10
      edm::ESHandle<SiStripNoises> noiseHandle_;
      iSetup.get<SiStripNoisesRcd>().get(noiseHandle_);
      SiStripNoises::Range detNoiseRange = noiseHandle_->getRange(DSViter->id());
      float noise = noiseHandle_->getNoise(leadingPosition, detNoiseRange);
      if (noise < 1.5)
        continue;
      if (leadingCharge >= 250 || noise >= 8 || leadingCharge / noise > 50 || leadingCharge / noise < 10)
        continue;
      // apply some correction to the leading charge, but only if it has not saturated.
      if (leadingCharge < 255) {
        // correct for modulethickness for TEC and TOB
        if ((((((DSViter->id()) >> 25) & 0x7f) == 0xd) || ((((DSViter->id()) >> 25) & 0x7f) == 0xe)) &&
            ((((DSViter->id()) >> 5) & 0x7) > 4))
          leadingCharge = uint8_t((leadingCharge * 0.64));
      }
      //code the time of flight == 0 in the digi
      SiStripRawDigi newSiStrip(leadingCharge);
      newds.push_back(newSiStrip);
    }
    //store into the detsetvector
    output.push_back(newds);
    LogDebug("produce") << "New edm::DetSet<SiStripRawDigi> added with fedkey = " << std::hex << std::setfill('0')
                        << std::setw(8) << connectionMap_[DSViter->id()] << std::dec;
  }
  // add the selected hits to the event.
  LogDebug("produce") << "Putting " << output.size() << " new hits in the event.";
  std::unique_ptr<edm::DetSetVector<SiStripRawDigi> > formatedOutput(new edm::DetSetVector<SiStripRawDigi>(output));
  iEvent.put(std::move(formatedOutput), "FineDelaySelection");
}
 
// ------------ method called once each job just before starting event loop  ------------
void SiStripFineDelayHit::beginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
  // Retrieve FED cabling object
  edm::ESHandle<SiStripFedCabling> cabling;
  iSetup.get<SiStripFedCablingRcd>().get(cabling);
  auto feds = cabling->fedIds();
  for (auto fedid = feds.begin(); fedid < feds.end(); ++fedid) {
    auto connections = cabling->fedConnections(*fedid);
    for (std::vector<FedChannelConnection>::const_iterator conn = connections.begin(); conn < connections.end();
         ++conn) {
      /*
       SiStripFedKey key(conn->fedId(),
                         SiStripFedKey::feUnit(conn->fedCh()),
             SiStripFedKey::feChan(conn->fedCh()));
       connectionMap_[conn->detId()] = key.key();
     */
      // the key is computed using an alternate formula for performance reasons.
      connectionMap_[conn->detId()] = ((conn->fedId() & sistrip::invalid_) << 16) | (conn->fedCh() & sistrip::invalid_);
    }
  }
}