L1TTrackMatch/plugins/L1TruthTrackFastJetProducer.cc

0001 ///////////////////////////////////////////////////////////////////////////
0002 //                                                                       //
0003 // Producer of TruTrkFastJet,                                            //
0004 // Cluster L1 tracks with truth info using fastjet                       //
0005 //                                                                       //
0006 // Created by: Claire Savard (Oct. 2023)                                 //
0007 //                                                                       //
0008 ///////////////////////////////////////////////////////////////////////////
0009
0010 // system include files
0011 #include <memory>
0012
0013 // user include files
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0023 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0024 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0025 #include "DataFormats/Math/interface/LorentzVector.h"
0026
0027 // L1 objects
0028 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0029 #include "DataFormats/L1TCorrelator/interface/TkJet.h"
0030 #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h"
0031 #include "DataFormats/L1Trigger/interface/Vertex.h"
0032
0033 // MC
0034 #include "SimDataFormats/Associations/interface/TTTrackAssociationMap.h"
0035 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0036
0037 // geometry
0038 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0039 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0040
0041 #include <fastjet/JetDefinition.hh>
0042
0043 #include <string>
0044 #include "TMath.h"
0045 #include "TH1.h"
0046
0047 using namespace l1t;
0048 using namespace edm;
0049 using namespace std;
0050
0051 //////////////////////////////
0052 //                          //
0053 //     CLASS DEFINITION     //
0054 //                          //
0055 //////////////////////////////
0056
0057 class L1TruthTrackFastJetProducer : public edm::stream::EDProducer<> {
0058 public:
0059   typedef TTTrack<Ref_Phase2TrackerDigi_> L1TTTrackType;
0060   typedef std::vector<L1TTTrackType> L1TTTrackCollectionType;
0061
0062   explicit L1TruthTrackFastJetProducer(const edm::ParameterSet&);
0063   ~L1TruthTrackFastJetProducer() override = default;
0064
0065   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0066
0067 private:
0068   void produce(edm::Event&, const edm::EventSetup&) override;
0069
0070   // track selection criteria
0071   const float trkZMax_;      // in [cm]
0072   const float trkPtMin_;     // in [GeV]
0073   const float trkEtaMax_;    // in [rad]
0074   const int trkNStubMin_;    // minimum number of stubs
0075   const int trkNPSStubMin_;  // minimum number of PS stubs
0076   const double coneSize_;    // Use anti-kt with this cone size
0077   const bool displaced_;     //use prompt/displaced tracks
0078
0079   const edm::EDGetTokenT<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > > trackToken_;
0080   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0081   const edm::EDGetTokenT<TTTrackAssociationMap<Ref_Phase2TrackerDigi_> > ttTrackMCTruthToken_;
0082 };
0083
0084 // constructor
0085 L1TruthTrackFastJetProducer::L1TruthTrackFastJetProducer(const edm::ParameterSet& iConfig)
0086     : trkZMax_((float)iConfig.getParameter<double>("trk_zMax")),
0087       trkPtMin_((float)iConfig.getParameter<double>("trk_ptMin")),
0088       trkEtaMax_((float)iConfig.getParameter<double>("trk_etaMax")),
0089       trkNStubMin_((int)iConfig.getParameter<int>("trk_nStubMin")),
0090       trkNPSStubMin_((int)iConfig.getParameter<int>("trk_nPSStubMin")),
0091       coneSize_((float)iConfig.getParameter<double>("coneSize")),
0092       displaced_(iConfig.getParameter<bool>("displaced")),
0093       trackToken_(consumes<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > >(
0094           iConfig.getParameter<edm::InputTag>("L1TrackInputTag"))),
0095       tTopoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>(edm::ESInputTag("", ""))),
0096       ttTrackMCTruthToken_(consumes<TTTrackAssociationMap<Ref_Phase2TrackerDigi_> >(
0097           iConfig.getParameter<edm::InputTag>("MCTruthTrackInputTag"))) {
0098   if (displaced_)
0099     produces<TkJetCollection>("L1TruthTrackFastJetsExtended");
0100   else
0101     produces<TkJetCollection>("L1TruthTrackFastJets");
0102 }
0103
0104 // producer
0105 void L1TruthTrackFastJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0106   std::unique_ptr<TkJetCollection> L1TrackFastJets(new TkJetCollection);
0107
0108   // L1 tracks
0109   edm::Handle<std::vector<TTTrack<Ref_Phase2TrackerDigi_> > > TTTrackHandle;
0110   iEvent.getByToken(trackToken_, TTTrackHandle);
0111   std::vector<TTTrack<Ref_Phase2TrackerDigi_> >::const_iterator iterL1Track;
0112
0113   // MC truth
0114   edm::Handle<TTTrackAssociationMap<Ref_Phase2TrackerDigi_> > MCTruthTTTrackHandle;
0115   iEvent.getByToken(ttTrackMCTruthToken_, MCTruthTTTrackHandle);
0116
0117   // Tracker Topology
0118   const TrackerTopology& tTopo = iSetup.getData(tTopoToken_);
0119
0120   fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, coneSize_);
0121   std::vector<fastjet::PseudoJet> JetInputs;
0122
0123   unsigned int this_l1track = 0;
0124   for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) {
0125     edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_> > l1track_ptr(TTTrackHandle, this_l1track);
0126     this_l1track++;
0127     std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_> >, TTStub<Ref_Phase2TrackerDigi_> > >
0128         theStubs = iterL1Track->getStubRefs();
0129
0130     // standard quality cuts
0131     if (std::abs(iterL1Track->z0()) > trkZMax_)
0132       continue;
0133     if (std::abs(iterL1Track->momentum().eta()) > trkEtaMax_)
0134       continue;
0135     if (iterL1Track->momentum().perp() < trkPtMin_)
0136       continue;
0137     int trk_nstub = (int)theStubs.size();
0138     if (trk_nstub < trkNStubMin_)
0139       continue;
0140
0141     int trk_nPS = 0;
0142     for (int istub = 0; istub < trk_nstub; istub++) {
0143       DetId detId(theStubs.at(istub)->getDetId());
0144       bool tmp_isPS = false;
0145       if (detId.det() == DetId::Detector::Tracker) {
0146         if (detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3)
0147           tmp_isPS = true;
0148         else if (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9)
0149           tmp_isPS = true;
0150       }
0151       if (tmp_isPS)
0152         trk_nPS++;
0153     }
0154     if (trk_nPS < trkNPSStubMin_)
0155       continue;
0156
0157     // check that trk is real and from hard interaction
0158     edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(l1track_ptr);
0159     if (my_tp.isNull())  // there is no tp match so the track is fake
0160       continue;
0161     int tp_eventid = my_tp->eventId().event();
0162     if (tp_eventid > 0)  // matched tp is from pileup
0163       continue;
0164
0165     fastjet::PseudoJet psuedoJet(iterL1Track->momentum().x(),
0166                                  iterL1Track->momentum().y(),
0167                                  iterL1Track->momentum().z(),
0168                                  iterL1Track->momentum().mag());
0169     JetInputs.push_back(psuedoJet);                     // input tracks for clustering
0170     JetInputs.back().set_user_index(this_l1track - 1);  // save track index in the collection
0171   }  // end loop over tracks
0172
0173   fastjet::ClusterSequence cs(JetInputs, jet_def);  // define the output jet collection
0174   std::vector<fastjet::PseudoJet> JetOutputs =
0175       fastjet::sorted_by_pt(cs.inclusive_jets(0));  // output jet collection, pT-ordered
0176
0177   for (unsigned int ijet = 0; ijet < JetOutputs.size(); ++ijet) {
0178     math::XYZTLorentzVector jetP4(
0179         JetOutputs[ijet].px(), JetOutputs[ijet].py(), JetOutputs[ijet].pz(), JetOutputs[ijet].modp());
0180     float sumpt = 0;
0181     float avgZ = 0;
0182     std::vector<edm::Ptr<L1TTTrackType> > L1TrackPtrs;
0183     std::vector<fastjet::PseudoJet> fjConstituents = fastjet::sorted_by_pt(cs.constituents(JetOutputs[ijet]));
0184
0185     for (unsigned int i = 0; i < fjConstituents.size(); ++i) {
0186       auto index = fjConstituents[i].user_index();
0187       edm::Ptr<L1TTTrackType> trkPtr(TTTrackHandle, index);
0188       L1TrackPtrs.push_back(trkPtr);  // L1Tracks in the jet
0189       sumpt = sumpt + trkPtr->momentum().perp();
0190       avgZ = avgZ + trkPtr->momentum().perp() * trkPtr->z0();
0191     }
0192     avgZ = avgZ / sumpt;
0193     edm::Ref<JetBxCollection> jetRef;
0194     TkJet trkJet(jetP4, jetRef, L1TrackPtrs, avgZ);
0195     L1TrackFastJets->push_back(trkJet);
0196   }  //end loop over Jet Outputs
0197
0198   if (displaced_)
0199     iEvent.put(std::move(L1TrackFastJets), "L1TruthTrackFastJetsExtended");
0200   else
0201     iEvent.put(std::move(L1TrackFastJets), "L1TruthTrackFastJets");
0202 }
0203
0204 void L1TruthTrackFastJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0205   // The following says we do not know what parameters are allowed so do no validation
0206   edm::ParameterSetDescription desc;
0207   desc.add<edm::InputTag>("L1TrackInputTag", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks"));
0208   desc.add<edm::InputTag>("MCTruthTrackInputTag", edm::InputTag("TTTrackAssociatorFromPixelDigis", "Level1TTTracks"));
0209   desc.add<double>("trk_zMax", 15.);
0210   desc.add<double>("trk_ptMin", 2.0);
0211   desc.add<double>("trk_etaMax", 2.4);
0212   desc.add<int>("trk_nStubMin", 4);
0213   desc.add<int>("trk_nPSStubMin", -1);
0214   desc.add<double>("coneSize", 0.4);
0215   desc.add<bool>("displaced", false);
0216   descriptions.add("l1tTruthTrackFastJets", desc);
0217 }
0218
0219 //define this as a plug-in
0220 DEFINE_FWK_MODULE(L1TruthTrackFastJetProducer);