101 void fillTrack(
int,
double,
double,
double,
double);
103 void fillEnergy(
int,
int,
double,
double,
double,
double,
double);
169 : verbosity_(iConfig.getUntrackedParameter<
int>(
"verbosity", 0)),
172 labelMuon_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"labelMuon")),
173 labelGenTrack_(iConfig.getUntrackedParameter<
edm::
InputTag>(
"labelTrack")),
174 theTrackQuality_(iConfig.getUntrackedParameter<
std::
string>(
"trackQuality",
"highPurity")),
175 minTrackP_(iConfig.getUntrackedParameter<double>(
"minTrackP", 1.0)),
176 maxTrackEta_(iConfig.getUntrackedParameter<double>(
"maxTrackEta", 2.5)),
177 tMinE_(iConfig.getUntrackedParameter<double>(
"timeMinCutECAL", -500.)),
178 tMaxE_(iConfig.getUntrackedParameter<double>(
"timeMaxCutECAL", 500.)),
179 tMinH_(iConfig.getUntrackedParameter<double>(
"timeMinCutHCAL", -500.)),
180 tMaxH_(iConfig.getUntrackedParameter<double>(
"timeMaxCutHCAL", 500.)),
181 isItAOD_(iConfig.getUntrackedParameter<
bool>(
"isItAOD",
false)),
182 vetoTrigger_(iConfig.getUntrackedParameter<
bool>(
"vetoTrigger",
false)),
183 doTree_(iConfig.getUntrackedParameter<
bool>(
"doTree",
false)),
184 vetoMuon_(iConfig.getUntrackedParameter<
bool>(
"vetoMuon",
true)),
185 vetoEcal_(iConfig.getUntrackedParameter<
bool>(
"vetoEcal",
false)),
186 cutMuon_(iConfig.getUntrackedParameter<double>(
"cutMuon", 0.1)),
187 cutEcal_(iConfig.getUntrackedParameter<double>(
"cutEcal", 2.0)),
188 cutRatio_(iConfig.getUntrackedParameter<double>(
"cutRatio", 0.90)),
189 puWeights_(iConfig.getUntrackedParameter<
std::
vector<double> >(
"puWeights")),
190 triggerEvent_(
edm::
InputTag(
"hltTriggerSummaryAOD",
"",
"HLT")),
191 theTriggerResultsLabel_(
edm::
InputTag(
"TriggerResults",
"",
"HLT")),
227 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
229 tok_topo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
230 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
232 tok_sevlv_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
250 double pBins[
nPBin_ + 1] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 9.0, 11.0, 15.0, 20.0, 25.0, 30.0, 40.0, 60.0, 100.0};
252 int pvBins[
nPVBin_ + 1] = {1, 2, 3, 5, 100};
265 std::vector<std::string>
trig;
267 std::vector<std::string> newNames = {
"HLT",
"PixelTracks_Multiplicity",
"HLT_Physics_",
"HLT_JetE",
"HLT_ZeroBias"};
270 desc.addUntracked<
int>(
"verbosity", 0);
271 desc.addUntracked<std::vector<std::string> >(
"triggers",
trig);
272 desc.addUntracked<std::vector<std::string> >(
"newNames", newNames);
276 desc.addUntracked<
double>(
"minTrackPt", 1.0);
277 desc.addUntracked<
double>(
"maxDxyPV", 0.02);
278 desc.addUntracked<
double>(
"maxDzPV", 0.02);
279 desc.addUntracked<
double>(
"maxChi2", 5.0);
280 desc.addUntracked<
double>(
"maxDpOverP", 0.1);
281 desc.addUntracked<
int>(
"minOuterHit", 4);
282 desc.addUntracked<
int>(
"minLayerCrossed", 8);
283 desc.addUntracked<
int>(
"maxInMiss", 0);
284 desc.addUntracked<
int>(
"maxOutMiss", 0);
285 desc.addUntracked<
double>(
"minTrackP", 1.0);
286 desc.addUntracked<
double>(
"maxTrackEta", 2.6);
287 desc.addUntracked<
double>(
"timeMinCutECAL", -500.0);
288 desc.addUntracked<
double>(
"timeMaxCutECAL", 500.0);
289 desc.addUntracked<
double>(
"timeMinCutHCAL", -500.0);
290 desc.addUntracked<
double>(
"timeMaxCutHCAL", 500.0);
291 desc.addUntracked<
bool>(
"isItAOD",
false);
292 desc.addUntracked<
bool>(
"vetoTrigger",
false);
293 desc.addUntracked<
bool>(
"doTree",
false);
294 desc.addUntracked<
bool>(
"vetoMuon",
true);
295 desc.addUntracked<
double>(
"cutMuon", 0.1);
296 desc.addUntracked<
bool>(
"vetoEcal",
false);
297 desc.addUntracked<
double>(
"cutEcal", 2.0);
298 desc.addUntracked<
double>(
"cutRatio", 0.9);
299 desc.addUntracked<std::vector<double> >(
"puWeights",
weights);
300 descriptions.
add(
"studyCaloResponse",
desc);
305 int counter0[1000] = {0};
306 int counter1[1000] = {0};
307 int counter2[1000] = {0};
308 int counter3[1000] = {0};
309 int counter4[1000] = {0};
310 int counter5[1000] = {0};
311 int counter6[1000] = {0};
312 int counter7[1000] = {0};
314 edm::LogVerbatim(
"IsoTrack") <<
"Event starts====================================";
315 int RunNo =
iEvent.id().run();
316 int EvtNo =
iEvent.id().event();
318 int Bunch =
iEvent.bunchCrossing();
320 std::vector<int> newAccept(
newNames_.size() + 1, 0);
322 edm::LogVerbatim(
"IsoTrack") <<
"RunNo " << RunNo <<
" EvtNo " << EvtNo <<
" Lumi " <<
Lumi <<
" Bunch " << Bunch;
330 if (!triggerEventHandle.
isValid()) {
344 const std::vector<std::string>& triggerNames_ =
triggerNames.triggerNames();
345 for (
unsigned int iHLT = 0; iHLT <
triggerResults->size(); iHLT++) {
357 if (ipos <= h_HLTAccept->GetNbinsX())
358 h_HLTAccept->GetXaxis()->SetBinLabel(ipos, newtriggerName.c_str());
361 edm::LogVerbatim(
"IsoTrack") <<
"Wrong trigger " << RunNo <<
" Event " << EvtNo <<
" Hlt " << iHLT;
375 if (newtriggerName.find(
trigNames_[
i]) != std::string::npos) {
380 triggerUse = newtriggerName;
389 if (newtriggerName.find(
newNames_[
i]) != std::string::npos) {
398 int iflg(0), indx(1);
400 iflg += (indx * newAccept[
i]);
407 edm::LogVerbatim(
"IsoTrack") <<
"Trigger check gives " <<
ok <<
" with " << triggerUse;
422 int ntrk(0), ngoodPV(0), nPV(-1);
423 int nvtxs = (
int)(recVtxs->size());
424 for (
int ind = 0; ind < nvtxs; ind++) {
425 if (!((*recVtxs)[ind].isFake()) && (*recVtxs)[ind].ndof() > 4)
442 edm::LogVerbatim(
"IsoTrack") <<
"Number of vertices: " << nvtxs <<
" Good " << ngoodPV <<
" Bin " << nPV
449 int npbin =
h_goodPV->FindBin(ngoodPV);
450 if (npbin > 0 && npbin <= (
int)(
puWeights_.size()))
467 double pt1 =
p.momentum().Rho();
468 double p1 =
p.momentum().R();
469 double eta1 =
p.momentum().Eta();
470 double phi1 =
p.momentum().Phi();
480 double phi0 = trk.phi();
483 double dPhi = phi0 - phi2;
500 reco::TrackCollection::const_iterator trkItr;
501 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); ++trkItr, ++ntrk) {
503 double pt1 = pTrack->
pt();
504 double p1 = pTrack->
p();
506 double phi1 = pTrack->
momentum().phi();
513 ++counter0[(
int)(
p1)];
518 std::vector<spr::propagatedTrackID> trkCaloDets;
520 std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
521 for (trkDetItr = trkCaloDets.begin(), ntrk = 0; trkDetItr != trkCaloDets.end(); trkDetItr++, ntrk++) {
522 const reco::Track* pTrack = &(*(trkDetItr->trkItr));
523 double pt1 = pTrack->
pt();
524 double p1 = pTrack->
p();
526 double phi1 = pTrack->
momentum().phi();
529 <<
" okEcal " << trkDetItr->okECAL;
532 bool vetoMuon(
false);
533 double chiGlobal(0),
dr(0);
534 bool goodGlob(
false);
536 if (muonEventHandle.
isValid()) {
537 for (reco::MuonCollection::const_iterator recMuon = muonEventHandle->begin();
538 recMuon != muonEventHandle->end();
540 if (((recMuon->isPFMuon()) && (recMuon->isGlobalMuon() || recMuon->isTrackerMuon())) &&
541 (recMuon->innerTrack()->validFraction() > 0.49) && (recMuon->innerTrack().isNonnull())) {
542 chiGlobal = ((recMuon->globalTrack().isNonnull()) ? recMuon->globalTrack()->normalizedChi2() : 999);
543 goodGlob = (recMuon->isGlobalMuon() && chiGlobal < 3 &&
544 recMuon->combinedQuality().chi2LocalPosition < 12 && recMuon->combinedQuality().trkKink < 20);
559 <<
":" << goodGlob <<
":" <<
dr;
562 double maxNearP31x31 =
572 std::pair<double, bool> e7x7P, e11x11P, e15x15P;
573 const DetId isoCell = trkDetItr->detIdECAL;
617 double maxNearHcalP7x7 =
620 double h3x3(0), h5x5(0), h7x7(0);
621 fillIsolation(0, maxNearP31x31, e11x11P.first, e15x15P.first);
623 edm::LogVerbatim(
"IsoTrack") <<
"Accepted Tracks reaching Ecal maxNearP31x31 " << maxNearP31x31 <<
" e11x11P "
624 << e11x11P.first <<
" e15x15P " << e15x15P.first <<
" okHCAL "
625 << trkDetItr->okHCAL;
628 if (trkDetItr->okHCAL) {
631 const DetId ClosestCell(trkDetItr->detIdHCAL);
676 double eByh = ((e11x11P.second) ? (e11x11P.first /
std::max(h3x3, 0.001)) : 0.0);
679 edm::LogVerbatim(
"IsoTrack") <<
"Tracks Reaching Hcal maxNearHcalP7x7/h5x5/h7x7 " << maxNearHcalP7x7 <<
"/"
680 << h5x5 <<
"/" << h7x7 <<
" eByh " << eByh <<
" notAnElec " << notAnElec;
699 if (maxNearP31x31 < 0 && notAnElec) {
702 if (maxNearHcalP7x7 < 0) {
705 if ((e11x11P.second) && (e15x15P.second) && (e15x15P.first - e11x11P.first) < 2.0) {
708 if (h7x7 - h5x5 < 2.0) {
721 ++counter7[(
int)(
p1)];
726 ++counter6[(
int)(
p1)];
731 ++counter5[(
int)(
p1)];
736 ++counter4[(
int)(
p1)];
741 ++counter3[(
int)(
p1)];
746 ++counter2[(
int)(
p1)];
751 ++counter1[(
int)(
p1)];
757 for (
int i = 0;
i < 1000; ++
i) {
781 h_nHLT =
fs_->
make<TH1I>(
"h_nHLT",
"size of trigger Names", 1000, 0, 1000);
782 h_HLTAccept =
fs_->
make<TH1I>(
"h_HLTAccept",
"HLT Accepts for all runs", 500, 0, 500);
783 for (
int i = 1;
i <= 500; ++
i)
785 h_nHLTvsRN =
fs_->
make<TH2I>(
"h_nHLTvsRN",
"size of trigger Names vs RunNo", 2168, 190949, 193116, 100, 400, 500);
786 h_HLTCorr =
fs_->
make<TH1I>(
"h_HLTCorr",
"Correlation among different paths", 100, 0, 100);
787 h_numberPV =
fs_->
make<TH1I>(
"h_numberPV",
"Number of Primary Vertex", 100, 0, 100);
788 h_goodPV =
fs_->
make<TH1I>(
"h_goodPV",
"Number of good Primary Vertex", 100, 0, 100);
789 h_goodRun =
fs_->
make<TH1I>(
"h_goodRun",
"Number of accepted events for Run", 4000, 190000, 1940000);
790 char hname[60], htit[200];
791 std::string CollectionNames[2] = {
"Reco",
"Propagated"};
792 for (
unsigned int i = 0;
i < 2;
i++) {
793 sprintf(hname,
"h_nTrk_%s", CollectionNames[
i].c_str());
794 sprintf(htit,
"Number of %s tracks", CollectionNames[
i].c_str());
798 "All",
"Quality",
"NoIso",
"okEcal",
"EcalCharIso",
"HcalCharIso",
"EcalNeutIso",
"HcalNeutIso"};
799 std::string particle[4] = {
"Electron",
"Pion",
"Kaon",
"Proton"};
800 for (
unsigned int i = 0;
i <=
nGen_ + 1;
i++) {
802 sprintf(hname,
"h_pt_%s", TrkNames[
i].c_str());
803 sprintf(htit,
"p_{T} of %s tracks", TrkNames[
i].c_str());
805 sprintf(hname,
"h_pt_%s_%d", TrkNames[7].c_str(),
i - 8);
806 sprintf(htit,
"p_{T} of %s tracks (PV=%d:%d)", TrkNames[7].c_str(),
pvBin_[
i - 8],
pvBin_[
i - 7] - 1);
808 sprintf(hname,
"h_pt_%s_%d", TrkNames[0].c_str(),
i -
nGen_);
809 sprintf(htit,
"p_{T} of %s Generator tracks", TrkNames[0].c_str());
811 sprintf(hname,
"h_pt_%s_%s", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
812 sprintf(htit,
"p_{T} of %s tracks (%s)", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
818 sprintf(hname,
"h_p_%s", TrkNames[
i].c_str());
819 sprintf(htit,
"Momentum of %s tracks", TrkNames[
i].c_str());
821 sprintf(hname,
"h_p_%s_%d", TrkNames[7].c_str(),
i - 8);
822 sprintf(htit,
"Momentum of %s tracks (PV=%d:%d)", TrkNames[7].c_str(),
pvBin_[
i - 8],
pvBin_[
i - 7] - 1);
824 sprintf(hname,
"h_p_%s_%d", TrkNames[0].c_str(),
i -
nGen_);
825 sprintf(htit,
"Momentum of %s Generator tracks", TrkNames[0].c_str());
827 sprintf(hname,
"h_p_%s_%s", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
828 sprintf(htit,
"Momentum of %s tracks (%s)", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
830 h_p[
i] =
fs_->
make<TH1D>(hname, htit, 400, 0, 200.0);
834 sprintf(hname,
"h_eta_%s", TrkNames[
i].c_str());
835 sprintf(htit,
"Eta of %s tracks", TrkNames[
i].c_str());
837 sprintf(hname,
"h_eta_%s_%d", TrkNames[7].c_str(),
i - 8);
838 sprintf(htit,
"Eta of %s tracks (PV=%d:%d)", TrkNames[7].c_str(),
pvBin_[
i - 8],
pvBin_[
i - 7] - 1);
840 sprintf(hname,
"h_eta_%s_%d", TrkNames[0].c_str(),
i -
nGen_);
841 sprintf(htit,
"Eta of %s Generator tracks", TrkNames[0].c_str());
843 sprintf(hname,
"h_eta_%s_%s", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
844 sprintf(htit,
"Eta of %s tracks (%s)", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
850 sprintf(hname,
"h_phi_%s", TrkNames[
i].c_str());
851 sprintf(htit,
"Phi of %s tracks", TrkNames[
i].c_str());
853 sprintf(hname,
"h_phi_%s_%d", TrkNames[7].c_str(),
i - 8);
854 sprintf(htit,
"Phi of %s tracks (PV=%d:%d)", TrkNames[7].c_str(),
pvBin_[
i - 8],
pvBin_[
i - 7] - 1);
856 sprintf(hname,
"h_phi_%s_%d", TrkNames[0].c_str(),
i -
nGen_);
857 sprintf(htit,
"Phi of %s Generator tracks", TrkNames[0].c_str());
859 sprintf(hname,
"h_phi_%s_%s", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
860 sprintf(htit,
"Phi of %s tracks (%s)", TrkNames[7].c_str(), particle[
i - 8 -
nPVBin_].c_str());
866 for (
unsigned int i = 0;
i < 2;
i++) {
867 sprintf(hname,
"h_maxNearP_%s", IsolationNames[
i].c_str());
868 sprintf(htit,
"Energy in ChargeIso region for %s", IsolationNames[
i].c_str());
872 sprintf(hname,
"h_ene1_%s", IsolationNames[
i].c_str());
873 sprintf(htit,
"Energy in smaller cone for %s", IsolationNames[
i].c_str());
877 sprintf(hname,
"h_ene2_%s", IsolationNames[
i].c_str());
878 sprintf(htit,
"Energy in bigger cone for %s", IsolationNames[
i].c_str());
882 sprintf(hname,
"h_ediff_%s", IsolationNames[
i].c_str());
883 sprintf(htit,
"Energy in NeutralIso region for %s", IsolationNames[
i].c_str());
888 "E_{7x7}",
"H_{3x3}",
"(E_{7x7}+H_{3x3})",
"E_{11x11}",
"H_{5x5}",
"{E_{11x11}+H_{5x5})"};
890 for (
int ip = 0; ip <
nPBin_; ++ip) {
891 for (
int ie = 0; ie <
nEtaBin_; ++ie) {
892 for (
int j = 0;
j < 6; ++
j) {
893 sprintf(hname,
"h_energy_%d_%d_%d_%d",
i, ip, ie,
j);
896 "%s/p (p=%4.1f:%4.1f; i#eta=%d:%d) for tracks with %s",
897 energyNames[
j].c_str(),
902 TrkNames[
i + 4].c_str());
905 "%s/p (p=%4.1f:%4.1f, i#eta=%d:%d, PV=%d:%d) for tracks with %s",
906 energyNames[
j].c_str(),
913 TrkNames[7].c_str());
916 "%s/p (p=%4.1f:%4.1f, i#eta=%d:%d %s) for tracks with %s",
917 energyNames[
j].c_str(),
923 TrkNames[7].c_str());
932 for (
int i = 0;
i < 8; ++
i) {
933 sprintf(hname,
"counter%d",
i);
934 sprintf(htit,
"Counter with cut %d",
i);
936 sprintf(hname,
"h_pTNew%d",
i);
937 sprintf(htit,
"Track momentum with cut %d",
i);
970 char hname[100], htit[400];
973 sprintf(hname,
"h_HLTAccepts_%i", iRun.
run());
974 sprintf(htit,
"HLT Accepts for Run No %i", iRun.
run());
975 TH1I* hnew =
fs_->
make<TH1I>(hname, htit, 500, 0, 500);
976 for (
int i = 1;
i <= 500; ++
i)
977 hnew->GetXaxis()->SetBinLabel(
i,
" ");
1026 int flag,
int ieta,
double p,
double enEcal1,
double enHcal1,
double enEcal2,
double enHcal2) {
1040 if (ip >= 0 && ie >= 0 && enEcal1 > 0.02 && enHcal1 > 0.1) {
1052 int length =
str.length();
1053 for (
int i = 0;
i < length - 2;
i++) {
1054 if (
str[
i] ==
'_' &&
str[
i + 1] ==
'v' && isdigit(
str.at(
i + 2))) {
1056 truncated_str =
str.substr(0,
z);
1059 return (truncated_str);
1067 reco::GenParticleCollection::const_iterator
p;
1068 double mindR(999.9);
1073 double dEta = pTrack->
eta() -
p->momentum().Eta();
1074 double phi1 = pTrack->
phi();
1075 double phi2 =
p->momentum().Phi();
1080 double dPhi = phi1 - phi2;