29 #include "TDirectory.h"
114 std::vector<double> &PixMaxP,
123 std::vector<bool> &Flags,
132 double &phiTriggered,
264 std::vector<math::XYZTLorentzVector>
vec_[3];
268 : hltPrescaleProvider_(iConfig, consumesCollector(), *this),
269 trigNames_(iConfig.getUntrackedParameter<
std::vector<
std::
string>>(
"Triggers")),
270 pixCandTag_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"pixCandTag")),
271 l1CandTag_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"l1CandTag")),
272 l2CandTag_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"l2CandTag")),
273 pixelTracksSources_(iConfig.getUntrackedParameter<
std::vector<
edm::
InputTag>>(
"pixelTracksSources")),
274 doL2L3_(iConfig.getUntrackedParameter<
bool>(
"doL2L3",
true)),
275 doTiming_(iConfig.getUntrackedParameter<
bool>(
"doTimingTree",
true)),
276 doMipCutTree_(iConfig.getUntrackedParameter<
bool>(
"doMipCutTree",
true)),
277 doTrkResTree_(iConfig.getUntrackedParameter<
bool>(
"doTrkResTree",
true)),
278 doChgIsolTree_(iConfig.getUntrackedParameter<
bool>(
"doChgIsolTree",
true)),
279 doStudyIsol_(iConfig.getUntrackedParameter<
bool>(
"doStudyIsol",
true)),
280 verbosity_(iConfig.getUntrackedParameter<
int>(
"verbosity", 0)),
281 pixelIsolationConeSizeAtEC_(iConfig.getUntrackedParameter<
std::vector<double>>(
"pixelIsolationConeSizeAtEC")),
282 minPTrackValue_(iConfig.getUntrackedParameter<double>(
"minPTrackValue")),
283 vtxCutSeed_(iConfig.getUntrackedParameter<double>(
"vertexCutSeed")),
284 vtxCutIsol_(iConfig.getUntrackedParameter<double>(
"vertexCutIsol")),
285 tauUnbiasCone_(iConfig.getUntrackedParameter<double>(
"tauUnbiasCone")),
286 prelimCone_(iConfig.getUntrackedParameter<double>(
"prelimCone")),
287 theTrackQuality_(iConfig.getUntrackedParameter<
std::
string>(
"trackQuality",
"highPurity")),
288 processName_(iConfig.getUntrackedParameter<
std::
string>(
"processName",
"HLT")),
289 dr_L1_(iConfig.getUntrackedParameter<double>(
"isolationL1", 1.0)),
290 a_coneR_(iConfig.getUntrackedParameter<double>(
"coneRadius", 34.98)),
291 a_charIsoR_(a_coneR_ + 28.9),
292 a_neutIsoR_(a_charIsoR_ * 0.726),
293 a_mipR_(iConfig.getUntrackedParameter<double>(
"coneRadiusMIP", 14.0)),
294 a_neutR1_(iConfig.getUntrackedParameter<double>(
"coneRadiusNeut1", 21.0)),
295 a_neutR2_(iConfig.getUntrackedParameter<double>(
"coneRadiusNeut2", 29.0)),
296 cutMip_(iConfig.getUntrackedParameter<double>(
"cutMIP", 1.0)),
297 cutCharge_(iConfig.getUntrackedParameter<double>(
"chargeIsolation", 2.0)),
298 cutNeutral_(iConfig.getUntrackedParameter<double>(
"neutralIsolation", 2.0)),
299 minRunNo_(iConfig.getUntrackedParameter<
int>(
"minRun")),
300 maxRunNo_(iConfig.getUntrackedParameter<
int>(
"maxRun")),
302 t_timeL2Prod(nullptr),
308 t_TrkselTkFlag(nullptr),
309 t_TrkqltyFlag(nullptr),
310 t_TrkMissFlag(nullptr),
311 t_TrkPVFlag(nullptr),
312 t_TrkNuIsolFlag(nullptr),
314 t_PixcandPt(nullptr),
315 t_PixcandEta(nullptr),
316 t_PixcandPhi(nullptr),
317 t_PixcandMaxP(nullptr),
318 t_PixTrkcandP(nullptr),
319 t_PixTrkcandPt(nullptr),
320 t_PixTrkcandEta(nullptr),
321 t_PixTrkcandPhi(nullptr),
322 t_PixTrkcandMaxP(nullptr),
323 t_PixTrkcandselTk(nullptr),
326 t_NFcandEta(nullptr),
327 t_NFcandPhi(nullptr),
328 t_NFcandEmip(nullptr),
329 t_NFTrkcandP(nullptr),
330 t_NFTrkcandPt(nullptr),
331 t_NFTrkcandEta(nullptr),
332 t_NFTrkcandPhi(nullptr),
333 t_NFTrkcandEmip(nullptr),
334 t_NFTrkMinDR(nullptr),
335 t_NFTrkMinDP1(nullptr),
336 t_NFTrkselTkFlag(nullptr),
337 t_NFTrkqltyFlag(nullptr),
338 t_NFTrkMissFlag(nullptr),
339 t_NFTrkPVFlag(nullptr),
340 t_NFTrkPropFlag(nullptr),
341 t_NFTrkChgIsoFlag(nullptr),
342 t_NFTrkNeuIsoFlag(nullptr),
343 t_NFTrkMipFlag(nullptr),
363 tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
365 tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel);
407 <<
"\t prelimCone " <<
prelimCone_ <<
"\t pixelIsolationConeSizeAtEC";
411 double pl[] = {20, 30, 40, 60, 80, 120};
412 for (
int i = 0;
i < 6; ++
i)
510 std::vector<std::string> triggers = {
"HLT_IsoTrackHB"};
512 std::vector<double> cones = {35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 63.9, 70.0};
514 desc.
addUntracked<std::vector<std::string>>(
"Triggers", triggers);
547 desc.
addUntracked<std::vector<double>>(
"pixelIsolationConeSizeAtEC", cones);
553 desc.
add<
unsigned int>(
"stageL1Trigger", 1);
554 descriptions.
add(
"isoTrigDefault", desc);
559 edm::LogVerbatim(
"IsoTrack") <<
"Event starts====================================";
561 int RunNo =
iEvent.id().run();
574 if (!triggerEventHandle.
isValid()) {
575 edm::LogWarning(
"IsoTrack") <<
"Error! Can't get the product hltTriggerSummaryAOD";
595 if (!
recVtxs_->empty() && !((*recVtxs_)[0].isFake())) {
612 for (reco::TrackCollection::const_iterator pit = iPixCol->begin(); pit != iPixCol->end(); pit++) {
622 for (
unsigned int ifilter = 0; ifilter <
triggerEvent.sizeFilters(); ++ifilter) {
623 std::string FilterNames[7] = {
"hltL1sL1SingleJet68",
624 "hltIsolPixelTrackL2FilterHE",
625 "ecalIsolPixelTrackFilterHE",
626 "hltIsolPixelTrackL3FilterHE",
627 "hltIsolPixelTrackL2FilterHB",
628 "ecalIsolPixelTrackFilterHB",
629 "hltIsolPixelTrackL3FilterHB"};
631 for (
int i = 0;
i < 7;
i++) {
632 if (
label == FilterNames[
i])
644 <<
iEvent.luminosityBlock() <<
" Bunch " <<
iEvent.bunchCrossing() <<
" mybxlumi "
647 edm::LogWarning(
"IsoTrack") <<
"Error! Can't get the product triggerResults";
651 std::vector<std::string>
modules;
655 const std::vector<std::string> &triggerNames_ =
triggerNames.triggerNames();
660 unsigned int triggerindx =
hltConfig.triggerIndex(triggerNames_[
i]);
661 const std::vector<std::string> &moduleLabels(
hltConfig.moduleLabels(triggerindx));
667 edm::LogVerbatim(
"IsoTrack") <<
"trigger that i want " << triggerNames_[
i] <<
" accept "
684 << triggerNames_[
i] <<
" accept " <<
hlt <<
" preL1 " << preL1 <<
" preHLT " << preHLT;
685 for (
int iv = 0; iv < 3; ++iv)
690 trigList_.insert(std::pair<unsigned int, unsigned int>(RunNo, 1));
694 for (
unsigned int ifilter = 0; ifilter <
triggerEvent.sizeFilters(); ++ifilter) {
695 std::vector<int>
Keys;
698 for (
unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
699 if (
label.find(moduleLabels[imodule]) != std::string::npos) {
702 for (
unsigned int ifiltrKey = 0; ifiltrKey <
triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
706 if (
label.find(
"L2Filter") != std::string::npos) {
707 vec_[1].push_back(v4);
708 }
else if (
label.find(
"L3Filter") != std::string::npos) {
709 vec_[2].push_back(v4);
711 vec_[0].push_back(v4);
716 <<
" phi " << TO.
phi() <<
" mass " << TO.
mass() <<
" Id " << TO.
id();
721 std::vector<reco::TrackCollection::const_iterator> goodTks;
727 reco::TrackCollection::const_iterator trkItr;
728 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++)
729 goodTks.push_back(trkItr);
769 for (
unsigned itrig = 0; itrig < triggerNames_.size(); itrig++) {
770 unsigned int triggerindx =
hltConfig.triggerIndex(triggerNames_[itrig]);
771 if (triggerindx >=
n)
772 edm::LogVerbatim(
"IsoTrack") << triggerNames_[itrig] <<
" " << triggerindx <<
" does not exist in"
773 <<
" the current menu";
775 edm::LogVerbatim(
"IsoTrack") << triggerNames_[itrig] <<
" " << triggerindx <<
" exists";
826 std::vector<double> &PixMaxP,
848 std::vector<bool> &Flags,
874 char hname[100], htit[100];
876 "L2Match",
"L2NoMatch",
"L3Match",
"L3NoMatch",
"HLTTrk",
"HLTGoodTrk",
877 "HLTIsoTrk",
"HLTMip",
"HLTSelect",
"nonHLTTrk",
"nonHLTGoodTrk",
"nonHLTIsoTrk",
878 "nonHLTMip",
"nonHLTSelect"};
894 t_TrkP =
new std::vector<double>();
957 t_ECone =
new std::vector<double>();
982 std::string FilterNames[7] = {
"hltL1sL1SingleJet68",
983 "hltIsolPixelTrackL2FilterHE",
984 "ecalIsolPixelTrackFilterHE",
985 "hltIsolPixelTrackL3FilterHE",
986 "hltIsolPixelTrackL2FilterHB",
987 "ecalIsolPixelTrackFilterHB",
988 "hltIsolPixelTrackL3FilterHB"};
989 for (
int i = 0;
i < 7;
i++) {
990 h_Filters->GetXaxis()->SetBinLabel(
i + 1, FilterNames[
i].c_str());
993 h_nHLT =
fs_->
make<TH1I>(
"h_nHLT",
"Size of trigger Names", 1000, 1, 1000);
994 h_HLT =
fs_->
make<TH1I>(
"h_HLT",
"HLT accept", 3, -1, 2);
997 h_Pre =
fs_->
make<TH1I>(
"h_Pre",
"Prescale", 3000, 0, 3000);
999 h_PreL1wt =
fs_->
make<TH1D>(
"h_PreL1wt",
"Weighted L1 Prescale", 500, 0, 500);
1000 h_PreHLTwt =
fs_->
make<TH1D>(
"h_PreHLTwt",
"Weighted HLT Prescale", 500, 0, 100);
1003 h_EnIn =
fs_->
make<TH1D>(
"h_EnInEcal",
"EnergyIn Ecal", 200, 0.0, 20.0);
1004 h_EnOut =
fs_->
make<TH1D>(
"h_EnOutEcal",
"EnergyOut Ecal", 200, 0.0, 20.0);
1006 fs_->
make<TH2D>(
"h_MipEnMatch",
"MipEn: HLT level vs Reco Level (Matched)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1008 "h_MipEnNoMatch",
"MipEn: HLT level vs Reco Level (No Match Found)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1011 h_nL3Objs =
fs_->
make<TH1I>(
"h_nL3Objs",
"Number of L3 objects", 10, 0, 10);
1020 "NewFilterRecoMatch",
1021 "NewFilterRecoNoMatch"};
1022 for (
int ipair = 0; ipair < 9; ipair++) {
1023 sprintf(hname,
"h_dEta%s", pairs[ipair].c_str());
1024 sprintf(htit,
"#Delta#eta for %s", pairs[ipair].c_str());
1025 h_dEta[ipair] =
fs_->
make<TH1D>(hname, htit, 200, -10.0, 10.0);
1026 h_dEta[ipair]->GetXaxis()->SetTitle(
"d#eta");
1028 sprintf(hname,
"h_dPhi%s", pairs[ipair].c_str());
1029 sprintf(htit,
"#Delta#phi for %s", pairs[ipair].c_str());
1030 h_dPhi[ipair] =
fs_->
make<TH1D>(hname, htit, 140, -7.0, 7.0);
1031 h_dPhi[ipair]->GetXaxis()->SetTitle(
"d#phi");
1033 sprintf(hname,
"h_dPt%s", pairs[ipair].c_str());
1034 sprintf(htit,
"#Delta dp_{T} for %s objects", pairs[ipair].c_str());
1035 h_dPt[ipair] =
fs_->
make<TH1D>(hname, htit, 400, -200.0, 200.0);
1036 h_dPt[ipair]->GetXaxis()->SetTitle(
"dp_{T} (GeV)");
1038 sprintf(hname,
"h_dP%s", pairs[ipair].c_str());
1039 sprintf(htit,
"#Delta p for %s objects", pairs[ipair].c_str());
1040 h_dP[ipair] =
fs_->
make<TH1D>(hname, htit, 400, -200.0, 200.0);
1041 h_dP[ipair]->GetXaxis()->SetTitle(
"dP (GeV)");
1043 sprintf(hname,
"h_dinvPt%s", pairs[ipair].c_str());
1044 sprintf(htit,
"#Delta (1/p_{T}) for %s objects", pairs[ipair].c_str());
1046 h_dinvPt[ipair]->GetXaxis()->SetTitle(
"d(1/p_{T})");
1047 sprintf(hname,
"h_mindR%s", pairs[ipair].c_str());
1048 sprintf(htit,
"min(#Delta R) for %s objects", pairs[ipair].c_str());
1050 h_mindR[ipair]->GetXaxis()->SetTitle(
"dR");
1053 for (
int lvl = 0; lvl < 2; lvl++) {
1054 sprintf(hname,
"h_dEtaL1%s",
levels[lvl + 1].c_str());
1055 sprintf(htit,
"#Delta#eta for L1 and %s objects",
levels[lvl + 1].c_str());
1058 sprintf(hname,
"h_dPhiL1%s",
levels[lvl + 1].c_str());
1059 sprintf(htit,
"#Delta#phi for L1 and %s objects",
levels[lvl + 1].c_str());
1062 sprintf(hname,
"h_dRL1%s",
levels[lvl + 1].c_str());
1063 sprintf(htit,
"#Delta R for L1 and %s objects",
levels[lvl + 1].c_str());
1068 int levmin = (
doL2L3_ ? 0 : 10);
1069 for (
int ilevel = levmin; ilevel < 20; ilevel++) {
1070 sprintf(hname,
"h_p%s",
levels[ilevel].c_str());
1071 sprintf(htit,
"p for %s objects",
levels[ilevel].c_str());
1072 h_p[ilevel] =
fs_->
make<TH1D>(hname, htit, 100, 0.0, 500.0);
1073 h_p[ilevel]->GetXaxis()->SetTitle(
"p (GeV)");
1075 sprintf(hname,
"h_pt%s",
levels[ilevel].c_str());
1076 sprintf(htit,
"p_{T} for %s objects",
levels[ilevel].c_str());
1077 h_pt[ilevel] =
fs_->
make<TH1D>(hname, htit, 100, 0.0, 500.0);
1078 h_pt[ilevel]->GetXaxis()->SetTitle(
"p_{T} (GeV)");
1080 sprintf(hname,
"h_eta%s",
levels[ilevel].c_str());
1081 sprintf(htit,
"#eta for %s objects",
levels[ilevel].c_str());
1082 h_eta[ilevel] =
fs_->
make<TH1D>(hname, htit, 100, -5.0, 5.0);
1083 h_eta[ilevel]->GetXaxis()->SetTitle(
"#eta");
1085 sprintf(hname,
"h_phi%s",
levels[ilevel].c_str());
1086 sprintf(htit,
"#phi for %s objects",
levels[ilevel].c_str());
1087 h_phi[ilevel] =
fs_->
make<TH1D>(hname, htit, 70, -3.5, 3.50);
1088 h_phi[ilevel]->GetXaxis()->SetTitle(
"#phi");
1093 for (
int icut = 0; icut < 2; icut++) {
1094 sprintf(hname,
"h_eMip%s",
cuts[icut].c_str());
1095 sprintf(htit,
"eMip for %s tracks",
cuts[icut].c_str());
1096 h_eMip[icut] =
fs_->
make<TH1D>(hname, htit, 200, 0.0, 10.0);
1097 h_eMip[icut]->GetXaxis()->SetTitle(
"E_{Mip} (GeV)");
1099 sprintf(hname,
"h_eMaxNearP%s",
cuts[icut].c_str());
1100 sprintf(htit,
"eMaxNearP for %s tracks",
cuts[icut].c_str());
1102 h_eMaxNearP[icut]->GetXaxis()->SetTitle(
"E_{MaxNearP} (GeV)");
1104 sprintf(hname,
"h_eNeutIso%s",
cuts[icut].c_str());
1105 sprintf(htit,
"eNeutIso for %s ",
cuts[icut].c_str());
1107 h_eNeutIso[icut]->GetXaxis()->SetTitle(
"E_{NeutIso} (GeV)");
1109 for (
int kcut = 0; kcut < 2; ++kcut) {
1110 for (
int lim = 0; lim < 5; ++lim) {
1111 sprintf(hname,
"h_etaCalibTracks%sCut%dLim%d",
cuts[icut].c_str(), kcut, lim);
1113 "#eta for %s isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1117 cuts2[kcut].c_str());
1121 sprintf(hname,
"h_etaMipTracks%sCut%dLim%d",
cuts[icut].c_str(), kcut, lim);
1123 "#eta for %s charge isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1127 cuts2[kcut].c_str());
1134 std::string ecut1[3] = {
"all",
"HLTMatched",
"HLTNotMatched"};
1136 int etac[48] = {-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, -13, -14, -15, -16,
1137 -17, -18, -19, -20, -21, -22, -23, -24, 1, 2, 3, 4, 5, 6, 7, 8,
1138 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24};
1139 for (
int icut = 0; icut < 6; icut++) {
1141 int i1 = (icut > 2 ? 1 : 0);
1142 int i2 = icut -
i1 * 3;
1143 for (
int kcut = 0; kcut < 48; kcut++) {
1144 for (
int lim = 0; lim < 5; ++lim) {
1145 sprintf(hname,
"h_eta%dEnHcal%s%s%d", etac[kcut], ecut1[
i2].c_str(), ecut2[
i1].c_str(), lim);
1147 "HCAL energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1153 h_eHcal[lim][icut][kcut] =
fs_->
make<TH1D>(hname, htit, 750, 0.0, 150.0);
1154 h_eHcal[lim][icut][kcut]->GetXaxis()->SetTitle(
"Energy (GeV)");
1155 sprintf(hname,
"h_eta%dEnCalo%s%s%d", etac[kcut], ecut1[
i2].c_str(), ecut2[
i1].c_str(), lim);
1157 "Calorimter energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1163 h_eCalo[lim][icut][kcut] =
fs_->
make<TH1D>(hname, htit, 750, 0.0, 150.0);
1164 h_eCalo[lim][icut][kcut]->GetXaxis()->SetTitle(
"Energy (GeV)");
1172 unsigned int preL1, preHLT;
1173 std::map<unsigned int, unsigned int>::iterator
itr;
1174 std::map<unsigned int, const std::pair<int, int>>::iterator itrPre;
1185 preL1 = (itrPre->second).
first;
1186 preHLT = (itrPre->second).
second;
1187 edm::LogVerbatim(
"IsoTrack") <<
itr->first <<
" " <<
itr->second <<
" " << itrPre->first <<
" " << preL1 <<
" "
1192 g_Pre->Fill(
itr->first, preL1 * preHLT);
1210 if (!trkCollection.
isValid()) {
1213 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1216 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1219 int nRH_eMipDR = 0, nNearTRKs = 0;
1220 std::vector<bool> selFlags;
1221 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++,
nTracks++) {
1222 double conehmaxNearP = 0, hCone = 0, eMipDR = 0.0;
1223 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1226 if (pTrack->
p() > 20) {
1231 trkDetItr->pointHCAL,
1232 trkDetItr->pointECAL,
1234 trkDetItr->directionECAL,
1239 oneCutParameters.
maxDzPV = 100;
1245 oneCutParameters.
maxDzPV = 100;
1253 <<
" nNearTRKs " << nNearTRKs;
1257 edm::LogVerbatim(
"IsoTrack") <<
"coneh " << conehmaxNearP <<
"ok " << trkDetItr->okECAL <<
" "
1258 << trkDetItr->okHCAL;
1262 trkDetItr->pointHCAL,
1263 trkDetItr->pointECAL,
1265 trkDetItr->directionECAL,
1270 trkDetItr->pointHCAL,
1271 trkDetItr->pointECAL,
1273 trkDetItr->directionECAL,
1275 double e_inCone = e2 -
e1;
1276 bool chgIsolFlag = (conehmaxNearP <
cutCharge_);
1277 bool mipFlag = (eMipDR <
cutMip_);
1279 bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1281 selFlags.push_back(selectTk);
1282 selFlags.push_back(qltyFlag);
1283 selFlags.push_back(qltyMissFlag);
1284 selFlags.push_back(qltyPVFlag);
1287 <<
" ; ok: " << trkDetItr->okECAL <<
"/" << trkDetItr->okHCAL
1288 <<
" ; chgiso: " << conehmaxNearP <<
"<" <<
cutCharge_ <<
"(" << chgIsolFlag
1291 if (chgIsolFlag && mipFlag && trkpropFlag) {
1292 double distFromHotCell = -99.0;
1293 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1295 std::vector<DetId> coneRecHitDetIds;
1298 trkDetItr->pointHCAL,
1299 trkDetItr->pointECAL,
1301 trkDetItr->directionHCAL,
1316 t_TrkP->push_back(pTrack->
p());
1337 edm::LogVerbatim(
"IsoTrack") <<
"size of Seeding TripletLayers hb/he " << layershb->
size() <<
"/"
1338 << layershe->
size();
1345 int nSeedHB = 0, nSeedHE = 0;
1354 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1355 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1356 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
1362 double ptTriggered = -10;
1363 double etaTriggered = -100;
1364 double phiTriggered = -100;
1365 for (
unsigned int p = 0;
p < l1tauobjref.size();
p++) {
1366 if (l1tauobjref[
p]->
pt() > ptTriggered) {
1367 ptTriggered = l1tauobjref[
p]->pt();
1368 phiTriggered = l1tauobjref[
p]->phi();
1369 etaTriggered = l1tauobjref[
p]->eta();
1372 for (
unsigned int p = 0;
p < l1jetobjref.size();
p++) {
1373 if (l1jetobjref[
p]->
pt() > ptTriggered) {
1374 ptTriggered = l1jetobjref[
p]->pt();
1375 phiTriggered = l1jetobjref[
p]->phi();
1376 etaTriggered = l1jetobjref[
p]->eta();
1379 for (
unsigned int p = 0;
p < l1forjetobjref.size();
p++) {
1380 if (l1forjetobjref[
p]->
pt() > ptTriggered) {
1381 ptTriggered = l1forjetobjref[
p]->pt();
1382 phiTriggered = l1forjetobjref[
p]->phi();
1383 etaTriggered = l1forjetobjref[
p]->eta();
1387 reco::VertexCollection::const_iterator vitSel;
1390 for (reco::VertexCollection::const_iterator vit = pVertHE->begin(); vit != pVertHE->end(); vit++) {
1406 reco::VertexCollection::const_iterator vitSel;
1408 bool vtxMatch(
false);
1409 for (reco::VertexCollection::const_iterator vit = pVertHB->begin(); vit != pVertHB->end(); vit++) {
1421 if (
R > 1.2 && vtxMatch)
1426 edm::LogVerbatim(
"IsoTrack") <<
"(HB/HE) nCand: " << nCandHB <<
"/" << nCandHE <<
"nSeed: " << nSeedHB <<
"/"
1442 if (!trkCollection.
isValid()) {
1445 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1448 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1451 edm::LogVerbatim(
"IsoTrack") <<
"Number of L2cands:" << L2cands->size() <<
" to be matched to something out of "
1452 << trkCaloDirections1.size() <<
" reco tracks";
1453 for (
unsigned int i = 0;
i < L2cands->size();
i++) {
1456 double enIn = candref->energyIn();
1457 h_EnIn->Fill(candref->energyIn());
1458 h_EnOut->Fill(candref->energyOut());
1460 candref->track()->px(), candref->track()->py(), candref->track()->pz(), candref->track()->p());
1462 edm::LogVerbatim(
"IsoTrack") <<
"HLT Level Candidate eta/phi/pt/enIn:" << candref->track()->eta() <<
"/"
1463 << candref->track()->phi() <<
"/" << candref->track()->pt() <<
"/"
1464 << candref->energyIn();
1466 double mindR = 999.9, mindP1 = 999.9, eMipDR = 0.0;
1467 std::vector<bool> selFlags;
1469 double conehmaxNearP = 0, hCone = 0;
1470 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++,
nTracks++) {
1471 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1473 double dr =
dR(v1, v2);
1474 double dp1 =
std::abs(1. / v1.r() - 1. / v2.r());
1476 edm::LogVerbatim(
"IsoTrack") <<
"This recotrack(eta/phi/pt) " << pTrack->
eta() <<
"/" << pTrack->
phi() <<
"/"
1477 << pTrack->
pt() <<
" has dr/dp= " <<
dr <<
"/" << dp1;
1482 int nRH_eMipDR = 0, nNearTRKs = 0;
1486 trkDetItr->pointHCAL,
1487 trkDetItr->pointECAL,
1489 trkDetItr->directionECAL,
1494 oneCutParameters.
maxDzPV = 100;
1500 oneCutParameters.
maxDzPV = 100;
1511 trkDetItr->pointHCAL,
1512 trkDetItr->pointECAL,
1514 trkDetItr->directionECAL,
1519 trkDetItr->pointHCAL,
1520 trkDetItr->pointECAL,
1522 trkDetItr->directionECAL,
1524 double e_inCone = e2 -
e1;
1525 bool chgIsolFlag = (conehmaxNearP <
cutCharge_);
1526 bool mipFlag = (eMipDR <
cutMip_);
1528 bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1530 selFlags.push_back(selectTk);
1531 selFlags.push_back(qltyFlag);
1532 selFlags.push_back(qltyMissFlag);
1533 selFlags.push_back(qltyPVFlag);
1534 selFlags.push_back(trkpropFlag);
1535 selFlags.push_back(chgIsolFlag);
1536 selFlags.push_back(neuIsolFlag);
1537 selFlags.push_back(mipFlag);
1538 double distFromHotCell = -99.0;
1539 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1541 std::vector<DetId> coneRecHitDetIds;
1544 trkDetItr->pointHCAL,
1545 trkDetItr->pointECAL,
1547 trkDetItr->directionHCAL,
1559 if (mindR >= 0.05) {
1572 std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1576 for (
int j = 0;
j < 3;
j++) {
1577 for (
unsigned int k = 0;
k <
vec_[
j].size();
k++) {
1585 double deta, dphi,
dr;
1587 for (
int lvl = 1; lvl < 3; lvl++) {
1588 for (
unsigned int i = 0;
i <
vec_[lvl].size();
i++) {
1593 edm::LogVerbatim(
"IsoTrack") <<
"lvl " << lvl <<
" i " <<
i <<
" deta " << deta <<
" dphi " << dphi <<
" dR "
1603 for (
unsigned int k = 0;
k <
vec_[2].size(); ++
k) {
1608 <<
vec_[2][
k].phi();
1609 for (
unsigned int j = 0;
j <
vec_[1].size();
j++) {
1613 mindRvec =
vec_[1][
j];
1630 edm::LogVerbatim(
"IsoTrack") <<
"Now Matching L3 track with reco: L3 Track (eta, phi) " <<
vec_[2][
k].eta() <<
":"
1631 <<
vec_[2][
k].phi() <<
" L2 Track " << mindRvec.eta() <<
":" << mindRvec.phi()
1633 reco::TrackCollection::const_iterator goodTk = trkCollection->end();
1634 if (trkCollection.
isValid()) {
1635 double mindP(9999.9);
1636 reco::TrackCollection::const_iterator trkItr;
1637 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++) {
1648 edm::LogVerbatim(
"IsoTrack") <<
"reco track: pt " << v4.pt() <<
" eta " << v4.eta() <<
" phi " << v4.phi()
1652 edm::LogVerbatim(
"IsoTrack") <<
"Now Matching at Reco level in step 1 DR: " << mindR <<
":" << mindP
1653 <<
" eta:phi " << mindRvec.eta() <<
":" << mindRvec.phi();
1654 if (mindR < 0.03 && mindP > 0.1) {
1655 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++) {
1667 edm::LogVerbatim(
"IsoTrack") <<
"Now Matching at Reco level in step 2 DR: " << mindR <<
":" << mindP
1668 <<
" eta:phi " << mindRvec.eta() <<
":" << mindRvec.phi();
1679 if (goodTk != trkCollection->end())
1680 goodTks.push_back(goodTk);
1686 std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1687 if (trkCollection.
isValid()) {
1691 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1694 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1700 <<
" energy " <<
hit->energy();
1706 <<
" energy " <<
hit->energy();
1714 unsigned int nTracks = 0, ngoodTk = 0, nselTk = 0;
1716 for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++,
nTracks++) {
1717 bool l3Track = (
std::find(goodTks.begin(), goodTks.end(), trkDetItr->trkItr) != goodTks.end());
1718 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1721 double eMipDR = 9999., e_inCone = 0, conehmaxNearP = 0, mindR = 999.9, hCone = 0;
1722 if (trkDetItr->okHCAL) {
1726 for (
unsigned k = 0;
k <
vec_[0].size(); ++
k) {
1732 edm::LogVerbatim(
"IsoTrack") <<
"Track ECAL " << trkDetItr->okECAL <<
" HCAL " << trkDetItr->okHCAL <<
" Flag "
1734 if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL) {
1736 int nRH_eMipDR = 0, nNearTRKs = 0;
1740 trkDetItr->pointHCAL,
1741 trkDetItr->pointECAL,
1743 trkDetItr->directionECAL,
1748 trkDetItr->pointHCAL,
1749 trkDetItr->pointECAL,
1751 trkDetItr->directionECAL,
1756 trkDetItr->pointHCAL,
1757 trkDetItr->pointECAL,
1759 trkDetItr->directionECAL,
1764 double distFromHotCell = -99.0;
1765 int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1767 std::vector<DetId> coneRecHitDetIds;
1770 trkDetItr->pointHCAL,
1771 trkDetItr->pointECAL,
1773 trkDetItr->directionHCAL,
1829 double &phiTriggered,
1834 edm::LogVerbatim(
"IsoTrack") <<
"Inside chgIsolation() with eta/phi Triggered: " << etaTriggered <<
"/"
1836 std::vector<double>
maxP;
1838 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1841 std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1848 std::vector<std::pair<unsigned int, std::pair<double, double>>> VecSeedsatEC;
1852 bool vtxMatch =
false;
1854 unsigned int ivSel =
recVtxs_->size();
1856 for (
unsigned int iv = 0; iv <
recVtxs_->size(); ++iv) {
1872 std::pair<double, double> seedCooAtEC;
1887 VecSeedsatEC.push_back(std::make_pair(iS, seedCooAtEC));
1891 for (
unsigned int l = 0;
l < VecSeedsatEC.size();
l++) {
1892 unsigned int iSeed = VecSeedsatEC[
l].first;
1898 for (
unsigned int j = 0;
j < VecSeedsatEC.size();
j++) {
1899 unsigned int iSurr = VecSeedsatEC[
j].first;
1900 if (iSeed != iSurr) {
1908 unsigned int ivSel =
recVtxs_->size();
1909 double minDZ2 = 100;
1910 for (
unsigned int iv = 0; iv <
recVtxs_->size(); ++iv) {
1920 VecSeedsatEC[
i].second.second,
1921 VecSeedsatEC[
j].second.first,
1922 VecSeedsatEC[
j].second.second);
1934 double conehmaxNearP = -1;
1935 bool selectTk =
false;
1936 double mindR = 999.9;
1939 for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++,
nTracks++) {
1941 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1943 double dr =
dR(v1, v2);
1959 std::vector<int> nGood(4, 0);
1960 if (trkCollection.
isValid()) {
1964 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1971 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1973 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1975 std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
1978 double ptTriggered(-10), etaTriggered(-100), phiTriggered(-100);
1979 for (
unsigned int p = 0;
p < l1tauobjref.size();
p++) {
1980 if (l1tauobjref[
p]->
pt() > ptTriggered) {
1981 ptTriggered = l1tauobjref[
p]->pt();
1982 phiTriggered = l1tauobjref[
p]->phi();
1983 etaTriggered = l1tauobjref[
p]->eta();
1986 for (
unsigned int p = 0;
p < l1jetobjref.size();
p++) {
1987 if (l1jetobjref[
p]->
pt() > ptTriggered) {
1988 ptTriggered = l1jetobjref[
p]->pt();
1989 phiTriggered = l1jetobjref[
p]->phi();
1990 etaTriggered = l1jetobjref[
p]->eta();
1993 for (
unsigned int p = 0;
p < l1forjetobjref.size();
p++) {
1994 if (l1forjetobjref[
p]->
pt() > ptTriggered) {
1995 ptTriggered = l1forjetobjref[
p]->pt();
1996 phiTriggered = l1forjetobjref[
p]->phi();
1997 etaTriggered = l1forjetobjref[
p]->eta();
2000 double pTriggered = ptTriggered * cosh(etaTriggered);
2002 ptTriggered *
cos(phiTriggered), ptTriggered *
sin(phiTriggered), pTriggered * tanh(etaTriggered), pTriggered);
2004 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
2006 for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++,
nTracks++) {
2007 const reco::Track *pTrack = &(*(trkDetItr->trkItr));
2010 double mindR =
dR(v4, pTrigger);
2012 edm::LogVerbatim(
"IsoTrack") <<
"Track ECAL " << trkDetItr->okECAL <<
" HCAL " << trkDetItr->okHCAL <<
" Flag "
2014 if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL && mindR > 1.0) {
2015 int nRH_eMipDR(0), nNearTRKs(0);
2019 trkDetItr->pointHCAL,
2020 trkDetItr->pointECAL,
2022 trkDetItr->directionECAL,
2024 double conehmaxNearP =
2026 if (conehmaxNearP < 2.0 && eMipDR < 1.0) {
2027 if (pTrack->
p() >= 20 && pTrack->
p() < 30) {
2029 }
else if (pTrack->
p() >= 30 && pTrack->
p() < 40) {
2031 }
else if (pTrack->
p() >= 40 && pTrack->
p() < 60) {
2033 }
else if (pTrack->
p() >= 60 && pTrack->
p() < 100) {
2041 for (
unsigned int ii = 0;
ii < nGood.size(); ++
ii)
2046 h_p[indx]->Fill(vec.r());
2047 h_pt[indx]->Fill(vec.pt());
2048 h_eta[indx]->Fill(vec.eta());
2049 h_phi[indx]->Fill(vec.phi());
2059 h_dEta[indx]->Fill(deta);
2060 h_dPhi[indx]->Fill(dphi);
2061 h_dPt[indx]->Fill(dpt);
2066 edm::LogVerbatim(
"IsoTrack") <<
"mindR for index " << indx <<
" is " <<
dr <<
" deta " << deta <<
" dphi " << dphi
2067 <<
" dpt " << dpt <<
" dinvpt " << dinvpt;
2072 h_eMip[indx]->Fill(eMipDR);
2076 for (
int lim = 0; lim < 5; ++lim) {
2097 if (
kk >= 0 && eMipDR > 0.01 && hCone > 1.0) {
2098 for (
int lim = 0; lim < 5; ++lim) {
2101 h_eCalo[lim][indx][
kk]->Fill(hCone + eMipDR);
2110 double phi1 =
vec1.phi();
2113 double phi2 =
vec2.phi();
2116 double dphi = phi1 - phi2;
2119 else if (dphi < -
M_PI)
2127 return std::sqrt(deta * deta + dphi * dphi);
2137 return ((1 /
vec1.pt()) - (1 /
vec2.pt()));
2142 for (
unsigned int k = 0;
k <
vec_[0].size(); ++
k) {
2153 return std::pair<double, double>(
eta,
phi);
2161 double Rcurv = 9999999;
2163 Rcurv =
pT * 33.3 * 100 / (
bfVal_ * 10);
2165 double ecDist =
zEE_;
2166 double ecRad =
rEB_;
2167 double theta = 2 * atan(
exp(-etaIP));
2171 if (fabs(etaIP) < 1.479) {
2172 if ((0.5 * ecRad / Rcurv) > 1) {
2177 double alpha1 = 2 * asin(0.5 * ecRad / Rcurv);
2180 zNew =
z * (Rcurv * alpha1) / ecRad + vtxZ;
2182 zNew = -
z * (Rcurv * alpha1) / ecRad + vtxZ;
2183 double zAbs = fabs(zNew);
2184 if (zAbs < ecDist) {
2185 etaEC = -
log(
tan(0.5 * atan(ecRad / zAbs)));
2188 if (zAbs > ecDist) {
2189 zAbs = (fabs(etaIP) / etaIP) * ecDist;
2190 double Zflight = fabs(zAbs - vtxZ);
2191 double alpha = (Zflight * ecRad) / (
z * Rcurv);
2192 double Rec = 2 * Rcurv *
sin(
alpha / 2);
2194 etaEC = -
log(
tan(0.5 * atan(Rec / ecDist)));
2198 zNew = (fabs(etaIP) / etaIP) * ecDist;
2199 double Zflight = fabs(zNew - vtxZ);
2200 double Rvirt = fabs(Zflight *
tan(
theta));
2201 double Rec = 2 * Rcurv *
sin(Rvirt / (2 * Rcurv));
2203 etaEC = -
log(
tan(0.5 * atan(Rec / ecDist)));
2215 std::pair<double, double> retVal(etaEC, phiEC);
2221 double theta1 = 2 * atan(
exp(-
eta1));
2222 double theta2 = 2 * atan(
exp(-
eta2));
2223 if (fabs(
eta1) < 1.479)
2225 else if (fabs(
eta1) > 1.479 && fabs(
eta1) < 7.0)