Skip to content

Commit

Permalink
Implement particle-level truth for PF
Browse files Browse the repository at this point in the history
  • Loading branch information
@jpata
jpata committed Mar 31, 2020
1 parent 887d27e commit a6d6ae1
Show file tree
Hide file tree
Showing 7 changed files with 1,277 additions and 50 deletions.
5 changes: 3 additions & 2 deletions RecoParticleFlow/PFProducer/python/simPFProducer_cfi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,9 @@
gsfTrackTimeErrorMap = cms.InputTag("gsfTrackTimeValueMapProducer:electronGsfTracksConfigurableFlatResolutionModelResolution"),
)

from Configuration.Eras.Modifier_phase2_timing_layer_cff import phase2_timing_layer
phase2_timing_layer.toModify(
from Configuration.Eras.Modifier_phase2_timing_layer_tile_cff import phase2_timing_layer_tile
from Configuration.Eras.Modifier_phase2_timing_layer_bar_cff import phase2_timing_layer_bar
(phase2_timing_layer_tile | phase2_timing_layer_bar).toModify(
simPFProducer,
trackTimeValueMap = cms.InputTag("tofPID:t0"),
trackTimeErrorMap = cms.InputTag("tofPID:sigmat0"),
Expand Down
116 changes: 68 additions & 48 deletions SimGeneral/CaloAnalysis/plugins/CaloTruthAccumulator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -206,6 +206,7 @@ class CaloTruthAccumulator : public DigiAccumulatorMixMod {
calo_particles m_caloParticles;
// geometry type (0 pre-TDR; 1 TDR)
int geometryType_;
bool doHGCAL;
};

/* Graph utility functions */
Expand Down Expand Up @@ -368,7 +369,8 @@ CaloTruthAccumulator::CaloTruthAccumulator(const edm::ParameterSet &config,
minEnergy_(config.getParameter<double>("MinEnergy")),
maxPseudoRapidity_(config.getParameter<double>("MaxPseudoRapidity")),
premixStage1_(config.getParameter<bool>("premixStage1")),
geometryType_(-1) {
geometryType_(-1),
doHGCAL(config.getParameter<bool>("doHGCAL")) {
producesCollector.produces<SimClusterCollection>("MergedCaloTruth");
producesCollector.produces<CaloParticleCollection>("MergedCaloTruth");
if (premixStage1_) {
Expand Down Expand Up @@ -400,34 +402,38 @@ void CaloTruthAccumulator::beginLuminosityBlock(edm::LuminosityBlock const &iLum
iSetup.get<CaloGeometryRecord>().get(geom);
const HGCalGeometry *eegeom = nullptr, *fhgeom = nullptr, *bhgeomnew = nullptr;
const HcalGeometry *bhgeom = nullptr;

eegeom = static_cast<const HGCalGeometry *>(
geom->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
// check if it's the new geometry
if (eegeom) {
geometryType_ = 1;
fhgeom = static_cast<const HGCalGeometry *>(
geom->getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
bhgeomnew = static_cast<const HGCalGeometry *>(
geom->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
} else {
geometryType_ = 0;
eegeom = static_cast<const HGCalGeometry *>(geom->getSubdetectorGeometry(DetId::Forward, HGCEE));
fhgeom = static_cast<const HGCalGeometry *>(geom->getSubdetectorGeometry(DetId::Forward, HGCHEF));
bhgeom = static_cast<const HcalGeometry *>(geom->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
}
hgtopo_[0] = &(eegeom->topology());
hgtopo_[1] = &(fhgeom->topology());
if (bhgeomnew)
hgtopo_[2] = &(bhgeomnew->topology());

for (unsigned i = 0; i < 3; ++i) {
if (hgtopo_[i])
hgddd_[i] = &(hgtopo_[i]->dddConstants());
bhgeom = static_cast<const HcalGeometry *>(geom->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));

if (doHGCAL) {
eegeom = static_cast<const HGCalGeometry *>(
geom->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
// check if it's the new geometry
if (eegeom) {
geometryType_ = 1;
fhgeom = static_cast<const HGCalGeometry *>(
geom->getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
bhgeomnew = static_cast<const HGCalGeometry *>(
geom->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
} else {
geometryType_ = 0;
eegeom = static_cast<const HGCalGeometry *>(geom->getSubdetectorGeometry(DetId::Forward, HGCEE));
fhgeom = static_cast<const HGCalGeometry *>(geom->getSubdetectorGeometry(DetId::Forward, HGCHEF));
bhgeom = static_cast<const HcalGeometry *>(geom->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
}
hgtopo_[0] = &(eegeom->topology());
hgtopo_[1] = &(fhgeom->topology());
if (bhgeomnew)
hgtopo_[2] = &(bhgeomnew->topology());

for (unsigned i = 0; i < 3; ++i) {
if (hgtopo_[i])
hgddd_[i] = &(hgtopo_[i]->dddConstants());
}
}

if (bhgeom)
if (bhgeom) {
hcddd_ = bhgeom->topology().dddConstants();
}
}

void CaloTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSetup const &setup) {
Expand Down Expand Up @@ -672,37 +678,51 @@ void CaloTruthAccumulator::fillSimHits(std::vector<std::pair<DetId, const PCaloH
edm::Handle<std::vector<PCaloHit>> hSimHits;
const bool isHcal = (collectionTag.instance().find("HcalHits") != std::string::npos);
event.getByLabel(collectionTag, hSimHits);

for (auto const &simHit : *hSimHits) {
DetId id(0);
const uint32_t simId = simHit.id();
if (geometryType_ == 1) {
// no test numbering in new geometry
id = simId;
} else if (isHcal) {
HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd_);
if (hid.subdet() == HcalEndcap)
id = hid;

//Relabel as necessary for HGCAL
if (doHGCAL) {
const uint32_t simId = simHit.id();
if (geometryType_ == 1) {
// no test numbering in new geometry
id = simId;
} else if (isHcal) {
HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd_);
if (hid.subdet() == HcalEndcap)
id = hid;
} else {
int subdet, layer, cell, sec, subsec, zp;
HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
const HGCalDDDConstants *ddd = hgddd_[subdet - 3];
std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo_[subdet - 3]->detectorType());
cell = recoLayerCell.first;
layer = recoLayerCell.second;
// skip simhits with bad barcodes or non-existant layers
if (layer == -1 || simHit.geantTrackId() == 0)
continue;
id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
}
} else {
int subdet, layer, cell, sec, subsec, zp;
HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
const HGCalDDDConstants *ddd = hgddd_[subdet - 3];
std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo_[subdet - 3]->detectorType());
cell = recoLayerCell.first;
layer = recoLayerCell.second;
// skip simhits with bad barcodes or non-existant layers
if (layer == -1 || simHit.geantTrackId() == 0)
continue;
id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
id = simHit.id();
//Relabel all HCAL hits
if (isHcal) {
HcalDetId hid = HcalHitRelabeller::relabel(simHit.id(), hcddd_);
id = hid;
}
}

if (DetId(0) == id)
if (id == DetId(0)) {
continue;
}
if (simHit.geantTrackId() == 0) {
continue;
}

uint32_t detId = id.rawId();
returnValue.emplace_back(id, &simHit);
simTrackDetIdEnergyMap[simHit.geantTrackId()][id.rawId()] += simHit.energy();

m_detIdToTotalSimEnergy[detId] += simHit.energy();
m_detIdToTotalSimEnergy[id.rawId()] += simHit.energy();
}
} // end of loop over InputTags
}
Expand Down
1 change: 1 addition & 0 deletions SimGeneral/MixingModule/python/caloTruthProducer_cfi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
MinEnergy = cms.double(0.5),
MaxPseudoRapidity = cms.double(5.0),
premixStage1 = cms.bool(False),
doHGCAL = cms.bool(True),
maximumPreviousBunchCrossing = cms.uint32(0),
maximumSubsequentBunchCrossing = cms.uint32(0),
simHitCollections = cms.PSet(
Expand Down
2 changes: 2 additions & 0 deletions Validation/RecoParticleFlow/plugins/BuildFile.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@
<use name="DataFormats/ParticleFlowReco"/>
<use name="DataFormats/ParticleFlowCandidate"/>
<use name="RecoParticleFlow/Benchmark"/>
<use name="CommonTools/UtilAlgos"/>
<use name="Geometry/HcalTowerAlgo"/>
<flags EDM_PLUGIN="1"/>
</library>
<export>
Expand Down
Loading

0 comments on commit a6d6ae1

Please sign in to comment.