Calo/src/FWECALCaloDataDetailViewBuilder.cc

0001 // FIXME - needed to set fixed eta-phi limits. Without the
0002 //         visible area may change widely depending on energy
0003 //         deposition availability
0004
0005 #include "TEveCaloData.h"
0006 #include "TEveViewer.h"
0007 #include "TEveCalo.h"
0008 #include "TAxis.h"
0009 #include "TMath.h"
0010 #include "THLimitsFinder.h"
0011 #include "TLatex.h"
0012
0013 #include "Fireworks/Core/interface/FWDetailViewBase.h"
0014 #include "Fireworks/Calo/interface/FWECALCaloDataDetailViewBuilder.h"
0015 #include "Fireworks/Core/interface/FWGeometry.h"
0016 #include "Fireworks/Core/interface/fw3dlego_xbins.h"
0017 #include "Fireworks/Core/interface/fwLog.h"
0018
0019 #include "DataFormats/FWLite/interface/Handle.h"
0020 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0021 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0022
0023 #include "TGeoMatrix.h"
0024 #include "TEveTrans.h"
0025
0026 #include <utility>
0027
0028 TEveCaloData* FWECALCaloDataDetailViewBuilder::buildCaloData(bool xyEE) {
0029   // get the hits from the event
0030
0031   edm::Handle<EcalRecHitCollection> handle_hits;
0032   const EcalRecHitCollection* hits = nullptr;
0033
0034   if (fabs(m_eta) < 1.5) {
0035     try {
0036       edm::InputTag tag("ecalRecHit", "EcalRecHitsEB");
0037       m_event->getByLabel(tag, handle_hits);
0038       if (handle_hits.isValid()) {
0039         hits = &*handle_hits;
0040       }
0041     } catch (...) {
0042       fwLog(fwlog::kWarning) << "FWECALCaloDataDetailViewBuilder::build():: Failed to access EcalRecHitsEB collection."
0043                              << std::endl;
0044     }
0045     if (!handle_hits.isValid()) {
0046       try {
0047         edm::InputTag tag("reducedEcalRecHitsEB");
0048         m_event->getByLabel(tag, handle_hits);
0049         if (handle_hits.isValid()) {
0050           hits = &*handle_hits;
0051         }
0052
0053       } catch (...) {
0054         fwLog(fwlog::kWarning)
0055             << "FWECALCaloDataDetailViewBuilder::build():: Failed to access reducedEcalRecHitsEB collection."
0056             << std::endl;
0057       }
0058     }
0059   } else {
0060     try {
0061       edm::InputTag tag("ecalRecHit", "EcalRecHitsEE");
0062       m_event->getByLabel(tag, handle_hits);
0063       if (handle_hits.isValid())
0064         hits = &*handle_hits;
0065     } catch (...) {
0066       fwLog(fwlog::kWarning) << "FWECALCaloDataDetailViewBuilder::build():: Failed to access ecalRecHitsEE collection."
0067                              << std::endl;
0068     }
0069
0070     if (!handle_hits.isValid()) {
0071       try {
0072         edm::InputTag tag("reducedEcalRecHitsEE");
0073         m_event->getByLabel(tag, handle_hits);
0074         if (handle_hits.isValid()) {
0075           hits = &*handle_hits;
0076         }
0077
0078       } catch (...) {
0079         fwLog(fwlog::kWarning)
0080             << "FWECALCaloDataDetailViewBuilder::build():: Failed to access reducedEcalRecHitsEE collection."
0081             << std::endl;
0082       }
0083     }
0084   }
0085
0086   // data
0087   TEveCaloDataVec* data = new TEveCaloDataVec(1 + m_colors.size());
0088   data->RefSliceInfo(0).Setup("hits (not clustered)", 0.0, m_defaultColor);
0089   for (size_t i = 0; i < m_colors.size(); ++i) {
0090     data->RefSliceInfo(i + 1).Setup("hits (clustered)", 0.0, m_colors[i]);
0091   }
0092
0093   if (handle_hits.isValid()) {
0094     fillData(hits, data, xyEE);
0095   }
0096
0097   // axis
0098   Double_t etaMin(0), etaMax(0), phiMin(0), phiMax(0);
0099   if (data->Empty()) {
0100     fwLog(fwlog::kWarning) << "FWECALCaloDataDetailViewBuilder::build():: No hits found in "
0101                            << Form("FWECALCaloDataDetailViewBuilder::build():: No hits found in eta[%f] phi[%f] region",
0102                                    m_eta,
0103                                    m_phi)
0104                            << ".\n";
0105
0106     // add dummy background
0107     float x = m_size * TMath::DegToRad();
0108     if (fabs(m_eta) < 1.5 || xyEE == false) {
0109       etaMin = m_eta - x;
0110       etaMax = m_eta + x;
0111       phiMin = m_phi - x;
0112       phiMax = m_phi + x;
0113       data->AddTower(etaMin, etaMax, phiMin, phiMax);
0114     } else {
0115       float theta = TEveCaloData::EtaToTheta(m_eta);
0116       float r = TMath::Tan(theta) * 290;
0117       phiMin = r * TMath::Cos(m_phi - x) - 300;
0118       phiMax = r * TMath::Cos(m_phi + x) + 300;
0119       etaMin = r * TMath::Sin(m_phi - x) - 300;
0120       etaMax = r * TMath::Sin(m_phi + x) + 300;
0121       data->AddTower(TMath::Min(etaMin, etaMax),
0122                      TMath::Max(etaMin, etaMax),
0123                      TMath::Min(phiMin, phiMax),
0124                      TMath::Max(phiMin, phiMax));
0125     }
0126     data->FillSlice(0, 0.1);
0127   }
0128
0129   TAxis* eta_axis = nullptr;
0130   TAxis* phi_axis = nullptr;
0131   data->GetEtaLimits(etaMin, etaMax);
0132   data->GetPhiLimits(phiMin, phiMax);
0133   //  printf("data rng %f %f %f %f\n",etaMin, etaMax, phiMin, phiMax );
0134
0135   if (fabs(m_eta) > 1.5 && xyEE) {
0136     eta_axis = new TAxis(10, etaMin, etaMax);
0137     phi_axis = new TAxis(10, phiMin, phiMax);
0138     eta_axis->SetTitle("X[cm]");
0139     phi_axis->SetTitle("Y[cm]");
0140     phi_axis->SetTitleSize(0.05);
0141     eta_axis->SetTitleSize(0.05);
0142   } else {
0143     std::vector<double> etaBinsWithinLimits;
0144     etaBinsWithinLimits.push_back(etaMin);
0145     for (unsigned int i = 0; i < 83; ++i)
0146       if (fw3dlego::xbins[i] > etaMin && fw3dlego::xbins[i] < etaMax)
0147         etaBinsWithinLimits.push_back(fw3dlego::xbins[i]);
0148     etaBinsWithinLimits.push_back(etaMax);
0149
0150     std::vector<double> phiBinsWithinLimits;
0151     phiBinsWithinLimits.push_back(phiMin);
0152     for (double phi = -M_PI; phi < M_PI; phi += M_PI / 36)
0153       if (phi > phiMin && phi < phiMax)
0154         phiBinsWithinLimits.push_back(phi);
0155     phiBinsWithinLimits.push_back(phiMax);
0156
0157     eta_axis = new TAxis((int)etaBinsWithinLimits.size() - 1, &etaBinsWithinLimits[0]);
0158     phi_axis = new TAxis((int)phiBinsWithinLimits.size() - 1, &phiBinsWithinLimits[0]);
0159
0160     eta_axis->SetTitleFont(122);
0161     eta_axis->SetTitle("h");
0162     eta_axis->SetTitleSize(0.07);
0163     phi_axis->SetTitleFont(122);
0164     phi_axis->SetTitle("f");
0165     phi_axis->SetTitleSize(0.07);
0166   }
0167   eta_axis->SetNdivisions(510);
0168   phi_axis->SetNdivisions(510);
0169   data->SetEtaBins(eta_axis);
0170   data->SetPhiBins(phi_axis);
0171   return data;
0172 }
0173
0174 //_______________________________________________________________
0175 TEveCaloLego* FWECALCaloDataDetailViewBuilder::build() {
0176   TEveCaloData* data = buildCaloData(true);
0177
0178   double etaMin, etaMax, phiMin, phiMax;
0179   data->GetEtaLimits(etaMin, etaMax);
0180   data->GetPhiLimits(phiMin, phiMax);
0181
0182   // lego
0183   TEveCaloLego* lego = new TEveCaloLego(data);
0184   lego->SetDrawNumberCellPixels(100);
0185   // scale and translate to real world coordinates
0186   lego->SetEta(etaMin, etaMax);
0187   lego->SetPhiWithRng((phiMin + phiMax) * 0.5, (phiMax - phiMin) * 0.5);  // phi range = 2* phiOffset
0188   Double_t legoScale = ((etaMax - etaMin) < (phiMax - phiMin)) ? (etaMax - etaMin) : (phiMax - phiMin);
0189   lego->InitMainTrans();
0190   lego->RefMainTrans().SetScale(legoScale, legoScale, legoScale * 0.5);
0191   lego->RefMainTrans().SetPos((etaMax + etaMin) * 0.5, (phiMax + phiMin) * 0.5, 0);
0192   lego->SetAutoRebin(kFALSE);
0193   lego->Set2DMode(TEveCaloLego::kValSizeOutline);
0194   lego->SetName("ECALDetail Lego");
0195   return lego;
0196 }
0197
0198 void FWECALCaloDataDetailViewBuilder::setColor(Color_t color, const std::vector<DetId>& detIds) {
0199   m_colors.push_back(color);
0200
0201   // get the slice for this group of detIds
0202   // note that the zeroth slice is the default one (all else)
0203   int slice = m_colors.size();
0204   // take a note of which slice these detids are going to go into
0205   for (size_t i = 0; i < detIds.size(); ++i)
0206     m_detIdsToColor[detIds[i]] = slice;
0207 }
0208
0209 void FWECALCaloDataDetailViewBuilder::showSuperCluster(const reco::SuperCluster& cluster, Color_t color) {
0210   std::vector<DetId> clusterDetIds;
0211   const std::vector<std::pair<DetId, float> >& hitsAndFractions = cluster.hitsAndFractions();
0212   for (size_t j = 0; j < hitsAndFractions.size(); ++j) {
0213     clusterDetIds.push_back(hitsAndFractions[j].first);
0214   }
0215
0216   setColor(color, clusterDetIds);
0217 }
0218
0219 void FWECALCaloDataDetailViewBuilder::showSuperClusters(Color_t color1, Color_t color2) {
0220   // get the superclusters from the event
0221   edm::Handle<reco::SuperClusterCollection> collection;
0222
0223   if (fabs(m_eta) < 1.5) {
0224     try {
0225       m_event->getByLabel(edm::InputTag("correctedHybridSuperClusters"), collection);
0226     } catch (...) {
0227       fwLog(fwlog::kWarning) << "no barrel superclusters are available" << std::endl;
0228     }
0229   } else {
0230     try {
0231       m_event->getByLabel(edm::InputTag("correctedMulti5x5SuperClustersWithPreshower"), collection);
0232     } catch (...) {
0233       fwLog(fwlog::kWarning) << "no endcap superclusters are available" << std::endl;
0234     }
0235   }
0236   if (collection.isValid()) {
0237     unsigned int colorIndex = 0;
0238     // sort clusters in eta so neighboring clusters have distinct colors
0239     reco::SuperClusterCollection sorted = *collection.product();
0240     std::sort(sorted.begin(), sorted.end(), superClusterEtaLess);
0241     for (size_t i = 0; i < sorted.size(); ++i) {
0242       if (!(fabs(sorted[i].eta() - m_eta) < (m_size * 0.0172) && fabs(sorted[i].phi() - m_phi) < (m_size * 0.0172)))
0243         continue;
0244
0245       if (colorIndex % 2 == 0)
0246         showSuperCluster(sorted[i], color1);
0247       else
0248         showSuperCluster(sorted[i], color2);
0249       ++colorIndex;
0250     }
0251   }
0252 }
0253
0254 void FWECALCaloDataDetailViewBuilder::fillData(const EcalRecHitCollection* hits, TEveCaloDataVec* data, bool xyEE) {
0255   const float barrelCR = m_size * 0.0172;  // barrel cell range
0256
0257   // loop on all the detids
0258   for (EcalRecHitCollection::const_iterator k = hits->begin(), kEnd = hits->end(); k != kEnd; ++k) {
0259     // get reco geometry
0260     double centerEta = 0;
0261     double centerPhi = 0;
0262     const float* points = m_geom->getCorners(k->id().rawId());
0263     if (points != nullptr) {
0264       TEveVector v;
0265       int j = 0;
0266       for (int i = 0; i < 8; ++i) {
0267         v += TEveVector(points[j], points[j + 1], points[j + 2]);
0268         j += 3;
0269       }
0270       centerEta = v.Eta();
0271       centerPhi = v.Phi();
0272     } else
0273       fwLog(fwlog::kInfo) << "cannot get geometry for DetId: " << k->id().rawId() << ". Ignored.\n";
0274
0275     double size = k->energy() / cosh(centerEta);
0276
0277     // check what slice to put in
0278     int slice = 0;
0279     std::map<DetId, int>::const_iterator itr = m_detIdsToColor.find(k->id());
0280     if (itr != m_detIdsToColor.end())
0281       slice = itr->second;
0282
0283     // if in the EB
0284     if (k->id().subdetId() == EcalBarrel || xyEE == false) {
0285       // do phi wrapping
0286       if (centerPhi > m_phi + M_PI)
0287         centerPhi -= 2 * M_PI;
0288       if (centerPhi < m_phi - M_PI)
0289         centerPhi += 2 * M_PI;
0290
0291       // check if the hit is in the window to be drawn
0292       if (!(fabs(centerEta - m_eta) < barrelCR && fabs(centerPhi - m_phi) < barrelCR))
0293         continue;
0294
0295       double minEta(10), maxEta(-10), minPhi(4), maxPhi(-4);
0296       if (points != nullptr) {
0297         // calorimeter crystalls have slightly non-symetrical form in eta-phi projection
0298         // so if we simply get the largest eta and phi, cells will overlap
0299         // therefore we get a smaller eta-phi range representing the inner square
0300         // we also should use only points from the inner face of the crystal, since
0301         // non-projecting direction of crystals leads to large shift in eta on outter
0302         // face.
0303         int j = 0;
0304         float eps = 0.005;
0305         for (unsigned int i = 0; i < 8; ++i) {
0306           TEveVector crystal(points[j], points[j + 1], points[j + 2]);
0307           j += 3;
0308           double eta = crystal.Eta();
0309           double phi = crystal.Phi();
0310           if (((k->id().subdetId() == EcalBarrel) && (crystal.Perp() > 135)) ||
0311               ((k->id().subdetId() == EcalEndcap) && (crystal.Perp() > 155)))
0312             continue;
0313           if (minEta - eta > eps)
0314             minEta = eta;
0315           if (eta - minEta > 0 && eta - minEta < eps)
0316             minEta = eta;
0317           if (eta - maxEta > eps)
0318             maxEta = eta;
0319           if (maxEta - eta > 0 && maxEta - eta < eps)
0320             maxEta = eta;
0321           if (minPhi - phi > eps)
0322             minPhi = phi;
0323           if (phi - minPhi > 0 && phi - minPhi < eps)
0324             minPhi = phi;
0325           if (phi - maxPhi > eps)
0326             maxPhi = phi;
0327           if (maxPhi - phi > 0 && maxPhi - phi < eps)
0328             maxPhi = phi;
0329         }
0330       } else {
0331         double delta = 0.0172 * 0.5;
0332         minEta = centerEta - delta;
0333         maxEta = centerEta + delta;
0334         minPhi = centerPhi - delta;
0335         maxPhi = centerPhi + delta;
0336       }
0337       if (minPhi >= (m_phi - barrelCR) && maxPhi <= (m_phi + barrelCR) && minEta >= (m_eta - barrelCR) &&
0338           maxEta <= (m_eta + barrelCR)) {
0339         // printf("add %f %f %f %f \n",minEta, maxEta, minPhi, maxPhi );
0340         data->AddTower(minEta, maxEta, minPhi, maxPhi);
0341         data->FillSlice(slice, size);
0342       }
0343       // otherwise in the EE
0344     } else if (k->id().subdetId() == EcalEndcap) {
0345       // check if the hit is in the window to be drawn
0346       double crystalSize = m_size * 0.0172;
0347       if (!(fabs(centerEta - m_eta) < (crystalSize) && fabs(centerPhi - m_phi) < (crystalSize)))
0348         continue;
0349
0350       if (points != nullptr) {
0351         double minX(9999), maxX(-9999), minY(9999), maxY(-9999);
0352         int j = 0;
0353         for (unsigned int i = 0; i < 8; ++i) {
0354           TEveVector crystal(points[j], points[j + 1], points[j + 2]);
0355           j += 3;
0356           double x = crystal.fX;
0357           double y = crystal.fY;
0358           if (fabs(crystal.fZ) > 330)
0359             continue;
0360           if (minX - x > 0.01)
0361             minX = x;
0362           if (x - maxX > 0.01)
0363             maxX = x;
0364           if (minY - y > 0.01)
0365             minY = y;
0366           if (y - maxY > 0.01)
0367             maxY = y;
0368         }
0369         data->AddTower(minX, maxX, minY, maxY);
0370         // printf("EE add %f %f %f %f \n",minX, maxX, minY, maxY );
0371       }
0372       data->FillSlice(slice, size);
0373     }
0374   }  // end loop on hits
0375
0376   data->DataChanged();
0377 }
0378
0379 double FWECALCaloDataDetailViewBuilder::makeLegend(
0380     double x0, double y0, Color_t clustered1, Color_t clustered2, Color_t supercluster) {
0381   Double_t fontsize = 0.07;
0382   TLatex* latex = new TLatex();
0383   Double_t x = x0;
0384   Double_t y = y0;
0385   Double_t boxH = 0.25 * fontsize;
0386   Double_t yStep = 0.04;
0387
0388   y -= yStep;
0389   latex->DrawLatex(x, y, "Energy types:");
0390   y -= yStep;
0391
0392   Double_t pos[4];
0393   pos[0] = x + 0.05;
0394   pos[2] = x + 0.20;
0395
0396   pos[1] = y;
0397   pos[3] = pos[1] + boxH;
0398   FWDetailViewBase::drawCanvasBox(pos, m_defaultColor);
0399   latex->DrawLatex(x + 0.25, y, "unclustered");
0400   y -= yStep;
0401   if (clustered1 < 0)
0402     return y;
0403
0404   pos[1] = y;
0405   pos[3] = pos[1] + boxH;
0406   FWDetailViewBase::drawCanvasBox(pos, clustered1);
0407   latex->DrawLatex(x + 0.25, y, "clustered");
0408   y -= yStep;
0409   if (clustered2 < 0)
0410     return y;
0411
0412   pos[1] = y;
0413   pos[3] = pos[1] + boxH;
0414   FWDetailViewBase::drawCanvasBox(pos, clustered2);
0415   latex->DrawLatex(x + 0.25, y, "clustered");
0416   y -= yStep;
0417   if (supercluster < 0)
0418     return y;
0419
0420   pos[1] = y;
0421   pos[3] = pos[1] + boxH;
0422   FWDetailViewBase::drawCanvasBox(pos, supercluster);
0423   latex->DrawLatex(x + 0.25, y, "super-cluster");
0424   y -= yStep;
0425
0426   return y;
0427 }