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

HGCal Reworked rechit scheme #13835

Merged
merged 5 commits into from Apr 6, 2016
Merged

HGCal Reworked rechit scheme #13835

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b8efe3b
working on tuning the digitizer
lgray Mar 22, 2016
4f73f0f
cross check of all conversion factors
lgray Mar 24, 2016
1963e4a
Uncalibrated RecHits are now in MIPs
lgray Mar 24, 2016
ce7b1a8
switch default to hex geom, use proper unit conversions
lgray Mar 24, 2016
4f5269c
address code review
lgray Apr 4, 2016
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
1 change: 0 additions & 1 deletion RecoLocalCalo/HGCalRecAlgos/interface/HGCalRecHitSimpleAlgo.h
View file
Open in desktop
Expand Up @@ -56,6 +56,5 @@ class HGCalRecHitSimpleAlgo : public HGCalRecHitAbsAlgo {
private:
float adcToGeVConstant_;
bool adcToGeVConstantIsSet_;

};
#endif
46 changes: 31 additions & 15 deletions RecoLocalCalo/HGCalRecAlgos/interface/HGCalUncalibRecHitRecWeightsAlgo.h
View file
Open in desktop
Expand Up @@ -17,20 +17,29 @@
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"

#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"

template<class C> class HGCalUncalibRecHitRecWeightsAlgo
{
public:
// destructor
virtual ~HGCalUncalibRecHitRecWeightsAlgo<C>() { };

virtual void set_isSiFESim(const bool isSiFE) { isSiFESim_ = isSiFE; }
void set_isSiFESim(const bool isSiFE) { isSiFESim_ = isSiFE; }

void set_ADCLSB(const double adclsb) { adcLSB_ = adclsb; }
void set_TDCLSB(const double tdclsb) { tdcLSB_ = tdclsb; }

virtual void set_ADCLSB(const double adclsb) { adcLSB_ = adclsb; }
virtual void set_TDCLSB(const double tdclsb) { tdcLSB_ = tdclsb; }
void set_toaLSBToNS(const double lsb2ns) { toaLSBToNS_ = lsb2ns; }

virtual void set_toaLSBToNS(const double lsb2ns) { toaLSBToNS_ = lsb2ns; }
void set_tdcOnsetfC(const double tdcOnset) { tdcOnsetfC_ = tdcOnset; }

virtual void set_tdcOnsetfC(const double tdcOnset) { tdcOnsetfC_ = tdcOnset; }
void set_fCPerMIP(const std::vector<double>& fCPerMIP) { fCPerMIP_ = fCPerMIP; }

void setGeometry(const HGCalGeometry* geom) {
if ( geom ) ddd_ = &(geom->topology().dddConstants());
else ddd_ = nullptr;
}

/// Compute parameters
virtual HGCUncalibratedRecHit makeRecHit( const C& dataFrame ) {
Expand All @@ -51,30 +60,37 @@ template<class C> class HGCalUncalibRecHitRecWeightsAlgo
// need to increment by one so TDC doesn't overlap with ADC last bin
amplitude_ = ( std::floor(tdcOnsetfC_/adcLSB_) + 1.0 )* adcLSB_ + double(sample.data()) * tdcLSB_;
jitter_ = double(sample.toa()) * toaLSBToNS_;
LogDebug("HGCUncalibratedRecHit") << "TDC+: set the energy to: " << amplitude_ << ' ' << sample.data()
LogDebug("HGCUncalibratedRecHit") << "TDC+: set the charge to: " << amplitude_ << ' ' << sample.data()
<< ' ' << tdcLSB_ << std::endl
<< "TDC+: set the jitter to: " << jitter_ << ' '
<< "TDC+: set the ToA to: " << jitter_ << ' '
<< sample.toa() << ' ' << toaLSBToNS_ << ' '
<< " flag=" << flag << std::endl;
}
else {
} else {
amplitude_ = double(sample.data()) * adcLSB_;
LogDebug("HGCUncalibratedRecHit") << "ADC+: set the energy to: " << amplitude_ << ' ' << sample.data()
LogDebug("HGCUncalibratedRecHit") << "ADC+: set the charge to: " << amplitude_ << ' ' << sample.data()
<< ' ' << adcLSB_ << ' ' << std::endl;
}
}
else {
} else {
amplitude_ = double(sample.data()) * adcLSB_;
LogDebug("HGCUncalibratedRecHit") << "ADC+: set the energy to: " << amplitude_ << ' ' << sample.data()
LogDebug("HGCUncalibratedRecHit") << "ADC+: set the charge to: " << amplitude_ << ' ' << sample.data()
<< ' ' << adcLSB_ << ' ' << std::endl;
}

int thickness = 1;
if( ddd_ != nullptr ) {
HGCalDetId hid(dataFrame.id());
thickness = ddd_->waferTypeL(hid.wafer());
}
amplitude_ = amplitude_/fCPerMIP_[thickness-1];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please check if fCPerMIP_[thickness-1] is zero before dividing by it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will check this in the constructor so it is guaranteed non-zero here.


LogDebug("HGCUncalibratedRecHit") << "Final uncalibrated amplitude : " << amplitude_ << std::endl;
return HGCUncalibratedRecHit( dataFrame.id(), amplitude_, pedestal_, jitter_, chi2_, flag);
}

private:
double adcLSB_, tdcLSB_, fCToMIP_, toaLSBToNS_, tdcOnsetfC_;
bool isSiFESim_;
double adcLSB_, tdcLSB_, toaLSBToNS_, tdcOnsetfC_;
bool isSiFESim_;
std::vector<double> fCPerMIP_;
const HGCalDDDConstants* ddd_;
};
#endif
1 change: 1 addition & 0 deletions RecoLocalCalo/HGCalRecProducers/BuildFile.xml
View file
Open in desktop
@@ -1,5 +1,6 @@
<use name="FWCore/MessageLogger"/>
<use name="FWCore/Framework"/>
<use name="Geometry/HGCalGeometry"/>
<use name="clhep"/>
<export>
<lib name="1"/>
Expand Down
1 change: 1 addition & 0 deletions RecoLocalCalo/HGCalRecProducers/plugins/BuildFile.xml
View file
Open in desktop
Expand Up @@ -8,6 +8,7 @@
<use name="CondFormats/DataRecord"/>
<use name="FWCore/MessageLogger"/>
<use name="FWCore/MessageService"/>
<use name="Geometry/HGCalGeometry"/>
<library file="*.cc" name="RecoLocalCaloHGCalRecProducersPlugins">
<flags EDM_PLUGIN="1"/>
</library>
22 changes: 21 additions & 1 deletion RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.cc
View file
Open in desktop
Expand Up @@ -13,10 +13,12 @@ HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet&ps) :
constexpr float keV2GeV = 1e-6;
// HGCee constants
HGCEE_keV2DIGI_ = ps.getParameter<double>("HGCEE_keV2DIGI");
HGCEE_fCPerMIP_ = ps.getParameter<std::vector<double> >("HGCEE_fCPerMIP");
hgceeUncalib2GeV_ = keV2GeV/HGCEE_keV2DIGI_;

// HGChef constants
HGCHEF_keV2DIGI_ = ps.getParameter<double>("HGCHEF_keV2DIGI");
HGCHEF_fCPerMIP_ = ps.getParameter<std::vector<double> >("HGCHEF_fCPerMIP");
hgchefUncalib2GeV_ = keV2GeV/HGCHEF_keV2DIGI_;

// HGCheb constants
Expand All @@ -25,6 +27,13 @@ HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet&ps) :
}

void HGCalRecHitWorkerSimple::set(const edm::EventSetup& es) {
edm::ESHandle<HGCalGeometry> hgceeGeoHandle;
edm::ESHandle<HGCalGeometry> hgchefGeoHandle;
es.get<IdealGeometryRecord>().get("HGCalEESensitive",hgceeGeoHandle);
es.get<IdealGeometryRecord>().get("HGCalHESiliconSensitive",hgchefGeoHandle);
ddds_[0] = &(hgceeGeoHandle->topology().dddConstants());
ddds_[1] = &(hgchefGeoHandle->topology().dddConstants());
ddds_[2] = nullptr;
}


Expand All @@ -34,13 +43,16 @@ HGCalRecHitWorkerSimple::run( const edm::Event & evt,
HGCRecHitCollection & result ) {
DetId detid=uncalibRH.id();
uint32_t recoFlag = 0;
const std::vector<double>* fCPerMIP = nullptr;

switch( detid.subdetId() ) {
case HGCEE:
rechitMaker_->setADCToGeVConstant(float(hgceeUncalib2GeV_) );
fCPerMIP = &HGCEE_fCPerMIP_;
break;
case HGCHEF:
rechitMaker_->setADCToGeVConstant(float(hgchefUncalib2GeV_) );
fCPerMIP = &HGCHEF_fCPerMIP_;
break;
case HGCHEB:
rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_) );
Expand All @@ -51,8 +63,16 @@ HGCalRecHitWorkerSimple::run( const edm::Event & evt,
}

// make the rechit and put in the output collection
if (recoFlag == 0) {
if (recoFlag == 0) {
HGCRecHit myrechit( rechitMaker_->makeRecHit(uncalibRH, 0) );
HGCalDetId hid(detid);
if( fCPerMIP != nullptr ) {
const int thk = ddds_[hid.subdetId()-3]->waferTypeL(hid.wafer());
// units out of rechit maker are MIP * (GeV/fC)
// so multiple
const double new_E = myrechit.energy()*(*fCPerMIP)[thk-1];
myrechit.setEnergy(new_E);
}
result.push_back(myrechit);
}

Expand Down
6 changes: 6 additions & 0 deletions RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.h
View file
Open in desktop
Expand Up @@ -13,6 +13,8 @@

#include "FWCore/Framework/interface/ESHandle.h"

#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"

class HGCalRecHitWorkerSimple : public HGCalRecHitWorkerBaseClass {
public:
HGCalRecHitWorkerSimple(const edm::ParameterSet&);
Expand All @@ -24,9 +26,13 @@ class HGCalRecHitWorkerSimple : public HGCalRecHitWorkerBaseClass {
protected:

double HGCEE_keV2DIGI_, hgceeUncalib2GeV_;
std::vector<double> HGCEE_fCPerMIP_;
double HGCHEF_keV2DIGI_, hgchefUncalib2GeV_;
std::vector<double> HGCHEF_fCPerMIP_;
double HGCHEB_keV2DIGI_, hgchebUncalib2GeV_;

std::array<const HGCalDDDConstants*, 3> ddds_;

std::vector<int> v_chstatus_;

std::vector<int> v_DB_reco_flags_;
Expand Down
16 changes: 15 additions & 1 deletion RecoLocalCalo/HGCalRecProducers/plugins/HGCalUncalibRecHitWorkerWeights.cc
View file
Open in desktop
Expand Up @@ -16,6 +16,7 @@ void configureIt(const edm::ParameterSet& conf,
constexpr char tdcSaturation[] = "tdcSaturation";
constexpr char tdcOnset[] = "tdcOnset";
constexpr char toaLSB_ns[] = "toaLSB_ns";
constexpr char fCPerMIP[] = "fCPerMIP";

if( conf.exists(isSiFE) ) {
maker.set_isSiFESim(conf.getParameter<bool>(isSiFE));
Expand Down Expand Up @@ -49,6 +50,13 @@ void configureIt(const edm::ParameterSet& conf,
} else {
maker.set_toaLSBToNS(-1.);
}

if( conf.exists(fCPerMIP) ) {
maker.set_fCPerMIP(conf.getParameter<std::vector<double> >(fCPerMIP));
} else {
maker.set_fCPerMIP(std::vector<double>({1.0}));
}

}

HGCalUncalibRecHitWorkerWeights::HGCalUncalibRecHitWorkerWeights(const edm::ParameterSet&ps) :
Expand All @@ -65,7 +73,13 @@ HGCalUncalibRecHitWorkerWeights::HGCalUncalibRecHitWorkerWeights(const edm::Para
void
HGCalUncalibRecHitWorkerWeights::set(const edm::EventSetup& es)
{

edm::ESHandle<HGCalGeometry> hgceeGeoHandle;
edm::ESHandle<HGCalGeometry> hgchefGeoHandle;
es.get<IdealGeometryRecord>().get("HGCalEESensitive",hgceeGeoHandle) ;
es.get<IdealGeometryRecord>().get("HGCalHESiliconSensitive",hgchefGeoHandle) ;
uncalibMaker_ee_.setGeometry(hgceeGeoHandle.product());
uncalibMaker_hef_.setGeometry(hgchefGeoHandle.product());
uncalibMaker_heb_.setGeometry(nullptr);
}


Expand Down
4 changes: 4 additions & 0 deletions RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py
View file
Open in desktop
@@ -1,6 +1,8 @@
import FWCore.ParameterSet.Config as cms

from SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi import *
from RecoLocalCalo.HGCalRecProducers.HGCalUncalibRecHit_cfi import *

# HGCAL rechit producer
HGCalRecHit = cms.EDProducer(
"HGCalRecHitProducer",
Expand All @@ -13,7 +15,9 @@

# digi constants
HGCEE_keV2DIGI = hgceeDigitizer.digiCfg.keV2fC,
HGCEE_fCPerMIP = HGCalUncalibRecHit.HGCEEConfig.fCPerMIP,
HGCHEF_keV2DIGI = hgchefrontDigitizer.digiCfg.keV2fC,
HGCHEF_fCPerMIP = HGCalUncalibRecHit.HGCHEFConfig.fCPerMIP,
HGCHEB_keV2DIGI = hgchebackDigitizer.digiCfg.keV2MIP,

# algo
Expand Down
6 changes: 4 additions & 2 deletions RecoLocalCalo/HGCalRecProducers/python/HGCalUncalibRecHit_cfi.py
View file
Open in desktop
Expand Up @@ -21,7 +21,8 @@
tdcNbits = hgceeDigitizer.digiCfg.feCfg.tdcNbits,
tdcSaturation = hgceeDigitizer.digiCfg.feCfg.tdcSaturation_fC,
tdcOnset = hgceeDigitizer.digiCfg.feCfg.tdcOnset_fC,
toaLSB_ns = hgceeDigitizer.digiCfg.feCfg.toaLSB_ns
toaLSB_ns = hgceeDigitizer.digiCfg.feCfg.toaLSB_ns,
fCPerMIP = cms.vdouble(1.25,2.57,3.88) #100um, 200um, 300um
),

HGCHEFConfig = cms.PSet(
Expand All @@ -33,7 +34,8 @@
tdcNbits = hgchefrontDigitizer.digiCfg.feCfg.tdcNbits,
tdcSaturation = hgchefrontDigitizer.digiCfg.feCfg.tdcSaturation_fC,
tdcOnset = hgchefrontDigitizer.digiCfg.feCfg.tdcOnset_fC,
toaLSB_ns = hgchefrontDigitizer.digiCfg.feCfg.toaLSB_ns
toaLSB_ns = hgchefrontDigitizer.digiCfg.feCfg.toaLSB_ns,
fCPerMIP = cms.vdouble(1.25,2.57,3.88) #100um, 200um, 300um
),

HGCHEBConfig = cms.PSet(
Expand Down
4 changes: 2 additions & 2 deletions RecoLocalCalo/HGCalRecProducers/test/testHGCalRecoLocal_cfg.py
View file
Open in desktop
Expand Up @@ -11,8 +11,8 @@
process.load('FWCore.MessageService.MessageLogger_cfi')
process.load('Configuration.EventContent.EventContent_cff')
process.load('SimGeneral.MixingModule.mixNoPU_cfi')
process.load('Configuration.Geometry.GeometryExtended2023Reco_cff')
process.load('Configuration.Geometry.GeometryExtended2023_cff')
process.load('Configuration.Geometry.GeometryExtended2023DevReco_cff')
process.load('Configuration.Geometry.GeometryExtended2023Dev_cff')
process.load('Configuration.StandardSequences.MagneticField_38T_PostLS1_cff')
process.load('Configuration.StandardSequences.Generator_cff')
process.load('IOMC.EventVertexGenerators.VtxSmearedGauss_cfi')
Expand Down
76 changes: 76 additions & 0 deletions RecoParticleFlow/PFClusterProducer/interface/PFRecHitQTests.h
View file
Open in desktop
Expand Up @@ -638,4 +638,80 @@ class PFRecHitQTestThresholdInMIPs : public PFRecHitQTestBase {
}
};

#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
class PFRecHitQTestThresholdInThicknessNormalizedMIPs : public PFRecHitQTestBase {
public:
PFRecHitQTestThresholdInThicknessNormalizedMIPs() :
geometryInstance_(""),
recHitEnergy_keV_(0.),
threshold_(0.),
mip_(0.),
recHitEnergyMultiplier_(0.) {
}

PFRecHitQTestThresholdInThicknessNormalizedMIPs(const edm::ParameterSet& iConfig) :
PFRecHitQTestBase(iConfig),
geometryInstance_(iConfig.getParameter<std::string>("geometryInstance")),
recHitEnergy_keV_(iConfig.getParameter<bool>("recHitEnergyIs_keV")),
threshold_(iConfig.getParameter<double>("thresholdInMIPs")),
mip_(iConfig.getParameter<double>("mipValueInkeV")),
recHitEnergyMultiplier_(iConfig.getParameter<double>("recHitEnergyMultiplier")) {
}

void beginEvent(const edm::Event& event,const edm::EventSetup& iSetup) {
edm::ESHandle<HGCalGeometry> geoHandle;
iSetup.get<IdealGeometryRecord>().get(geometryInstance_,geoHandle);
ddd_ = &(geoHandle->topology().dddConstants());
}

bool test(reco::PFRecHit& hit,const EcalRecHit& rh,bool& clean) {
throw cms::Exception("WrongDetector")
<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
return false;
}
bool test(reco::PFRecHit& hit,const HBHERecHit& rh,bool& clean) {
throw cms::Exception("WrongDetector")
<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
return false;
}

bool test(reco::PFRecHit& hit,const HFRecHit& rh,bool& clean) {
throw cms::Exception("WrongDetector")
<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
return false;
}
bool test(reco::PFRecHit& hit,const HORecHit& rh,bool& clean) {
throw cms::Exception("WrongDetector")
<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
return false;
}

bool test(reco::PFRecHit& hit,const CaloTower& rh,bool& clean) {
throw cms::Exception("WrongDetector")
<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
return false;
}

bool test(reco::PFRecHit& hit,const HGCRecHit& rh,bool& clean) {
const double newE = ( recHitEnergy_keV_ ?
1.0e-6*rh.energy()*recHitEnergyMultiplier_ :
rh.energy()*recHitEnergyMultiplier_ );
const int wafer = HGCalDetId(rh.detid()).wafer();
const float mult = (float) ddd_->waferTypeL(wafer); // 1 for 100um, 2 for 200um, 3 for 300um
hit.setEnergy(newE);
return pass(hit,mult);
}

protected:
const std::string geometryInstance_;
const bool recHitEnergy_keV_;
const double threshold_,mip_,recHitEnergyMultiplier_;
const HGCalDDDConstants* ddd_;

bool pass(const reco::PFRecHit& hit, const float mult) {
const double hitValueInMIPs = 1e6*hit.energy()/(mult*mip_);
return hitValueInMIPs > threshold_;
}
};

#endif