Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:08:35 

0001 // -*- C++ -*-
0002 //
0003 // Package:    SiStripFineDelayHit
0004 // Class:      SiStripFineDelayHit
0005 //
0006 /**\class SiStripFineDelayHit SiStripFineDelayHit.cc DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayHit.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Christophe DELAERE
0015 //         Created:  Fri Nov 17 10:52:42 CET 2006
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <utility>
0022 #include <vector>
0023 #include <algorithm>
0024 
0025 // user include files
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/Framework/interface/ConsumesCollector.h"
0030 #include "FWCore/Utilities/interface/InputTag.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "DataFormats/Common/interface/Ref.h"
0033 #include "DataFormats/DetId/interface/DetId.h"
0034 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0035 #include "DataFormats/TrackReco/interface/Track.h"
0036 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0037 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0038 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0039 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0040 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0041 #include "DataFormats/Candidate/interface/Candidate.h"
0042 #include "DataFormats/SiStripCommon/interface/ConstantsForRunType.h"
0043 #include "DataFormats/SiStripCommon/interface/SiStripFedKey.h"
0044 #include "CondFormats/SiStripObjects/interface/FedChannelConnection.h"
0045 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0046 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0047 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0048 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0049 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0050 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0051 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0052 #include "Geometry/CommonTopologies/interface/Topology.h"
0053 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0054 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0055 #include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayHit.h"
0056 #include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTLA.h"
0057 #include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTOF.h"
0058 
0059 // ROOT includes
0060 #include "TMath.h"
0061 
0062 //
0063 // constructors and destructor
0064 //
0065 SiStripFineDelayHit::SiStripFineDelayHit(const edm::ParameterSet& iConfig) : event_(nullptr) {
0066   //register your products
0067   produces<edm::DetSetVector<SiStripRawDigi> >("FineDelaySelection");
0068   //now do what ever other initialization is needed
0069   anglefinder_ = new SiStripFineDelayTLA(iConfig, consumesCollector());
0070   cosmic_ = iConfig.getParameter<bool>("cosmic");
0071   field_ = iConfig.getParameter<bool>("MagneticField");
0072   maxAngle_ = iConfig.getParameter<double>("MaxTrackAngle");
0073   minTrackP2_ = iConfig.getParameter<double>("MinTrackMomentum") * iConfig.getParameter<double>("MinTrackMomentum");
0074   maxClusterDistance_ = iConfig.getParameter<double>("MaxClusterDistance");
0075   /*
0076    clusterLabel_ = iConfig.getParameter<edm::InputTag>("ClustersLabel");
0077    trackLabel_ = iConfig.getParameter<edm::InputTag>("TracksLabel");
0078    seedLabel_  = iConfig.getParameter<edm::InputTag>("SeedsLabel");
0079    inputModuleLabel_ = iConfig.getParameter<edm::InputTag>( "InputModuleLabel" ) ;
0080    digiLabel_ = iConfig.getParameter<edm::InputTag>("DigiLabel");
0081    */
0082   clustersToken_ =
0083       consumes<edmNew::DetSetVector<SiStripCluster> >(iConfig.getParameter<edm::InputTag>("ClustersLabel"));
0084   trackToken_ = consumes<std::vector<Trajectory> >(iConfig.getParameter<edm::InputTag>("TracksLabel"));
0085   trackCollectionToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("TracksLabel"));
0086   seedcollToken_ = consumes<TrajectorySeedCollection>(iConfig.getParameter<edm::InputTag>("SeedsLabel"));
0087   inputModuleToken_ = consumes<SiStripEventSummary>(iConfig.getParameter<edm::InputTag>("InputModuleLabel"));
0088   digiToken_ = consumes<edm::DetSetVector<SiStripDigi> >(iConfig.getParameter<edm::InputTag>("DigiLabel"));
0089 
0090   homeMadeClusters_ = iConfig.getParameter<bool>("NoClustering");
0091   explorationWindow_ = iConfig.getParameter<uint32_t>("ExplorationWindow");
0092   noTracking_ = iConfig.getParameter<bool>("NoTracking");
0093   mode_ = 0;
0094 
0095   tkGeomToken_ = esConsumes();
0096   tTopoToken_ = esConsumes();
0097   fedCablingToken_ = esConsumes<edm::Transition::BeginRun>();
0098   noiseToken_ = esConsumes();
0099 }
0100 
0101 SiStripFineDelayHit::~SiStripFineDelayHit() {
0102   // do anything here that needs to be done at desctruction time
0103   // (e.g. close files, deallocate resources etc.)
0104   delete anglefinder_;
0105 }
0106 
0107 //
0108 // member functions
0109 //
0110 SiStripFineDelayHit::DeviceMask SiStripFineDelayHit::deviceMask(const StripSubdetector::SubDetector subdet,
0111                                                                 const int substructure,
0112                                                                 const TrackerTopology* tkrTopo) {
0113   uint32_t rootDetId = 0;
0114   uint32_t maskDetId = 0;
0115 
0116   switch (subdet) {
0117     case StripSubdetector::TIB: {
0118       rootDetId = tkrTopo->tibDetId(substructure, 0, 0, 0, 0, 0).rawId();
0119       maskDetId = tkrTopo->tibDetId(15, 0, 0, 0, 0, 0).rawId();
0120       break;
0121     }
0122     case StripSubdetector::TID: {
0123       rootDetId = tkrTopo->tidDetId(substructure > 0 ? 2 : 1, abs(substructure), 0, 0, 0, 0).rawId();
0124       maskDetId = tkrTopo->tidDetId(3, 15, 0, 0, 0, 0).rawId();
0125       break;
0126     }
0127     case StripSubdetector::TOB: {
0128       rootDetId = tkrTopo->tobDetId(substructure, 0, 0, 0, 0).rawId();
0129       maskDetId = tkrTopo->tobDetId(15, 0, 0, 0, 0).rawId();
0130       break;
0131     }
0132     case StripSubdetector::TEC: {
0133       rootDetId = tkrTopo->tecDetId(substructure > 0 ? 2 : 1, abs(substructure), 0, 0, 0, 0, 0).rawId();
0134       maskDetId = tkrTopo->tecDetId(3, 15, 0, 0, 0, 0, 0).rawId();
0135       break;
0136     }
0137     default:
0138       break;
0139   }
0140   return std::make_pair(maskDetId, rootDetId);
0141 }
0142 
0143 std::vector<std::pair<uint32_t, std::pair<double, double> > > SiStripFineDelayHit::detId(
0144     const TrackerGeometry& tracker,
0145     const TrackerTopology* tkrTopo,
0146     const reco::Track* tk,
0147     const std::vector<Trajectory>& trajVec,
0148     const StripSubdetector::SubDetector subdet,
0149     const int substructure) {
0150   if (substructure == 0xff)
0151     return detId(tracker, tkrTopo, tk, trajVec, 0, 0);
0152   // first determine the root detId we are looking for
0153   DeviceMask mask = deviceMask(subdet, substructure, tkrTopo);
0154   // then call the method that loops on recHits
0155   return detId(tracker, tkrTopo, tk, trajVec, mask.first, mask.second);
0156 }
0157 
0158 std::vector<std::pair<uint32_t, std::pair<double, double> > > SiStripFineDelayHit::detId(
0159     const TrackerGeometry& tracker,
0160     const TrackerTopology* tkrTopo,
0161     const reco::Track* tk,
0162     const std::vector<Trajectory>& trajVec,
0163     const uint32_t& maskDetId,
0164     const uint32_t& rootDetId) {
0165   bool onDisk = ((maskDetId == tkrTopo->tidDetId(3, 15, 0, 0, 0, 0).rawId()) ||
0166                  (maskDetId == tkrTopo->tecDetId(3, 15, 0, 0, 0, 0, 0).rawId()));
0167   std::vector<std::pair<uint32_t, std::pair<double, double> > > result;
0168   std::vector<uint32_t> usedDetids;
0169   // now loop on recHits to find the right detId plus the track local angle
0170   std::vector<std::pair<std::pair<DetId, LocalPoint>, float> > hitangle;
0171   if (!cosmic_) {
0172     // use trajectories in event.
0173     // we have first to find the right trajectory for the considered track.
0174     for (std::vector<Trajectory>::const_iterator traj = trajVec.begin(); traj < trajVec.end(); ++traj) {
0175       if (((traj->lastMeasurement().recHit()->geographicalId().rawId() ==
0176             (*(tk->recHitsEnd() - 1))->geographicalId().rawId()) &&
0177            (traj->lastMeasurement().recHit()->localPosition().x() == (*(tk->recHitsEnd() - 1))->localPosition().x())) ||
0178           ((traj->firstMeasurement().recHit()->geographicalId().rawId() ==
0179             (*(tk->recHitsEnd() - 1))->geographicalId().rawId()) &&
0180            (traj->firstMeasurement().recHit()->localPosition().x() == (*(tk->recHitsEnd() - 1))->localPosition().x()))) {
0181         hitangle = anglefinder_->findtrackangle(*traj);
0182         break;
0183       }
0184     }
0185   } else {
0186     edm::Handle<TrajectorySeedCollection> seedcoll;
0187     //    event_->getByLabel(seedLabel_,seedcoll);
0188     event_->getByToken(seedcollToken_, seedcoll);
0189     // use trajectories in event.
0190     hitangle = anglefinder_->findtrackangle(trajVec);
0191   }
0192   LogDebug("DetId") << "number of hits for the track: " << hitangle.size();
0193   std::vector<std::pair<std::pair<DetId, LocalPoint>, float> >::iterator iter;
0194   // select the interesting DetIds, based on the ID and TLA
0195   for (iter = hitangle.begin(); iter != hitangle.end(); iter++) {
0196     // check the detId.
0197     // if substructure was 0xff, then maskDetId and rootDetId == 0
0198     // this implies all detids are accepted. (also if maskDetId=rootDetId=0 explicitely).
0199     // That "unusual" mode of operation allows to analyze also Latency scans
0200     LogDebug("DetId") << "check the detid: " << std::hex << (iter->first.first.rawId()) << " vs " << rootDetId
0201                       << " with a mask of " << maskDetId << std::dec << std::endl;
0202 
0203     if (((iter->first.first.rawId() & maskDetId) != rootDetId))
0204       continue;
0205     if (std::find(usedDetids.begin(), usedDetids.end(), iter->first.first.rawId()) != usedDetids.end())
0206       continue;
0207     // check the local angle (extended to the equivalent angle correction)
0208     LogDebug("DetId") << "check the angle: " << fabs((iter->second));
0209     if (1 - fabs(fabs(iter->second) - 1) < cos(maxAngle_ / 180. * TMath::Pi()))
0210       continue;
0211     // returns the detid + the time of flight to there
0212     std::pair<uint32_t, std::pair<double, double> > el;
0213     std::pair<double, double> subel;
0214     el.first = iter->first.first.rawId();
0215     // here, we compute the TOF.
0216     // For cosmics, some track parameters are missing. Parameters are recomputed.
0217     // for our calculation, the track momemtum at any point is enough:
0218     // only used without B field or for the sign of Pz.
0219     double trackParameters[5];
0220     for (int i = 0; i < 5; i++)
0221       trackParameters[i] = tk->parameters()[i];
0222     if (cosmic_)
0223       SiStripFineDelayTOF::trackParameters(*tk, trackParameters);
0224     double hit[3];
0225     const GeomDetUnit* det(tracker.idToDetUnit(iter->first.first));
0226     Surface::GlobalPoint gp = det->surface().toGlobal(iter->first.second);
0227     hit[0] = gp.x();
0228     hit[1] = gp.y();
0229     hit[2] = gp.z();
0230     double phit[3];
0231     phit[0] = tk->momentum().x();
0232     phit[1] = tk->momentum().y();
0233     phit[2] = tk->momentum().z();
0234     subel.first = SiStripFineDelayTOF::timeOfFlight(cosmic_, field_, trackParameters, hit, phit, onDisk);
0235     subel.second = iter->second;
0236     el.second = subel;
0237     // returns the detid + TOF
0238     result.push_back(el);
0239     usedDetids.push_back(el.first);
0240   }
0241   return result;
0242 }
0243 
0244 bool SiStripFineDelayHit::rechit(reco::Track* tk, uint32_t det_id) {
0245   for (trackingRecHit_iterator it = tk->recHitsBegin(); it != tk->recHitsEnd(); it++)
0246     if ((*it)->geographicalId().rawId() == det_id) {
0247       return (*it)->isValid();
0248       break;
0249     }
0250   return false;
0251 }
0252 
0253 // VI January 2012: FIXME
0254 // do not understand what is going on here: each hit has a cluster: by definition will be the closest!
0255 std::pair<const SiStripCluster*, double> SiStripFineDelayHit::closestCluster(
0256     const TrackerGeometry& tracker,
0257     const reco::Track* tk,
0258     const uint32_t& det_id,
0259     const edmNew::DetSetVector<SiStripCluster>& clusters,
0260     const edm::DetSetVector<SiStripDigi>& hits) {
0261   std::pair<const SiStripCluster*, double> result(nullptr, 0.);
0262   double hitStrip = -1;
0263   int nstrips = -1;
0264   // localize the crossing point of the track on the module
0265   for (trackingRecHit_iterator it = tk->recHitsBegin(); it != tk->recHitsEnd(); it++) {
0266     LogDebug("closestCluster") << "(*it)->geographicalId().rawId() vs det_id" << (*it)->geographicalId().rawId() << " "
0267                                << det_id;
0268     //handle the mono rechits
0269     if ((*it)->geographicalId().rawId() == det_id) {
0270       if (!(*it)->isValid())
0271         continue;
0272       LogDebug("closestCluster") << " using the single mono hit";
0273       LocalPoint lp = (*it)->localPosition();
0274       const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker.idToDet((*it)->geographicalId()));
0275       MeasurementPoint p = gdu->topology().measurementPosition(lp);
0276       hitStrip = p.x();
0277       nstrips = (dynamic_cast<const StripTopology*>(&(gdu->topology())))->nstrips();
0278       break;
0279     }
0280     /* FIXME: local position is not there anymore...
0281     //handle stereo part of matched hits
0282     //one could try to cast to SiStripMatchedRecHit2D but it is faster to look at the detid
0283     else if((det_id - (*it)->geographicalId().rawId())==1) {
0284       const SiStripMatchedRecHit2D* hit2D = dynamic_cast<const SiStripMatchedRecHit2D*>(&(**it));
0285       if(!hit2D) continue; // this is a security that should never trigger
0286       const SiStripRecHit2D* stereo = hit2D->stereoHit();
0287       if(!stereo) continue; // this is a security that should never trigger
0288       if(!stereo->isValid()) continue;
0289       LogDebug("closestCluster") << " using the stereo hit";
0290       LocalPoint lp = stereo->localPosition();
0291       const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker.idToDet(stereo->geographicalId()));
0292       MeasurementPoint p = gdu->topology().measurementPosition(lp);
0293       hitStrip = p.x();
0294       nstrips = (dynamic_cast<const StripTopology*>(&(gdu->topology())))->nstrips();
0295       break;
0296     }
0297     //handle mono part of matched hits
0298     //one could try to cast to SiStripMatchedRecHit2D but it is faster to look at the detid
0299     else if((det_id - (*it)->geographicalId().rawId())==2) {
0300       const SiStripMatchedRecHit2D* hit2D = dynamic_cast<const SiStripMatchedRecHit2D*>(&(**it));
0301       if(!hit2D) continue; // this is a security that should never trigger
0302       const SiStripRecHit2D* mono = hit2D->monoHit();
0303       if(!mono) continue; // this is a security that should never trigger
0304       if(!mono->isValid()) continue;
0305       LogDebug("closestCluster") << " using the mono hit";
0306       LocalPoint lp = mono->localPosition();
0307       const GeomDetUnit* gdu = static_cast<const GeomDetUnit*>(tracker.idToDet(mono->geographicalId()));
0308       MeasurementPoint p = gdu->topology().measurementPosition(lp);
0309       hitStrip = p.x();
0310       nstrips = (dynamic_cast<const StripTopology*>(&(gdu->topology())))->nstrips();
0311       break;
0312     }
0313     */
0314   }
0315   LogDebug("closestCluster") << " hit strip = " << hitStrip;
0316   if (hitStrip < 0)
0317     return result;
0318   if (homeMadeClusters_) {
0319     // take the list of digis on the module
0320     for (edm::DetSetVector<SiStripDigi>::const_iterator DSViter = hits.begin(); DSViter != hits.end(); DSViter++) {
0321       if (DSViter->id == det_id) {
0322         // loop from hitstrip-n to hitstrip+n (explorationWindow_) and select the highest strip
0323         int minStrip = int(round(hitStrip)) - explorationWindow_;
0324         minStrip = minStrip < 0 ? 0 : minStrip;
0325         int maxStrip = int(round(hitStrip)) + explorationWindow_ + 1;
0326         maxStrip = maxStrip >= nstrips ? nstrips - 1 : maxStrip;
0327         edm::DetSet<SiStripDigi>::const_iterator rangeStart = DSViter->end();
0328         edm::DetSet<SiStripDigi>::const_iterator rangeStop = DSViter->end();
0329         for (edm::DetSet<SiStripDigi>::const_iterator digiIt = DSViter->begin(); digiIt != DSViter->end(); ++digiIt) {
0330           if (digiIt->strip() >= minStrip && rangeStart == DSViter->end())
0331             rangeStart = digiIt;
0332           if (digiIt->strip() <= maxStrip)
0333             rangeStop = digiIt;
0334         }
0335         if (rangeStart != DSViter->end()) {
0336           if (rangeStop != DSViter->end())
0337             ++rangeStop;
0338           // build a fake cluster
0339           LogDebug("closestCluster") << "build a fake cluster ";
0340           SiStripCluster* newCluster =
0341               new SiStripCluster(SiStripCluster::SiStripDigiRange(rangeStart, rangeStop));  // /!\ ownership transfered
0342           result.first = newCluster;
0343           result.second = fabs(newCluster->barycenter() - hitStrip);
0344         }
0345         break;
0346       }
0347     }
0348   } else {
0349     // loop on the detsetvector<cluster> to find the right one
0350     for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = clusters.begin(); DSViter != clusters.end();
0351          DSViter++)
0352       if (DSViter->id() == det_id) {
0353         LogDebug("closestCluster") << " detset with the right detid. ";
0354         edmNew::DetSet<SiStripCluster>::const_iterator begin = DSViter->begin();
0355         edmNew::DetSet<SiStripCluster>::const_iterator end = DSViter->end();
0356         //find the cluster close to the hitStrip
0357         result.second = 1000.;
0358         for (edmNew::DetSet<SiStripCluster>::const_iterator iter = begin; iter != end; ++iter) {
0359           double dist = fabs(iter->barycenter() - hitStrip);
0360           if (dist < result.second) {
0361             result.second = dist;
0362             result.first = &(*iter);
0363           }
0364         }
0365         break;
0366       }
0367   }
0368   return result;
0369 }
0370 
0371 // ------------ method called to produce the data  ------------
0372 void SiStripFineDelayHit::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0373   using namespace edm;
0374   // Retrieve commissioning information from "event summary"
0375   edm::Handle<SiStripEventSummary> runsummary;
0376   //   iEvent.getByLabel( inputModuleLabel_, runsummary );
0377   iEvent.getByToken(inputModuleToken_, runsummary);
0378   if (runsummary->runType() == sistrip::APV_LATENCY)
0379     mode_ = 2;  // LatencyScan
0380   else if (runsummary->runType() == sistrip::FINE_DELAY)
0381     mode_ = 1;  // DelayScan
0382   else {
0383     mode_ = 0;  //unknown
0384     return;
0385   }
0386 
0387   if (noTracking_) {
0388     produceNoTracking(iEvent, iSetup);
0389     return;
0390   }
0391   event_ = &iEvent;
0392   // container for the selected hits
0393   std::vector<edm::DetSet<SiStripRawDigi> > output;
0394   output.reserve(100);
0395   // access the tracks
0396   edm::Handle<reco::TrackCollection> trackCollection;
0397   //   iEvent.getByLabel(trackLabel_,trackCollection);
0398   iEvent.getByToken(trackCollectionToken_, trackCollection);
0399   const reco::TrackCollection* tracks = trackCollection.product();
0400   const auto& tracker = iSetup.getData(tkGeomToken_);
0401   if (!tracks->empty()) {
0402     anglefinder_->init(iEvent, iSetup);
0403     LogDebug("produce") << "Found " << tracks->size() << " tracks.";
0404     // look at the hits if one needs them
0405     edm::Handle<edm::DetSetVector<SiStripDigi> > hits;
0406     const edm::DetSetVector<SiStripDigi>* hitSet = nullptr;
0407     if (homeMadeClusters_) {
0408       //       iEvent.getByLabel(digiLabel_,hits);
0409       iEvent.getByToken(digiToken_, hits);
0410       hitSet = hits.product();
0411     }
0412     // look at the clusters
0413     edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters;
0414     //     iEvent.getByLabel(clusterLabel_, clusters);
0415     iEvent.getByToken(clustersToken_, clusters);
0416     const edmNew::DetSetVector<SiStripCluster>* clusterSet = clusters.product();
0417     // look at the trajectories if they are in the event
0418     std::vector<Trajectory> trajVec;
0419     edm::Handle<std::vector<Trajectory> > TrajectoryCollection;
0420     //     iEvent.getByLabel(trackLabel_,TrajectoryCollection);
0421     iEvent.getByToken(trackToken_, TrajectoryCollection);
0422     trajVec = *(TrajectoryCollection.product());
0423     // Get TrackerTopology
0424     const auto tTopo = &iSetup.getData(tTopoToken_);
0425     // loop on tracks
0426     for (reco::TrackCollection::const_iterator itrack = tracks->begin(); itrack < tracks->end(); itrack++) {
0427       // first check the track Pt
0428       if ((itrack->px() * itrack->px() + itrack->py() * itrack->py() + itrack->pz() * itrack->pz()) < minTrackP2_)
0429         continue;
0430       // check that we have something in the layer we are interested in
0431       std::vector<std::pair<uint32_t, std::pair<double, double> > > intersections;
0432       if (mode_ == 1) {
0433         // Retrieve and decode commissioning information from "event summary"
0434         edm::Handle<SiStripEventSummary> summary;
0435         //         iEvent.getByLabel( inputModuleLabel_, summary );
0436         iEvent.getByToken(inputModuleToken_, summary);
0437         uint32_t layerCode = (const_cast<SiStripEventSummary*>(summary.product())->layerScanned()) >> 16;
0438         StripSubdetector::SubDetector subdet = StripSubdetector::TIB;
0439         if (((layerCode >> 6) & 0x3) == 0)
0440           subdet = StripSubdetector::TIB;
0441         else if (((layerCode >> 6) & 0x3) == 1)
0442           subdet = StripSubdetector::TOB;
0443         else if (((layerCode >> 6) & 0x3) == 2)
0444           subdet = StripSubdetector::TID;
0445         else if (((layerCode >> 6) & 0x3) == 3)
0446           subdet = StripSubdetector::TEC;
0447         int32_t layerIdx = (layerCode & 0xF) * (((layerCode >> 4) & 0x3) ? -1 : 1);
0448         intersections = detId(tracker, tTopo, &(*itrack), trajVec, subdet, layerIdx);
0449       } else {
0450         // for latency scans, no layer is specified -> no cut on detid
0451         intersections = detId(tracker, tTopo, &(*itrack), trajVec);
0452       }
0453       LogDebug("produce") << "  Found " << intersections.size() << " interesting intersections." << std::endl;
0454       for (std::vector<std::pair<uint32_t, std::pair<double, double> > >::iterator it = intersections.begin();
0455            it < intersections.end();
0456            it++) {
0457         std::pair<const SiStripCluster*, double> candidateCluster =
0458             closestCluster(tracker, &(*itrack), it->first, *clusterSet, *hitSet);
0459         if (candidateCluster.first) {
0460           LogDebug("produce") << "    Found a cluster." << std::endl;
0461           // cut on the distance
0462           if (candidateCluster.second > maxClusterDistance_)
0463             continue;
0464           LogDebug("produce") << "    The cluster is close enough." << std::endl;
0465           // build the rawdigi corresponding to the leading strip and save it
0466           // here, only the leading strip is retained. All other rawdigis in the module are set to 0.
0467           const auto& amplitudes = candidateCluster.first->amplitudes();
0468           uint8_t leadingCharge = 0;
0469           uint8_t leadingStrip = candidateCluster.first->firstStrip();
0470           uint8_t leadingPosition = 0;
0471           for (auto amplit = amplitudes.begin(); amplit < amplitudes.end(); amplit++, leadingStrip++) {
0472             if (leadingCharge < *amplit) {
0473               leadingCharge = *amplit;
0474               leadingPosition = leadingStrip;
0475             }
0476           }
0477 
0478           // look for an existing detset
0479           std::vector<edm::DetSet<SiStripRawDigi> >::iterator newdsit = output.begin();
0480           for (; newdsit != output.end() && newdsit->detId() != connectionMap_[it->first]; ++newdsit) {
0481           }
0482           // if there is no detset yet, create it.
0483           if (newdsit == output.end()) {
0484             edm::DetSet<SiStripRawDigi> newds(connectionMap_[it->first]);
0485             output.push_back(newds);
0486             newdsit = output.end() - 1;
0487           }
0488 
0489           LogDebug("produce") << " New Hit...   TOF:" << it->second.first << ", charge: " << int(leadingCharge)
0490                               << " at " << int(leadingPosition) << "." << std::endl
0491                               << "Angular correction: " << it->second.second << " giving a final value of "
0492                               << int(leadingCharge * fabs(it->second.second))
0493                               << " for fed key = " << connectionMap_[it->first] << " (detid=" << it->first << ")";
0494           // apply corrections to the leading charge, but only if it has not saturated.
0495           if (leadingCharge < 255) {
0496             // correct the leading charge for the crossing angle
0497             leadingCharge = uint8_t(leadingCharge * fabs(it->second.second));
0498             // correct for module thickness for TEC and TOB
0499             if ((((it->first >> 25) & 0x7f) == 0xd) ||
0500                 ((((it->first >> 25) & 0x7f) == 0xe) && (((it->first >> 5) & 0x7) > 4)))
0501               leadingCharge = uint8_t((leadingCharge * 0.64));
0502           }
0503           //code the time of flight in the digi
0504           unsigned int tof = abs(int(round(it->second.first * 10)));
0505           tof = tof > 255 ? 255 : tof;
0506           SiStripRawDigi newSiStrip(leadingCharge + (tof << 8));
0507           newdsit->push_back(newSiStrip);
0508           LogDebug("produce") << "New edm::DetSet<SiStripRawDigi> added.";
0509         }
0510         if (homeMadeClusters_)
0511           delete candidateCluster.first;  // we are owner of home-made clusters
0512       }
0513     }
0514   }
0515   // add the selected hits to the event.
0516   LogDebug("produce") << "Putting " << output.size() << " new hits in the event.";
0517   std::unique_ptr<edm::DetSetVector<SiStripRawDigi> > formatedOutput(new edm::DetSetVector<SiStripRawDigi>(output));
0518   iEvent.put(std::move(formatedOutput), "FineDelaySelection");
0519 }
0520 
0521 // Simple solution when tracking is not available/ not working
0522 void SiStripFineDelayHit::produceNoTracking(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0523   event_ = &iEvent;
0524   // Get TrackerTopology
0525   const auto tTopo = &iSetup.getData(tTopoToken_);
0526   // container for the selected hits
0527   std::vector<edm::DetSet<SiStripRawDigi> > output;
0528   output.reserve(100);
0529   // Retrieve and decode commissioning information from "event summary"
0530   edm::Handle<SiStripEventSummary> summary;
0531   //   iEvent.getByLabel( inputModuleLabel_, summary );
0532   iEvent.getByToken(inputModuleToken_, summary);
0533   uint32_t layerCode = (const_cast<SiStripEventSummary*>(summary.product())->layerScanned()) >> 16;
0534   StripSubdetector::SubDetector subdet = StripSubdetector::TIB;
0535   if (((layerCode >> 6) & 0x3) == 0)
0536     subdet = StripSubdetector::TIB;
0537   else if (((layerCode >> 6) & 0x3) == 1)
0538     subdet = StripSubdetector::TOB;
0539   else if (((layerCode >> 6) & 0x3) == 2)
0540     subdet = StripSubdetector::TID;
0541   else if (((layerCode >> 6) & 0x3) == 3)
0542     subdet = StripSubdetector::TEC;
0543   int32_t layerIdx = (layerCode & 0xF) * (((layerCode >> 4) & 0x3) ? -1 : 1);
0544   DeviceMask mask = deviceMask(subdet, layerIdx, tTopo);
0545   // look at the clusters
0546   edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters;
0547   //   iEvent.getByLabel(clusterLabel_,clusters);
0548   iEvent.getByToken(clustersToken_, clusters);
0549   for (edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter = clusters->begin(); DSViter != clusters->end();
0550        DSViter++) {
0551     // check that we are in the layer of interest
0552     if (mode_ == 1 && ((DSViter->id() & mask.first) != mask.second))
0553       continue;
0554     // iterate over clusters
0555     edmNew::DetSet<SiStripCluster>::const_iterator begin = DSViter->begin();
0556     edmNew::DetSet<SiStripCluster>::const_iterator end = DSViter->end();
0557     edm::DetSet<SiStripRawDigi> newds(connectionMap_[DSViter->id()]);
0558     for (edmNew::DetSet<SiStripCluster>::const_iterator iter = begin; iter != end; ++iter) {
0559       // build the rawdigi corresponding to the leading strip and save it
0560       // here, only the leading strip is retained. All other rawdigis in the module are set to 0.
0561       auto const& amplitudes = iter->amplitudes();
0562       uint8_t leadingCharge = 0;
0563       uint8_t leadingStrip = iter->firstStrip();
0564       uint8_t leadingPosition = 0;
0565       for (auto amplit = amplitudes.begin(); amplit < amplitudes.end(); amplit++, leadingStrip++) {
0566         if (leadingCharge < *amplit) {
0567           leadingCharge = *amplit;
0568           leadingPosition = leadingStrip;
0569         }
0570       }
0571       // apply some sanity cuts. This is needed since we don't use tracking to clean clusters
0572       // 1.5< noise <8
0573       // charge<250
0574       // 50 > s/n > 10
0575       const auto& noises = iSetup.getData(noiseToken_);
0576       SiStripNoises::Range detNoiseRange = noises.getRange(DSViter->id());
0577       float noise = noises.getNoise(leadingPosition, detNoiseRange);
0578       if (noise < 1.5)
0579         continue;
0580       if (leadingCharge >= 250 || noise >= 8 || leadingCharge / noise > 50 || leadingCharge / noise < 10)
0581         continue;
0582       // apply some correction to the leading charge, but only if it has not saturated.
0583       if (leadingCharge < 255) {
0584         // correct for modulethickness for TEC and TOB
0585         if ((((((DSViter->id()) >> 25) & 0x7f) == 0xd) || ((((DSViter->id()) >> 25) & 0x7f) == 0xe)) &&
0586             ((((DSViter->id()) >> 5) & 0x7) > 4))
0587           leadingCharge = uint8_t((leadingCharge * 0.64));
0588       }
0589       //code the time of flight == 0 in the digi
0590       SiStripRawDigi newSiStrip(leadingCharge);
0591       newds.push_back(newSiStrip);
0592     }
0593     //store into the detsetvector
0594     output.push_back(newds);
0595     LogDebug("produce") << "New edm::DetSet<SiStripRawDigi> added with fedkey = " << std::hex << std::setfill('0')
0596                         << std::setw(8) << connectionMap_[DSViter->id()] << std::dec;
0597   }
0598   // add the selected hits to the event.
0599   LogDebug("produce") << "Putting " << output.size() << " new hits in the event.";
0600   std::unique_ptr<edm::DetSetVector<SiStripRawDigi> > formatedOutput(new edm::DetSetVector<SiStripRawDigi>(output));
0601   iEvent.put(std::move(formatedOutput), "FineDelaySelection");
0602 }
0603 
0604 // ------------ method called once each job just before starting event loop  ------------
0605 void SiStripFineDelayHit::beginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0606   // Retrieve FED cabling object
0607   const auto& cabling = iSetup.getData(fedCablingToken_);
0608   auto feds = cabling.fedIds();
0609   for (auto fedid = feds.begin(); fedid < feds.end(); ++fedid) {
0610     auto connections = cabling.fedConnections(*fedid);
0611     for (std::vector<FedChannelConnection>::const_iterator conn = connections.begin(); conn < connections.end();
0612          ++conn) {
0613       /*
0614        SiStripFedKey key(conn->fedId(),
0615                          SiStripFedKey::feUnit(conn->fedCh()),
0616              SiStripFedKey::feChan(conn->fedCh()));
0617        connectionMap_[conn->detId()] = key.key();
0618      */
0619       // the key is computed using an alternate formula for performance reasons.
0620       connectionMap_[conn->detId()] = ((conn->fedId() & sistrip::invalid_) << 16) | (conn->fedCh() & sistrip::invalid_);
0621     }
0622   }
0623 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!