File indexing completed on 2023-10-25 09:46:16
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "Fireworks/Calo/interface/FWHeatmapProxyBuilderTemplate.h"
0003 #include "Fireworks/Core/interface/FWEventItem.h"
0004 #include "Fireworks/Core/interface/FWGeometry.h"
0005 #include "Fireworks/Core/interface/BuilderUtils.h"
0006 #include "DataFormats/HGCalReco/interface/Trackster.h"
0007 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0008 #include "DataFormats/Common/interface/ValueMap.h"
0009 #include "DataFormats/DetId/interface/DetId.h"
0010 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0011 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0012
0013 #include "TEveBoxSet.h"
0014 #include "TGeoSphere.h"
0015 #include "TGeoTube.h"
0016 #include "TEveGeoShape.h"
0017 #include "TEveStraightLineSet.h"
0018
0019 #include <cmath>
0020
0021 class FWTracksterLayersProxyBuilder : public FWHeatmapProxyBuilderTemplate<ticl::Trackster> {
0022 public:
0023 FWTracksterLayersProxyBuilder(void) {}
0024 ~FWTracksterLayersProxyBuilder(void) override {}
0025
0026 REGISTER_PROXYBUILDER_METHODS();
0027
0028 FWTracksterLayersProxyBuilder(const FWTracksterLayersProxyBuilder &) = delete;
0029 const FWTracksterLayersProxyBuilder &operator=(const FWTracksterLayersProxyBuilder &) = delete;
0030
0031 private:
0032 edm::Handle<edm::ValueMap<std::pair<float, float>>> TimeValueMapHandle_;
0033 edm::Handle<std::vector<reco::CaloCluster>> layerClustersHandle_;
0034 double timeLowerBound_, timeUpperBound_;
0035 long layer_;
0036 double saturation_energy_;
0037 bool heatmap_;
0038 bool z_plus_;
0039 bool z_minus_;
0040 bool enableTimeFilter_;
0041 bool enablePositionLines_;
0042 bool enableEdges_;
0043 double displayMode_;
0044 double proportionalityFactor_;
0045
0046 void setItem(const FWEventItem *iItem) override;
0047
0048 void build(const FWEventItem *iItem, TEveElementList *product, const FWViewContext *vc) override;
0049 void build(const ticl::Trackster &iData,
0050 unsigned int iIndex,
0051 TEveElement &oItemHolder,
0052 const FWViewContext *) override;
0053 };
0054
0055 void FWTracksterLayersProxyBuilder::setItem(const FWEventItem *iItem) {
0056 FWHeatmapProxyBuilderTemplate::setItem(iItem);
0057 if (iItem) {
0058 iItem->getConfig()->assertParam("EnablePositionLines", false);
0059 iItem->getConfig()->assertParam("EnableEdges", false);
0060 iItem->getConfig()->assertParam("EnableTimeFilter", false);
0061 iItem->getConfig()->assertParam("TimeLowerBound(ns)", 0.01, 0.0, 75.0);
0062 iItem->getConfig()->assertParam("TimeUpperBound(ns)", 0.01, 0.0, 75.0);
0063 iItem->getConfig()->assertParam("DisplayMode", 0.0, 0.0, 5.0);
0064 iItem->getConfig()->assertParam("ProportionalityFactor", 1.0, 0.0, 1.0);
0065 }
0066 }
0067
0068 void FWTracksterLayersProxyBuilder::build(const FWEventItem *iItem, TEveElementList *product, const FWViewContext *vc) {
0069 iItem->getEvent()->getByLabel(edm::InputTag("hgcalLayerClusters", "timeLayerCluster"), TimeValueMapHandle_);
0070 iItem->getEvent()->getByLabel(edm::InputTag("hgcalLayerClusters"), layerClustersHandle_);
0071 if (TimeValueMapHandle_.isValid()) {
0072 timeLowerBound_ = item()->getConfig()->value<double>("TimeLowerBound(ns)");
0073 timeUpperBound_ = item()->getConfig()->value<double>("TimeUpperBound(ns)");
0074 if (timeLowerBound_ > timeUpperBound_) {
0075 edm::LogWarning("InvalidParameters")
0076 << "lower time bound is larger than upper time bound. Maybe opposite is desired?";
0077 }
0078 } else {
0079 iItem->getEvent()->getByLabel(edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster"), TimeValueMapHandle_);
0080 edm::LogWarning("DataNotFound|InvalidData")
0081 << __FILE__ << ":" << __LINE__
0082 << " couldn't locate 'hgcalLayerClusters:timeLayerCluster' ValueMap in input file. Trying to access "
0083 "'hgcalMergeLayerClusters:timeLayerClusters' ValueMap";
0084 if (!TimeValueMapHandle_.isValid()) {
0085 edm::LogWarning("DataNotFound|InvalidData")
0086 << __FILE__ << ":" << __LINE__
0087 << " couldn't locate 'hgcalMergeLayerClusters:timeLayerCluster' ValueMap in input file.";
0088 }
0089 }
0090
0091 if (!layerClustersHandle_.isValid()) {
0092 iItem->getEvent()->getByLabel(edm::InputTag("hgcalMergeLayerClusters"), layerClustersHandle_);
0093 edm::LogWarning("DataNotFound|InvalidData")
0094 << __FILE__ << ":" << __LINE__
0095 << " couldn't locate 'hgcalLayerClusters' collection "
0096 "in input file. Trying to access 'hgcalMergeLayerClusters' collection.";
0097 if (!layerClustersHandle_.isValid()) {
0098 edm::LogWarning("DataNotFound|InvalidData")
0099 << __FILE__ << ":" << __LINE__ << " couldn't locate 'hgcalMergeLayerClusters' collection in input file.";
0100 }
0101 }
0102
0103 layer_ = item()->getConfig()->value<long>("Layer");
0104 saturation_energy_ = item()->getConfig()->value<double>("EnergyCutOff");
0105 heatmap_ = item()->getConfig()->value<bool>("Heatmap");
0106 z_plus_ = item()->getConfig()->value<bool>("Z+");
0107 z_minus_ = item()->getConfig()->value<bool>("Z-");
0108 enableTimeFilter_ = item()->getConfig()->value<bool>("EnableTimeFilter");
0109 enablePositionLines_ = item()->getConfig()->value<bool>("EnablePositionLines");
0110 enableEdges_ = item()->getConfig()->value<bool>("EnableEdges");
0111 displayMode_ = item()->getConfig()->value<double>("DisplayMode");
0112 proportionalityFactor_ = item()->getConfig()->value<double>("ProportionalityFactor");
0113
0114 FWHeatmapProxyBuilderTemplate::build(iItem, product, vc);
0115 }
0116
0117 void FWTracksterLayersProxyBuilder::build(const ticl::Trackster &iData,
0118 unsigned int iIndex,
0119 TEveElement &oItemHolder,
0120 const FWViewContext *) {
0121 if (enableTimeFilter_ && TimeValueMapHandle_.isValid()) {
0122 const float time = TimeValueMapHandle_->get(iIndex).first;
0123 if (time < timeLowerBound_ || time > timeUpperBound_)
0124 return;
0125 }
0126
0127 const ticl::Trackster &trackster = iData;
0128 const size_t N = trackster.vertices().size();
0129 const std::vector<reco::CaloCluster> &layerClusters = *layerClustersHandle_;
0130 TEveStraightLineSet *position_marker = nullptr;
0131
0132 if (enablePositionLines_) {
0133 position_marker = new TEveStraightLineSet;
0134 position_marker->SetLineWidth(2);
0135 position_marker->SetLineColor(kWhite);
0136 }
0137
0138 for (size_t i = 0; i < N; ++i) {
0139 const reco::CaloCluster layerCluster = layerClusters[trackster.vertices(i)];
0140 const math::XYZPoint &position = layerCluster.position();
0141 const size_t nHits = layerCluster.size();
0142 const double energy = layerCluster.energy();
0143 float radius = 0;
0144 auto detIdOnLayer = layerCluster.seed();
0145
0146 const auto *parameters = item()->getGeom()->getParameters(detIdOnLayer);
0147 const int layer = parameters[1];
0148 const int zside = parameters[2];
0149 const bool isSilicon = parameters[3];
0150
0151 auto const z_selection_is_on = z_plus_ ^ z_minus_;
0152 auto const z_plus_selection_ok = z_plus_ && (zside == 1);
0153 auto const z_minus_selection_ok = z_minus_ && (zside == -1);
0154 if (!z_minus_ && !z_plus_)
0155 continue;
0156 if (z_selection_is_on && !(z_plus_selection_ok || z_minus_selection_ok))
0157 continue;
0158
0159 if (layer_ > 0 && (layer != layer_))
0160 continue;
0161
0162 if (displayMode_ == 0) {
0163 radius = sqrt(nHits);
0164 } else if (displayMode_ == 1) {
0165 radius = nHits;
0166 } else if (displayMode_ == 2) {
0167 radius = energy;
0168 } else if (displayMode_ == 3) {
0169 radius = energy / nHits;
0170 } else if (displayMode_ == 4) {
0171 float area = 0;
0172 if (!isSilicon) {
0173 const bool isFine = (HGCScintillatorDetId(layerCluster.seed()).type() == 0);
0174 float dphi = (isFine) ? 1.0 * M_PI / 180. : 1.25 * M_PI / 180.;
0175 int ir = HGCScintillatorDetId(layerCluster.seed()).iradiusAbs();
0176 float dr = (isFine) ? (0.0484 * static_cast<float>(ir) + 2.1) : (0.075 * static_cast<float>(ir) + 2.0);
0177 float r = (isFine) ? (0.0239 * static_cast<float>(pow(ir, 2)) + 2.02 * static_cast<float>(ir) + 119.6)
0178 : (0.0367 * static_cast<float>(pow(ir, 2)) + 1.7 * static_cast<float>(ir) + 90.7);
0179 area = r * dr * dphi;
0180 } else {
0181 const bool isFine = (HGCSiliconDetId(layerCluster.seed()).type() == 0);
0182 float side = (isFine) ? 0.465 : 0.698;
0183 area = pow(side, 2) * 3 * sqrt(3) / 2;
0184 }
0185 radius = sqrt(nHits * area) / M_PI;
0186 }
0187
0188 auto *eveCircle = new TEveGeoShape("Circle");
0189 auto tube = new TGeoTube(0., proportionalityFactor_ * radius, 0.1);
0190 eveCircle->SetShape(tube);
0191 eveCircle->InitMainTrans();
0192 eveCircle->RefMainTrans().Move3PF(position.x(), position.y(), position.z());
0193 setupAddElement(eveCircle, &oItemHolder);
0194
0195 if (heatmap_) {
0196 const float normalized_energy = fmin(energy / saturation_energy_, 1.0f);
0197 const uint8_t colorFactor = gradient_steps * normalized_energy;
0198 eveCircle->SetFillColor(
0199 TColor::GetColor(gradient[0][colorFactor], gradient[1][colorFactor], gradient[2][colorFactor]));
0200 } else {
0201 eveCircle->SetMainColor(item()->modelInfo(iIndex).displayProperties().color());
0202 eveCircle->SetMainTransparency(item()->defaultDisplayProperties().transparency());
0203 }
0204
0205
0206 const float crossScale = 1.0f + fmin(energy, 5.0f);
0207 if (enablePositionLines_) {
0208 auto const &pos = layerCluster.position();
0209 const float position_crossScale = crossScale * 0.5;
0210 position_marker->AddLine(
0211 pos.x() - position_crossScale, pos.y(), pos.z(), pos.x() + position_crossScale, pos.y(), pos.z());
0212 position_marker->AddLine(
0213 pos.x(), pos.y() - position_crossScale, pos.z(), pos.x(), pos.y() + position_crossScale, pos.z());
0214 }
0215 }
0216
0217 if (enablePositionLines_)
0218 oItemHolder.AddElement(position_marker);
0219
0220 if (enableEdges_) {
0221 auto &edges = trackster.edges();
0222
0223 TEveStraightLineSet *adjacent_marker = new TEveStraightLineSet;
0224 adjacent_marker->SetLineWidth(2);
0225 adjacent_marker->SetLineColor(kYellow);
0226
0227 TEveStraightLineSet *non_adjacent_marker = new TEveStraightLineSet;
0228 non_adjacent_marker->SetLineWidth(2);
0229 non_adjacent_marker->SetLineColor(kRed);
0230
0231 for (auto edge : edges) {
0232 auto doublet = std::make_pair(layerClusters[edge[0]], layerClusters[edge[1]]);
0233
0234 int layerIn = item()->getGeom()->getParameters(doublet.first.seed())[1];
0235 int layerOut = item()->getGeom()->getParameters(doublet.second.seed())[1];
0236
0237 const bool isAdjacent = std::abs(layerOut - layerIn) == 1;
0238
0239
0240 if (layer_ == 0 || fabs(layerIn - layer_) == 0 || fabs(layerOut - layer_) == 0) {
0241 if (isAdjacent)
0242 adjacent_marker->AddLine(doublet.first.x(),
0243 doublet.first.y(),
0244 doublet.first.z(),
0245 doublet.second.x(),
0246 doublet.second.y(),
0247 doublet.second.z());
0248 else
0249 non_adjacent_marker->AddLine(doublet.first.x(),
0250 doublet.first.y(),
0251 doublet.first.z(),
0252 doublet.second.x(),
0253 doublet.second.y(),
0254 doublet.second.z());
0255 }
0256 }
0257 oItemHolder.AddElement(adjacent_marker);
0258 oItemHolder.AddElement(non_adjacent_marker);
0259 }
0260 }
0261
0262 REGISTER_FWPROXYBUILDER(FWTracksterLayersProxyBuilder, ticl::Trackster, "Trackster layers", FWViewType::kISpyBit);
