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

Last change on this file since 233 was 233, checked in by barrand, 17 years ago
  • Property svn:eol-style set to native
File size: 14.1 KB
Line 
1// this :
2#include <G4Lab/PhysicsTableAccessor.h>
3
4// Inventor :
5#include <Inventor/SbString.h>
6#include <Inventor/SbName.h>
7#include <Inventor/nodes/SoSeparator.h>
8
9#ifdef WIN32
10#undef pascal // Clash between windef.h and Geant4/SystemOfUnits.hh
11#endif
12
13// Geant4 :
14#include <G4ParticleTable.hh>
15#include <G4PhysicsTable.hh>
16#include <G4ParticleDefinition.hh>
17#include <G4Element.hh>
18#include <G4VProcess.hh>
19#include <G4ProcessVector.hh>
20#include <G4ProcessManager.hh>
21
22// Lib :
23#include <Slash/Core/ISession.h>
24#include <Slash/Data/IIterator.h>
25#include <Lib/Out.h>
26#include <Lib/Value.h>
27#include <Lib/smanip.h>
28#include <Lib/fmanip.h>
29#include <Lib/dirmanip.h>
30#include <Lib/Debug.h>
31
32// AIDA :
33#include <AIDA/IAnalysisFactory.h>
34#include <AIDA/ITreeFactory.h>
35#include <AIDA/ITree.h>
36#include <AIDA/IPlotterFactory.h>
37#include <AIDA/IPlotter.h>
38#include <AIDA/IPlotterRegion.h>
39#include <AIDA/IFunctionFactory.h>
40#include <AIDA/IFunction.h>
41
42namespace G4Lab {
43class PhysicsTable {
44public:
45  PhysicsTable(const std::string& aTableName,
46               const std::string& aProcessName,
47               const std::string& aParticleName,
48               G4PhysicsTable* aTable,
49               G4VProcess* aProcess,
50               G4ParticleDefinition* aParticle)
51    :fTableName(aTableName)
52    ,fProcessName(aProcessName)
53    ,fParticleName(aParticleName)
54    ,fTable(aTable)
55    ,fProcess(aProcess)
56    ,fParticle(aParticle){}
57  virtual ~PhysicsTable() {
58    delete fTable;
59  }
60public:
61  std::string fTableName;
62  std::string fProcessName;
63  std::string fParticleName;
64  G4PhysicsTable* fTable;
65  G4VProcess* fProcess;
66  G4ParticleDefinition* fParticle;
67  Lib::Debug fDebug;
68};
69}
70
71//////////////////////////////////////////////////////////////////////////////
72G4Lab::PhysicsTableAccessor::PhysicsTableAccessor(
73 Slash::Core::ISession& aSession
74,AIDA::IAnalysisFactory* aAIDA
75)
76:Lib::BaseAccessor(aSession.printer())
77,fSession(aSession)
78,fType("PhysicsTable")
79,fAIDA(aAIDA)
80//////////////////////////////////////////////////////////////////////////////
81//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
82{
83  addProperty("table",Lib::Property::STRING);
84  addProperty("process",Lib::Property::STRING);
85  addProperty("particle",Lib::Property::STRING);
86  addProperty("number",Lib::Property::INTEGER);
87  addProperty("name",Lib::Property::STRING);
88}
89//////////////////////////////////////////////////////////////////////////////
90G4Lab::PhysicsTableAccessor::~PhysicsTableAccessor(
91) 
92//////////////////////////////////////////////////////////////////////////////
93//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
94{
95  for(unsigned int index=0;index<fTables.size();index++) delete fTables[index];
96}
97//////////////////////////////////////////////////////////////////////////////
98std::string G4Lab::PhysicsTableAccessor::name(
99) const
100//////////////////////////////////////////////////////////////////////////////
101//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
102{
103  return fType;
104}
105namespace G4Lab {
106  class PhysicsTableIterator : public virtual Slash::Data::IIterator {
107  public: //Slash::Data::IIterator
108    virtual Slash::Data::IAccessor::Data object() {
109      if(fIterator==fVector.end()) return 0;
110      return *fIterator;
111    }
112    virtual void next() { ++fIterator; }
113    virtual void* tag() { return 0;}
114  public:
115    PhysicsTableIterator(std::vector<G4Lab::PhysicsTable*>& aVector)
116      :fVector(aVector) {
117      fIterator = fVector.begin();
118    }
119    virtual ~PhysicsTableIterator() {}
120  private:
121    std::vector<G4Lab::PhysicsTable*>& fVector;
122    std::vector<G4Lab::PhysicsTable*>::iterator fIterator;
123  };
124}
125//////////////////////////////////////////////////////////////////////////////
126Slash::Data::IIterator* G4Lab::PhysicsTableAccessor::iterator(
127) 
128//////////////////////////////////////////////////////////////////////////////
129//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
130{
131  if(!fTables.size() ) {
132    G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
133    if(particleTable) {
134      G4ParticleTable::G4PTblDicIterator* particleTableIterator =
135        particleTable->GetIterator();
136      if(particleTableIterator) {
137        particleTableIterator->reset();
138        while( (*particleTableIterator)() ){
139          G4ParticleDefinition* particle = particleTableIterator->value();
140          if(particle) {
141            G4String particleName = particle->GetParticleName();
142            //Lib::Out out(fSession.printer());
143            //out << "debug : Particle : " << particleName << Lib::endl;
144            // Processes for this particle :
145            G4ProcessManager* processManager = particle->GetProcessManager();
146            if(processManager) {
147              G4ProcessVector* processList = processManager->GetProcessList();
148              if(processList) {
149                int number = processList->entries();
150                for (int index=0;index<number;index++){
151                  G4VProcess* process = (*processList)(index);
152                  std::vector<std::string> names;
153                  std::vector<G4PhysicsTable*> tables;
154                  if(findProcessTables(particle,
155                                       process,
156                                       names,
157                                       tables)) {
158                    for(unsigned int i=0;i<tables.size();i++) {
159                      std::vector<std::string> words;
160                      Lib::smanip::words(names[i],".",words);
161                      if(words.size()==3) {
162                        fTables.push_back(new PhysicsTable(words[0],
163                                                           words[1],
164                                                           words[2],
165                                                           tables[i],
166                                                           process,
167                                                           particle));
168                      }
169                    }
170                  }
171                }
172              }
173            }
174          }
175        }
176      }
177    }
178  }
179
180  return new PhysicsTableIterator(fTables);
181}
182//////////////////////////////////////////////////////////////////////////////
183Slash::Core::IValue* G4Lab::PhysicsTableAccessor::findValue(
184 Slash::Data::IAccessor::Data aData
185,const std::string& aName
186,void*
187) 
188//////////////////////////////////////////////////////////////////////////////
189//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
190{
191  PhysicsTable* obj = (PhysicsTable*)aData;
192  if(aName=="table") {
193    return new Lib::Value(obj->fTableName);
194  } else if(aName=="process") {
195    return new Lib::Value(obj->fProcessName);
196  } else if(aName=="particle") {
197    return new Lib::Value(obj->fParticleName);
198  } else if(aName=="number") {
199    return new Lib::Value((int)obj->fTable->size());
200  } else if(aName=="name") {
201    std::string name = 
202      obj->fTableName + "." + obj->fProcessName + "." + obj->fParticleName;
203    return new Lib::Value(name);
204  } else {
205    return new Lib::Value();
206  }
207}
208//////////////////////////////////////////////////////////////////////////////
209void G4Lab::PhysicsTableAccessor::beginVisualize(
210) 
211//////////////////////////////////////////////////////////////////////////////
212//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
213{
214  if(!fAIDA) {
215    Lib::Out out(fSession.printer());
216    out << "G4Lab::PhysicsTableAccessor::beginVisualize :"
217        << " AIDA not found. Can't then find an AIDA plotter." 
218        << Lib::endl;
219  }
220}
221//////////////////////////////////////////////////////////////////////////////
222void G4Lab::PhysicsTableAccessor::visualize(
223 Slash::Data::IAccessor::Data aData
224,void*
225) 
226//////////////////////////////////////////////////////////////////////////////
227//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
228{
229  PhysicsTable* obj = (PhysicsTable*)aData;
230
231  std::string value;
232
233  bool superpose = true;
234  if(fSession.parameterValue("modeling.superpose",value))
235    if(!Lib::smanip::tobool(value,superpose)) superpose = true;
236
237  if(!fAIDA) return;
238   
239  AIDA::ITreeFactory* tf = fAIDA->createTreeFactory();
240  if(!tf) return;
241  AIDA::ITree* memoryTree = tf->create();
242  delete tf;
243  if(!memoryTree) return;
244
245  AIDA::IFunctionFactory* functionFactory = 
246    fAIDA->createFunctionFactory(*memoryTree);
247  if(!functionFactory) {
248    Lib::Out out(fSession.printer());
249    out << "Can't create an function factory." << Lib::endl;
250    return;
251  }
252
253  AIDA::IPlotterFactory* plotterFactory = fAIDA->createPlotterFactory();
254  if(!plotterFactory) {
255    Lib::Out out(fSession.printer());
256    out << "Can't create a plotter factory." << Lib::endl;
257    delete functionFactory;
258    return;
259  }
260
261  AIDA::IPlotter* plotter = plotterFactory->create();
262  if(!plotter) {
263    Lib::Out out(fSession.printer());
264    out << "Can't create a plotter factory." << Lib::endl;
265    delete functionFactory;
266    delete plotterFactory;
267    return;
268  }
269
270  int number = obj->fTable->entries();
271  if(!superpose) {
272    if(number>0) plotter->createRegions(1,number,0);
273    else plotter->createRegions(1,1,0);
274  }
275
276  for (int count=0;count<number;count++){
277    G4PhysicsVector* physicsVector = (*(obj->fTable))[count];
278    size_t n = physicsVector->GetVectorLength();
279    if(n>=2) {
280      std::string name;
281      Lib::smanip::printf(name,1024,"%s.%s.%s.%d",
282                    obj->fTableName.c_str(),
283                    obj->fProcessName.c_str(),
284                    obj->fParticleName.c_str(),
285                    count);
286      std::vector<double> params;
287      params.push_back(1); // Dimension of the grid (here one).
288      params.push_back(n); // Number of entries for first grid axis.
289      size_t i;
290      for (i=0;i<n;i++) params.push_back(i); // Xs.
291      for (i=0;i<n;i++) params.push_back((*physicsVector)[i]); //Values.
292      AIDA::IFunction* function = 
293        functionFactory->createFunctionByName(name,"Grid1D");
294      if(!function) {
295        Lib::Out out(fSession.printer());
296        out << "Can't create \"" << name << "\" function." << Lib::endl;
297        break;
298      }
299      function->setParameters(params);
300     
301      AIDA::IPlotterRegion& region = plotter->currentRegion();
302
303      region.plot(*function);
304      region.setParameter("plotter.xAxisAutomated","FALSE");
305      std::string value1;
306      Lib::smanip::printf(value1,32,"%g",0.);
307      region.setParameter("plotter.xAxisMinimum",value1);
308      std::string value2;
309      Lib::smanip::printf(value2,32,"%g",(double)(n-1));
310      region.setParameter("plotter.xAxisMaximum",value2);
311      std::string value3;
312      Lib::smanip::printf(value3,32,"%d",4 * n);
313      region.setParameter("plotter.xNumberOfPoints",value3);
314      /*
315      fSession.out().println("debug : min : %s, max : %s, steps : %s\n.",
316                             value1.c_str(),value2.c_str(),value3.c_str());
317      */
318
319      if(!superpose) plotter->next();
320     
321    }
322  }
323 
324  delete plotter;
325  delete functionFactory;
326  delete plotterFactory;
327   
328}
329//////////////////////////////////////////////////////////////////////////////
330void G4Lab::PhysicsTableAccessor::endVisualize(
331) 
332//////////////////////////////////////////////////////////////////////////////
333//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
334{
335}
336//////////////////////////////////////////////////////////////////////////////
337bool G4Lab::PhysicsTableAccessor::findProcessTables(
338 G4ParticleDefinition* aParticle
339,G4VProcess* aProcess
340,std::vector<std::string>& aTables
341,std::vector<G4PhysicsTable*>& aPhysicsTables
342) 
343//////////////////////////////////////////////////////////////////////////////
344//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
345{
346  aTables.clear();
347  aPhysicsTables.clear();
348
349  std::string tmpName;
350  if(!Lib::fmanip::tmpname(tmpName,".","G4Lab_tmp","")) {
351    Lib::Out out(fSession.printer());
352    out << "Can't build a temporary directory name." << Lib::endl;
353    return false;
354  }
355
356  if(!Lib::dirmanip::create(tmpName)) {
357    Lib::Out out(fSession.printer());
358    out << "Can't create \"" << tmpName << "\" directory." << Lib::endl;
359    return false;
360  }
361
362#ifdef WIN32 //binary mode out on Windows.
363  G4bool ascii = true;
364#else
365  G4bool ascii = false;
366#endif
367  aProcess->StorePhysicsTable(aParticle,tmpName,ascii);
368
369  // Analyse the content of the directory :
370  std::vector<std::string> files;
371  if(!Lib::dirmanip::entries(tmpName,files,false)) {
372    Lib::Out out(fSession.printer());
373    out << "Can't get files in \"" << tmpName << "\" directory." << Lib::endl;
374  } else {
375    // Analyse files :
376    Lib::smanip::remove(files,".");
377    Lib::smanip::remove(files,"..");
378    for(unsigned int index=0;index<files.size();index++) {
379      std::vector<std::string> words;
380      Lib::smanip::words(files[index],".",words);
381      //Lib::Out out(fSession.printer());
382      //out << "debug : Analyse \"" << files[index] << "\" ." << Lib::endl;
383      if(words.size()==4) {
384        std::string table = words[0];
385        std::string process = words[1];
386        std::string particle = words[2];
387        //Lib::smanip::replace(aTables[i],".dat","");
388        std::string filename = 
389          aProcess->GetPhysicsTableFileName(aParticle,
390                                            tmpName,
391                                            table,
392                                            ascii);
393        G4PhysicsTable* physicsTable = 
394          new G4PhysicsTable(G4Element::GetNumberOfElements());
395        if(physicsTable) {
396          G4bool status = physicsTable->RetrievePhysicsTable(filename,ascii);
397          if(!status){
398            Lib::Out out(fSession.printer());
399            out << "Unable to read \"" << filename << "\" file." << Lib::endl;
400            delete physicsTable;
401            physicsTable = 0;
402          }
403        }
404        if(physicsTable) {
405          aTables.push_back(table+"."+process+"."+particle);
406          aPhysicsTables.push_back(physicsTable);
407        }
408      }
409    }
410  }
411
412  // Remove the temporary directory :
413  if(!Lib::dirmanip::remove(tmpName)) {
414    Lib::Out out(fSession.printer());
415    out << "Can't remove \"" << tmpName << "\" directory." << Lib::endl;
416    aTables.clear();
417    aPhysicsTables.clear();
418    return false;
419  }
420
421  return true;
422}
Note: See TracBrowser for help on using the repository browser.