Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-15 23:40:44

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