File indexing completed on 2024-04-06 11:59:17
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <unordered_map>
0022
0023
0024 #include "TROOT.h"
0025 #include "TH1.h"
0026 #include "TH2.h"
0027 #include "TSystem.h"
0028 #include "TFile.h"
0029 #include "TProfile.h"
0030 #include "TDirectory.h"
0031 #include "TTree.h"
0032 #include "TMath.h"
0033
0034
0035 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0036 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0037 #include "Calibration/IsolatedParticles/interface/eCone.h"
0038 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0039 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0040
0041 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0042
0043 #include "DataFormats/Common/interface/RefToBase.h"
0044 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0045 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0046 #include "DataFormats/Math/interface/LorentzVector.h"
0047 #include "DataFormats/Math/interface/Point3D.h"
0048 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0049 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0050 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0053 #include "DataFormats/TrackReco/interface/TrackBase.h"
0054
0055 #include "DataFormats/TrackReco/interface/HitPattern.h"
0056
0057 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0058 #include "DataFormats/VertexReco/interface/Vertex.h"
0059 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0060
0061 #include "DataFormats/Common/interface/TriggerResults.h"
0062 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0063 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0064 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0065 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0066
0067 #include "FWCore/Common/interface/TriggerNames.h"
0068 #include "FWCore/Framework/interface/Frameworkfwd.h"
0069 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0070 #include "FWCore/Framework/interface/Event.h"
0071 #include "FWCore/Framework/interface/MakerMacros.h"
0072 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0073 #include "FWCore/ServiceRegistry/interface/Service.h"
0074 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0075 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0076
0077 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0078 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0079 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0080 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0081 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0082 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0083 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0084
0085 #include "HLTrigger/HLTcore/interface/HLTPrescaleProvider.h"
0086 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0087
0088 #include "MagneticField/Engine/interface/MagneticField.h"
0089 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0090
0091 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0092 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0093
0094 #include "TrackingTools/TransientTrackingRecHit/interface/SeedingLayerSetsHits.h"
0095
0096 class IsoTrig : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0097 public:
0098 explicit IsoTrig(const edm::ParameterSet &);
0099 ~IsoTrig() override;
0100
0101 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0102
0103 private:
0104 void analyze(const edm::Event &, const edm::EventSetup &) override;
0105 void beginJob() override;
0106 void endJob() override;
0107 void beginRun(edm::Run const &, edm::EventSetup const &) override;
0108 void endRun(edm::Run const &, edm::EventSetup const &) override {}
0109
0110 void clearMipCutTreeVectors();
0111 void clearChgIsolnTreeVectors();
0112 void pushChgIsolnTreeVecs(math::XYZTLorentzVector &Pixcand,
0113 math::XYZTLorentzVector &Trkcand,
0114 std::vector<double> &PixMaxP,
0115 double &TrkMaxP,
0116 bool &selTk);
0117 void pushMipCutTreeVecs(math::XYZTLorentzVector &NFcand,
0118 math::XYZTLorentzVector &Trkcand,
0119 double &EmipNFcand,
0120 double &EmipTrkcand,
0121 double &mindR,
0122 double &mindP1,
0123 std::vector<bool> &Flags,
0124 double hCone);
0125 void StudyTrkEbyP(edm::Handle<reco::TrackCollection> &trkCollection);
0126 void studyTiming(const edm::Event &theEvent);
0127 void studyMipCut(edm::Handle<reco::TrackCollection> &trkCollection,
0128 edm::Handle<reco::IsolatedPixelTrackCandidateCollection> &L2cands);
0129 void studyTrigger(edm::Handle<reco::TrackCollection> &, std::vector<reco::TrackCollection::const_iterator> &);
0130 void studyIsolation(edm::Handle<reco::TrackCollection> &, std::vector<reco::TrackCollection::const_iterator> &);
0131 void chgIsolation(double &etaTriggered,
0132 double &phiTriggered,
0133 edm::Handle<reco::TrackCollection> &trkCollection,
0134 const edm::Event &theEvent);
0135 void getGoodTracks(const edm::Event &, edm::Handle<reco::TrackCollection> &);
0136 void fillHist(int, math::XYZTLorentzVector &);
0137 void fillDifferences(int, math::XYZTLorentzVector &, math::XYZTLorentzVector &, bool);
0138 void fillCuts(int, double, double, double, math::XYZTLorentzVector &, int, bool);
0139 void fillEnergy(int, int, double, double, math::XYZTLorentzVector &);
0140 double dEta(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0141 double dPhi(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0142 double dR(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0143 double dPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0144 double dP(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0145 double dinvPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0146 std::pair<double, double> etaPhiTrigger();
0147 std::pair<double, double> GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ);
0148 double getDistInCM(double eta1, double phi1, double eta2, double phi2);
0149
0150
0151 HLTPrescaleProvider hltPrescaleProvider_;
0152 const std::vector<std::string> trigNames_;
0153 const edm::InputTag pixCandTag_, l1CandTag_, l2CandTag_;
0154 const std::vector<edm::InputTag> pixelTracksSources_;
0155 const bool doL2L3_, doTiming_, doMipCutTree_;
0156 const bool doTrkResTree_, doChgIsolTree_, doStudyIsol_;
0157 const int verbosity_;
0158 const std::vector<double> pixelIsolationConeSizeAtEC_;
0159 const double minPTrackValue_, vtxCutSeed_, vtxCutIsol_;
0160 const double tauUnbiasCone_, prelimCone_;
0161 std::string theTrackQuality_;
0162 const std::string processName_;
0163 double rEB_, zEE_, bfVal_;
0164 spr::trackSelectionParameters selectionParameters_;
0165 const double dr_L1_, a_coneR_, a_charIsoR_, a_neutIsoR_;
0166 const double a_mipR_, a_neutR1_, a_neutR2_, cutMip_;
0167 const double cutCharge_, cutNeutral_;
0168 const int minRunNo_, maxRunNo_;
0169 edm::EDGetTokenT<LumiDetails> tok_lumi_;
0170 edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0171 edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0172 edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hlt_;
0173 edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0174 edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0175 edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0176 edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0177 edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0178 edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0179 edm::EDGetTokenT<reco::VertexCollection> tok_verthb_, tok_verthe_;
0180 edm::EDGetTokenT<SeedingLayerSetsHits> tok_SeedingLayerHB_;
0181 edm::EDGetTokenT<SeedingLayerSetsHits> tok_SeedingLayerHE_;
0182 edm::EDGetTokenT<SiPixelRecHitCollection> tok_SiPixelRecHits_;
0183 edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_pixtk_;
0184 edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_l1cand_;
0185 edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_l2cand_;
0186 std::vector<edm::EDGetTokenT<reco::TrackCollection>> tok_pixtks_;
0187
0188 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0189 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0190
0191 std::vector<reco::TrackRef> pixelTrackRefsHB_, pixelTrackRefsHE_;
0192 edm::Handle<HBHERecHitCollection> hbhe_;
0193 edm::Handle<EcalRecHitCollection> barrelRecHitsHandle_;
0194 edm::Handle<EcalRecHitCollection> endcapRecHitsHandle_;
0195 edm::Handle<reco::BeamSpot> beamSpotH_;
0196 edm::Handle<reco::VertexCollection> recVtxs_;
0197
0198 const MagneticField *bField_;
0199 const CaloGeometry *geo_;
0200 math::XYZPoint leadPV_;
0201
0202 std::unordered_map<unsigned int, unsigned int> trigList_;
0203 std::unordered_map<unsigned int, const std::pair<double, double>> trigPreList_;
0204 bool changed_;
0205 double pLimits_[6];
0206 edm::Service<TFileService> fs_;
0207 TTree *MipCutTree_, *ChgIsolnTree_, *TrkResTree_, *TimingTree_;
0208 std::vector<double> *t_timeL2Prod;
0209 std::vector<int> *t_nPixCand;
0210 std::vector<int> *t_nPixSeed;
0211 std::vector<int> *t_nGoodTk;
0212
0213 std::vector<double> *t_TrkhCone;
0214 std::vector<double> *t_TrkP;
0215 std::vector<bool> *t_TrkselTkFlag;
0216 std::vector<bool> *t_TrkqltyFlag;
0217 std::vector<bool> *t_TrkMissFlag;
0218 std::vector<bool> *t_TrkPVFlag;
0219 std::vector<bool> *t_TrkNuIsolFlag;
0220
0221 std::vector<double> *t_PixcandP;
0222 std::vector<double> *t_PixcandPt;
0223 std::vector<double> *t_PixcandEta;
0224 std::vector<double> *t_PixcandPhi;
0225 std::vector<std::vector<double>> *t_PixcandMaxP;
0226 std::vector<double> *t_PixTrkcandP;
0227 std::vector<double> *t_PixTrkcandPt;
0228 std::vector<double> *t_PixTrkcandEta;
0229 std::vector<double> *t_PixTrkcandPhi;
0230 std::vector<double> *t_PixTrkcandMaxP;
0231 std::vector<bool> *t_PixTrkcandselTk;
0232
0233 std::vector<double> *t_NFcandP;
0234 std::vector<double> *t_NFcandPt;
0235 std::vector<double> *t_NFcandEta;
0236 std::vector<double> *t_NFcandPhi;
0237 std::vector<double> *t_NFcandEmip;
0238 std::vector<double> *t_NFTrkcandP;
0239 std::vector<double> *t_NFTrkcandPt;
0240 std::vector<double> *t_NFTrkcandEta;
0241 std::vector<double> *t_NFTrkcandPhi;
0242 std::vector<double> *t_NFTrkcandEmip;
0243 std::vector<double> *t_NFTrkMinDR;
0244 std::vector<double> *t_NFTrkMinDP1;
0245 std::vector<bool> *t_NFTrkselTkFlag;
0246 std::vector<bool> *t_NFTrkqltyFlag;
0247 std::vector<bool> *t_NFTrkMissFlag;
0248 std::vector<bool> *t_NFTrkPVFlag;
0249 std::vector<bool> *t_NFTrkPropFlag;
0250 std::vector<bool> *t_NFTrkChgIsoFlag;
0251 std::vector<bool> *t_NFTrkNeuIsoFlag;
0252 std::vector<bool> *t_NFTrkMipFlag;
0253 std::vector<double> *t_ECone;
0254
0255 TH1D *h_EnIn, *h_EnOut;
0256 TH2D *h_MipEnMatch, *h_MipEnNoMatch;
0257 TH1I *h_nHLT, *h_HLT, *h_PreL1, *h_PreHLT;
0258 TH1I *h_Pre, *h_nL3Objs, *h_Filters;
0259 TH1D *h_PreL1wt, *h_PreHLTwt, *h_L1ObjEnergy;
0260 TH1D *h_p[20], *h_pt[20], *h_eta[20], *h_phi[20];
0261 TH1D *h_dEtaL1[2], *h_dPhiL1[2], *h_dRL1[2];
0262 TH1D *h_dEta[9], *h_dPhi[9], *h_dPt[9], *h_dP[9];
0263 TH1D *h_dinvPt[9], *h_mindR[9], *h_eMip[2];
0264 TH1D *h_eMaxNearP[2], *h_eNeutIso[2];
0265 TH1D *h_etaCalibTracks[5][2][2], *h_etaMipTracks[5][2][2];
0266 TH1D *h_eHcal[5][6][48], *h_eCalo[5][6][48];
0267 TH1D *g_Pre, *g_PreL1, *g_PreHLT;
0268 TH1I *g_Accepts;
0269 std::vector<math::XYZTLorentzVector> vec_[3];
0270 };
0271
0272 IsoTrig::IsoTrig(const edm::ParameterSet &iConfig)
0273 : hltPrescaleProvider_(iConfig, consumesCollector(), *this),
0274 trigNames_(iConfig.getUntrackedParameter<std::vector<std::string>>("Triggers")),
0275 pixCandTag_(iConfig.getUntrackedParameter<edm::InputTag>("pixCandTag")),
0276 l1CandTag_(iConfig.getUntrackedParameter<edm::InputTag>("l1CandTag")),
0277 l2CandTag_(iConfig.getUntrackedParameter<edm::InputTag>("l2CandTag")),
0278 pixelTracksSources_(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("pixelTracksSources")),
0279 doL2L3_(iConfig.getUntrackedParameter<bool>("doL2L3", true)),
0280 doTiming_(iConfig.getUntrackedParameter<bool>("doTimingTree", true)),
0281 doMipCutTree_(iConfig.getUntrackedParameter<bool>("doMipCutTree", true)),
0282 doTrkResTree_(iConfig.getUntrackedParameter<bool>("doTrkResTree", true)),
0283 doChgIsolTree_(iConfig.getUntrackedParameter<bool>("doChgIsolTree", true)),
0284 doStudyIsol_(iConfig.getUntrackedParameter<bool>("doStudyIsol", true)),
0285 verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0286 pixelIsolationConeSizeAtEC_(iConfig.getUntrackedParameter<std::vector<double>>("pixelIsolationConeSizeAtEC")),
0287 minPTrackValue_(iConfig.getUntrackedParameter<double>("minPTrackValue")),
0288 vtxCutSeed_(iConfig.getUntrackedParameter<double>("vertexCutSeed")),
0289 vtxCutIsol_(iConfig.getUntrackedParameter<double>("vertexCutIsol")),
0290 tauUnbiasCone_(iConfig.getUntrackedParameter<double>("tauUnbiasCone")),
0291 prelimCone_(iConfig.getUntrackedParameter<double>("prelimCone")),
0292 theTrackQuality_(iConfig.getUntrackedParameter<std::string>("trackQuality", "highPurity")),
0293 processName_(iConfig.getUntrackedParameter<std::string>("processName", "HLT")),
0294 dr_L1_(iConfig.getUntrackedParameter<double>("isolationL1", 1.0)),
0295 a_coneR_(iConfig.getUntrackedParameter<double>("coneRadius", 34.98)),
0296 a_charIsoR_(a_coneR_ + 28.9),
0297 a_neutIsoR_(a_charIsoR_ * 0.726),
0298 a_mipR_(iConfig.getUntrackedParameter<double>("coneRadiusMIP", 14.0)),
0299 a_neutR1_(iConfig.getUntrackedParameter<double>("coneRadiusNeut1", 21.0)),
0300 a_neutR2_(iConfig.getUntrackedParameter<double>("coneRadiusNeut2", 29.0)),
0301 cutMip_(iConfig.getUntrackedParameter<double>("cutMIP", 1.0)),
0302 cutCharge_(iConfig.getUntrackedParameter<double>("chargeIsolation", 2.0)),
0303 cutNeutral_(iConfig.getUntrackedParameter<double>("neutralIsolation", 2.0)),
0304 minRunNo_(iConfig.getUntrackedParameter<int>("minRun")),
0305 maxRunNo_(iConfig.getUntrackedParameter<int>("maxRun")),
0306 changed_(false),
0307 t_timeL2Prod(nullptr),
0308 t_nPixCand(nullptr),
0309 t_nPixSeed(nullptr),
0310 t_nGoodTk(nullptr),
0311 t_TrkhCone(nullptr),
0312 t_TrkP(nullptr),
0313 t_TrkselTkFlag(nullptr),
0314 t_TrkqltyFlag(nullptr),
0315 t_TrkMissFlag(nullptr),
0316 t_TrkPVFlag(nullptr),
0317 t_TrkNuIsolFlag(nullptr),
0318 t_PixcandP(nullptr),
0319 t_PixcandPt(nullptr),
0320 t_PixcandEta(nullptr),
0321 t_PixcandPhi(nullptr),
0322 t_PixcandMaxP(nullptr),
0323 t_PixTrkcandP(nullptr),
0324 t_PixTrkcandPt(nullptr),
0325 t_PixTrkcandEta(nullptr),
0326 t_PixTrkcandPhi(nullptr),
0327 t_PixTrkcandMaxP(nullptr),
0328 t_PixTrkcandselTk(nullptr),
0329 t_NFcandP(nullptr),
0330 t_NFcandPt(nullptr),
0331 t_NFcandEta(nullptr),
0332 t_NFcandPhi(nullptr),
0333 t_NFcandEmip(nullptr),
0334 t_NFTrkcandP(nullptr),
0335 t_NFTrkcandPt(nullptr),
0336 t_NFTrkcandEta(nullptr),
0337 t_NFTrkcandPhi(nullptr),
0338 t_NFTrkcandEmip(nullptr),
0339 t_NFTrkMinDR(nullptr),
0340 t_NFTrkMinDP1(nullptr),
0341 t_NFTrkselTkFlag(nullptr),
0342 t_NFTrkqltyFlag(nullptr),
0343 t_NFTrkMissFlag(nullptr),
0344 t_NFTrkPVFlag(nullptr),
0345 t_NFTrkPropFlag(nullptr),
0346 t_NFTrkChgIsoFlag(nullptr),
0347 t_NFTrkNeuIsoFlag(nullptr),
0348 t_NFTrkMipFlag(nullptr),
0349 t_ECone(nullptr) {
0350 usesResource(TFileService::kSharedResource);
0351
0352
0353 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0354 selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("minTrackPt", 10.0);
0355 selectionParameters_.minQuality = trackQuality_;
0356 selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("maxDxyPV", 0.2);
0357 selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("maxDzPV", 5.0);
0358 selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("maxChi2", 5.0);
0359 selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("maxDpOverP", 0.1);
0360 selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("minOuterHit", 4);
0361 selectionParameters_.minLayerCrossed = iConfig.getUntrackedParameter<int>("minLayerCrossed", 8);
0362 selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("maxInMiss", 0);
0363 selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("maxOutMiss", 0);
0364
0365
0366 tok_lumi_ = consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"));
0367 edm::InputTag triggerEvent_("hltTriggerSummaryAOD", "", processName_);
0368 tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
0369 edm::InputTag theTriggerResultsLabel("TriggerResults", "", processName_);
0370 tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel);
0371 tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
0372 tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0373 tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0374 tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0375 tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0376 tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0377 tok_pixtk_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(pixCandTag_);
0378 tok_l1cand_ = consumes<trigger::TriggerFilterObjectWithRefs>(l1CandTag_);
0379 tok_l2cand_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(l2CandTag_);
0380 if (doTiming_) {
0381 tok_verthb_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
0382 tok_verthe_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
0383 tok_hlt_ = consumes<trigger::TriggerFilterObjectWithRefs>(edm::InputTag("hltL1sL1SingleJet68"));
0384 tok_SeedingLayerHB_ = consumes<SeedingLayerSetsHits>(edm::InputTag("hltPixelLayerTripletsHITHB"));
0385 tok_SeedingLayerHE_ = consumes<SeedingLayerSetsHits>(edm::InputTag("hltPixelLayerTripletsHITHE"));
0386 tok_SiPixelRecHits_ = consumes<SiPixelRecHitCollection>(edm::InputTag("hltSiPixelRecHits"));
0387 }
0388 if (doChgIsolTree_) {
0389 for (unsigned int k = 0; k < pixelTracksSources_.size(); ++k) {
0390
0391
0392 tok_pixtks_.push_back(consumes<reco::TrackCollection>(pixelTracksSources_[k]));
0393 }
0394 }
0395 if (verbosity_ >= 0) {
0396 edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0397 << "\t minPt " << selectionParameters_.minPt << "\t theTrackQuality "
0398 << theTrackQuality_ << "\t minQuality " << selectionParameters_.minQuality
0399 << "\t maxDxyPV " << selectionParameters_.maxDxyPV << "\t maxDzPV "
0400 << selectionParameters_.maxDzPV << "\t maxChi2 " << selectionParameters_.maxChi2
0401 << "\t maxDpOverP " << selectionParameters_.maxDpOverP << "\t minOuterHit "
0402 << selectionParameters_.minOuterHit << "\t minLayerCrossed "
0403 << selectionParameters_.minLayerCrossed << "\t maxInMiss "
0404 << selectionParameters_.maxInMiss << "\t maxOutMiss "
0405 << selectionParameters_.maxOutMiss << "\t a_coneR " << a_coneR_ << "\t a_charIsoR "
0406 << a_charIsoR_ << "\t a_neutIsoR " << a_neutIsoR_ << "\t a_mipR " << a_mipR_
0407 << "\t a_neutR " << a_neutR1_ << ":" << a_neutR2_ << "\t cuts (MIP " << cutMip_
0408 << " : Charge " << cutCharge_ << " : Neutral " << cutNeutral_ << ")";
0409 edm::LogVerbatim("IsoTrack") << "Charge Isolation parameters:"
0410 << "\t minPTrackValue " << minPTrackValue_ << "\t vtxCutSeed " << vtxCutSeed_
0411 << "\t vtxCutIsol " << vtxCutIsol_ << "\t tauUnbiasCone " << tauUnbiasCone_
0412 << "\t prelimCone " << prelimCone_ << "\t pixelIsolationConeSizeAtEC";
0413 for (unsigned int k = 0; k < pixelIsolationConeSizeAtEC_.size(); ++k)
0414 edm::LogVerbatim("IsoTrack") << "[" << k << "] " << pixelIsolationConeSizeAtEC_[k];
0415 }
0416 double pl[] = {20, 30, 40, 60, 80, 120};
0417 for (int i = 0; i < 6; ++i)
0418 pLimits_[i] = pl[i];
0419 rEB_ = 123.8;
0420 zEE_ = 317.0;
0421
0422 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0423 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0424 }
0425
0426 IsoTrig::~IsoTrig() {
0427
0428
0429 if (t_timeL2Prod)
0430 delete t_timeL2Prod;
0431 if (t_nPixCand)
0432 delete t_nPixCand;
0433 if (t_nPixSeed)
0434 delete t_nPixSeed;
0435 if (t_nGoodTk)
0436 delete t_nGoodTk;
0437 if (t_TrkhCone)
0438 delete t_TrkhCone;
0439 if (t_TrkP)
0440 delete t_TrkP;
0441 if (t_TrkselTkFlag)
0442 delete t_TrkselTkFlag;
0443 if (t_TrkqltyFlag)
0444 delete t_TrkqltyFlag;
0445 if (t_TrkMissFlag)
0446 delete t_TrkMissFlag;
0447 if (t_TrkPVFlag)
0448 delete t_TrkPVFlag;
0449 if (t_TrkNuIsolFlag)
0450 delete t_TrkNuIsolFlag;
0451 if (t_PixcandP)
0452 delete t_PixcandP;
0453 if (t_PixcandPt)
0454 delete t_PixcandPt;
0455 if (t_PixcandEta)
0456 delete t_PixcandEta;
0457 if (t_PixcandPhi)
0458 delete t_PixcandPhi;
0459 if (t_PixcandMaxP)
0460 delete t_PixcandMaxP;
0461 if (t_PixTrkcandP)
0462 delete t_PixTrkcandP;
0463 if (t_PixTrkcandPt)
0464 delete t_PixTrkcandPt;
0465 if (t_PixTrkcandEta)
0466 delete t_PixTrkcandEta;
0467 if (t_PixTrkcandPhi)
0468 delete t_PixTrkcandPhi;
0469 if (t_PixTrkcandMaxP)
0470 delete t_PixTrkcandMaxP;
0471 if (t_PixTrkcandselTk)
0472 delete t_PixTrkcandselTk;
0473 if (t_NFcandP)
0474 delete t_NFcandP;
0475 if (t_NFcandPt)
0476 delete t_NFcandPt;
0477 if (t_NFcandEta)
0478 delete t_NFcandEta;
0479 if (t_NFcandPhi)
0480 delete t_NFcandPhi;
0481 if (t_NFcandEmip)
0482 delete t_NFcandEmip;
0483 if (t_NFTrkcandP)
0484 delete t_NFTrkcandP;
0485 if (t_NFTrkcandPt)
0486 delete t_NFTrkcandPt;
0487 if (t_NFTrkcandEta)
0488 delete t_NFTrkcandEta;
0489 if (t_NFTrkcandPhi)
0490 delete t_NFTrkcandPhi;
0491 if (t_NFTrkcandEmip)
0492 delete t_NFTrkcandEmip;
0493 if (t_NFTrkMinDR)
0494 delete t_NFTrkMinDR;
0495 if (t_NFTrkMinDP1)
0496 delete t_NFTrkMinDP1;
0497 if (t_NFTrkselTkFlag)
0498 delete t_NFTrkselTkFlag;
0499 if (t_NFTrkqltyFlag)
0500 delete t_NFTrkqltyFlag;
0501 if (t_NFTrkMissFlag)
0502 delete t_NFTrkMissFlag;
0503 if (t_NFTrkPVFlag)
0504 delete t_NFTrkPVFlag;
0505 if (t_NFTrkPropFlag)
0506 delete t_NFTrkPropFlag;
0507 if (t_NFTrkChgIsoFlag)
0508 delete t_NFTrkChgIsoFlag;
0509 if (t_NFTrkNeuIsoFlag)
0510 delete t_NFTrkNeuIsoFlag;
0511 if (t_NFTrkMipFlag)
0512 delete t_NFTrkMipFlag;
0513 if (t_ECone)
0514 delete t_ECone;
0515 }
0516
0517 void IsoTrig::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0518 std::vector<std::string> triggers = {"HLT_IsoTrackHB"};
0519 std::vector<edm::InputTag> tags = {edm::InputTag("hltHITPixelTracksHB"), edm::InputTag("hltHITPixelTracksHE")};
0520 std::vector<double> cones = {35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 63.9, 70.0};
0521 edm::ParameterSetDescription desc;
0522 desc.addUntracked<std::vector<std::string>>("Triggers", triggers);
0523 desc.addUntracked<edm::InputTag>("pixCandTag", edm::InputTag(" "));
0524 desc.addUntracked<edm::InputTag>("l1CandTag", edm::InputTag("hltL1sV0SingleJet60"));
0525 desc.addUntracked<edm::InputTag>("l2CandTag", edm::InputTag("isolEcalPixelTrackProd"));
0526 desc.addUntracked<bool>("doL2L3", false);
0527 desc.addUntracked<bool>("doTimingTree", false);
0528 desc.addUntracked<bool>("doMipCutTree", false);
0529 desc.addUntracked<bool>("doTrkResTree", true);
0530 desc.addUntracked<bool>("doChgIsolTree", false);
0531 desc.addUntracked<bool>("doStudyIsol", false);
0532 desc.addUntracked<int>("verbosity", 0);
0533 desc.addUntracked<std::string>("processName", "HLT");
0534 desc.addUntracked<std::string>("trackQuality", "highPurity");
0535 desc.addUntracked<double>("minTrackPt", 10.0);
0536 desc.addUntracked<double>("maxDxyPV", 0.02);
0537 desc.addUntracked<double>("maxDzPV", 0.02);
0538 desc.addUntracked<double>("maxChi2", 5.0);
0539 desc.addUntracked<double>("maxDpOverP", 0.1);
0540 desc.addUntracked<int>("minOuterHit", 4);
0541 desc.addUntracked<int>("minLayerCrossed", 8);
0542 desc.addUntracked<int>("maxInMiss", 0);
0543 desc.addUntracked<int>("maxOutMiss", 0);
0544 desc.addUntracked<double>("isolationL1", 1.0);
0545 desc.addUntracked<double>("coneRadius", 34.98);
0546 desc.addUntracked<double>("coneRadiusMIP", 14.0);
0547 desc.addUntracked<double>("coneRadiusNeut1", 21.0);
0548 desc.addUntracked<double>("coneRadiusNeut2", 29.0);
0549 desc.addUntracked<double>("cutMIP", 1.0);
0550 desc.addUntracked<double>("chargeIsolation", 2.0);
0551 desc.addUntracked<double>("neutralIsolation", 2.0);
0552 desc.addUntracked<int>("minRun", 190456);
0553 desc.addUntracked<int>("maxRun", 203002);
0554 desc.addUntracked<std::vector<edm::InputTag>>("pixelTracksSources", tags);
0555 desc.addUntracked<std::vector<double>>("pixelIsolationConeSizeAtEC", cones);
0556 desc.addUntracked<double>("minPTrackValue", 0.0);
0557 desc.addUntracked<double>("vertexCutSeed", 101.0);
0558 desc.addUntracked<double>("vertexCutIsol", 101.0);
0559 desc.addUntracked<double>("tauUnbiasCone", 1.2);
0560 desc.addUntracked<double>("prelimCone", 1.0);
0561 desc.add<unsigned int>("stageL1Trigger", 1);
0562 descriptions.add("isoTrigDefault", desc);
0563 }
0564
0565 void IsoTrig::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0566 if (verbosity_ % 10 > 0)
0567 edm::LogVerbatim("IsoTrack") << "Event starts====================================";
0568
0569 int RunNo = iEvent.id().run();
0570
0571 HLTConfigProvider const &hltConfig = hltPrescaleProvider_.hltConfigProvider();
0572
0573 bField_ = &iSetup.getData(tok_magField_);
0574 geo_ = &iSetup.getData(tok_geom_);
0575 GlobalVector BField = bField_->inTesla(GlobalPoint(0, 0, 0));
0576 bfVal_ = BField.mag();
0577
0578 trigger::TriggerEvent triggerEvent;
0579 edm::Handle<trigger::TriggerEvent> triggerEventHandle;
0580 iEvent.getByToken(tok_trigEvt_, triggerEventHandle);
0581 if (!triggerEventHandle.isValid()) {
0582 edm::LogWarning("IsoTrack") << "Error! Can't get the product hltTriggerSummaryAOD";
0583
0584 } else {
0585 triggerEvent = *(triggerEventHandle.product());
0586 }
0587 const trigger::TriggerObjectCollection &TOC(triggerEvent.getObjects());
0588
0589 edm::Handle<edm::TriggerResults> triggerResults;
0590 iEvent.getByToken(tok_trigRes_, triggerResults);
0591
0592 edm::Handle<reco::TrackCollection> trkCollection;
0593 iEvent.getByToken(tok_genTrack_, trkCollection);
0594
0595 iEvent.getByToken(tok_EB_, barrelRecHitsHandle_);
0596 iEvent.getByToken(tok_EE_, endcapRecHitsHandle_);
0597
0598 iEvent.getByToken(tok_hbhe_, hbhe_);
0599
0600 iEvent.getByToken(tok_recVtx_, recVtxs_);
0601 iEvent.getByToken(tok_bs_, beamSpotH_);
0602 if (!recVtxs_->empty() && !((*recVtxs_)[0].isFake())) {
0603 leadPV_ = math::XYZPoint((*recVtxs_)[0].x(), (*recVtxs_)[0].y(), (*recVtxs_)[0].z());
0604 } else if (beamSpotH_.isValid()) {
0605 leadPV_ = beamSpotH_->position();
0606 }
0607
0608 if ((verbosity_ / 100) % 10 > 0) {
0609 edm::LogVerbatim("IsoTrack") << "Primary Vertex " << leadPV_;
0610 if (beamSpotH_.isValid())
0611 edm::LogVerbatim("IsoTrack") << "Beam Spot " << beamSpotH_->position();
0612 }
0613 pixelTrackRefsHE_.clear();
0614 pixelTrackRefsHB_.clear();
0615 for (unsigned int iPix = 0; iPix < pixelTracksSources_.size(); iPix++) {
0616 edm::Handle<reco::TrackCollection> iPixCol;
0617 iEvent.getByToken(tok_pixtks_[iPix], iPixCol);
0618 if (iPixCol.isValid()) {
0619 for (reco::TrackCollection::const_iterator pit = iPixCol->begin(); pit != iPixCol->end(); pit++) {
0620 if (iPix == 0)
0621 pixelTrackRefsHB_.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0622 pixelTrackRefsHE_.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0623 }
0624 }
0625 }
0626 if (doTiming_)
0627 getGoodTracks(iEvent, trkCollection);
0628
0629 for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0630 std::string FilterNames[7] = {"hltL1sL1SingleJet68",
0631 "hltIsolPixelTrackL2FilterHE",
0632 "ecalIsolPixelTrackFilterHE",
0633 "hltIsolPixelTrackL3FilterHE",
0634 "hltIsolPixelTrackL2FilterHB",
0635 "ecalIsolPixelTrackFilterHB",
0636 "hltIsolPixelTrackL3FilterHB"};
0637 std::string label = triggerEvent.filterTag(ifilter).label();
0638 for (int i = 0; i < 7; i++) {
0639 if (label == FilterNames[i])
0640 h_Filters->Fill(i);
0641 }
0642 }
0643 edm::InputTag lumiProducer("LumiProducer", "", "RECO");
0644 edm::Handle<LumiDetails> Lumid;
0645 iEvent.getLuminosityBlock().getByToken(tok_lumi_, Lumid);
0646 float mybxlumi = -1;
0647 if (Lumid.isValid())
0648 mybxlumi = Lumid->lumiValue(LumiDetails::kOCC1, iEvent.bunchCrossing()) * 6.37;
0649 if (verbosity_ % 10 > 0)
0650 edm::LogVerbatim("IsoTrack") << "RunNo " << RunNo << " EvtNo " << iEvent.id().event() << " Lumi "
0651 << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing() << " mybxlumi "
0652 << mybxlumi;
0653 if (!triggerResults.isValid()) {
0654 edm::LogWarning("IsoTrack") << "Error! Can't get the product triggerResults";
0655
0656
0657 } else {
0658 std::vector<std::string> modules;
0659 h_nHLT->Fill(triggerResults->size());
0660 const edm::TriggerNames &triggerNames = iEvent.triggerNames(*triggerResults);
0661
0662 const std::vector<std::string> &triggerNames_ = triggerNames.triggerNames();
0663 if (verbosity_ % 10 > 1)
0664 edm::LogVerbatim("IsoTrack") << "number of HLTs " << triggerNames_.size();
0665 double preL1(-1), preHLT(-1), prescale(-1);
0666 int hlt(-1);
0667 for (unsigned int i = 0; i < triggerResults->size(); i++) {
0668 unsigned int triggerindx = hltConfig.triggerIndex(triggerNames_[i]);
0669 const std::vector<std::string> &moduleLabels(hltConfig.moduleLabels(triggerindx));
0670
0671 for (unsigned int in = 0; in < trigNames_.size(); ++in) {
0672
0673 if (triggerNames_[i].find(trigNames_[in]) != std::string::npos) {
0674 if (verbosity_ % 10 > 0)
0675 edm::LogVerbatim("IsoTrack") << "trigger that i want " << triggerNames_[i] << " accept "
0676 << triggerResults->accept(i);
0677 hlt = triggerResults->accept(i);
0678 h_HLT->Fill(hlt);
0679
0680 if (hlt > 0) {
0681 edm::Handle<reco::IsolatedPixelTrackCandidateCollection> Pixcands;
0682 iEvent.getByToken(tok_pixtk_, Pixcands);
0683 edm::Handle<trigger::TriggerFilterObjectWithRefs> L1cands;
0684 iEvent.getByToken(tok_l1cand_, L1cands);
0685
0686 auto const prescales = hltPrescaleProvider_.prescaleValues<double>(iEvent, iSetup, triggerNames_[i]);
0687 preL1 = prescales.first;
0688 preHLT = prescales.second;
0689 prescale = preL1 * preHLT;
0690 if (verbosity_ % 10 > 0)
0691 edm::LogVerbatim("IsoTrack")
0692 << triggerNames_[i] << " accept " << hlt << " preL1 " << preL1 << " preHLT " << preHLT;
0693 for (int iv = 0; iv < 3; ++iv)
0694 vec_[iv].clear();
0695 if (trigList_.find(RunNo) != trigList_.end()) {
0696 trigList_[RunNo] += 1;
0697 } else {
0698 trigList_.insert({RunNo, 1});
0699 trigPreList_.insert({RunNo, prescales});
0700 }
0701
0702 for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0703 std::vector<int> Keys;
0704 std::string label = triggerEvent.filterTag(ifilter).label();
0705
0706 for (unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
0707 if (label.find(moduleLabels[imodule]) != std::string::npos) {
0708 if (verbosity_ % 10 > 0)
0709 edm::LogVerbatim("IsoTrack") << "FILTERNAME " << label;
0710 for (unsigned int ifiltrKey = 0; ifiltrKey < triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
0711 Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
0712 const trigger::TriggerObject &TO(TOC[Keys[ifiltrKey]]);
0713 math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
0714 if (label.find("L2Filter") != std::string::npos) {
0715 vec_[1].push_back(v4);
0716 } else if (label.find("L3Filter") != std::string::npos) {
0717 vec_[2].push_back(v4);
0718 } else {
0719 vec_[0].push_back(v4);
0720 h_L1ObjEnergy->Fill(TO.energy());
0721 }
0722 if (verbosity_ % 10 > 0)
0723 edm::LogVerbatim("IsoTrack") << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta()
0724 << " phi " << TO.phi() << " mass " << TO.mass() << " Id " << TO.id();
0725 }
0726 }
0727 }
0728 }
0729 std::vector<reco::TrackCollection::const_iterator> goodTks;
0730 if (doL2L3_) {
0731 h_nL3Objs->Fill(vec_[2].size());
0732 studyTrigger(trkCollection, goodTks);
0733 } else {
0734 if (trkCollection.isValid()) {
0735 reco::TrackCollection::const_iterator trkItr;
0736 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++)
0737 goodTks.push_back(trkItr);
0738 }
0739 }
0740
0741 if (doStudyIsol_)
0742 studyIsolation(trkCollection, goodTks);
0743 if (doTrkResTree_)
0744 StudyTrkEbyP(trkCollection);
0745
0746 std::pair<double, double> etaphi = etaPhiTrigger();
0747 edm::Handle<reco::IsolatedPixelTrackCandidateCollection> L2cands;
0748 iEvent.getByToken(tok_l2cand_, L2cands);
0749 if (!L2cands.isValid()) {
0750 if (verbosity_ % 10 > 0)
0751 edm::LogVerbatim("IsoTrack") << " trigCand is not valid ";
0752 } else {
0753 if (doMipCutTree_)
0754 studyMipCut(trkCollection, L2cands);
0755 }
0756 if (!pixelTracksSources_.empty())
0757 if (doChgIsolTree_ && !pixelTrackRefsHE_.empty())
0758 chgIsolation(etaphi.first, etaphi.second, trkCollection, iEvent);
0759 }
0760 break;
0761 }
0762 }
0763 }
0764 h_PreL1->Fill(preL1);
0765 h_PreHLT->Fill(preHLT);
0766 h_Pre->Fill(prescale);
0767 h_PreL1wt->Fill(preL1, mybxlumi);
0768 h_PreHLTwt->Fill(preHLT, mybxlumi);
0769
0770
0771
0772 if (changed_) {
0773 changed_ = false;
0774 if ((verbosity_ / 10) % 10 > 1) {
0775 edm::LogVerbatim("IsoTrack") << "New trigger menu found !!!";
0776 const unsigned int n(hltConfig.size());
0777 for (unsigned itrig = 0; itrig < triggerNames_.size(); itrig++) {
0778 unsigned int triggerindx = hltConfig.triggerIndex(triggerNames_[itrig]);
0779 if (triggerindx >= n)
0780 edm::LogVerbatim("IsoTrack") << triggerNames_[itrig] << " " << triggerindx << " does not exist in"
0781 << " the current menu";
0782 else
0783 edm::LogVerbatim("IsoTrack") << triggerNames_[itrig] << " " << triggerindx << " exists";
0784 }
0785 }
0786 }
0787 }
0788 if (doTiming_)
0789 studyTiming(iEvent);
0790 }
0791
0792 void IsoTrig::clearChgIsolnTreeVectors() {
0793 t_PixcandP->clear();
0794 t_PixcandPt->clear();
0795 t_PixcandEta->clear();
0796 t_PixcandPhi->clear();
0797 for (unsigned int i = 0; i < t_PixcandMaxP->size(); i++)
0798 t_PixcandMaxP[i].clear();
0799 t_PixcandMaxP->clear();
0800 t_PixTrkcandP->clear();
0801 t_PixTrkcandPt->clear();
0802 t_PixTrkcandEta->clear();
0803 t_PixTrkcandPhi->clear();
0804 t_PixTrkcandMaxP->clear();
0805 t_PixTrkcandselTk->clear();
0806 }
0807
0808 void IsoTrig::clearMipCutTreeVectors() {
0809 t_NFcandP->clear();
0810 t_NFcandPt->clear();
0811 t_NFcandEta->clear();
0812 t_NFcandPhi->clear();
0813 t_NFcandEmip->clear();
0814 t_NFTrkcandP->clear();
0815 t_NFTrkcandPt->clear();
0816 t_NFTrkcandEta->clear();
0817 t_NFTrkcandPhi->clear();
0818 t_NFTrkcandEmip->clear();
0819 t_NFTrkMinDR->clear();
0820 t_NFTrkMinDP1->clear();
0821 t_NFTrkselTkFlag->clear();
0822 t_NFTrkqltyFlag->clear();
0823 t_NFTrkMissFlag->clear();
0824 t_NFTrkPVFlag->clear();
0825 t_NFTrkPropFlag->clear();
0826 t_NFTrkChgIsoFlag->clear();
0827 t_NFTrkNeuIsoFlag->clear();
0828 t_NFTrkMipFlag->clear();
0829 t_ECone->clear();
0830 }
0831
0832 void IsoTrig::pushChgIsolnTreeVecs(math::XYZTLorentzVector &Pixcand,
0833 math::XYZTLorentzVector &Trkcand,
0834 std::vector<double> &PixMaxP,
0835 double &TrkMaxP,
0836 bool &selTk) {
0837 t_PixcandP->push_back(Pixcand.r());
0838 t_PixcandPt->push_back(Pixcand.pt());
0839 t_PixcandEta->push_back(Pixcand.eta());
0840 t_PixcandPhi->push_back(Pixcand.phi());
0841 t_PixcandMaxP->push_back(PixMaxP);
0842 t_PixTrkcandP->push_back(Trkcand.r());
0843 t_PixTrkcandPt->push_back(Trkcand.pt());
0844 t_PixTrkcandEta->push_back(Trkcand.eta());
0845 t_PixTrkcandPhi->push_back(Trkcand.phi());
0846 t_PixTrkcandMaxP->push_back(TrkMaxP);
0847 t_PixTrkcandselTk->push_back(selTk);
0848 }
0849
0850 void IsoTrig::pushMipCutTreeVecs(math::XYZTLorentzVector &NFcand,
0851 math::XYZTLorentzVector &Trkcand,
0852 double &EmipNFcand,
0853 double &EmipTrkcand,
0854 double &mindR,
0855 double &mindP1,
0856 std::vector<bool> &Flags,
0857 double hCone) {
0858 t_NFcandP->push_back(NFcand.r());
0859 t_NFcandPt->push_back(NFcand.pt());
0860 t_NFcandEta->push_back(NFcand.eta());
0861 t_NFcandPhi->push_back(NFcand.phi());
0862 t_NFcandEmip->push_back(EmipNFcand);
0863 t_NFTrkcandP->push_back(Trkcand.r());
0864 t_NFTrkcandPt->push_back(Trkcand.pt());
0865 t_NFTrkcandEta->push_back(Trkcand.eta());
0866 t_NFTrkcandPhi->push_back(Trkcand.phi());
0867 t_NFTrkcandEmip->push_back(EmipTrkcand);
0868 t_NFTrkMinDR->push_back(mindR);
0869 t_NFTrkMinDP1->push_back(mindP1);
0870 t_NFTrkselTkFlag->push_back(Flags[0]);
0871 t_NFTrkqltyFlag->push_back(Flags[1]);
0872 t_NFTrkMissFlag->push_back(Flags[2]);
0873 t_NFTrkPVFlag->push_back(Flags[3]);
0874 t_NFTrkPropFlag->push_back(Flags[4]);
0875 t_NFTrkChgIsoFlag->push_back(Flags[5]);
0876 t_NFTrkNeuIsoFlag->push_back(Flags[6]);
0877 t_NFTrkMipFlag->push_back(Flags[7]);
0878 t_ECone->push_back(hCone);
0879 }
0880
0881 void IsoTrig::beginJob() {
0882 char hname[100], htit[100];
0883 std::string levels[20] = {"L1", "L2", "L3", "Reco", "RecoMatch", "RecoNoMatch",
0884 "L2Match", "L2NoMatch", "L3Match", "L3NoMatch", "HLTTrk", "HLTGoodTrk",
0885 "HLTIsoTrk", "HLTMip", "HLTSelect", "nonHLTTrk", "nonHLTGoodTrk", "nonHLTIsoTrk",
0886 "nonHLTMip", "nonHLTSelect"};
0887 if (doTiming_) {
0888 TimingTree_ = fs_->make<TTree>("TimingTree", "TimingTree");
0889 t_timeL2Prod = new std::vector<double>();
0890 t_nPixCand = new std::vector<int>();
0891 t_nPixSeed = new std::vector<int>();
0892 t_nGoodTk = new std::vector<int>();
0893
0894 TimingTree_->Branch("t_timeL2Prod", "std::vector<double>", &t_timeL2Prod);
0895 TimingTree_->Branch("t_nPixCand", "std::vector<int>", &t_nPixCand);
0896 TimingTree_->Branch("t_nPixSeed", "std::vector<int>", &t_nPixSeed);
0897 TimingTree_->Branch("t_nGoodTk", "std::vector<int>", &t_nGoodTk);
0898 }
0899 if (doTrkResTree_) {
0900 TrkResTree_ = fs_->make<TTree>("TrkRestree", "TrkResTree");
0901 t_TrkhCone = new std::vector<double>();
0902 t_TrkP = new std::vector<double>();
0903 t_TrkselTkFlag = new std::vector<bool>();
0904 t_TrkqltyFlag = new std::vector<bool>();
0905 t_TrkMissFlag = new std::vector<bool>();
0906 t_TrkPVFlag = new std::vector<bool>();
0907 t_TrkNuIsolFlag = new std::vector<bool>();
0908
0909 TrkResTree_->Branch("t_TrkhCone", "std::vector<double>", &t_TrkhCone);
0910 TrkResTree_->Branch("t_TrkP", "std::vector<double>", &t_TrkP);
0911 TrkResTree_->Branch("t_TrkselTkFlag", "std::vector<bool>", &t_TrkselTkFlag);
0912 TrkResTree_->Branch("t_TrkqltyFlag", "std::vector<bool>", &t_TrkqltyFlag);
0913 TrkResTree_->Branch("t_TrkMissFlag", "std::vector<bool>", &t_TrkMissFlag);
0914 TrkResTree_->Branch("t_TrkPVFlag", "std::vector<bool>", &t_TrkPVFlag);
0915 TrkResTree_->Branch("t_TrkNuIsolFlag", "std::vector<bool>", &t_TrkNuIsolFlag);
0916 }
0917 if (doChgIsolTree_) {
0918 ChgIsolnTree_ = fs_->make<TTree>("ChgIsolnTree", "ChgIsolnTree");
0919 t_PixcandP = new std::vector<double>();
0920 t_PixcandPt = new std::vector<double>();
0921 t_PixcandEta = new std::vector<double>();
0922 t_PixcandPhi = new std::vector<double>();
0923 t_PixcandMaxP = new std::vector<std::vector<double>>();
0924 t_PixTrkcandP = new std::vector<double>();
0925 t_PixTrkcandPt = new std::vector<double>();
0926 t_PixTrkcandEta = new std::vector<double>();
0927 t_PixTrkcandPhi = new std::vector<double>();
0928 t_PixTrkcandMaxP = new std::vector<double>();
0929 t_PixTrkcandselTk = new std::vector<bool>();
0930
0931 ChgIsolnTree_->Branch("t_PixcandP", "std::vector<double>", &t_PixcandP);
0932 ChgIsolnTree_->Branch("t_PixcandPt", "std::vector<double>", &t_PixcandPt);
0933 ChgIsolnTree_->Branch("t_PixcandEta", "std::vector<double>", &t_PixcandEta);
0934 ChgIsolnTree_->Branch("t_PixcandPhi", "std::vector<double>", &t_PixcandPhi);
0935 ChgIsolnTree_->Branch("t_PixcandMaxP", "std::vector<std::vector<double> >", &t_PixcandMaxP);
0936 ChgIsolnTree_->Branch("t_PixTrkcandP", "std::vector<double>", &t_PixTrkcandP);
0937 ChgIsolnTree_->Branch("t_PixTrkcandPt", "std::vector<double>", &t_PixTrkcandPt);
0938 ChgIsolnTree_->Branch("t_PixTrkcandEta", "std::vector<double>", &t_PixTrkcandEta);
0939 ChgIsolnTree_->Branch("t_PixTrkcandPhi", "std::vector<double>", &t_PixTrkcandPhi);
0940 ChgIsolnTree_->Branch("t_PixTrkcandMaxP", "std::vector<double>", &t_PixTrkcandMaxP);
0941 ChgIsolnTree_->Branch("t_PixTrkcandselTk", "std::vector<bool>", &t_PixTrkcandselTk);
0942 }
0943 if (doMipCutTree_) {
0944 MipCutTree_ = fs_->make<TTree>("MipCutTree", "MipCutTree");
0945 t_NFcandP = new std::vector<double>();
0946 t_NFcandPt = new std::vector<double>();
0947 t_NFcandEta = new std::vector<double>();
0948 t_NFcandPhi = new std::vector<double>();
0949 t_NFcandEmip = new std::vector<double>();
0950 t_NFTrkcandP = new std::vector<double>();
0951 t_NFTrkcandPt = new std::vector<double>();
0952 t_NFTrkcandEta = new std::vector<double>();
0953 t_NFTrkcandPhi = new std::vector<double>();
0954 t_NFTrkcandEmip = new std::vector<double>();
0955 t_NFTrkMinDR = new std::vector<double>();
0956 t_NFTrkMinDP1 = new std::vector<double>();
0957 t_NFTrkselTkFlag = new std::vector<bool>();
0958 t_NFTrkqltyFlag = new std::vector<bool>();
0959 t_NFTrkMissFlag = new std::vector<bool>();
0960 t_NFTrkPVFlag = new std::vector<bool>();
0961 t_NFTrkPropFlag = new std::vector<bool>();
0962 t_NFTrkChgIsoFlag = new std::vector<bool>();
0963 t_NFTrkNeuIsoFlag = new std::vector<bool>();
0964 t_NFTrkMipFlag = new std::vector<bool>();
0965 t_ECone = new std::vector<double>();
0966
0967 MipCutTree_->Branch("t_NFcandP", "std::vector<double>", &t_NFcandP);
0968 MipCutTree_->Branch("t_NFcandPt", "std::vector<double>", &t_NFcandPt);
0969 MipCutTree_->Branch("t_NFcandEta", "std::vector<double>", &t_NFcandEta);
0970 MipCutTree_->Branch("t_NFcandPhi", "std::vector<double>", &t_NFcandPhi);
0971 MipCutTree_->Branch("t_NFcandEmip", "std::vector<double>", &t_NFcandEmip);
0972 MipCutTree_->Branch("t_NFTrkcandP", "std::vector<double>", &t_NFTrkcandP);
0973 MipCutTree_->Branch("t_NFTrkcandPt", "std::vector<double>", &t_NFTrkcandPt);
0974 MipCutTree_->Branch("t_NFTrkcandEta", "std::vector<double>", &t_NFTrkcandEta);
0975 MipCutTree_->Branch("t_NFTrkcandPhi", "std::vector<double>", &t_NFTrkcandPhi);
0976 MipCutTree_->Branch("t_NFTrkcandEmip", "std::vector<double>", &t_NFTrkcandEmip);
0977 MipCutTree_->Branch("t_NFTrkMinDR", "std::vector<double>", &t_NFTrkMinDR);
0978 MipCutTree_->Branch("t_NFTrkMinDP1", "std::vector<double>", &t_NFTrkMinDP1);
0979 MipCutTree_->Branch("t_NFTrkselTkFlag", "std::vector<bool>", &t_NFTrkselTkFlag);
0980 MipCutTree_->Branch("t_NFTrkqltyFlag", "std::vector<bool>", &t_NFTrkqltyFlag);
0981 MipCutTree_->Branch("t_NFTrkMissFlag", "std::vector<bool>", &t_NFTrkMissFlag);
0982 MipCutTree_->Branch("t_NFTrkPVFlag", "std::vector<bool>", &t_NFTrkPVFlag);
0983 MipCutTree_->Branch("t_NFTrkPropFlag", "std::vector<bool>", &t_NFTrkPropFlag);
0984 MipCutTree_->Branch("t_NFTrkChgIsoFlag", "std::vector<bool>", &t_NFTrkChgIsoFlag);
0985 MipCutTree_->Branch("t_NFTrkNeuIsoFlag", "std::vector<bool>", &t_NFTrkNeuIsoFlag);
0986 MipCutTree_->Branch("t_NFTrkMipFlag", "std::vector<bool>", &t_NFTrkMipFlag);
0987 MipCutTree_->Branch("t_ECone", "std::vector<double>", &t_ECone);
0988 }
0989 h_Filters = fs_->make<TH1I>("h_Filters", "Filter Accepts", 10, 0, 10);
0990 std::string FilterNames[7] = {"hltL1sL1SingleJet68",
0991 "hltIsolPixelTrackL2FilterHE",
0992 "ecalIsolPixelTrackFilterHE",
0993 "hltIsolPixelTrackL3FilterHE",
0994 "hltIsolPixelTrackL2FilterHB",
0995 "ecalIsolPixelTrackFilterHB",
0996 "hltIsolPixelTrackL3FilterHB"};
0997 for (int i = 0; i < 7; i++) {
0998 h_Filters->GetXaxis()->SetBinLabel(i + 1, FilterNames[i].c_str());
0999 }
1000
1001 h_nHLT = fs_->make<TH1I>("h_nHLT", "Size of trigger Names", 1000, 1, 1000);
1002 h_HLT = fs_->make<TH1I>("h_HLT", "HLT accept", 3, -1, 2);
1003 h_PreL1 = fs_->make<TH1I>("h_PreL1", "L1 Prescale", 500, 0, 500);
1004 h_PreHLT = fs_->make<TH1I>("h_PreHLT", "HLT Prescale", 50, 0, 50);
1005 h_Pre = fs_->make<TH1I>("h_Pre", "Prescale", 3000, 0, 3000);
1006
1007 h_PreL1wt = fs_->make<TH1D>("h_PreL1wt", "Weighted L1 Prescale", 500, 0, 500);
1008 h_PreHLTwt = fs_->make<TH1D>("h_PreHLTwt", "Weighted HLT Prescale", 500, 0, 100);
1009 h_L1ObjEnergy = fs_->make<TH1D>("h_L1ObjEnergy", "Energy of L1Object", 500, 0.0, 500.0);
1010
1011 h_EnIn = fs_->make<TH1D>("h_EnInEcal", "EnergyIn Ecal", 200, 0.0, 20.0);
1012 h_EnOut = fs_->make<TH1D>("h_EnOutEcal", "EnergyOut Ecal", 200, 0.0, 20.0);
1013 h_MipEnMatch =
1014 fs_->make<TH2D>("h_MipEnMatch", "MipEn: HLT level vs Reco Level (Matched)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1015 h_MipEnNoMatch = fs_->make<TH2D>(
1016 "h_MipEnNoMatch", "MipEn: HLT level vs Reco Level (No Match Found)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1017
1018 if (doL2L3_) {
1019 h_nL3Objs = fs_->make<TH1I>("h_nL3Objs", "Number of L3 objects", 10, 0, 10);
1020
1021 std::string pairs[9] = {"L2L3",
1022 "L2L3Match",
1023 "L2L3NoMatch",
1024 "L3Reco",
1025 "L3RecoMatch",
1026 "L3RecoNoMatch",
1027 "NewFilterReco",
1028 "NewFilterRecoMatch",
1029 "NewFilterRecoNoMatch"};
1030 for (int ipair = 0; ipair < 9; ipair++) {
1031 sprintf(hname, "h_dEta%s", pairs[ipair].c_str());
1032 sprintf(htit, "#Delta#eta for %s", pairs[ipair].c_str());
1033 h_dEta[ipair] = fs_->make<TH1D>(hname, htit, 200, -10.0, 10.0);
1034 h_dEta[ipair]->GetXaxis()->SetTitle("d#eta");
1035
1036 sprintf(hname, "h_dPhi%s", pairs[ipair].c_str());
1037 sprintf(htit, "#Delta#phi for %s", pairs[ipair].c_str());
1038 h_dPhi[ipair] = fs_->make<TH1D>(hname, htit, 140, -7.0, 7.0);
1039 h_dPhi[ipair]->GetXaxis()->SetTitle("d#phi");
1040
1041 sprintf(hname, "h_dPt%s", pairs[ipair].c_str());
1042 sprintf(htit, "#Delta dp_{T} for %s objects", pairs[ipair].c_str());
1043 h_dPt[ipair] = fs_->make<TH1D>(hname, htit, 400, -200.0, 200.0);
1044 h_dPt[ipair]->GetXaxis()->SetTitle("dp_{T} (GeV)");
1045
1046 sprintf(hname, "h_dP%s", pairs[ipair].c_str());
1047 sprintf(htit, "#Delta p for %s objects", pairs[ipair].c_str());
1048 h_dP[ipair] = fs_->make<TH1D>(hname, htit, 400, -200.0, 200.0);
1049 h_dP[ipair]->GetXaxis()->SetTitle("dP (GeV)");
1050
1051 sprintf(hname, "h_dinvPt%s", pairs[ipair].c_str());
1052 sprintf(htit, "#Delta (1/p_{T}) for %s objects", pairs[ipair].c_str());
1053 h_dinvPt[ipair] = fs_->make<TH1D>(hname, htit, 500, -0.4, 0.1);
1054 h_dinvPt[ipair]->GetXaxis()->SetTitle("d(1/p_{T})");
1055 sprintf(hname, "h_mindR%s", pairs[ipair].c_str());
1056 sprintf(htit, "min(#Delta R) for %s objects", pairs[ipair].c_str());
1057 h_mindR[ipair] = fs_->make<TH1D>(hname, htit, 500, 0.0, 1.0);
1058 h_mindR[ipair]->GetXaxis()->SetTitle("dR");
1059 }
1060
1061 for (int lvl = 0; lvl < 2; lvl++) {
1062 sprintf(hname, "h_dEtaL1%s", levels[lvl + 1].c_str());
1063 sprintf(htit, "#Delta#eta for L1 and %s objects", levels[lvl + 1].c_str());
1064 h_dEtaL1[lvl] = fs_->make<TH1D>(hname, htit, 400, -10.0, 10.0);
1065
1066 sprintf(hname, "h_dPhiL1%s", levels[lvl + 1].c_str());
1067 sprintf(htit, "#Delta#phi for L1 and %s objects", levels[lvl + 1].c_str());
1068 h_dPhiL1[lvl] = fs_->make<TH1D>(hname, htit, 280, -7.0, 7.0);
1069
1070 sprintf(hname, "h_dRL1%s", levels[lvl + 1].c_str());
1071 sprintf(htit, "#Delta R for L1 and %s objects", levels[lvl + 1].c_str());
1072 h_dRL1[lvl] = fs_->make<TH1D>(hname, htit, 100, 0.0, 10.0);
1073 }
1074 }
1075
1076 int levmin = (doL2L3_ ? 0 : 10);
1077 for (int ilevel = levmin; ilevel < 20; ilevel++) {
1078 sprintf(hname, "h_p%s", levels[ilevel].c_str());
1079 sprintf(htit, "p for %s objects", levels[ilevel].c_str());
1080 h_p[ilevel] = fs_->make<TH1D>(hname, htit, 100, 0.0, 500.0);
1081 h_p[ilevel]->GetXaxis()->SetTitle("p (GeV)");
1082
1083 sprintf(hname, "h_pt%s", levels[ilevel].c_str());
1084 sprintf(htit, "p_{T} for %s objects", levels[ilevel].c_str());
1085 h_pt[ilevel] = fs_->make<TH1D>(hname, htit, 100, 0.0, 500.0);
1086 h_pt[ilevel]->GetXaxis()->SetTitle("p_{T} (GeV)");
1087
1088 sprintf(hname, "h_eta%s", levels[ilevel].c_str());
1089 sprintf(htit, "#eta for %s objects", levels[ilevel].c_str());
1090 h_eta[ilevel] = fs_->make<TH1D>(hname, htit, 100, -5.0, 5.0);
1091 h_eta[ilevel]->GetXaxis()->SetTitle("#eta");
1092
1093 sprintf(hname, "h_phi%s", levels[ilevel].c_str());
1094 sprintf(htit, "#phi for %s objects", levels[ilevel].c_str());
1095 h_phi[ilevel] = fs_->make<TH1D>(hname, htit, 70, -3.5, 3.50);
1096 h_phi[ilevel]->GetXaxis()->SetTitle("#phi");
1097 }
1098
1099 std::string cuts[2] = {"HLTMatched", "HLTNotMatched"};
1100 std::string cuts2[2] = {"All", "Away from L1"};
1101 for (int icut = 0; icut < 2; icut++) {
1102 sprintf(hname, "h_eMip%s", cuts[icut].c_str());
1103 sprintf(htit, "eMip for %s tracks", cuts[icut].c_str());
1104 h_eMip[icut] = fs_->make<TH1D>(hname, htit, 200, 0.0, 10.0);
1105 h_eMip[icut]->GetXaxis()->SetTitle("E_{Mip} (GeV)");
1106
1107 sprintf(hname, "h_eMaxNearP%s", cuts[icut].c_str());
1108 sprintf(htit, "eMaxNearP for %s tracks", cuts[icut].c_str());
1109 h_eMaxNearP[icut] = fs_->make<TH1D>(hname, htit, 240, -2.0, 10.0);
1110 h_eMaxNearP[icut]->GetXaxis()->SetTitle("E_{MaxNearP} (GeV)");
1111
1112 sprintf(hname, "h_eNeutIso%s", cuts[icut].c_str());
1113 sprintf(htit, "eNeutIso for %s ", cuts[icut].c_str());
1114 h_eNeutIso[icut] = fs_->make<TH1D>(hname, htit, 200, 0.0, 10.0);
1115 h_eNeutIso[icut]->GetXaxis()->SetTitle("E_{NeutIso} (GeV)");
1116
1117 for (int kcut = 0; kcut < 2; ++kcut) {
1118 for (int lim = 0; lim < 5; ++lim) {
1119 sprintf(hname, "h_etaCalibTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
1120 sprintf(htit,
1121 "#eta for %s isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1122 cuts[icut].c_str(),
1123 pLimits_[lim],
1124 pLimits_[lim + 1],
1125 cuts2[kcut].c_str());
1126 h_etaCalibTracks[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 60, -30.0, 30.0);
1127 h_etaCalibTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
1128
1129 sprintf(hname, "h_etaMipTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
1130 sprintf(htit,
1131 "#eta for %s charge isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1132 cuts[icut].c_str(),
1133 pLimits_[lim],
1134 pLimits_[lim + 1],
1135 cuts2[kcut].c_str());
1136 h_etaMipTracks[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 60, -30.0, 30.0);
1137 h_etaMipTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
1138 }
1139 }
1140 }
1141
1142 std::string ecut1[3] = {"all", "HLTMatched", "HLTNotMatched"};
1143 std::string ecut2[2] = {"without", "with"};
1144 int etac[48] = {-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, -13, -14, -15, -16,
1145 -17, -18, -19, -20, -21, -22, -23, -24, 1, 2, 3, 4, 5, 6, 7, 8,
1146 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24};
1147 for (int icut = 0; icut < 6; icut++) {
1148
1149 int i1 = (icut > 2 ? 1 : 0);
1150 int i2 = icut - i1 * 3;
1151 for (int kcut = 0; kcut < 48; kcut++) {
1152 for (int lim = 0; lim < 5; ++lim) {
1153 sprintf(hname, "h_eta%dEnHcal%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
1154 sprintf(htit,
1155 "HCAL energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1156 etac[kcut],
1157 ecut1[i2].c_str(),
1158 pLimits_[lim],
1159 pLimits_[lim + 1],
1160 ecut2[i1].c_str());
1161 h_eHcal[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 750, 0.0, 150.0);
1162 h_eHcal[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
1163 sprintf(hname, "h_eta%dEnCalo%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
1164 sprintf(htit,
1165 "Calorimter energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1166 etac[kcut],
1167 ecut1[i2].c_str(),
1168 pLimits_[lim],
1169 pLimits_[lim + 1],
1170 ecut2[i1].c_str());
1171 h_eCalo[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 750, 0.0, 150.0);
1172 h_eCalo[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
1173 }
1174 }
1175 }
1176 }
1177
1178
1179 void IsoTrig::endJob() {
1180 edm::LogWarning("IsoTrack") << trigNames_.size() << "Triggers were run. RunNo vs HLT accepts for";
1181 for (unsigned int i = 0; i < trigNames_.size(); ++i)
1182 edm::LogWarning("IsoTrack") << "[" << i << "]: " << trigNames_[i];
1183 unsigned int n = maxRunNo_ - minRunNo_ + 1;
1184 g_Pre = fs_->make<TH1D>("h_PrevsRN", "PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1185 g_PreL1 = fs_->make<TH1D>("h_PreL1vsRN", "L1 PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1186 g_PreHLT = fs_->make<TH1D>("h_PreHLTvsRN", "HLT PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1187 g_Accepts = fs_->make<TH1I>("h_HLTAcceptsvsRN", "HLT Accepts Vs Run Number", n, minRunNo_, maxRunNo_);
1188
1189 for (auto const &[runNum, nAccept] : trigList_) {
1190 auto const &triggerPrescales = trigPreList_[runNum];
1191 auto const preL1 = triggerPrescales.first;
1192 auto const preHLT = triggerPrescales.second;
1193 edm::LogVerbatim("IsoTrack") << runNum << " " << nAccept << " " << preL1 << " " << preHLT;
1194 g_Accepts->Fill(runNum, nAccept);
1195 g_PreL1->Fill(runNum, preL1);
1196 g_PreHLT->Fill(runNum, preHLT);
1197 g_Pre->Fill(runNum, preL1 * preHLT);
1198 }
1199 }
1200
1201 void IsoTrig::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
1202 edm::LogWarning("IsoTrack") << "Run " << iRun.run() << " hltconfig.init "
1203 << hltPrescaleProvider_.init(iRun, iSetup, processName_, changed_);
1204 }
1205
1206 void IsoTrig::StudyTrkEbyP(edm::Handle<reco::TrackCollection> &trkCollection) {
1207 t_TrkselTkFlag->clear();
1208 t_TrkqltyFlag->clear();
1209 t_TrkMissFlag->clear();
1210 t_TrkPVFlag->clear();
1211 t_TrkNuIsolFlag->clear();
1212 t_TrkhCone->clear();
1213 t_TrkP->clear();
1214
1215 if (!trkCollection.isValid()) {
1216 edm::LogVerbatim("IsoTrack") << "trkCollection.isValid is false";
1217 } else {
1218 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1219 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1220 spr::propagateCALO(
1221 trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections1, ((verbosity_ / 100) % 10 > 2));
1222 unsigned int nTracks = 0;
1223 int nRH_eMipDR = 0, nNearTRKs = 0;
1224 std::vector<bool> selFlags;
1225 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++, nTracks++) {
1226 double conehmaxNearP = 0, hCone = 0, eMipDR = 0.0;
1227 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1228 if (verbosity_ % 10 > 0)
1229 edm::LogVerbatim("IsoTrack") << "track no. " << nTracks << " p(): " << pTrack->p();
1230 if (pTrack->p() > 20) {
1231 math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1232 eMipDR = spr::eCone_ecal(geo_,
1233 barrelRecHitsHandle_,
1234 endcapRecHitsHandle_,
1235 trkDetItr->pointHCAL,
1236 trkDetItr->pointECAL,
1237 a_mipR_,
1238 trkDetItr->directionECAL,
1239 nRH_eMipDR);
1240 bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1241 spr::trackSelectionParameters oneCutParameters = selectionParameters_;
1242 oneCutParameters.maxDxyPV = 10;
1243 oneCutParameters.maxDzPV = 100;
1244 oneCutParameters.maxInMiss = 2;
1245 oneCutParameters.maxOutMiss = 2;
1246 bool qltyFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1247 oneCutParameters = selectionParameters_;
1248 oneCutParameters.maxDxyPV = 10;
1249 oneCutParameters.maxDzPV = 100;
1250 bool qltyMissFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1251 oneCutParameters = selectionParameters_;
1252 oneCutParameters.maxInMiss = 2;
1253 oneCutParameters.maxOutMiss = 2;
1254 bool qltyPVFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1255 if ((verbosity_ / 1000) % 10 > 1)
1256 edm::LogVerbatim("IsoTrack") << "sel " << selectTk << "ntracks " << nTracks << " a_charIsoR " << a_charIsoR_
1257 << " nNearTRKs " << nNearTRKs;
1258 conehmaxNearP = spr::chargeIsolationCone(
1259 nTracks, trkCaloDirections1, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1260 if ((verbosity_ / 1000) % 10 > 1)
1261 edm::LogVerbatim("IsoTrack") << "coneh " << conehmaxNearP << "ok " << trkDetItr->okECAL << " "
1262 << trkDetItr->okHCAL;
1263 double e1 = spr::eCone_ecal(geo_,
1264 barrelRecHitsHandle_,
1265 endcapRecHitsHandle_,
1266 trkDetItr->pointHCAL,
1267 trkDetItr->pointECAL,
1268 a_neutR1_,
1269 trkDetItr->directionECAL,
1270 nRH_eMipDR);
1271 double e2 = spr::eCone_ecal(geo_,
1272 barrelRecHitsHandle_,
1273 endcapRecHitsHandle_,
1274 trkDetItr->pointHCAL,
1275 trkDetItr->pointECAL,
1276 a_neutR2_,
1277 trkDetItr->directionECAL,
1278 nRH_eMipDR);
1279 double e_inCone = e2 - e1;
1280 bool chgIsolFlag = (conehmaxNearP < cutCharge_);
1281 bool mipFlag = (eMipDR < cutMip_);
1282 bool neuIsolFlag = (e_inCone < cutNeutral_);
1283 bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1284 selFlags.clear();
1285 selFlags.push_back(selectTk);
1286 selFlags.push_back(qltyFlag);
1287 selFlags.push_back(qltyMissFlag);
1288 selFlags.push_back(qltyPVFlag);
1289 if (verbosity_ % 10 > 0)
1290 edm::LogVerbatim("IsoTrack") << "emip: " << eMipDR << "<" << cutMip_ << "(" << mipFlag << ")"
1291 << " ; ok: " << trkDetItr->okECAL << "/" << trkDetItr->okHCAL
1292 << " ; chgiso: " << conehmaxNearP << "<" << cutCharge_ << "(" << chgIsolFlag
1293 << ")";
1294
1295 if (chgIsolFlag && mipFlag && trkpropFlag) {
1296 double distFromHotCell = -99.0;
1297 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1298 GlobalPoint gposHotCell(0., 0., 0.);
1299 std::vector<DetId> coneRecHitDetIds;
1300 hCone = spr::eCone_hcal(geo_,
1301 hbhe_,
1302 trkDetItr->pointHCAL,
1303 trkDetItr->pointECAL,
1304 a_coneR_,
1305 trkDetItr->directionHCAL,
1306 nRecHitsCone,
1307 coneRecHitDetIds,
1308 distFromHotCell,
1309 ietaHotCell,
1310 iphiHotCell,
1311 gposHotCell,
1312 -1);
1313
1314 t_TrkselTkFlag->push_back(selFlags[0]);
1315 t_TrkqltyFlag->push_back(selFlags[1]);
1316 t_TrkMissFlag->push_back(selFlags[2]);
1317 t_TrkPVFlag->push_back(selFlags[3]);
1318 t_TrkNuIsolFlag->push_back(neuIsolFlag);
1319 t_TrkhCone->push_back(hCone);
1320 t_TrkP->push_back(pTrack->p());
1321 }
1322 }
1323 }
1324 if (verbosity_ % 10 > 0)
1325 edm::LogVerbatim("IsoTrack") << "Filling " << t_TrkP->size() << " tracks in TrkRestree out of " << nTracks;
1326 }
1327 TrkResTree_->Fill();
1328 }
1329
1330 void IsoTrig::studyTiming(const edm::Event &theEvent) {
1331 t_timeL2Prod->clear();
1332 t_nPixCand->clear();
1333 t_nPixSeed->clear();
1334
1335 if (verbosity_ % 10 > 0) {
1336 edm::Handle<SeedingLayerSetsHits> hblayers, helayers;
1337 theEvent.getByToken(tok_SeedingLayerHB_, hblayers);
1338 theEvent.getByToken(tok_SeedingLayerHE_, helayers);
1339 const SeedingLayerSetsHits *layershb = hblayers.product();
1340 const SeedingLayerSetsHits *layershe = helayers.product();
1341 edm::LogVerbatim("IsoTrack") << "size of Seeding TripletLayers hb/he " << layershb->size() << "/"
1342 << layershe->size();
1343 edm::Handle<SiPixelRecHitCollection> rchts;
1344 theEvent.getByToken(tok_SiPixelRecHits_, rchts);
1345 const SiPixelRecHitCollection *rechits = rchts.product();
1346 edm::LogVerbatim("IsoTrack") << "size of SiPixelRechits " << rechits->size();
1347 }
1348 int nCandHB = pixelTrackRefsHB_.size(), nCandHE = pixelTrackRefsHE_.size();
1349 int nSeedHB = 0, nSeedHE = 0;
1350
1351 if (nCandHE > 0) {
1352 edm::Handle<reco::VertexCollection> pVertHB, pVertHE;
1353 theEvent.getByToken(tok_verthb_, pVertHB);
1354 theEvent.getByToken(tok_verthe_, pVertHE);
1355 edm::Handle<trigger::TriggerFilterObjectWithRefs> l1trigobj;
1356 theEvent.getByToken(tok_l1cand_, l1trigobj);
1357
1358 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1359 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1360 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
1361
1362 l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
1363 l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
1364 l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
1365
1366 double ptTriggered = -10;
1367 double etaTriggered = -100;
1368 double phiTriggered = -100;
1369 for (unsigned int p = 0; p < l1tauobjref.size(); p++) {
1370 if (l1tauobjref[p]->pt() > ptTriggered) {
1371 ptTriggered = l1tauobjref[p]->pt();
1372 phiTriggered = l1tauobjref[p]->phi();
1373 etaTriggered = l1tauobjref[p]->eta();
1374 }
1375 }
1376 for (unsigned int p = 0; p < l1jetobjref.size(); p++) {
1377 if (l1jetobjref[p]->pt() > ptTriggered) {
1378 ptTriggered = l1jetobjref[p]->pt();
1379 phiTriggered = l1jetobjref[p]->phi();
1380 etaTriggered = l1jetobjref[p]->eta();
1381 }
1382 }
1383 for (unsigned int p = 0; p < l1forjetobjref.size(); p++) {
1384 if (l1forjetobjref[p]->pt() > ptTriggered) {
1385 ptTriggered = l1forjetobjref[p]->pt();
1386 phiTriggered = l1forjetobjref[p]->phi();
1387 etaTriggered = l1forjetobjref[p]->eta();
1388 }
1389 }
1390 for (unsigned iS = 0; iS < pixelTrackRefsHE_.size(); iS++) {
1391 reco::VertexCollection::const_iterator vitSel;
1392 double minDZ = 100;
1393 bool vtxMatch;
1394 for (reco::VertexCollection::const_iterator vit = pVertHE->begin(); vit != pVertHE->end(); vit++) {
1395 if (fabs(pixelTrackRefsHE_[iS]->dz(vit->position())) < minDZ) {
1396 minDZ = fabs(pixelTrackRefsHE_[iS]->dz(vit->position()));
1397 vitSel = vit;
1398 }
1399 }
1400
1401 if ((fabs(pixelTrackRefsHE_[iS]->dxy(vitSel->position())) < vtxCutSeed_) || (minDZ == 100))
1402 vtxMatch = true;
1403
1404
1405 double R = deltaR(etaTriggered, phiTriggered, pixelTrackRefsHE_[iS]->eta(), pixelTrackRefsHE_[iS]->phi());
1406 if (R > tauUnbiasCone_ && vtxMatch)
1407 nSeedHE++;
1408 }
1409 for (unsigned iS = 0; iS < pixelTrackRefsHB_.size(); iS++) {
1410 reco::VertexCollection::const_iterator vitSel;
1411 double minDZ = 100;
1412 bool vtxMatch(false);
1413 for (reco::VertexCollection::const_iterator vit = pVertHB->begin(); vit != pVertHB->end(); vit++) {
1414 if (fabs(pixelTrackRefsHB_[iS]->dz(vit->position())) < minDZ) {
1415 minDZ = fabs(pixelTrackRefsHB_[iS]->dz(vit->position()));
1416 vitSel = vit;
1417 }
1418 }
1419
1420 if ((fabs(pixelTrackRefsHB_[iS]->dxy(vitSel->position())) < 101.0) || (minDZ == 100))
1421 vtxMatch = true;
1422
1423
1424 double R = deltaR(etaTriggered, phiTriggered, pixelTrackRefsHB_[iS]->eta(), pixelTrackRefsHB_[iS]->phi());
1425 if (R > 1.2 && vtxMatch)
1426 nSeedHB++;
1427 }
1428
1429 if (verbosity_ % 10 > 0) {
1430 edm::LogVerbatim("IsoTrack") << "(HB/HE) nCand: " << nCandHB << "/" << nCandHE << "nSeed: " << nSeedHB << "/"
1431 << nSeedHE;
1432 }
1433 }
1434 t_nPixSeed->push_back(nSeedHB);
1435 t_nPixSeed->push_back(nSeedHE);
1436 t_nPixCand->push_back(nCandHB);
1437 t_nPixCand->push_back(nCandHE);
1438
1439 TimingTree_->Fill();
1440 }
1441 void IsoTrig::studyMipCut(edm::Handle<reco::TrackCollection> &trkCollection,
1442 edm::Handle<reco::IsolatedPixelTrackCandidateCollection> &L2cands) {
1443 clearMipCutTreeVectors();
1444 if (verbosity_ % 10 > 0)
1445 edm::LogVerbatim("IsoTrack") << "inside studymipcut";
1446 if (!trkCollection.isValid()) {
1447 edm::LogWarning("IsoTrack") << "trkCollection.isValid is false";
1448 } else {
1449 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1450 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1451 spr::propagateCALO(
1452 trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections1, ((verbosity_ / 100) % 10 > 2));
1453 if (verbosity_ % 10 > 0)
1454 edm::LogVerbatim("IsoTrack") << "Number of L2cands:" << L2cands->size() << " to be matched to something out of "
1455 << trkCaloDirections1.size() << " reco tracks";
1456 for (unsigned int i = 0; i < L2cands->size(); i++) {
1457 edm::Ref<reco::IsolatedPixelTrackCandidateCollection> candref =
1458 edm::Ref<reco::IsolatedPixelTrackCandidateCollection>(L2cands, i);
1459 double enIn = candref->energyIn();
1460 h_EnIn->Fill(candref->energyIn());
1461 h_EnOut->Fill(candref->energyOut());
1462 math::XYZTLorentzVector v1(
1463 candref->track()->px(), candref->track()->py(), candref->track()->pz(), candref->track()->p());
1464 if (verbosity_ % 10 > 1)
1465 edm::LogVerbatim("IsoTrack") << "HLT Level Candidate eta/phi/pt/enIn:" << candref->track()->eta() << "/"
1466 << candref->track()->phi() << "/" << candref->track()->pt() << "/"
1467 << candref->energyIn();
1468 math::XYZTLorentzVector mindRvec;
1469 double mindR = 999.9, mindP1 = 999.9, eMipDR = 0.0;
1470 std::vector<bool> selFlags;
1471 unsigned int nTracks = 0;
1472 double conehmaxNearP = 0, hCone = 0;
1473 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++, nTracks++) {
1474 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1475 math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1476 double dr = dR(v1, v2);
1477 double dp1 = std::abs(1. / v1.r() - 1. / v2.r());
1478 if (verbosity_ % 1000 > 0)
1479 edm::LogVerbatim("IsoTrack") << "This recotrack(eta/phi/pt) " << pTrack->eta() << "/" << pTrack->phi() << "/"
1480 << pTrack->pt() << " has dr/dp= " << dr << "/" << dp1;
1481 if (dr < mindR) {
1482 mindR = dr;
1483 mindP1 = dp1;
1484 mindRvec = v2;
1485 int nRH_eMipDR = 0, nNearTRKs = 0;
1486 eMipDR = spr::eCone_ecal(geo_,
1487 barrelRecHitsHandle_,
1488 endcapRecHitsHandle_,
1489 trkDetItr->pointHCAL,
1490 trkDetItr->pointECAL,
1491 a_mipR_,
1492 trkDetItr->directionECAL,
1493 nRH_eMipDR);
1494 bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1495 spr::trackSelectionParameters oneCutParameters = selectionParameters_;
1496 oneCutParameters.maxDxyPV = 10;
1497 oneCutParameters.maxDzPV = 100;
1498 oneCutParameters.maxInMiss = 2;
1499 oneCutParameters.maxOutMiss = 2;
1500 bool qltyFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1501 oneCutParameters = selectionParameters_;
1502 oneCutParameters.maxDxyPV = 10;
1503 oneCutParameters.maxDzPV = 100;
1504 bool qltyMissFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1505 oneCutParameters = selectionParameters_;
1506 oneCutParameters.maxInMiss = 2;
1507 oneCutParameters.maxOutMiss = 2;
1508 bool qltyPVFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1509 conehmaxNearP = spr::chargeIsolationCone(
1510 nTracks, trkCaloDirections1, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1511 double e1 = spr::eCone_ecal(geo_,
1512 barrelRecHitsHandle_,
1513 endcapRecHitsHandle_,
1514 trkDetItr->pointHCAL,
1515 trkDetItr->pointECAL,
1516 a_neutR1_,
1517 trkDetItr->directionECAL,
1518 nRH_eMipDR);
1519 double e2 = spr::eCone_ecal(geo_,
1520 barrelRecHitsHandle_,
1521 endcapRecHitsHandle_,
1522 trkDetItr->pointHCAL,
1523 trkDetItr->pointECAL,
1524 a_neutR2_,
1525 trkDetItr->directionECAL,
1526 nRH_eMipDR);
1527 double e_inCone = e2 - e1;
1528 bool chgIsolFlag = (conehmaxNearP < cutCharge_);
1529 bool mipFlag = (eMipDR < cutMip_);
1530 bool neuIsolFlag = (e_inCone < cutNeutral_);
1531 bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1532 selFlags.clear();
1533 selFlags.push_back(selectTk);
1534 selFlags.push_back(qltyFlag);
1535 selFlags.push_back(qltyMissFlag);
1536 selFlags.push_back(qltyPVFlag);
1537 selFlags.push_back(trkpropFlag);
1538 selFlags.push_back(chgIsolFlag);
1539 selFlags.push_back(neuIsolFlag);
1540 selFlags.push_back(mipFlag);
1541 double distFromHotCell = -99.0;
1542 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1543 GlobalPoint gposHotCell(0., 0., 0.);
1544 std::vector<DetId> coneRecHitDetIds;
1545 hCone = spr::eCone_hcal(geo_,
1546 hbhe_,
1547 trkDetItr->pointHCAL,
1548 trkDetItr->pointECAL,
1549 a_coneR_,
1550 trkDetItr->directionHCAL,
1551 nRecHitsCone,
1552 coneRecHitDetIds,
1553 distFromHotCell,
1554 ietaHotCell,
1555 iphiHotCell,
1556 gposHotCell,
1557 -1);
1558 }
1559 }
1560 pushMipCutTreeVecs(v1, mindRvec, enIn, eMipDR, mindR, mindP1, selFlags, hCone);
1561 fillDifferences(6, v1, mindRvec, (verbosity_ % 10 > 0));
1562 if (mindR >= 0.05) {
1563 fillDifferences(8, v1, mindRvec, (verbosity_ % 10 > 0));
1564 h_MipEnNoMatch->Fill(candref->energyIn(), eMipDR);
1565 } else {
1566 fillDifferences(7, v1, mindRvec, (verbosity_ % 10 > 0));
1567 h_MipEnMatch->Fill(candref->energyIn(), eMipDR);
1568 }
1569 }
1570 }
1571 MipCutTree_->Fill();
1572 }
1573
1574 void IsoTrig::studyTrigger(edm::Handle<reco::TrackCollection> &trkCollection,
1575 std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1576 if (verbosity_ % 10 > 0)
1577 edm::LogVerbatim("IsoTrack") << "Inside StudyTrigger";
1578
1579 for (int j = 0; j < 3; j++) {
1580 for (unsigned int k = 0; k < vec_[j].size(); k++) {
1581 if (verbosity_ % 10 > 0)
1582 edm::LogVerbatim("IsoTrack") << "vec[" << j << "][" << k << "] pt " << vec_[j][k].pt() << " eta "
1583 << vec_[j][k].eta() << " phi " << vec_[j][k].phi();
1584 fillHist(j, vec_[j][k]);
1585 }
1586 }
1587
1588 double deta, dphi, dr;
1589
1590 for (int lvl = 1; lvl < 3; lvl++) {
1591 for (unsigned int i = 0; i < vec_[lvl].size(); i++) {
1592 deta = dEta(vec_[0][0], vec_[lvl][i]);
1593 dphi = dPhi(vec_[0][0], vec_[lvl][i]);
1594 dr = dR(vec_[0][0], vec_[lvl][i]);
1595 if (verbosity_ % 10 > 1)
1596 edm::LogVerbatim("IsoTrack") << "lvl " << lvl << " i " << i << " deta " << deta << " dphi " << dphi << " dR "
1597 << dr;
1598 h_dEtaL1[lvl - 1]->Fill(deta);
1599 h_dPhiL1[lvl - 1]->Fill(dphi);
1600 h_dRL1[lvl - 1]->Fill(std::sqrt(dr));
1601 }
1602 }
1603
1604 math::XYZTLorentzVector mindRvec;
1605 double mindR;
1606 for (unsigned int k = 0; k < vec_[2].size(); ++k) {
1607
1608 mindR = 999.9;
1609 if (verbosity_ % 10 > 1)
1610 edm::LogVerbatim("IsoTrack") << "L3obj: pt " << vec_[2][k].pt() << " eta " << vec_[2][k].eta() << " phi "
1611 << vec_[2][k].phi();
1612 for (unsigned int j = 0; j < vec_[1].size(); j++) {
1613 dr = dR(vec_[2][k], vec_[1][j]);
1614 if (dr < mindR) {
1615 mindR = dr;
1616 mindRvec = vec_[1][j];
1617 }
1618 }
1619 fillDifferences(0, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1620 if (mindR < 0.03) {
1621 fillDifferences(1, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1622 fillHist(6, mindRvec);
1623 fillHist(8, vec_[2][k]);
1624 } else {
1625 fillDifferences(2, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1626 fillHist(7, mindRvec);
1627 fillHist(9, vec_[2][k]);
1628 }
1629
1630
1631 mindR = 999.9;
1632 if (verbosity_ % 10 > 0)
1633 edm::LogVerbatim("IsoTrack") << "Now Matching L3 track with reco: L3 Track (eta, phi) " << vec_[2][k].eta() << ":"
1634 << vec_[2][k].phi() << " L2 Track " << mindRvec.eta() << ":" << mindRvec.phi()
1635 << " dR " << mindR;
1636 reco::TrackCollection::const_iterator goodTk = trkCollection->end();
1637 if (trkCollection.isValid()) {
1638 double mindP(9999.9);
1639 reco::TrackCollection::const_iterator trkItr;
1640 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++) {
1641 math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(), trkItr->pz(), trkItr->p());
1642 double deltaR = dR(v4, vec_[2][k]);
1643 double dp = std::abs(v4.r() / vec_[2][k].r() - 1.0);
1644 if (deltaR < mindR) {
1645 mindR = deltaR;
1646 mindP = dp;
1647 mindRvec = v4;
1648 goodTk = trkItr;
1649 }
1650 if ((verbosity_ / 10) % 10 > 1 && deltaR < 1.0)
1651 edm::LogVerbatim("IsoTrack") << "reco track: pt " << v4.pt() << " eta " << v4.eta() << " phi " << v4.phi()
1652 << " DR " << deltaR;
1653 }
1654 if (verbosity_ % 10 > 0)
1655 edm::LogVerbatim("IsoTrack") << "Now Matching at Reco level in step 1 DR: " << mindR << ":" << mindP
1656 << " eta:phi " << mindRvec.eta() << ":" << mindRvec.phi();
1657 if (mindR < 0.03 && mindP > 0.1) {
1658 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++) {
1659 math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(), trkItr->pz(), trkItr->p());
1660 double deltaR = dR(v4, vec_[2][k]);
1661 double dp = std::abs(v4.r() / vec_[2][k].r() - 1.0);
1662 if (dp < mindP && deltaR < 0.03) {
1663 mindR = deltaR;
1664 mindP = dp;
1665 mindRvec = v4;
1666 goodTk = trkItr;
1667 }
1668 }
1669 if (verbosity_ % 10 > 0)
1670 edm::LogVerbatim("IsoTrack") << "Now Matching at Reco level in step 2 DR: " << mindR << ":" << mindP
1671 << " eta:phi " << mindRvec.eta() << ":" << mindRvec.phi();
1672 }
1673 fillDifferences(3, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1674 fillHist(3, mindRvec);
1675 if (mindR < 0.03) {
1676 fillDifferences(4, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1677 fillHist(4, mindRvec);
1678 } else {
1679 fillDifferences(5, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1680 fillHist(5, mindRvec);
1681 }
1682 if (goodTk != trkCollection->end())
1683 goodTks.push_back(goodTk);
1684 }
1685 }
1686 }
1687
1688 void IsoTrig::studyIsolation(edm::Handle<reco::TrackCollection> &trkCollection,
1689 std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1690 if (trkCollection.isValid()) {
1691 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1692 spr::propagateCALO(
1693 trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections, ((verbosity_ / 100) % 10 > 2));
1694
1695 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1696 if ((verbosity_ / 1000) % 10 > 1) {
1697 edm::LogVerbatim("IsoTrack") << "n of barrelRecHitsHandle " << barrelRecHitsHandle_->size();
1698 for (EcalRecHitCollection::const_iterator hit = barrelRecHitsHandle_->begin(); hit != barrelRecHitsHandle_->end();
1699 ++hit) {
1700 edm::LogVerbatim("IsoTrack") << "hit : detid(ieta,iphi) " << (EBDetId)(hit->id()) << " time " << hit->time()
1701 << " energy " << hit->energy();
1702 }
1703 edm::LogVerbatim("IsoTrack") << "n of endcapRecHitsHandle " << endcapRecHitsHandle_->size();
1704 for (EcalRecHitCollection::const_iterator hit = endcapRecHitsHandle_->begin(); hit != endcapRecHitsHandle_->end();
1705 ++hit) {
1706 edm::LogVerbatim("IsoTrack") << "hit : detid(ieta,iphi) " << (EEDetId)(hit->id()) << " time " << hit->time()
1707 << " energy " << hit->energy();
1708 }
1709 edm::LogVerbatim("IsoTrack") << "n of hbhe " << hbhe_->size();
1710 for (HBHERecHitCollection::const_iterator hit = hbhe_->begin(); hit != hbhe_->end(); ++hit) {
1711 edm::LogVerbatim("IsoTrack") << "hit : detid(ieta,iphi) " << hit->id() << " time " << hit->time() << " energy "
1712 << hit->energy();
1713 }
1714 }
1715 unsigned int nTracks = 0;
1716 int ieta = 999;
1717 for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++, nTracks++) {
1718 bool l3Track = (std::find(goodTks.begin(), goodTks.end(), trkDetItr->trkItr) != goodTks.end());
1719 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1720 math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1721 bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1722 double eMipDR = 9999., e_inCone = 0, conehmaxNearP = 0, mindR = 999.9, hCone = 0;
1723 if (trkDetItr->okHCAL) {
1724 HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
1725 ieta = detId.ieta();
1726 }
1727 for (unsigned k = 0; k < vec_[0].size(); ++k) {
1728 double deltaR = dR(v4, vec_[0][k]);
1729 if (deltaR < mindR)
1730 mindR = deltaR;
1731 }
1732 if ((verbosity_ / 100) % 10 > 1)
1733 edm::LogVerbatim("IsoTrack") << "Track ECAL " << trkDetItr->okECAL << " HCAL " << trkDetItr->okHCAL << " Flag "
1734 << selectTk;
1735 if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL) {
1736 int nRH_eMipDR = 0, nNearTRKs = 0;
1737 double e1 = spr::eCone_ecal(geo_,
1738 barrelRecHitsHandle_,
1739 endcapRecHitsHandle_,
1740 trkDetItr->pointHCAL,
1741 trkDetItr->pointECAL,
1742 a_neutR1_,
1743 trkDetItr->directionECAL,
1744 nRH_eMipDR);
1745 double e2 = spr::eCone_ecal(geo_,
1746 barrelRecHitsHandle_,
1747 endcapRecHitsHandle_,
1748 trkDetItr->pointHCAL,
1749 trkDetItr->pointECAL,
1750 a_neutR2_,
1751 trkDetItr->directionECAL,
1752 nRH_eMipDR);
1753 eMipDR = spr::eCone_ecal(geo_,
1754 barrelRecHitsHandle_,
1755 endcapRecHitsHandle_,
1756 trkDetItr->pointHCAL,
1757 trkDetItr->pointECAL,
1758 a_mipR_,
1759 trkDetItr->directionECAL,
1760 nRH_eMipDR);
1761 conehmaxNearP =
1762 spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1763 e_inCone = e2 - e1;
1764 double distFromHotCell = -99.0;
1765 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1766 GlobalPoint gposHotCell(0., 0., 0.);
1767 std::vector<DetId> coneRecHitDetIds;
1768 hCone = spr::eCone_hcal(geo_,
1769 hbhe_,
1770 trkDetItr->pointHCAL,
1771 trkDetItr->pointECAL,
1772 a_coneR_,
1773 trkDetItr->directionHCAL,
1774 nRecHitsCone,
1775 coneRecHitDetIds,
1776 distFromHotCell,
1777 ietaHotCell,
1778 iphiHotCell,
1779 gposHotCell,
1780 -1);
1781 }
1782 if (l3Track) {
1783 fillHist(10, v4);
1784 if (selectTk) {
1785 fillHist(11, v4);
1786 fillCuts(0, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR > dr_L1_));
1787 if (conehmaxNearP < cutCharge_) {
1788 fillHist(12, v4);
1789 if (eMipDR < cutMip_) {
1790 fillHist(13, v4);
1791 fillEnergy(1, ieta, hCone, eMipDR, v4);
1792 fillEnergy(0, ieta, hCone, eMipDR, v4);
1793 if (e_inCone < cutNeutral_) {
1794 fillHist(14, v4);
1795 fillEnergy(4, ieta, hCone, eMipDR, v4);
1796 fillEnergy(3, ieta, hCone, eMipDR, v4);
1797 }
1798 }
1799 }
1800 }
1801 } else if (doL2L3_) {
1802 fillHist(15, v4);
1803 if (selectTk) {
1804 fillHist(16, v4);
1805 fillCuts(1, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR > dr_L1_));
1806 if (conehmaxNearP < cutCharge_) {
1807 fillHist(17, v4);
1808 if (eMipDR < cutMip_) {
1809 fillHist(18, v4);
1810 fillEnergy(2, ieta, hCone, eMipDR, v4);
1811 fillEnergy(0, ieta, hCone, eMipDR, v4);
1812 if (e_inCone < cutNeutral_) {
1813 fillHist(19, v4);
1814 fillEnergy(5, ieta, hCone, eMipDR, v4);
1815 fillEnergy(3, ieta, hCone, eMipDR, v4);
1816 }
1817 }
1818 }
1819 }
1820 }
1821 }
1822 }
1823 }
1824
1825 void IsoTrig::chgIsolation(double &etaTriggered,
1826 double &phiTriggered,
1827 edm::Handle<reco::TrackCollection> &trkCollection,
1828 const edm::Event &theEvent) {
1829 clearChgIsolnTreeVectors();
1830 if (verbosity_ % 10 > 0)
1831 edm::LogVerbatim("IsoTrack") << "Inside chgIsolation() with eta/phi Triggered: " << etaTriggered << "/"
1832 << phiTriggered;
1833 std::vector<double> maxP;
1834
1835 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1836 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1837 spr::propagateCALO(trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections1, ((verbosity_ / 100) % 10 > 2));
1838 if (verbosity_ % 10 > 0)
1839 edm::LogVerbatim("IsoTrack") << "Propagated TrkCollection";
1840 for (unsigned int k = 0; k < pixelIsolationConeSizeAtEC_.size(); ++k)
1841 maxP.push_back(0);
1842 unsigned i = pixelTrackRefsHE_.size();
1843 std::vector<std::pair<unsigned int, std::pair<double, double>>> VecSeedsatEC;
1844
1845 for (unsigned iS = 0; iS < pixelTrackRefsHE_.size(); iS++) {
1846 if (pixelTrackRefsHE_[iS]->p() > minPTrackValue_) {
1847 bool vtxMatch = false;
1848
1849 unsigned int ivSel = recVtxs_->size();
1850 double minDZ = 100;
1851 for (unsigned int iv = 0; iv < recVtxs_->size(); ++iv) {
1852 if (fabs(pixelTrackRefsHE_[iS]->dz((*recVtxs_)[iv].position())) < minDZ) {
1853 minDZ = fabs(pixelTrackRefsHE_[iS]->dz((*recVtxs_)[iv].position()));
1854 ivSel = iv;
1855 }
1856 }
1857
1858 if (ivSel == recVtxs_->size()) {
1859 vtxMatch = true;
1860 } else if (fabs(pixelTrackRefsHE_[iS]->dxy((*recVtxs_)[ivSel].position())) < vtxCutSeed_) {
1861 vtxMatch = true;
1862 }
1863
1864 double R = deltaR(etaTriggered, phiTriggered, pixelTrackRefsHE_[iS]->eta(), pixelTrackRefsHE_[iS]->phi());
1865 if (R > tauUnbiasCone_ && vtxMatch) {
1866
1867 std::pair<double, double> seedCooAtEC;
1868
1869 if (minDZ != 100)
1870 seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefsHE_[iS]->eta(),
1871 pixelTrackRefsHE_[iS]->phi(),
1872 pixelTrackRefsHE_[iS]->pt(),
1873 pixelTrackRefsHE_[iS]->charge(),
1874 (*recVtxs_)[ivSel].z());
1875
1876 else
1877 seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefsHE_[iS]->eta(),
1878 pixelTrackRefsHE_[iS]->phi(),
1879 pixelTrackRefsHE_[iS]->pt(),
1880 pixelTrackRefsHE_[iS]->charge(),
1881 0);
1882 VecSeedsatEC.push_back(std::make_pair(iS, seedCooAtEC));
1883 }
1884 }
1885 }
1886 for (unsigned int l = 0; l < VecSeedsatEC.size(); l++) {
1887 unsigned int iSeed = VecSeedsatEC[l].first;
1888 math::XYZTLorentzVector v1(pixelTrackRefsHE_[iSeed]->px(),
1889 pixelTrackRefsHE_[iSeed]->py(),
1890 pixelTrackRefsHE_[iSeed]->pz(),
1891 pixelTrackRefsHE_[iSeed]->p());
1892
1893 for (unsigned int j = 0; j < VecSeedsatEC.size(); j++) {
1894 unsigned int iSurr = VecSeedsatEC[j].first;
1895 if (iSeed != iSurr) {
1896
1897
1898
1899 if (deltaR(pixelTrackRefsHE_[iSeed]->eta(),
1900 pixelTrackRefsHE_[iSeed]->phi(),
1901 pixelTrackRefsHE_[iSurr]->eta(),
1902 pixelTrackRefsHE_[iSurr]->phi()) < prelimCone_) {
1903 unsigned int ivSel = recVtxs_->size();
1904 double minDZ2 = 100;
1905 for (unsigned int iv = 0; iv < recVtxs_->size(); ++iv) {
1906 if (fabs(pixelTrackRefsHE_[iSurr]->dz((*recVtxs_)[iv].position())) < minDZ2) {
1907 minDZ2 = fabs(pixelTrackRefsHE_[iSurr]->dz((*recVtxs_)[iv].position()));
1908 ivSel = iv;
1909 }
1910 }
1911
1912 if (minDZ2 == 100 || fabs(pixelTrackRefsHE_[iSurr]->dxy((*recVtxs_)[ivSel].position())) < vtxCutIsol_) {
1913
1914 double dist = getDistInCM(VecSeedsatEC[i].second.first,
1915 VecSeedsatEC[i].second.second,
1916 VecSeedsatEC[j].second.first,
1917 VecSeedsatEC[j].second.second);
1918 for (unsigned int k = 0; k < pixelIsolationConeSizeAtEC_.size(); ++k) {
1919 if (dist < pixelIsolationConeSizeAtEC_[k]) {
1920 if (pixelTrackRefsHE_[iSurr]->p() > maxP[k])
1921 maxP[k] = pixelTrackRefsHE_[iSurr]->p();
1922 }
1923 }
1924 }
1925 }
1926 }
1927 }
1928
1929 double conehmaxNearP = -1;
1930 bool selectTk = false;
1931 double mindR = 999.9;
1932 int nTracks = 0;
1933 math::XYZTLorentzVector mindRvec;
1934 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++, nTracks++) {
1935 int nNearTRKs = 0;
1936 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1937 math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1938 double dr = dR(v1, v2);
1939 if (dr < mindR) {
1940 selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1941 conehmaxNearP = spr::chargeIsolationCone(
1942 nTracks, trkCaloDirections1, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1943 mindR = dr;
1944 mindRvec = v2;
1945 }
1946 }
1947 pushChgIsolnTreeVecs(v1, mindRvec, maxP, conehmaxNearP, selectTk);
1948 }
1949 ChgIsolnTree_->Fill();
1950 }
1951
1952 void IsoTrig::getGoodTracks(const edm::Event &iEvent, edm::Handle<reco::TrackCollection> &trkCollection) {
1953 t_nGoodTk->clear();
1954 std::vector<int> nGood(4, 0);
1955 if (trkCollection.isValid()) {
1956 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1957 spr::propagateCALO(
1958 trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections, ((verbosity_ / 100) % 10 > 2));
1959
1960
1961 edm::Handle<trigger::TriggerFilterObjectWithRefs> l1trigobj;
1962 iEvent.getByToken(tok_l1cand_, l1trigobj);
1963
1964 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1965 l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
1966 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1967 l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
1968 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
1969 l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
1970
1971 double ptTriggered(-10), etaTriggered(-100), phiTriggered(-100);
1972 for (unsigned int p = 0; p < l1tauobjref.size(); p++) {
1973 if (l1tauobjref[p]->pt() > ptTriggered) {
1974 ptTriggered = l1tauobjref[p]->pt();
1975 phiTriggered = l1tauobjref[p]->phi();
1976 etaTriggered = l1tauobjref[p]->eta();
1977 }
1978 }
1979 for (unsigned int p = 0; p < l1jetobjref.size(); p++) {
1980 if (l1jetobjref[p]->pt() > ptTriggered) {
1981 ptTriggered = l1jetobjref[p]->pt();
1982 phiTriggered = l1jetobjref[p]->phi();
1983 etaTriggered = l1jetobjref[p]->eta();
1984 }
1985 }
1986 for (unsigned int p = 0; p < l1forjetobjref.size(); p++) {
1987 if (l1forjetobjref[p]->pt() > ptTriggered) {
1988 ptTriggered = l1forjetobjref[p]->pt();
1989 phiTriggered = l1forjetobjref[p]->phi();
1990 etaTriggered = l1forjetobjref[p]->eta();
1991 }
1992 }
1993 double pTriggered = ptTriggered * cosh(etaTriggered);
1994 math::XYZTLorentzVector pTrigger(
1995 ptTriggered * cos(phiTriggered), ptTriggered * sin(phiTriggered), pTriggered * tanh(etaTriggered), pTriggered);
1996
1997 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1998 unsigned int nTracks(0);
1999 for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++, nTracks++) {
2000 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
2001 math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
2002 bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
2003 double mindR = dR(v4, pTrigger);
2004 if ((verbosity_ / 100) % 10 > 1)
2005 edm::LogVerbatim("IsoTrack") << "Track ECAL " << trkDetItr->okECAL << " HCAL " << trkDetItr->okHCAL << " Flag "
2006 << selectTk;
2007 if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL && mindR > 1.0) {
2008 int nRH_eMipDR(0), nNearTRKs(0);
2009 double eMipDR = spr::eCone_ecal(geo_,
2010 barrelRecHitsHandle_,
2011 endcapRecHitsHandle_,
2012 trkDetItr->pointHCAL,
2013 trkDetItr->pointECAL,
2014 a_mipR_,
2015 trkDetItr->directionECAL,
2016 nRH_eMipDR);
2017 double conehmaxNearP =
2018 spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
2019 if (conehmaxNearP < 2.0 && eMipDR < 1.0) {
2020 if (pTrack->p() >= 20 && pTrack->p() < 30) {
2021 ++nGood[0];
2022 } else if (pTrack->p() >= 30 && pTrack->p() < 40) {
2023 ++nGood[1];
2024 } else if (pTrack->p() >= 40 && pTrack->p() < 60) {
2025 ++nGood[2];
2026 } else if (pTrack->p() >= 60 && pTrack->p() < 100) {
2027 ++nGood[3];
2028 }
2029 }
2030 }
2031 }
2032 }
2033
2034 for (unsigned int ii = 0; ii < nGood.size(); ++ii)
2035 t_nGoodTk->push_back(nGood[ii]);
2036 }
2037
2038 void IsoTrig::fillHist(int indx, math::XYZTLorentzVector &vec) {
2039 h_p[indx]->Fill(vec.r());
2040 h_pt[indx]->Fill(vec.pt());
2041 h_eta[indx]->Fill(vec.eta());
2042 h_phi[indx]->Fill(vec.phi());
2043 }
2044
2045 void IsoTrig::fillDifferences(int indx, math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2, bool debug) {
2046 double dr = dR(vec1, vec2);
2047 double deta = dEta(vec1, vec2);
2048 double dphi = dPhi(vec1, vec2);
2049 double dpt = dPt(vec1, vec2);
2050 double dp = dP(vec1, vec2);
2051 double dinvpt = dinvPt(vec1, vec2);
2052 h_dEta[indx]->Fill(deta);
2053 h_dPhi[indx]->Fill(dphi);
2054 h_dPt[indx]->Fill(dpt);
2055 h_dP[indx]->Fill(dp);
2056 h_dinvPt[indx]->Fill(dinvpt);
2057 h_mindR[indx]->Fill(dr);
2058 if (debug)
2059 edm::LogVerbatim("IsoTrack") << "mindR for index " << indx << " is " << dr << " deta " << deta << " dphi " << dphi
2060 << " dpt " << dpt << " dinvpt " << dinvpt;
2061 }
2062
2063 void IsoTrig::fillCuts(
2064 int indx, double eMipDR, double conehmaxNearP, double e_inCone, math::XYZTLorentzVector &vec, int ieta, bool cut) {
2065 h_eMip[indx]->Fill(eMipDR);
2066 h_eMaxNearP[indx]->Fill(conehmaxNearP);
2067 h_eNeutIso[indx]->Fill(e_inCone);
2068 if ((conehmaxNearP < cutCharge_) && (eMipDR < cutMip_)) {
2069 for (int lim = 0; lim < 5; ++lim) {
2070 if ((vec.r() > pLimits_[lim]) && (vec.r() <= pLimits_[lim + 1])) {
2071 h_etaMipTracks[lim][indx][0]->Fill((double)(ieta));
2072 if (cut)
2073 h_etaMipTracks[lim][indx][1]->Fill((double)(ieta));
2074 if (e_inCone < cutNeutral_) {
2075 h_etaCalibTracks[lim][indx][0]->Fill((double)(ieta));
2076 if (cut)
2077 h_etaCalibTracks[lim][indx][1]->Fill((double)(ieta));
2078 }
2079 }
2080 }
2081 }
2082 }
2083
2084 void IsoTrig::fillEnergy(int indx, int ieta, double hCone, double eMipDR, math::XYZTLorentzVector &vec) {
2085 int kk(-1);
2086 if (ieta > 0 && ieta < 25)
2087 kk = 23 + ieta;
2088 else if (ieta > -25 && ieta < 0)
2089 kk = -(ieta + 1);
2090 if (kk >= 0 && eMipDR > 0.01 && hCone > 1.0) {
2091 for (int lim = 0; lim < 5; ++lim) {
2092 if ((vec.r() > pLimits_[lim]) && (vec.r() <= pLimits_[lim + 1])) {
2093 h_eHcal[lim][indx][kk]->Fill(hCone);
2094 h_eCalo[lim][indx][kk]->Fill(hCone + eMipDR);
2095 }
2096 }
2097 }
2098 }
2099
2100 double IsoTrig::dEta(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) { return (vec1.eta() - vec2.eta()); }
2101
2102 double IsoTrig::dPhi(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2103 double phi1 = vec1.phi();
2104 if (phi1 < 0)
2105 phi1 += 2.0 * M_PI;
2106 double phi2 = vec2.phi();
2107 if (phi2 < 0)
2108 phi2 += 2.0 * M_PI;
2109 double dphi = phi1 - phi2;
2110 if (dphi > M_PI)
2111 dphi -= 2. * M_PI;
2112 else if (dphi < -M_PI)
2113 dphi += 2. * M_PI;
2114 return dphi;
2115 }
2116
2117 double IsoTrig::dR(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2118 double deta = dEta(vec1, vec2);
2119 double dphi = dPhi(vec1, vec2);
2120 return std::sqrt(deta * deta + dphi * dphi);
2121 }
2122
2123 double IsoTrig::dPt(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) { return (vec1.pt() - vec2.pt()); }
2124
2125 double IsoTrig::dP(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2126 return (std::abs(vec1.r() - vec2.r()));
2127 }
2128
2129 double IsoTrig::dinvPt(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2130 return ((1 / vec1.pt()) - (1 / vec2.pt()));
2131 }
2132
2133 std::pair<double, double> IsoTrig::etaPhiTrigger() {
2134 double eta(0), phi(0), ptmax(0);
2135 for (unsigned int k = 0; k < vec_[0].size(); ++k) {
2136 if (k == 0) {
2137 eta = vec_[0][k].eta();
2138 phi = vec_[0][k].phi();
2139 ptmax = vec_[0][k].pt();
2140 } else if (vec_[0][k].pt() > ptmax) {
2141 eta = vec_[0][k].eta();
2142 phi = vec_[0][k].phi();
2143 ptmax = vec_[0][k].pt();
2144 }
2145 }
2146 return std::pair<double, double>(eta, phi);
2147 }
2148
2149 std::pair<double, double> IsoTrig::GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ) {
2150 double deltaPhi = 0;
2151 double etaEC = 100;
2152 double phiEC = 100;
2153
2154 double Rcurv = 9999999;
2155 if (bfVal_ != 0)
2156 Rcurv = pT * 33.3 * 100 / (bfVal_ * 10);
2157
2158 double ecDist = zEE_;
2159 double ecRad = rEB_;
2160 double theta = 2 * atan(exp(-etaIP));
2161 double zNew = 0;
2162 if (theta > 0.5 * M_PI)
2163 theta = M_PI - theta;
2164 if (fabs(etaIP) < 1.479) {
2165 if ((0.5 * ecRad / Rcurv) > 1) {
2166 etaEC = 10000;
2167 deltaPhi = 0;
2168 } else {
2169 deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
2170 double alpha1 = 2 * asin(0.5 * ecRad / Rcurv);
2171 double z = ecRad / tan(theta);
2172 if (etaIP > 0)
2173 zNew = z * (Rcurv * alpha1) / ecRad + vtxZ;
2174 else
2175 zNew = -z * (Rcurv * alpha1) / ecRad + vtxZ;
2176 double zAbs = fabs(zNew);
2177 if (zAbs < ecDist) {
2178 etaEC = -log(tan(0.5 * atan(ecRad / zAbs)));
2179 deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
2180 }
2181 if (zAbs > ecDist) {
2182 zAbs = (fabs(etaIP) / etaIP) * ecDist;
2183 double Zflight = fabs(zAbs - vtxZ);
2184 double alpha = (Zflight * ecRad) / (z * Rcurv);
2185 double Rec = 2 * Rcurv * sin(alpha / 2);
2186 deltaPhi = -charge * alpha / 2;
2187 etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
2188 }
2189 }
2190 } else {
2191 zNew = (fabs(etaIP) / etaIP) * ecDist;
2192 double Zflight = fabs(zNew - vtxZ);
2193 double Rvirt = fabs(Zflight * tan(theta));
2194 double Rec = 2 * Rcurv * sin(Rvirt / (2 * Rcurv));
2195 deltaPhi = -(charge) * (Rvirt / (2 * Rcurv));
2196 etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
2197 }
2198
2199 if (zNew < 0)
2200 etaEC = -etaEC;
2201 phiEC = phiIP + deltaPhi;
2202
2203 if (phiEC < -M_PI)
2204 phiEC += 2 * M_PI;
2205 if (phiEC > M_PI)
2206 phiEC -= 2 * M_PI;
2207
2208 std::pair<double, double> retVal(etaEC, phiEC);
2209 return retVal;
2210 }
2211
2212 double IsoTrig::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
2213 double Rec;
2214 double theta1 = 2 * atan(exp(-eta1));
2215 double theta2 = 2 * atan(exp(-eta2));
2216 if (fabs(eta1) < 1.479)
2217 Rec = rEB_;
2218 else if (fabs(eta1) > 1.479 && fabs(eta1) < 7.0)
2219 Rec = tan(theta1) * zEE_;
2220 else
2221 return 1000;
2222
2223
2224 double angle =
2225 acos((sin(theta1) * sin(theta2) * (sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2)) + cos(theta1) * cos(theta2)));
2226 if (angle < 0.5 * M_PI)
2227 return fabs((Rec / sin(theta1)) * tan(angle));
2228 else
2229 return 1000;
2230 }
2231
2232
2233 DEFINE_FWK_MODULE(IsoTrig);
