DetectorView.cxx

Go to the documentation of this file.

#include <iostream>
#include <vector>
#include "TCanvas.h"
#include "TPad.h"
#include "TH3F.h"
#include "TLine.h"
#include "TBox.h"
#include "TPolyLine3D.h"
#include "RawData/RawDigit.h"
#include "RawData/DAQHeader.h"
#include "EventDataModel/EventHandle.h"
#include "EventDisplayBase/IoModule.h"
#include "EventDisplayBase/View3D.h"
#include "EventDisplay/DetectorView.h"
#include "EventDisplay/DetectorViewOption.h"
#include "EventDisplay/Style.h"
#include "Geometry/Geometry.h"
#include "Geometry/PlaneGeo.h"
#include "Geometry/CellGeo.h"
#include "RecoBase/PlaneCluster.h"
#include "RecoBase/TrackBase.h"
#include "Simulation/MCTruth.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleList.h"

using namespace evd;
using namespace recobase;

DetectorView::DetectorView(TGMainFrame* mf) : evdb::Canvas(mf) 
{
  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 DetectorView::DetectorView(), using default fardet geometry" << std::endl;
  }
  if(header.size() > 0) det = header[0]->DetId();
  geo::Geometry *geom = &geo::Geometry::Instance(evt.Header().Run(),det);

  evdb::Canvas::fCanvas->cd();

  f3DviewPad = new TPad("fYviewPad","Yview",0.0,0.0,0.75,0.9);
  f3DviewPad->Draw();

  f3DviewPad->cd();
  f3DviewHisto = new TH3F("f3DviewHisto",";z (cm);x (cm); y(cm)",
                          200,0.,geom->DetLength(),
                          32,-geom->DetHalfWidth(), geom->DetHalfWidth(),
                          32,-geom->DetHalfHeight(),geom->DetHalfHeight());
  f3DviewHisto->Draw();

  f3Dview = new evdb::View3D();

  fOption = &DetectorViewOption::Instance();
}

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

DetectorView::~DetectorView() { }

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

void DetectorView::DrawMCTruth(edm::EventHandle& evt)
{
  double zmin = 1e20, zmax = -1e20;
  double xmin = 1e20, xmax = -1e20;
  double ymin = 1e20, ymax = -1e20;
  
  std::vector<const sim::MCTruth*> mctruth;
  
  try { evt.MC().Get(".",mctruth); }
  catch (edm::Exception e) { }
  
  for (int i=0; i<mctruth[0]->NParticles(); ++i) {
    const TParticle& p = mctruth[0]->GetParticle(i);
    
    double xyz0[3];
    xyz0[0] = p.Vx();
    xyz0[1] = p.Vy();
    xyz0[2] = p.Vz();
    
    double xyz1[3];
    // Incident particles
    if (p.GetStatusCode()==0) {
      xyz1[0] = xyz0[0] - p.Px()/2E-3;
      xyz1[1] = xyz0[1] - p.Py()/2E-3;
      xyz1[2] = xyz0[2] - p.Pz()/2E-3;
    }
    // Particles to be tracked
    else if (p.GetStatusCode()==1) {
      xyz1[0] = xyz0[0] + p.Px()/2E-3;
      xyz1[1] = xyz0[1] + p.Py()/2E-3;
      xyz1[2] = xyz0[2] + p.Pz()/2E-3;
    }
    // Intermediate states
    else continue;
    
    // Find limits of display
    if (xyz0[0]<xmin) xmin = xyz0[0];
    if (xyz1[0]<xmin) xmin = xyz1[0];
    if (xyz0[0]>xmax) xmax = xyz0[0];
    if (xyz1[0]>xmax) xmax = xyz1[0];
    if (xyz0[1]<ymin) ymin = xyz0[1];
    if (xyz1[1]<ymin) ymin = xyz1[1];
    if (xyz0[1]>ymax) ymax = xyz0[1];
    if (xyz1[1]>ymax) ymax = xyz1[1];
    if (xyz0[2]<zmin) zmin = xyz0[2];
    if (xyz1[2]<zmin) zmin = xyz1[2];
    if (xyz0[2]>zmax) zmax = xyz0[2];
    if (xyz1[2]>zmax) zmax = xyz1[2];
    
    int ic = Style::ColorFromPDG(p.GetPdgCode());
    int iw = Style::LineWidthFromPDG(p.GetPdgCode());
    int is = Style::LineStyleFromPDG(p.GetPdgCode());

    TPolyLine3D& line = f3Dview->AddPolyLine3D(100,ic,iw,is);
    double dz = 2.;
    int npnt = (int)((xyz1[2]-xyz0[2])/dz); // 2 cm step sizes
    double dxdz = p.Px()/p.Pz();
    double dydz = p.Py()/p.Pz();
    double x, y, z;
    for (int i=0; i<npnt; ++i) {
      x = xyz0[0] + i*dz*dxdz;
      y = xyz0[1] + i*dz*dydz;
      z = xyz0[2] + i*dz;
      line.SetPoint(i,z,x,y);
    }
  }

  std::vector<const sim::ParticleList*> particleList;
  try {
    evt.DetSim().Get(".",particleList);
  }
  catch (edm::Exception e) { }
  
  std::cerr << "Interactions in this record=" 
            << particleList.size()<<std::endl;
  
  for (unsigned int iev=0; iev<particleList.size(); iev++) {
    std::vector<sim::Particle> particles=particleList[iev]->fParticles;
    std::cerr <<"Event "<<iev<<" has "<<particles.size()
              <<" particles."<<std::endl;
    for (unsigned int i=0; i<particles.size(); ++i) {
      std::cerr << particles[i].PdgCode()<<" ";
      if ( (i>0 && i%10==0) || i==particles.size()-1) std::cerr << std::endl;
    }
  }

  // Force the x and y views to have the same vertical scale
  if (xmax-xmin > ymax-ymin) {
    double tmp = 0.5*((xmax-xmin)-(ymax-ymin));
    ymax += tmp;
    ymin -= tmp;
  }
  else {
    double tmp = 0.5*((ymax-ymin)-(xmax-xmin));
    xmax += tmp;
    xmin -= tmp;
  }
  xmin -= 100;  xmax += 100;
  ymin -= 100;  ymax += 100;
  zmin -= 100;  zmax += 100;
  
  f3DviewHisto->GetXaxis()->SetRangeUser(zmin,zmax);
  f3DviewHisto->GetYaxis()->SetRangeUser(xmin,xmax);
  f3DviewHisto->GetZaxis()->SetRangeUser(ymin,ymax);

}

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

void DetectorView::Draw(const char* /*opt*/) 
{
  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 DetectorView::Draw(), using default fardet geometry" << std::endl;
  }
  if(header.size() > 0) det = header[0]->DetId();
  fGeo = &geo::Geometry::Instance(evt.Header().Run(),det);
  
  f3Dview->Clear();

  if (fOption->fDrawMCTruthVectors)
    DrawMCTruth(evt);

  f3DviewPad->cd();
  f3DviewHisto->Draw();
  f3Dview->Draw();
  
  evdb::Canvas::fCanvas->Update();
}

---