Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:21:57

0001 //////////////////////////
0002 //  Producer by Anders  //
0003 //     and Emmanuele    //
0004 //    july 2012 @ CU    //
0005 //////////////////////////
0006 
0007 ////////////////////
0008 // FRAMEWORK HEADERS
0009 #include "FWCore/PluginManager/interface/ModuleDef.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 //
0012 #include "FWCore/Framework/interface/one/EDProducer.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/ParameterSet/interface/FileInPath.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 //
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 
0023 ///////////////////////
0024 // DATA FORMATS HEADERS
0025 #include "DataFormats/Common/interface/Handle.h"
0026 #include "DataFormats/Common/interface/Ref.h"
0027 //
0028 #include "DataFormats/DetId/interface/DetId.h"
0029 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
0030 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
0031 #include "DataFormats/Common/interface/DetSetVector.h"
0032 #include "DataFormats/L1TrackTrigger/interface/TTBV.h"
0033 //
0034 #include "L1Trigger/TrackFindingTracklet/interface/SLHCEvent.h"
0035 #include "L1Trigger/TrackFindingTracklet/interface/L1TStub.h"
0036 
0037 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0038 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0039 #include "SimDataFormats/Track/interface/SimTrack.h"
0040 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0041 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0042 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0043 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0044 #include "SimTracker/TrackTriggerAssociation/interface/TTClusterAssociationMap.h"
0045 //
0046 #include "DataFormats/Math/interface/LorentzVector.h"
0047 #include "DataFormats/Math/interface/Vector3D.h"
0048 //
0049 #include "DataFormats/L1TrackTrigger/interface/TTCluster.h"
0050 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
0051 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0052 #include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h"
0053 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0054 #include "DataFormats/L1TrackTrigger/interface/TTDTC.h"
0055 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0056 //
0057 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0058 #include "DataFormats/Candidate/interface/Candidate.h"
0059 //
0060 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0061 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0062 #include "DataFormats/SiPixelDetId/interface/PixelChannelIdentifier.h"
0063 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0064 
0065 ////////////////////////////
0066 // DETECTOR GEOMETRY HEADERS
0067 #include "MagneticField/Engine/interface/MagneticField.h"
0068 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0069 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0070 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0071 #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h"
0072 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyBuilder.h"
0073 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0074 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0075 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0076 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0077 //
0078 #include "Geometry/Records/interface/StackedTrackerGeometryRecord.h"
0079 
0080 ///////////////
0081 // Tracklet emulation
0082 #include "L1Trigger/TrackFindingTracklet/interface/Settings.h"
0083 #include "L1Trigger/TrackFindingTracklet/interface/Sector.h"
0084 #include "L1Trigger/TrackFindingTracklet/interface/Track.h"
0085 #include "L1Trigger/TrackFindingTracklet/interface/TrackletEventProcessor.h"
0086 #include "L1Trigger/TrackFindingTracklet/interface/ChannelAssignment.h"
0087 #include "L1Trigger/TrackFindingTracklet/interface/Tracklet.h"
0088 #include "L1Trigger/TrackFindingTracklet/interface/Residual.h"
0089 #include "L1Trigger/TrackFindingTracklet/interface/Stub.h"
0090 #include "L1Trigger/TrackFindingTracklet/interface/StubKiller.h"
0091 #include "L1Trigger/TrackFindingTracklet/interface/StubStreamData.h"
0092 #include "L1Trigger/TrackFindingTracklet/interface/HitPatternHelper.h"
0093 
0094 ////////////////
0095 // PHYSICS TOOLS
0096 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0097 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0098 
0099 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
0100 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0101 
0102 #include "L1Trigger/TrackTrigger/interface/StubPtConsistency.h"
0103 #include "L1Trigger/TrackTrigger/interface/L1TrackQuality.h"
0104 
0105 //////////////
0106 // STD HEADERS
0107 #include <memory>
0108 #include <string>
0109 #include <iostream>
0110 #include <fstream>
0111 
0112 //////////////
0113 // NAMESPACES
0114 using namespace edm;
0115 using namespace std;
0116 using namespace tt;
0117 using namespace trklet;
0118 
0119 //////////////////////////////
0120 //                          //
0121 //     CLASS DEFINITION     //
0122 //                          //
0123 //////////////////////////////
0124 
0125 /////////////////////////////////////
0126 // this class is needed to make a map
0127 // between different types of stubs
0128 struct L1TStubCompare {
0129 public:
0130   bool operator()(const trklet::L1TStub& a, const trklet::L1TStub& b) const {
0131     if (a.x() != b.x())
0132       return (b.x() > a.x());
0133     else if (a.y() != b.y())
0134       return (b.y() > a.y());
0135     else if (a.z() != b.z())
0136       return (a.z() > b.z());
0137     else
0138       return a.bend() > b.bend();
0139   }
0140 };
0141 
0142 class L1FPGATrackProducer : public edm::one::EDProducer<edm::one::WatchRuns> {
0143 public:
0144   /// Constructor/destructor
0145   explicit L1FPGATrackProducer(const edm::ParameterSet& iConfig);
0146   ~L1FPGATrackProducer() override;
0147 
0148 private:
0149   int eventnum;
0150 
0151   /// Containers of parameters passed by python configuration file
0152   edm::ParameterSet config;
0153 
0154   bool readMoreMcTruth_;
0155 
0156   /// File path for configuration files
0157   edm::FileInPath fitPatternFile;
0158   edm::FileInPath memoryModulesFile;
0159   edm::FileInPath processingModulesFile;
0160   edm::FileInPath wiresFile;
0161 
0162   edm::FileInPath tableTEDFile;
0163   edm::FileInPath tableTREFile;
0164 
0165   string asciiEventOutName_;
0166   std::ofstream asciiEventOut_;
0167 
0168   // settings containing various constants for the tracklet processing
0169   trklet::Settings settings_;
0170 
0171   // event processor for the tracklet track finding
0172   trklet::TrackletEventProcessor eventProcessor;
0173 
0174   // used to "kill" stubs from a selected area of the detector
0175   StubKiller* stubKiller_;
0176   int failScenario_;
0177 
0178   unsigned int nHelixPar_;
0179   bool extended_;
0180   bool reduced_;
0181 
0182   bool trackQuality_;
0183   std::unique_ptr<L1TrackQuality> trackQualityModel_;
0184 
0185   std::map<string, vector<int>> dtclayerdisk;
0186 
0187   edm::InputTag MCTruthClusterInputTag;
0188   edm::InputTag MCTruthStubInputTag;
0189   edm::InputTag TrackingParticleInputTag;
0190 
0191   const edm::EDGetTokenT<reco::BeamSpot> getTokenBS_;
0192   const edm::EDGetTokenT<TTDTC> getTokenDTC_;
0193   edm::EDGetTokenT<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>> getTokenTTClusterMCTruth_;
0194   edm::EDGetTokenT<std::vector<TrackingParticle>> getTokenTrackingParticle_;
0195 
0196   // ED output token for clock and bit accurate tracks
0197   const edm::EDPutTokenT<Streams> putTokenTracks_;
0198   // ED output token for clock and bit accurate stubs
0199   const edm::EDPutTokenT<StreamsStub> putTokenStubs_;
0200   // ChannelAssignment token
0201   const ESGetToken<ChannelAssignment, ChannelAssignmentRcd> esGetTokenChannelAssignment_;
0202   // helper class to assign tracks to channel
0203   const ChannelAssignment* channelAssignment_;
0204 
0205   // helper class to store DTC configuration
0206   const Setup* setup_;
0207   // helper class to store configuration needed by HitPatternHelper
0208   const hph::Setup* setupHPH_;
0209 
0210   // Setup token
0211   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> esGetTokenBfield_;
0212   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> esGetTokenTGeom_;
0213   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> esGetTokenTTopo_;
0214   const edm::ESGetToken<tt::Setup, tt::SetupRcd> esGetToken_;
0215   const edm::ESGetToken<hph::Setup, hph::SetupRcd> esGetTokenHPH_;
0216 
0217   /// ///////////////// ///
0218   /// MANDATORY METHODS ///
0219   void beginRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
0220   void endRun(edm::Run const&, edm::EventSetup const&) override;
0221   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0222 };
0223 
0224 //////////////
0225 // CONSTRUCTOR
0226 L1FPGATrackProducer::L1FPGATrackProducer(edm::ParameterSet const& iConfig)
0227     : config(iConfig),
0228       readMoreMcTruth_(iConfig.getParameter<bool>("readMoreMcTruth")),
0229       MCTruthClusterInputTag(readMoreMcTruth_ ? config.getParameter<edm::InputTag>("MCTruthClusterInputTag")
0230                                               : edm::InputTag()),
0231       MCTruthStubInputTag(readMoreMcTruth_ ? config.getParameter<edm::InputTag>("MCTruthStubInputTag")
0232                                            : edm::InputTag()),
0233       TrackingParticleInputTag(readMoreMcTruth_ ? iConfig.getParameter<edm::InputTag>("TrackingParticleInputTag")
0234                                                 : edm::InputTag()),
0235       // book ED products
0236       getTokenBS_(consumes<reco::BeamSpot>(config.getParameter<edm::InputTag>("BeamSpotSource"))),
0237       getTokenDTC_(consumes<TTDTC>(edm::InputTag(iConfig.getParameter<edm::InputTag>("InputTagTTDTC")))),
0238       // book ED output token for clock and bit accurate tracks
0239       putTokenTracks_(produces<Streams>("Level1TTTracks")),
0240       // book ED output token for clock and bit accurate stubs
0241       putTokenStubs_(produces<StreamsStub>("Level1TTTracks")),
0242       // book ES products
0243       esGetTokenChannelAssignment_(esConsumes<ChannelAssignment, ChannelAssignmentRcd, Transition::BeginRun>()),
0244       esGetTokenBfield_(esConsumes<edm::Transition::BeginRun>()),
0245       esGetTokenTGeom_(esConsumes()),
0246       esGetTokenTTopo_(esConsumes()),
0247       esGetToken_(esConsumes<tt::Setup, tt::SetupRcd, edm::Transition::BeginRun>()),
0248       esGetTokenHPH_(esConsumes<hph::Setup, hph::SetupRcd, edm::Transition::BeginRun>()) {
0249   if (readMoreMcTruth_) {
0250     getTokenTTClusterMCTruth_ = consumes<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>>(MCTruthClusterInputTag);
0251     getTokenTrackingParticle_ = consumes<std::vector<TrackingParticle>>(TrackingParticleInputTag);
0252   }
0253 
0254   produces<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>("Level1TTTracks").setBranchAlias("Level1TTTracks");
0255 
0256   asciiEventOutName_ = iConfig.getUntrackedParameter<string>("asciiFileName", "");
0257 
0258   fitPatternFile = iConfig.getParameter<edm::FileInPath>("fitPatternFile");
0259   processingModulesFile = iConfig.getParameter<edm::FileInPath>("processingModulesFile");
0260   memoryModulesFile = iConfig.getParameter<edm::FileInPath>("memoryModulesFile");
0261   wiresFile = iConfig.getParameter<edm::FileInPath>("wiresFile");
0262 
0263   failScenario_ = iConfig.getUntrackedParameter<int>("FailScenario", 0);
0264 
0265   extended_ = iConfig.getParameter<bool>("Extended");
0266   reduced_ = iConfig.getParameter<bool>("Reduced");
0267   nHelixPar_ = iConfig.getParameter<unsigned int>("Hnpar");
0268 
0269   if (extended_) {
0270     tableTEDFile = iConfig.getParameter<edm::FileInPath>("tableTEDFile");
0271     tableTREFile = iConfig.getParameter<edm::FileInPath>("tableTREFile");
0272   }
0273 
0274   // initial ES products
0275   channelAssignment_ = nullptr;
0276   setup_ = nullptr;
0277 
0278   // --------------------------------------------------------------------------------
0279   // set options in Settings based on inputs from configuration files
0280   // --------------------------------------------------------------------------------
0281 
0282   settings_.setExtended(extended_);
0283   settings_.setReduced(reduced_);
0284   settings_.setNHelixPar(nHelixPar_);
0285 
0286   settings_.setFitPatternFile(fitPatternFile.fullPath());
0287   settings_.setProcessingModulesFile(processingModulesFile.fullPath());
0288   settings_.setMemoryModulesFile(memoryModulesFile.fullPath());
0289   settings_.setWiresFile(wiresFile.fullPath());
0290 
0291   settings_.setFakefit(iConfig.getParameter<bool>("Fakefit"));
0292   settings_.setStoreTrackBuilderOutput(iConfig.getParameter<bool>("StoreTrackBuilderOutput"));
0293   settings_.setRemovalType(iConfig.getParameter<string>("RemovalType"));
0294   settings_.setDoMultipleMatches(iConfig.getParameter<bool>("DoMultipleMatches"));
0295 
0296   if (extended_) {
0297     settings_.setTableTEDFile(tableTEDFile.fullPath());
0298     settings_.setTableTREFile(tableTREFile.fullPath());
0299 
0300     //FIXME: The TED and TRE tables are currently disabled by default, so we
0301     //need to allow for the additional tracklets that will eventually be
0302     //removed by these tables, once they are finalized
0303     settings_.setNbitstrackletindex(15);
0304   }
0305 
0306   eventnum = 0;
0307   if (not asciiEventOutName_.empty()) {
0308     asciiEventOut_.open(asciiEventOutName_.c_str());
0309   }
0310 
0311   if (settings_.debugTracklet()) {
0312     edm::LogVerbatim("Tracklet") << "fit pattern :     " << fitPatternFile.fullPath()
0313                                  << "\n process modules : " << processingModulesFile.fullPath()
0314                                  << "\n memory modules :  " << memoryModulesFile.fullPath()
0315                                  << "\n wires          :  " << wiresFile.fullPath();
0316     if (extended_) {
0317       edm::LogVerbatim("Tracklet") << "table_TED    :  " << tableTEDFile.fullPath()
0318                                    << "\n table_TRE    :  " << tableTREFile.fullPath();
0319     }
0320   }
0321 
0322   trackQuality_ = iConfig.getParameter<bool>("TrackQuality");
0323   if (trackQuality_) {
0324     trackQualityModel_ = std::make_unique<L1TrackQuality>(iConfig.getParameter<edm::ParameterSet>("TrackQualityPSet"));
0325   }
0326   if (settings_.storeTrackBuilderOutput() && (settings_.doMultipleMatches() || !settings_.removalType().empty())) {
0327     cms::Exception exception("ConfigurationNotSupported.");
0328     exception.addContext("L1FPGATrackProducer::produce");
0329     if (settings_.doMultipleMatches())
0330       exception << "Storing of TrackBuilder output does not support doMultipleMatches.";
0331     if (!settings_.removalType().empty())
0332       exception << "Storing of TrackBuilder output does not support duplicate removal.";
0333     throw exception;
0334   }
0335 }
0336 
0337 /////////////
0338 // DESTRUCTOR
0339 L1FPGATrackProducer::~L1FPGATrackProducer() {
0340   if (asciiEventOut_.is_open()) {
0341     asciiEventOut_.close();
0342   }
0343 }
0344 
0345 ///////END RUN
0346 //
0347 void L1FPGATrackProducer::endRun(const edm::Run& run, const edm::EventSetup& iSetup) {}
0348 
0349 ////////////
0350 // BEGIN JOB
0351 void L1FPGATrackProducer::beginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0352   ////////////////////////
0353   // GET MAGNETIC FIELD //
0354   const MagneticField* theMagneticField = &iSetup.getData(esGetTokenBfield_);
0355   double mMagneticFieldStrength = theMagneticField->inTesla(GlobalPoint(0, 0, 0)).z();
0356   settings_.setBfield(mMagneticFieldStrength);
0357 
0358   setup_ = &iSetup.getData(esGetToken_);
0359 
0360   settings_.passSetup(setup_);
0361 
0362   setupHPH_ = &iSetup.getData(esGetTokenHPH_);
0363   // Tracklet pattern reco output channel info.
0364   channelAssignment_ = &iSetup.getData(esGetTokenChannelAssignment_);
0365   // initialize the tracklet event processing (this sets all the processing & memory modules, wiring, etc)
0366   eventProcessor.init(settings_, setup_);
0367 }
0368 
0369 //////////
0370 // PRODUCE
0371 void L1FPGATrackProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0372   typedef std::map<trklet::L1TStub,
0373                    edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>,
0374                    L1TStubCompare>
0375       stubMapType;
0376   typedef std::map<unsigned int,
0377                    edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
0378       stubIndexMapType;
0379   typedef edm::Ref<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_>>, TTCluster<Ref_Phase2TrackerDigi_>>
0380       TTClusterRef;
0381 
0382   /// Prepare output
0383   auto L1TkTracksForOutput = std::make_unique<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>();
0384 
0385   stubMapType stubMap;
0386   stubIndexMapType stubIndexMap;
0387 
0388   ////////////
0389   // GET BS //
0390   edm::Handle<reco::BeamSpot> beamSpotHandle;
0391   iEvent.getByToken(getTokenBS_, beamSpotHandle);
0392   math::XYZPoint bsPosition = beamSpotHandle->position();
0393 
0394   eventnum++;
0395   trklet::SLHCEvent ev;
0396   ev.setEventNum(eventnum);
0397   ev.setIP(bsPosition.x(), bsPosition.y());
0398 
0399   // tracking particles
0400   edm::Handle<std::vector<TrackingParticle>> TrackingParticleHandle;
0401   if (readMoreMcTruth_)
0402     iEvent.getByToken(getTokenTrackingParticle_, TrackingParticleHandle);
0403 
0404   // tracker topology
0405   const TrackerTopology* const tTopo = &iSetup.getData(esGetTokenTTopo_);
0406   const TrackerGeometry* const theTrackerGeom = &iSetup.getData(esGetTokenTGeom_);
0407 
0408   // check killing stubs for detector degradation studies
0409   // if failType = 0, StubKiller does not kill any modules
0410   int failType = 0;
0411   if (failScenario_ < 0 || failScenario_ > 9) {
0412     edm::LogVerbatim("Tracklet") << "Invalid fail scenario! Ignoring input";
0413   } else
0414     failType = failScenario_;
0415 
0416   stubKiller_ = new StubKiller();
0417   stubKiller_->initialise(failType, tTopo, theTrackerGeom);
0418 
0419   ////////////////////////
0420   // GET THE PRIMITIVES //
0421   edm::Handle<TTDTC> handleDTC;
0422   iEvent.getByToken<TTDTC>(getTokenDTC_, handleDTC);
0423 
0424   // must be defined for code to compile, even if it's not used unless readMoreMcTruth_ is true
0425   map<edm::Ptr<TrackingParticle>, int> translateTP;
0426 
0427   // MC truth association maps
0428   edm::Handle<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>> MCTruthTTClusterHandle;
0429   if (readMoreMcTruth_) {
0430     iEvent.getByToken(getTokenTTClusterMCTruth_, MCTruthTTClusterHandle);
0431 
0432     ////////////////////////////////////////////////
0433     /// LOOP OVER TRACKING PARTICLES & GET SIMTRACKS
0434 
0435     int ntps = 1;  //count from 1 ; 0 will mean invalid
0436 
0437     int this_tp = 0;
0438     if (readMoreMcTruth_) {
0439       for (const auto& iterTP : *TrackingParticleHandle) {
0440         edm::Ptr<TrackingParticle> tp_ptr(TrackingParticleHandle, this_tp);
0441         this_tp++;
0442 
0443         // only keep TPs producing a cluster
0444         if (MCTruthTTClusterHandle->findTTClusterRefs(tp_ptr).empty())
0445           continue;
0446 
0447         if (iterTP.g4Tracks().empty()) {
0448           continue;
0449         }
0450 
0451         int sim_eventid = iterTP.g4Tracks().at(0).eventId().event();
0452         int sim_type = iterTP.pdgId();
0453         float sim_pt = iterTP.pt();
0454         float sim_eta = iterTP.eta();
0455         float sim_phi = iterTP.phi();
0456 
0457         float vx = iterTP.vertex().x();
0458         float vy = iterTP.vertex().y();
0459         float vz = iterTP.vertex().z();
0460 
0461         if (sim_pt < 1.0 || std::abs(vz) > 100.0 || hypot(vx, vy) > 50.0)
0462           continue;
0463 
0464         ev.addL1SimTrack(sim_eventid, ntps, sim_type, sim_pt, sim_eta, sim_phi, vx, vy, vz);
0465 
0466         translateTP[tp_ptr] = ntps;
0467         ntps++;
0468 
0469       }  //end loop over TPs
0470     }
0471 
0472   }  // end if (readMoreMcTruth_)
0473 
0474   /////////////////////////////////
0475   /// READ DTC STUB INFORMATION ///
0476   /////////////////////////////////
0477 
0478   // Process stubs in each region and channel within that tracking region
0479   unsigned int theStubIndex = 0;
0480   for (const int& region : handleDTC->tfpRegions()) {
0481     for (const int& channel : handleDTC->tfpChannels()) {
0482       // Get the DTC name & ID from the channel
0483       unsigned int atcaSlot = channel % 12;
0484       string dtcname = settings_.slotToDTCname(atcaSlot);
0485       if (channel % 24 >= 12)
0486         dtcname = "neg" + dtcname;
0487       dtcname += (channel < 24) ? "_A" : "_B";  // which detector region
0488       int dtcId = setup_->dtcId(region, channel);
0489 
0490       // Get the stubs from the DTC
0491       const tt::StreamStub& streamFromDTC{handleDTC->stream(region, channel)};
0492 
0493       // Prepare the DTC stubs for the IR
0494       for (size_t stubIndex = 0; stubIndex < streamFromDTC.size(); ++stubIndex) {
0495         const tt::FrameStub& stub{streamFromDTC[stubIndex]};
0496         const TTStubRef& stubRef = stub.first;
0497 
0498         if (stubRef.isNull())
0499           continue;
0500 
0501         const GlobalPoint& ttPos = setup_->stubPos(stubRef);
0502 
0503         //Get the 2 bits for the layercode
0504         string layerword = stub.second.to_string().substr(61, 2);
0505         unsigned int layercode = 2 * (layerword[0] - '0') + layerword[1] - '0';
0506         assert(layercode < 4);
0507 
0508         //translation from the two bit layercode to the layer/disk number of each of the
0509         //12 channels (dtcs)
0510         // FIX: take this from DTC cabling map.
0511         static const int layerdisktab[12][4] = {{0, 6, 8, 10},
0512                                                 {0, 7, 9, -1},
0513                                                 {1, 7, -1, -1},
0514                                                 {6, 8, 10, -1},
0515                                                 {2, 7, -1, -1},
0516                                                 {2, 9, -1, -1},
0517                                                 {3, 4, -1, -1},
0518                                                 {4, -1, -1, -1},
0519                                                 {5, -1, -1, -1},
0520                                                 {5, 8, -1, -1},
0521                                                 {6, 9, -1, -1},
0522                                                 {7, 10, -1, -1}};
0523 
0524         int layerdisk = layerdisktab[channel % 12][layercode];
0525         assert(layerdisk != -1);
0526 
0527         //Get the 36 bit word - skip the lowest 3 buts (status and layer code)
0528         constexpr int DTCLinkWordSize = 64;
0529         constexpr int StubWordSize = 36;
0530         constexpr int LayerandStatusCodeSize = 3;
0531         string stubword =
0532             stub.second.to_string().substr(DTCLinkWordSize - StubWordSize - LayerandStatusCodeSize, StubWordSize);
0533         string stubwordhex = "";
0534 
0535         //Loop over the 9 words in the 36 bit stub word
0536         for (unsigned int i = 0; i < 9; i++) {
0537           bitset<4> bits(stubword.substr(i * 4, 4));
0538           ulong val = bits.to_ulong();
0539           stubwordhex += ((val < 10) ? ('0' + val) : ('A' + val - 10));
0540         }
0541 
0542         /// Get the Inner and Outer TTCluster
0543         edm::Ref<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_>>, TTCluster<Ref_Phase2TrackerDigi_>>
0544             innerCluster = stub.first->clusterRef(0);
0545         edm::Ref<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_>>, TTCluster<Ref_Phase2TrackerDigi_>>
0546             outerCluster = stub.first->clusterRef(1);
0547 
0548         // -----------------------------------------------------
0549         // check module orientation, if flipped, need to store that information for track fit
0550         // -----------------------------------------------------
0551 
0552         const DetId innerDetId = innerCluster->getDetId();
0553         const GeomDetUnit* det_inner = theTrackerGeom->idToDetUnit(innerDetId);
0554         const auto* theGeomDet_inner = dynamic_cast<const PixelGeomDetUnit*>(det_inner);
0555         const PixelTopology* topol_inner = dynamic_cast<const PixelTopology*>(&(theGeomDet_inner->specificTopology()));
0556 
0557         MeasurementPoint coords_inner = innerCluster->findAverageLocalCoordinatesCentered();
0558         LocalPoint clustlp_inner = topol_inner->localPosition(coords_inner);
0559         GlobalPoint posStub_inner = theGeomDet_inner->surface().toGlobal(clustlp_inner);
0560 
0561         const DetId outerDetId = outerCluster->getDetId();
0562         const GeomDetUnit* det_outer = theTrackerGeom->idToDetUnit(outerDetId);
0563         const auto* theGeomDet_outer = dynamic_cast<const PixelGeomDetUnit*>(det_outer);
0564         const PixelTopology* topol_outer = dynamic_cast<const PixelTopology*>(&(theGeomDet_outer->specificTopology()));
0565 
0566         MeasurementPoint coords_outer = outerCluster->findAverageLocalCoordinatesCentered();
0567         LocalPoint clustlp_outer = topol_outer->localPosition(coords_outer);
0568         GlobalPoint posStub_outer = theGeomDet_outer->surface().toGlobal(clustlp_outer);
0569 
0570         bool isFlipped = (posStub_outer.mag() < posStub_inner.mag());
0571 
0572         vector<int> assocTPs;
0573 
0574         for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over both clusters that make up stub.
0575 
0576           const TTClusterRef& ttClusterRef = stubRef->clusterRef(iClus);
0577 
0578           // Now identify all TP's contributing to either cluster in stub.
0579           if (readMoreMcTruth_) {
0580             vector<edm::Ptr<TrackingParticle>> vecTpPtr =
0581                 MCTruthTTClusterHandle->findTrackingParticlePtrs(ttClusterRef);
0582 
0583             for (const edm::Ptr<TrackingParticle>& tpPtr : vecTpPtr) {
0584               if (translateTP.find(tpPtr) != translateTP.end()) {
0585                 if (iClus == 0) {
0586                   assocTPs.push_back(translateTP.at(tpPtr));
0587                 } else {
0588                   assocTPs.push_back(-translateTP.at(tpPtr));
0589                 }
0590                 // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0591               } else {
0592                 assocTPs.push_back(0);
0593               }
0594             }
0595           }
0596         }
0597 
0598         double stubbend = stubRef->bendFE();  //stubRef->rawBend()
0599         if (ttPos.z() < -120) {
0600           stubbend = -stubbend;
0601         }
0602 
0603         bool barrel = (layerdisk < N_LAYER);
0604         // See  https://github.com/cms-sw/cmssw/tree/master/Geometry/TrackerNumberingBuilder
0605         enum TypeBarrel { nonBarrel = 0, tiltedMinus = 1, tiltedPlus = 2, flat = 3 };
0606         const TypeBarrel type = static_cast<TypeBarrel>(tTopo->tobSide(innerDetId));
0607         bool tiltedBarrel = barrel && (type == tiltedMinus || type == tiltedPlus);
0608         unsigned int tiltedRingId = 0;
0609         // Tilted module ring no. (Increasing 1 to 12 as |z| increases).
0610         if (tiltedBarrel) {
0611           tiltedRingId = tTopo->tobRod(innerDetId);
0612           if (type == tiltedMinus) {
0613             unsigned int layp1 = 1 + layerdisk;  // Setup counts from 1
0614             unsigned int nTilted = setup_->numTiltedLayerRing(layp1);
0615             tiltedRingId = 1 + nTilted - tiltedRingId;
0616           }
0617         }
0618         // Endcap module ring number (1-15) in endcap disks.
0619         unsigned int endcapRingId = barrel ? 0 : tTopo->tidRing(innerDetId);
0620 
0621         const unsigned int intDetId = innerDetId.rawId();
0622 
0623         // check killing stubs for detector degredation studies
0624         const TTStub<Ref_Phase2TrackerDigi_>* theStub = &(*stubRef);
0625         bool killThisStub = stubKiller_->killStub(theStub);
0626         if (!killThisStub) {
0627           ev.addStub(dtcname,
0628                      region,
0629                      layerdisk,
0630                      stubwordhex,
0631                      setup_->psModule(dtcId),
0632                      isFlipped,
0633                      tiltedBarrel,
0634                      tiltedRingId,
0635                      endcapRingId,
0636                      intDetId,
0637                      ttPos.x(),
0638                      ttPos.y(),
0639                      ttPos.z(),
0640                      stubbend,
0641                      stubRef->innerClusterPosition(),
0642                      assocTPs,
0643                      theStubIndex);
0644 
0645           const trklet::L1TStub& lastStub = ev.lastStub();
0646           stubMap[lastStub] = stubRef;
0647           stubIndexMap[lastStub.uniqueIndex()] = stub.first;
0648           theStubIndex++;
0649         }
0650       }
0651     }
0652   }
0653 
0654   //////////////////////////
0655   // NOW RUN THE L1 tracking
0656 
0657   if (!asciiEventOutName_.empty()) {
0658     ev.write(asciiEventOut_);
0659   }
0660 
0661   const std::vector<trklet::Track>& tracks = eventProcessor.tracks();
0662 
0663   // max number of projection layers
0664   const unsigned int maxNumProjectionLayers = channelAssignment_->maxNumProjectionLayers();
0665   // number of track channels
0666   const unsigned int numStreamsTrack = N_SECTOR * channelAssignment_->numChannelsTrack();
0667   // number of stub channels
0668   const unsigned int numStreamsStub = N_SECTOR * channelAssignment_->numChannelsStub();
0669   // number of seeding layers
0670   const unsigned int numSeedingLayers = channelAssignment_->numSeedingLayers();
0671   // max number of stub channel per track
0672   const unsigned int numStubChannel = maxNumProjectionLayers + numSeedingLayers;
0673   // number of stub channels if all seed types streams padded to have same number of stub channels (for coding simplicity)
0674   const unsigned int numStreamsStubRaw = numStreamsTrack * numStubChannel;
0675 
0676   // Streams formatted to allow this code to run outside CMSSW.
0677   vector<vector<string>> streamsTrackRaw(numStreamsTrack);
0678   vector<vector<StubStreamData>> streamsStubRaw(numStreamsStubRaw);
0679 
0680   // this performs the actual tracklet event processing
0681   eventProcessor.event(ev, streamsTrackRaw, streamsStubRaw);
0682 
0683   for (const auto& track : tracks) {
0684     if (track.duplicate())
0685       continue;
0686 
0687     // this is where we create the TTTrack object
0688     double tmp_rinv = track.rinv(settings_);
0689     double tmp_phi = track.phi0(settings_);
0690     double tmp_tanL = track.tanL(settings_);
0691     double tmp_z0 = track.z0(settings_);
0692     double tmp_d0 = track.d0(settings_);
0693     double tmp_chi2rphi = track.chisqrphi();
0694     double tmp_chi2rz = track.chisqrz();
0695     unsigned int tmp_hit = track.hitpattern();
0696 
0697     TTTrack<Ref_Phase2TrackerDigi_> aTrack(tmp_rinv,
0698                                            tmp_phi,
0699                                            tmp_tanL,
0700                                            tmp_z0,
0701                                            tmp_d0,
0702                                            tmp_chi2rphi,
0703                                            tmp_chi2rz,
0704                                            0,
0705                                            0,
0706                                            0,
0707                                            tmp_hit,
0708                                            settings_.nHelixPar(),
0709                                            settings_.bfield());
0710 
0711     unsigned int trksector = track.sector();
0712     unsigned int trkseed = (unsigned int)abs(track.seed());
0713 
0714     aTrack.setPhiSector(trksector);
0715     aTrack.setTrackSeedType(trkseed);
0716 
0717     const vector<trklet::L1TStub>& stubptrs = track.stubs();
0718     vector<trklet::L1TStub> stubs;
0719 
0720     stubs.reserve(stubptrs.size());
0721     for (const auto& stubptr : stubptrs) {
0722       stubs.push_back(stubptr);
0723     }
0724 
0725     int countStubs = 0;
0726     stubMapType::const_iterator it;
0727     stubIndexMapType::const_iterator itIndex;
0728     for (const auto& itstubs : stubs) {
0729       itIndex = stubIndexMap.find(itstubs.uniqueIndex());
0730       if (itIndex != stubIndexMap.end()) {
0731         aTrack.addStubRef(itIndex->second);
0732         countStubs = countStubs + 1;
0733       } else {
0734         // could not find stub in stub map
0735       }
0736     }
0737 
0738     // pt consistency
0739     aTrack.setStubPtConsistency(
0740         StubPtConsistency::getConsistency(aTrack, theTrackerGeom, tTopo, settings_.bfield(), settings_.nHelixPar()));
0741 
0742     if (trackQuality_) {
0743       trackQualityModel_->setL1TrackQuality(aTrack);
0744     }
0745 
0746     //    hph::HitPatternHelper hph(setupHPH_, tmp_hit, tmp_tanL, tmp_z0);
0747     //    if (trackQuality_) {
0748     //      trackQualityModel_->setBonusFeatures(hph.bonusFeatures());
0749     //    }
0750 
0751     // set TTTrack word
0752     aTrack.setTrackWordBits();
0753     // test track word
0754     //aTrack.testTrackWordBits();
0755 
0756     L1TkTracksForOutput->push_back(aTrack);
0757   }
0758 
0759   iEvent.put(std::move(L1TkTracksForOutput), "Level1TTTracks");
0760 
0761   // produce clock and bit accurate stream output tracks and stubs.
0762   // from end of tracklet pattern recognition.
0763   // Convertion here is from stream format that allows this code to run
0764   // outside CMSSW to the EDProduct one.
0765   Streams streamsTrack(numStreamsTrack);
0766   StreamsStub streamsStub(numStreamsStub);
0767 
0768   for (unsigned int chanTrk = 0; chanTrk < numStreamsTrack; chanTrk++) {
0769     for (unsigned int itk = 0; itk < streamsTrackRaw[chanTrk].size(); itk++) {
0770       std::string bitsTrk = streamsTrackRaw[chanTrk][itk];
0771       int iSeed = chanTrk % channelAssignment_->numChannelsTrack();  // seed type
0772       streamsTrack[chanTrk].emplace_back(bitsTrk);
0773 
0774       const unsigned int chanStubOffsetIn = chanTrk * numStubChannel;
0775       const unsigned int chanStubOffsetOut = channelAssignment_->offsetStub(chanTrk);
0776       const unsigned int numProjLayers = channelAssignment_->numProjectionLayers(iSeed);
0777       TTBV hitMap(0, numProjLayers + numSeedingLayers);
0778       // remove padding from stub stream
0779       for (unsigned int iproj = 0; iproj < numStubChannel; iproj++) {
0780         // FW current has one (perhaps invalid) stub per layer per track.
0781         const StubStreamData& stubdata = streamsStubRaw[chanStubOffsetIn + iproj][itk];
0782         const L1TStub& stub = stubdata.stub();
0783         if (!stubdata.valid())
0784           continue;
0785         const TTStubRef& ttStubRef = stubMap[stub];
0786         const int seedType = stubdata.iSeed();
0787         const int layerId = setup_->layerId(ttStubRef);
0788         const int channelId = channelAssignment_->channelId(seedType, layerId);
0789         hitMap.set(channelId);
0790         streamsStub[chanStubOffsetOut + channelId].emplace_back(ttStubRef, stubdata.dataBits());
0791       }
0792       for (int layerId : hitMap.ids(false)) {  // invalid stubs
0793         streamsStub[chanStubOffsetOut + layerId].emplace_back(tt::FrameStub());
0794       }
0795     }
0796   }
0797 
0798   iEvent.emplace(putTokenTracks_, std::move(streamsTrack));
0799   iEvent.emplace(putTokenStubs_, std::move(streamsStub));
0800 
0801 }  /// End of produce()
0802 
0803 // ///////////////////////////
0804 // // DEFINE THIS AS A PLUG-IN
0805 DEFINE_FWK_MODULE(L1FPGATrackProducer);