source: snovis/trunk/source/G4Lab/cxx/PhysicsTableAccessor.cxx @ 288

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

// Inventor :
#include <Inventor/SbString.h>
#include <Inventor/SbName.h>
#include <Inventor/nodes/SoSeparator.h>

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

// Geant4 :
#include <G4ParticleTable.hh>
#include <G4PhysicsTable.hh>
#include <G4ParticleDefinition.hh>
#include <G4Element.hh>
#include <G4VProcess.hh>
#include <G4ProcessVector.hh>
#include <G4ProcessManager.hh>

// Lib :
#include <Slash/Core/ISession.h>
#include <Slash/Data/IIterator.h>
#include <Lib/Out.h>
#include <Lib/sout.h>
#include <Lib/Value.h>
#include <Lib/smanip.h>
#include <Lib/fmanip.h>
#include <Lib/dirmanip.h>
#include <Lib/Cast.h>
#include <Lib/Debug.h>

// AIDA :
#include <AIDA/IAnalysisFactory.h>
#include <AIDA/ITreeFactory.h>
#include <AIDA/ITree.h>
#include <AIDA/IPlotterFactory.h>
#include <AIDA/IPlotter.h>
#include <AIDA/IPlotterRegion.h>
#include <AIDA/IFunctionFactory.h>
#include <AIDA/IFunction.h>

// G4Lab :
#include <G4Lab/Interfaces/IGeant4Manager.h>

namespace G4Lab {
class PhysicsTable {
public:
  PhysicsTable(const std::string& aTableName,
               const std::string& aProcessName,
               const std::string& aParticleName,
               G4PhysicsTable* aTable,
               G4VProcess* aProcess,
               G4ParticleDefinition* aParticle)
    :fTableName(aTableName)
    ,fProcessName(aProcessName)
    ,fParticleName(aParticleName)
    ,fTable(aTable)
    ,fProcess(aProcess)
    ,fParticle(aParticle){}
  virtual ~PhysicsTable() {
    delete fTable;
  }
public:
  std::string fTableName;
  std::string fProcessName;
  std::string fParticleName;
  G4PhysicsTable* fTable;
  G4VProcess* fProcess;
  G4ParticleDefinition* fParticle;
  Lib::Debug fDebug;
};
}

//////////////////////////////////////////////////////////////////////////////
G4Lab::PhysicsTableAccessor::PhysicsTableAccessor(
 Slash::Core::ISession& aSession
,IGeant4Manager& aManager
,AIDA::IAnalysisFactory* aAIDA
)
:Lib::BaseAccessor(aSession.printer())
,fSession(aSession)
,fManager(aManager)
,fType("PhysicsTable")
,fAIDA(aAIDA)
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  addProperty("table",Lib::Property::STRING);
  addProperty("process",Lib::Property::STRING);
  addProperty("particle",Lib::Property::STRING);
  addProperty("number",Lib::Property::INTEGER);
  addProperty("name",Lib::Property::STRING);
}
//////////////////////////////////////////////////////////////////////////////
G4Lab::PhysicsTableAccessor::~PhysicsTableAccessor(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  for(unsigned int index=0;index<fTables.size();index++) delete fTables[index];
}
//////////////////////////////////////////////////////////////////////////////
void* G4Lab::PhysicsTableAccessor::cast(
 const std::string& aClass
) const 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if_Lib_SCast(Slash::Data::IVisualizer)
  else return Lib::BaseAccessor::cast(aClass);
}
//////////////////////////////////////////////////////////////////////////////
std::string G4Lab::PhysicsTableAccessor::name(
) const
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  return fType;
}
namespace G4Lab {
  class PhysicsTableIterator : 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:
    PhysicsTableIterator(std::vector<G4Lab::PhysicsTable*>& aVector)
      :fVector(aVector) {
      fIterator = fVector.begin();
    }
    virtual ~PhysicsTableIterator() {}
  private:
    std::vector<G4Lab::PhysicsTable*>& fVector;
    std::vector<G4Lab::PhysicsTable*>::iterator fIterator;
  };
}
//////////////////////////////////////////////////////////////////////////////
Slash::Data::IIterator* G4Lab::PhysicsTableAccessor::iterator(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(!fTables.size()) {
    if(!buildProcessTables()) return 0;
  }
  return new PhysicsTableIterator(fTables);
}
//////////////////////////////////////////////////////////////////////////////
Slash::Core::IValue* G4Lab::PhysicsTableAccessor::findValue(
 Slash::Data::IAccessor::Data aData
,const std::string& aName
,void*
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  PhysicsTable* obj = (PhysicsTable*)aData;
  if(aName=="table") {
    return new Lib::Value(obj->fTableName);
  } else if(aName=="process") {
    return new Lib::Value(obj->fProcessName);
  } else if(aName=="particle") {
    return new Lib::Value(obj->fParticleName);
  } else if(aName=="number") {
    return new Lib::Value((int)obj->fTable->size());
  } else if(aName=="name") {
    std::string name = 
      obj->fTableName + "." + obj->fProcessName + "." + obj->fParticleName;
    return new Lib::Value(name);
  } else {
    return new Lib::Value();
  }
}
//////////////////////////////////////////////////////////////////////////////
void G4Lab::PhysicsTableAccessor::beginVisualize(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(!fAIDA) {
    Lib::Out out(fSession.printer());
    out << "G4Lab::PhysicsTableAccessor::beginVisualize :"
        << " AIDA not found. Can't then find an AIDA plotter." 
        << Lib::endl;
  }
}
//////////////////////////////////////////////////////////////////////////////
void G4Lab::PhysicsTableAccessor::visualize(
 Slash::Data::IAccessor::Data aData
,void*
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  PhysicsTable* obj = (PhysicsTable*)aData;

  std::string value;

  bool superpose = true;
  if(fSession.parameterValue("modeling.superpose",value))
    if(!Lib::smanip::tobool(value,superpose)) superpose = true;

  if(!fAIDA) return;
    
  AIDA::ITreeFactory* tf = fAIDA->createTreeFactory();
  if(!tf) return;
  AIDA::ITree* memoryTree = tf->create();
  delete tf;
  if(!memoryTree) return;

  AIDA::IFunctionFactory* functionFactory = 
    fAIDA->createFunctionFactory(*memoryTree);
  if(!functionFactory) {
    Lib::Out out(fSession.printer());
    out << "Can't create an function factory." << Lib::endl;
    return;
  }

  AIDA::IPlotterFactory* plotterFactory = fAIDA->createPlotterFactory();
  if(!plotterFactory) {
    Lib::Out out(fSession.printer());
    out << "Can't create a plotter factory." << Lib::endl;
    delete functionFactory;
    return;
  }

  AIDA::IPlotter* plotter = plotterFactory->create();
  if(!plotter) {
    Lib::Out out(fSession.printer());
    out << "Can't create a plotter factory." << Lib::endl;
    delete functionFactory;
    delete plotterFactory;
    return;
  }

  int number = obj->fTable->entries();
  if(!superpose) {
    if(number>0) plotter->createRegions(1,number,0);
    else plotter->createRegions(1,1,0);
  }

  for (int count=0;count<number;count++){
    G4PhysicsVector* physicsVector = (*(obj->fTable))[count];
    size_t n = physicsVector->GetVectorLength();
    if(n<2) {
      Lib::Out out(fSession.printer());
      out << "PhysicsVector with " << int(n) << " entries. Can't plot." 
          << Lib::endl;
      continue;
    }

    std::string name;
    Lib::smanip::printf(name,1024,"%s.%s.%s.%d",
                  obj->fTableName.c_str(),
                  obj->fProcessName.c_str(),
                  obj->fParticleName.c_str(),
                  count);
    std::vector<double> params;
    params.push_back(1); // Dimension of the grid (here one).
    params.push_back(n); // Number of entries for first grid axis.
    size_t i;
    for (i=0;i<n;i++) params.push_back(i); // Xs.
    for (i=0;i<n;i++) params.push_back((*physicsVector)[i]); //Values.
    AIDA::IFunction* function = 
      functionFactory->createFunctionByName(name,"Grid1D");
    if(!function) {
      Lib::Out out(fSession.printer());
      out << "Can't create \"" << name << "\" function." << Lib::endl;
      break;
    }
    function->setParameters(params);
      
    AIDA::IPlotterRegion& region = plotter->currentRegion();

    region.plot(*function);
    region.setParameter("plotter.xAxisAutomated","FALSE");
    std::string value1;
    Lib::smanip::printf(value1,32,"%g",0.);
    region.setParameter("plotter.xAxisMinimum",value1);
    std::string value2;
    Lib::smanip::printf(value2,32,"%g",(double)(n-1));
    region.setParameter("plotter.xAxisMaximum",value2);
    //std::string value3;
    //Lib::smanip::printf(value3,32,"%d",4 * n);
    //region.setParameter("plotter.xNumberOfPoints",value3);
    /*
    fSession.out().println("debug : min : %s, max : %s, steps : %s\n.",
                           value1.c_str(),value2.c_str(),value3.c_str());
    */
    if(!superpose) plotter->next();
      
  }
  
  delete plotter;
  delete functionFactory;
  delete plotterFactory;
    
}
//////////////////////////////////////////////////////////////////////////////
void G4Lab::PhysicsTableAccessor::endVisualize(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
}
//////////////////////////////////////////////////////////////////////////////
bool G4Lab::PhysicsTableAccessor::buildProcessTables(
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  if(fTables.size()) return true; //done.

  if(!fManager.isRunning()) {
    Lib::Out out(fSession.printer());
    out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
        << " It is needed to have started a run to get physics tables." 
        << Lib::endl;
    return false;
  }

  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  if(!particleTable) {
    Lib::Out out(fSession.printer());
    out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
        << " Can't get particle table." 
        << Lib::endl;
    return false;
  }

  G4ParticleTable::G4PTblDicIterator* particleTableIterator = 
    particleTable->GetIterator();
  if(!particleTableIterator) {
    Lib::Out out(fSession.printer());
    out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
        << " Can't get particle table iterator." 
        << Lib::endl;
    return false;
  }

  particleTableIterator->reset();
  while( (*particleTableIterator)() ){
    G4ParticleDefinition* particle = particleTableIterator->value();
    if(!particle) {
      Lib::Out out(fSession.printer());
      out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
          << " Can't get particle table particle." 
          << Lib::endl;
      return false;
    }
    G4String particleName = particle->GetParticleName();
    //Lib::Out out(fSession.printer());
    //out << "debug : Particle : " << particleName << Lib::endl;

    // Processes for this particle :
    G4ProcessManager* processManager = particle->GetProcessManager();
    if(!processManager) {
      Lib::Out out(fSession.printer());
      out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
          << " Can't get process manager of particle " 
          << particleName << "."
          << Lib::endl;
      return false;
    }

    G4ProcessVector* processList = processManager->GetProcessList();
    if(!processList) {
      Lib::Out out(fSession.printer());
      out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
          << " Can't get process list for process manager of particle " 
          << particleName << "."
          << Lib::endl;
      return false;
    }

    int number = processList->entries();
    for (int index=0;index<number;index++){
      G4VProcess* process = (*processList)(index);
      std::vector<std::string> names;
      std::vector<G4PhysicsTable*> tables;
      if(!findProcessTables(particle,process,names,tables)) {
        return false;
      }

      for(unsigned int i=0;i<tables.size();i++) {
        std::vector<std::string> words;
        Lib::smanip::words(names[i],".",words);
        if(words.size()!=3) {
          Lib::Out out(fSession.printer());
          out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
              << " Three words separted with a dot expected in " 
              << names[i] << " for particle " << particleName << "."
              << Lib::endl;
          return false;
        }
        fTables.push_back(new PhysicsTable(words[0],
                                           words[1],
                                           words[2],
                                           tables[i],
                                           process,
                                           particle));
      }
    }
  }

  return true;
}
//////////////////////////////////////////////////////////////////////////////
bool G4Lab::PhysicsTableAccessor::findProcessTables(
 G4ParticleDefinition* aParticle
,G4VProcess* aProcess
,std::vector<std::string>& aTables
,std::vector<G4PhysicsTable*>& aPhysicsTables
) 
//////////////////////////////////////////////////////////////////////////////
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
{
  aTables.clear();
  aPhysicsTables.clear();

  std::string tmpName;
  if(!Lib::fmanip::tmpname(tmpName,".","G4Lab_tmp","")) {
    Lib::Out out(fSession.printer());
    out << "Can't build a temporary directory name." << Lib::endl;
    return false;
  }

  if(!Lib::dirmanip::create(tmpName)) {
    Lib::Out out(fSession.printer());
    out << "Can't create \"" << tmpName << "\" directory." << Lib::endl;
    return false;
  }

#ifdef WIN32 //binary mode out on Windows.
  G4bool ascii = true;
#else
  G4bool ascii = false;
  //G4bool ascii = true;
#endif

  if(!aProcess->StorePhysicsTable(aParticle,tmpName,ascii)) {
    Lib::Out out(fSession.printer());
    out << "Can't store physics table for particle "
        << aParticle->GetParticleName()
        << " and process " << aProcess->GetProcessName() << "."
        << Lib::endl;
    //return false;
  } else {

    // Analyse the content of the directory :
    std::vector<std::string> files;
    if(!Lib::dirmanip::entries(tmpName,files,false)) {
      Lib::Out out(fSession.printer());
      out << "Can't get files in \"" << tmpName << "\" directory." 
          << Lib::endl;
    } else {
      // Analyse files :
      Lib::smanip::remove(files,".");
      Lib::smanip::remove(files,"..");
/*
      if(!files.size()) {
        Lib::Out out(fSession.printer());
        out << "No files in \"" << tmpName << "\" directory for particle " 
            << aParticle->GetParticleName() 
            << " and process " << aProcess->GetProcessName() << "."
            << Lib::endl;
        //return false;
      }
*/
      for(unsigned int index=0;index<files.size();index++) {
        std::vector<std::string> words;
        Lib::smanip::words(files[index],".",words);
        //Lib::Out out(fSession.printer());
        //out << "debug : Analyse \"" << files[index] << "\" ." << Lib::endl;
        if(words.size()!=4) {
          Lib::Out out(fSession.printer());
          out << "Bad syntax in file name " << Lib::sout(files[index])
              << " for particle " << aParticle->GetParticleName() 
              << " and process " << aProcess->GetProcessName() << "."
              << Lib::endl;
          //aTables.clear();
          //aPhysicsTables.clear();
          //return false;
          continue;
        }
        std::string table = words[0];
        std::string process = words[1];
        std::string particle = words[2];
        //Lib::smanip::replace(aTables[i],".dat","");
        std::string filename = 
          aProcess->GetPhysicsTableFileName(aParticle,
                                            tmpName,
                                            table,
                                            ascii);
        G4PhysicsTable* physicsTable = 
          new G4PhysicsTable(G4Element::GetNumberOfElements());
        G4bool status = physicsTable->RetrievePhysicsTable(filename,ascii);
        if(!status){
          Lib::Out out(fSession.printer());
          out << "Unable to read \"" << filename << "\" file." << Lib::endl;
          delete physicsTable;
          physicsTable = 0;
        }
        if(physicsTable) {
          aTables.push_back(table+"."+process+"."+particle);
          aPhysicsTables.push_back(physicsTable);
        }
      }
    }

  }

  // Remove the temporary directory :
  if(!Lib::dirmanip::remove(tmpName)) {
    Lib::Out out(fSession.printer());
    out << "Can't remove \"" << tmpName << "\" directory." << Lib::endl;
    aTables.clear();
    aPhysicsTables.clear();
    return false;
  }

  return true;
}