29 #include "TDirectory.h"
113 std::vector<double> &PixMaxP,
122 std::vector<bool> &Flags,
131 double &phiTriggered,
267 std::vector<math::XYZTLorentzVector>
vec_[3];
271 : hltPrescaleProvider_(iConfig, consumesCollector(), *this),
273 pixCandTag_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"pixCandTag")),
274 l1CandTag_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"l1CandTag")),
275 l2CandTag_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"l2CandTag")),
276 pixelTracksSources_(iConfig.getUntrackedParameter<
std::
vector<
edm::
InputTag>>(
"pixelTracksSources")),
277 doL2L3_(iConfig.getUntrackedParameter<
bool>(
"doL2L3",
true)),
278 doTiming_(iConfig.getUntrackedParameter<
bool>(
"doTimingTree",
true)),
279 doMipCutTree_(iConfig.getUntrackedParameter<
bool>(
"doMipCutTree",
true)),
280 doTrkResTree_(iConfig.getUntrackedParameter<
bool>(
"doTrkResTree",
true)),
281 doChgIsolTree_(iConfig.getUntrackedParameter<
bool>(
"doChgIsolTree",
true)),
282 doStudyIsol_(iConfig.getUntrackedParameter<
bool>(
"doStudyIsol",
true)),
283 verbosity_(iConfig.getUntrackedParameter<
int>(
"verbosity", 0)),
284 pixelIsolationConeSizeAtEC_(iConfig.getUntrackedParameter<
std::
vector<double>>(
"pixelIsolationConeSizeAtEC")),
285 minPTrackValue_(iConfig.getUntrackedParameter<double>(
"minPTrackValue")),
286 vtxCutSeed_(iConfig.getUntrackedParameter<double>(
"vertexCutSeed")),
287 vtxCutIsol_(iConfig.getUntrackedParameter<double>(
"vertexCutIsol")),
288 tauUnbiasCone_(iConfig.getUntrackedParameter<double>(
"tauUnbiasCone")),
289 prelimCone_(iConfig.getUntrackedParameter<double>(
"prelimCone")),
290 theTrackQuality_(iConfig.getUntrackedParameter<
std::
string>(
"trackQuality",
"highPurity")),
291 processName_(iConfig.getUntrackedParameter<
std::
string>(
"processName",
"HLT")),
292 dr_L1_(iConfig.getUntrackedParameter<double>(
"isolationL1", 1.0)),
293 a_coneR_(iConfig.getUntrackedParameter<double>(
"coneRadius", 34.98)),
294 a_charIsoR_(a_coneR_ + 28.9),
295 a_neutIsoR_(a_charIsoR_ * 0.726),
296 a_mipR_(iConfig.getUntrackedParameter<double>(
"coneRadiusMIP", 14.0)),
297 a_neutR1_(iConfig.getUntrackedParameter<double>(
"coneRadiusNeut1", 21.0)),
298 a_neutR2_(iConfig.getUntrackedParameter<double>(
"coneRadiusNeut2", 29.0)),
299 cutMip_(iConfig.getUntrackedParameter<double>(
"cutMIP", 1.0)),
300 cutCharge_(iConfig.getUntrackedParameter<double>(
"chargeIsolation", 2.0)),
301 cutNeutral_(iConfig.getUntrackedParameter<double>(
"neutralIsolation", 2.0)),
302 minRunNo_(iConfig.getUntrackedParameter<
int>(
"minRun")),
303 maxRunNo_(iConfig.getUntrackedParameter<
int>(
"maxRun")),
305 t_timeL2Prod(nullptr),
311 t_TrkselTkFlag(nullptr),
312 t_TrkqltyFlag(nullptr),
313 t_TrkMissFlag(nullptr),
314 t_TrkPVFlag(nullptr),
315 t_TrkNuIsolFlag(nullptr),
317 t_PixcandPt(nullptr),
318 t_PixcandEta(nullptr),
319 t_PixcandPhi(nullptr),
320 t_PixcandMaxP(nullptr),
321 t_PixTrkcandP(nullptr),
322 t_PixTrkcandPt(nullptr),
323 t_PixTrkcandEta(nullptr),
324 t_PixTrkcandPhi(nullptr),
325 t_PixTrkcandMaxP(nullptr),
326 t_PixTrkcandselTk(nullptr),
329 t_NFcandEta(nullptr),
330 t_NFcandPhi(nullptr),
331 t_NFcandEmip(nullptr),
332 t_NFTrkcandP(nullptr),
333 t_NFTrkcandPt(nullptr),
334 t_NFTrkcandEta(nullptr),
335 t_NFTrkcandPhi(nullptr),
336 t_NFTrkcandEmip(nullptr),
337 t_NFTrkMinDR(nullptr),
338 t_NFTrkMinDP1(nullptr),
339 t_NFTrkselTkFlag(nullptr),
340 t_NFTrkqltyFlag(nullptr),
341 t_NFTrkMissFlag(nullptr),
342 t_NFTrkPVFlag(nullptr),
343 t_NFTrkPropFlag(nullptr),
344 t_NFTrkChgIsoFlag(nullptr),
345 t_NFTrkNeuIsoFlag(nullptr),
346 t_NFTrkMipFlag(nullptr),
366 tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
368 tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel);
410 <<
"\t prelimCone " <<
prelimCone_ <<
"\t pixelIsolationConeSizeAtEC";
414 double pl[] = {20, 30, 40, 60, 80, 120};
415 for (
int i = 0;
i < 6; ++
i)
420 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
421 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
516 std::vector<std::string> triggers = {
"HLT_IsoTrackHB"};
518 std::vector<double> cones = {35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 63.9, 70.0};
520 desc.addUntracked<std::vector<std::string>>(
"Triggers", triggers);
524 desc.addUntracked<
bool>(
"doL2L3",
false);
525 desc.addUntracked<
bool>(
"doTimingTree",
false);
526 desc.addUntracked<
bool>(
"doMipCutTree",
false);
527 desc.addUntracked<
bool>(
"doTrkResTree",
true);
528 desc.addUntracked<
bool>(
"doChgIsolTree",
false);
529 desc.addUntracked<
bool>(
"doStudyIsol",
false);
530 desc.addUntracked<
int>(
"verbosity", 0);
533 desc.addUntracked<
double>(
"minTrackPt", 10.0);
534 desc.addUntracked<
double>(
"maxDxyPV", 0.02);
535 desc.addUntracked<
double>(
"maxDzPV", 0.02);
536 desc.addUntracked<
double>(
"maxChi2", 5.0);
537 desc.addUntracked<
double>(
"maxDpOverP", 0.1);
538 desc.addUntracked<
int>(
"minOuterHit", 4);
539 desc.addUntracked<
int>(
"minLayerCrossed", 8);
540 desc.addUntracked<
int>(
"maxInMiss", 0);
541 desc.addUntracked<
int>(
"maxOutMiss", 0);
542 desc.addUntracked<
double>(
"isolationL1", 1.0);
543 desc.addUntracked<
double>(
"coneRadius", 34.98);
544 desc.addUntracked<
double>(
"coneRadiusMIP", 14.0);
545 desc.addUntracked<
double>(
"coneRadiusNeut1", 21.0);
546 desc.addUntracked<
double>(
"coneRadiusNeut2", 29.0);
547 desc.addUntracked<
double>(
"cutMIP", 1.0);
548 desc.addUntracked<
double>(
"chargeIsolation", 2.0);
549 desc.addUntracked<
double>(
"neutralIsolation", 2.0);
550 desc.addUntracked<
int>(
"minRun", 190456);
551 desc.addUntracked<
int>(
"maxRun", 203002);
552 desc.addUntracked<std::vector<edm::InputTag>>(
"pixelTracksSources",
tags);
553 desc.addUntracked<std::vector<double>>(
"pixelIsolationConeSizeAtEC", cones);
554 desc.addUntracked<
double>(
"minPTrackValue", 0.0);
555 desc.addUntracked<
double>(
"vertexCutSeed", 101.0);
556 desc.addUntracked<
double>(
"vertexCutIsol", 101.0);
557 desc.addUntracked<
double>(
"tauUnbiasCone", 1.2);
558 desc.addUntracked<
double>(
"prelimCone", 1.0);
559 desc.add<
unsigned int>(
"stageL1Trigger", 1);
560 descriptions.
add(
"isoTrigDefault",
desc);
565 edm::LogVerbatim(
"IsoTrack") <<
"Event starts====================================";
567 int RunNo =
iEvent.id().run();
579 if (!triggerEventHandle.
isValid()) {
580 edm::LogWarning(
"IsoTrack") <<
"Error! Can't get the product hltTriggerSummaryAOD";
600 if (!
recVtxs_->empty() && !((*recVtxs_)[0].isFake())) {
617 for (reco::TrackCollection::const_iterator pit = iPixCol->begin(); pit != iPixCol->end(); pit++) {
627 for (
unsigned int ifilter = 0; ifilter <
triggerEvent.sizeFilters(); ++ifilter) {
628 std::string FilterNames[7] = {
"hltL1sL1SingleJet68",
629 "hltIsolPixelTrackL2FilterHE",
630 "ecalIsolPixelTrackFilterHE",
631 "hltIsolPixelTrackL3FilterHE",
632 "hltIsolPixelTrackL2FilterHB",
633 "ecalIsolPixelTrackFilterHB",
634 "hltIsolPixelTrackL3FilterHB"};
636 for (
int i = 0;
i < 7;
i++) {
637 if (
label == FilterNames[
i])
649 <<
iEvent.luminosityBlock() <<
" Bunch " <<
iEvent.bunchCrossing() <<
" mybxlumi "
652 edm::LogWarning(
"IsoTrack") <<
"Error! Can't get the product triggerResults";
656 std::vector<std::string>
modules;
660 const std::vector<std::string> &triggerNames_ =
triggerNames.triggerNames();
665 unsigned int triggerindx =
hltConfig.triggerIndex(triggerNames_[
i]);
666 const std::vector<std::string> &moduleLabels(
hltConfig.moduleLabels(triggerindx));
672 edm::LogVerbatim(
"IsoTrack") <<
"trigger that i want " << triggerNames_[
i] <<
" accept "
689 << triggerNames_[
i] <<
" accept " <<
hlt <<
" preL1 " << preL1 <<
" preHLT " << preHLT;
690 for (
int iv = 0;
iv < 3; ++
iv)
695 trigList_.insert(std::pair<unsigned int, unsigned int>(RunNo, 1));
699 for (
unsigned int ifilter = 0; ifilter <
triggerEvent.sizeFilters(); ++ifilter) {
700 std::vector<int>
Keys;
703 for (
unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
704 if (
label.find(moduleLabels[imodule]) != std::string::npos) {
707 for (
unsigned int ifiltrKey = 0; ifiltrKey <
triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
711 if (
label.find(
"L2Filter") != std::string::npos) {
712 vec_[1].push_back(v4);
713 }
else if (
label.find(
"L3Filter") != std::string::npos) {
714 vec_[2].push_back(v4);
716 vec_[0].push_back(v4);
721 <<
" phi " << TO.
phi() <<
" mass " << TO.
mass() <<
" Id " << TO.
id();
726 std::vector<reco::TrackCollection::const_iterator> goodTks;
732 reco::TrackCollection::const_iterator trkItr;
733 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++)
734 goodTks.push_back(trkItr);
774 for (
unsigned itrig = 0; itrig < triggerNames_.size(); itrig++) {
775 unsigned int triggerindx =
hltConfig.triggerIndex(triggerNames_[itrig]);
776 if (triggerindx >=
n)
777 edm::LogVerbatim(
"IsoTrack") << triggerNames_[itrig] <<
" " << triggerindx <<
" does not exist in"
778 <<
" the current menu";
780 edm::LogVerbatim(
"IsoTrack") << triggerNames_[itrig] <<
" " << triggerindx <<
" exists";
831 std::vector<double> &PixMaxP,
853 std::vector<bool> &Flags,
879 char hname[100], htit[100];
881 "L2Match",
"L2NoMatch",
"L3Match",
"L3NoMatch",
"HLTTrk",
"HLTGoodTrk",
882 "HLTIsoTrk",
"HLTMip",
"HLTSelect",
"nonHLTTrk",
"nonHLTGoodTrk",
"nonHLTIsoTrk",
883 "nonHLTMip",
"nonHLTSelect"};
899 t_TrkP =
new std::vector<double>();
962 t_ECone =
new std::vector<double>();
987 std::string FilterNames[7] = {
"hltL1sL1SingleJet68",
988 "hltIsolPixelTrackL2FilterHE",
989 "ecalIsolPixelTrackFilterHE",
990 "hltIsolPixelTrackL3FilterHE",
991 "hltIsolPixelTrackL2FilterHB",
992 "ecalIsolPixelTrackFilterHB",
993 "hltIsolPixelTrackL3FilterHB"};
994 for (
int i = 0;
i < 7;
i++) {
995 h_Filters->GetXaxis()->SetBinLabel(
i + 1, FilterNames[
i].c_str());
998 h_nHLT =
fs_->
make<TH1I>(
"h_nHLT",
"Size of trigger Names", 1000, 1, 1000);
999 h_HLT =
fs_->
make<TH1I>(
"h_HLT",
"HLT accept", 3, -1, 2);
1000 h_PreL1 =
fs_->
make<TH1I>(
"h_PreL1",
"L1 Prescale", 500, 0, 500);
1002 h_Pre =
fs_->
make<TH1I>(
"h_Pre",
"Prescale", 3000, 0, 3000);
1004 h_PreL1wt =
fs_->
make<TH1D>(
"h_PreL1wt",
"Weighted L1 Prescale", 500, 0, 500);
1005 h_PreHLTwt =
fs_->
make<TH1D>(
"h_PreHLTwt",
"Weighted HLT Prescale", 500, 0, 100);
1008 h_EnIn =
fs_->
make<TH1D>(
"h_EnInEcal",
"EnergyIn Ecal", 200, 0.0, 20.0);
1009 h_EnOut =
fs_->
make<TH1D>(
"h_EnOutEcal",
"EnergyOut Ecal", 200, 0.0, 20.0);
1011 fs_->
make<TH2D>(
"h_MipEnMatch",
"MipEn: HLT level vs Reco Level (Matched)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1013 "h_MipEnNoMatch",
"MipEn: HLT level vs Reco Level (No Match Found)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1016 h_nL3Objs =
fs_->
make<TH1I>(
"h_nL3Objs",
"Number of L3 objects", 10, 0, 10);
1025 "NewFilterRecoMatch",
1026 "NewFilterRecoNoMatch"};
1027 for (
int ipair = 0; ipair < 9; ipair++) {
1028 sprintf(hname,
"h_dEta%s", pairs[ipair].c_str());
1029 sprintf(htit,
"#Delta#eta for %s", pairs[ipair].c_str());
1030 h_dEta[ipair] =
fs_->
make<TH1D>(hname, htit, 200, -10.0, 10.0);
1031 h_dEta[ipair]->GetXaxis()->SetTitle(
"d#eta");
1033 sprintf(hname,
"h_dPhi%s", pairs[ipair].c_str());
1034 sprintf(htit,
"#Delta#phi for %s", pairs[ipair].c_str());
1035 h_dPhi[ipair] =
fs_->
make<TH1D>(hname, htit, 140, -7.0, 7.0);
1036 h_dPhi[ipair]->GetXaxis()->SetTitle(
"d#phi");
1038 sprintf(hname,
"h_dPt%s", pairs[ipair].c_str());
1039 sprintf(htit,
"#Delta dp_{T} for %s objects", pairs[ipair].c_str());
1040 h_dPt[ipair] =
fs_->
make<TH1D>(hname, htit, 400, -200.0, 200.0);
1041 h_dPt[ipair]->GetXaxis()->SetTitle(
"dp_{T} (GeV)");
1043 sprintf(hname,
"h_dP%s", pairs[ipair].c_str());
1044 sprintf(htit,
"#Delta p for %s objects", pairs[ipair].c_str());
1045 h_dP[ipair] =
fs_->
make<TH1D>(hname, htit, 400, -200.0, 200.0);
1046 h_dP[ipair]->GetXaxis()->SetTitle(
"dP (GeV)");
1048 sprintf(hname,
"h_dinvPt%s", pairs[ipair].c_str());
1049 sprintf(htit,
"#Delta (1/p_{T}) for %s objects", pairs[ipair].c_str());
1051 h_dinvPt[ipair]->GetXaxis()->SetTitle(
"d(1/p_{T})");
1052 sprintf(hname,
"h_mindR%s", pairs[ipair].c_str());
1053 sprintf(htit,
"min(#Delta R) for %s objects", pairs[ipair].c_str());
1055 h_mindR[ipair]->GetXaxis()->SetTitle(
"dR");
1058 for (
int lvl = 0; lvl < 2; lvl++) {
1059 sprintf(hname,
"h_dEtaL1%s",
levels[lvl + 1].c_str());
1060 sprintf(htit,
"#Delta#eta for L1 and %s objects",
levels[lvl + 1].c_str());
1063 sprintf(hname,
"h_dPhiL1%s",
levels[lvl + 1].c_str());
1064 sprintf(htit,
"#Delta#phi for L1 and %s objects",
levels[lvl + 1].c_str());
1067 sprintf(hname,
"h_dRL1%s",
levels[lvl + 1].c_str());
1068 sprintf(htit,
"#Delta R for L1 and %s objects",
levels[lvl + 1].c_str());
1073 int levmin = (
doL2L3_ ? 0 : 10);
1074 for (
int ilevel = levmin; ilevel < 20; ilevel++) {
1075 sprintf(hname,
"h_p%s",
levels[ilevel].c_str());
1076 sprintf(htit,
"p for %s objects",
levels[ilevel].c_str());
1077 h_p[ilevel] =
fs_->
make<TH1D>(hname, htit, 100, 0.0, 500.0);
1078 h_p[ilevel]->GetXaxis()->SetTitle(
"p (GeV)");
1080 sprintf(hname,
"h_pt%s",
levels[ilevel].c_str());
1081 sprintf(htit,
"p_{T} for %s objects",
levels[ilevel].c_str());
1082 h_pt[ilevel] =
fs_->
make<TH1D>(hname, htit, 100, 0.0, 500.0);
1083 h_pt[ilevel]->GetXaxis()->SetTitle(
"p_{T} (GeV)");
1085 sprintf(hname,
"h_eta%s",
levels[ilevel].c_str());
1086 sprintf(htit,
"#eta for %s objects",
levels[ilevel].c_str());
1087 h_eta[ilevel] =
fs_->
make<TH1D>(hname, htit, 100, -5.0, 5.0);
1088 h_eta[ilevel]->GetXaxis()->SetTitle(
"#eta");
1090 sprintf(hname,
"h_phi%s",
levels[ilevel].c_str());
1091 sprintf(htit,
"#phi for %s objects",
levels[ilevel].c_str());
1092 h_phi[ilevel] =
fs_->
make<TH1D>(hname, htit, 70, -3.5, 3.50);
1093 h_phi[ilevel]->GetXaxis()->SetTitle(
"#phi");
1098 for (
int icut = 0; icut < 2; icut++) {
1099 sprintf(hname,
"h_eMip%s",
cuts[icut].c_str());
1100 sprintf(htit,
"eMip for %s tracks",
cuts[icut].c_str());
1101 h_eMip[icut] =
fs_->
make<TH1D>(hname, htit, 200, 0.0, 10.0);
1102 h_eMip[icut]->GetXaxis()->SetTitle(
"E_{Mip} (GeV)");
1104 sprintf(hname,
"h_eMaxNearP%s",
cuts[icut].c_str());
1105 sprintf(htit,
"eMaxNearP for %s tracks",
cuts[icut].c_str());
1107 h_eMaxNearP[icut]->GetXaxis()->SetTitle(
"E_{MaxNearP} (GeV)");
1109 sprintf(hname,
"h_eNeutIso%s",
cuts[icut].c_str());
1110 sprintf(htit,
"eNeutIso for %s ",
cuts[icut].c_str());
1112 h_eNeutIso[icut]->GetXaxis()->SetTitle(
"E_{NeutIso} (GeV)");
1114 for (
int kcut = 0; kcut < 2; ++kcut) {
1115 for (
int lim = 0; lim < 5; ++lim) {
1116 sprintf(hname,
"h_etaCalibTracks%sCut%dLim%d",
cuts[icut].c_str(), kcut, lim);
1118 "#eta for %s isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1122 cuts2[kcut].c_str());
1126 sprintf(hname,
"h_etaMipTracks%sCut%dLim%d",
cuts[icut].c_str(), kcut, lim);
1128 "#eta for %s charge isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1132 cuts2[kcut].c_str());
1139 std::string ecut1[3] = {
"all",
"HLTMatched",
"HLTNotMatched"};
1141 int etac[48] = {-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, -13, -14, -15, -16,
1142 -17, -18, -19, -20, -21, -22, -23, -24, 1, 2, 3, 4, 5, 6, 7, 8,
1143 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24};
1144 for (
int icut = 0; icut < 6; icut++) {
1146 int i1 = (icut > 2 ? 1 : 0);
1147 int i2 = icut -
i1 * 3;
1148 for (
int kcut = 0; kcut < 48; kcut++) {
1149 for (
int lim = 0; lim < 5; ++lim) {
1150 sprintf(hname,
"h_eta%dEnHcal%s%s%d", etac[kcut], ecut1[
i2].c_str(), ecut2[
i1].c_str(), lim);
1152 "HCAL energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1158 h_eHcal[lim][icut][kcut] =
fs_->
make<TH1D>(hname, htit, 750, 0.0, 150.0);
1159 h_eHcal[lim][icut][kcut]->GetXaxis()->SetTitle(
"Energy (GeV)");
1160 sprintf(hname,
"h_eta%dEnCalo%s%s%d", etac[kcut], ecut1[
i2].c_str(), ecut2[
i1].c_str(), lim);
1162 "Calorimter energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1168 h_eCalo[lim][icut][kcut] =
fs_->
make<TH1D>(hname, htit, 750, 0.0, 150.0);
1169 h_eCalo[lim][icut][kcut]->GetXaxis()->SetTitle(
"Energy (GeV)");
1177 unsigned int preL1, preHLT;
1178 std::map<unsigned int, unsigned int>::iterator itr;
1179 std::map<unsigned int, const std::pair<int, int>>::iterator itrPre;
1190 preL1 = (itrPre->second).
first;
1191 preHLT = (itrPre->second).
second;
1192 edm::LogVerbatim(
"IsoTrack") << itr->first <<
" " << itr->second <<
" " << itrPre->first <<
" " << preL1 <<
" "
1194 g_Accepts->Fill(itr->first, itr->second);
1195 g_PreL1->Fill(itr->first, preL1);
1196 g_PreHLT->Fill(itr->first, preHLT);
1197 g_Pre->Fill(itr->first, preL1 * preHLT);
1215 if (!trkCollection.
isValid()) {
1218 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1219 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
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));
1230 if (pTrack->
p() > 20) {
1235 trkDetItr->pointHCAL,
1236 trkDetItr->pointECAL,
1238 trkDetItr->directionECAL,
1243 oneCutParameters.
maxDzPV = 100;
1249 oneCutParameters.
maxDzPV = 100;
1257 <<
" nNearTRKs " << nNearTRKs;
1261 edm::LogVerbatim(
"IsoTrack") <<
"coneh " << conehmaxNearP <<
"ok " << trkDetItr->okECAL <<
" "
1262 << trkDetItr->okHCAL;
1266 trkDetItr->pointHCAL,
1267 trkDetItr->pointECAL,
1269 trkDetItr->directionECAL,
1274 trkDetItr->pointHCAL,
1275 trkDetItr->pointECAL,
1277 trkDetItr->directionECAL,
1279 double e_inCone = e2 -
e1;
1280 bool chgIsolFlag = (conehmaxNearP <
cutCharge_);
1281 bool mipFlag = (eMipDR <
cutMip_);
1283 bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1285 selFlags.push_back(selectTk);
1286 selFlags.push_back(qltyFlag);
1287 selFlags.push_back(qltyMissFlag);
1288 selFlags.push_back(qltyPVFlag);
1291 <<
" ; ok: " << trkDetItr->okECAL <<
"/" << trkDetItr->okHCAL
1292 <<
" ; chgiso: " << conehmaxNearP <<
"<" <<
cutCharge_ <<
"(" << chgIsolFlag
1295 if (chgIsolFlag && mipFlag && trkpropFlag) {
1296 double distFromHotCell = -99.0;
1297 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1299 std::vector<DetId> coneRecHitDetIds;
1302 trkDetItr->pointHCAL,
1303 trkDetItr->pointECAL,
1305 trkDetItr->directionHCAL,
1320 t_TrkP->push_back(pTrack->
p());
1341 edm::LogVerbatim(
"IsoTrack") <<
"size of Seeding TripletLayers hb/he " << layershb->
size() <<
"/"
1342 << layershe->
size();
1349 int nSeedHB = 0, nSeedHE = 0;
1358 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1359 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1360 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
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();
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();
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();
1391 reco::VertexCollection::const_iterator vitSel;
1394 for (reco::VertexCollection::const_iterator vit = pVertHE->begin(); vit != pVertHE->end(); vit++) {
1410 reco::VertexCollection::const_iterator vitSel;
1412 bool vtxMatch(
false);
1413 for (reco::VertexCollection::const_iterator vit = pVertHB->begin(); vit != pVertHB->end(); vit++) {
1425 if (
R > 1.2 && vtxMatch)
1430 edm::LogVerbatim(
"IsoTrack") <<
"(HB/HE) nCand: " << nCandHB <<
"/" << nCandHE <<
"nSeed: " << nSeedHB <<
"/"
1446 if (!trkCollection.
isValid()) {
1449 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1450 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
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++) {
1459 double enIn = candref->energyIn();
1460 h_EnIn->Fill(candref->energyIn());
1461 h_EnOut->Fill(candref->energyOut());
1463 candref->track()->px(), candref->track()->py(), candref->track()->pz(), candref->track()->p());
1465 edm::LogVerbatim(
"IsoTrack") <<
"HLT Level Candidate eta/phi/pt/enIn:" << candref->track()->eta() <<
"/"
1466 << candref->track()->phi() <<
"/" << candref->track()->pt() <<
"/"
1467 << candref->energyIn();
1469 double mindR = 999.9, mindP1 = 999.9, eMipDR = 0.0;
1470 std::vector<bool> selFlags;
1472 double conehmaxNearP = 0, hCone = 0;
1473 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++,
nTracks++) {
1474 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1476 double dr =
dR(v1, v2);
1477 double dp1 =
std::abs(1. / v1.r() - 1. / v2.r());
1479 edm::LogVerbatim(
"IsoTrack") <<
"This recotrack(eta/phi/pt) " << pTrack->
eta() <<
"/" << pTrack->
phi() <<
"/"
1480 << pTrack->
pt() <<
" has dr/dp= " <<
dr <<
"/" << dp1;
1485 int nRH_eMipDR = 0, nNearTRKs = 0;
1489 trkDetItr->pointHCAL,
1490 trkDetItr->pointECAL,
1492 trkDetItr->directionECAL,
1497 oneCutParameters.
maxDzPV = 100;
1503 oneCutParameters.
maxDzPV = 100;
1514 trkDetItr->pointHCAL,
1515 trkDetItr->pointECAL,
1517 trkDetItr->directionECAL,
1522 trkDetItr->pointHCAL,
1523 trkDetItr->pointECAL,
1525 trkDetItr->directionECAL,
1527 double e_inCone = e2 -
e1;
1528 bool chgIsolFlag = (conehmaxNearP <
cutCharge_);
1529 bool mipFlag = (eMipDR <
cutMip_);
1531 bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
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;
1544 std::vector<DetId> coneRecHitDetIds;
1547 trkDetItr->pointHCAL,
1548 trkDetItr->pointECAL,
1550 trkDetItr->directionHCAL,
1562 if (mindR >= 0.05) {
1575 std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1579 for (
int j = 0;
j < 3;
j++) {
1580 for (
unsigned int k = 0;
k <
vec_[
j].size();
k++) {
1588 double deta, dphi,
dr;
1590 for (
int lvl = 1; lvl < 3; lvl++) {
1591 for (
unsigned int i = 0;
i <
vec_[lvl].size();
i++) {
1596 edm::LogVerbatim(
"IsoTrack") <<
"lvl " << lvl <<
" i " <<
i <<
" deta " << deta <<
" dphi " << dphi <<
" dR "
1606 for (
unsigned int k = 0;
k <
vec_[2].size(); ++
k) {
1611 <<
vec_[2][
k].phi();
1612 for (
unsigned int j = 0;
j <
vec_[1].size();
j++) {
1616 mindRvec =
vec_[1][
j];
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()
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++) {
1651 edm::LogVerbatim(
"IsoTrack") <<
"reco track: pt " << v4.pt() <<
" eta " << v4.eta() <<
" phi " << v4.phi()
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++) {
1670 edm::LogVerbatim(
"IsoTrack") <<
"Now Matching at Reco level in step 2 DR: " << mindR <<
":" << mindP
1671 <<
" eta:phi " << mindRvec.eta() <<
":" << mindRvec.phi();
1682 if (goodTk != trkCollection->end())
1683 goodTks.push_back(goodTk);
1689 std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1690 if (trkCollection.
isValid()) {
1691 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1695 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1701 <<
" energy " <<
hit->energy();
1707 <<
" energy " <<
hit->energy();
1715 unsigned int nTracks = 0, ngoodTk = 0, nselTk = 0;
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));
1722 double eMipDR = 9999., e_inCone = 0, conehmaxNearP = 0, mindR = 999.9, hCone = 0;
1723 if (trkDetItr->okHCAL) {
1727 for (
unsigned k = 0;
k <
vec_[0].size(); ++
k) {
1733 edm::LogVerbatim(
"IsoTrack") <<
"Track ECAL " << trkDetItr->okECAL <<
" HCAL " << trkDetItr->okHCAL <<
" Flag "
1735 if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL) {
1737 int nRH_eMipDR = 0, nNearTRKs = 0;
1741 trkDetItr->pointHCAL,
1742 trkDetItr->pointECAL,
1744 trkDetItr->directionECAL,
1749 trkDetItr->pointHCAL,
1750 trkDetItr->pointECAL,
1752 trkDetItr->directionECAL,
1757 trkDetItr->pointHCAL,
1758 trkDetItr->pointECAL,
1760 trkDetItr->directionECAL,
1765 double distFromHotCell = -99.0;
1766 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1768 std::vector<DetId> coneRecHitDetIds;
1771 trkDetItr->pointHCAL,
1772 trkDetItr->pointECAL,
1774 trkDetItr->directionHCAL,
1830 double &phiTriggered,
1835 edm::LogVerbatim(
"IsoTrack") <<
"Inside chgIsolation() with eta/phi Triggered: " << etaTriggered <<
"/"
1837 std::vector<double>
maxP;
1839 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1840 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1847 std::vector<std::pair<unsigned int, std::pair<double, double>>> VecSeedsatEC;
1851 bool vtxMatch =
false;
1853 unsigned int ivSel =
recVtxs_->size();
1871 std::pair<double, double> seedCooAtEC;
1886 VecSeedsatEC.push_back(std::make_pair(iS, seedCooAtEC));
1890 for (
unsigned int l = 0;
l < VecSeedsatEC.size();
l++) {
1891 unsigned int iSeed = VecSeedsatEC[
l].first;
1897 for (
unsigned int j = 0;
j < VecSeedsatEC.size();
j++) {
1898 unsigned int iSurr = VecSeedsatEC[
j].first;
1899 if (iSeed != iSurr) {
1907 unsigned int ivSel =
recVtxs_->size();
1908 double minDZ2 = 100;
1919 VecSeedsatEC[
i].second.second,
1920 VecSeedsatEC[
j].second.first,
1921 VecSeedsatEC[
j].second.second);
1933 double conehmaxNearP = -1;
1934 bool selectTk =
false;
1935 double mindR = 999.9;
1938 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++,
nTracks++) {
1940 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1942 double dr =
dR(v1, v2);
1958 std::vector<int> nGood(4, 0);
1959 if (trkCollection.
isValid()) {
1960 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1968 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1970 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1972 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
1975 double ptTriggered(-10), etaTriggered(-100), phiTriggered(-100);
1976 for (
unsigned int p = 0;
p < l1tauobjref.size();
p++) {
1977 if (l1tauobjref[
p]->
pt() > ptTriggered) {
1978 ptTriggered = l1tauobjref[
p]->pt();
1979 phiTriggered = l1tauobjref[
p]->phi();
1980 etaTriggered = l1tauobjref[
p]->eta();
1983 for (
unsigned int p = 0;
p < l1jetobjref.size();
p++) {
1984 if (l1jetobjref[
p]->
pt() > ptTriggered) {
1985 ptTriggered = l1jetobjref[
p]->pt();
1986 phiTriggered = l1jetobjref[
p]->phi();
1987 etaTriggered = l1jetobjref[
p]->eta();
1990 for (
unsigned int p = 0;
p < l1forjetobjref.size();
p++) {
1991 if (l1forjetobjref[
p]->
pt() > ptTriggered) {
1992 ptTriggered = l1forjetobjref[
p]->pt();
1993 phiTriggered = l1forjetobjref[
p]->phi();
1994 etaTriggered = l1forjetobjref[
p]->eta();
1997 double pTriggered = ptTriggered * cosh(etaTriggered);
1999 ptTriggered *
cos(phiTriggered), ptTriggered *
sin(phiTriggered), pTriggered * tanh(etaTriggered), pTriggered);
2001 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
2003 for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++,
nTracks++) {
2004 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
2007 double mindR =
dR(v4, pTrigger);
2009 edm::LogVerbatim(
"IsoTrack") <<
"Track ECAL " << trkDetItr->okECAL <<
" HCAL " << trkDetItr->okHCAL <<
" Flag "
2011 if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL && mindR > 1.0) {
2012 int nRH_eMipDR(0), nNearTRKs(0);
2016 trkDetItr->pointHCAL,
2017 trkDetItr->pointECAL,
2019 trkDetItr->directionECAL,
2021 double conehmaxNearP =
2023 if (conehmaxNearP < 2.0 && eMipDR < 1.0) {
2024 if (pTrack->
p() >= 20 && pTrack->
p() < 30) {
2026 }
else if (pTrack->
p() >= 30 && pTrack->
p() < 40) {
2028 }
else if (pTrack->
p() >= 40 && pTrack->
p() < 60) {
2030 }
else if (pTrack->
p() >= 60 && pTrack->
p() < 100) {
2038 for (
unsigned int ii = 0;
ii < nGood.size(); ++
ii)
2043 h_p[indx]->Fill(vec.r());
2044 h_pt[indx]->Fill(vec.pt());
2045 h_eta[indx]->Fill(vec.eta());
2046 h_phi[indx]->Fill(vec.phi());
2056 h_dEta[indx]->Fill(deta);
2057 h_dPhi[indx]->Fill(dphi);
2058 h_dPt[indx]->Fill(dpt);
2063 edm::LogVerbatim(
"IsoTrack") <<
"mindR for index " << indx <<
" is " <<
dr <<
" deta " << deta <<
" dphi " << dphi
2064 <<
" dpt " << dpt <<
" dinvpt " << dinvpt;
2069 h_eMip[indx]->Fill(eMipDR);
2073 for (
int lim = 0; lim < 5; ++lim) {
2094 if (
kk >= 0 && eMipDR > 0.01 && hCone > 1.0) {
2095 for (
int lim = 0; lim < 5; ++lim) {
2098 h_eCalo[lim][indx][
kk]->Fill(hCone + eMipDR);
2107 double phi1 =
vec1.phi();
2110 double phi2 =
vec2.phi();
2113 double dphi = phi1 - phi2;
2116 else if (dphi < -
M_PI)
2124 return std::sqrt(deta * deta + dphi * dphi);
2134 return ((1 /
vec1.pt()) - (1 /
vec2.pt()));
2139 for (
unsigned int k = 0;
k <
vec_[0].size(); ++
k) {
2150 return std::pair<double, double>(
eta,
phi);
2158 double Rcurv = 9999999;
2160 Rcurv =
pT * 33.3 * 100 / (
bfVal_ * 10);
2162 double ecDist =
zEE_;
2163 double ecRad =
rEB_;
2164 double theta = 2 * atan(
exp(-etaIP));
2168 if (fabs(etaIP) < 1.479) {
2169 if ((0.5 * ecRad / Rcurv) > 1) {
2174 double alpha1 = 2 * asin(0.5 * ecRad / Rcurv);
2177 zNew =
z * (Rcurv * alpha1) / ecRad + vtxZ;
2179 zNew = -
z * (Rcurv * alpha1) / ecRad + vtxZ;
2180 double zAbs = fabs(zNew);
2181 if (zAbs < ecDist) {
2182 etaEC = -
log(
tan(0.5 * atan(ecRad / zAbs)));
2185 if (zAbs > ecDist) {
2186 zAbs = (fabs(etaIP) / etaIP) * ecDist;
2187 double Zflight = fabs(zAbs - vtxZ);
2188 double alpha = (Zflight * ecRad) / (
z * Rcurv);
2189 double Rec = 2 * Rcurv *
sin(
alpha / 2);
2191 etaEC = -
log(
tan(0.5 * atan(Rec / ecDist)));
2195 zNew = (fabs(etaIP) / etaIP) * ecDist;
2196 double Zflight = fabs(zNew - vtxZ);
2197 double Rvirt = fabs(Zflight *
tan(
theta));
2198 double Rec = 2 * Rcurv *
sin(Rvirt / (2 * Rcurv));
2200 etaEC = -
log(
tan(0.5 * atan(Rec / ecDist)));
2212 std::pair<double, double> retVal(etaEC, phiEC);
2218 double theta1 = 2 * atan(
exp(-
eta1));
2219 double theta2 = 2 * atan(
exp(-
eta2));
2220 if (fabs(
eta1) < 1.479)
2222 else if (fabs(
eta1) > 1.479 && fabs(
eta1) < 7.0)