JMarkovPathDisplayer.cc

Go to the documentation of this file.

/*******************************************************************************
 * \file 3D display of generated photon paths

Not complete yet. It should be updated to use the path weights properly.

*******************************************************************************/

// C++ standard library
#include<iostream>
#include<sstream>
#include<iomanip>
#include<vector>
#include<cmath>
#include<cstdlib>

// ROOT
#include "TRandom3.h"
#include "TFile.h"
#include "TPolyLine3D.h"
#include "TPolyMarker3D.h"
#include "TAxis3D.h"
#include "TView.h"
#include "TView3D.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TH3F.h"

// JPP
#include "Jeep/JParser.hh"
#include "JMarkov/JPhotonPath.hh"
#include "JMarkov/JPhotonPathReader.hh"
using namespace std ;
using namespace JPP;

// this little line removes the make_field macro defined in JParser.hh
// so that we can call the
// JPARSER::make_field(T,string) function ourselves
#undef make_field
using namespace JPP;

int main( int argc, char** argv ) {
  cout << "JMarkovPathDisplayer" << endl 
       << "Written by Martijn Jongen" << endl 
       << endl ; 
  cout << "Type '" << argv[0] << " -h!' to display the command-line options." << endl ;
  cout << endl ;

  vector<string> ifnames ;
  string ofname = "out.root" ;
  int npaths_to_draw ;

  try {
    JParser<string> zap ; // this argument parser can handle strings    
    zap["f"] = make_field(ifnames,"input file name (binary file containing JPhotonPaths)") ;
    zap["o"] = make_field(ofname,"output file name") ;
    zap["n"] = make_field(npaths_to_draw,"max number of paths to draw") = 200 ;

    if (zap.read(argc, argv) != 0) {
      return 1 ;
    }
  }
  catch(const exception &error) {
    // ignore exceptions
  }

  // read the paths from the file
  vector<JMARKOV::JPhotonPath> paths ;
  cout << "Reading files." << endl ;
  int nread_total = 0 ;
  for( vector<string>::iterator it=ifnames.begin() ; it!=ifnames.end() ; ++it ) {
    cout << "Reading '" << *it << "'." << endl ;
    // read the paths from the file
    JMARKOV::JPhotonPathReader reader ;
    reader.open(it->c_str()) ;
    if( !reader.is_open() ) {
      cerr << "FATAL ERROR: unable to open input file '" << *it << "'." << endl ;
      exit(1) ;
    }
    int nread = 0 ;
    JMARKOV::JPhotonPath* p = NULL ;
    while( reader.hasNext() && nread_total<=npaths_to_draw ) {
      p = reader.next() ;
      paths.push_back(*p) ;
      ++nread ;
      ++nread_total ;
    }
    cout << "--> read " << nread << " paths." << endl ;
    reader.close() ;
  }
  cout << endl ;
  cout << "Done reading paths." << endl ;
  cout << endl ;

  if( nread_total == 0 ) {
    cerr << "FATAL ERROR: could not read any JPhotonPaths from the file(s)." << endl ;
    exit(1) ;
  }

  cout << "Done reading files. Read " << nread_total << " paths from it." << endl ;
  cout << endl ;

  // variables to find out how large our view has to be
  double xmin    = 1.0/0.0  ;
  double xmax    = -1.0/0.0 ;
  double ymin    = 1.0/0.0  ;
  double ymax    = -1.0/0.0 ;
  double zmin    = 1.0/0.0  ;
  double zmax    = -1.0/0.0 ;

  // location of the source and target
  JGEOMETRY3D::JPosition3D xsource ;
  JGEOMETRY3D::JPosition3D xtarget ;
  
  TRandom3* ran = new TRandom3(0) ;

  // limit number of paths to draw
  npaths_to_draw = min(npaths_to_draw,nread_total) ;

  vector<TPolyLine3D*> to_draw ;
  for( int i=0 ; i<npaths_to_draw ; ++i ) {
    int nvert = paths[i].size() ; // number of vertices

    // create graphics object TPolyLine that will be used to draw the path
    TPolyLine3D* pl = new TPolyLine3D(nvert) ;
    pl->SetLineWidth(1) ;
    pl->SetLineColor(kCyan) ;
    to_draw.push_back(pl) ;

    for( int j=0 ; j<nvert ; ++j ) {
      double x = paths[i][j].getX() ;
      double y = paths[i][j].getY() ;
      double z = paths[i][j].getZ() ;
      if( x>xmax ) xmax = x ;
      if( y>ymax ) ymax = y ;
      if( x<xmin ) xmin = x ;
      if( y<ymin ) ymin = y ;
      if( z<zmin )  zmin = z  ;
      if( z>zmax )  zmax = z  ;
      pl->SetPoint( j, x, y, z ) ;
    }
    xsource = paths[i][0] ;
    xtarget = paths[i][nvert-1] ;
  }

  // adjust x,y, z min and max such that the length of the x, y and z
  // ranges are equal
  double range = 0 ;
  if( xmax-xmin > range ) range = xmax-xmin ;
  if( ymax-ymin > range ) range = ymax-ymin ;
  if( zmax-zmin > range ) range = zmax-zmin ;
  range *= 1.1 ;
  double extra ;
  // adjust x
  extra = 0.5*(range-xmax-xmin) ;
  xmax += extra ;
  xmin -= extra ;
  // adjust y
  extra = 0.5*(range-ymax-ymin) ;
  ymax += extra ;
  ymin -= extra ;
  // adjust z
  extra = 0.5*(range-zmax-zmin) ;
  zmax += extra ;
  zmin -= extra ;

  cout << "DIMENSIONS: " << endl
       << "x = [ " << xmin << " , " << xmax << " ] [m]" << endl 
       << "y = [ " << ymin << " , " << xmax << " ] [m]" << endl 
       << "z = [ " << zmin << " , " << xmax << " ] [m]" << endl ;
  cout << endl ;

  // prepare TH3F to show the local path density
  cout << "Creating and filling TH3F" << endl ;
  int nbinsh3 = 200 ;
  TH3F* h3 = new TH3F("hPathDensity_Z",
                      "Local path density;X [m];Y [m];Z [m]",
                      nbinsh3,xmin,xmax,
                      nbinsh3,ymin,ymax,
                      nbinsh3,zmin,zmax ) ;
  h3->SetOption("colz") ;
  for( Int_t zbin=1 ; zbin<=nbinsh3 ; ++zbin ) {
    double z = h3->GetZaxis()->GetBinCenter(zbin) ;
    for( int i=0 ; i<nread_total ; ++i ) {
      vector<JGEOMETRY3D::JPosition3D> coords = paths[i].getPointsWithZ(z) ;
      for( vector<JGEOMETRY3D::JPosition3D>::iterator it=coords.begin() ; it!=coords.end() ; ++it ) {
        h3->Fill( it->getX(), it->getY(), it->getZ() ) ;
      }
    }
  }
  cout << "TH3F filled!" << endl ;
  cout << endl ;

  // prepare another TH3F to show the local path density
  cout << "Creating and filling TH3F" << endl ;
  TH3F* h3X = new TH3F("hPathDensity_X",
                      "Local path density;X [m];Y [m];Z [m]",
                      nbinsh3,xmin,xmax,
                      nbinsh3,ymin,ymax,
                      nbinsh3,zmin,zmax ) ;
  h3X->SetOption("colz") ;
  for( Int_t xbin=1 ; xbin<=nbinsh3 ; ++xbin ) {
    double x = h3X->GetXaxis()->GetBinCenter(xbin) ;
    for( int i=0 ; i<nread_total ; ++i ) {
      vector<JGEOMETRY3D::JPosition3D> coords = paths[i].getPointsWithX(x) ;
      for( vector<JGEOMETRY3D::JPosition3D>::iterator it=coords.begin() ; it!=coords.end() ; ++it ) {
        h3X->Fill( it->getX(), it->getY(), it->getZ() ) ;
      }
    }
  }
  cout << "TH3F filled!" << endl ;
  cout << endl ;

  // prepare another TH3F to show the local path density
  /*cout << "Creating and filling TH3F" << endl ;
  nbinsh3 = 80 ;
  TH3F* h3sphere = new TH3F("hPathDensity_sphere",
                      "Local path density;X [m];Y [m];Z [m]",
                      nbinsh3,xmin,xmax,
                      nbinsh3,ymin,ymax,
                      nbinsh3,zmin,zmax ) ;
  h3sphere->SetOption("colz") ;
  for( Int_t xbin=1 ; xbin<=nbinsh3 ; ++xbin ) {
    double x = h3sphere->GetXaxis()->GetBinCenter(xbin) ;
    cout << "x = " << x << endl ;
    for( Int_t ybin=1 ; ybin<=nbinsh3 ; ++ybin ) {
      double y = h3sphere->GetYaxis()->GetBinCenter(ybin) ;
      for( Int_t zbin=1 ; zbin<=nbinsh3 ; ++zbin ) {
        double z = h3sphere->GetZaxis()->GetBinCenter(zbin) ;
        Int_t bin = h3sphere->GetBin(xbin,ybin,zbin) ;
        for( int i=0 ; i<nread_total ; ++i ) {
          JPosition3D pos(x,y,z) ;
          if ( paths[i].hitsSphere( pos, range/nbinsh3 ) ) {
            h3sphere->Fill( x,y,z ) ;
          }
        }
      }
    }
  }
  cout << "TH3F filled!" << endl ;
  cout << endl ;*/

  // make a canvas that shows the 3D paths
  double csize = 2048 ;
  TCanvas* cDisplay = new TCanvas("cDisplay","Display",10,10,csize,csize) ;
  TPad* pad1 = new TPad("pad1","pad1",0,0,1,1,kBlue+4) ;
  pad1->Draw() ;
  pad1->cd() ;
  TView3D* view = (TView3D*)TView::CreateView(1) ;
  //double scale = 1.0/sqrt(2.0) ; // scale to use more of the visible area
  view->SetRange(xmin,ymin,zmin,
                 xmax,ymax,zmax   ) ;
  //TAxis3D rulers ;
  //rulers.Draw();

  // draw source and target position
  TPolyMarker3D* pm = new TPolyMarker3D(2,1) ;
  pm->SetPoint(0, xsource.getX(), xsource.getY(), xsource.getZ() ) ;
  pm->SetPoint(1, xtarget.getX(), xtarget.getY(), xtarget.getZ() ) ;
  pm->SetMarkerSize(2) ;
  pm->SetMarkerColor(kRed) ;
  pm->SetMarkerStyle(20) ;
  pm->Draw() ;

  for( vector<TPolyLine3D*>::iterator it=to_draw.begin() ; it!=to_draw.end() ; ++it ) {
    (*it)->Draw() ;
  }

  cDisplay->SaveAs("display.png") ;
  cDisplay->SaveAs("display.pdf") ;
  cout << "Output in 'display.png'." << endl ;
  cout << endl ;

  /*{
    // delete the file if it is there (otherwise there will be complaints)
    gSystem->Unlink("myGif.gif") ;
    const int nrot = 40 ;
    for( int i=0 ; i<nrot ; ++i ) {
      view->RotateView( i*360.0/nrot,90 ) ;
      if( i==nrot-1 ) {
        cDisplay->Print("myGif.gif++10++") ;
      } else {
        cDisplay->Print("myGif.gif+10") ;
      }
    }
    
    }
\author mjongen
  */
  TFile* fout = new TFile(ofname.c_str(),"recreate") ;
  cDisplay->Write() ;
  h3->Write() ;
  h3X->Write() ;
  //h3sphere->Write() ;
  fout->Close() ;
  cout << "Output written to '" << ofname << "'." << endl ;
  cout << endl ;
  cout << "Done!" << endl ;
}