Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Run2-hcx148 Correct the access of Geometry information for Hcal #20349

Merged
merged 4 commits into from Sep 11, 2017
Merged

Run2-hcx148 Correct the access of Geometry information for Hcal #20349

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5d973c2
Correct the access of Geometry information for Hcal
Sep 1, 2017
5feb420
After CodeChecker Reports
Sep 1, 2017
4ee6363
Take care of code checker comments
Sep 2, 2017
06cb1b0
Take care of code checker comments
Sep 2, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
14 changes: 2 additions & 12 deletions JetMETCorrections/MinBias/BuildFile.xml
View file
Open in desktop
@@ -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
Open in desktop
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
Open in desktop
@@ -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);
}
}
}

}
}