Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-13 23:03:09

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     edm::LogWarning("DataNotFound|InvalidData") << "couldn't locate 'timeLayerCluster' ValueMap in root file.";
0080   }
0081 
0082   if (!layerClustersHandle_.isValid()) {
0083     edm::LogWarning("DataNotFound|InvalidData") << "couldn't locate 'timeLayerCluster' ValueMap in root file.";
0084   }
0085 
0086   layer_ = item()->getConfig()->value<long>("Layer");
0087   saturation_energy_ = item()->getConfig()->value<double>("EnergyCutOff");
0088   heatmap_ = item()->getConfig()->value<bool>("Heatmap");
0089   z_plus_ = item()->getConfig()->value<bool>("Z+");
0090   z_minus_ = item()->getConfig()->value<bool>("Z-");
0091   enableTimeFilter_ = item()->getConfig()->value<bool>("EnableTimeFilter");
0092   enablePositionLines_ = item()->getConfig()->value<bool>("EnablePositionLines");
0093   enableEdges_ = item()->getConfig()->value<bool>("EnableEdges");
0094   displayMode_ = item()->getConfig()->value<double>("DisplayMode");
0095   proportionalityFactor_ = item()->getConfig()->value<double>("ProportionalityFactor");
0096 
0097   FWHeatmapProxyBuilderTemplate::build(iItem, product, vc);
0098 }
0099 
0100 void FWTracksterLayersProxyBuilder::build(const ticl::Trackster &iData,
0101                                           unsigned int iIndex,
0102                                           TEveElement &oItemHolder,
0103                                           const FWViewContext *) {
0104   if (enableTimeFilter_ && TimeValueMapHandle_.isValid()) {
0105     const float time = TimeValueMapHandle_->get(iIndex).first;
0106     if (time < timeLowerBound_ || time > timeUpperBound_)
0107       return;
0108   }
0109 
0110   const ticl::Trackster &trackster = iData;
0111   const size_t N = trackster.vertices().size();
0112   const std::vector<reco::CaloCluster> &layerClusters = *layerClustersHandle_;
0113   TEveStraightLineSet *position_marker = nullptr;
0114 
0115   if (enablePositionLines_) {
0116     position_marker = new TEveStraightLineSet;
0117     position_marker->SetLineWidth(2);
0118     position_marker->SetLineColor(kWhite);
0119   }
0120 
0121   for (size_t i = 0; i < N; ++i) {
0122     const reco::CaloCluster layerCluster = layerClusters[trackster.vertices(i)];
0123     const math::XYZPoint &position = layerCluster.position();
0124     const size_t nHits = layerCluster.size();
0125     const double energy = layerCluster.energy();
0126     float radius = 0;
0127     auto detIdOnLayer = layerCluster.seed();
0128 
0129     const auto *parameters = item()->getGeom()->getParameters(detIdOnLayer);
0130     const int layer = parameters[1];
0131     const int zside = parameters[2];
0132     const bool isSilicon = parameters[3];
0133 
0134     auto const z_selection_is_on = z_plus_ ^ z_minus_;
0135     auto const z_plus_selection_ok = z_plus_ && (zside == 1);
0136     auto const z_minus_selection_ok = z_minus_ && (zside == -1);
0137     if (!z_minus_ && !z_plus_)
0138       continue;
0139     if (z_selection_is_on && !(z_plus_selection_ok || z_minus_selection_ok))
0140       continue;
0141 
0142     if (layer_ > 0 && (layer != layer_))
0143       continue;
0144 
0145     if (displayMode_ == 0) {
0146       radius = sqrt(nHits);
0147     } else if (displayMode_ == 1) {
0148       radius = nHits;
0149     } else if (displayMode_ == 2) {
0150       radius = energy;
0151     } else if (displayMode_ == 3) {
0152       radius = energy / nHits;
0153     } else if (displayMode_ == 4) {
0154       float area = 0;
0155       if (!isSilicon) {
0156         const bool isFine = (HGCScintillatorDetId(layerCluster.seed()).type() == 0);
0157         float dphi = (isFine) ? 1.0 * M_PI / 180. : 1.25 * M_PI / 180.;
0158         int ir = HGCScintillatorDetId(layerCluster.seed()).iradiusAbs();
0159         float dr = (isFine) ? (0.0484 * static_cast<float>(ir) + 2.1) : (0.075 * static_cast<float>(ir) + 2.0);
0160         float r = (isFine) ? (0.0239 * static_cast<float>(pow(ir, 2)) + 2.02 * static_cast<float>(ir) + 119.6)
0161                            : (0.0367 * static_cast<float>(pow(ir, 2)) + 1.7 * static_cast<float>(ir) + 90.7);
0162         area = r * dr * dphi;
0163       } else {
0164         const bool isFine = (HGCSiliconDetId(layerCluster.seed()).type() == 0);
0165         float side = (isFine) ? 0.465 : 0.698;
0166         area = pow(side, 2) * 3 * sqrt(3) / 2;
0167       }
0168       radius = sqrt(nHits * area) / M_PI;
0169     }
0170 
0171     auto *eveCircle = new TEveGeoShape("Circle");
0172     auto tube = new TGeoTube(0., proportionalityFactor_ * radius, 0.1);
0173     eveCircle->SetShape(tube);
0174     eveCircle->InitMainTrans();
0175     eveCircle->RefMainTrans().Move3PF(position.x(), position.y(), position.z());
0176     setupAddElement(eveCircle, &oItemHolder);
0177     // Apply heatmap color coding **after** the call to setupAddElement, that will internally setup the color.
0178     if (heatmap_) {
0179       const float normalized_energy = fmin(energy / saturation_energy_, 1.0f);
0180       const uint8_t colorFactor = gradient_steps * normalized_energy;
0181       eveCircle->SetFillColor(
0182           TColor::GetColor(gradient[0][colorFactor], gradient[1][colorFactor], gradient[2][colorFactor]));
0183     } else {
0184       eveCircle->SetMainColor(item()->modelInfo(iIndex).displayProperties().color());
0185       eveCircle->SetMainTransparency(item()->defaultDisplayProperties().transparency());
0186     }
0187 
0188     // seed and cluster position
0189     const float crossScale = 1.0f + fmin(energy, 5.0f);
0190     if (enablePositionLines_) {
0191       auto const &pos = layerCluster.position();
0192       const float position_crossScale = crossScale * 0.5;
0193       position_marker->AddLine(
0194           pos.x() - position_crossScale, pos.y(), pos.z(), pos.x() + position_crossScale, pos.y(), pos.z());
0195       position_marker->AddLine(
0196           pos.x(), pos.y() - position_crossScale, pos.z(), pos.x(), pos.y() + position_crossScale, pos.z());
0197     }
0198   }
0199 
0200   if (enablePositionLines_)
0201     oItemHolder.AddElement(position_marker);
0202 
0203   if (enableEdges_) {
0204     auto &edges = trackster.edges();
0205 
0206     TEveStraightLineSet *adjacent_marker = new TEveStraightLineSet;
0207     adjacent_marker->SetLineWidth(2);
0208     adjacent_marker->SetLineColor(kYellow);
0209 
0210     TEveStraightLineSet *non_adjacent_marker = new TEveStraightLineSet;
0211     non_adjacent_marker->SetLineWidth(2);
0212     non_adjacent_marker->SetLineColor(kRed);
0213 
0214     for (auto edge : edges) {
0215       auto doublet = std::make_pair(layerClusters[edge[0]], layerClusters[edge[1]]);
0216 
0217       int layerIn = item()->getGeom()->getParameters(doublet.first.seed())[1];
0218       int layerOut = item()->getGeom()->getParameters(doublet.second.seed())[1];
0219 
0220       const bool isAdjacent = std::abs(layerOut - layerIn) == 1;
0221 
0222       // draw 3D cross
0223       if (layer_ == 0 || fabs(layerIn - layer_) == 0 || fabs(layerOut - layer_) == 0) {
0224         if (isAdjacent)
0225           adjacent_marker->AddLine(doublet.first.x(),
0226                                    doublet.first.y(),
0227                                    doublet.first.z(),
0228                                    doublet.second.x(),
0229                                    doublet.second.y(),
0230                                    doublet.second.z());
0231         else
0232           non_adjacent_marker->AddLine(doublet.first.x(),
0233                                        doublet.first.y(),
0234                                        doublet.first.z(),
0235                                        doublet.second.x(),
0236                                        doublet.second.y(),
0237                                        doublet.second.z());
0238       }
0239     }
0240     oItemHolder.AddElement(adjacent_marker);
0241     oItemHolder.AddElement(non_adjacent_marker);
0242   }
0243 }
0244 
0245 REGISTER_FWPROXYBUILDER(FWTracksterLayersProxyBuilder, ticl::Trackster, "Trackster layers", FWViewType::kISpyBit);