Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:14:01 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!