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

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

// Inventor :
#include <Inventor/nodes/SoSeparator.h>
#include <Inventor/nodes/SoLightModel.h>
#include <Inventor/nodes/SoDrawStyle.h>
#include <Inventor/nodes/SoCoordinate3.h>
#include <Inventor/nodes/SoIndexedLineSet.h>
#include <Inventor/engines/SoTimeCounter.h>
#include <Inventor/engines/SoCalculator.h>

// HEPVis :
#include <HEPVis/nodes/SoHighlightMaterial.h>

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

// Geant4 :
#include <G4RunManager.hh>
#include <G4VTrajectoryPoint.hh>
#include <G4AttValue.hh>
#include <G4UnitsTable.hh>

// Slash :
#include <Slash/Core/ISession.h>
#include <Slash/Data/IIterator.h>

// Lib :
#include <Lib/Debug.h>
#include <Lib/Out.h>
#include <Lib/Value.h>
#include <Lib/smanip.h>
#include <Lib/sout.h>

// G4Lab :
#include <G4Lab/Interfaces/IGeant4Trajectory.h>
#include <G4Lab/SoG4Trajectories.h>
#include <G4Lab/BestUnit.h>

//////////////////////////////////////////////////////////////////////////////
G4Lab::TrajectoryAccessor::TrajectoryAccessor(
 Slash::Core::ISession& aSession
,G4RunManager* aRunManager
)
:OnX::InventorAccessor(aSession)
,fType("Trajectory")
,fRunManager(aRunManager)
,fAttDefs(0)
,fModeling("immediate_all")
,fIGeant4Trajectory(false)
,fSeparator(0)
,fCoordinate(0)
,fIndex(0)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  Lib::Debug::increment("G4Lab::TrajectoryAccessor");
}
//////////////////////////////////////////////////////////////////////////////
G4Lab::TrajectoryAccessor::~TrajectoryAccessor(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  Lib::Debug::decrement("G4Lab::TrajectoryAccessor");
}
//////////////////////////////////////////////////////////////////////////////
std::string G4Lab::TrajectoryAccessor::name(
) const
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  return fType;
}
namespace G4Lab {
  class TrajectoryIterator : public virtual Slash::Data::IIterator {
  public: //Slash::Data::IIterator
    virtual Slash::Data::IAccessor::Data object() {
      if(fIterator==fVector.end()) return 0;
      return *fIterator;
    }
    virtual void next() { ++fIterator; }
    virtual void* tag() { return 0;}
  public:
    TrajectoryIterator(G4TrajectoryContainer& aVector)
      :fVector(*(aVector.GetVector())) {
      fIterator = fVector.begin();
    }
    virtual ~TrajectoryIterator() {}
  private:
    TrajectoryVector& fVector;
    TrajectoryVector::iterator fIterator;
  };
}
//////////////////////////////////////////////////////////////////////////////
Slash::Data::IIterator* G4Lab::TrajectoryAccessor::iterator(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  // Set properties here since a user G4VTrajectory may
  // had been declared very lately.
  clearProperties();
  // G4VTrajectory properties :
  addProperty("particle",Lib::Property::STRING);
  addProperty("track",Lib::Property::INTEGER);
  addProperty("parent",Lib::Property::INTEGER);
  addProperty("charge",Lib::Property::DOUBLE);
  addProperty("pdg",Lib::Property::INTEGER);
  addProperty("px",Lib::Property::DOUBLE);
  addProperty("py",Lib::Property::DOUBLE);
  addProperty("pz",Lib::Property::DOUBLE);

  fIGeant4Trajectory = isIGeant4Trajectory();
  if(fIGeant4Trajectory) {
    // G4Lab::Trajectory properties :
    addProperty("energy",Lib::Property::DOUBLE);
    addProperty("totalEnergy",Lib::Property::DOUBLE);
    addProperty("globalTime",Lib::Property::DOUBLE);
    addProperty("process",Lib::Property::STRING);
    addProperty("stopX",Lib::Property::DOUBLE);
    addProperty("stopY",Lib::Property::DOUBLE);
    addProperty("stopZ",Lib::Property::DOUBLE);
    addProperty("stopPV",Lib::Property::STRING);
  }
  //addProperty("id",Lib::Property::POINTER);
  addAttDefsProperties();

  if(!fRunManager) {
    Lib::Out out(fSession.printer());
    out << "No G4RunManager." << Lib::endl;
    return 0;
  }
  const G4Event* event = fRunManager->GetCurrentEvent();
  if(!event) {
    Lib::Out out(fSession.printer());
    out << "No event." << Lib::endl;
    return 0;
  }
  G4TrajectoryContainer* trajectoryContainer = 
    event->GetTrajectoryContainer();
  if(!trajectoryContainer) {
    Lib::Out out(fSession.printer());
    out << "No trajectory container." << Lib::endl;
    return 0;
  }

  return new TrajectoryIterator(*trajectoryContainer);
}
//////////////////////////////////////////////////////////////////////////////
Slash::Core::IValue* G4Lab::TrajectoryAccessor::findValue(
 Slash::Data::IAccessor::Data aData
,const std::string& aName
,void* aTag
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  G4VTrajectory* obj = (G4VTrajectory*)aData;
  IGeant4Trajectory* iobj = dynamic_cast<IGeant4Trajectory*>(obj);

  //if(aName=="id") {
    //return new Lib::Value((void*)obj);
  //} else 
  // G4VTrajectory data :
  if(aName=="particle") {
    return new Lib::Value(obj->GetParticleName());
  } else if(aName=="track") {
    return new Lib::Value(obj->GetTrackID());
  } else if(aName=="parent") {
    return new Lib::Value(obj->GetParentID());
  } else if(aName=="charge") {
    return new Lib::Value(obj->GetCharge());
  } else if(aName=="pdg") {
    return new Lib::Value(obj->GetPDGEncoding());
  } else if(aName=="px") {
    return new Lib::Value(obj->GetInitialMomentum().x());
  } else if(aName=="py") {
    return new Lib::Value(obj->GetInitialMomentum().y());
  } else if(aName=="pz") {
    return new Lib::Value(obj->GetInitialMomentum().z());
  } else if(iobj) {
    if(aName=="energy") {
      return new Lib::Value(iobj->kineticEnergy());
    } else if(aName=="totalEnergy") {
      return new Lib::Value(iobj->totalEnergy());
    } else if(aName=="globalTime") {
      return new Lib::Value(iobj->globalTime());
    } else if(aName=="process") {
      return new Lib::Value(iobj->creatorProcessName());
    } else if(aName=="stopX") {
      return new Lib::Value(iobj->stoppingPoint()[0]);
    } else if(aName=="stopY") {
      return new Lib::Value(iobj->stoppingPoint()[1]);
    } else if(aName=="stopZ") {
      return new Lib::Value(iobj->stoppingPoint()[2]);
    } else if(aName=="stopPV") {
      G4VPhysicalVolume* pv = iobj->stoppingPhysicalVolume();
      if(!pv) return new Lib::Value();
      return new Lib::Value(pv->GetName());
    }
  }
  return findAttDefsValue(aData,aName,aTag);
}
//////////////////////////////////////////////////////////////////////////////
void G4Lab::TrajectoryAccessor::beginVisualize(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  OnX::InventorAccessor::beginVisualize();

  // modeling.trajectories:immediate_all
  // modeling.trajectories:immediate_time
  // modeling.trajectories:pickable
  if(!fSession.parameterValue("modeling.trajectories",fModeling))
    fModeling = "immediate_all";

  //fSession.out().println("debug : modeling \"%s\".",fModeling.c_str());

  if(fModeling=="pickable") {

    fSeparator = new SoSeparator;
    fCoordinate = new SoCoordinate3;
    fSeparator->addChild(fCoordinate);
    fIndex = 0;

  } else if( (fModeling=="immediate_all") || (fModeling=="immediate_time") ) {
    // Scene graph :
    SoSeparator* separator = new SoSeparator;
    separator->setName("sceneGraph");
    
    separator->addChild(fSoGC.getHighlightMaterial());
    separator->addChild(fSoGC.getDrawStyle());
    separator->addChild(fSoGC.getLightModel());

    SoG4Trajectories* soG4Trajectories = new SoG4Trajectories(fRunManager);
    separator->addChild(soG4Trajectories);

    if(fModeling=="immediate_time") {
      soG4Trajectories->model.setValue(SoG4Trajectories::TIMED);
      //soG4Trajectories->timeSteps.setValue(1000); 

      SoCalculator* calculator = new SoCalculator();
      calculator->a.connectFrom(&soG4Trajectories->timeSteps);
      calculator->expression.set1Value(0,"oa=2.0/a");

      SoTimeCounter* timeCounter = new SoTimeCounter;
#ifdef WIN32
#undef max
#endif
      timeCounter->max.connectFrom(&soG4Trajectories->timeSteps);
      timeCounter->frequency.connectFrom(&calculator->oa);
      timeCounter->reset.connectFrom(&soG4Trajectories->timeIndex);

      soG4Trajectories->timeIndex.connectFrom(&timeCounter->output);

    }

    fSoRegion->doIt(SbAddNode(separator,"dynamicScene"));
  }
}
//////////////////////////////////////////////////////////////////////////////
void G4Lab::TrajectoryAccessor::visualize(
 Slash::Data::IAccessor::Data aData
,void*
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(!fSeparator) return; // not in "pickable" mode.

  G4VTrajectory* obj = (G4VTrajectory*)aData;
  if(!obj) return;

  int pointn = obj->GetPointEntries();
  if(pointn<=0) return;
  
  SbVec3f* points = new SbVec3f[pointn];
  int32_t* coordIndex = new int32_t[pointn+1];

  for (int i = 0; i < pointn ; i++) {
    G4VTrajectoryPoint* tp = obj->GetPoint(i); //?
    G4ThreeVector pos = tp->GetPosition();
    points[i].setValue((float)pos.x(),(float)pos.y(),(float)pos.z());
    coordIndex[i] = fIndex+i;
  }
  coordIndex[pointn] = SO_END_LINE_INDEX;
  
  //Get a style corresponding to Trajectory type.
 {std::string s = obj->GetParticleName();
  Lib::smanip::strip(s);
  std::string style = "Trajectory("+s+")";
  if(!isStyle(style)) {
    style = "Trajectory_"+s; //Backward compatibility.
    if(!isStyle(style)) {
      style = "Trajectory"; //Default.
    }
  }
  fillSoGC(style);}

  // Build name (for picking) :
  std::string s;
  Lib::smanip::printf(s,128,"Trajectory/0x%lx",(unsigned long)obj);
  SbName name(s.c_str());

  // Scene graph :
  SoSeparator* separator = new SoSeparator;
  separator->setName("sceneGraph");
  fSeparator->addChild(separator);
  
  separator->addChild(fSoGC.getHighlightMaterial());
  separator->addChild(fSoGC.getDrawStyle());
  separator->addChild(fSoGC.getLightModel());

  fCoordinate->point.setValues(fIndex,pointn,points);
  fIndex += pointn;

  SoIndexedLineSet* lineSet = new SoIndexedLineSet;
  lineSet->coordIndex.setValues(0,pointn+1,coordIndex);
  lineSet->setName(name);
  separator->addChild(lineSet);

  delete [] coordIndex;
  delete [] points;
}
//////////////////////////////////////////////////////////////////////////////
void G4Lab::TrajectoryAccessor::endVisualize(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(fSeparator) {
    if(fSeparator->getNumChildren()==1) {
      fSeparator->unref();
    } else {
      fSoRegion->doIt(SbAddNode(fSeparator,"dynamicScene"));
    }
    fSeparator = 0;
    fCoordinate = 0;
    fIndex = 0;
  }
}
//////////////////////////////////////////////////////////////////////////////
bool G4Lab::TrajectoryAccessor::isIGeant4Trajectory(
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(!fRunManager) return false;
  const G4Event* event = fRunManager->GetCurrentEvent();
  if(!event) return false;
  G4TrajectoryContainer* trajectoryContainer = 
    event->GetTrajectoryContainer();
  if(!trajectoryContainer) return false;
  int number = trajectoryContainer->entries();
  if(number<=0) return false;
  G4VTrajectory* vt = (*trajectoryContainer)[(size_t)0];
  IGeant4Trajectory* trajectory = dynamic_cast<IGeant4Trajectory*>(vt);
  return (trajectory?true:false);
}
//////////////////////////////////////////////////////////////////////////////
void G4Lab::TrajectoryAccessor::addAttDefsProperties(
)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  fAttDefs = 0;
  if(!fRunManager) return;
  const G4Event* event = fRunManager->GetCurrentEvent();
  if(!event) return;
  G4TrajectoryContainer* trajectoryContainer = 
    event->GetTrajectoryContainer();
  if(!trajectoryContainer) return;
  int number = trajectoryContainer->entries();
  if(number<=0) return;
  G4VTrajectory* vt = (*trajectoryContainer)[(size_t)0];

  // AttDefs properties :
  fAttDefs = (std::map<G4String,G4AttDef>*)vt->GetAttDefs();
  if(!fAttDefs) return;

  std::map<G4String,G4AttDef>::const_iterator it;
  for(it=fAttDefs->begin();it!=fAttDefs->end();++it) {
    const std::string& sname = (*it).second.GetName();
    const std::string& stype = (*it).second.GetValueType();    

    if(sname=="Color") continue; //Do not put some AttDef as property.

    //printf("debug : add %s %s\n",sname.c_str(),stype.c_str());

    if( (stype=="int")   ||
        (stype=="G4int") ) {
      addProperty((*it).first,Lib::Property::INTEGER);
    } else if( (stype=="double")   ||
               (stype=="float")    ||
               (stype=="G4double") ||
               (stype=="G4float")  ) {
      addProperty((*it).first,Lib::Property::DOUBLE);
    } else if( (stype=="std::string")   ||
               (stype=="G4String")  ) {
      addProperty((*it).first,Lib::Property::STRING);
    } else if( (stype=="std::vector<double>")   ||
               (stype=="std::vector<float>")    ||
               (stype=="G4ThreeVector")         ||
               (stype=="std::vector<G4double>") ||
               (stype=="std::vector<G4float>")  ) {
      addProperty((*it).first,Lib::Property::VECTOR_DOUBLE);
    } else {
      addProperty((*it).first,Lib::Property::STRING);
    }
  }
}
//////////////////////////////////////////////////////////////////////////////
Slash::Core::IValue* G4Lab::TrajectoryAccessor::findAttDefsValue(
 Slash::Data::IAccessor::Data aData
,const std::string& aName
,void*
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  //printf("debug : find %s (%ld)...\n",aName.c_str(),fAttDefs);

  if(!fAttDefs) return new Lib::Value();

  G4VTrajectory* obj = (G4VTrajectory*)aData;
  std::vector<G4AttValue>* vec = obj->CreateAttValues();
  if(!vec) {
    Lib::Out out(printer());
    out << "G4Lab::TrajectoryAccessor::findAttDefsValue :" 
        << " problem to get AttValues fot trajectory."
        << Lib::endl;
    return new Lib::Value();
  }
  
  unsigned int number = vec->size();
  for(unsigned int index=0;index<number;index++) {
    const G4AttValue& val = (*vec)[index];
    if(aName==val.GetName()) {
      const std::string& stype = (*fAttDefs)[val.GetName()].GetValueType();
      if( (stype=="int")   ||
          (stype=="G4int") ) {
        int v;
        if(!Lib::smanip::toint(val.GetValue(),v)) {
          Lib::Out out(printer());
          out << "G4Lab::TrajectoryAccessor::findAttDefsValue :" 
              << " " << Lib::sout(val.GetValue()) << " not a G4int." 
              << Lib::endl;
          delete vec;
          return new Lib::Value();
        }
        delete vec;
        return new Lib::Value(v);
      } else if( (stype=="double")   ||
                 (stype=="float")    ||
                 (stype=="G4double") ||
                 (stype=="G4float")  ) {
        double v = 0;
        if(!Lib::smanip::todouble(val.GetValue(),v)) {
          std::vector<double> vs;
          if(!BestUnit::toDoubles(printer(),val.GetValue(),vs)) {
            delete vec;
            return new Lib::Value();
          }
          if(vs.size()!=1) {
            Lib::Out out(printer());
            out << "G4Lab::TrajectoryAccessor::findAttDefsValue :" 
                << " " << Lib::sout(val.GetValue()) << " not a double." 
                << Lib::endl;
            delete vec;
            return new Lib::Value();
          }
          v = vs[0];
        }
        delete vec;
        return new Lib::Value(v);
      } else if( (stype=="std::string")   ||
                 (stype=="G4String")  ) {
        delete vec;
        return new Lib::Value(val.GetValue());
      } else if(stype=="G4ThreeVector") {
        std::vector<double> vs;
        if(!BestUnit::toDoubles(printer(),val.GetValue(),vs)) {
          delete vec;
          return new Lib::Value();
        }
        if(vs.size()!=3) {
          Lib::Out out(printer());
          out << "G4Lab::TrajectoryAccessor::findAttDefsValue :" 
              << " " << Lib::sout(val.GetValue()) << " not a vector." 
              << Lib::endl;
          delete vec;
          return new Lib::Value();
        }
        delete vec;
        return new Lib::Value(vs);
      } else if( (stype=="std::vector<float>")    ||
                 (stype=="std::vector<G4float>")  ) {
        void* p;
        if(!Lib::smanip::topointer(val.GetValue(),p)) {
          Lib::Out out(printer());
          out << "G4Lab::TrajectoryAccessor::findAttDefsValue :" 
              << " " << Lib::sout(val.GetValue()) << " not a pointer." 
              << Lib::endl;
          delete vec;
          return new Lib::Value();
        }
        std::vector<float>* vvec = (std::vector<float>*)p;
        unsigned int vn = vvec->size();
        std::vector<double> array(vn);
        for(unsigned int vi=0;vi<vn;vi++) array[vi] = double((*vvec)[vi]);
        delete vec;
        return new Lib::Value(array);
      } else if( (stype=="std::vector<double>")   ||
                 (stype=="std::vector<G4double>") ){
        void* p;
        if(!Lib::smanip::topointer(val.GetValue(),p)) {
          Lib::Out out(printer());
          out << "G4Lab::TrajectoryAccessor::findAttDefsValue :" 
              << " " << Lib::sout(val.GetValue()) << " not a pointer." 
              << Lib::endl;
          delete vec;
          return new Lib::Value();
        }
        std::vector<double>* vvec = (std::vector<double>*)p;
        delete vec;
        return new Lib::Value(*vvec);
      } else {
        delete vec;
        return new Lib::Value(val.GetValue());
      }
    }
  }
  Lib::Out out(printer());
  out << "G4Lab::TrajectoryAccessor::findAttDefsValue :" 
      << " AttValue not found for property " << Lib::sout(aName) << "." 
      << Lib::endl;
  delete vec;
  return new Lib::Value();
}