13 #include "TDirectory.h"
15 #include "TLorentzVector.h"
16 #include "TInterpreter.h"
100 std::array<int, 3>
fillTree(std::vector<math::XYZTLorentzVector>& vecL1,
101 std::vector<math::XYZTLorentzVector>& vecL3,
103 std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
104 std::vector<spr::propagatedTrackID>& trkCaloDets,
123 const std::vector<DetId>& ids,
124 std::vector<double>& edet,
126 std::vector<unsigned int>* detIds,
127 std::vector<double>* hitEnergies);
197 : trigNames_(iConfig.getParameter<
std::vector<
std::
string>>(
"triggers")),
198 theTrackQuality_(iConfig.getParameter<
std::
string>(
"trackQuality")),
199 processName_(iConfig.getParameter<
std::
string>(
"processName")),
200 l1Filter_(iConfig.getParameter<
std::
string>(
"l1Filter")),
201 l2Filter_(iConfig.getParameter<
std::
string>(
"l2Filter")),
202 l3Filter_(iConfig.getParameter<
std::
string>(
"l3Filter")),
203 a_coneR_(iConfig.getParameter<double>(
"coneRadius")),
204 a_mipR_(iConfig.getParameter<double>(
"coneRadiusMIP")),
205 a_mipR2_(iConfig.getParameter<double>(
"coneRadiusMIP2")),
206 a_mipR3_(iConfig.getParameter<double>(
"coneRadiusMIP3")),
207 a_mipR4_(iConfig.getParameter<double>(
"coneRadiusMIP4")),
208 a_mipR5_(iConfig.getParameter<double>(
"coneRadiusMIP5")),
209 pTrackMin_(iConfig.getParameter<double>(
"minimumTrackP")),
210 eEcalMax_(iConfig.getParameter<double>(
"maximumEcalEnergy")),
211 maxRestrictionP_(iConfig.getParameter<double>(
"maxTrackP")),
212 slopeRestrictionP_(iConfig.getParameter<double>(
"slopeTrackP")),
213 hcalScale_(iConfig.getUntrackedParameter<double>(
"hHcalScale", 1.0)),
214 eIsolate1_(iConfig.getParameter<double>(
"isolationEnergyTight")),
215 eIsolate2_(iConfig.getParameter<double>(
"isolationEnergyLoose")),
216 pTrackLow_(iConfig.getParameter<double>(
"momentumLow")),
217 pTrackHigh_(iConfig.getParameter<double>(
"momentumHigh")),
218 prescaleLow_(iConfig.getParameter<
int>(
"prescaleLow")),
219 prescaleHigh_(iConfig.getParameter<
int>(
"prescaleHigh")),
220 useRaw_(iConfig.getUntrackedParameter<
int>(
"useRaw", 0)),
221 dataType_(iConfig.getUntrackedParameter<
int>(
"dataType", 0)),
222 mode_(iConfig.getUntrackedParameter<
int>(
"outMode", 11)),
223 ignoreTrigger_(iConfig.getUntrackedParameter<
bool>(
"ignoreTriggers",
false)),
224 useL1Trigger_(iConfig.getUntrackedParameter<
bool>(
"useL1Trigger",
false)),
225 unCorrect_(iConfig.getUntrackedParameter<
bool>(
"unCorrect",
false)),
226 collapseDepth_(iConfig.getUntrackedParameter<
bool>(
"collapseDepth",
false)),
227 hitEthrEB_(iConfig.getParameter<double>(
"EBHitEnergyThreshold")),
228 hitEthrEE0_(iConfig.getParameter<double>(
"EEHitEnergyThreshold0")),
229 hitEthrEE1_(iConfig.getParameter<double>(
"EEHitEnergyThreshold1")),
230 hitEthrEE2_(iConfig.getParameter<double>(
"EEHitEnergyThreshold2")),
231 hitEthrEE3_(iConfig.getParameter<double>(
"EEHitEnergyThreshold3")),
232 hitEthrEELo_(iConfig.getParameter<double>(
"EEHitEnergyThresholdLow")),
233 hitEthrEEHi_(iConfig.getParameter<double>(
"EEHitEnergyThresholdHigh")),
234 triggerEvent_(iConfig.getParameter<
edm::
InputTag>(
"labelTriggerEvent")),
235 theTriggerResultsLabel_(iConfig.getParameter<
edm::
InputTag>(
"labelTriggerResult")),
236 labelGenTrack_(iConfig.getParameter<
std::
string>(
"labelTrack")),
237 labelRecVtx_(iConfig.getParameter<
std::
string>(
"labelVertex")),
238 labelEB_(iConfig.getParameter<
std::
string>(
"labelEBRecHit")),
239 labelEE_(iConfig.getParameter<
std::
string>(
"labelEERecHit")),
240 labelHBHE_(iConfig.getParameter<
std::
string>(
"labelHBHERecHit")),
241 labelTower_(iConfig.getParameter<
std::
string>(
"labelCaloTower")),
242 l1TrigName_(iConfig.getUntrackedParameter<
std::
string>(
"l1TrigName",
"L1_SingleJet60")),
243 oldID_(iConfig.getUntrackedParameter<
std::vector<
int>>(
"oldID")),
244 newDepth_(iConfig.getUntrackedParameter<
std::vector<
int>>(
"newDepth")),
245 hep17_(iConfig.getUntrackedParameter<
bool>(
"hep17")),
253 const double isolationRadius(28.9), innerR(10.0), outerR(30.0);
279 for (
unsigned int k = 0;
k <
oldID_.size(); ++
k) {
289 tok_bs_ = consumes<reco::BeamSpot>(labelBS);
294 tok_alg_ = consumes<BXVector<GlobalAlgBlk>>(algTag);
296 if (modnam.empty()) {
319 <<
"Parameters read from config file \n"
343 for (
unsigned int k = 0;
k <
oldID_.size(); ++
k)
347 for (
int i = 0;
i < 10;
i++)
349 for (
int i = 0;
i < 8; ++
i)
357 unsigned int k1(0), k2(0);
374 <<
" Luminosity " <<
iEvent.luminosityBlock() <<
" Bunch " <<
iEvent.bunchCrossing();
405 respCorrs->
setTopo(theHBHETopology);
415 if (!trkCollection.
isValid()) {
434 if (recVtxs.
isValid() && !(recVtxs->empty())) {
436 for (
unsigned int k = 0;
k < recVtxs->size(); ++
k) {
437 if (!((*recVtxs)[
k].isFake()) && ((*recVtxs)[
k].ndof() > 4)) {
457 if (!barrelRecHitsHandle.
isValid()) {
463 if (!endcapRecHitsHandle.
isValid()) {
469 if (!
hbhe.isValid()) {
477 std::vector<spr::propagatedTrackDirection> trkCaloDirections;
479 std::vector<spr::propagatedTrackID> trkCaloDets;
481 std::vector<math::XYZTLorentzVector> vecL1, vecL3;
506 for (
const auto& decision : finalDecisions) {
507 if (decision.first.find(
l1TrigName_) != std::string::npos) {
514 <<
" from a list of " << finalDecisions.size() <<
" decisions";
524 for (
unsigned int iHLT = 0; iHLT <
triggerResults->size(); iHLT++) {
545 std::array<int, 3> ntksave{{0, 0, 0}};
573 if (!triggerEventHandle.
isValid()) {
580 std::vector<std::string>
modules;
583 for (
unsigned int iHLT = 0; iHLT <
triggerResults->size(); iHLT++) {
588 std::vector<math::XYZTLorentzVector> vecL2;
592 for (
unsigned int ifilter = 0; ifilter <
triggerEvent.sizeFilters(); ++ifilter) {
593 std::vector<int>
Keys;
596 for (
unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
597 if (
label.find(moduleLabels[imodule]) != std::string::npos) {
601 for (
unsigned int ifiltrKey = 0; ifiltrKey <
triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
614 <<
"key " << ifiltrKey <<
" : pt " << TO.
pt() <<
" eta " << TO.
eta() <<
" phi " << TO.
phi()
615 <<
" mass " << TO.
mass() <<
" Id " << TO.
id();
620 <<
"sizes " << vecL1.size() <<
":" << vecL2.size() <<
":" << vecL3.size();
628 for (
unsigned int i = 0;
i < vecL2.size();
i++) {
629 double dr =
dR(vecL1[0], vecL2[
i]);
635 mindRvec1 = vecL2[
i];
639 edm::LogVerbatim(
"HcalIsoTrack") <<
"L2 object closest to L1 " << mindRvec1 <<
" at Dr " << mindR1;
642 if (!vecL1.empty()) {
649 if (!vecL3.empty()) {
693 tree =
fs->
make<TTree>(
"CalibTree",
"CalibTree");
695 tree->Branch(
"t_Run", &
t_Run,
"t_Run/I");
710 tree->Branch(
"t_p", &
t_p,
"t_p/D");
711 tree->Branch(
"t_pt", &
t_pt,
"t_pt/D");
712 tree->Branch(
"t_phi", &
t_phi,
"t_phi/D");
734 t_DetIds =
new std::vector<unsigned int>();
735 t_DetIds1 =
new std::vector<unsigned int>();
736 t_DetIds3 =
new std::vector<unsigned int>();
741 tree->Branch(
"t_DetIds",
"std::vector<unsigned int>", &
t_DetIds);
744 tree->Branch(
"t_DetIds1",
"std::vector<unsigned int>", &
t_DetIds1);
745 tree->Branch(
"t_DetIds3",
"std::vector<unsigned int>", &
t_DetIds3);
749 tree2 =
fs->
make<TTree>(
"EventInfo",
"Event Information");
781 <<
" init flag " <<
flag <<
" change flag " << changed_;
789 for (
unsigned itrig = 0; itrig <
trigNames_.size(); itrig++) {
791 if (triggerindx >=
n) {
793 <<
"the current menu";
811 std::vector<std::string>
trig = {
"HLT_PFJet40",
821 desc.
add<std::vector<std::string>>(
"triggers",
trig);
828 desc.
add<
double>(
"minTrackPt", 1.0);
829 desc.
add<
double>(
"maxDxyPV", 0.02);
830 desc.
add<
double>(
"maxDzPV", 0.02);
831 desc.
add<
double>(
"maxChi2", 5.0);
832 desc.
add<
double>(
"maxDpOverP", 0.1);
833 desc.
add<
int>(
"minOuterHit", 4);
834 desc.
add<
int>(
"minLayerCrossed", 8);
835 desc.
add<
int>(
"maxInMiss", 0);
836 desc.
add<
int>(
"maxOutMiss", 0);
838 desc.
add<
double>(
"minimumTrackP", 20.0);
839 desc.
add<
double>(
"coneRadius", 34.98);
841 desc.
add<
double>(
"coneRadiusMIP", 14.0);
842 desc.
add<
double>(
"coneRadiusMIP2", 18.0);
843 desc.
add<
double>(
"coneRadiusMIP3", 20.0);
844 desc.
add<
double>(
"coneRadiusMIP4", 22.0);
845 desc.
add<
double>(
"coneRadiusMIP5", 24.0);
846 desc.
add<
double>(
"maximumEcalEnergy", 2.0);
848 desc.
add<
double>(
"maxTrackP", 8.0);
849 desc.
add<
double>(
"slopeTrackP", 0.05090504066);
850 desc.
add<
double>(
"isolationEnergyTight", 2.0);
851 desc.
add<
double>(
"isolationEnergyLoose", 10.0);
853 desc.
add<
double>(
"EBHitEnergyThreshold", 0.08);
854 desc.
add<
double>(
"EEHitEnergyThreshold0", 0.30);
855 desc.
add<
double>(
"EEHitEnergyThreshold1", 0.00);
856 desc.
add<
double>(
"EEHitEnergyThreshold2", 0.00);
857 desc.
add<
double>(
"EEHitEnergyThreshold3", 0.00);
858 desc.
add<
double>(
"EEHitEnergyThresholdLow", 0.30);
859 desc.
add<
double>(
"EEHitEnergyThresholdHigh", 0.30);
861 desc.
add<
double>(
"momentumLow", 40.0);
862 desc.
add<
double>(
"momentumHigh", 60.0);
863 desc.
add<
int>(
"prescaleLow", 1);
864 desc.
add<
int>(
"prescaleHigh", 1);
890 std::vector<int>
dummy;
894 descriptions.
add(
"HcalIsoTrkAnalyzer", desc);
898 std::vector<math::XYZTLorentzVector>& vecL3,
900 std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
901 std::vector<spr::propagatedTrackID>& trkCaloDets,
913 int nSave(0),
nLoose(0), nTight(0);
915 std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
916 unsigned int nTracks(0), nselTracks(0);
917 t_nTrk = trkCaloDirections.size();
919 for (trkDetItr = trkCaloDirections.begin(),
nTracks = 0; trkDetItr != trkCaloDirections.end();
921 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
925 << pTrack->
eta() <<
"|" << pTrack->
phi() <<
"|" << pTrack->
p();
928 for (
unsigned int k = 0;
k < vecL3.size(); ++
k) {
929 double dr =
dR(vecL3[
k], v4);
934 t_mindR1 = (!vecL1.empty()) ?
dR(vecL1[0], v4) : 999;
939 if (trkDetItr->okHCAL) {
943 if (
t_p > 40.0 &&
t_p <= 60.0)
950 oneCutParameters.
maxDzPV = 100;
953 bool qltyFlag =
spr::goodTrack(pTrack, leadPV, oneCutParameters,
false);
956 oneCutParameters.
maxDzPV = 100;
968 edm::LogVerbatim(
"HcalIsoTrack") <<
"qltyFlag|okECAL|okHCAL : " << qltyFlag <<
"|" << trkDetItr->okECAL <<
"|"
969 << trkDetItr->okHCAL <<
" eIsolation " << eIsolation;
971 t_qltyFlag = (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL);
976 std::vector<DetId> eIds;
977 std::vector<double> eHit;
981 trkDetItr->pointHCAL,
982 trkDetItr->pointECAL,
984 trkDetItr->directionECAL,
988 for (
unsigned int k = 0;
k < eIds.size(); ++
k) {
998 std::vector<DetId> eIds2;
999 std::vector<double> eHit2;
1001 barrelRecHitsHandle,
1002 endcapRecHitsHandle,
1003 trkDetItr->pointHCAL,
1004 trkDetItr->pointECAL,
1006 trkDetItr->directionECAL,
1010 for (
unsigned int k = 0;
k < eIds2.size(); ++
k) {
1020 std::vector<DetId> eIds3;
1021 std::vector<double> eHit3;
1023 barrelRecHitsHandle,
1024 endcapRecHitsHandle,
1025 trkDetItr->pointHCAL,
1026 trkDetItr->pointECAL,
1028 trkDetItr->directionECAL,
1032 for (
unsigned int k = 0;
k < eIds3.size(); ++
k) {
1042 std::vector<DetId> eIds4;
1043 std::vector<double> eHit4;
1045 barrelRecHitsHandle,
1046 endcapRecHitsHandle,
1047 trkDetItr->pointHCAL,
1048 trkDetItr->pointECAL,
1050 trkDetItr->directionECAL,
1054 for (
unsigned int k = 0;
k < eIds4.size(); ++
k) {
1064 std::vector<DetId> eIds5;
1065 std::vector<double> eHit5;
1067 barrelRecHitsHandle,
1068 endcapRecHitsHandle,
1069 trkDetItr->pointHCAL,
1070 trkDetItr->pointECAL,
1072 trkDetItr->directionECAL,
1076 for (
unsigned int k = 0;
k < eIds5.size(); ++
k) {
1087 const DetId cellE(trkDetItr->detIdECAL);
1089 barrelRecHitsHandle,
1090 endcapRecHitsHandle,
1102 barrelRecHitsHandle,
1103 endcapRecHitsHandle,
1114 if (e11x11P.second && e15x15P.second) {
1115 t_eAnnular = (e15x15P.first - e11x11P.first);
1117 t_eAnnular = -(e15x15P.first - e11x11P.first);
1120 const DetId cellH(trkDetItr->detIdHCAL);
1122 theHBHETopology, cellH,
hbhe, 2, 2,
false,
true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1124 theHBHETopology, cellH,
hbhe, 3, 3,
false,
true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1128 <<
" Annular E (Ecal) " << e11x11P.first <<
":" << e15x15P.first <<
":"
1139 int nRecHits(-999), nRecHits1(-999), nRecHits3(-999);
1140 std::vector<DetId> ids, ids1, ids3;
1141 std::vector<double> edet0, edet1, edet3;
1144 trkDetItr->pointHCAL,
1145 trkDetItr->pointECAL,
1147 trkDetItr->directionHCAL,
1153 for (
unsigned k = 0;
k < ids.size(); ++
k)
1161 trkDetItr->pointHCAL,
1162 trkDetItr->pointECAL,
1164 trkDetItr->directionHCAL,
1170 for (
unsigned k = 0;
k < ids1.size(); ++
k)
1178 trkDetItr->pointHCAL,
1179 trkDetItr->pointECAL,
1181 trkDetItr->directionHCAL,
1187 for (
unsigned k = 0;
k < ids3.size(); ++
k)
1198 << pTrack->
eta() <<
"|" <<
t_phi <<
"|" <<
t_p <<
" Generator Level p "
1203 for (
unsigned int ll = 0; ll <
t_DetIds->size(); ll++) {
1207 for (
unsigned int ll = 0; ll <
t_DetIds1->size(); ll++) {
1211 for (
unsigned int ll = 0; ll <
t_DetIds3->size(); ll++) {
1261 std::array<int, 3>
i3{{nSave,
nLoose, nTight}};
1273 double mindR(999.9);
1278 pmom =
p.momentum().R();
1286 std::vector<double> sumPFNallSMDQH2;
1292 for (
const auto& pf_it : (*
tower)) {
1297 hadder += pf_it.hadEt();
1299 sumPFNallSMDQH2.emplace_back(hadder);
1304 std::sort(sumPFNallSMDQH2.begin(), sumPFNallSMDQH2.end());
1305 if (sumPFNallSMDQH2.size() % 2)
1306 evt_smdq = sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 1) / 2];
1308 evt_smdq = (sumPFNallSMDQH2[sumPFNallSMDQH2.size() / 2] + sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 2) / 2]) / 2.;
1334 for (
unsigned int k = 0;
k <
oldID_.size(); ++
k) {
1344 const std::vector<DetId>& ids,
1345 std::vector<double>& edet,
1347 std::vector<unsigned int>* detIds,
1348 std::vector<double>* hitEnergies) {
1351 for (
unsigned int k = 0;
k < ids.size(); ++
k) {
1358 for (
const auto& en : edet)
1362 edm::LogWarning(
"HcalIsoTrack") <<
"Check inconsistent energies: " << indx <<
" " << eHcal <<
":" << ehcal
1363 <<
" from " << ids.size() <<
" cells";
1367 std::map<HcalDetId, double> hitMap;
1368 for (
unsigned int k = 0;
k < ids.size(); ++
k) {
1370 auto itr = hitMap.find(
id);
1371 if (
itr == hitMap.end()) {
1372 hitMap[
id] = edet[
k];
1374 (
itr->second) += edet[
k];
1377 detIds->reserve(hitMap.size());
1378 hitEnergies->reserve(hitMap.size());
1379 for (
const auto&
hit : hitMap) {
1380 detIds->emplace_back(
hit.first.rawId());
1381 hitEnergies->emplace_back(
hit.second);
1384 detIds->reserve(ids.size());
1385 hitEnergies->reserve(ids.size());
1386 for (
unsigned int k = 0;
k < ids.size(); ++
k) {
1387 detIds->emplace_back(ids[
k].rawId());
1388 hitEnergies->emplace_back(edet[
k]);
1392 edm::LogVerbatim(
"HcalIsoTrack") <<
"Input to storeEnergy with " << ids.size() <<
" cells";
1393 for (
unsigned int k = 0;
k < ids.size(); ++
k)
1395 edm::LogVerbatim(
"HcalIsoTrack") <<
"Output of storeEnergy with " << detIds->size() <<
" cells and Etot " << eHcal;
1396 for (
unsigned int k = 0;
k < detIds->size(); ++
k)