Back to home page

Project CMSSW displayed by LXR

 
 

    


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;                   // stop default
0029   const FWTracksterLayersProxyBuilder &operator=(const FWTracksterLayersProxyBuilder &) = delete;  // stop default
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     // Apply heatmap color coding **after** the call to setupAddElement, that will internally setup the color.
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     // seed and cluster position
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       // draw 3D cross
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);