Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-14 01:44:02

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