Display3D.cxx

Go to the documentation of this file.

#include "EventDisplay/Display3D.h"
#include <vector>
#include "TCanvas.h"
#include "TView3D.h"
#include "TGLViewer.h"
#include "TMarker3DBox.h"
#include "TPolyLine3D.h"
#include "TPolyMarker3D.h"
#include "TView.h"
#include "TH1F.h"
#include "EventDisplayBase/IoModule.h"
#include "EventDisplayBase/View3D.h"
#include "EventDataModel/edm.h"
#include "Geometry/geo.h"
#include "RawData/RawDigit.h"
#include "RawData/DAQHeader.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/FLSHit.h"
#include "EventDisplay/Style.h"
#include "EventDisplay/ColorScale.h"
using namespace evd;

Display3D::Display3D(TGMainFrame* mf) : 
  evdb::Canvas(mf),
  fQscale(new ColorScale(0.0, 512.0, kRainbow, kSqrt)),
  fPad3D(0),
  fTimePad(0),
  fGLViewer3D(0)
{
  evdb::Canvas::fCanvas->cd();

  fTimePad = new TPad("fTimePad","Hits vs. Time",0.7,0.7,1.0,1.0,0);
  fTimePad->Draw();
  
  fChargePad = new TPad("fChargePad","Hits vs. Charge",0.7,0.4,1.0,0.7);
  fChargePad->SetLogx();
  fChargePad->SetLogy();
  fChargePad->Draw();
  
  fPad3D = new TPad("fPad3D","NOvA Event Display",0.0,0.0,0.7,0.7,0);
  fPad3D->Draw();

  fView = new evdb::View3D();
  
  fXlo = 0.0;
  fXhi = 0.0;
  fYlo = 0.0;
  fYhi = 0.0;
  fZlo = 0.0;
  fZhi = 0.0;
  fTlo = -10.0;
  fThi = +20.0;

  this->BookHistograms();
}

//......................................................................

void Display3D::BookHistograms()
{
  fRawTHisto = new TH1F("fRawTHisto",";usec;hits",300,-10.0,20.0);
  this->StdHistoOpts(fRawTHisto);

  fMCTHisto = new TH1F("fMCTHisto",";usec;hits",300,-10.0,20.0);
  this->StdHistoOpts(fMCTHisto);
  fMCTHisto->SetLineColor(2);

  fRawQHisto = new TH1F("fRawQHisto",";'MeV';hits",200,0.0,100.0);
  this->StdHistoOpts(fRawQHisto);
  
  fMCQHisto = new TH1F("fMCQHisto",";MeV;hits",200,0.0,100.0);
  this->StdHistoOpts(fMCQHisto);
  fMCQHisto->SetLineColor(2);
}

//......................................................................

void Display3D::StdHistoOpts(TH1F* h)
{
  h->SetLabelFont(42,"X");
  h->SetLabelFont(42,"Y");
  h->SetLabelSize(0.06,"X");
  h->SetLabelSize(0.06,"X");
  
  h->SetTitleFont(52,"X");
  h->SetTitleFont(52,"Y");
  h->SetTitleSize(0.08,"X");
  h->SetTitleSize(0.08,"Y");
  h->SetTitleOffset(0.55,"X");
  h->SetTitleOffset(0.45,"Y");
  
  h->GetXaxis()->SetNdivisions(510);
  h->GetYaxis()->SetNdivisions(505);
}

//......................................................................

Display3D::~Display3D() 
{
  delete fView;  fView  = 0;
  delete fPad3D; fPad3D = 0;
}

//......................................................................

void Display3D::Draw(const char* /*opt*/) 
{

  if (fXlo == fXhi && fYlo == fYhi && fZhi == fZlo) {
    edm::EventHandle& evt  = evdb::IoModule::Instance()->GetEvent();
    std::vector<const rawdata::DAQHeader*> header;
    short int det = rawdata::kFar;
    try{ evt.DAQ().Get("./",header); }
    catch(edm::Exception e){
      std::cerr << "Error retrieving daq header, while looking "<<
        "in Display3D::Draw()" << std::endl;
    }
    if(header.size() > 0) det = header[0]->DetId();

    geo::Geometry&    geom = geo::Geometry::Instance(evt.Header().Run(),det);
    fXlo = -geom.DetHalfWidth();
    fXhi =  geom.DetHalfWidth();
    fYlo = -geom.DetHalfHeight();
    fYhi =  geom.DetHalfHeight();
    fZlo =  0.0;
    fZhi =  geom.DetLength();
  }
  
  // fBomb control how far out to shift the detector faces from their
  // actual positions. The factor should be made adjustable. TODO
  double fBombFac = 0.05;
  fBomb = 
    (fXhi-fXlo)>(fYhi-fYhi) ? fBombFac*(fXhi-fXlo):fBombFac*(fYhi-fYlo);

  fView->Clear();
  fRawTHisto->Reset();
  fMCTHisto-> Reset();
  fRawQHisto->Reset();
  fMCQHisto-> Reset();

  this->DrawOutline();
  this->DrawMCHits();
  this->DrawRawData();

  int irep = 0;
  fPad3D->Clear();
  fPad3D->cd();
  fView->Draw();
  if (fPad3D->GetView()==0) {
    double sz = 1600.0;
    double rmin[] = {-sz,-sz,2500.0-sz};
    double rmax[] = { sz, sz,2500.0+sz};
    TView3D* v = new TView3D(1,rmin,rmax);
    v->SetPerspective();
    v->SetView(226.0,226.0,135.0,irep);
    fPad3D->SetView(v); // ROOT takes ownership of object *v
  }
  if (fGLViewer3D==0) {
    fGLViewer3D = (TGLViewer*)fPad3D->GetViewer3D("ogl");
    fGLViewer3D->SetClearColor(kWhite);
    fGLViewer3D->PreferLocalFrame();
    fGLViewer3D->SetResetCamerasOnUpdate(0);
  }
  else {
    fPad3D->GetViewer3D()->PadPaint(fPad3D);
  }
  fPad3D->Update();

  fTimePad->Clear();
  fTimePad->cd();
  this->SetGoodMax(fRawTHisto,fMCTHisto);
  fRawTHisto->Draw();
  fMCTHisto->Draw("same");
  fTimePad->Update();

  fChargePad->Clear();
  fChargePad->cd();
  this->SetGoodMax(fRawQHisto,fMCQHisto);
  fRawQHisto->Draw();
  fMCQHisto->Draw("same");
  fChargePad->Update();
}

//......................................................................

void Display3D::SetGoodMax(TH1F* h1, TH1F* h2) 
{
  if (h2->GetMaximum()>h1->GetMaximum()) {
    h1->SetMaximum(1.05*h2->GetMaximum());
  }
}

//......................................................................
    
void Display3D::DrawOutline() 
{
  int i, is;
  int cont;
  double x, y, z;
  double x0 = 0.5*(fXlo+fXhi);
  double y0 = 0.5*(fYlo+fYhi);
  double z0 = 0.5*(fZlo+fZhi);
  double dx = 100.0; // cm
  double dy = 100.0; // cm
  double dz = 100.0; // cm
  int c = kBlue-10;
  int w = 1;
  int s = 1;
  // Grid lines marking x-tick marks
  for (cont=1, i=0; cont==1 ; ++i) {
    cont = 0;
    for (is=-1; is<=1; is += 2) {
      x = x0+is*i*dx;
      if (x>=fXlo && x<=fXhi) {
        TPolyLine3D& p = fView->AddPolyLine3D(0,c,w,s);
        p.SetPoint(0,x,fYhi+fBomb,fZlo);
        p.SetPoint(1,x,fYhi+fBomb,fZhi);
        cont = 1;
      }
    }
  }
  // Grid lines marking y-tick marks
  for (cont=1, i=0; cont==1 ; ++i) {
    cont = 0;
    for (is=-1; is<=1; is += 2) {
      y = y0+is*i*dy;
      if (y>=fYlo && y<=fYhi) {
        TPolyLine3D& p = fView->AddPolyLine3D(0,c,w,s);
        p.SetPoint(0,fXhi+fBomb,y,fZlo);
        p.SetPoint(1,fXhi+fBomb,y,fZhi);
        cont = 1;
      }
    }
  }
  // Grid lines marking z-tick marks
  for (cont=1, i=0; cont==1 ; ++i) {
    cont = 0;
    for (is=-1; is<=1; is += 2) {
      z = z0+is*i*dz;
      if (z>=fZlo && z<=fZhi) {
        TPolyLine3D& p1 = fView->AddPolyLine3D(0,c,w,s);
        p1.SetPoint(0,fXhi+fBomb,fYlo,z);
        p1.SetPoint(1,fXhi+fBomb,fYhi,z);

        TPolyLine3D& p2 = fView->AddPolyLine3D(0,c,w,s);
        p2.SetPoint(0,fXlo,fYhi+fBomb,z);
        p2.SetPoint(1,fXhi,fYhi+fBomb,z);
        cont = 1;
      }
    }
  }
  // Complete the box on the otherside
  TPolyLine3D& p1 = fView->AddPolyLine3D(0,c,w,s);
  p1.SetPoint(0,fXlo-fBomb,fYlo,fZlo);
  p1.SetPoint(1,fXlo-fBomb,fYhi,fZlo);
  p1.SetPoint(2,fXlo-fBomb,fYhi,fZhi);
  p1.SetPoint(3,fXlo-fBomb,fYlo,fZhi);
  p1.SetPoint(4,fXlo-fBomb,fYlo,fZlo);

  TPolyLine3D& p3 = fView->AddPolyLine3D(0,c,w,s);
  p3.SetPoint(0,fXhi+fBomb,fYlo,fZlo);
  p3.SetPoint(1,fXhi+fBomb,fYhi,fZlo);
  p3.SetPoint(2,fXhi+fBomb,fYhi,fZhi);
  p3.SetPoint(3,fXhi+fBomb,fYlo,fZhi);
  p3.SetPoint(4,fXhi+fBomb,fYlo,fZlo);

  TPolyLine3D& p2 = fView->AddPolyLine3D(0,c,w,s);
  p2.SetPoint(0,fXlo,fYlo-fBomb,fZlo);
  p2.SetPoint(1,fXhi,fYlo-fBomb,fZlo);
  p2.SetPoint(2,fXhi,fYlo-fBomb,fZhi);
  p2.SetPoint(3,fXlo,fYlo-fBomb,fZhi);
  p2.SetPoint(4,fXlo,fYlo-fBomb,fZlo);

  TPolyLine3D& p4 = fView->AddPolyLine3D(0,c,w,s);
  p4.SetPoint(0,fXlo,fYhi+fBomb,fZlo);
  p4.SetPoint(1,fXhi,fYhi+fBomb,fZlo);
  p4.SetPoint(2,fXhi,fYhi+fBomb,fZhi);
  p4.SetPoint(3,fXlo,fYhi+fBomb,fZhi);
  p4.SetPoint(4,fXlo,fYhi+fBomb,fZlo);
}

//......................................................................

void Display3D::DrawMCHits() 
{
  edm::EventHandle& evt  = evdb::IoModule::Instance()->GetEvent();
  std::vector<const rawdata::DAQHeader*> header;
  short int det = rawdata::kFar;
  try{ evt.DAQ().Get("./",header); }
  catch(edm::Exception e){
    std::cerr << "Error retrieving daq header, while looking "<<
      "in Display3D::DrawMCHits(), using default fardet geometry" << std::endl;
  }
  if(header.size() > 0) det = header[0]->DetId();

  geo::Geometry&    geom = geo::Geometry::Instance(evt.Header().Run(),det);

  std::vector<const sim::FLSHitList*> hitlist;
  try         { evt.DetSim().Get(".",hitlist); }
  catch (...) { return; }

  for (unsigned int i=0; i<hitlist.size(); ++i) {
    for (unsigned int j=0; j<hitlist[i]->fHits.size(); ++j) {
      const sim::FLSHit& hit = hitlist[i]->fHits[j];
      
      // Scale units to MeV and usec
      fMCQHisto->Fill(hit.fEdep*1.0E3);
      fMCTHisto->Fill(hit.fT*1.0E6);
      
      int c;
      double sz;
      double lpos[3], wpos[3];

      lpos[0] = 0.5*(hit.fEntryX+hit.fExitX);
      lpos[1] = 0.5*(hit.fEntryY+hit.fExitY);
      lpos[2] = 0.5*(hit.fEntryZ+hit.fExitZ);
      geom.Plane(hit.fPlaneId).Cell(hit.fCellId).LocalToWorld(lpos,wpos);
      
      // scale the marker size so that volume is proportional to
      // charge
      sz = pow(1.0E0+1.E5*hit.fEdep,0.33);
      c  = Style::ColorFromPDG(hit.fPDG);
      TPolyMarker3D& m = fView->AddPolyMarker3D(1,c,24,sz);
      m.SetPoint(0, wpos[0], wpos[1], wpos[2]);
    }
  }
}

//......................................................................

void Display3D::DrawRawData() 
{
  edm::EventHandle& evt  = evdb::IoModule::Instance()->GetEvent();
  std::vector<const rawdata::DAQHeader*> header;
  short int det = rawdata::kFar;
  try{ evt.DAQ().Get("./",header); }
  catch(edm::Exception e){
    std::cerr << "Error retrieving daq header, while looking "<<
      "in Display3D::RawData(), using default fardet geometry" << std::endl;
  }
  if(header.size() > 0) det = header[0]->DetId();
  geo::Geometry&    geom = geo::Geometry::Instance(evt.Header().Run(),det);
  
  std::vector<const rawdata::RawDigit*> rd;
  try         { evt.Raw().Get(".",rd); }
  catch (...) { return; }

  double pos[3];
  double dx, dy, dz;
  for (unsigned int i=0; i<rd.size(); ++i) {
    const rawdata::RawDigit& r = (*rd[i]);
    const geo::PlaneGeo&     p = geom.Plane(r.GetPlane());
    const geo::CellGeo&      c = p.Cell(r.GetCell());
    
    c.GetCenter(pos);
    for (int j=0; j<r.NADC(); ++j) {
      fRawTHisto->Fill(r.TDC(j));
      fRawQHisto->Fill(0.1*r.ADC(j));
      if (p.View()==geo::kX) {
        dx = c.HalfW();
        dy = 100.0*r.ADC(j)/1024.0;
        dz = c.HalfD();
        TMarker3DBox& m =
          fView->AddMarker3DBox(pos[0],fYhi+fBomb,pos[2], dx, dy, dz);
        m.SetLineColor(fQscale->GetColor(r.ADC(j)));
      }
      else if (p.View()==geo::kY) {
        dx = 100.0*r.ADC(j)/1024.0;
        dy = c.HalfW();
        dz = c.HalfD();
        TMarker3DBox& m =
          fView->AddMarker3DBox(fXhi+fBomb,pos[1],pos[2], dx, dy, dz);
        m.SetLineColor(fQscale->GetColor(r.ADC(j)));
      }
      else abort();
    }
  }
}

---