FWECALDetailViewBuilder.cc

Go to the documentation of this file.

// FIXME - needed to set fixed eta-phi limits. Without the
//         visible area may change widely depending on energy
//         deposition availability

#include "TEveCaloData.h"
#include "TEveViewer.h"
#include "TEveCalo.h"
#include "TAxis.h"
#include "THLimitsFinder.h"
#include "TLatex.h"

#include "Fireworks/Core/interface/FWDetailViewBase.h"
#include "Fireworks/Calo/interface/FWECALDetailViewBuilder.h"
#include "Fireworks/Core/interface/FWGeometry.h"
#include "Fireworks/Core/interface/fw3dlego_xbins.h"
#include "Fireworks/Core/interface/fwLog.h"

#include "DataFormats/FWLite/interface/Handle.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"

#include "TGeoMatrix.h"
#include "TEveTrans.h"

#include <utility>

TEveCaloLego* FWECALDetailViewBuilder::build()
{
   // get the hits from the event

   edm::Handle<EcalRecHitCollection> handle_hits;
   const EcalRecHitCollection *hits = 0;

   if (fabs(m_eta) < 1.5)
   {
      try
      {
         edm::InputTag tag("ecalRecHit", "EcalRecHitsEB");
         m_event->getByLabel(tag, handle_hits);
     if (handle_hits.isValid())
        hits = &*handle_hits;
      }
      catch (...)
      {
         fwLog(fwlog::kWarning) <<"no barrel ECAL rechits are available, "
            "showing crystal location but not energy" << std::endl;
      }
   }
   else
   {
      try
      {
         edm::InputTag tag("ecalRecHit", "EcalRecHitsEE");
         m_event->getByLabel(tag, handle_hits);
     if (handle_hits.isValid())
        hits = &*handle_hits;
      }
      catch (...)
      {
         fwLog(fwlog::kWarning) <<"no endcap ECAL rechits are available, "
            "showing crystal location but not energy" << std::endl;
      }
   }
     
   // data
   TEveCaloDataVec* data = new TEveCaloDataVec( 1 + m_colors.size() );
   data->RefSliceInfo(0).Setup("hits (not clustered)", 0.0, m_defaultColor );
   for( size_t i = 0; i < m_colors.size(); ++i )
   {
      data->RefSliceInfo(i + 1).Setup( "hits (clustered)", 0.0, m_colors[i] );
   }

   if( handle_hits.isValid() ) 
      // fill
      fillData( hits, data );

   // axis
   Double_t etaMin(0), etaMax(0), phiMin(0), phiMax(0);
   // it's hard to define properly visible area in X-Y,
   // so we rely on auto limits
   data->GetEtaLimits(etaMin, etaMax);
   data->GetPhiLimits(phiMin, phiMax);
   Double_t bl, bh, bw;
   Int_t bn, n = 20;
   THLimitsFinder::Optimize(etaMin, etaMax, n, bl, bh, bn, bw);
   data->SetEtaBins( new TAxis(bn, bl, bh));
   THLimitsFinder::Optimize(phiMin, phiMax, n, bl, bh, bn, bw);
   data->SetPhiBins( new TAxis(bn, bl, bh));

   // make tower grid
   std::vector<double> etaBinsWithinLimits;
   etaBinsWithinLimits.push_back(etaMin);
   for (unsigned int i=0; i<83; ++i)
      if ( fw3dlego::xbins[i] > etaMin && fw3dlego::xbins[i] < etaMax )
         etaBinsWithinLimits.push_back(fw3dlego::xbins[i]);
   etaBinsWithinLimits.push_back(etaMax);
   Double_t* eta_bins = new Double_t[etaBinsWithinLimits.size()];
   for (unsigned int i=0; i<etaBinsWithinLimits.size(); ++i)
      eta_bins[i] = etaBinsWithinLimits[i];

   std::vector<double> phiBinsWithinLimits;
   phiBinsWithinLimits.push_back(phiMin);
   for ( double phi = -M_PI; phi < M_PI; phi += M_PI/36 )
      if ( phi > phiMin && phi < phiMax ) // it's stupid, I know, but I'm lazy right now
         phiBinsWithinLimits.push_back(phi);
   phiBinsWithinLimits.push_back(phiMax);
   Double_t* phi_bins = new Double_t[phiBinsWithinLimits.size()];
   for (unsigned int i=0; i<phiBinsWithinLimits.size(); ++i)
      phi_bins[i] = phiBinsWithinLimits[i];
   if (fabs(m_eta) > 1.5) {
      data->GetEtaBins()->SetTitle("X[cm]");
      data->GetPhiBins()->SetTitle("Y[cm]");
      data->GetPhiBins()->SetTitleSize(0.03);
      data->GetEtaBins()->SetTitleSize(0.03);
   } else {
      data->SetEtaBins(new TAxis(etaBinsWithinLimits.size()-1,eta_bins));
      data->SetPhiBins(new TAxis(phiBinsWithinLimits.size()-1,phi_bins));
      data->GetEtaBins()->SetTitleFont(122);
      data->GetEtaBins()->SetTitle("h");
      data->GetPhiBins()->SetTitleFont(122);
      data->GetPhiBins()->SetTitle("f");
      data->GetPhiBins()->SetTitleSize(0.05);
      data->GetEtaBins()->SetTitleSize(0.05);
   }
   delete [] eta_bins;
   delete [] phi_bins;

   // lego
   TEveCaloLego *lego = new TEveCaloLego(data);
   lego->SetDrawNumberCellPixels(100);
   // scale and translate to real world coordinates
   lego->SetEta(etaMin, etaMax);
   lego->SetPhiWithRng((phiMin+phiMax)*0.5, (phiMax-phiMin)*0.5); // phi range = 2* phiOffset
   Double_t legoScale = ((etaMax - etaMin) < (phiMax - phiMin)) ? (etaMax - etaMin) : (phiMax - phiMin);
   lego->InitMainTrans();
   lego->RefMainTrans().SetScale(legoScale, legoScale, legoScale*0.5);
   lego->RefMainTrans().SetPos((etaMax+etaMin)*0.5, (phiMax+phiMin)*0.5, 0);
   lego->SetAutoRebin(kFALSE);
   lego->Set2DMode(TEveCaloLego::kValSizeOutline);
   lego->SetName("ECALDetail Lego");
   return lego;

}

void FWECALDetailViewBuilder::setColor(Color_t color, const std::vector<DetId> &detIds)
{

   m_colors.push_back(color);

   // get the slice for this group of detIds
   // note that the zeroth slice is the default one (all else)
   int slice = m_colors.size();
   // take a note of which slice these detids are going to go into
   for (size_t i = 0; i < detIds.size(); ++i)
      m_detIdsToColor[detIds[i]] = slice;
}

void
FWECALDetailViewBuilder::showSuperCluster( const reco::SuperCluster &cluster, Color_t color )
{
   std::vector<DetId> clusterDetIds;
   const std::vector<std::pair<DetId, float> > &hitsAndFractions = cluster.hitsAndFractions();
   for (size_t j = 0; j < hitsAndFractions.size(); ++j)
   {
      clusterDetIds.push_back(hitsAndFractions[j].first);
   }

   setColor( color, clusterDetIds );
}

void
FWECALDetailViewBuilder::showSuperClusters( Color_t color1, Color_t color2 )
{
   // get the superclusters from the event
   edm::Handle<reco::SuperClusterCollection> collection;

   if( fabs( m_eta ) < 1.5 ) {
      try {
         m_event->getByLabel(edm::InputTag("correctedHybridSuperClusters"), collection);
      }
      catch (...)
      {
         fwLog(fwlog::kWarning) <<"no barrel superclusters are available" << std::endl;
      }
   } else {
      try {
         m_event->getByLabel(edm::InputTag("correctedMulti5x5SuperClustersWithPreshower"), collection);
      }
      catch (...)
      {
         fwLog(fwlog::kWarning) <<"no endcap superclusters are available" << std::endl;
      }
   }
   if( collection.isValid() )
   {
      unsigned int colorIndex = 0;
      // sort clusters in eta so neighboring clusters have distinct colors
      reco::SuperClusterCollection sorted = *collection.product();
      std::sort( sorted.begin(), sorted.end(), superClusterEtaLess );
      for( size_t i = 0; i < sorted.size(); ++i )
      {
     if( !(fabs(sorted[i].eta() - m_eta) < (m_size*0.0172)
           && fabs(sorted[i].phi() - m_phi) < (m_size*0.0172)) )
       continue;

     if( colorIndex %2 == 0 )
        showSuperCluster( sorted[i], color1 );
     else
        showSuperCluster( sorted[i], color2 );
     ++colorIndex;
      }
   }
}

void
FWECALDetailViewBuilder::fillData( const EcalRecHitCollection *hits,
                   TEveCaloDataVec *data )
{
   const float barrelCR = m_size*0.0172; // barrel cell range

   // loop on all the detids
   for( EcalRecHitCollection::const_iterator k = hits->begin(), kEnd = hits->end();
    k != kEnd; ++k )
   {
      // get reco geometry
      double centerEta = 0;
      double centerPhi = 0;
      const float* points = m_geom->getCorners( k->id().rawId());
      if( points != 0 )
      {
     TEveVector v;
     int j = 0;
     for( int i = 0; i < 8; ++i )
     {   
        v += TEveVector( points[j], points[j + 1], points[j + 2] );
        j +=3;
     }
     centerEta = v.Eta();
     centerPhi = v.Phi();
      }
      else
     fwLog( fwlog::kInfo ) << "cannot get geometry for DetId: "<< k->id().rawId() << ". Ignored.\n";
      
      double size = k->energy() / cosh( centerEta );

      // check what slice to put in
      int slice = 0;
      std::map<DetId, int>::const_iterator itr = m_detIdsToColor.find(k->id());
      if (itr != m_detIdsToColor.end()) slice = itr->second;

      // if in the EB
      if( k->id().subdetId() == EcalBarrel )
      {
     // do phi wrapping
     if( centerPhi > m_phi + M_PI) centerPhi -= 2 * M_PI;
     if( centerPhi < m_phi - M_PI) centerPhi += 2 * M_PI;

     // check if the hit is in the window to be drawn
     if( !( fabs( centerEta - m_eta ) < barrelCR
        && fabs( centerPhi - m_phi ) < barrelCR )) continue;

     double minEta(10), maxEta(-10), minPhi(4), maxPhi(-4);
     if( points != 0 )
     {
        // calorimeter crystalls have slightly non-symetrical form in eta-phi projection
        // so if we simply get the largest eta and phi, cells will overlap
        // therefore we get a smaller eta-phi range representing the inner square
        // we also should use only points from the inner face of the crystal, since
        // non-projecting direction of crystals leads to large shift in eta on outter
        // face.
        int j = 0;
        for( unsigned int i = 0; i < 8; ++i )
        {
           TEveVector crystal( points[j], points[j + 1], points[j + 2] );
           j += 3;
           double eta = crystal.Eta();
           double phi = crystal.Phi();
           if ( crystal.Perp() > 135 ) continue;
           if ( minEta - eta > 0.01) minEta = eta;
           if ( eta - minEta > 0 && eta - minEta < 0.01 ) minEta = eta;
           if ( eta - maxEta > 0.01) maxEta = eta;
           if ( maxEta - eta > 0 && maxEta - eta < 0.01 ) maxEta = eta;
           if ( minPhi - phi > 0.01) minPhi = phi;
           if ( phi - minPhi > 0 && phi - minPhi < 0.01 ) minPhi = phi;
           if ( phi - maxPhi > 0.01) maxPhi = phi;
           if ( maxPhi - phi > 0 && maxPhi - phi < 0.01 ) maxPhi = phi;
        }
     }
     else 
     {
        double delta = 0.0172 * 0.5;
        minEta = centerEta - delta;
        maxEta = centerEta + delta;
        minPhi = centerPhi - delta;
        maxPhi = centerPhi + delta;
     }
     if( minPhi >= ( m_phi - barrelCR ) && maxPhi <= ( m_phi + barrelCR ) &&
         minEta >= ( m_eta - barrelCR ) && maxEta <= ( m_eta + barrelCR ))
     {
        data->AddTower( minEta, maxEta, minPhi, maxPhi );
        data->FillSlice( slice, size );
     }
     // otherwise in the EE
      }
      else if( k->id().subdetId() == EcalEndcap )
      {
     // check if the hit is in the window to be drawn
     double crystalSize = m_size * 0.0172;
     if( !( fabs( centerEta - m_eta ) < ( crystalSize )
        && fabs( centerPhi - m_phi ) < ( crystalSize )))
        continue;

     if( points != 0 )
     {
        double minX(9999), maxX(-9999), minY(9999), maxY(-9999);
        int j = 0;
        for( unsigned int i = 0; i < 8; ++i )
        {
           TEveVector crystal( points[j], points[j + 1], points[j + 2] );
           j += 3;
           double x = crystal.fX;
           double y = crystal.fY;
           if( fabs( crystal.fZ ) > 330 ) continue;
           if( minX - x > 0.01 ) minX = x;
           if( x - maxX > 0.01 ) maxX = x;
           if( minY - y > 0.01 ) minY = y;
           if( y - maxY > 0.01 ) maxY = y;
        }
        data->AddTower( minX, maxX, minY, maxY );
     }
     data->FillSlice( slice, size );
      }
   } // end loop on hits

   data->DataChanged();
 }

double
FWECALDetailViewBuilder::makeLegend( double x0, double y0,
                                     Color_t clustered1, Color_t clustered2,
                                     Color_t supercluster
                                     )
{
   Double_t fontsize = 0.07;
   TLatex* latex = new TLatex();
   Double_t x = x0;
   Double_t y = y0;
   Double_t boxH = 0.25*fontsize;
   Double_t yStep = 0.04;

   y -= yStep;
   latex->DrawLatex(x, y, "Energy types:");
   y -= yStep;

   Double_t pos[4];
   pos[0] = x+0.05;
   pos[2] = x+0.20;

   pos[1] = y; pos[3] = pos[1] + boxH;
   FWDetailViewBase::drawCanvasBox(pos, m_defaultColor);
   latex->DrawLatex(x+0.25, y, "unclustered");
   y -= yStep;
   if (clustered1<0) return y;

   pos[1] = y; pos[3] = pos[1] + boxH;
   FWDetailViewBase::drawCanvasBox(pos, clustered1);
   latex->DrawLatex(x+0.25, y, "clustered");
   y -= yStep;
   if (clustered2<0) return y;

   pos[1] = y; pos[3] = pos[1] + boxH;
   FWDetailViewBase::drawCanvasBox(pos, clustered2);
   latex->DrawLatex(x+0.25, y, "clustered");
   y -= yStep;
   if (supercluster<0) return y;

   pos[1] = y; pos[3] = pos[1] + boxH;
   FWDetailViewBase::drawCanvasBox(pos, supercluster);
   latex->DrawLatex(x+0.25, y, "super-cluster");
   y -= yStep;

   return y;
}