TrackFindingTracklet/test/L1TrackNtupleMaker.cc

0001 //////////////////////////////////////////////////////////////////////
0002 //                                                                  //
0003 //  Analyzer for making mini-ntuple for L1 track performance plots  //
0004 //                                                                  //
0005 //////////////////////////////////////////////////////////////////////
0006
0007 ////////////////////
0008 // FRAMEWORK HEADERS
0009 #include "FWCore/PluginManager/interface/ModuleDef.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/Run.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/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020
0021 ///////////////////////
0022 // DATA FORMATS HEADERS
0023 #include "DataFormats/Common/interface/Handle.h"
0024 #include "DataFormats/Common/interface/Ref.h"
0025
0026 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0027 #include "DataFormats/L1TrackTrigger/interface/TTCluster.h"
0028 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
0029 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
0030 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0031 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0032 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0033 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0034 #include "SimDataFormats/Associations/interface/TTClusterAssociationMap.h"
0035 #include "SimDataFormats/Associations/interface/TTStubAssociationMap.h"
0036 #include "SimDataFormats/Associations/interface/TTTrackAssociationMap.h"
0037 #include "Geometry/Records/interface/StackedTrackerGeometryRecord.h"
0038
0039 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0040 #include "DataFormats/JetReco/interface/GenJet.h"
0041
0042 ////////////////////////////
0043 // DETECTOR GEOMETRY HEADERS
0044 #include "MagneticField/Engine/interface/MagneticField.h"
0045 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0046 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0047 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0048 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0049 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0050
0051 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0052 #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h"
0053 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyBuilder.h"
0054 #include "Geometry/Records/interface/StackedTrackerGeometryRecord.h"
0055
0056 ////////////////
0057 // PHYSICS TOOLS
0058 #include "L1Trigger/TrackTrigger/interface/Setup.h"
0059 #include "L1Trigger/TrackerTFP/interface/LayerEncoding.h"
0060 #include "L1Trigger/TrackFindingTracklet/interface/HitPatternHelper.h"
0061 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0062 #include "CLHEP/Units/PhysicalConstants.h"
0063
0064 ///////////////
0065 // ROOT HEADERS
0066 #include <TROOT.h>
0067 #include <TCanvas.h>
0068 #include <TTree.h>
0069 #include <TFile.h>
0070 #include <TMath.h>
0071 #include <TF1.h>
0072 #include <TH2F.h>
0073 #include <TH1F.h>
0074
0075 //////////////
0076 // STD HEADERS
0077 #include <memory>
0078 #include <string>
0079 #include <iostream>
0080
0081 //////////////
0082 // NAMESPACES
0083 using namespace std;
0084 using namespace edm;
0085
0086 //////////////////////////////
0087 //                          //
0088 //     CLASS DEFINITION     //
0089 //                          //
0090 //////////////////////////////
0091
0092 class L1TrackNtupleMaker : public one::EDAnalyzer<one::WatchRuns, one::SharedResources> {
0093 public:
0094   // Constructor/destructor
0095   explicit L1TrackNtupleMaker(const edm::ParameterSet& iConfig);
0096   ~L1TrackNtupleMaker() override;
0097
0098   // Mandatory methods
0099   void beginJob() override;
0100   void endJob() override;
0101   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0102   void beginRun(const Run& iEvent, const EventSetup& iSetup) override {}
0103   void endRun(const Run& iEvent, const EventSetup& iSetup) override {}
0104
0105 protected:
0106 private:
0107   //-----------------------------------------------------------------------------------------------
0108   // Containers of parameters passed by python configuration file
0109   edm::ParameterSet config;
0110
0111   int MyProcess;       // 11/13/211 for single electrons/muons/pions, 6/15 for pions from ttbar/taus, 1 for inclusive
0112   bool DebugMode;      // lots of debug printout statements
0113   bool SaveAllTracks;  // store in ntuples not only truth-matched tracks but ALL tracks
0114   bool SaveStubs;      // option to save also stubs in the ntuples (makes them large...)
0115   int L1Tk_nPar;       // use 4 or 5 parameter track fit?
0116   int TP_minNStub;  // require TPs to have >= minNStub (defining efficiency denominator) (==0 means to only require >= 1 cluster)
0117   int TP_minNStubLayer;  // require TPs to have stubs in >= minNStubLayer layers/disks (defining efficiency denominator)
0118   double TP_minPt;       // save TPs with pt > minPt
0119   double TP_maxEta;      // save TPs with |eta| < maxEta
0120   double TP_maxZ0;       // save TPs with |z0| < maxZ0
0121   int L1Tk_minNStub;     // require L1 tracks to have >= minNStub (this is mostly for tracklet purposes)
0122
0123   bool TrackingInJets;  // do tracking in jets?
0124
0125   edm::InputTag L1TrackInputTag;       // L1 track collection
0126   edm::InputTag MCTruthTrackInputTag;  // MC truth collection
0127   edm::InputTag MCTruthClusterInputTag;
0128   edm::InputTag L1StubInputTag;
0129   edm::InputTag MCTruthStubInputTag;
0130   edm::InputTag TrackingParticleInputTag;
0131   edm::InputTag TrackingVertexInputTag;
0132   edm::InputTag GenJetInputTag;
0133
0134   edm::EDGetTokenT<edmNew::DetSetVector<TTCluster<Ref_Phase2TrackerDigi_>>> ttClusterToken_;
0135   edm::EDGetTokenT<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>> ttStubToken_;
0136   edm::EDGetTokenT<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>> ttClusterMCTruthToken_;
0137   edm::EDGetTokenT<TTStubAssociationMap<Ref_Phase2TrackerDigi_>> ttStubMCTruthToken_;
0138
0139   edm::EDGetTokenT<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>> ttTrackToken_;
0140   edm::EDGetTokenT<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>> ttTrackMCTruthToken_;
0141
0142   edm::EDGetTokenT<std::vector<TrackingParticle>> TrackingParticleToken_;
0143   edm::EDGetTokenT<std::vector<TrackingVertex>> TrackingVertexToken_;
0144
0145   edm::EDGetTokenT<std::vector<reco::GenJet>> GenJetToken_;
0146
0147   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> getTokenTrackerGeom_;
0148   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> getTokenTrackerTopo_;
0149   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> getTokenBField_;
0150   edm::ESGetToken<hph::Setup, hph::SetupRcd> getTokenHPHSetup_;
0151   edm::ESGetToken<tt::Setup, tt::SetupRcd> getTokenSetup_;
0152   edm::ESGetToken<trackerTFP::LayerEncoding, trackerTFP::DataFormatsRcd> getTokenLayerEncoding_;
0153   //-----------------------------------------------------------------------------------------------
0154   // tree & branches for mini-ntuple
0155
0156   bool available_;  // ROOT file for histograms is open.
0157
0158   TTree* eventTree;
0159
0160   // all L1 tracks
0161   std::vector<float>* m_trk_pt;
0162   std::vector<float>* m_trk_eta;
0163   std::vector<float>* m_trk_phi;
0164   std::vector<float>* m_trk_d0;  // (filled if L1Tk_nPar==5, else 999)
0165   std::vector<float>* m_trk_z0;
0166   std::vector<float>* m_trk_chi2;
0167   std::vector<float>* m_trk_chi2_dof;
0168   std::vector<float>* m_trk_chi2rphi;
0169   std::vector<float>* m_trk_chi2rphi_dof;
0170   std::vector<float>* m_trk_chi2rz;
0171   std::vector<float>* m_trk_chi2rz_dof;
0172   std::vector<float>* m_trk_bendchi2;
0173   std::vector<int>* m_trk_nstub;
0174   std::vector<int>* m_trk_lhits;
0175   std::vector<int>* m_trk_dhits;
0176   std::vector<int>* m_trk_seed;
0177   std::vector<int>* m_trk_hitpattern;
0178   std::vector<int>* m_trk_lhits_hitpattern;  // 6-digit hit mask (barrel layer only) dervied from hitpattern
0179   std::vector<int>* m_trk_dhits_hitpattern;  // disk only
0180   std::vector<int>* m_trk_nPSstub_hitpattern;
0181   std::vector<int>* m_trk_n2Sstub_hitpattern;
0182   std::vector<int>* m_trk_nLostPSstub_hitpattern;
0183   std::vector<int>* m_trk_nLost2Sstub_hitpattern;
0184   std::vector<int>* m_trk_nLoststub_V1_hitpattern;  // Same as the definiton of "nlaymiss_interior" in TrackQuality.cc
0185   std::vector<int>* m_trk_nLoststub_V2_hitpattern;  // A tighter version of "nlaymiss_interior"
0186   std::vector<int>* m_trk_charge;
0187   std::vector<unsigned int>* m_trk_phiSector;
0188   std::vector<int>* m_trk_etaSector;
0189   std::vector<int>* m_trk_genuine;
0190   std::vector<int>* m_trk_loose;
0191   std::vector<int>* m_trk_unknown;
0192   std::vector<int>* m_trk_combinatoric;
0193   std::vector<int>* m_trk_fake;  //0 fake, 1 track from primary interaction, 2 secondary track
0194   std::vector<float>* m_trk_MVA1;
0195   std::vector<int>* m_trk_matchtp_pdgid;
0196   std::vector<float>* m_trk_matchtp_pt;
0197   std::vector<float>* m_trk_matchtp_eta;
0198   std::vector<float>* m_trk_matchtp_phi;
0199   std::vector<float>* m_trk_matchtp_z0;
0200   std::vector<float>* m_trk_matchtp_lxy;
0201   std::vector<float>* m_trk_matchtp_d0;
0202   std::vector<int>* m_trk_injet;          //is the track within dR<0.4 of a genjet with pt > 30 GeV?
0203   std::vector<int>* m_trk_injet_highpt;   //is the track within dR<0.4 of a genjet with pt > 100 GeV?
0204   std::vector<int>* m_trk_injet_vhighpt;  //is the track within dR<0.4 of a genjet with pt > 200 GeV?
0205   std::vector<std::vector<int>>* m_trk_layers;
0206
0207   // all tracking particles
0208   std::vector<float>* m_tp_pt;
0209   std::vector<float>* m_tp_eta;
0210   std::vector<float>* m_tp_phi;
0211   std::vector<float>* m_tp_lxy;
0212   std::vector<float>* m_tp_d0;
0213   std::vector<float>* m_tp_z0;
0214   std::vector<float>* m_tp_d0_prod;
0215   std::vector<float>* m_tp_z0_prod;
0216   std::vector<int>* m_tp_pdgid;
0217   std::vector<int>* m_tp_nmatch;
0218   std::vector<int>* m_tp_nstub;
0219   std::vector<int>* m_tp_eventid;
0220   std::vector<int>* m_tp_charge;
0221   std::vector<int>* m_tp_injet;
0222   std::vector<int>* m_tp_injet_highpt;
0223   std::vector<int>* m_tp_injet_vhighpt;
0224
0225   // *L1 track* properties if m_tp_nmatch > 0
0226   std::vector<float>* m_matchtrk_pt;
0227   std::vector<float>* m_matchtrk_eta;
0228   std::vector<float>* m_matchtrk_phi;
0229   std::vector<float>* m_matchtrk_d0;  //this variable is only filled if L1Tk_nPar==5
0230   std::vector<float>* m_matchtrk_z0;
0231   std::vector<float>* m_matchtrk_chi2;
0232   std::vector<float>* m_matchtrk_chi2_dof;
0233   std::vector<float>* m_matchtrk_chi2rphi;
0234   std::vector<float>* m_matchtrk_chi2rphi_dof;
0235   std::vector<float>* m_matchtrk_chi2rz;
0236   std::vector<float>* m_matchtrk_chi2rz_dof;
0237   std::vector<float>* m_matchtrk_bendchi2;
0238   std::vector<float>* m_matchtrk_MVA1;
0239   std::vector<int>* m_matchtrk_nstub;
0240   std::vector<int>* m_matchtrk_lhits;
0241   std::vector<int>* m_matchtrk_dhits;
0242   std::vector<int>* m_matchtrk_seed;
0243   std::vector<int>* m_matchtrk_hitpattern;
0244   std::vector<int>* m_matchtrk_charge;
0245   std::vector<int>* m_matchtrk_injet;
0246   std::vector<int>* m_matchtrk_injet_highpt;
0247   std::vector<int>* m_matchtrk_injet_vhighpt;
0248
0249   // ALL stubs
0250   std::vector<float>* m_allstub_x;
0251   std::vector<float>* m_allstub_y;
0252   std::vector<float>* m_allstub_z;
0253
0254   std::vector<int>* m_allstub_isBarrel;  // stub is in barrel (1) or in disk (0)
0255   std::vector<int>* m_allstub_layer;
0256   std::vector<int>* m_allstub_isPSmodule;
0257   std::vector<int>* m_allstub_isTiltedBarrel;
0258
0259   std::vector<float>* m_allstub_trigDisplace;
0260   std::vector<float>* m_allstub_trigOffset;
0261   std::vector<float>* m_allstub_trigPos;
0262   std::vector<float>* m_allstub_trigBend;
0263
0264   // stub associated with tracking particle ?
0265   std::vector<int>* m_allstub_matchTP_pdgid;  // -999 if not matched
0266   std::vector<float>* m_allstub_matchTP_pt;   // -999 if not matched
0267   std::vector<float>* m_allstub_matchTP_eta;  // -999 if not matched
0268   std::vector<float>* m_allstub_matchTP_phi;  // -999 if not matched
0269
0270   std::vector<int>* m_allstub_genuine;
0271
0272   // track jet variables (for each gen jet, store the sum of pt of TPs / tracks inside jet cone)
0273   std::vector<float>* m_jet_eta;
0274   std::vector<float>* m_jet_phi;
0275   std::vector<float>* m_jet_pt;
0276   std::vector<float>* m_jet_tp_sumpt;
0277   std::vector<float>* m_jet_trk_sumpt;
0278   std::vector<float>* m_jet_matchtrk_sumpt;
0279 };
0280
0281 //////////////////////////////////
0282 //                              //
0283 //     CLASS IMPLEMENTATION     //
0284 //                              //
0285 //////////////////////////////////
0286
0287 //////////////
0288 // CONSTRUCTOR
0289 L1TrackNtupleMaker::L1TrackNtupleMaker(edm::ParameterSet const& iConfig) : config(iConfig) {
0290   usesResource("TFileService");
0291   MyProcess = iConfig.getParameter<int>("MyProcess");
0292   DebugMode = iConfig.getParameter<bool>("DebugMode");
0293   SaveAllTracks = iConfig.getParameter<bool>("SaveAllTracks");
0294   SaveStubs = iConfig.getParameter<bool>("SaveStubs");
0295   L1Tk_nPar = iConfig.getParameter<int>("L1Tk_nPar");
0296   TP_minNStub = iConfig.getParameter<int>("TP_minNStub");
0297   TP_minNStubLayer = iConfig.getParameter<int>("TP_minNStubLayer");
0298   TP_minPt = iConfig.getParameter<double>("TP_minPt");
0299   TP_maxEta = iConfig.getParameter<double>("TP_maxEta");
0300   TP_maxZ0 = iConfig.getParameter<double>("TP_maxZ0");
0301   L1TrackInputTag = iConfig.getParameter<edm::InputTag>("L1TrackInputTag");
0302   MCTruthTrackInputTag = iConfig.getParameter<edm::InputTag>("MCTruthTrackInputTag");
0303   L1Tk_minNStub = iConfig.getParameter<int>("L1Tk_minNStub");
0304
0305   TrackingInJets = iConfig.getParameter<bool>("TrackingInJets");
0306
0307   L1StubInputTag = iConfig.getParameter<edm::InputTag>("L1StubInputTag");
0308   MCTruthClusterInputTag = iConfig.getParameter<edm::InputTag>("MCTruthClusterInputTag");
0309   MCTruthStubInputTag = iConfig.getParameter<edm::InputTag>("MCTruthStubInputTag");
0310   TrackingParticleInputTag = iConfig.getParameter<edm::InputTag>("TrackingParticleInputTag");
0311   TrackingVertexInputTag = iConfig.getParameter<edm::InputTag>("TrackingVertexInputTag");
0312   GenJetInputTag = iConfig.getParameter<edm::InputTag>("GenJetInputTag");
0313
0314   ttTrackToken_ = consumes<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>>(L1TrackInputTag);
0315   ttTrackMCTruthToken_ = consumes<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>>(MCTruthTrackInputTag);
0316   ttStubToken_ = consumes<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>>(L1StubInputTag);
0317   ttClusterMCTruthToken_ = consumes<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>>(MCTruthClusterInputTag);
0318   ttStubMCTruthToken_ = consumes<TTStubAssociationMap<Ref_Phase2TrackerDigi_>>(MCTruthStubInputTag);
0319
0320   TrackingParticleToken_ = consumes<std::vector<TrackingParticle>>(TrackingParticleInputTag);
0321   TrackingVertexToken_ = consumes<std::vector<TrackingVertex>>(TrackingVertexInputTag);
0322   GenJetToken_ = consumes<std::vector<reco::GenJet>>(GenJetInputTag);
0323
0324   getTokenTrackerGeom_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0325   getTokenTrackerTopo_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0326   getTokenBField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0327   getTokenHPHSetup_ = esConsumes<hph::Setup, hph::SetupRcd>();
0328   getTokenSetup_ = esConsumes<tt::Setup, tt::SetupRcd>();
0329   getTokenLayerEncoding_ = esConsumes<trackerTFP::LayerEncoding, trackerTFP::DataFormatsRcd>();
0330 }
0331
0332 /////////////
0333 // DESTRUCTOR
0334 L1TrackNtupleMaker::~L1TrackNtupleMaker() {}
0335
0336 //////////
0337 // END JOB
0338 void L1TrackNtupleMaker::endJob() {
0339   // things to be done at the exit of the event Loop
0340   edm::LogVerbatim("Tracklet") << "L1TrackNtupleMaker::endJob";
0341
0342   // clean up
0343   delete m_trk_pt;
0344   delete m_trk_eta;
0345   delete m_trk_phi;
0346   delete m_trk_z0;
0347   delete m_trk_d0;
0348   delete m_trk_chi2;
0349   delete m_trk_chi2_dof;
0350   delete m_trk_chi2rphi;
0351   delete m_trk_chi2rphi_dof;
0352   delete m_trk_chi2rz;
0353   delete m_trk_chi2rz_dof;
0354   delete m_trk_bendchi2;
0355   delete m_trk_nstub;
0356   delete m_trk_lhits;
0357   delete m_trk_dhits;
0358   delete m_trk_seed;
0359   delete m_trk_hitpattern;
0360   delete m_trk_lhits_hitpattern;
0361   delete m_trk_dhits_hitpattern;
0362   delete m_trk_nPSstub_hitpattern;
0363   delete m_trk_n2Sstub_hitpattern;
0364   delete m_trk_nLostPSstub_hitpattern;
0365   delete m_trk_nLost2Sstub_hitpattern;
0366   delete m_trk_nLoststub_V1_hitpattern;
0367   delete m_trk_nLoststub_V2_hitpattern;
0368   delete m_trk_charge;
0369   delete m_trk_phiSector;
0370   delete m_trk_etaSector;
0371   delete m_trk_genuine;
0372   delete m_trk_loose;
0373   delete m_trk_unknown;
0374   delete m_trk_combinatoric;
0375   delete m_trk_fake;
0376   delete m_trk_MVA1;
0377   delete m_trk_matchtp_pdgid;
0378   delete m_trk_matchtp_pt;
0379   delete m_trk_matchtp_eta;
0380   delete m_trk_matchtp_phi;
0381   delete m_trk_matchtp_z0;
0382   delete m_trk_matchtp_lxy;
0383   delete m_trk_matchtp_d0;
0384   delete m_trk_injet;
0385   delete m_trk_injet_highpt;
0386   delete m_trk_injet_vhighpt;
0387
0388   delete m_tp_pt;
0389   delete m_tp_eta;
0390   delete m_tp_phi;
0391   delete m_tp_lxy;
0392   delete m_tp_d0;
0393   delete m_tp_z0;
0394   delete m_tp_d0_prod;
0395   delete m_tp_z0_prod;
0396   delete m_tp_pdgid;
0397   delete m_tp_nmatch;
0398   delete m_tp_nstub;
0399   delete m_tp_eventid;
0400   delete m_tp_charge;
0401   delete m_tp_injet;
0402   delete m_tp_injet_highpt;
0403   delete m_tp_injet_vhighpt;
0404
0405   delete m_matchtrk_pt;
0406   delete m_matchtrk_eta;
0407   delete m_matchtrk_phi;
0408   delete m_matchtrk_z0;
0409   delete m_matchtrk_d0;
0410   delete m_matchtrk_chi2;
0411   delete m_matchtrk_chi2_dof;
0412   delete m_matchtrk_chi2rphi;
0413   delete m_matchtrk_chi2rphi_dof;
0414   delete m_matchtrk_chi2rz;
0415   delete m_matchtrk_chi2rz_dof;
0416   delete m_matchtrk_bendchi2;
0417   delete m_matchtrk_MVA1;
0418   delete m_matchtrk_nstub;
0419   delete m_matchtrk_dhits;
0420   delete m_matchtrk_lhits;
0421   delete m_matchtrk_seed;
0422   delete m_matchtrk_hitpattern;
0423   delete m_matchtrk_charge;
0424   delete m_matchtrk_injet;
0425   delete m_matchtrk_injet_highpt;
0426   delete m_matchtrk_injet_vhighpt;
0427
0428   delete m_allstub_x;
0429   delete m_allstub_y;
0430   delete m_allstub_z;
0431   delete m_allstub_isBarrel;
0432   delete m_allstub_layer;
0433   delete m_allstub_isPSmodule;
0434   delete m_allstub_trigDisplace;
0435   delete m_allstub_trigOffset;
0436   delete m_allstub_trigPos;
0437   delete m_allstub_trigBend;
0438   delete m_allstub_matchTP_pdgid;
0439   delete m_allstub_matchTP_pt;
0440   delete m_allstub_matchTP_eta;
0441   delete m_allstub_matchTP_phi;
0442   delete m_allstub_genuine;
0443
0444   delete m_jet_eta;
0445   delete m_jet_phi;
0446   delete m_jet_pt;
0447   delete m_jet_tp_sumpt;
0448   delete m_jet_trk_sumpt;
0449   delete m_jet_matchtrk_sumpt;
0450 }
0451
0452 ////////////
0453 // BEGIN JOB
0454 void L1TrackNtupleMaker::beginJob() {
0455   // things to be done before entering the event Loop
0456   edm::LogVerbatim("Tracklet") << "L1TrackNtupleMaker::beginJob";
0457
0458   //-----------------------------------------------------------------------------------------------
0459   // book histograms / make ntuple
0460   edm::Service<TFileService> fs;
0461   available_ = fs.isAvailable();
0462   if (not available_)
0463     return;  // No ROOT file open.
0464
0465   // initilize
0466   m_trk_pt = new std::vector<float>;
0467   m_trk_eta = new std::vector<float>;
0468   m_trk_phi = new std::vector<float>;
0469   m_trk_z0 = new std::vector<float>;
0470   m_trk_d0 = new std::vector<float>;
0471   m_trk_chi2 = new std::vector<float>;
0472   m_trk_chi2_dof = new std::vector<float>;
0473   m_trk_chi2rphi = new std::vector<float>;
0474   m_trk_chi2rphi_dof = new std::vector<float>;
0475   m_trk_chi2rz = new std::vector<float>;
0476   m_trk_chi2rz_dof = new std::vector<float>;
0477   m_trk_bendchi2 = new std::vector<float>;
0478   m_trk_nstub = new std::vector<int>;
0479   m_trk_lhits = new std::vector<int>;
0480   m_trk_dhits = new std::vector<int>;
0481   m_trk_seed = new std::vector<int>;
0482   m_trk_hitpattern = new std::vector<int>;
0483   m_trk_lhits_hitpattern = new std::vector<int>;
0484   m_trk_dhits_hitpattern = new std::vector<int>;
0485   m_trk_nPSstub_hitpattern = new std::vector<int>;
0486   m_trk_n2Sstub_hitpattern = new std::vector<int>;
0487   m_trk_nLostPSstub_hitpattern = new std::vector<int>;
0488   m_trk_nLost2Sstub_hitpattern = new std::vector<int>;
0489   m_trk_nLoststub_V1_hitpattern = new std::vector<int>;
0490   m_trk_nLoststub_V2_hitpattern = new std::vector<int>;
0491   m_trk_charge = new std::vector<int>;
0492   m_trk_phiSector = new std::vector<unsigned int>;
0493   m_trk_etaSector = new std::vector<int>;
0494   m_trk_genuine = new std::vector<int>;
0495   m_trk_loose = new std::vector<int>;
0496   m_trk_unknown = new std::vector<int>;
0497   m_trk_combinatoric = new std::vector<int>;
0498   m_trk_fake = new std::vector<int>;
0499   m_trk_MVA1 = new std::vector<float>;
0500   m_trk_matchtp_pdgid = new std::vector<int>;
0501   m_trk_matchtp_pt = new std::vector<float>;
0502   m_trk_matchtp_eta = new std::vector<float>;
0503   m_trk_matchtp_phi = new std::vector<float>;
0504   m_trk_matchtp_z0 = new std::vector<float>;
0505   m_trk_matchtp_lxy = new std::vector<float>;
0506   m_trk_matchtp_d0 = new std::vector<float>;
0507   m_trk_injet = new std::vector<int>;
0508   m_trk_injet_highpt = new std::vector<int>;
0509   m_trk_injet_vhighpt = new std::vector<int>;
0510   m_trk_layers = new std::vector<std::vector<int>>;
0511
0512   m_tp_pt = new std::vector<float>;
0513   m_tp_eta = new std::vector<float>;
0514   m_tp_phi = new std::vector<float>;
0515   m_tp_lxy = new std::vector<float>;
0516   m_tp_d0 = new std::vector<float>;
0517   m_tp_z0 = new std::vector<float>;
0518   m_tp_d0_prod = new std::vector<float>;
0519   m_tp_z0_prod = new std::vector<float>;
0520   m_tp_pdgid = new std::vector<int>;
0521   m_tp_nmatch = new std::vector<int>;
0522   m_tp_nstub = new std::vector<int>;
0523   m_tp_eventid = new std::vector<int>;
0524   m_tp_charge = new std::vector<int>;
0525   m_tp_injet = new std::vector<int>;
0526   m_tp_injet_highpt = new std::vector<int>;
0527   m_tp_injet_vhighpt = new std::vector<int>;
0528
0529   m_matchtrk_pt = new std::vector<float>;
0530   m_matchtrk_eta = new std::vector<float>;
0531   m_matchtrk_phi = new std::vector<float>;
0532   m_matchtrk_z0 = new std::vector<float>;
0533   m_matchtrk_d0 = new std::vector<float>;
0534   m_matchtrk_chi2 = new std::vector<float>;
0535   m_matchtrk_chi2_dof = new std::vector<float>;
0536   m_matchtrk_chi2rphi = new std::vector<float>;
0537   m_matchtrk_chi2rphi_dof = new std::vector<float>;
0538   m_matchtrk_chi2rz = new std::vector<float>;
0539   m_matchtrk_chi2rz_dof = new std::vector<float>;
0540   m_matchtrk_bendchi2 = new std::vector<float>;
0541   m_matchtrk_MVA1 = new std::vector<float>;
0542   m_matchtrk_nstub = new std::vector<int>;
0543   m_matchtrk_dhits = new std::vector<int>;
0544   m_matchtrk_lhits = new std::vector<int>;
0545   m_matchtrk_seed = new std::vector<int>;
0546   m_matchtrk_hitpattern = new std::vector<int>;
0547   m_matchtrk_charge = new std::vector<int>;
0548   m_matchtrk_injet = new std::vector<int>;
0549   m_matchtrk_injet_highpt = new std::vector<int>;
0550   m_matchtrk_injet_vhighpt = new std::vector<int>;
0551
0552   m_allstub_x = new std::vector<float>;
0553   m_allstub_y = new std::vector<float>;
0554   m_allstub_z = new std::vector<float>;
0555
0556   m_allstub_isBarrel = new std::vector<int>;
0557   m_allstub_layer = new std::vector<int>;
0558   m_allstub_isPSmodule = new std::vector<int>;
0559   m_allstub_isTiltedBarrel = new std::vector<int>;
0560   m_allstub_trigDisplace = new std::vector<float>;
0561   m_allstub_trigOffset = new std::vector<float>;
0562   m_allstub_trigPos = new std::vector<float>;
0563   m_allstub_trigBend = new std::vector<float>;
0564
0565   m_allstub_matchTP_pdgid = new std::vector<int>;
0566   m_allstub_matchTP_pt = new std::vector<float>;
0567   m_allstub_matchTP_eta = new std::vector<float>;
0568   m_allstub_matchTP_phi = new std::vector<float>;
0569
0570   m_allstub_genuine = new std::vector<int>;
0571
0572   m_jet_eta = new std::vector<float>;
0573   m_jet_phi = new std::vector<float>;
0574   m_jet_pt = new std::vector<float>;
0575   m_jet_tp_sumpt = new std::vector<float>;
0576   m_jet_trk_sumpt = new std::vector<float>;
0577   m_jet_matchtrk_sumpt = new std::vector<float>;
0578
0579   // ntuple
0580   eventTree = fs->make<TTree>("eventTree", "Event tree");
0581
0582   if (SaveAllTracks) {
0583     eventTree->Branch("trk_pt", &m_trk_pt);
0584     eventTree->Branch("trk_eta", &m_trk_eta);
0585     eventTree->Branch("trk_phi", &m_trk_phi);
0586     eventTree->Branch("trk_d0", &m_trk_d0);
0587     eventTree->Branch("trk_z0", &m_trk_z0);
0588     eventTree->Branch("trk_chi2", &m_trk_chi2);
0589     eventTree->Branch("trk_chi2_dof", &m_trk_chi2_dof);
0590     eventTree->Branch("trk_chi2rphi", &m_trk_chi2rphi);
0591     eventTree->Branch("trk_chi2rphi_dof", &m_trk_chi2rphi_dof);
0592     eventTree->Branch("trk_chi2rz", &m_trk_chi2rz);
0593     eventTree->Branch("trk_chi2rz_dof", &m_trk_chi2rz_dof);
0594     eventTree->Branch("trk_bendchi2", &m_trk_bendchi2);
0595     eventTree->Branch("trk_nstub", &m_trk_nstub);
0596     eventTree->Branch("trk_lhits", &m_trk_lhits);
0597     eventTree->Branch("trk_dhits", &m_trk_dhits);
0598     eventTree->Branch("trk_seed", &m_trk_seed);
0599     eventTree->Branch("trk_hitpattern", &m_trk_hitpattern);
0600     eventTree->Branch("trk_lhits_hitpattern", &m_trk_lhits_hitpattern);
0601     eventTree->Branch("trk_dhits_hitpattern", &m_trk_dhits_hitpattern);
0602     eventTree->Branch("trk_nPSstub_hitpattern", &m_trk_nPSstub_hitpattern);
0603     eventTree->Branch("trk_n2Sstub_hitpattern", &m_trk_n2Sstub_hitpattern);
0604     eventTree->Branch("trk_nLostPSstub_hitpattern", &m_trk_nLostPSstub_hitpattern);
0605     eventTree->Branch("trk_nLost2Sstub_hitpattern", &m_trk_nLost2Sstub_hitpattern);
0606     eventTree->Branch("trk_nLoststub_V1_hitpattern", &m_trk_nLoststub_V1_hitpattern);
0607     eventTree->Branch("trk_nLoststub_V2_hitpattern", &m_trk_nLoststub_V2_hitpattern);
0608     eventTree->Branch("trk_charge", &m_trk_charge);
0609     eventTree->Branch("trk_phiSector", &m_trk_phiSector);
0610     eventTree->Branch("trk_etaSector", &m_trk_etaSector);
0611     eventTree->Branch("trk_genuine", &m_trk_genuine);
0612     eventTree->Branch("trk_loose", &m_trk_loose);
0613     eventTree->Branch("trk_unknown", &m_trk_unknown);
0614     eventTree->Branch("trk_combinatoric", &m_trk_combinatoric);
0615     eventTree->Branch("trk_fake", &m_trk_fake);
0616     eventTree->Branch("trk_MVA1", &m_trk_MVA1);
0617     eventTree->Branch("trk_matchtp_pdgid", &m_trk_matchtp_pdgid);
0618     eventTree->Branch("trk_matchtp_pt", &m_trk_matchtp_pt);
0619     eventTree->Branch("trk_matchtp_eta", &m_trk_matchtp_eta);
0620     eventTree->Branch("trk_matchtp_phi", &m_trk_matchtp_phi);
0621     eventTree->Branch("trk_matchtp_z0", &m_trk_matchtp_z0);
0622     eventTree->Branch("trk_matchtp_lxy", &m_trk_matchtp_lxy);
0623     eventTree->Branch("trk_matchtp_d0", &m_trk_matchtp_d0);
0624     if (TrackingInJets) {
0625       eventTree->Branch("trk_injet", &m_trk_injet);
0626       eventTree->Branch("trk_injet_highpt", &m_trk_injet_highpt);
0627       eventTree->Branch("trk_injet_vhighpt", &m_trk_injet_vhighpt);
0628     }
0629     eventTree->Branch("m_trk_layers", &m_trk_layers);
0630   }
0631
0632   eventTree->Branch("tp_pt", &m_tp_pt);
0633   eventTree->Branch("tp_eta", &m_tp_eta);
0634   eventTree->Branch("tp_phi", &m_tp_phi);
0635   eventTree->Branch("tp_lxy", &m_tp_lxy);
0636   eventTree->Branch("tp_d0", &m_tp_d0);
0637   eventTree->Branch("tp_z0", &m_tp_z0);
0638   eventTree->Branch("tp_d0_prod", &m_tp_d0_prod);
0639   eventTree->Branch("tp_z0_prod", &m_tp_z0_prod);
0640   eventTree->Branch("tp_pdgid", &m_tp_pdgid);
0641   eventTree->Branch("tp_nmatch", &m_tp_nmatch);
0642   eventTree->Branch("tp_nstub", &m_tp_nstub);
0643   eventTree->Branch("tp_eventid", &m_tp_eventid);
0644   eventTree->Branch("tp_charge", &m_tp_charge);
0645   if (TrackingInJets) {
0646     eventTree->Branch("tp_injet", &m_tp_injet);
0647     eventTree->Branch("tp_injet_highpt", &m_tp_injet_highpt);
0648     eventTree->Branch("tp_injet_vhighpt", &m_tp_injet_vhighpt);
0649   }
0650
0651   eventTree->Branch("matchtrk_pt", &m_matchtrk_pt);
0652   eventTree->Branch("matchtrk_eta", &m_matchtrk_eta);
0653   eventTree->Branch("matchtrk_phi", &m_matchtrk_phi);
0654   eventTree->Branch("matchtrk_z0", &m_matchtrk_z0);
0655   eventTree->Branch("matchtrk_d0", &m_matchtrk_d0);
0656   eventTree->Branch("matchtrk_chi2", &m_matchtrk_chi2);
0657   eventTree->Branch("matchtrk_chi2_dof", &m_matchtrk_chi2_dof);
0658   eventTree->Branch("matchtrk_chi2rphi", &m_matchtrk_chi2rphi);
0659   eventTree->Branch("matchtrk_chi2rphi_dof", &m_matchtrk_chi2rphi_dof);
0660   eventTree->Branch("matchtrk_chi2rz", &m_matchtrk_chi2rz);
0661   eventTree->Branch("matchtrk_chi2rz_dof", &m_matchtrk_chi2rz_dof);
0662   eventTree->Branch("matchtrk_bendchi2", &m_matchtrk_bendchi2);
0663   eventTree->Branch("matchtrk_MVA1", &m_matchtrk_MVA1);
0664   eventTree->Branch("matchtrk_nstub", &m_matchtrk_nstub);
0665   eventTree->Branch("matchtrk_lhits", &m_matchtrk_lhits);
0666   eventTree->Branch("matchtrk_dhits", &m_matchtrk_dhits);
0667   eventTree->Branch("matchtrk_seed", &m_matchtrk_seed);
0668   eventTree->Branch("matchtrk_hitpattern", &m_matchtrk_hitpattern);
0669   eventTree->Branch("matchtrk_charge", &m_matchtrk_charge);
0670   if (TrackingInJets) {
0671     eventTree->Branch("matchtrk_injet", &m_matchtrk_injet);
0672     eventTree->Branch("matchtrk_injet_highpt", &m_matchtrk_injet_highpt);
0673     eventTree->Branch("matchtrk_injet_vhighpt", &m_matchtrk_injet_vhighpt);
0674   }
0675
0676   if (SaveStubs) {
0677     eventTree->Branch("allstub_x", &m_allstub_x);
0678     eventTree->Branch("allstub_y", &m_allstub_y);
0679     eventTree->Branch("allstub_z", &m_allstub_z);
0680
0681     eventTree->Branch("allstub_isBarrel", &m_allstub_isBarrel);
0682     eventTree->Branch("allstub_layer", &m_allstub_layer);
0683     eventTree->Branch("allstub_isPSmodule", &m_allstub_isPSmodule);
0684     eventTree->Branch("allstub_isTiltedBarrel", &m_allstub_isTiltedBarrel);
0685
0686     eventTree->Branch("allstub_trigDisplace", &m_allstub_trigDisplace);
0687     eventTree->Branch("allstub_trigOffset", &m_allstub_trigOffset);
0688     eventTree->Branch("allstub_trigPos", &m_allstub_trigPos);
0689     eventTree->Branch("allstub_trigBend", &m_allstub_trigBend);
0690
0691     eventTree->Branch("allstub_matchTP_pdgid", &m_allstub_matchTP_pdgid);
0692     eventTree->Branch("allstub_matchTP_pt", &m_allstub_matchTP_pt);
0693     eventTree->Branch("allstub_matchTP_eta", &m_allstub_matchTP_eta);
0694     eventTree->Branch("allstub_matchTP_phi", &m_allstub_matchTP_phi);
0695
0696     eventTree->Branch("allstub_genuine", &m_allstub_genuine);
0697   }
0698
0699   if (TrackingInJets) {
0700     eventTree->Branch("jet_eta", &m_jet_eta);
0701     eventTree->Branch("jet_phi", &m_jet_phi);
0702     eventTree->Branch("jet_pt", &m_jet_pt);
0703     eventTree->Branch("jet_tp_sumpt", &m_jet_tp_sumpt);
0704     eventTree->Branch("jet_trk_sumpt", &m_jet_trk_sumpt);
0705     eventTree->Branch("jet_matchtrk_sumpt", &m_jet_matchtrk_sumpt);
0706   }
0707 }
0708
0709 //////////
0710 // ANALYZE
0711 void L1TrackNtupleMaker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0712   if (not available_)
0713     return;  // No ROOT file open.
0714
0715   if (!(MyProcess == 13 || MyProcess == 11 || MyProcess == 211 || MyProcess == 6 || MyProcess == 15 ||
0716         MyProcess == 1)) {
0717     edm::LogVerbatim("Tracklet") << "The specified MyProcess is invalid! Exiting...";
0718     return;
0719   }
0720
0721   if (!(L1Tk_nPar == 4 || L1Tk_nPar == 5)) {
0722     edm::LogVerbatim("Tracklet") << "Invalid number of track parameters, specified L1Tk_nPar == " << L1Tk_nPar
0723                                  << " but only 4/5 are valid options! Exiting...";
0724     return;
0725   }
0726
0727   // clear variables
0728   if (SaveAllTracks) {
0729     m_trk_pt->clear();
0730     m_trk_eta->clear();
0731     m_trk_phi->clear();
0732     m_trk_d0->clear();
0733     m_trk_z0->clear();
0734     m_trk_chi2->clear();
0735     m_trk_chi2_dof->clear();
0736     m_trk_chi2rphi->clear();
0737     m_trk_chi2rphi_dof->clear();
0738     m_trk_chi2rz->clear();
0739     m_trk_chi2rz_dof->clear();
0740     m_trk_bendchi2->clear();
0741     m_trk_nstub->clear();
0742     m_trk_lhits->clear();
0743     m_trk_dhits->clear();
0744     m_trk_seed->clear();
0745     m_trk_hitpattern->clear();
0746     m_trk_lhits_hitpattern->clear();
0747     m_trk_dhits_hitpattern->clear();
0748     m_trk_nPSstub_hitpattern->clear();
0749     m_trk_n2Sstub_hitpattern->clear();
0750     m_trk_nLostPSstub_hitpattern->clear();
0751     m_trk_nLost2Sstub_hitpattern->clear();
0752     m_trk_nLoststub_V1_hitpattern->clear();
0753     m_trk_nLoststub_V2_hitpattern->clear();
0754     m_trk_charge->clear();
0755     m_trk_phiSector->clear();
0756     m_trk_etaSector->clear();
0757     m_trk_genuine->clear();
0758     m_trk_loose->clear();
0759     m_trk_unknown->clear();
0760     m_trk_combinatoric->clear();
0761     m_trk_fake->clear();
0762     m_trk_MVA1->clear();
0763     m_trk_matchtp_pdgid->clear();
0764     m_trk_matchtp_pt->clear();
0765     m_trk_matchtp_eta->clear();
0766     m_trk_matchtp_phi->clear();
0767     m_trk_matchtp_z0->clear();
0768     m_trk_matchtp_lxy->clear();
0769     m_trk_matchtp_d0->clear();
0770     m_trk_injet->clear();
0771     m_trk_injet_highpt->clear();
0772     m_trk_injet_vhighpt->clear();
0773     m_trk_layers->clear();
0774   }
0775
0776   m_tp_pt->clear();
0777   m_tp_eta->clear();
0778   m_tp_phi->clear();
0779   m_tp_lxy->clear();
0780   m_tp_d0->clear();
0781   m_tp_z0->clear();
0782   m_tp_d0_prod->clear();
0783   m_tp_z0_prod->clear();
0784   m_tp_pdgid->clear();
0785   m_tp_nmatch->clear();
0786   m_tp_nstub->clear();
0787   m_tp_eventid->clear();
0788   m_tp_charge->clear();
0789   m_tp_injet->clear();
0790   m_tp_injet_highpt->clear();
0791   m_tp_injet_vhighpt->clear();
0792
0793   m_matchtrk_pt->clear();
0794   m_matchtrk_eta->clear();
0795   m_matchtrk_phi->clear();
0796   m_matchtrk_z0->clear();
0797   m_matchtrk_d0->clear();
0798   m_matchtrk_chi2->clear();
0799   m_matchtrk_chi2_dof->clear();
0800   m_matchtrk_chi2rphi->clear();
0801   m_matchtrk_chi2rphi_dof->clear();
0802   m_matchtrk_chi2rz->clear();
0803   m_matchtrk_chi2rz_dof->clear();
0804   m_matchtrk_bendchi2->clear();
0805   m_matchtrk_MVA1->clear();
0806   m_matchtrk_nstub->clear();
0807   m_matchtrk_lhits->clear();
0808   m_matchtrk_dhits->clear();
0809   m_matchtrk_seed->clear();
0810   m_matchtrk_hitpattern->clear();
0811   m_matchtrk_charge->clear();
0812   m_matchtrk_injet->clear();
0813   m_matchtrk_injet_highpt->clear();
0814   m_matchtrk_injet_vhighpt->clear();
0815
0816   if (SaveStubs) {
0817     m_allstub_x->clear();
0818     m_allstub_y->clear();
0819     m_allstub_z->clear();
0820
0821     m_allstub_isBarrel->clear();
0822     m_allstub_layer->clear();
0823     m_allstub_isPSmodule->clear();
0824     m_allstub_isTiltedBarrel->clear();
0825
0826     m_allstub_trigDisplace->clear();
0827     m_allstub_trigOffset->clear();
0828     m_allstub_trigPos->clear();
0829     m_allstub_trigBend->clear();
0830
0831     m_allstub_matchTP_pdgid->clear();
0832     m_allstub_matchTP_pt->clear();
0833     m_allstub_matchTP_eta->clear();
0834     m_allstub_matchTP_phi->clear();
0835
0836     m_allstub_genuine->clear();
0837   }
0838
0839   m_jet_eta->clear();
0840   m_jet_phi->clear();
0841   m_jet_pt->clear();
0842   m_jet_tp_sumpt->clear();
0843   m_jet_trk_sumpt->clear();
0844   m_jet_matchtrk_sumpt->clear();
0845
0846   // -----------------------------------------------------------------------------------------------
0847   // retrieve various containers
0848   // -----------------------------------------------------------------------------------------------
0849
0850   // L1 tracks
0851   edm::Handle<std::vector<TTTrack<Ref_Phase2TrackerDigi_>>> TTTrackHandle;
0852   iEvent.getByToken(ttTrackToken_, TTTrackHandle);
0853
0854   // L1 stubs
0855   edm::Handle<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>> TTStubHandle;
0856   if (SaveStubs)
0857     iEvent.getByToken(ttStubToken_, TTStubHandle);
0858
0859   // MC truth association maps
0860   edm::Handle<TTClusterAssociationMap<Ref_Phase2TrackerDigi_>> MCTruthTTClusterHandle;
0861   iEvent.getByToken(ttClusterMCTruthToken_, MCTruthTTClusterHandle);
0862   edm::Handle<TTStubAssociationMap<Ref_Phase2TrackerDigi_>> MCTruthTTStubHandle;
0863   iEvent.getByToken(ttStubMCTruthToken_, MCTruthTTStubHandle);
0864   edm::Handle<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>> MCTruthTTTrackHandle;
0865   iEvent.getByToken(ttTrackMCTruthToken_, MCTruthTTTrackHandle);
0866
0867   // tracking particles
0868   edm::Handle<std::vector<TrackingParticle>> TrackingParticleHandle;
0869   edm::Handle<std::vector<TrackingVertex>> TrackingVertexHandle;
0870   iEvent.getByToken(TrackingParticleToken_, TrackingParticleHandle);
0871   //iEvent.getByToken(TrackingVertexToken_, TrackingVertexHandle);
0872
0873   // -----------------------------------------------------------------------------------------------
0874   // more for TTStubs
0875   edm::ESHandle<TrackerGeometry> tGeomHandle = iSetup.getHandle(getTokenTrackerGeom_);
0876
0877   edm::ESHandle<TrackerTopology> tTopoHandle = iSetup.getHandle(getTokenTrackerTopo_);
0878
0879   edm::ESHandle<MagneticField> bFieldHandle = iSetup.getHandle(getTokenBField_);
0880
0881   edm::ESHandle<hph::Setup> hphHandle = iSetup.getHandle(getTokenHPHSetup_);
0882   edm::ESHandle<tt::Setup> handleSetup = iSetup.getHandle(getTokenSetup_);
0883   edm::ESHandle<trackerTFP::LayerEncoding> handleLayerEncoding = iSetup.getHandle(getTokenLayerEncoding_);
0884
0885   const TrackerTopology* const tTopo = tTopoHandle.product();
0886   const TrackerGeometry* const theTrackerGeom = tGeomHandle.product();
0887   const hph::Setup* hphSetup = hphHandle.product();
0888   const tt::Setup* setup = handleSetup.product();
0889   const trackerTFP::LayerEncoding* layerEncoding = handleLayerEncoding.product();
0890
0891   // ----------------------------------------------------------------------------------------------
0892   // loop over L1 stubs
0893   // ----------------------------------------------------------------------------------------------
0894
0895   if (SaveStubs) {
0896     for (auto gd = theTrackerGeom->dets().begin(); gd != theTrackerGeom->dets().end(); gd++) {
0897       DetId detid = (*gd)->geographicalId();
0898       if (detid.subdetId() != StripSubdetector::TOB && detid.subdetId() != StripSubdetector::TID)
0899         continue;
0900       if (!tTopo->isLower(detid))
0901         continue;                              // loop on the stacks: choose the lower arbitrarily
0902       DetId stackDetid = tTopo->stack(detid);  // Stub module detid
0903
0904       if (TTStubHandle->find(stackDetid) == TTStubHandle->end())
0905         continue;
0906
0907       // Get the DetSets of the Clusters
0908       edmNew::DetSet<TTStub<Ref_Phase2TrackerDigi_>> stubs = (*TTStubHandle)[stackDetid];
0909       const GeomDetUnit* det0 = theTrackerGeom->idToDetUnit(detid);
0910       const auto* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(det0);
0911       const PixelTopology* topol = dynamic_cast<const PixelTopology*>(&(theGeomDet->specificTopology()));
0912
0913       // loop over stubs
0914       for (auto stubIter = stubs.begin(); stubIter != stubs.end(); ++stubIter) {
0915         edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>> tempStubPtr =
0916             edmNew::makeRefTo(TTStubHandle, stubIter);
0917
0918         int isBarrel = 0;
0919         int layer = -999999;
0920         if (detid.subdetId() == StripSubdetector::TOB) {
0921           isBarrel = 1;
0922           layer = static_cast<int>(tTopo->layer(detid));
0923         } else if (detid.subdetId() == StripSubdetector::TID) {
0924           isBarrel = 0;
0925           layer = static_cast<int>(tTopo->layer(detid));
0926         } else {
0927           edm::LogVerbatim("Tracklet") << "WARNING -- neither TOB or TID stub, shouldn't happen...";
0928           layer = -1;
0929         }
0930
0931         int isPSmodule = 0;
0932         if (topol->nrows() == 960)
0933           isPSmodule = 1;
0934
0935         const unsigned int tobSide = tTopo->tobSide(detid);  // nonBarrel = 0, tiltedMinus = 1, tiltedPlus = 2, flat = 3
0936         int isTiltedBarrel = 0;
0937         if (isBarrel == 1 && (tobSide == 1 || tobSide == 2))
0938           isTiltedBarrel = 1;
0939
0940         MeasurementPoint coords = tempStubPtr->clusterRef(0)->findAverageLocalCoordinatesCentered();
0941         LocalPoint clustlp = topol->localPosition(coords);
0942         GlobalPoint posStub = theGeomDet->surface().toGlobal(clustlp);
0943
0944         double tmp_stub_x = posStub.x();
0945         double tmp_stub_y = posStub.y();
0946         double tmp_stub_z = posStub.z();
0947
0948         float trigDisplace = tempStubPtr->rawBend();
0949         float trigOffset = tempStubPtr->bendOffset();
0950         float trigPos = tempStubPtr->innerClusterPosition();
0951         float trigBend = tempStubPtr->bendFE();
0952
0953         m_allstub_x->push_back(tmp_stub_x);
0954         m_allstub_y->push_back(tmp_stub_y);
0955         m_allstub_z->push_back(tmp_stub_z);
0956
0957         m_allstub_isBarrel->push_back(isBarrel);
0958         m_allstub_layer->push_back(layer);
0959         m_allstub_isPSmodule->push_back(isPSmodule);
0960         m_allstub_isTiltedBarrel->push_back(isTiltedBarrel);
0961
0962         m_allstub_trigDisplace->push_back(trigDisplace);
0963         m_allstub_trigOffset->push_back(trigOffset);
0964         m_allstub_trigPos->push_back(trigPos);
0965         m_allstub_trigBend->push_back(trigBend);
0966
0967         // matched to tracking particle?
0968         edm::Ptr<TrackingParticle> my_tp = MCTruthTTStubHandle->findTrackingParticlePtr(tempStubPtr);
0969
0970         int myTP_pdgid = -999;
0971         float myTP_pt = -999;
0972         float myTP_eta = -999;
0973         float myTP_phi = -999;
0974
0975         if (my_tp.isNull() == false) {
0976           int tmp_eventid = my_tp->eventId().event();
0977
0978           if (tmp_eventid > 0)
0979             continue;  // this means stub from pileup track
0980
0981           myTP_pdgid = my_tp->pdgId();
0982           myTP_pt = my_tp->p4().pt();
0983           myTP_eta = my_tp->p4().eta();
0984           myTP_phi = my_tp->p4().phi();
0985         }
0986
0987         m_allstub_matchTP_pdgid->push_back(myTP_pdgid);
0988         m_allstub_matchTP_pt->push_back(myTP_pt);
0989         m_allstub_matchTP_eta->push_back(myTP_eta);
0990         m_allstub_matchTP_phi->push_back(myTP_phi);
0991
0992         int tmp_stub_genuine = 0;
0993         if (MCTruthTTStubHandle->isGenuine(tempStubPtr))
0994           tmp_stub_genuine = 1;
0995
0996         m_allstub_genuine->push_back(tmp_stub_genuine);
0997       }
0998     }
0999   }
1000
1001   // ----------------------------------------------------------------------------------------------
1002   // tracking in jets
1003   // ----------------------------------------------------------------------------------------------
1004
1005   std::vector<math::XYZTLorentzVector> v_jets;
1006   std::vector<int> v_jets_highpt;
1007   std::vector<int> v_jets_vhighpt;
1008
1009   if (TrackingInJets) {
1010     // gen jets
1011     if (DebugMode)
1012       edm::LogVerbatim("Tracklet") << "get genjets";
1013     edm::Handle<std::vector<reco::GenJet>> GenJetHandle;
1014     iEvent.getByToken(GenJetToken_, GenJetHandle);
1015
1016     if (GenJetHandle.isValid()) {
1017       if (DebugMode)
1018         edm::LogVerbatim("Tracklet") << "loop over genjets";
1019       std::vector<reco::GenJet>::const_iterator iterGenJet;
1020       for (iterGenJet = GenJetHandle->begin(); iterGenJet != GenJetHandle->end(); ++iterGenJet) {
1021         reco::GenJet myJet = reco::GenJet(*iterGenJet);
1022
1023         if (myJet.pt() < 30.0)
1024           continue;
1025         if (std::abs(myJet.eta()) > 2.5)
1026           continue;
1027
1028         if (DebugMode)
1029           edm::LogVerbatim("Tracklet") << "genjet pt = " << myJet.pt() << ", eta = " << myJet.eta();
1030
1031         bool ishighpt = false;
1032         bool isveryhighpt = false;
1033         if (myJet.pt() > 100.0)
1034           ishighpt = true;
1035         if (myJet.pt() > 200.0)
1036           isveryhighpt = true;
1037
1038         const math::XYZTLorentzVector& jetP4 = myJet.p4();
1039         v_jets.push_back(jetP4);
1040         if (ishighpt)
1041           v_jets_highpt.push_back(1);
1042         else
1043           v_jets_highpt.push_back(0);
1044         if (isveryhighpt)
1045           v_jets_vhighpt.push_back(1);
1046         else
1047           v_jets_vhighpt.push_back(0);
1048
1049       }  // end loop over genjets
1050     }  // end isValid
1051
1052   }  // end TrackingInJets
1053
1054   const int NJETS = 10;
1055   float jets_tp_sumpt[NJETS] = {0};        //sum pt of TPs with dR<0.4 of jet
1056   float jets_matchtrk_sumpt[NJETS] = {0};  //sum pt of tracks matched to TP with dR<0.4 of jet
1057   float jets_trk_sumpt[NJETS] = {0};       //sum pt of all tracks with dR<0.4 of jet
1058
1059   // ----------------------------------------------------------------------------------------------
1060   // loop over L1 tracks
1061   // ----------------------------------------------------------------------------------------------
1062
1063   if (SaveAllTracks) {
1064     if (DebugMode) {
1065       edm::LogVerbatim("Tracklet") << "\n Loop over L1 tracks!";
1066       edm::LogVerbatim("Tracklet") << "\n Looking at " << L1Tk_nPar << "-parameter tracks!";
1067     }
1068
1069     int this_l1track = 0;
1070     std::vector<TTTrack<Ref_Phase2TrackerDigi_>>::const_iterator iterL1Track;
1071     for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) {
1072       edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_>> l1track_ptr(TTTrackHandle, this_l1track);
1073       this_l1track++;
1074
1075       float tmp_trk_pt = iterL1Track->momentum().perp();
1076       float tmp_trk_eta = iterL1Track->momentum().eta();
1077       float tmp_trk_phi = iterL1Track->momentum().phi();
1078       float tmp_trk_z0 = iterL1Track->z0();  //cm
1079       float tmp_trk_tanL = iterL1Track->tanL();
1080       int tmp_trk_charge = (int)TMath::Sign(1, iterL1Track->rInv());
1081       bool usingNewKF = hphSetup->useNewKF();
1082       if (usingNewKF) {
1083         // Skip crazy tracks to avoid crash (as NewKF applies no cuts to kill them).
1084         constexpr float crazy_z0_cut = 30.;  // Cut to kill any crazy tracks found by New KF (which applies no cuts)
1085         if (fabs(tmp_trk_z0) > crazy_z0_cut)
1086           continue;
1087       }
1088
1089       int tmp_trk_hitpattern = 0;
1090       tmp_trk_hitpattern = (int)iterL1Track->hitPattern();
1091       hph::HitPatternHelper hph(hphSetup, tmp_trk_hitpattern, tmp_trk_tanL, tmp_trk_z0);
1092       std::vector<int> hitpattern_expanded_binary = hph.binary();
1093       int tmp_trk_lhits_hitpattern = 0;
1094       int tmp_trk_dhits_hitpattern = 0;
1095       for (int i = 0; i < (int)hitpattern_expanded_binary.size(); i++) {
1096         if (hitpattern_expanded_binary[i]) {
1097           if (i < 6) {
1098             tmp_trk_lhits_hitpattern += pow(10, i);
1099           } else {
1100             tmp_trk_dhits_hitpattern += pow(10, i - 6);
1101           }
1102         }
1103       }
1104       int tmp_trk_nPSstub_hitpattern = hph.numPS();
1105       int tmp_trk_n2Sstub_hitpattern = hph.num2S();
1106       int tmp_trk_nLostPSstub_hitpattern = hph.numMissingPS();
1107       int tmp_trk_nLost2Sstub_hitpattern = hph.numMissing2S();
1108       int tmp_trk_nLoststub_V1_hitpattern = hph.numMissingInterior1();
1109       int tmp_trk_nLoststub_V2_hitpattern = hph.numMissingInterior2();
1110
1111       float tmp_trk_d0 = -999;
1112       if (L1Tk_nPar == 5) {
1113         float tmp_trk_x0 = iterL1Track->POCA().x();
1114         float tmp_trk_y0 = iterL1Track->POCA().y();
1115         tmp_trk_d0 = tmp_trk_x0 * sin(tmp_trk_phi) - tmp_trk_y0 * cos(tmp_trk_phi);
1116       }
1117
1118       float tmp_trk_chi2 = iterL1Track->chi2();
1119       float tmp_trk_chi2rphi = iterL1Track->chi2XY();
1120       float tmp_trk_chi2rz = iterL1Track->chi2Z();
1121       float tmp_trk_bendchi2 = iterL1Track->stubPtConsistency();
1122       float tmp_trk_MVA1 = iterL1Track->trkMVA1();
1123
1124       std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
1125           stubRefs = iterL1Track->getStubRefs();
1126       int tmp_trk_nstub = (int)stubRefs.size();
1127       int ndof = 2 * tmp_trk_nstub - L1Tk_nPar;
1128       int ndofrphi = tmp_trk_nstub - L1Tk_nPar + 2;
1129       int ndofrz = tmp_trk_nstub - 2;
1130       float tmp_trk_chi2_dof = (float)tmp_trk_chi2 / ndof;
1131       float tmp_trk_chi2rphi_dof = (float)tmp_trk_chi2rphi / ndofrphi;
1132       float tmp_trk_chi2rz_dof = (float)tmp_trk_chi2rz / ndofrz;
1133
1134       int tmp_trk_seed = 0;
1135       tmp_trk_seed = (int)iterL1Track->trackSeedType();
1136
1137       unsigned int tmp_trk_phiSector = iterL1Track->phiSector();
1138       int tmp_trk_etaSector = hph.etaSector();
1139
1140       // ----------------------------------------------------------------------------------------------
1141       // loop over stubs on tracks
1142
1143       //float tmp_trk_bend_chi2 = 0;
1144       int tmp_trk_dhits = 0;
1145       int tmp_trk_lhits = 0;
1146
1147       if (true) {
1148         // loop over stubs
1149         for (int is = 0; is < tmp_trk_nstub; is++) {
1150           //detID of stub
1151           DetId detIdStub = theTrackerGeom->idToDet((stubRefs.at(is)->clusterRef(0))->getDetId())->geographicalId();
1152
1153           MeasurementPoint coords = stubRefs.at(is)->clusterRef(0)->findAverageLocalCoordinatesCentered();
1154           const GeomDet* theGeomDet = theTrackerGeom->idToDet(detIdStub);
1155           Global3DPoint posStub = theGeomDet->surface().toGlobal(theGeomDet->topology().localPosition(coords));
1156
1157           double x = posStub.x();
1158           double y = posStub.y();
1159           double z = posStub.z();
1160
1161           int layer = -999999;
1162           if (detIdStub.subdetId() == StripSubdetector::TOB) {
1163             layer = static_cast<int>(tTopo->layer(detIdStub));
1164             if (DebugMode)
1165               edm::LogVerbatim("Tracklet")
1166                   << "   stub in layer " << layer << " at position x y z = " << x << " " << y << " " << z;
1167             tmp_trk_lhits += pow(10, layer - 1);
1168           } else if (detIdStub.subdetId() == StripSubdetector::TID) {
1169             layer = static_cast<int>(tTopo->layer(detIdStub));
1170             if (DebugMode)
1171               edm::LogVerbatim("Tracklet")
1172                   << "   stub in disk " << layer << " at position x y z = " << x << " " << y << " " << z;
1173             tmp_trk_dhits += pow(10, layer - 1);
1174           }
1175
1176         }  //end loop over stubs
1177       }
1178       // ----------------------------------------------------------------------------------------------
1179
1180       int tmp_trk_genuine = 0;
1181       int tmp_trk_loose = 0;
1182       int tmp_trk_unknown = 0;
1183       int tmp_trk_combinatoric = 0;
1184       if (MCTruthTTTrackHandle->isLooselyGenuine(l1track_ptr))
1185         tmp_trk_loose = 1;
1186       if (MCTruthTTTrackHandle->isGenuine(l1track_ptr))
1187         tmp_trk_genuine = 1;
1188       if (MCTruthTTTrackHandle->isUnknown(l1track_ptr))
1189         tmp_trk_unknown = 1;
1190       if (MCTruthTTTrackHandle->isCombinatoric(l1track_ptr))
1191         tmp_trk_combinatoric = 1;
1192
1193       if (DebugMode) {
1194         edm::LogVerbatim("Tracklet") << "L1 track,"
1195                                      << " pt: " << tmp_trk_pt << " eta: " << tmp_trk_eta << " phi: " << tmp_trk_phi
1196                                      << " z0: " << tmp_trk_z0 << " chi2: " << tmp_trk_chi2
1197                                      << " chi2rphi: " << tmp_trk_chi2rphi << " chi2rz: " << tmp_trk_chi2rz
1198                                      << " nstub: " << tmp_trk_nstub;
1199         if (tmp_trk_genuine)
1200           edm::LogVerbatim("Tracklet") << "    (is genuine)";
1201         if (tmp_trk_unknown)
1202           edm::LogVerbatim("Tracklet") << "    (is unknown)";
1203         if (tmp_trk_combinatoric)
1204           edm::LogVerbatim("Tracklet") << "    (is combinatoric)";
1205       }
1206
1207       m_trk_pt->push_back(tmp_trk_pt);
1208       m_trk_eta->push_back(tmp_trk_eta);
1209       m_trk_phi->push_back(tmp_trk_phi);
1210       m_trk_z0->push_back(tmp_trk_z0);
1211       if (L1Tk_nPar == 5)
1212         m_trk_d0->push_back(tmp_trk_d0);
1213       else
1214         m_trk_d0->push_back(999.);
1215       m_trk_chi2->push_back(tmp_trk_chi2);
1216       m_trk_chi2_dof->push_back(tmp_trk_chi2_dof);
1217       m_trk_chi2rphi->push_back(tmp_trk_chi2rphi);
1218       m_trk_chi2rphi_dof->push_back(tmp_trk_chi2rphi_dof);
1219       m_trk_chi2rz->push_back(tmp_trk_chi2rz);
1220       m_trk_chi2rz_dof->push_back(tmp_trk_chi2rz_dof);
1221       m_trk_bendchi2->push_back(tmp_trk_bendchi2);
1222       m_trk_MVA1->push_back(tmp_trk_MVA1);
1223       m_trk_nstub->push_back(tmp_trk_nstub);
1224       m_trk_dhits->push_back(tmp_trk_dhits);
1225       m_trk_lhits->push_back(tmp_trk_lhits);
1226       m_trk_seed->push_back(tmp_trk_seed);
1227       m_trk_hitpattern->push_back(tmp_trk_hitpattern);
1228       m_trk_lhits_hitpattern->push_back(tmp_trk_lhits_hitpattern);
1229       m_trk_dhits_hitpattern->push_back(tmp_trk_dhits_hitpattern);
1230       m_trk_nPSstub_hitpattern->push_back(tmp_trk_nPSstub_hitpattern);
1231       m_trk_n2Sstub_hitpattern->push_back(tmp_trk_n2Sstub_hitpattern);
1232       m_trk_nLostPSstub_hitpattern->push_back(tmp_trk_nLostPSstub_hitpattern);
1233       m_trk_nLost2Sstub_hitpattern->push_back(tmp_trk_nLost2Sstub_hitpattern);
1234       m_trk_nLoststub_V1_hitpattern->push_back(tmp_trk_nLoststub_V1_hitpattern);
1235       m_trk_nLoststub_V2_hitpattern->push_back(tmp_trk_nLoststub_V2_hitpattern);
1236       m_trk_charge->push_back(tmp_trk_charge);
1237       m_trk_phiSector->push_back(tmp_trk_phiSector);
1238       m_trk_etaSector->push_back(tmp_trk_etaSector);
1239       m_trk_genuine->push_back(tmp_trk_genuine);
1240       m_trk_loose->push_back(tmp_trk_loose);
1241       m_trk_unknown->push_back(tmp_trk_unknown);
1242       m_trk_combinatoric->push_back(tmp_trk_combinatoric);
1243
1244       // ----------------------------------------------------------------------------------------------
1245       // for studying the fake rate
1246       // ----------------------------------------------------------------------------------------------
1247
1248       edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(l1track_ptr);
1249
1250       int myFake = 0;
1251
1252       int tmp_matchtp_pdgid = -999;
1253       float tmp_matchtp_pt = -999;
1254       float tmp_matchtp_eta = -999;
1255       float tmp_matchtp_phi = -999;
1256       float tmp_matchtp_z0 = -999;
1257       float tmp_matchtp_lxy = -999;
1258       float tmp_matchtp_d0 = -999;
1259
1260       if (my_tp.isNull())
1261         myFake = 0;
1262       else {
1263         int tmp_eventid = my_tp->eventId().event();
1264
1265         if (tmp_eventid > 0)
1266           myFake = 2;
1267         else
1268           myFake = 1;
1269
1270         tmp_matchtp_pdgid = my_tp->pdgId();
1271         tmp_matchtp_pt = my_tp->pt();
1272         tmp_matchtp_eta = my_tp->eta();
1273         tmp_matchtp_phi = my_tp->phi();
1274
1275         float tmp_matchtp_vz = my_tp->vz();
1276         float tmp_matchtp_vx = my_tp->vx();
1277         float tmp_matchtp_vy = my_tp->vy();
1278         tmp_matchtp_lxy = sqrt(tmp_matchtp_vx * tmp_matchtp_vx + tmp_matchtp_vy * tmp_matchtp_vy);
1279
1280         // ----------------------------------------------------------------------------------------------
1281         // get d0/z0 propagated back to the IP
1282
1283         float tmp_matchtp_t = 1.0 / tan(2.0 * atan(exp(-tmp_matchtp_eta)));
1284
1285         float delx = -tmp_matchtp_vx;
1286         float dely = -tmp_matchtp_vy;
1287
1288         float b_field = bFieldHandle.product()->inTesla(GlobalPoint(0, 0, 0)).z();
1289         float c_converted = CLHEP::c_light / 1.0E5;
1290         float r2_inv = my_tp->charge() * c_converted * b_field / tmp_matchtp_pt / 2.0;
1291
1292         float tmp_matchtp_x0p = delx - (1. / (2. * r2_inv) * sin(tmp_matchtp_phi));
1293         float tmp_matchtp_y0p = dely + (1. / (2. * r2_inv) * cos(tmp_matchtp_phi));
1294         float tmp_matchtp_rp = sqrt(tmp_matchtp_x0p * tmp_matchtp_x0p + tmp_matchtp_y0p * tmp_matchtp_y0p);
1295         tmp_matchtp_d0 = my_tp->charge() * tmp_matchtp_rp - (1. / (2. * r2_inv));
1296
1297         static double pi = M_PI;
1298         float delphi = tmp_matchtp_phi - atan2(-r2_inv * tmp_matchtp_x0p, r2_inv * tmp_matchtp_y0p);
1299         if (delphi < -pi)
1300           delphi += 2.0 * pi;
1301         if (delphi > pi)
1302           delphi -= 2.0 * pi;
1303         tmp_matchtp_z0 = tmp_matchtp_vz + tmp_matchtp_t * delphi / (2.0 * r2_inv);
1304         // ----------------------------------------------------------------------------------------------
1305
1306         if (DebugMode) {
1307           edm::LogVerbatim("Tracklet") << "TP matched to track has pt = " << my_tp->p4().pt()
1308                                        << " eta = " << my_tp->momentum().eta() << " phi = " << my_tp->momentum().phi()
1309                                        << " z0 = " << my_tp->vertex().z() << " pdgid = " << my_tp->pdgId()
1310                                        << " lxy = " << tmp_matchtp_lxy;
1311         }
1312       }
1313
1314       m_trk_fake->push_back(myFake);
1315
1316       m_trk_matchtp_pdgid->push_back(tmp_matchtp_pdgid);
1317       m_trk_matchtp_pt->push_back(tmp_matchtp_pt);
1318       m_trk_matchtp_eta->push_back(tmp_matchtp_eta);
1319       m_trk_matchtp_phi->push_back(tmp_matchtp_phi);
1320       m_trk_matchtp_z0->push_back(tmp_matchtp_z0);
1321       m_trk_matchtp_lxy->push_back(tmp_matchtp_lxy);
1322       m_trk_matchtp_d0->push_back(tmp_matchtp_d0);
1323
1324       // ----------------------------------------------------------------------------------------------
1325       // for tracking in jets
1326       // ----------------------------------------------------------------------------------------------
1327
1328       if (TrackingInJets) {
1329         if (DebugMode)
1330           edm::LogVerbatim("Tracklet") << "doing tracking in jets now";
1331
1332         int InJet = 0;
1333         int InJetHighpt = 0;
1334         int InJetVeryHighpt = 0;
1335
1336         for (int ij = 0; ij < (int)v_jets.size(); ij++) {
1337           float deta = tmp_trk_eta - (v_jets.at(ij)).eta();
1338           float dphi = tmp_trk_phi - (v_jets.at(ij)).phi();
1339           while (dphi > 3.14159)
1340             dphi = std::abs(2 * 3.14159 - dphi);
1341           float dR = sqrt(deta * deta + dphi * dphi);
1342
1343           if (dR < 0.4) {
1344             InJet = 1;
1345             if (v_jets_highpt.at(ij) == 1)
1346               InJetHighpt = 1;
1347             if (v_jets_vhighpt.at(ij) == 1)
1348               InJetVeryHighpt = 1;
1349             if (ij < NJETS)
1350               jets_trk_sumpt[ij] += tmp_trk_pt;
1351           }
1352         }
1353
1354         m_trk_injet->push_back(InJet);
1355         m_trk_injet_highpt->push_back(InJetHighpt);
1356         m_trk_injet_vhighpt->push_back(InJetVeryHighpt);
1357
1358       }  //end tracking in jets
1359
1360       // layer encoding
1361       const TTBV hitPattern((int)iterL1Track->hitPattern(), setup->numLayers());
1362       const double zT = iterL1Track->z0() + setup->chosenRofZ() * iterL1Track->tanL();
1363       const vector<int>& le = layerEncoding->layerEncoding(zT);
1364       vector<int> layers;
1365       layers.reserve(hitPattern.size());
1366       for (int layer : hitPattern.ids())
1367         layers.push_back(le[layer]);
1368       m_trk_layers->push_back(layers);
1369
1370     }  //end track loop
1371
1372   }  //end if SaveAllTracks
1373
1374   // ----------------------------------------------------------------------------------------------
1375   // loop over tracking particles
1376   // ----------------------------------------------------------------------------------------------
1377
1378   if (DebugMode)
1379     edm::LogVerbatim("Tracklet") << "\n Loop over tracking particles!";
1380
1381   int this_tp = 0;
1382   std::vector<TrackingParticle>::const_iterator iterTP;
1383   for (iterTP = TrackingParticleHandle->begin(); iterTP != TrackingParticleHandle->end(); ++iterTP) {
1384     edm::Ptr<TrackingParticle> tp_ptr(TrackingParticleHandle, this_tp);
1385     this_tp++;
1386
1387     int tmp_eventid = iterTP->eventId().event();
1388     if (MyProcess != 1 && tmp_eventid > 0)
1389       continue;  //only care about tracking particles from the primary interaction (except for MyProcess==1, i.e. looking at all TPs)
1390
1391     float tmp_tp_pt = iterTP->pt();
1392     float tmp_tp_charge = iterTP->charge();
1393     float tmp_tp_eta = iterTP->eta();
1394
1395     if (tmp_tp_pt < TP_minPt)  // Save CPU by applying these cuts here.
1396       continue;
1397     if (tmp_tp_charge == 0.)
1398       continue;
1399     if (std::abs(tmp_tp_eta) > TP_maxEta)
1400       continue;
1401
1402     float tmp_tp_phi = iterTP->phi();
1403     float tmp_tp_vz = iterTP->vz();
1404     float tmp_tp_vx = iterTP->vx();
1405     float tmp_tp_vy = iterTP->vy();
1406     int tmp_tp_pdgid = iterTP->pdgId();
1407     float tmp_tp_z0_prod = tmp_tp_vz;
1408     float tmp_tp_d0_prod = tmp_tp_vx * sin(tmp_tp_phi) - tmp_tp_vy * cos(tmp_tp_phi);
1409
1410     // ----------------------------------------------------------------------------------------------
1411     // get d0/z0 propagated back to the IP
1412
1413     float tmp_tp_t = 1.0 / tan(2.0 * atan(exp(-tmp_tp_eta)));
1414
1415     float delx = -tmp_tp_vx;
1416     float dely = -tmp_tp_vy;
1417
1418     float b_field = bFieldHandle.product()->inTesla(GlobalPoint(0, 0, 0)).z();
1419     float c_converted = CLHEP::c_light / 1.0E5;
1420     float r2_inv = tmp_tp_charge * c_converted * b_field / tmp_tp_pt / 2.0;
1421
1422     float tmp_tp_x0p = delx - (1. / (2. * r2_inv) * sin(tmp_tp_phi));
1423     float tmp_tp_y0p = dely + (1. / (2. * r2_inv) * cos(tmp_tp_phi));
1424     float tmp_tp_rp = sqrt(tmp_tp_x0p * tmp_tp_x0p + tmp_tp_y0p * tmp_tp_y0p);
1425     float tmp_tp_d0 = tmp_tp_charge * tmp_tp_rp - (1. / (2. * r2_inv));
1426
1427     static double pi = M_PI;
1428     float delphi = tmp_tp_phi - atan2(-r2_inv * tmp_tp_x0p, r2_inv * tmp_tp_y0p);
1429     if (delphi < -pi)
1430       delphi += 2.0 * pi;
1431     if (delphi > pi)
1432       delphi -= 2.0 * pi;
1433     float tmp_tp_z0 = tmp_tp_vz + tmp_tp_t * delphi / (2.0 * r2_inv);
1434     // ----------------------------------------------------------------------------------------------
1435
1436     if (MyProcess == 13 && abs(tmp_tp_pdgid) != 13)
1437       continue;
1438     if (MyProcess == 11 && abs(tmp_tp_pdgid) != 11)
1439       continue;
1440     if ((MyProcess == 6 || MyProcess == 15 || MyProcess == 211) && abs(tmp_tp_pdgid) != 211)
1441       continue;
1442
1443     if (std::abs(tmp_tp_z0) > TP_maxZ0)
1444       continue;
1445
1446     // for pions in ttbar, only consider TPs coming from near the IP!
1447     float lxy = sqrt(tmp_tp_vx * tmp_tp_vx + tmp_tp_vy * tmp_tp_vy);
1448     float tmp_tp_lxy = lxy;
1449     if (MyProcess == 6 && (lxy > 1.0))
1450       continue;
1451
1452     if (DebugMode)
1453       edm::LogVerbatim("Tracklet") << "Tracking particle, pt: " << tmp_tp_pt << " eta: " << tmp_tp_eta
1454                                    << " phi: " << tmp_tp_phi << " z0: " << tmp_tp_z0 << " d0: " << tmp_tp_d0
1455                                    << " z_prod: " << tmp_tp_z0_prod << " d_prod: " << tmp_tp_d0_prod
1456                                    << " pdgid: " << tmp_tp_pdgid << " eventID: " << iterTP->eventId().event()
1457                                    << " ttclusters " << MCTruthTTClusterHandle->findTTClusterRefs(tp_ptr).size()
1458                                    << " ttstubs " << MCTruthTTStubHandle->findTTStubRefs(tp_ptr).size() << " tttracks "
1459                                    << MCTruthTTTrackHandle->findTTTrackPtrs(tp_ptr).size();
1460
1461     // ----------------------------------------------------------------------------------------------
1462     // only consider TPs associated with >= 1 cluster, or >= X stubs, or have stubs in >= X layers (configurable options)
1463
1464     if (MCTruthTTClusterHandle->findTTClusterRefs(tp_ptr).empty()) {
1465       if (DebugMode)
1466         edm::LogVerbatim("Tracklet") << "No matching TTClusters for TP, continuing...";
1467       continue;
1468     }
1469
1470     std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
1471         theStubRefs = MCTruthTTStubHandle->findTTStubRefs(tp_ptr);
1472     int nStubTP = (int)theStubRefs.size();
1473
1474     // how many layers/disks have stubs?
1475     int hasStubInLayer[11] = {0};
1476     for (auto& theStubRef : theStubRefs) {
1477       DetId detid(theStubRef->getDetId());
1478
1479       int layer = -1;
1480       if (detid.subdetId() == StripSubdetector::TOB) {
1481         layer = static_cast<int>(tTopo->layer(detid)) - 1;  //fill in array as entries 0-5
1482       } else if (detid.subdetId() == StripSubdetector::TID) {
1483         layer = static_cast<int>(tTopo->layer(detid)) + 5;  //fill in array as entries 6-10
1484       }
1485
1486       //bool isPS = (theTrackerGeom->getDetectorType(detid)==TrackerGeometry::ModuleType::Ph2PSP);
1487
1488       //treat genuine stubs separately (==2 is genuine, ==1 is not)
1489       if (MCTruthTTStubHandle->findTrackingParticlePtr(theStubRef).isNull() && hasStubInLayer[layer] < 2)
1490         hasStubInLayer[layer] = 1;
1491       else
1492         hasStubInLayer[layer] = 2;
1493     }
1494
1495     int nStubLayerTP = 0;
1496     int nStubLayerTP_g = 0;
1497     for (int isum : hasStubInLayer) {
1498       if (isum >= 1)
1499         nStubLayerTP += 1;
1500       if (isum == 2)
1501         nStubLayerTP_g += 1;
1502     }
1503
1504     if (DebugMode)
1505       edm::LogVerbatim("Tracklet") << "TP is associated with " << nStubTP << " stubs, and has stubs in " << nStubLayerTP
1506                                    << " different layers/disks, and has GENUINE stubs in " << nStubLayerTP_g
1507                                    << " layers ";
1508
1509     if (TP_minNStub > 0) {
1510       if (DebugMode)
1511         edm::LogVerbatim("Tracklet") << "Only consider TPs with >= " << TP_minNStub << " stubs";
1512       if (nStubTP < TP_minNStub) {
1513         if (DebugMode)
1514           edm::LogVerbatim("Tracklet") << "TP fails minimum nbr stubs requirement! Continuing...";
1515         continue;
1516       }
1517     }
1518     if (TP_minNStubLayer > 0) {
1519       if (DebugMode)
1520         edm::LogVerbatim("Tracklet") << "Only consider TPs with stubs in >= " << TP_minNStubLayer << " layers/disks";
1521       if (nStubLayerTP < TP_minNStubLayer) {
1522         if (DebugMode)
1523           edm::LogVerbatim("Tracklet") << "TP fails stubs in minimum nbr of layers/disks requirement! Continuing...";
1524         continue;
1525       }
1526     }
1527
1528     // ----------------------------------------------------------------------------------------------
1529     // look for L1 tracks matched to the tracking particle
1530
1531     std::vector<edm::Ptr<TTTrack<Ref_Phase2TrackerDigi_>>> matchedTracks =
1532         MCTruthTTTrackHandle->findTTTrackPtrs(tp_ptr);
1533
1534     int nMatch = 0;
1535     int i_track = -1;
1536     float i_chi2dof = 99999;
1537
1538     if (!matchedTracks.empty()) {
1539       if (DebugMode && (matchedTracks.size() > 1))
1540         edm::LogVerbatim("Tracklet") << "TrackingParticle has more than one matched L1 track!";
1541
1542       // ----------------------------------------------------------------------------------------------
1543       // loop over matched L1 tracks
1544       // here, "match" means tracks that can be associated to a TrackingParticle with at least one hit of at least one of its clusters
1545       // https://twiki.cern.ch/twiki/bin/viewauth/CMS/SLHCTrackerTriggerSWTools#MC_truth_for_TTTrack
1546
1547       for (int it = 0; it < (int)matchedTracks.size(); it++) {
1548         bool tmp_trk_genuine = false;
1549         bool tmp_trk_loosegenuine = false;
1550         if (MCTruthTTTrackHandle->isGenuine(matchedTracks.at(it)))
1551           tmp_trk_genuine = true;
1552         if (MCTruthTTTrackHandle->isLooselyGenuine(matchedTracks.at(it)))
1553           tmp_trk_loosegenuine = true;
1554         if (!tmp_trk_loosegenuine)
1555           continue;
1556
1557         if (DebugMode) {
1558           if (MCTruthTTTrackHandle->findTrackingParticlePtr(matchedTracks.at(it)).isNull()) {
1559             edm::LogVerbatim("Tracklet") << "track matched to TP is NOT uniquely matched to a TP";
1560           } else {
1561             edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(matchedTracks.at(it));
1562             edm::LogVerbatim("Tracklet") << "TP matched to track matched to TP ... tp pt = " << my_tp->p4().pt()
1563                                          << " eta = " << my_tp->momentum().eta() << " phi = " << my_tp->momentum().phi()
1564                                          << " z0 = " << my_tp->vertex().z();
1565           }
1566           edm::LogVerbatim("Tracklet") << "   ... matched L1 track has pt = " << matchedTracks.at(it)->momentum().perp()
1567                                        << " eta = " << matchedTracks.at(it)->momentum().eta()
1568                                        << " phi = " << matchedTracks.at(it)->momentum().phi()
1569                                        << " chi2 = " << matchedTracks.at(it)->chi2()
1570                                        << " consistency = " << matchedTracks.at(it)->stubPtConsistency()
1571                                        << " z0 = " << matchedTracks.at(it)->z0()
1572                                        << " nstub = " << matchedTracks.at(it)->getStubRefs().size();
1573           if (tmp_trk_genuine)
1574             edm::LogVerbatim("Tracklet") << "    (genuine!) ";
1575           if (tmp_trk_loosegenuine)
1576             edm::LogVerbatim("Tracklet") << "    (loose genuine!) ";
1577         }
1578
1579         // ----------------------------------------------------------------------------------------------
1580         // further require L1 track to be (loosely) genuine, that there is only one TP matched to the track
1581         // + have >= L1Tk_minNStub stubs for it to be a valid match (only relevant is your track collection
1582         // e.g. stores 3-stub tracks but at plot level you require >= 4 stubs (--> tracklet case)
1583
1584         std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
1585             stubRefs = matchedTracks.at(it)->getStubRefs();
1586         int tmp_trk_nstub = stubRefs.size();
1587
1588         if (tmp_trk_nstub < L1Tk_minNStub)
1589           continue;
1590
1591         /*
1592     // PS stubs
1593     int tmp_trk_nPSstub = 0;
1594     for (int is=0; is<tmp_trk_nstub; is++) {
1595       DetId detIdStub = theTrackerGeom->idToDet( (stubRefs.at(is)->clusterRef(0))->getDetId() )->geographicalId();
1596       DetId stackDetid = tTopo->stack(detIdStub);
1597       bool isPS = (theTrackerGeom->getDetectorType(stackDetid)==TrackerGeometry::ModuleType::Ph2PSP);
1598       if (isPS) tmp_trk_nPSstub++;
1599     }
1600     */
1601
1602         float dmatch_pt = 999;
1603         float dmatch_eta = 999;
1604         float dmatch_phi = 999;
1605         int match_id = 999;
1606
1607         edm::Ptr<TrackingParticle> my_tp = MCTruthTTTrackHandle->findTrackingParticlePtr(matchedTracks.at(it));
1608         dmatch_pt = std::abs(my_tp->p4().pt() - tmp_tp_pt);
1609         dmatch_eta = std::abs(my_tp->p4().eta() - tmp_tp_eta);
1610         dmatch_phi = std::abs(my_tp->p4().phi() - tmp_tp_phi);
1611         match_id = my_tp->pdgId();
1612
1613         float tmp_trk_chi2dof = (matchedTracks.at(it)->chi2()) / (2 * tmp_trk_nstub - L1Tk_nPar);
1614
1615         // ensure that track is uniquely matched to the TP we are looking at!
1616         if (dmatch_pt < 0.1 && dmatch_eta < 0.1 && dmatch_phi < 0.1 && tmp_tp_pdgid == match_id && tmp_trk_genuine) {
1617           nMatch++;
1618           if (i_track < 0 || tmp_trk_chi2dof < i_chi2dof) {
1619             i_track = it;
1620             i_chi2dof = tmp_trk_chi2dof;
1621           }
1622         }
1623
1624       }  // end loop over matched L1 tracks
1625
1626     }  // end has at least 1 matched L1 track
1627     // ----------------------------------------------------------------------------------------------
1628
1629     float tmp_matchtrk_pt = -999;
1630     float tmp_matchtrk_eta = -999;
1631     float tmp_matchtrk_phi = -999;
1632     float tmp_matchtrk_z0 = -999;
1633     float tmp_matchtrk_d0 = -999;
1634     float tmp_matchtrk_chi2 = -999;
1635     float tmp_matchtrk_chi2_dof = -999;
1636     float tmp_matchtrk_chi2rphi = -999;
1637     float tmp_matchtrk_chi2rphi_dof = -999;
1638     float tmp_matchtrk_chi2rz = -999;
1639     float tmp_matchtrk_chi2rz_dof = -999;
1640     float tmp_matchtrk_bendchi2 = -999;
1641     float tmp_matchtrk_MVA1 = -999;
1642     int tmp_matchtrk_charge = -999;
1643     int tmp_matchtrk_nstub = -999;
1644     int tmp_matchtrk_dhits = -999;
1645     int tmp_matchtrk_lhits = -999;
1646     int tmp_matchtrk_seed = -999;
1647     int tmp_matchtrk_hitpattern = -999;
1648
1649     if (nMatch > 1 && DebugMode)
1650       edm::LogVerbatim("Tracklet") << "WARNING *** 2 or more matches to genuine L1 tracks ***";
1651
1652     if (nMatch > 0) {
1653       tmp_matchtrk_pt = matchedTracks.at(i_track)->momentum().perp();
1654       tmp_matchtrk_charge = (int)TMath::Sign(1, matchedTracks.at(i_track)->rInv());
1655       tmp_matchtrk_eta = matchedTracks.at(i_track)->momentum().eta();
1656       tmp_matchtrk_phi = matchedTracks.at(i_track)->momentum().phi();
1657       tmp_matchtrk_z0 = matchedTracks.at(i_track)->z0();
1658
1659       if (L1Tk_nPar == 5) {
1660         float tmp_matchtrk_x0 = matchedTracks.at(i_track)->POCA().x();
1661         float tmp_matchtrk_y0 = matchedTracks.at(i_track)->POCA().y();
1662         tmp_matchtrk_d0 = tmp_matchtrk_x0 * sin(tmp_matchtrk_phi) - tmp_matchtrk_y0 * cos(tmp_matchtrk_phi);
1663       }
1664
1665       tmp_matchtrk_chi2 = matchedTracks.at(i_track)->chi2();
1666       tmp_matchtrk_chi2rphi = matchedTracks.at(i_track)->chi2XY();
1667       tmp_matchtrk_chi2rz = matchedTracks.at(i_track)->chi2Z();
1668       tmp_matchtrk_bendchi2 = matchedTracks.at(i_track)->stubPtConsistency();
1669       tmp_matchtrk_MVA1 = matchedTracks.at(i_track)->trkMVA1();
1670       tmp_matchtrk_nstub = (int)matchedTracks.at(i_track)->getStubRefs().size();
1671       tmp_matchtrk_seed = (int)matchedTracks.at(i_track)->trackSeedType();
1672       tmp_matchtrk_hitpattern = (int)matchedTracks.at(i_track)->hitPattern();
1673
1674       int ndof = 2 * tmp_matchtrk_nstub - L1Tk_nPar;
1675       int ndofrphi = tmp_matchtrk_nstub - L1Tk_nPar + 2;
1676       int ndofrz = tmp_matchtrk_nstub - 2;
1677       tmp_matchtrk_chi2_dof = (float)tmp_matchtrk_chi2 / ndof;
1678       tmp_matchtrk_chi2rphi_dof = (float)tmp_matchtrk_chi2rphi / ndofrphi;
1679       tmp_matchtrk_chi2rz_dof = (float)tmp_matchtrk_chi2rz / ndofrz;
1680
1681       // ------------------------------------------------------------------------------------------
1682
1683       //float tmp_matchtrk_bend_chi2 = 0;
1684
1685       tmp_matchtrk_dhits = 0;
1686       tmp_matchtrk_lhits = 0;
1687
1688       std::vector<edm::Ref<edmNew::DetSetVector<TTStub<Ref_Phase2TrackerDigi_>>, TTStub<Ref_Phase2TrackerDigi_>>>
1689           stubRefs = matchedTracks.at(i_track)->getStubRefs();
1690       int tmp_nstub = stubRefs.size();
1691
1692       for (int is = 0; is < tmp_nstub; is++) {
1693         DetId detIdStub = theTrackerGeom->idToDet((stubRefs.at(is)->clusterRef(0))->getDetId())->geographicalId();
1694         /*
1695     MeasurementPoint coords = stubRefs.at(is)->clusterRef(0)->findAverageLocalCoordinatesCentered();
1696     const GeomDet* theGeomDet = theTrackerGeom->idToDet(detIdStub);
1697     Global3DPoint posStub = theGeomDet->surface().toGlobal( theGeomDet->topology().localPosition(coords) );
1698     */
1699
1700         int layer = -999999;
1701         if (detIdStub.subdetId() == StripSubdetector::TOB) {
1702           layer = static_cast<int>(tTopo->layer(detIdStub));
1703           tmp_matchtrk_lhits += pow(10, layer - 1);
1704         } else if (detIdStub.subdetId() == StripSubdetector::TID) {
1705           layer = static_cast<int>(tTopo->layer(detIdStub));
1706           tmp_matchtrk_dhits += pow(10, layer - 1);
1707         }
1708
1709         // ------------------------------------------------------------------------------------------
1710       }
1711     }
1712
1713     m_tp_pt->push_back(tmp_tp_pt);
1714     m_tp_eta->push_back(tmp_tp_eta);
1715     m_tp_phi->push_back(tmp_tp_phi);
1716     m_tp_lxy->push_back(tmp_tp_lxy);
1717     m_tp_z0->push_back(tmp_tp_z0);
1718     m_tp_d0->push_back(tmp_tp_d0);
1719     m_tp_z0_prod->push_back(tmp_tp_z0_prod);
1720     m_tp_d0_prod->push_back(tmp_tp_d0_prod);
1721     m_tp_pdgid->push_back(tmp_tp_pdgid);
1722     m_tp_nmatch->push_back(nMatch);
1723     m_tp_nstub->push_back(nStubTP);
1724     m_tp_eventid->push_back(tmp_eventid);
1725     m_tp_charge->push_back(tmp_tp_charge);
1726
1727     m_matchtrk_pt->push_back(tmp_matchtrk_pt);
1728     m_matchtrk_eta->push_back(tmp_matchtrk_eta);
1729     m_matchtrk_phi->push_back(tmp_matchtrk_phi);
1730     m_matchtrk_z0->push_back(tmp_matchtrk_z0);
1731     m_matchtrk_d0->push_back(tmp_matchtrk_d0);
1732     m_matchtrk_chi2->push_back(tmp_matchtrk_chi2);
1733     m_matchtrk_chi2rphi->push_back(tmp_matchtrk_chi2rphi);
1734     m_matchtrk_chi2rz->push_back(tmp_matchtrk_chi2rz);
1735     m_matchtrk_bendchi2->push_back(tmp_matchtrk_bendchi2);
1736     m_matchtrk_MVA1->push_back(tmp_matchtrk_MVA1);
1737     m_matchtrk_nstub->push_back(tmp_matchtrk_nstub);
1738     m_matchtrk_dhits->push_back(tmp_matchtrk_dhits);
1739     m_matchtrk_lhits->push_back(tmp_matchtrk_lhits);
1740     m_matchtrk_seed->push_back(tmp_matchtrk_seed);
1741     m_matchtrk_hitpattern->push_back(tmp_matchtrk_hitpattern);
1742     m_matchtrk_charge->push_back(tmp_matchtrk_charge);
1743     m_matchtrk_chi2_dof->push_back(tmp_matchtrk_chi2_dof);
1744     m_matchtrk_chi2rphi_dof->push_back(tmp_matchtrk_chi2rphi_dof);
1745     m_matchtrk_chi2rz_dof->push_back(tmp_matchtrk_chi2rz_dof);
1746
1747     // ----------------------------------------------------------------------------------------------
1748     // for tracking in jets
1749     // ----------------------------------------------------------------------------------------------
1750
1751     if (TrackingInJets) {
1752       if (DebugMode)
1753         edm::LogVerbatim("Tracklet") << "check if TP/matched track is within jet";
1754
1755       int tp_InJet = 0;
1756       int matchtrk_InJet = 0;
1757       int tp_InJetHighpt = 0;
1758       int matchtrk_InJetHighpt = 0;
1759       int tp_InJetVeryHighpt = 0;
1760       int matchtrk_InJetVeryHighpt = 0;
1761
1762       for (int ij = 0; ij < (int)v_jets.size(); ij++) {
1763         float deta = tmp_tp_eta - (v_jets.at(ij)).eta();
1764         float dphi = tmp_tp_phi - (v_jets.at(ij)).phi();
1765         while (dphi > 3.14159)
1766           dphi = std::abs(2 * 3.14159 - dphi);
1767         float dR = sqrt(deta * deta + dphi * dphi);
1768         if (dR < 0.4) {
1769           tp_InJet = 1;
1770           if (v_jets_highpt.at(ij) == 1)
1771             tp_InJetHighpt = 1;
1772           if (v_jets_vhighpt.at(ij) == 1)
1773             tp_InJetVeryHighpt = 1;
1774           if (ij < NJETS)
1775             jets_tp_sumpt[ij] += tmp_tp_pt;
1776         }
1777
1778         if (nMatch > 0) {
1779           deta = tmp_matchtrk_eta - (v_jets.at(ij)).eta();
1780           dphi = tmp_matchtrk_phi - (v_jets.at(ij)).phi();
1781           while (dphi > 3.14159)
1782             dphi = std::abs(2 * 3.14159 - dphi);
1783           dR = sqrt(deta * deta + dphi * dphi);
1784           if (dR < 0.4) {
1785             matchtrk_InJet = 1;
1786             if (v_jets_highpt.at(ij) == 1)
1787               matchtrk_InJetHighpt = 1;
1788             if (v_jets_vhighpt.at(ij) == 1)
1789               matchtrk_InJetVeryHighpt = 1;
1790             if (ij < NJETS)
1791               jets_matchtrk_sumpt[ij] += tmp_matchtrk_pt;
1792           }
1793         }
1794       }
1795
1796       m_tp_injet->push_back(tp_InJet);
1797       m_tp_injet_highpt->push_back(tp_InJetHighpt);
1798       m_tp_injet_vhighpt->push_back(tp_InJetVeryHighpt);
1799       m_matchtrk_injet->push_back(matchtrk_InJet);
1800       m_matchtrk_injet_highpt->push_back(matchtrk_InJetHighpt);
1801       m_matchtrk_injet_vhighpt->push_back(matchtrk_InJetVeryHighpt);
1802
1803     }  //end TrackingInJets
1804
1805   }  //end loop tracking particles
1806
1807   if (TrackingInJets) {
1808     for (int ij = 0; ij < (int)v_jets.size(); ij++) {
1809       if (ij < NJETS) {
1810         m_jet_eta->push_back((v_jets.at(ij)).eta());
1811         m_jet_phi->push_back((v_jets.at(ij)).phi());
1812         m_jet_pt->push_back((v_jets.at(ij)).pt());
1813         m_jet_tp_sumpt->push_back(jets_tp_sumpt[ij]);
1814         m_jet_trk_sumpt->push_back(jets_trk_sumpt[ij]);
1815         m_jet_matchtrk_sumpt->push_back(jets_matchtrk_sumpt[ij]);
1816       }
1817     }
1818   }
1819
1820   eventTree->Fill();
1821
1822 }  // end of analyze()
1823
1824 ///////////////////////////
1825 // DEFINE THIS AS A PLUG-IN
1826 DEFINE_FWK_MODULE(L1TrackNtupleMaker);