Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-06-03 00:12:21

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