L1TTrackMatch/plugins/L1FastTrackingJetProducer.cc

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 // Producer of L1FastTrackingJets
0004 //
0005 // FastTracking Jets: Jets created by running the FastJet clustering algorithm on L1 tracks that have been matched to a Tracking Particle
0006 // Author: G.Karathanasis , CU Boulder
0007 ///////////////////////////////////////////////////////////////////////////
0008
0009 // system include files
0010 #include <memory>
0011
0012 // user include files
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/stream/EDProducer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0022 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0023 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0024 #include "DataFormats/Math/interface/LorentzVector.h"
0025
0026 // L1 objects
0027 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0028 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0029 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0030 #include "DataFormats/L1Trigger/interface/Vertex.h"
0031
0032 // geometry
0033 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0034 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0035
0036 //mc
0037 #include "SimDataFormats/Associations/interface/TTTrackAssociationMap.h"
0038 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0039
0040 #include <fastjet/JetDefinition.hh>
0041
0042 #include <string>
0043 #include "TMath.h"
0044 #include "TH1.h"
0045
0046 using namespace l1t;
0047 using namespace edm;
0048 using namespace std;
0049
0050 //////////////////////////////
0051 //                          //
0052 //     CLASS DEFINITION     //
0053 //                          //
0054 //////////////////////////////
0055
0056 class L1FastTrackingJetProducer : public edm::stream::EDProducer<> {
0057 public:
0058   typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
0059   typedef std::vector<L1TTTrackType> L1TTTrackCollectionType;
0060
0061   explicit L1FastTrackingJetProducer(const edm::ParameterSet&);
0062   ~L1FastTrackingJetProducer() override;
0063
0064   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0065
0066 private:
0067   virtual void beginJob();
0068   void produce(edm::Event&, const edm::EventSetup&) override;
0069   virtual void endJob();
0070
0071   // track selection criteria
0072   const float trkZMax_;          // in [cm]
0073   const float trkChi2dofMax_;    // maximum track chi2dof
0074   const double trkBendChi2Max_;  // maximum track bendchi2
0075   const float trkPtMin_;         // in [GeV]
0076   const float trkEtaMax_;        // in [rad]
0077   const int trkNStubMin_;        // minimum number of stubs
0078   const int trkNPSStubMin_;      // minimum number of PS stubs
0079   const double deltaZ0Cut_;      // save with |L1z-z0| < maxZ0
0080   const double coneSize_;        // Use anti-kt with this cone size
0081   const bool doTightChi2_;
0082   const float trkPtTightChi2_;
0083   const float trkChi2dofTightChi2_;
0084   const bool displaced_;  //use prompt/displaced tracks
0085   bool selectTrkMatchGenTight_;
0086   bool selectTrkMatchGenLoose_;
0087   bool selectTrkMatchGenOrPU_;
0088
0089   const edm::EDGetTokenT<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>> trackToken_;
0090   edm::EDGetTokenT<std::vector<l1t::Vertex>> pvToken_;
0091   const edm::EDGetTokenT<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>> genToken_;
0092   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0093   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tGeomToken_;
0094 };
0095
0096 // constructor
0097 L1FastTrackingJetProducer::L1FastTrackingJetProducer(const edm::ParameterSet& iConfig)
0098     : trkZMax_((float)iConfig.getParameter<double>("trk_zMax")),
0099       trkChi2dofMax_((float)iConfig.getParameter<double>("trk_chi2dofMax")),
0100       trkBendChi2Max_(iConfig.getParameter<double>("trk_bendChi2Max")),
0101       trkPtMin_((float)iConfig.getParameter<double>("trk_ptMin")),
0102       trkEtaMax_((float)iConfig.getParameter<double>("trk_etaMax")),
0103       trkNStubMin_((int)iConfig.getParameter<int>("trk_nStubMin")),
0104       trkNPSStubMin_((int)iConfig.getParameter<int>("trk_nPSStubMin")),
0105       deltaZ0Cut_((float)iConfig.getParameter<double>("deltaZ0Cut")),
0106       coneSize_((float)iConfig.getParameter<double>("coneSize")),
0107       doTightChi2_(iConfig.getParameter<bool>("doTightChi2")),
0108       trkPtTightChi2_((float)iConfig.getParameter<double>("trk_ptTightChi2")),
0109       trkChi2dofTightChi2_((float)iConfig.getParameter<double>("trk_chi2dofTightChi2")),
0110       displaced_(iConfig.getParameter<bool>("displaced")),
0111       selectTrkMatchGenTight_(iConfig.getParameter<bool>("selectTrkMatchGenTight")),
0112       selectTrkMatchGenLoose_(iConfig.getParameter<bool>("selectTrkMatchGenLoose")),
0113       selectTrkMatchGenOrPU_(iConfig.getParameter<bool>("selectTrkMatchGenOrPU")),
0114       trackToken_(consumes<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>(
0115           iConfig.getParameter<edm::InputTag>("L1TrackInputTag"))),
0116       pvToken_(consumes<std::vector<l1t::Vertex>>(iConfig.getParameter<edm::InputTag>("L1PrimaryVertexTag"))),
0117       genToken_(
0118           consumes<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>>(iConfig.getParameter<edm::InputTag>("GenInfo"))),
0119       tTopoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>(edm::ESInputTag("", ""))),
0120       tGeomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>(edm::ESInputTag("", ""))) {
0121   if (displaced_)
0122     produces<TkJetCollection>("L1FastTrackingJetsExtended");
0123   else
0124     produces<TkJetCollection>("L1FastTrackingJets");
0125 }
0126
0127 // destructor
0128 L1FastTrackingJetProducer::~L1FastTrackingJetProducer() {}
0129
0130 // producer
0131 void L1FastTrackingJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0132   std::unique_ptr<TkJetCollection> L1FastTrackingJets(new TkJetCollection);
0133
0134   // L1 tracks
0135   edm::Handle<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>> TTTrackHandle;
0136   iEvent.getByToken(trackToken_, TTTrackHandle);
0137   std::vector<TTTrack<Ref_Phase2TrackerDigi_>>::const_iterator iterL1Track;
0138
0139   // Gen
0140   edm::Handle<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>> MCTrkAssociation;
0141   if (!iEvent.isRealData())
0142     iEvent.getByToken(genToken_, MCTrkAssociation);
0143
0144   // Tracker Topology
0145   const TrackerTopology& tTopo = iSetup.getData(tTopoToken_);
0146   //const TrackerGeometry& tGeom = iSetup.getData(tGeomToken_);
0147
0148   edm::Handle<std::vector<l1t::Vertex>> L1VertexHandle;
0149   iEvent.getByToken(pvToken_, L1VertexHandle);
0150
0151   fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, coneSize_);
0152   std::vector<fastjet::PseudoJet> JetInputs;
0153
0154   float recoVtx = L1VertexHandle->begin()->z0();
0155   unsigned int this_l1track = 0;
0156   for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) {
0157     this_l1track++;
0158     float trk_pt = iterL1Track->momentum().perp();
0159     float trk_z0 = iterL1Track->z0();
0160     float trk_chi2dof = iterL1Track->chi2Red();
0161     float trk_bendchi2 = iterL1Track->stubPtConsistency();
0162     std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
0163         theStubs = iterL1Track->getStubRefs();
0164     int trk_nstub = (int)theStubs.size();
0165
0166     if (std::abs(trk_z0) > trkZMax_)
0167       continue;
0168     if (std::abs(iterL1Track->momentum().eta()) > trkEtaMax_)
0169       continue;
0170     if (trk_pt < trkPtMin_)
0171       continue;
0172     if (trk_nstub < trkNStubMin_)
0173       continue;
0174     if (trk_chi2dof > trkChi2dofMax_)
0175       continue;
0176     if (trk_bendchi2 > trkBendChi2Max_)
0177       continue;
0178     if (doTightChi2_ && (trk_pt > trkPtTightChi2_ && trk_chi2dof > trkChi2dofTightChi2_))
0179       continue;
0180
0181     int trk_nPS = 0;
0182     for (int istub = 0; istub < trk_nstub; istub++) {
0183       DetId detId(theStubs.at(istub)->getDetId());
0184       bool tmp_isPS = false;
0185       if (detId.det() == DetId::Detector::Tracker) {
0186         if (detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3)
0187           tmp_isPS = true;
0188         else if (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9)
0189           tmp_isPS = true;
0190       }
0191       if (tmp_isPS)
0192         trk_nPS++;
0193     }
0194     if (trk_nPS < trkNPSStubMin_)
0195       continue;
0196     if (std::abs(recoVtx - trk_z0) > deltaZ0Cut_)
0197       continue;
0198     if (!iEvent.isRealData()) {
0199       edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_>> trk_ptr(TTTrackHandle, this_l1track);
0200
0201       if (!(MCTrkAssociation->isGenuine(trk_ptr)) && selectTrkMatchGenTight_)
0202         continue;
0203       if (!(MCTrkAssociation->isLooselyGenuine(trk_ptr) || MCTrkAssociation->isGenuine(trk_ptr)) &&
0204           selectTrkMatchGenLoose_)
0205         continue;
0206       if (!(MCTrkAssociation->isLooselyGenuine(trk_ptr) || MCTrkAssociation->isGenuine(trk_ptr) ||
0207             MCTrkAssociation->isCombinatoric(trk_ptr)) &&
0208           selectTrkMatchGenOrPU_)
0209         continue;
0210     }
0211
0212     fastjet::PseudoJet psuedoJet(iterL1Track->momentum().x(),
0213                                  iterL1Track->momentum().y(),
0214                                  iterL1Track->momentum().z(),
0215                                  iterL1Track->momentum().mag());
0216     JetInputs.push_back(psuedoJet);                     // input tracks for clustering
0217     JetInputs.back().set_user_index(this_l1track - 1);  // save track index in the collection
0218   }  // end loop over tracks
0219
0220   fastjet::ClusterSequence cs(JetInputs, jet_def);  // define the output jet collection
0221   std::vector<fastjet::PseudoJet> JetOutputs =
0222       fastjet::sorted_by_pt(cs.inclusive_jets(0));  // output jet collection, pT-ordered
0223
0224   for (unsigned int ijet = 0; ijet < JetOutputs.size(); ++ijet) {
0225     math::XYZTLorentzVector jetP4(
0226         JetOutputs[ijet].px(), JetOutputs[ijet].py(), JetOutputs[ijet].pz(), JetOutputs[ijet].modp());
0227     float sumpt = 0;
0228     float avgZ = 0;
0229     std::vector<edm::Ptr<L1TTTrackType>> L1TrackPtrs;
0230     std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(cs.constituents(JetOutputs[ijet]));
0231
0232     for (unsigned int i = 0; i < fjConstituents.size(); ++i) {
0233       auto index = fjConstituents[i].user_index();
0234       edm::Ptr<L1TTTrackType> trkPtr(TTTrackHandle, index);
0235       L1TrackPtrs.push_back(trkPtr);  // L1Tracks in the jet
0236       sumpt = sumpt + trkPtr->momentum().perp();
0237       avgZ = avgZ + trkPtr->momentum().perp() * trkPtr->z0();
0238     }
0239     avgZ = avgZ / sumpt;
0240     edm::Ref<JetBxCollection> jetRef;
0241     TkJet trkJet(jetP4, jetRef, L1TrackPtrs, avgZ);
0242     L1FastTrackingJets->push_back(trkJet);
0243   }  //end loop over Jet Outputs
0244
0245   if (displaced_)
0246     iEvent.put(std::move(L1FastTrackingJets), "L1FastTrackingJetsExtended");
0247   else
0248     iEvent.put(std::move(L1FastTrackingJets), "L1FastTrackingJets");
0249 }
0250
0251 void L1FastTrackingJetProducer::beginJob() {}
0252
0253 void L1FastTrackingJetProducer::endJob() {}
0254
0255 void L1FastTrackingJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0256   {
0257     // L1FastTrackingJets
0258     edm::ParameterSetDescription desc;
0259     desc.add<edm::InputTag>("L1TrackInputTag", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks"));
0260     desc.add<std::string>("L1PrimaryVertexTag", "L1Vertices");
0261     desc.add<edm::InputTag>("GenInfo", edm::InputTag("TTTrackAssociatorFromPixelDigis", "Level1TTTracks"));
0262     desc.add<double>("trk_zMax", 15.0);
0263     desc.add<double>("trk_chi2dofMax", 10.0);
0264     desc.add<double>("trk_bendChi2Max", 2.2);
0265     desc.add<double>("trk_ptMin", 2.0);
0266     desc.add<double>("trk_etaMax", 2.5);
0267     desc.add<int>("trk_nStubMin", 4);
0268     desc.add<int>("trk_nPSStubMin", -1);
0269     desc.add<double>("deltaZ0Cut", 0.5);
0270     desc.add<bool>("doTightChi2", true);
0271     desc.add<double>("trk_ptTightChi2", 20.0);
0272     desc.add<double>("trk_chi2dofTightChi2", 5.0);
0273     desc.add<double>("coneSize", 0.4);
0274     desc.add<bool>("displaced", false);
0275     desc.add<bool>("selectTrkMatchGenTight", true);
0276     desc.add<bool>("selectTrkMatchGenLoose", false);
0277     desc.add<bool>("selectTrkMatchGenOrPU", false);
0278     descriptions.add("L1FastTrackingJets", desc);
0279     // or use the following to generate the label from the module's C++ type
0280     //descriptions.addWithDefaultLabel(desc);
0281   }
0282   /*{
0283     // L1FastTrackingJetsExtended
0284     desc.add<double>("trk_bendChi2Max", 2.4);
0285     desc.add<double>("trk_ptMin", 3.0);
0286     desc.add<double>("trk_etaMax", 2.5);
0287     desc.add<int>("trk_nStubMin", 4);
0288     desc.add<int>("trk_nPSStubMin", -1);
0289     desc.add<double>("deltaZ0Cut", 3.0);
0290     desc.add<bool>("doTightChi2", true);
0291     desc.add<double>("trk_ptTightChi2", 20.0);
0292     desc.add<double>("trk_chi2dofTightChi2", 5.0);
0293     desc.add<double>("coneSize", 0.4);
0294     desc.add<bool>("displaced", true);
0295     desc.add<bool>("selectTrkMatchGenTight", true);
0296     desc.add<bool>("selectTrkMatchGenLoose", false);
0297     desc.add<bool>("selectTrkMatchGenOrPU", false);
0298     descriptions.add("L1FastTrackingJetsExtended", desc);
0299     // or use the following to generate the label from the module's C++ type
0300     //descriptions.addWithDefaultLabel(desc);
0301   }*/
0302 }
0303
0304 //define this as a plug-in
0305 DEFINE_FWK_MODULE(L1FastTrackingJetProducer);