source: snovis/trunk/source/G4Lab/cxx/SoG4Trajectories.cxx

// this :
#include <G4Lab/SoG4Trajectories.h>

// Inventor :
#include <Inventor/SoPickedPoint.h>
#include <Inventor/errors/SoDebugError.h>
#include <Inventor/actions/SoGLRenderAction.h>
#include <Inventor/actions/SoRayPickAction.h>
#include <Inventor/actions/SoCallbackAction.h>
#include <Inventor/actions/SoCallbackAction.h>
#include <Inventor/sensors/SoFieldSensor.h>
#include <Inventor/elements/SoLazyElement.h>
#include <Inventor/bundles/SoMaterialBundle.h>

// HEPVis :
#include <HEPVis/actions/SoAlternateRepAction.h>

#ifdef WIN32
#undef pascal // Clash between windef.h and Geant4/SystemOfnits.hh
#endif

// Geant4 :
#include <G4RunManager.hh>

// G4Lab :
#include <G4Lab/Trajectory.h>

//WIN32 : have HEPVis/SbGL.h after G4SteppingManager.h.
// Else the windows.h will clash with std::max.
#include <HEPVis/SbGL.h> 

#define MINIMUM(a,b) ((a)<(b)?a:b)
#define MAXIMUM(a,b) ((a)>(b)?a:b)

SO_NODE_SOURCE(SoG4Trajectories)
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::initClass(
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  SO_NODE_INIT_CLASS(SoG4Trajectories,SoShape,"Shape");
}
//////////////////////////////////////////////////////////////////////////////
SoG4Trajectories::SoG4Trajectories(
)
:fRunManager(0)
,fContainer(0)
,fSensor(0)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  SO_NODE_CONSTRUCTOR(SoG4Trajectories);
  SO_NODE_ADD_FIELD(model,(ALL));
  SO_NODE_ADD_FIELD(timeStart,(0));
  SO_NODE_ADD_FIELD(timeInterval,(0.05F));
  SO_NODE_ADD_FIELD(timeSteps,(100));
  SO_NODE_ADD_FIELD(timeIndex,(0));

  SO_NODE_ADD_FIELD(timeMin,(0));
  SO_NODE_ADD_FIELD(timeMax,(0));
  SO_NODE_ADD_FIELD(verbose,(FALSE));
  
  SO_NODE_ADD_FIELD(alternateRep,(NULL));

  SO_NODE_DEFINE_ENUM_VALUE(Model,ALL);
  SO_NODE_DEFINE_ENUM_VALUE(Model,TIMED);
  
  SO_NODE_SET_SF_ENUM_TYPE(model,Model);

  fSensor = new SoFieldSensor(sensorCB,this);
  fSensor->attach(&timeInterval);

  timeSteps.addAuditor(&timeInterval,SoNotRec::FIELD);
  timeStart.addAuditor(&timeInterval,SoNotRec::FIELD);

}
//////////////////////////////////////////////////////////////////////////////
SoG4Trajectories::SoG4Trajectories(
 G4RunManager* aRunManager
)
:fRunManager(aRunManager)
,fContainer(0)
,fSensor(0)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  SO_NODE_CONSTRUCTOR(SoG4Trajectories);
  SO_NODE_ADD_FIELD(model,(ALL));
  SO_NODE_ADD_FIELD(timeStart,(0));
  SO_NODE_ADD_FIELD(timeInterval,(0.05F));
  SO_NODE_ADD_FIELD(timeSteps,(100));
  SO_NODE_ADD_FIELD(timeIndex,(0));

  SO_NODE_ADD_FIELD(timeMin,(0));
  SO_NODE_ADD_FIELD(timeMax,(0));
  SO_NODE_ADD_FIELD(verbose,(FALSE));
  
  SO_NODE_ADD_FIELD(alternateRep,(NULL));

  SO_NODE_DEFINE_ENUM_VALUE(Model,ALL);
  SO_NODE_DEFINE_ENUM_VALUE(Model,TIMED);
  
  SO_NODE_SET_SF_ENUM_TYPE(model,Model);

  fSensor = new SoFieldSensor(sensorCB,this);
  fSensor->attach(&timeInterval);

  timeSteps.addAuditor(&timeInterval,SoNotRec::FIELD);
  timeStart.addAuditor(&timeInterval,SoNotRec::FIELD);

}
//////////////////////////////////////////////////////////////////////////////
SoG4Trajectories::~SoG4Trajectories(
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  fSensor->detach();
  timeSteps.removeAuditor(&timeInterval,SoNotRec::FIELD);
  timeStart.removeAuditor(&timeInterval,SoNotRec::FIELD);
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::GLRender (
 SoGLRenderAction* aAction
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(!fRunManager) return;
  const G4Event* event = fRunManager->GetCurrentEvent();
  if(!event) return;
  G4TrajectoryContainer* trajectoryContainer = event->GetTrajectoryContainer();
  if(!trajectoryContainer) return;

  SoState* state = aAction->getState();
  state->push();

  SoLazyElement::setLightModel(state, SoLazyElement::BASE_COLOR);

  glPushAttrib(GL_ENABLE_BIT);
  glPushAttrib(GL_CURRENT_BIT);

  SoMaterialBundle mb(aAction);
  mb.sendFirst();

  const SbColor& color = SoLazyElement::getDiffuse(aAction->getState(),0);
  float red,green,blue;
  color.getValue(red,green,blue);
  glDisable(GL_LIGHTING);
  glColor3f(red,green,blue);

  if((SoG4Trajectories::Model)model.getValue()==SoG4Trajectories::TIMED) {

    if(fContainer!=trajectoryContainer) { // Event had changed.
      timeOrder(trajectoryContainer);
      fContainer = trajectoryContainer;
    }

    int timen = fTimeOrderedPoints.size();
    int ti = timeIndex.getValue();
    if((ti>=0)&&(ti<timen)) {

      const std::vector<PointIdentifier>& ids = fTimeOrderedPoints[ti];
      int idn = ids.size();
      if(verbose.getValue()==TRUE)
        SoDebugError::postInfo("SoG4RunManager::GLRender",
                               "debug : render slice %d with %d points.",
                               ti,idn);

      for(int t=0;t<=ti;t++) {
        const std::vector<PointIdentifier>& ids = fTimeOrderedPoints[t];
        int idn = ids.size();
            
        // First pass to draw points :
        glBegin(GL_POINTS);
        int i;
        for(i=0;i<idn;i++) {
          int itraj = ids[i].fTrajectory;
          int ipoint = ids[i].fPoint;
          G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)itraj];
          int pointn = trajectory->GetPointEntries();
          if( (pointn==1) || (ipoint==pointn-1) ) {
            G4VTrajectoryPoint* tp = trajectory->GetPoint(ipoint);
            G4ThreeVector pos = tp->GetPosition();
            glVertex3d(pos.x(),pos.y(),pos.z());
          }
        }
        glEnd();
        
        // Draw segments :
        glBegin(GL_LINES);
        for(i=0;i<idn;i++) {
          int itraj = ids[i].fTrajectory;
          int ipoint = ids[i].fPoint;
          G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)itraj];
          int pointn = trajectory->GetPointEntries();
          if( (pointn==1) || (ipoint==pointn-1) ) {
          } else {
            G4VTrajectoryPoint* tp = trajectory->GetPoint(ipoint);
            G4ThreeVector pos = tp->GetPosition();
            glVertex3d(pos.x(),pos.y(),pos.z());
            tp = trajectory->GetPoint(ipoint+1);
            pos = tp->GetPosition();
            glVertex3d(pos.x(),pos.y(),pos.z());
          }
        }
        glEnd();
      }
    }

  } else {

    int number = trajectoryContainer->entries();
    for(int index=0;index<number;index++) {
      G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)index];
      int pointn = trajectory->GetPointEntries();
      if(pointn==1) {
        glBegin(GL_POINTS);
        G4VTrajectoryPoint* tp = trajectory->GetPoint(0);
        G4ThreeVector pos = tp->GetPosition();
        glVertex3d(pos.x(),pos.y(),pos.z());
        glEnd();
      } else {
        glBegin(GL_LINE_STRIP);
        for (int i = 0; i < pointn ; i++) {
          G4VTrajectoryPoint* tp = trajectory->GetPoint(i);
          G4ThreeVector pos = tp->GetPosition();
          glVertex3d(pos.x(),pos.y(),pos.z());
        }
        glEnd();
      }
    }

  }

  glPopAttrib(); //GL_CURRENT_BIT
  glPopAttrib(); //GL_ENABLE_BIT

  state->pop();
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::rayPick (
 SoRayPickAction* aAction
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(!shouldRayPick(aAction)) return;

  aAction->setObjectSpace();

  if(fRunManager) {
    const G4Event* event = fRunManager->GetCurrentEvent();
    if(event) {
      G4TrajectoryContainer* trajectoryContainer = 
        event->GetTrajectoryContainer();
      if(trajectoryContainer) {
        int number = trajectoryContainer->entries();
        for(int index=0;index<number;index++) {
          G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)index];
          int pointn = trajectory->GetPointEntries();
          if(!pointn) continue;
          G4VTrajectoryPoint* tp = trajectory->GetPoint(0);
          G4ThreeVector pos = tp->GetPosition();
          SbVec3f v0;
          v0.setValue((float)pos.x(),(float)pos.y(),(float)pos.z());
          for (int i = 1; i < pointn ; i++) {
            tp = trajectory->GetPoint(i);
            pos = tp->GetPosition();
            SbVec3f v1((float)pos.x(),(float)pos.y(),(float)pos.z());
            SbVec3f isect;
            SbBool hit = aAction->intersect(v0, v1, isect);
            if(hit && aAction->isBetweenPlanes(isect)) {
              SoPickedPoint* pp = aAction->addIntersection(isect);
              if(pp) {
                pp->setMaterialIndex(0);
                pp->setObjectNormal(SbVec3f(0.0f, 0.0f, 1.0f));
              }
            }
            v0 = v1;
          }
        }
      }
    }
  }
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::computeBBox (
 SoAction* //aAction
,SbBox3f& aBox
,SbVec3f& aCenter
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(fRunManager) {
    const G4Event* event = fRunManager->GetCurrentEvent();
    if(event) {
      G4TrajectoryContainer* trajectoryContainer = 
        event->GetTrajectoryContainer();
      if(trajectoryContainer) {
        int number = trajectoryContainer->entries();
        double xmin = 0;
        double ymin = 0;
        double zmin = 0;
        double xmax = 0;
        double ymax = 0;
        double zmax = 0;
        bool first = true;
        for(int index=0;index<number;index++) {
          G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)index];
          int pointn = trajectory->GetPointEntries();
          for (int i = 0; i < pointn ; i++) {
            G4VTrajectoryPoint* tp = trajectory->GetPoint(i);
            G4ThreeVector pos = tp->GetPosition();
            if(first) {
              xmax = xmin = pos.x();
              ymax = ymin = pos.y();
              zmax = zmin = pos.z();
              first = false;
            } else {
              xmin = MINIMUM(xmin,pos.x());
              ymin = MINIMUM(ymin,pos.y());
              zmin = MINIMUM(zmin,pos.z());
              xmax = MAXIMUM(xmax,pos.x());
              ymax = MAXIMUM(ymax,pos.y());
              zmax = MAXIMUM(zmax,pos.z());
            }
          }
        }
        if(!first) {
          aBox.setBounds(SbVec3f((float)xmin,(float)ymin,(float)zmin),
                         SbVec3f((float)xmax,(float)ymax,(float)zmax));
        }
      }
    }
  }
  aCenter = aBox.getCenter();
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::generatePrimitives(
 SoAction* //aAction
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::flush(
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  glFlush();
  glFinish();
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::timeOrder(
 G4TrajectoryContainer* aContainer
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  fTimeOrderedPoints.clear();

  if((SoG4Trajectories::Model)model.getValue()
     !=SoG4Trajectories::TIMED) return;

  if(!aContainer) return;

  int number = aContainer->entries();
  if(number<=0) return;

  // Compute min, max time :
  double minTime = 0;
  double maxTime = 0;
  int numberOfPoints = 0;
  for(int index=0;index<number;index++) {
    IGeant4Trajectory* trajectory = 
      dynamic_cast<IGeant4Trajectory*>((*aContainer)[(size_t)index]);
    if(!trajectory) {
      fTimeOrderedPoints.clear();
      return;
    }
    unsigned int pointn = trajectory->pointEntries();
    for (unsigned int i = 0; i < pointn ; i++) {
      double globalTime = trajectory->pointGlobalTime(i);
      if(numberOfPoints==0) {
        minTime = globalTime;
        maxTime = globalTime;
      } else {
        minTime = MINIMUM(minTime,globalTime);
        maxTime = MAXIMUM(maxTime,globalTime);
      }
      numberOfPoints++;
    }
  }    



  SbBool flag = enableNotify(FALSE);
  timeMin.setValue((float)minTime);
  timeMax.setValue((float)maxTime);
  enableNotify(flag);

  // Construct the time ordered point list :
  int timen = timeSteps.getValue();
  if(timen>0) {
    //double deltaTime = (maxTime - minTime)/timen;
    double deltaTime = timeInterval.getValue();
    int count = 0;
    for(int t=0;t<=timen;t++) {
      fTimeOrderedPoints.push_back(std::vector<PointIdentifier>());
      std::vector<PointIdentifier>& ids = fTimeOrderedPoints[t];
      double tmin = timeStart.getValue() + t * deltaTime;
      double tmax = tmin + deltaTime;
      for(int index=0;index<number;index++) {
        IGeant4Trajectory* trajectory = 
          dynamic_cast<IGeant4Trajectory*>((*aContainer)[(size_t)index]);
        unsigned int pointn = trajectory->pointEntries();
        for (unsigned int i = 0; i < pointn ; i++) {
          double globalTime = trajectory->pointGlobalTime(i);
          if((tmin<=globalTime) && (globalTime<tmax)) {
            ids.push_back(PointIdentifier(index,i));
            count++;
          }
        }
      }
    }

    /*    
    if(count!=numberOfPoints) {
      SoDebugError::postWarning("SoG4RunManager::GLRender",
        "point count problem ; %d points, but %d ordered.",
        numberOfPoints,count);
    }
    */
  }

  if(verbose.getValue()==TRUE)
    SoDebugError::postInfo("SoG4RunManager::GLRender",
    "points : %d, minTime : %g, maxTime : %g, timeStart : %g, timeInterval : %g, timeSteps : %d.",
    numberOfPoints,minTime,maxTime,
    timeStart.getValue(),timeInterval.getValue(),timeSteps.getValue());

}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::sensorCB (
 void* aData
,SoSensor* //aSensor
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  SoG4Trajectories* This = (SoG4Trajectories*)aData;

  if(!This->fRunManager) return;
  const G4Event* event = This->fRunManager->GetCurrentEvent();
  if(!event) return;

  G4TrajectoryContainer* trajectoryContainer = 
    event->GetTrajectoryContainer();
  if(!trajectoryContainer) return;

  This->timeOrder(trajectoryContainer);

}

#include <Inventor/nodes/SoSeparator.h>
#include <Inventor/nodes/SoCoordinate3.h>
//#include <Inventor/nodes/SoLineSet.h>
#include <Inventor/nodes/SoIndexedLineSet.h>
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::generateAlternateRep(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(!fRunManager) return;
  const G4Event* event = fRunManager->GetCurrentEvent();
  if(!event) return;
  G4TrajectoryContainer* trajectoryContainer = event->GetTrajectoryContainer();
  if(!trajectoryContainer) return;

  /* With SoLineSets :
  int number = trajectoryContainer->entries();
  for(int index=0;index<number;index++) {
    G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)index];
    int pointn = trajectory->GetPointEntries();
    if(pointn>=1) {
      SbVec3f* points = new SbVec3f[pointn];
      for (int i = 0; i < pointn ; i++) {
        G4VTrajectoryPoint* tp = trajectory->GetPoint(i);
        G4ThreeVector pos = tp->GetPosition();
        points[i].setValue((float)pos.x(),(float)pos.y(),(float)pos.z());
      }

      // Scene graph :
      SoSeparator* separator = new SoSeparator;

      SoCoordinate3* coordinate3 = new SoCoordinate3;
      coordinate3->point.setValues(0,pointn,points);
      delete [] points;
      separator->addChild(coordinate3);

      SoLineSet* lineSet = new SoLineSet;
      lineSet->numVertices.setValues(0,1,&pointn);
      separator->addChild(lineSet);

      alternateRep.setValue(separator);
    }
  } */

  // With one SoIndexedLineSet :
  int npoint = 0;
  int ncoord = 0;
 {int number = trajectoryContainer->entries();
  for(int index=0;index<number;index++) {
    G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)index];
    int pointn = trajectory->GetPointEntries();
    if(pointn>=1) {
      npoint += pointn;
      ncoord += pointn + 1;
    }
  }}
  if(npoint<=0) return;

  SbVec3f* points = new SbVec3f[npoint];
  int32_t* coords = new int32_t[ncoord];

  int ipoint = 0;
  int icoord = 0;

  int number = trajectoryContainer->entries();
  for(int index=0;index<number;index++) {
    G4VTrajectory* trajectory = (*trajectoryContainer)[(size_t)index];
    int pointn = trajectory->GetPointEntries();
    if(pointn>=1) {
      for (int i = 0; i < pointn ; i++) {
        G4VTrajectoryPoint* tp = trajectory->GetPoint(i);
        G4ThreeVector pos = tp->GetPosition();
        points[ipoint].setValue((float)pos.x(),(float)pos.y(),(float)pos.z());
        coords[icoord] = ipoint;
        ipoint++;
        icoord++;
      }
      coords[icoord] = SO_END_LINE_INDEX;
      icoord++;
    }
  }

  SoSeparator* separator = new SoSeparator;

  SoCoordinate3* coordinate3 = new SoCoordinate3;
  coordinate3->point.setValues(0,ipoint,points);
  delete [] points;
  separator->addChild(coordinate3);

  SoIndexedLineSet* indexedLineSet = new SoIndexedLineSet;
  indexedLineSet->coordIndex.setValues(0,icoord,coords);
  delete [] coords;
  separator->addChild(indexedLineSet);

  alternateRep.setValue(separator);
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::clearAlternateRep(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  alternateRep.setValue(NULL);
}
//////////////////////////////////////////////////////////////////////////////
void SoG4Trajectories::doAction(
 SoAction* aAction
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  SO_ALTERNATEREP_DO_ACTION(aAction)
  SoShape::doAction(aAction);
}