Skip to content

Commit

Permalink
Merge pull request #20349 from bsunanda/Run2-hcx148
Browse files Browse the repository at this point in the history 
Run2-hcx148 Correct the access of Geometry information for Hcal
  • Loading branch information
@cmsbuild
cmsbuild committed Sep 11, 2017
2 parents 6caa3ff + 06cb1b0 commit 69441fc
Show file tree
Hide file tree
Showing 3 changed files with 128 additions and 175 deletions.
14 changes: 2 additions & 12 deletions JetMETCorrections/MinBias/BuildFile.xml
View file
@@ -1,15 +1,5 @@
<use name="boost"/>
<use name="DataFormats/CaloTowers"/>
<use name="DataFormats/HcalRecHit"/>
<use name="DataFormats/EcalDetId"/>
<use name="DataFormats/DetId"/>
<use name="DataFormats/RecoCandidate"/>
<use name="Geometry/CaloTopology"/>
<use name="Geometry/CaloGeometry"/>
<use name="FWCore/Framework"/>
<use name="FWCore/ParameterSet"/>
<use name="FWCore/PluginManager"/>
<use name="Geometry/Records"/>
<use name="clhep"/>
<use name="Geometry/HcalTowerAlgo"/>
<use name="CommonTools/UtilAlgos"/>
<use name="root"/>
<flags EDM_PLUGIN="1"/>
34 changes: 12 additions & 22 deletions JetMETCorrections/MinBias/interface/MinBias.h
View file
Expand Up @@ -13,17 +13,10 @@
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
#include "DataFormats/DetId/interface/DetId.h"

#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"

#include "TFile.h"
#include "TTree.h"


Expand All @@ -42,11 +35,11 @@ namespace cms
class MinBias : public edm::EDAnalyzer {
public:
explicit MinBias(const edm::ParameterSet&);
~MinBias();

virtual void analyze(const edm::Event&, const edm::EventSetup&);
virtual void beginJob() ;
virtual void endJob() ;
void analyze(const edm::Event&, const edm::EventSetup&) override;
void beginJob() override;
void beginRun(edm::Run const&, edm::EventSetup const&) override;
void endJob() override;

private:
// ----------member data ---------------------------
Expand All @@ -57,24 +50,21 @@ class MinBias : public edm::EDAnalyzer {
edm::EDGetTokenT<HORecHitCollection> hoToken_;
edm::EDGetTokenT<HFRecHitCollection> hfToken_;
// stuff for histogramms
// output file name with histograms
std::string fOutputFileName ;
bool allowMissingInputs_;
//
TFile* hOutputFile ;
// TH1D* hCalo1[8000], *hCalo2;
TTree * myTree;
TTree * myTree_;
//
int mydet, mysubd, depth, iphi, ieta, ievent;
int mydet, mysubd, depth, iphi, ieta;
float phi,eta;
float mom1,mom2,mom3,mom4,occup;
const CaloGeometry* geo;
const CaloGeometry* geo_;
// counters
std::map<DetId,double> theFillDetMap0;
std::map<DetId,double> theFillDetMap1;
std::map<DetId,double> theFillDetMap2;
std::map<DetId,double> theFillDetMap3;
std::map<DetId,double> theFillDetMap4;
std::map<DetId,double> theFillDetMap0_;
std::map<DetId,double> theFillDetMap1_;
std::map<DetId,double> theFillDetMap2_;
std::map<DetId,double> theFillDetMap3_;
std::map<DetId,double> theFillDetMap4_;

};
}
255 changes: 114 additions & 141 deletions JetMETCorrections/MinBias/src/MinBias.cc
View file
@@ -1,100 +1,101 @@
// user include files
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "JetMETCorrections/MinBias/interface/MinBias.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"

using namespace std;
namespace cms
{
MinBias::MinBias(const edm::ParameterSet& iConfig)
{
namespace cms {
MinBias::MinBias(const edm::ParameterSet& iConfig) : geo_(nullptr) {
// get names of modules, producing object collections
hbheLabel_= iConfig.getParameter<std::string>("hbheInput");
hoLabel_= iConfig.getParameter<std::string>("hoInput");
hfLabel_= iConfig.getParameter<std::string>("hfInput");
hoLabel_ = iConfig.getParameter<std::string>("hoInput");
hfLabel_ = iConfig.getParameter<std::string>("hfInput");
hbheToken_= mayConsume<HBHERecHitCollection>(edm::InputTag(hbheLabel_));
hoToken_=mayConsume<HORecHitCollection>(edm::InputTag(hoLabel_));
hfToken_=mayConsume<HFRecHitCollection>(edm::InputTag(hfLabel_));
hoToken_ = mayConsume<HORecHitCollection>(edm::InputTag(hoLabel_));
hfToken_ = mayConsume<HFRecHitCollection>(edm::InputTag(hfLabel_));
allowMissingInputs_=iConfig.getUntrackedParameter<bool>("AllowMissingInputs",false);
// get name of output file with histogramms
fOutputFileName = iConfig.getUntrackedParameter<std::string>("HistOutFile");

}


MinBias::~MinBias()
{

// do anything here that needs to be done at desctruction time
// (e.g. close files, deallocate resources etc.)
void MinBias::beginJob() {
edm::Service<TFileService> fs;
myTree_ = fs->make<TTree>("RecJet","RecJet Tree");
myTree_->Branch("mydet", &mydet, "mydet/I");
myTree_->Branch("mysubd", &mysubd, "mysubd/I");
myTree_->Branch("depth", &depth, "depth/I");
myTree_->Branch("ieta", &ieta, "ieta/I");
myTree_->Branch("iphi", &iphi, "iphi/I");
myTree_->Branch("eta", &eta, "eta/F");
myTree_->Branch("phi", &phi, "phi/F");
myTree_->Branch("mom1", &mom1, "mom1/F");
myTree_->Branch("mom2", &mom2, "mom2/F");
myTree_->Branch("mom3", &mom3, "mom3/F");
myTree_->Branch("mom4", &mom4, "mom4/F");

}

void MinBias::beginJob()
{
hOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
myTree = new TTree("RecJet","RecJet Tree");
myTree->Branch("mydet", &mydet, "mydet/I");
myTree->Branch("mysubd", &mysubd, "mysubd/I");
myTree->Branch("depth", &depth, "depth/I");
myTree->Branch("ieta", &ieta, "ieta/I");
myTree->Branch("iphi", &iphi, "iphi/I");
myTree->Branch("eta", &eta, "eta/F");
myTree->Branch("phi", &phi, "phi/F");
myTree->Branch("mom1", &mom1, "mom1/F");
myTree->Branch("mom2", &mom2, "mom2/F");
myTree->Branch("mom3", &mom3, "mom3/F");
myTree->Branch("mom4", &mom4, "mom4/F");

ievent = 0;

void MinBias::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
edm::ESHandle<CaloGeometry> pG;
iSetup.get<CaloGeometryRecord>().get(pG);
geo_ = pG.product();
std::vector<DetId> did = geo_->getValidDetIds();

for (auto const& id : did) {
if( (id).det() == DetId::Hcal ) {
theFillDetMap0_[id] = 0.;
theFillDetMap1_[id] = 0.;
theFillDetMap2_[id] = 0.;
theFillDetMap3_[id] = 0.;
theFillDetMap4_[id] = 0.;
}
}
}

void MinBias::endJob()
{

const std::vector<DetId>& did = geo->getSubdetectorGeometry( DetId::Hcal, 1 )->getValidDetIds() ;
int i=0;
for(std::vector<DetId>::const_iterator id=did.begin(); id != did.end(); id++)
{
// if( (*id).det() == DetId::Hcal ) {
GlobalPoint pos = geo->getPosition(*id);
mydet = ((*id).rawId()>>28)&0xF;
mysubd = ((*id).rawId()>>25)&0x7;
depth = HcalDetId(*id).depth();
ieta = HcalDetId(*id).ieta();
iphi = HcalDetId(*id).iphi();
phi = pos.phi();
eta = pos.eta();
if ( theFillDetMap0[*id] > 0. )
{
mom1 = theFillDetMap1[*id]/theFillDetMap0[*id];
mom2 = theFillDetMap2[*id]/theFillDetMap0[*id]-(mom1*mom1);
mom3 = theFillDetMap3[*id]/theFillDetMap0[*id]-3.*mom1*theFillDetMap2[*id]/theFillDetMap0[*id]+
2.*pow(mom2,3);
mom4 = (theFillDetMap4[*id]-4.*mom1*theFillDetMap3[*id]+6.*pow(mom1,2)*theFillDetMap2[*id])/theFillDetMap0[*id]-3.*pow(mom1,4);

} else
{
mom1 = 0.; mom2 = 0.; mom3 = 0.; mom4 = 0.;
}
std::cout<<" Detector "<<(*id).rawId()<<" mydet "<<mydet<<" "<<mysubd<<" "<<depth<<" "<<
HcalDetId(*id).subdet()<<" "<<ieta<<" "<<iphi<<" "<<pos.eta()<<" "<<pos.phi()<<std::endl;
std::cout<<" Energy "<<mom1<<" "<<mom2<<std::endl;
myTree->Fill();
i++;
// }
}
std::cout<<" The number of CaloDet records "<<did.size()<<std::endl;
std::cout<<" The number of Hcal records "<<i<<std::endl;


std::cout << "===== Start writing user histograms =====" << std::endl;
hOutputFile->SetCompressionLevel(2);
hOutputFile->cd();
myTree->Write();
hOutputFile->Close() ;
std::cout << "===== End writing user histograms =======" << std::endl;
void MinBias::endJob() {
const HcalGeometry* hgeo = (const HcalGeometry*)(geo_->getSubdetectorGeometry(DetId::Hcal,1));
const std::vector<DetId>& did = hgeo->getValidDetIds() ;
int i=0;
for (const auto& id : did) {
// if( id.det() == DetId::Hcal ) {
GlobalPoint pos = hgeo->getPosition(id);
mydet = (int)(id.det());
mysubd = (id.subdetId());
depth = HcalDetId(id).depth();
ieta = HcalDetId(id).ieta();
iphi = HcalDetId(id).iphi();
phi = pos.phi();
eta = pos.eta();
if ( theFillDetMap0_[id] > 0. ) {
mom1 = theFillDetMap1_[id]/theFillDetMap0_[id];
mom2 = theFillDetMap2_[id]/theFillDetMap0_[id]-(mom1*mom1);
mom3 = theFillDetMap3_[id]/theFillDetMap0_[id]-3.*mom1*theFillDetMap2_[id]/theFillDetMap0_[id]+
2.*pow(mom2,3);
mom4 = (theFillDetMap4_[id]-4.*mom1*theFillDetMap3_[id]+6.*pow(mom1,2)*theFillDetMap2_[id])/theFillDetMap0_[id]-3.*pow(mom1,4);

} else {
mom1 = 0.; mom2 = 0.; mom3 = 0.; mom4 = 0.;
}
edm::LogWarning("MinBias")<<" Detector "<< id.rawId()<<" mydet "<<mydet
<<" "<<mysubd<<" "<<depth<<" "<<ieta<<" "
<<iphi<<" "<<pos.eta()<<" "<<pos.phi();
edm::LogWarning("MinBias")<<" Energy "<<mom1<<" "<<mom2<<std::endl;
myTree_->Fill();
i++;
// }
}
edm::LogWarning("MinBias")<<" The number of CaloDet records "<<did.size();
edm::LogWarning("MinBias")<<" The number of Hcal records "<<i<<std::endl;

/*
std::cout << "===== Start writing user histograms =====" << std::endl;
hOutputFile->SetCompressionLevel(2);
hOutputFile->cd();
myTree->Write();
hOutputFile->Close() ;
std::cout << "===== End writing user histograms =======" << std::endl;
*/
}


Expand All @@ -103,30 +104,7 @@ void MinBias::endJob()
//

// ------------ method called to produce the data ------------
void
MinBias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{

using namespace edm;
if(ievent == 0 ){
edm::ESHandle<CaloGeometry> pG;
iSetup.get<CaloGeometryRecord>().get(pG);
geo = pG.product();
std::vector<DetId> did = geo->getValidDetIds();

for(std::vector<DetId>::const_iterator id=did.begin(); id != did.end(); id++)
{
if( (*id).det() == DetId::Hcal ) {
theFillDetMap0[*id] = 0.;
theFillDetMap1[*id] = 0.;
theFillDetMap2[*id] = 0.;
theFillDetMap3[*id] = 0.;
theFillDetMap4[*id] = 0.;
}
}
}


void MinBias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup){

if (!hbheLabel_.empty()) {
edm::Handle<HBHERecHitCollection> hbhe;
Expand All @@ -137,20 +115,19 @@ MinBias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
*hbhe; // will throw the proper exception
}
} else {
for(HBHERecHitCollection::const_iterator hbheItr = (*hbhe).begin();
hbheItr != (*hbhe).end(); ++hbheItr)
{
DetId id = (hbheItr)->detid();
if( (*hbheItr).energy() > 0. ) std::cout<<" Energy = "<<(*hbheItr).energy()<<std::endl;
theFillDetMap0[id] = theFillDetMap0[id]+ 1.;
theFillDetMap1[id] = theFillDetMap1[id]+(*hbheItr).energy();
theFillDetMap2[id] = theFillDetMap2[id]+pow((*hbheItr).energy(),2);
theFillDetMap3[id] = theFillDetMap3[id]+pow((*hbheItr).energy(),3);
theFillDetMap4[id] = theFillDetMap4[id]+pow((*hbheItr).energy(),4);
}
for (auto const& hbheItr : (HBHERecHitCollection)(*hbhe)) {
DetId id = (hbheItr).detid();
if (hbheItr.energy() > 0. )
edm::LogWarning("MinBias")<<" Energy = "<<hbheItr.energy();
theFillDetMap0_[id] += 1.;
theFillDetMap1_[id] += hbheItr.energy();
theFillDetMap2_[id] += pow(hbheItr.energy(),2);
theFillDetMap3_[id] += pow(hbheItr.energy(),3);
theFillDetMap4_[id] += pow(hbheItr.energy(),4);
}
}
}

if (!hoLabel_.empty()) {
edm::Handle<HORecHitCollection> ho;
iEvent.getByToken(hoToken_,ho);
Expand All @@ -159,20 +136,18 @@ MinBias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
if (!allowMissingInputs_) {
*ho; // will throw the proper exception
}
} else {
for(HORecHitCollection::const_iterator hoItr = (*ho).begin();
hoItr != (*ho).end(); ++hoItr)
{
DetId id = (hoItr)->detid();
theFillDetMap0[id] = theFillDetMap0[id]+ 1.;
theFillDetMap1[id] = theFillDetMap1[id]+(*hoItr).energy();
theFillDetMap2[id] = theFillDetMap2[id]+pow((*hoItr).energy(),2);
theFillDetMap3[id] = theFillDetMap3[id]+pow((*hoItr).energy(),3);
theFillDetMap4[id] = theFillDetMap4[id]+pow((*hoItr).energy(),4);
}
} else {
for (auto const& hoItr : (HORecHitCollection)(*ho)) {
DetId id = hoItr.detid();
theFillDetMap0_[id] += 1.;
theFillDetMap1_[id] += hoItr.energy();
theFillDetMap2_[id] += pow(hoItr.energy(),2);
theFillDetMap3_[id] += pow(hoItr.energy(),3);
theFillDetMap4_[id] += pow(hoItr.energy(),4);
}
}
}

if (!hfLabel_.empty()) {
edm::Handle<HFRecHitCollection> hf;
iEvent.getByToken(hfToken_,hf);
Expand All @@ -181,19 +156,17 @@ MinBias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
if (!allowMissingInputs_) {
*hf; // will throw the proper exception
}
} else {
for(HFRecHitCollection::const_iterator hfItr = (*hf).begin();
hfItr != (*hf).end(); ++hfItr)
{
DetId id = (hfItr)->detid();
theFillDetMap0[id] = theFillDetMap0[id]+ 1.;
theFillDetMap1[id] = theFillDetMap1[id]+(*hfItr).energy();
theFillDetMap2[id] = theFillDetMap2[id]+pow((*hfItr).energy(),2);
theFillDetMap3[id] = theFillDetMap3[id]+pow((*hfItr).energy(),3);
theFillDetMap4[id] = theFillDetMap4[id]+pow((*hfItr).energy(),4);
}
} else {
for (auto const hfItr : (HFRecHitCollection)(*hf)) {
DetId id = hfItr.detid();
theFillDetMap0_[id] += 1.;
theFillDetMap1_[id] += hfItr.energy();
theFillDetMap2_[id] += pow(hfItr.energy(),2);
theFillDetMap3_[id] += pow(hfItr.energy(),3);
theFillDetMap4_[id] += pow(hfItr.energy(),4);
}
}
}

}
}

0 comments on commit 69441fc

Please sign in to comment.