Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #11045 from kkiesel/fastsimECALresponseCorrection76X
Fastsim ECAL response correction (2)
- Loading branch information
Showing
12 changed files
with
511 additions
and
7 deletions.
There are no files selected for viewing
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
45 changes: 45 additions & 0 deletions
45
FastSimulation/Calorimetry/interface/KKCorrectionFactors.h
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,45 @@ | ||
// -*- C++ -*- | ||
// | ||
// Package: FastSimulation/Calorimetry | ||
// Class: KKCorrectionFactors | ||
// | ||
/**\class KKCorrectionFactorsr | ||
Description: Returns scale from an TH3F histogram in a root file | ||
*/ | ||
// | ||
// Original Author: Maximilian Knut Kiesel | ||
// Created: Fri, 07 Aug 2015 | ||
// | ||
// | ||
|
||
|
||
#ifndef CALORESPONSE_H | ||
#define CALORESPONSE_H | ||
|
||
#include "FWCore/Utilities/interface/Exception.h" | ||
#include "FWCore/ParameterSet/interface/ParameterSet.h" | ||
|
||
#include <TH3F.h> | ||
#include <TFile.h> | ||
#include <TROOT.h> | ||
|
||
|
||
class KKCorrectionFactors { | ||
|
||
public: | ||
/* Constructor: pset must contain two strings: "fileName" is the name of the | ||
* file in which the TH3F named "histogramName" is saved. | ||
*/ | ||
KKCorrectionFactors( const edm::ParameterSet& pset ); | ||
~KKCorrectionFactors() { delete h3_; }; | ||
|
||
float getScale( float genEnergy, float genEta, float simEnergy ) const; | ||
|
||
private: | ||
// histogram which contains the scales | ||
TH3F* h3_; | ||
bool interpolate3D_; | ||
|
||
}; | ||
|
||
#endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,12 @@ | ||
import FWCore.ParameterSet.Config as cms | ||
|
||
ECALScaleBlock = cms.PSet( | ||
ECALResponseScaling = cms.PSet( | ||
fileName = cms.untracked.string("FastSimulation/Calorimetry/data/scaleECALFastsim.root"), | ||
histogramName = cms.untracked.string("scaleVsEVsEta"), | ||
interpolate3D = cms.untracked.bool(False) | ||
) | ||
) | ||
|
||
|
||
|
39 changes: 39 additions & 0 deletions
39
FastSimulation/Calorimetry/python/TreeWriterForEcalCorrection_cfg.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,39 @@ | ||
import FWCore.ParameterSet.Config as cms | ||
process = cms.Process("ANA") | ||
|
||
process.load("FWCore.MessageService.MessageLogger_cfi") | ||
process.MessageLogger.cerr.FwkReport.reportEvery = 10000 | ||
process.MessageLogger.cerr.default.limit = 100000000 | ||
|
||
|
||
process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) | ||
|
||
from python.fastFileList import filelist | ||
#from python.fullFileList import filelist | ||
|
||
process.source = cms.Source("PoolSource", | ||
fileNames = cms.untracked.vstring( | ||
#'/store/user/kiesel/SingleElectronGun_E80_fast/SingleElectronGun_E80_fast/a88c9ccbd07595a7d33ae9d4d6a917a0/out_2_1_X2L.root', | ||
#'/store/user/kiesel/SingleElectronGun_E80_full/SingleElectronGun_E80_full/1d73d2566f2f721e3a7146e7729d743b/out_48_1_kkW.root' | ||
filelist | ||
) | ||
) | ||
|
||
outFileName = "out_tree.root" | ||
|
||
# guess the output name from the 1. input name | ||
fn0 = process.source.fileNames[0] | ||
if "full" in fn0 or "Full" in fn0: | ||
outFileName = "full_tree.root" | ||
if "fast" in fn0 or "Fast" in fn0: | ||
outFileName = "fast_tree.root" | ||
|
||
process.TFileService = cms.Service("TFileService", | ||
fileName = cms.string( outFileName ) | ||
) | ||
|
||
process.load("Configuration.StandardSequences.Analysis_cff") | ||
|
||
process.treeWriterForEcalCorrection = cms.EDAnalyzer('TreeWriterForEcalCorrection') | ||
process.p = cms.Path( process.treeWriterForEcalCorrection ) | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,71 @@ | ||
#include "FastSimulation/Calorimetry/interface/KKCorrectionFactors.h" | ||
|
||
KKCorrectionFactors::KKCorrectionFactors( const edm::ParameterSet& pset ) | ||
{ | ||
|
||
interpolate3D_ = pset.exists("interpolate3D") && pset.getUntrackedParameter<bool>("interpolate3D"); | ||
|
||
// Get the filename and histogram name | ||
std::string fileName = pset.getUntrackedParameter<std::string>("fileName"); | ||
std::string histogramName = pset.getUntrackedParameter<std::string>("histogramName"); | ||
|
||
// Read histo | ||
edm::FileInPath myDataFile( fileName ); | ||
TFile * myFile = TFile::Open( myDataFile.fullPath().c_str() ); | ||
|
||
gROOT->cd(); // create histogram in memory | ||
auto obj = myFile->Get( histogramName.c_str() ); | ||
|
||
// Check if histogram exists in file | ||
if(!obj) { | ||
throw cms::Exception( "FileReadError" ) | ||
<< "Histogram \"" << histogramName | ||
<< "\" not found in file \"" << fileName | ||
<< "\".\n"; | ||
} | ||
h3_ = new TH3F( *((TH3F*)obj) ); | ||
|
||
delete myFile; | ||
|
||
} | ||
|
||
|
||
float KKCorrectionFactors::getScale( float genE, float genEta, float simE ) const | ||
{ | ||
|
||
float r = simE / genE; | ||
float scale = 1.; | ||
|
||
if( interpolate3D_ | ||
// TH3::Interpolate can only interpolate inside the bondaries of the histogram | ||
&& genE > h3_->GetXaxis()->GetXmin() | ||
&& genE < h3_->GetXaxis()->GetXmax() | ||
&& genEta > h3_->GetYaxis()->GetXmin() | ||
&& genEta < h3_->GetYaxis()->GetXmax() | ||
&& r < h3_->GetZaxis()->GetXmax() | ||
&& r > h3_->GetZaxis()->GetXmax() ) { | ||
|
||
scale = h3_->Interpolate( genE, genEta, r ); | ||
|
||
} else { // interpolation in r is mandatory | ||
|
||
int binE = h3_->GetXaxis()->FindFixBin( genE ); | ||
int binEta = h3_->GetYaxis()->FindFixBin( genEta ); | ||
|
||
// find the two bins which are closest to the actual value | ||
auto binWidthR = h3_->GetZaxis()->GetBinWidth(0); | ||
int binRup = h3_->GetZaxis()->FindFixBin( r + binWidthR/2. ); | ||
int binRdn = h3_->GetZaxis()->FindFixBin( r - binWidthR/2. ); | ||
|
||
auto scaleUp = h3_->GetBinContent( binE, binEta, binRup ); | ||
auto scaleDn = h3_->GetBinContent( binE, binEta, binRdn ); | ||
|
||
// make a linear extrapolation between neighbour bins if they are not zero | ||
auto Rdn = h3_->GetZaxis()->GetBinCenter( binRdn ); | ||
scale = scaleDn + (scaleUp-scaleDn) * ( r - Rdn ) / binWidthR; | ||
|
||
} | ||
return scale; | ||
|
||
} | ||
|
112 changes: 112 additions & 0 deletions
112
FastSimulation/Calorimetry/src/TreeWriterForEcalCorrection.cc
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,112 @@ | ||
// system include files | ||
#include <memory> | ||
|
||
// user include files | ||
#include "FWCore/Framework/interface/Frameworkfwd.h" | ||
#include "FWCore/Framework/interface/EDAnalyzer.h" | ||
|
||
#include "FWCore/Framework/interface/Event.h" | ||
#include "FWCore/Framework/interface/MakerMacros.h" | ||
|
||
#include "FWCore/ParameterSet/interface/ParameterSet.h" | ||
|
||
#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" | ||
#include "DataFormats/HepMCCandidate/interface/GenParticle.h" | ||
|
||
#include "FWCore/Framework/interface/ESHandle.h" | ||
|
||
#include "FWCore/ServiceRegistry/interface/Service.h" | ||
#include "CommonTools/UtilAlgos/interface/TFileService.h" | ||
|
||
#include "TTree.h" | ||
|
||
|
||
class TreeWriterForEcalCorrection : public edm::EDAnalyzer { | ||
public: | ||
explicit TreeWriterForEcalCorrection(const edm::ParameterSet&); | ||
~TreeWriterForEcalCorrection(){}; | ||
|
||
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); | ||
|
||
|
||
private: | ||
virtual void analyze(const edm::Event&, const edm::EventSetup&) override; | ||
|
||
edm::Service<TFileService> file; | ||
TTree* tree; | ||
float tree_e, tree_eta, tree_response; | ||
}; | ||
|
||
TreeWriterForEcalCorrection::TreeWriterForEcalCorrection(const edm::ParameterSet& iConfig) | ||
{ | ||
tree = file->make<TTree>("responseTree", "same info as 3dhisto"); | ||
tree->Branch( "e", &tree_e, "e/F"); | ||
tree->Branch( "eta", &tree_eta, "eta/F"); | ||
tree->Branch( "r", &tree_response, "r/F"); | ||
} | ||
|
||
void | ||
TreeWriterForEcalCorrection::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) | ||
{ | ||
// get generated particles | ||
edm::Handle<reco::GenParticleCollection> GenParticles; | ||
iEvent.getByLabel("genParticles", "", GenParticles); | ||
|
||
// As this module is intended for single particle guns, there should be | ||
// exactly one generated particle. | ||
if( GenParticles->size() != 1 ) { | ||
throw cms::Exception("MismatchedInputFiles") | ||
<< "Intended for particle guns only\n"; | ||
} | ||
|
||
// I assume here that the tracker simulation is disabled and no vertex | ||
// smearing is done, so that the generated particles is exactly the same | ||
// particle as the particle hitting the entrace of the ECAL. | ||
reco::GenParticle gen = GenParticles->at(0); | ||
float genE = gen.energy(); | ||
float genEta = gen.eta(); | ||
// genEta is positive for my photongun sample per definition. | ||
// If not, you could take the absolute value here. | ||
|
||
// get sim hits | ||
edm::Handle<edm::PCaloHitContainer> SimHitsEB; | ||
edm::Handle<edm::PCaloHitContainer> SimHitsEE; | ||
edm::Handle<edm::PCaloHitContainer> SimHitsES; | ||
|
||
// Finds out automatically, if this is fullsim or fastsim | ||
bool isFastSim = iEvent.getByLabel( "famosSimHits", "EcalHitsEB", SimHitsEB ); | ||
if( isFastSim ) { | ||
iEvent.getByLabel( "famosSimHits", "EcalHitsEE", SimHitsEE ); | ||
iEvent.getByLabel( "famosSimHits", "EcalHitsES", SimHitsES ); | ||
} else { | ||
iEvent.getByLabel( "g4SimHits", "EcalHitsEB", SimHitsEB ); | ||
iEvent.getByLabel( "g4SimHits", "EcalHitsEE", SimHitsEE ); | ||
iEvent.getByLabel( "g4SimHits", "EcalHitsES", SimHitsES ); | ||
} | ||
|
||
// merge them into one single vector | ||
auto SimHits = *SimHitsEB; | ||
SimHits.insert( SimHits.end(), SimHitsEE->begin(), SimHitsEE->end() ); | ||
SimHits.insert( SimHits.end(), SimHitsES->begin(), SimHitsES->end() ); | ||
|
||
// As we only had one generated particle (and hopefully no pileup), | ||
// the total energy is due to the generated particle only | ||
float energyTotal = 0; | ||
for( auto const& Hit : SimHits ) { | ||
energyTotal += Hit.energy(); | ||
} | ||
|
||
tree_e = genE; | ||
tree_eta = genEta; | ||
tree_response = energyTotal/genE; | ||
tree->Fill(); | ||
|
||
} | ||
|
||
void | ||
TreeWriterForEcalCorrection::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { | ||
edm::ParameterSetDescription desc; | ||
descriptions.add( "ecalScaleFactorCalculator", desc ); | ||
} | ||
|
||
DEFINE_FWK_MODULE(TreeWriterForEcalCorrection); |
Oops, something went wrong.