Skip to content

Commit

Permalink
Update isotrack code needed for validation of simulation
Browse files Browse the repository at this point in the history
  • Loading branch information
Sunanda committed Sep 2, 2016
1 parent b92f0ad commit 8c128f5
Show file tree
Hide file tree
Showing 4 changed files with 176 additions and 122 deletions.
113 changes: 66 additions & 47 deletions Calibration/IsolatedParticles/plugins/StudyHLT.cc
View file
Expand Up @@ -57,6 +57,7 @@ StudyHLT::StudyHLT(const edm::ParameterSet& iConfig) : nRun(0) {
tMaxH_ = iConfig.getUntrackedParameter<double>("TimeMaxCutHCAL", 500.);
isItAOD_ = iConfig.getUntrackedParameter<bool>("IsItAOD", false);
doTree_ = iConfig.getUntrackedParameter<bool>("DoTree", false);
puWeights_ = iConfig.getUntrackedParameter<std::vector<double> >("PUWeights");
triggerEvent_ = edm::InputTag("hltTriggerSummaryAOD","","HLT");
theTriggerResultsLabel_ = edm::InputTag("TriggerResults","","HLT");

Expand All @@ -75,6 +76,7 @@ StudyHLT::StudyHLT(const edm::ParameterSet& iConfig) : nRun(0) {
tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
}
tok_ew_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));

edm::LogInfo("IsoTrack") << "Verbosity " << verbosity_ << " with "
<< trigNames_.size() << " triggers:";
Expand Down Expand Up @@ -126,16 +128,15 @@ void StudyHLT::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup)
int Lumi = iEvent.luminosityBlock();
int Bunch = iEvent.bunchCrossing();

edm::Handle<LumiDetails> Lumid;
iEvent.getLuminosityBlock().getByToken(tok_lumi,Lumid);

std::string newNames[5]={"HLT","PixelTracks_Multiplicity","HLT_Physics_","HLT_JetE","HLT_ZeroBias"};
int newAccept[5];
for (int i=0; i<5; ++i) newAccept[i] = 0;
float mybxlumi=-1;
if (Lumid.isValid())
mybxlumi=Lumid->lumiValue(LumiDetails::kOCC1,iEvent.bunchCrossing())*6.37;

/*
edm::Handle<LumiDetails> Lumid;
iEvent.getLuminosityBlock().getByToken(tok_lumi,Lumid);
if (Lumid.isValid()) mybxlumi=Lumid->lumiValue(LumiDetails::kOCC1,iEvent.bunchCrossing())*6.37;
*/
if (verbosity_ > 0)
edm::LogInfo("IsoTrack") << "RunNo " << RunNo << " EvtNo " << EvtNo
<< " Lumi " << Lumi << " Bunch " << Bunch
Expand All @@ -146,9 +147,7 @@ void StudyHLT::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup)
iEvent.getByToken(tok_trigEvt,triggerEventHandle);

bool ok(false);
if (trigNames_.size() < 1) {
ok = true;
} else if (!triggerEventHandle.isValid()) {
if (!triggerEventHandle.isValid()) {
edm::LogWarning("IsoTrack") << "Error! Can't get the product "
<< triggerEvent_.label();
} else {
Expand Down Expand Up @@ -192,22 +191,26 @@ void StudyHLT::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup)
h_HLTAccepts[nRun]->Fill(iHLT+1);
h_HLTAccept->Fill(ipos);
}
for (unsigned int i=0; i<trigNames_.size(); ++i) {
if (newtriggerName.find(trigNames_[i].c_str())!=std::string::npos) {
if (verbosity_%10 > 0)
edm::LogInfo("IsoTrack") << newtriggerName;
if (hlt > 0) {
ok = true;
tr_TrigName.push_back(newtriggerName);
if (trigNames_.size() < 1) {
ok = true;
} else {
for (unsigned int i=0; i<trigNames_.size(); ++i) {
if (newtriggerName.find(trigNames_[i].c_str())!=std::string::npos) {
if (verbosity_%10 > 0)
edm::LogInfo("IsoTrack") << newtriggerName;
if (hlt > 0) {
ok = true;
tr_TrigName.push_back(newtriggerName);
}
}
}
}
for (int i=0; i<5; ++i) {
if (newtriggerName.find(newNames[i].c_str())!=std::string::npos) {
if (verbosity_%10 > 0)
edm::LogInfo("IsoTrack") << "[" << i << "] " << newNames[i]
<< " : " << newtriggerName;
if (hlt > 0) newAccept[i] = 1;
for (int i=0; i<5; ++i) {
if (newtriggerName.find(newNames[i].c_str())!=std::string::npos) {
if (verbosity_%10 > 0)
edm::LogInfo("IsoTrack") << "[" << i << "] " << newNames[i]
<< " : " << newtriggerName;
if (hlt > 0) newAccept[i] = 1;
}
}
}
}
Expand Down Expand Up @@ -246,22 +249,37 @@ void StudyHLT::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup)

edm::Handle<reco::VertexCollection> recVtxs;
iEvent.getByToken(tok_recVtx_,recVtxs);
int ntrk(0), ngoodPV(0), nPV(-1);
for (unsigned int ind=0; ind<recVtxs->size(); ind++) {
int ntrk(0), ngoodPV(0), nPV(-1);
int nvtxs = (int)(recVtxs->size());
for (int ind=0; ind<nvtxs; ind++) {
if (!((*recVtxs)[ind].isFake()) && (*recVtxs)[ind].ndof() > 4) ngoodPV++;
}
for (int i=0; i<nPVBin; ++i) {
if (ngoodPV >= pvBin[i] && ngoodPV < pvBin[i+1]) {
nPV = i; break;
}
}

tr_eventWeight = 1.0;
edm::Handle<GenEventInfoProduct> genEventInfo;
iEvent.getByToken(tok_ew_, genEventInfo);
if (genEventInfo.isValid()) tr_eventWeight = genEventInfo->weight();

if ((verbosity_/10)%10 > 0)
edm::LogInfo("IsoTrack") << "Number of vertices: " << recVtxs->size()
<< " Good " << ngoodPV << " Bin " << nPV;
h_numberPV->Fill((int)(recVtxs->size()));
h_goodPV->Fill(ngoodPV);
edm::LogInfo("IsoTrack") << "Number of vertices: " << nvtxs
<< " Good " << ngoodPV << " Bin " << nPV
<< " Event weight " << tr_eventWeight;
h_numberPV->Fill(nvtxs,tr_eventWeight);
h_goodPV->Fill(ngoodPV,tr_eventWeight);
tr_goodPV = ngoodPV;

if (puWeights_.size() > 0) {
int npbin = h_goodPV->FindBin(ngoodPV);
if (npbin > 0 && npbin <= (int)(puWeights_.size()))
tr_eventWeight *= puWeights_[npbin-1];
else
tr_eventWeight = 0;
}
edm::Handle<reco::TrackCollection> trkCollection;
iEvent.getByToken(tok_genTrack_, trkCollection);
reco::TrackCollection::const_iterator trkItr;
Expand All @@ -275,7 +293,7 @@ void StudyHLT::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup)
fillTrack(0, pt1,p1,eta1,phi1);
if (quality) fillTrack(1, pt1,p1,eta1,phi1);
}
h_ntrk[0]->Fill(ntrk);
h_ntrk[0]->Fill(ntrk,tr_eventWeight);

std::vector<spr::propagatedTrackID> trkCaloDets;
spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDets, ((verbosity_/100)%10 > 0));
Expand Down Expand Up @@ -371,7 +389,7 @@ void StudyHLT::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup)
}
}
}
h_ntrk[1]->Fill(ntrk);
h_ntrk[1]->Fill(ntrk,tr_eventWeight);
if (tr_TrkPt.size() > 0 && doTree_) tree_->Fill();
}
firstEvent_ = false;
Expand Down Expand Up @@ -485,8 +503,9 @@ void StudyHLT::beginJob() {
// Now the tree
if (doTree_) {
tree_ = fs_->make<TTree>("testTree", "new HLT Tree");
tree_->Branch("tr_goodRun", &tr_goodRun, "tr_goodRun/I");
tree_->Branch("tr_goodPV", &tr_goodPV, "tr_goodPV/I");
tree_->Branch("tr_goodRun", &tr_goodRun, "tr_goodRun/I");
tree_->Branch("tr_goodPV", &tr_goodPV, "tr_goodPV/I");
tree_->Branch("tr_eventWeight", &tr_eventWeight, "tr_eventWeight/D");
tree_->Branch("tr_tr_TrigName", &tr_TrigName);
tree_->Branch("tr_TrkPt", &tr_TrkPt);
tree_->Branch("tr_TrkP", &tr_TrkP);
Expand Down Expand Up @@ -544,17 +563,17 @@ void StudyHLT::clear() {
}

void StudyHLT::fillTrack(int i, double pt, double p, double eta, double phi){
h_pt[i]->Fill(pt);
h_p[i]->Fill(p);
h_eta[i]->Fill(eta);
h_phi[i]->Fill(phi);
h_pt[i]->Fill(pt,tr_eventWeight);
h_p[i]->Fill(p,tr_eventWeight);
h_eta[i]->Fill(eta,tr_eventWeight);
h_phi[i]->Fill(phi,tr_eventWeight);
}

void StudyHLT::fillIsolation(int i, double emaxnearP, double eneutIso1, double eneutIso2){
h_maxNearP[i]->Fill(emaxnearP);
h_ene1[i]->Fill(eneutIso1);
h_ene2[i]->Fill(eneutIso2);
h_ediff[i]->Fill(eneutIso2-eneutIso1);
h_maxNearP[i]->Fill(emaxnearP,tr_eventWeight);
h_ene1[i]->Fill(eneutIso1,tr_eventWeight);
h_ene2[i]->Fill(eneutIso2,tr_eventWeight);
h_ediff[i]->Fill(eneutIso2-eneutIso1,tr_eventWeight);
}

void StudyHLT::fillEnergy(int flag, int ieta, double p, double enEcal1,
Expand All @@ -567,12 +586,12 @@ void StudyHLT::fillEnergy(int flag, int ieta, double p, double enEcal1,
if (ieta >= etaBin[i] && ieta < etaBin[i+1]) { ie = i; break; }
}
if (ip >= 0 && ie >= 0 && enEcal1 > 0.02 && enHcal1 > 0.1) {
h_energy[flag][ip][ie][0]->Fill(enEcal1/p);
h_energy[flag][ip][ie][1]->Fill(enHcal1/p);
h_energy[flag][ip][ie][2]->Fill((enEcal1+enHcal1)/p);
h_energy[flag][ip][ie][3]->Fill(enEcal2/p);
h_energy[flag][ip][ie][4]->Fill(enHcal2/p);
h_energy[flag][ip][ie][5]->Fill((enEcal2+enHcal2)/p);
h_energy[flag][ip][ie][0]->Fill(enEcal1/p,tr_eventWeight);
h_energy[flag][ip][ie][1]->Fill(enHcal1/p,tr_eventWeight);
h_energy[flag][ip][ie][2]->Fill((enEcal1+enHcal1)/p,tr_eventWeight);
h_energy[flag][ip][ie][3]->Fill(enEcal2/p,tr_eventWeight);
h_energy[flag][ip][ie][4]->Fill(enHcal2/p,tr_eventWeight);
h_energy[flag][ip][ie][5]->Fill((enEcal2+enHcal2)/p,tr_eventWeight);
}
}

Expand Down
4 changes: 4 additions & 0 deletions Calibration/IsolatedParticles/plugins/StudyHLT.h
View file
Expand Up @@ -30,6 +30,7 @@
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

class StudyHLT : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::SharedResources> {

Expand Down Expand Up @@ -61,6 +62,7 @@ class StudyHLT : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::Share
spr::trackSelectionParameters selectionParameters_;
std::vector<std::string> trigNames_, HLTNames_;
std::string theTrackQuality_;
std::vector<double> puWeights_;
double minTrackP_, maxTrackEta_;
double tMinE_, tMaxE_, tMinH_, tMaxH_;
bool isItAOD_, changed_, firstEvent_, doTree_;
Expand All @@ -75,6 +77,7 @@ class StudyHLT : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::Share
edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;

TH1I *h_nHLT, *h_HLTAccept, *h_HLTCorr, *h_numberPV;
TH1I *h_goodPV, *h_goodRun;
Expand All @@ -89,6 +92,7 @@ class StudyHLT : public edm::one::EDAnalyzer<edm::one::WatchRuns,edm::one::Share
int nRun, etaBin[nEtaBin+1], pvBin[nPVBin+1];
double pBin[nPBin+1];
int tr_goodPV, tr_goodRun;
double tr_eventWeight;
std::vector<std::string> tr_TrigName;
std::vector<double> tr_TrkPt, tr_TrkP, tr_TrkEta, tr_TrkPhi;
std::vector<double> tr_MaxNearP31X31, tr_MaxNearHcalP7x7;
Expand Down
3 changes: 2 additions & 1 deletion Calibration/IsolatedParticles/python/studyHLT_cfi.py
View file
Expand Up @@ -20,5 +20,6 @@
TimeMinCutHCAL = cms.untracked.double(-500.0),
TimeMaxCutHCAL = cms.untracked.double(500.0),
IsItAOD = cms.untracked.bool(False),
DoTree = cms.untracked.bool(True),
DoTree = cms.untracked.bool(False),
PUWeights = cms.untracked.vdouble([]),
)

0 comments on commit 8c128f5

Please sign in to comment.