Calo/src/FWECALDetailViewBuilder.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 "TEvePointSet.h"
0008 #include "TEveCalo.h"
0009 #include "TEveCompound.h"
0010 #include "TAxis.h"
0011 #include "TMath.h"
0012 #include "THLimitsFinder.h"
0013 #include "TLatex.h"
0014
0015 #include "Fireworks/Core/interface/FWDetailViewBase.h"
0016 #include "Fireworks/Calo/interface/FWECALDetailViewBuilder.h"
0017 #include "Fireworks/Calo/interface/FWBoxRecHit.h"
0018 #include "Fireworks/Core/interface/FWGeometry.h"
0019 #include "Fireworks/Core/interface/fw3dlego_xbins.h"
0020 #include "Fireworks/Core/interface/fwLog.h"
0021 #include "Fireworks/Core/interface/Context.h"
0022
0023 #include "DataFormats/FWLite/interface/Handle.h"
0024 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0025 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0026
0027 #include "TGeoMatrix.h"
0028 #include "TEveTrans.h"
0029
0030 #include <utility>
0031
0032 FWECALDetailViewBuilder::FWECALDetailViewBuilder(
0033     const edm::EventBase *event, const FWGeometry *geom, float eta, float phi, int size, Color_t defaultColor)
0034
0035     : m_event(event),
0036       m_geom(geom),
0037       m_eta(eta),
0038       m_phi(phi),
0039       m_size(size),
0040       m_defaultColor(defaultColor),
0041       m_towerList(nullptr) {}
0042
0043 TEveCaloData *FWECALDetailViewBuilder::buildCaloData(bool) {
0044   // get the hits from the event
0045
0046   // data
0047   TEveCaloDataVec *data = new TEveCaloDataVec(1);
0048   data->SetWrapTwoPi(false);
0049   data->RefSliceInfo(0).Setup("hits (not clustered)", 0.0, m_defaultColor);
0050
0051   fillData(data);
0052
0053   // axis
0054   float etaMin = m_eta - sizeRad();
0055   float etaMax = m_eta + sizeRad();
0056   float phiMin = m_phi - sizeRad();
0057   float phiMax = m_phi + sizeRad();
0058
0059   data->AddTower(m_eta - sizeRad(), m_eta + sizeRad(), m_phi - sizeRad(), m_phi + sizeRad());
0060
0061   data->FillSlice(0, 0.1);
0062
0063   TAxis *eta_axis = nullptr;
0064   TAxis *phi_axis = nullptr;
0065
0066   //  printf("data rng %f %f %f %f\n",etaMin, etaMax, phiMin, phiMax );
0067   std::vector<double> etaBinsWithinLimits;
0068   etaBinsWithinLimits.push_back(etaMin);
0069   for (unsigned int i = 0; i < 83; ++i)
0070     if (fw3dlego::xbins[i] > etaMin && fw3dlego::xbins[i] < etaMax)
0071       etaBinsWithinLimits.push_back(fw3dlego::xbins[i]);
0072   etaBinsWithinLimits.push_back(etaMax);
0073
0074   std::vector<double> phiBinsWithinLimits;
0075   phiBinsWithinLimits.push_back(phiMin);
0076   for (double phi = -M_PI; phi < M_PI; phi += M_PI / 36)
0077     if (phi > phiMin && phi < phiMax)
0078       phiBinsWithinLimits.push_back(phi);
0079   phiBinsWithinLimits.push_back(phiMax);
0080
0081   eta_axis = new TAxis((int)etaBinsWithinLimits.size() - 1, &etaBinsWithinLimits[0]);
0082   phi_axis = new TAxis((int)phiBinsWithinLimits.size() - 1, &phiBinsWithinLimits[0]);
0083
0084   eta_axis->SetTitleFont(122);
0085   eta_axis->SetTitle("h");
0086   eta_axis->SetTitleSize(0.07);
0087   phi_axis->SetTitleFont(122);
0088   phi_axis->SetTitle("f");
0089   phi_axis->SetTitleSize(0.07);
0090
0091   eta_axis->SetNdivisions(510);
0092   phi_axis->SetNdivisions(510);
0093   data->SetEtaBins(eta_axis);
0094   data->SetPhiBins(phi_axis);
0095   return data;
0096 }
0097
0098 //_______________________________________________________________
0099 TEveCaloLego *FWECALDetailViewBuilder::build() {
0100   // axis
0101   float etaMin = m_eta - sizeRad();
0102   float etaMax = m_eta + sizeRad();
0103   float phiMin = m_phi - sizeRad();
0104   float phiMax = m_phi + sizeRad();
0105
0106   m_towerList = new TEveElementList("TowerHolder");
0107   TEveCaloData *data = buildCaloData(true);
0108
0109   // lego
0110   TEveCaloLego *lego = new TEveCaloLego();
0111   lego->SetData(data);
0112   lego->AddElement(m_towerList);
0113   lego->SetAutoRange(false);
0114   lego->SetDrawNumberCellPixels(100);
0115   // scale and translate to real world coordinates
0116   lego->SetEta(etaMin, etaMax);
0117   lego->SetPhiWithRng((phiMin + phiMax) * 0.5, (phiMax - phiMin) * 0.5);  // phi range = 2* phiOffset
0118   Double_t legoScale = sizeRad() * 2;
0119   lego->InitMainTrans();
0120   lego->RefMainTrans().SetScale(legoScale, legoScale, legoScale * 0.5);
0121   lego->RefMainTrans().SetPos(m_eta, m_phi, -0.01);
0122   lego->SetAutoRebin(kFALSE);
0123   lego->SetName("ECALDetail Lego");
0124
0125   // cut & paste from FWLegoViewBase
0126   lego->SetScaleAbs(true);
0127   lego->SetHasFixedHeightIn2DMode(true);
0128   lego->SetFixedHeightValIn2DMode(0.001);
0129
0130   TEvePointSet *ps = new TEvePointSet("origin");
0131   ps->SetNextPoint(m_eta, m_phi, 0.01);
0132   ps->SetMarkerSize(0.05);
0133   ps->SetMarkerStyle(2);
0134   ps->SetMainColor(kGreen);
0135   ps->SetMarkerColor(kGreen);
0136   lego->AddElement(ps);
0137
0138   return lego;
0139 }
0140
0141 void FWECALDetailViewBuilder::setColor(Color_t color, const std::vector<DetId> &detIds) {
0142   for (size_t i = 0; i < detIds.size(); ++i)
0143     m_detIdsToColor[detIds[i]] = color;
0144 }
0145
0146 void FWECALDetailViewBuilder::showSuperCluster(const reco::SuperCluster &cluster, Color_t color) {
0147   std::vector<DetId> clusterDetIds;
0148   const std::vector<std::pair<DetId, float> > &hitsAndFractions = cluster.hitsAndFractions();
0149   for (size_t j = 0; j < hitsAndFractions.size(); ++j) {
0150     clusterDetIds.push_back(hitsAndFractions[j].first);
0151   }
0152
0153   setColor(color, clusterDetIds);
0154 }
0155
0156 void FWECALDetailViewBuilder::showSuperClusters(Color_t color1, Color_t color2) {
0157   // get the superclusters from the event
0158   edm::Handle<reco::SuperClusterCollection> collection;
0159
0160   if (fabs(m_eta) < 1.5) {
0161     try {
0162       m_event->getByLabel(edm::InputTag("correctedHybridSuperClusters"), collection);
0163     } catch (...) {
0164       fwLog(fwlog::kWarning) << "no barrel superclusters are available" << std::endl;
0165     }
0166   } else {
0167     try {
0168       m_event->getByLabel(edm::InputTag("correctedMulti5x5SuperClustersWithPreshower"), collection);
0169     } catch (...) {
0170       fwLog(fwlog::kWarning) << "no endcap superclusters are available" << std::endl;
0171     }
0172   }
0173   if (collection.isValid()) {
0174     unsigned int colorIndex = 0;
0175     // sort clusters in eta so neighboring clusters have distinct colors
0176     reco::SuperClusterCollection sorted = *collection.product();
0177     std::sort(sorted.begin(), sorted.end(), superClusterEtaLess);
0178     for (size_t i = 0; i < sorted.size(); ++i) {
0179       if (!(fabs(sorted[i].eta() - m_eta) < sizeRad() && fabs(sorted[i].phi() - m_phi) < sizeRad()))
0180         continue;
0181
0182       if (colorIndex % 2 == 0)
0183         showSuperCluster(sorted[i], color1);
0184       else
0185         showSuperCluster(sorted[i], color2);
0186       ++colorIndex;
0187     }
0188   }
0189 }
0190
0191 namespace {
0192   float calculateEt(const TEveVector &centre, float e) {
0193     TEveVector vec = centre;
0194     float et;
0195
0196     vec.Normalize();
0197     vec *= e;
0198     et = vec.Perp();
0199
0200     return et;
0201   }
0202
0203 }  // namespace
0204 //------------------------------------------------------------------
0205 void FWECALDetailViewBuilder::fillEtaPhi(const EcalRecHitCollection *hits, TEveCaloDataVec *data) {
0206   // printf("filletaphi \n");
0207   const float area = sizeRad();  // barrel cell range, AMT this is available in context
0208
0209   double eta1 = m_eta - area;
0210   double eta2 = m_eta + area;
0211   double phi1 = m_phi - area;
0212   double phi2 = m_phi + area;
0213
0214   std::vector<FWBoxRecHit *> boxes;
0215   for (EcalRecHitCollection::const_iterator hitIt = hits->begin(); hitIt != hits->end(); ++hitIt) {
0216     const float *corners = m_geom->getCorners(hitIt->detid());
0217     float energy, et;
0218     std::vector<TEveVector> etaphiCorners(8);
0219
0220     if (corners == nullptr)
0221       continue;
0222
0223     for (int i = 0; i < 4; ++i) {
0224       TEveVector cv = TEveVector(corners[i * 3], corners[i * 3 + 1], corners[i * 3 + 2]);
0225       etaphiCorners[i].fX = cv.Eta();  // Conversion of rechit X/Y values for plotting in Eta/Phi
0226       etaphiCorners[i].fY = cv.Phi();
0227       etaphiCorners[i].fZ = 0.0;
0228
0229       etaphiCorners[i + 4].fX =
0230           etaphiCorners[i].fX;  // Top can simply be plotted exactly over the top of the bottom face
0231       etaphiCorners[i + 4].fY = etaphiCorners[i].fY;
0232       etaphiCorners[i + 4].fZ = 0.001;
0233       //  printf("%f %f %d \n",  etaphiCorners[i].fX, etaphiCorners[i].fY, i);
0234     }
0235
0236     TEveVector center;
0237     for (int i = 0; i < 4; ++i)
0238       center += etaphiCorners[i];
0239     center *= 1.f / 4.f;
0240
0241     if (center.fX < eta1 || center.fX > eta2)
0242       continue;
0243     if (center.fY < phi1 || center.fY > phi2)
0244       continue;
0245
0246     // Stop phi wrap
0247     float dPhi1 = etaphiCorners[2].fY - etaphiCorners[1].fY;
0248     float dPhi2 = etaphiCorners[3].fY - etaphiCorners[0].fY;
0249     float dPhi3 = etaphiCorners[1].fY - etaphiCorners[2].fY;
0250     float dPhi4 = etaphiCorners[0].fY - etaphiCorners[3].fY;
0251
0252     if (dPhi1 > 1)
0253       etaphiCorners[2].fY = etaphiCorners[2].fY - (2 * TMath::Pi());
0254     if (dPhi2 > 1)
0255       etaphiCorners[3].fY = etaphiCorners[3].fY - (2 * TMath::Pi());
0256     if (dPhi3 > 1)
0257       etaphiCorners[2].fY = etaphiCorners[2].fY + (2 * TMath::Pi());
0258     if (dPhi4 > 1)
0259       etaphiCorners[3].fY = etaphiCorners[3].fY + (2 * TMath::Pi());
0260
0261     energy = hitIt->energy();
0262     et = calculateEt(center, energy);
0263     Color_t bcolor = m_defaultColor;
0264     std::map<DetId, int>::const_iterator itr = m_detIdsToColor.find(hitIt->id());
0265     if (itr != m_detIdsToColor.end())
0266       bcolor = itr->second;
0267
0268     m_boxes.push_back(new FWBoxRecHit(etaphiCorners, m_towerList, energy, et));
0269     TEveElement::List_i pIt = m_boxes.back()->getTower()->BeginParents();
0270     TEveCompound *comp = dynamic_cast<TEveCompound *>(*pIt);
0271     comp->SetMainColor(bcolor);
0272     m_boxes.back()->getTower()->SetPickable(true);
0273     m_boxes.back()->getTower()->SetElementTitle(Form("rawId = %d, et = %f", hitIt->id().rawId(), et));
0274   }  // loop hits
0275 }
0276
0277 //---------------------------------------------------------------------------------------
0278
0279 void FWECALDetailViewBuilder::fillData(TEveCaloDataVec *data) {
0280   {  // barrel
0281     const EcalRecHitCollection *hitsEB = nullptr;
0282     edm::Handle<EcalRecHitCollection> handle_hitsEB;
0283
0284     // RECO
0285     try {
0286       edm::InputTag tag("ecalRecHit", "EcalRecHitsEB");
0287       m_event->getByLabel(tag, handle_hitsEB);
0288       if (handle_hitsEB.isValid()) {
0289         hitsEB = &*handle_hitsEB;
0290       }
0291     } catch (...) {
0292       fwLog(fwlog::kWarning) << "FWECALDetailViewBuilder::fillData():: Failed to access EcalRecHitsEB collection."
0293                              << std::endl;
0294     }
0295
0296     // AOD
0297     if (!handle_hitsEB.isValid()) {
0298       try {
0299         edm::InputTag tag("reducedEcalRecHitsEB");
0300         m_event->getByLabel(tag, handle_hitsEB);
0301         if (handle_hitsEB.isValid()) {
0302           hitsEB = &*handle_hitsEB;
0303         }
0304
0305       } catch (...) {
0306         fwLog(fwlog::kWarning)
0307             << "FWECALDetailViewBuilder::filData():: Failed to access reducedEcalRecHitsEB collection." << std::endl;
0308       }
0309     }
0310
0311     // MINIAOD
0312     if (!handle_hitsEB.isValid()) {
0313       try {
0314         edm::InputTag tag("reducedEgamma", "reducedEBRecHits");
0315         m_event->getByLabel(tag, handle_hitsEB);
0316         if (handle_hitsEB.isValid()) {
0317           hitsEB = &*handle_hitsEB;
0318         }
0319
0320       } catch (...) {
0321         fwLog(fwlog::kWarning) << "FWECALDetailViewBuilder::filData():: Failed to access reducedEgamma collection."
0322                                << std::endl;
0323       }
0324     }
0325
0326     if (handle_hitsEB.isValid()) {
0327       fillEtaPhi(hitsEB, data);
0328     }
0329   }
0330
0331   {  // endcap
0332
0333     const EcalRecHitCollection *hitsEE = nullptr;
0334     edm::Handle<EcalRecHitCollection> handle_hitsEE;
0335
0336     // RECO
0337     try {
0338       edm::InputTag tag("ecalRecHit", "EcalRecHitsEE");
0339       m_event->getByLabel(tag, handle_hitsEE);
0340       if (handle_hitsEE.isValid())
0341         hitsEE = &*handle_hitsEE;
0342     } catch (...) {
0343       fwLog(fwlog::kWarning) << "FWECALDetailViewBuilder::fillData():: Failed to access ecalRecHitsEE collection."
0344                              << std::endl;
0345     }
0346
0347     // AOD
0348     if (!handle_hitsEE.isValid()) {
0349       try {
0350         edm::InputTag tag("reducedEcalRecHitsEE");
0351         m_event->getByLabel(tag, handle_hitsEE);
0352         if (handle_hitsEE.isValid()) {
0353           hitsEE = &*handle_hitsEE;
0354         }
0355
0356       } catch (...) {
0357         fwLog(fwlog::kWarning)
0358             << "FWECALDetailViewBuilder::fillData():: Failed to access reducedEcalRecHitsEE collection." << std::endl;
0359       }
0360
0361       // MINIAOD
0362       if (!handle_hitsEE.isValid()) {
0363         try {
0364           edm::InputTag tag("reducedEgamma", "reducedEERecHits");
0365           m_event->getByLabel(tag, handle_hitsEE);
0366           if (handle_hitsEE.isValid()) {
0367             hitsEE = &*handle_hitsEE;
0368           }
0369
0370         } catch (...) {
0371           fwLog(fwlog::kWarning)
0372               << "FWECALDetailViewBuilder::fillData():: Failed to access reducedEcalRecHitsEE collection." << std::endl;
0373         }
0374       }
0375     }
0376
0377     if (handle_hitsEE.isValid()) {
0378       fillEtaPhi(hitsEE, data);
0379     }
0380   }
0381
0382   if (m_boxes.empty())
0383     return;
0384
0385   bool plotEt = true;
0386   float maxEnergy = 0;
0387   int maxEnergyIdx = 0;
0388   // get max energy in EE and EB
0389
0390   int cnt = 0;
0391   for (auto &i : m_boxes) {
0392     if (i->getEnergy(plotEt) > maxEnergy) {
0393       maxEnergy = i->getEnergy(plotEt);
0394       maxEnergyIdx = cnt;
0395     }
0396     cnt++;
0397   }
0398
0399   m_boxes[maxEnergyIdx]->setIsTallest();
0400
0401   // AMT ... max size can be an external parameter
0402   float scale = 0.3 / maxEnergy;
0403   for (auto &i : m_boxes) {
0404     i->updateScale(scale, log(maxEnergy + 1), plotEt);
0405     i->getTower()->SetDrawFrame(true);
0406   }
0407   data->DataChanged();
0408 }
0409
0410 double FWECALDetailViewBuilder::makeLegend(
0411     double x0, double y0, Color_t clustered1, Color_t clustered2, Color_t supercluster) {
0412   Double_t fontsize = 0.07;
0413   TLatex *latex = new TLatex();
0414   Double_t x = x0;
0415   Double_t y = y0;
0416   Double_t boxH = 0.25 * fontsize;
0417   Double_t yStep = 0.04;
0418
0419   y -= yStep;
0420   latex->DrawLatex(x, y, "Energy types:");
0421   y -= yStep;
0422
0423   Double_t pos[4];
0424   pos[0] = x + 0.05;
0425   pos[2] = x + 0.20;
0426
0427   pos[1] = y;
0428   pos[3] = pos[1] + boxH;
0429   FWDetailViewBase::drawCanvasBox(pos, m_defaultColor);
0430   latex->DrawLatex(x + 0.25, y, "unclustered");
0431   y -= yStep;
0432   if (clustered1 < 0)
0433     return y;
0434
0435   pos[1] = y;
0436   pos[3] = pos[1] + boxH;
0437   FWDetailViewBase::drawCanvasBox(pos, clustered1);
0438   latex->DrawLatex(x + 0.25, y, "clustered");
0439   y -= yStep;
0440   if (clustered2 < 0)
0441     return y;
0442
0443   pos[1] = y;
0444   pos[3] = pos[1] + boxH;
0445   FWDetailViewBase::drawCanvasBox(pos, clustered2);
0446   latex->DrawLatex(x + 0.25, y, "clustered");
0447   y -= yStep;
0448   if (supercluster < 0)
0449     return y;
0450
0451   pos[1] = y;
0452   pos[3] = pos[1] + boxH;
0453   FWDetailViewBase::drawCanvasBox(pos, supercluster);
0454   latex->DrawLatex(x + 0.25, y, "super-cluster");
0455   y -= yStep;
0456
0457   return y;
0458 }
0459 //______________________________________________________________________________
0460
0461 float FWECALDetailViewBuilder::sizeRad() const {
0462   float rs = m_size * TMath::DegToRad();
0463   return rs;
0464 }