Skip to content

Commit

Permalink
Merge pull request #16350 from bsunanda/Run2-alca70
Browse files Browse the repository at this point in the history 
bsunanda:Run2-alca70 Generator level information usage for closure tests
  • Loading branch information
@cmsbuild
cmsbuild committed Oct 26, 2016
2 parents ba926f4 + 2c7f20e commit f35b4d6
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 74 deletions.
54 changes: 29 additions & 25 deletions Calibration/HcalCalibAlgos/test/CalibMonitor.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
// Usage:
// .L CalibMonitor.C+g
// CalibMonitor c1(fname, dirname, dupFileName, outFileName, prefix,
// flag, numb, dataMC);
// flag, numb, dataMC, useGen);
// c1.Loop();
// c1.SavePlot(histFileName,mode);
//
Expand Down Expand Up @@ -34,6 +34,8 @@
// for loose/tight selection)
// numb (int) = number of eta bins (42 for -21:21)
// dataMC (bool) = true/false for data/MC
// useGen (bool) = false/true to use generator level momentum
// or reconstruction level momentum
//
// histFileName (std::string)= name of the file containing saved histograms
// mode (bool) = true/false if the hitogram file to be opened
Expand Down Expand Up @@ -139,7 +141,7 @@ public :
CalibMonitor(std::string fname, std::string dirname,
std::string dupFileName, std::string outTxtFileName,
std::string prefix="", int flag=0, int numb=42,
bool datMC=true);
bool datMC=true, bool useGen=false);
virtual ~CalibMonitor();
virtual Int_t Cut(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry);
Expand All @@ -157,7 +159,7 @@ private:
static const unsigned int npbin=5, kp50=2;
std::string fname_, dirnm_, prefix_, outTxtFileName_;
int flag_, numb_;
bool dataMC_, plotStandard_, flexibleSelect_;
bool dataMC_, plotStandard_, flexibleSelect_, useGen_;
double log16by24_;
std::vector<Long64_t> entries_;
std::vector<double> etas_, ps_, dl1_;
Expand All @@ -171,12 +173,12 @@ private:

CalibMonitor::CalibMonitor(std::string fname, std::string dirnm,
std::string dupFileName, std::string outTxtFileName,
std::string prefix, int flag, int numb,
bool datMC) : fname_(fname), dirnm_(dirnm),
prefix_(prefix),
outTxtFileName_(outTxtFileName),
flag_(flag), numb_(numb),
dataMC_(datMC) {
std::string prefix, int flag, int numb, bool dataMC,
bool useGen) : fname_(fname), dirnm_(dirnm),
prefix_(prefix),
outTxtFileName_(outTxtFileName),
flag_(flag), numb_(numb),
dataMC_(dataMC), useGen_(useGen) {
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree

Expand Down Expand Up @@ -532,8 +534,9 @@ void CalibMonitor::Loop() {

// if (Cut(ientry) < 0) continue;
int kp(-1), jp(-1);
double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
for (unsigned int k=1; k<ps_.size(); ++k ) {
if (t_p >= ps_[k-1] && t_p < ps_[k]) {
if (pmom >= ps_[k-1] && pmom < ps_[k]) {
kp = k - 1; break;
}
}
Expand All @@ -560,29 +563,29 @@ void CalibMonitor::Loop() {
<< kv1 << ":" << kd << ":" << kd1 << ":" << jp
<< std::endl;
if (plotStandard_) {
h_p[0]->Fill(t_p,t_EventWeight);
h_p[0]->Fill(pmom,t_EventWeight);
h_eta[0]->Fill(t_ieta,t_EventWeight);
if (kp >= 0) h_eta0[kp]->Fill(t_ieta,t_EventWeight);
}
double cut = (t_p > 20) ? 2.0 : 0.0;
double rcut= (t_p > 20) ? 0.25: 0.1;
double cut = (pmom > 20) ? 2.0 : 0.0;
double rcut= (pmom > 20) ? 0.25: 0.1;

// Some Standard plots for control
if (plotStandard_) {
if (t_qltyFlag) {
h_p[1]->Fill(t_p,t_EventWeight);
h_p[1]->Fill(pmom,t_EventWeight);
h_eta[1]->Fill(t_ieta,t_EventWeight);
if (kp >= 0) h_eta1[kp]->Fill(t_ieta,t_EventWeight);
if (t_selectTk) {
h_p[2]->Fill(t_p,t_EventWeight);
h_p[2]->Fill(pmom,t_EventWeight);
h_eta[2]->Fill(t_ieta,t_EventWeight);
if (kp >= 0) h_eta2[kp]->Fill(t_ieta,t_EventWeight);
if (t_hmaxNearP < cut) {
h_p[3]->Fill(t_p,t_EventWeight);
h_p[3]->Fill(pmom,t_EventWeight);
h_eta[3]->Fill(t_ieta,t_EventWeight);
if (kp >= 0) h_eta3[kp]->Fill(t_ieta,t_EventWeight);
if (t_eMipDR < 1.0) {
h_p[4]->Fill(t_p,t_EventWeight);
h_p[4]->Fill(pmom,t_EventWeight);
h_eta[4]->Fill(t_ieta,t_EventWeight);
if (kp >= 0) {
h_eta4[kp]->Fill(t_ieta,t_EventWeight);
Expand All @@ -598,9 +601,9 @@ void CalibMonitor::Loop() {
// Selection of good track and energy measured in Hcal
double rat(1.0), eHcal(t_eHcal);
bool goodTk = GoodTrack(eHcal, cut, debug);
if (t_p > 0) rat = (eHcal/(t_p-t_eMipDR));
if (pmom > 0) rat = (eHcal/(pmom-t_eMipDR));
if (debug)
std::cout << "Entry " << jentry << " p|eHcal|ratio " << t_p << "|"
std::cout << "Entry " << jentry << " p|eHcal|ratio " << pmom << "|"
<< t_eHcal << "|" << eHcal << "|" << rat << "|" << kp << "|"
<< kv << "|" << jp << " Cuts " << t_qltyFlag << "|"
<< t_selectTk << "|" << (t_hmaxNearP < cut) << "|"
Expand Down Expand Up @@ -632,7 +635,7 @@ void CalibMonitor::Loop() {
h_etaF[kp][0]->Fill(rat,t_EventWeight);
}
}
if (t_p > 20.0) {
if (pmom > 20.0) {
if (plotStandard_) {
h_etaX[kp1][kv]->Fill(eta,rat,t_EventWeight);
h_etaX[kp1][kv1]->Fill(eta,rat,t_EventWeight);
Expand All @@ -646,12 +649,12 @@ void CalibMonitor::Loop() {
}
}
}
if (t_p > 20.0) {
if (pmom > 20.0) {
kount++;
if (((flag_/100)%10) != 0) {
good++;
fileout << good << " " << jentry << " " << t_Run << " "
<< t_Event << " " << t_ieta << " " << t_p << std::endl;
<< t_Event << " " << t_ieta << " " << pmom << std::endl;
}
}
}
Expand All @@ -671,14 +674,15 @@ void CalibMonitor::Loop() {
bool CalibMonitor::GoodTrack(double& eHcal, double &cut, bool debug) {

bool select(true);
double pmom = (useGen_ && (t_gentrackP>0)) ? t_gentrackP : t_p;
if (debug) std::cout << "GoodTrack input " << eHcal << ":" << cut;
if (flexibleSelect_) {
double eta = (t_ieta > 0) ? t_ieta : -t_ieta;
cut = 2.0*exp(eta*log16by24_);
select = ((t_qltyFlag) && (t_hmaxNearP < cut) && (t_eMipDR < 1.0));
double fac = (t_p > 0) ?
(1.0 - 0.375 * (eHcal/t_p) *
((t_eHcalDelta/t_p) - 0.45*(t_eHcalDelta/t_p)*(t_eHcalDelta/t_p))) : 1.;
double fac = (pmom > 0) ?
(1.0 - 0.375 * (eHcal/pmom) *
((t_eHcalDelta/pmom) - 0.45*(t_eHcalDelta/pmom)*(t_eHcalDelta/pmom))) : 1.;
eHcal *= fac;
} else {
select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < cut) &&
Expand Down
72 changes: 56 additions & 16 deletions Calibration/HcalCalibAlgos/test/HcalIsoTrkAnalyzer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,6 @@
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

//Triggers
#include "DataFormats/Common/interface/TriggerResults.h"
Expand All @@ -34,6 +33,11 @@
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/Math/interface/deltaR.h"

//Generator information
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

Expand Down Expand Up @@ -89,8 +93,11 @@ class HcalIsoTrkAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::
const CaloGeometry* geo,
edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
edm::Handle<HBHERecHitCollection>& hbhe);
edm::Handle<HBHERecHitCollection>& hbhe,
edm::Handle<reco::GenParticleCollection>& genParticles);
double dR(math::XYZTLorentzVector&, math::XYZTLorentzVector&);
double trackP(const reco::Track* ,
const edm::Handle<reco::GenParticleCollection>&);

unsigned int nRun_;
edm::Service<TFileService> fs;
Expand All @@ -107,21 +114,23 @@ class HcalIsoTrkAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::
std::string theTrackQuality_, processName_, labelHBHE_;
std::string l1Filter_, l2Filter_, l3Filter_;
bool ignoreTrigger_, useRaw_;
edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;
edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;

TTree *tree, *tree2;
int t_Run, t_Event, t_ieta, t_goodPV, t_DataType;
double t_EventWeight, t_l1pt, t_l1eta, t_l1phi;
double t_l3pt, t_l3eta, t_l3phi, t_p, t_mindR1;
double t_mindR2, t_eMipDR, t_eHcal, t_eHcalDelta, t_hmaxNearP;
double t_mindR2, t_eMipDR, t_eHcal, t_eHcalDelta;
double t_hmaxNearP, t_gentrackP;
bool t_selectTk,t_qltyFlag,t_qltyMissFlag,t_qltyPVFlag;
std::vector<unsigned int> *t_DetIds;
std::vector<double> *t_HitEnergies, pbin;
Expand Down Expand Up @@ -185,6 +194,7 @@ HcalIsoTrkAnalyzer::HcalIsoTrkAnalyzer(const edm::ParameterSet& iConfig) :
tok_bs_ = consumes<reco::BeamSpot>(labelBS);
tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
tok_ew_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
tok_parts_ = consumes<reco::GenParticleCollection>(edm::InputTag("genParticles"));
if (modnam == "") {
tok_recVtx_ = consumes<reco::VertexCollection>(labelRecVtx_);
tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit",labelEB_));
Expand Down Expand Up @@ -262,6 +272,10 @@ void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const
edm::ESHandle<CaloGeometry> pG;
iSetup.get<CaloGeometryRecord>().get(pG);
const CaloGeometry* geo = pG.product();

//=== genParticle information
edm::Handle<reco::GenParticleCollection> genParticles;
iEvent.getByToken(tok_parts_, genParticles);

bool okC(true);
//Get track collection
Expand Down Expand Up @@ -345,7 +359,8 @@ void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const
t_l1pt = t_l1eta = t_l1phi = 0;
t_l3pt = t_l3eta = t_l3phi = 0;
t_TracksSaved = fillTree(vecL1, vecL3, leadPV, trkCaloDirections, geo,
barrelRecHitsHandle, endcapRecHitsHandle, hbhe);
barrelRecHitsHandle, endcapRecHitsHandle, hbhe,
genParticles);
} else {
trigger::TriggerEvent triggerEvent;
edm::Handle<trigger::TriggerEvent> triggerEventHandle;
Expand Down Expand Up @@ -460,7 +475,8 @@ void HcalIsoTrkAnalyzer::analyze(edm::Event const& iEvent, edm::EventSetup const
// Now fill in the tree for each selected track
if (!done) {
ntksave = fillTree(vecL1, vecL3, leadPV, trkCaloDirections, geo,
barrelRecHitsHandle,endcapRecHitsHandle, hbhe);
barrelRecHitsHandle,endcapRecHitsHandle, hbhe,
genParticles);
done = true;
}
t_TracksSaved += ntksave;
Expand Down Expand Up @@ -499,6 +515,7 @@ void HcalIsoTrkAnalyzer::beginJob() {
tree->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
tree->Branch("t_qltyMissFlag",&t_qltyMissFlag,"t_qltyMissFlag/O");
tree->Branch("t_qltyPVFlag", &t_qltyPVFlag, "t_qltyPVFlag/O)");
tree->Branch("t_gentrackP", &t_gentrackP, "t_gentrackP/D");

t_DetIds = new std::vector<unsigned int>();
t_HitEnergies = new std::vector<double>();
Expand Down Expand Up @@ -576,7 +593,8 @@ int HcalIsoTrkAnalyzer::fillTree(std::vector< math::XYZTLorentzVector>& vecL1,
const CaloGeometry* geo,
edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
edm::Handle<HBHERecHitCollection>& hbhe) {
edm::Handle<HBHERecHitCollection>& hbhe,
edm::Handle<reco::GenParticleCollection>& genParticles) {

int nsave = 0;
//Loop over tracks
Expand Down Expand Up @@ -674,11 +692,15 @@ int HcalIsoTrkAnalyzer::fillTree(std::vector< math::XYZTLorentzVector>& vecL1,
<< t_eHcal << ":" << ehcal
<< " from " << t_DetIds->size()
<< " cells\n";

t_gentrackP = trackP(pTrack, genParticles);

#ifdef DebugLog
edm::LogInfo("HcalIsoTrack") << "This track : " << nTracks
<< " (pt|eta|phi|p) :" << pTrack->pt()
<< "|" << pTrack->eta() << "|"
<< pTrack->phi() << "|" << t_p;
<< pTrack->phi() << "|" << t_p
<< " Generator Level p " << t_gentrackP;
edm::LogInfo("HcalIsoTrack") << "e_MIP " << t_eMipDR
<< " Chg Isolation " << t_hmaxNearP
<< " eHcal" << t_eHcal << ":" << t_eHcalDelta
Expand Down Expand Up @@ -712,5 +734,23 @@ double HcalIsoTrkAnalyzer::dR(math::XYZTLorentzVector& vec1, math::XYZTLorentzVe
return reco::deltaR(vec1.eta(),vec1.phi(),vec2.eta(),vec2.phi());
}

double HcalIsoTrkAnalyzer::trackP(const reco::Track* pTrack,
const edm::Handle<reco::GenParticleCollection>& genParticles) {

double pmom = -1.0;
if (genParticles.isValid()) {
reco::GenParticleCollection::const_iterator p;
double mindR(999.9);
for (p=genParticles->begin(); p!=genParticles->end(); ++p) {
double dR = reco::deltaR(pTrack->eta(), pTrack->phi(),
p->momentum().Eta(), p->momentum().Phi());
if (dR < mindR) {
mindR = dR; pmom = p->momentum().R();
}
}
}
return pmom;
}

//define this as a plug-in
DEFINE_FWK_MODULE(HcalIsoTrkAnalyzer);

0 comments on commit f35b4d6

Please sign in to comment.