Fastsim ECAL response correction (2)

FastSimulation/Calorimetry/interface/CalorimetryManager.h

@@ -14,6 +14,8 @@
 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
 #include "FastSimulation/CaloHitMakers/interface/PreshowerHitMaker.h"
 
+#include "FastSimulation/Calorimetry/interface/KKCorrectionFactors.h"
+
 // For the uint32_t
 //#include <boost/cstdint.hpp>
 #include <map>
@@ -186,5 +188,8 @@ class CalorimetryManager{
   bool useShowerLibrary;
   bool useCorrectionSL;
   FastHFShowerLibrary *theHFShowerLibrary;
+
+  std::unique_ptr<KKCorrectionFactors> ecalCorrection;
+
 };
 #endif

FastSimulation/Calorimetry/interface/KKCorrectionFactors.h

@@ -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

FastSimulation/Calorimetry/python/Calorimetry_cff.py

@@ -7,8 +7,12 @@
 from FastSimulation.Calorimetry.HcalResponse_cfi import *
 from FastSimulation.Calorimetry.HSParameters_cfi import *
 #from FastSimulation.Configuration.CommonInputs_cff import *
+
+from FastSimulation.Calorimetry.ECALResponse_cfi import *
+
 FamosCalorimetryBlock = cms.PSet(
     Calorimetry = cms.PSet(
+        #ECALScaleBlock, # comment out to disable scaling
         HSParameterBlock,
         HCALResponseBlock,
         ECAL = cms.PSet(

FastSimulation/Calorimetry/python/ECALResponse_cfi.py

@@ -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)
+    )
+)
+
+
+

FastSimulation/Calorimetry/python/TreeWriterForEcalCorrection_cfg.py

@@ -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 )
+

FastSimulation/Calorimetry/src/CalorimetryManager.cc

@@ -53,6 +53,7 @@
 //#include "DataFormats/EcalDetId/interface/EcalDetId.h"
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/Utilities/interface/EDMException.h"
 
 //ROOT headers
 #include "TROOT.h"
@@ -139,8 +140,11 @@ CalorimetryManager::CalorimetryManager(FSimEvent * aSimEvent,
        fastMuHCAL.getParameter<bool>("EnergyLoss") || 
        fastMuHCAL.getParameter<bool>("MultipleScattering") )
     theMuonHcalEffects = new MaterialEffects(fastMuHCAL);
-
-
+
+  if( fastCalo.exists("ECALResponseScaling") ) {
+    ecalCorrection = std::unique_ptr<KKCorrectionFactors>( new KKCorrectionFactors( fastCalo.getParameter<edm::ParameterSet>("ECALResponseScaling") ) );
+  }
+
 }
 
 void CalorimetryManager::clean()
@@ -458,10 +462,19 @@ void CalorimetryManager::EMShowerSimulation(const FSimTrack& myTrack,
   theShower.setHcal(&myHcalHitMaker);
   theShower.compute();
   //myHistos->fill("h502", myPart->eta(),myGrid.totalX0());
-
+
+  // calculate the total simulated energy for this particle
+  float simE = 0;
+  for( const auto& mapIterator : myGrid.getHits() ) {
+    simE += mapIterator.second;
+  }
+
+  auto scale = ecalCorrection ? ecalCorrection->getScale( myTrack.ecalEntrance().e(),
+     std::abs( myTrack.ecalEntrance().eta() ), simE ) : 1.;
+
   // Save the hits !
-  updateECAL(myGrid.getHits(),onEcal,myTrack.id());
-  
+  updateECAL( myGrid.getHits(), onEcal,myTrack.id(), scale );
+
   // Now fill the HCAL hits
   updateHCAL(myHcalHitMaker.getHits(),myTrack.id());
 
@@ -474,8 +487,6 @@ void CalorimetryManager::EMShowerSimulation(const FSimTrack& myTrack,
 
 }
 
-
-
 void CalorimetryManager::reconstructHCAL(const FSimTrack& myTrack,
                                          RandomEngineAndDistribution const* random)
 {

FastSimulation/Calorimetry/src/KKCorrectionFactors.cc

@@ -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;
+
+}
+

FastSimulation/Calorimetry/src/TreeWriterForEcalCorrection.cc

@@ -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);