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

Last change on this file since 288 was 288, checked in by barrand, 18 years ago
  • Property svn:eol-style set to native
File size: 17.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/sout.h>
27#include <Lib/Value.h>
28#include <Lib/smanip.h>
29#include <Lib/fmanip.h>
30#include <Lib/dirmanip.h>
31#include <Lib/Cast.h>
32#include <Lib/Debug.h>
33
34// AIDA :
35#include <AIDA/IAnalysisFactory.h>
36#include <AIDA/ITreeFactory.h>
37#include <AIDA/ITree.h>
38#include <AIDA/IPlotterFactory.h>
39#include <AIDA/IPlotter.h>
40#include <AIDA/IPlotterRegion.h>
41#include <AIDA/IFunctionFactory.h>
42#include <AIDA/IFunction.h>
43
44// G4Lab :
45#include <G4Lab/Interfaces/IGeant4Manager.h>
46
47namespace G4Lab {
48class PhysicsTable {
49public:
50 PhysicsTable(const std::string& aTableName,
51 const std::string& aProcessName,
52 const std::string& aParticleName,
53 G4PhysicsTable* aTable,
54 G4VProcess* aProcess,
55 G4ParticleDefinition* aParticle)
56 :fTableName(aTableName)
57 ,fProcessName(aProcessName)
58 ,fParticleName(aParticleName)
59 ,fTable(aTable)
60 ,fProcess(aProcess)
61 ,fParticle(aParticle){}
62 virtual ~PhysicsTable() {
63 delete fTable;
64 }
65public:
66 std::string fTableName;
67 std::string fProcessName;
68 std::string fParticleName;
69 G4PhysicsTable* fTable;
70 G4VProcess* fProcess;
71 G4ParticleDefinition* fParticle;
72 Lib::Debug fDebug;
73};
74}
75
76//////////////////////////////////////////////////////////////////////////////
77G4Lab::PhysicsTableAccessor::PhysicsTableAccessor(
78 Slash::Core::ISession& aSession
79,IGeant4Manager& aManager
80,AIDA::IAnalysisFactory* aAIDA
81)
82:Lib::BaseAccessor(aSession.printer())
83,fSession(aSession)
84,fManager(aManager)
85,fType("PhysicsTable")
86,fAIDA(aAIDA)
87//////////////////////////////////////////////////////////////////////////////
88//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
89{
90 addProperty("table",Lib::Property::STRING);
91 addProperty("process",Lib::Property::STRING);
92 addProperty("particle",Lib::Property::STRING);
93 addProperty("number",Lib::Property::INTEGER);
94 addProperty("name",Lib::Property::STRING);
95}
96//////////////////////////////////////////////////////////////////////////////
97G4Lab::PhysicsTableAccessor::~PhysicsTableAccessor(
98)
99//////////////////////////////////////////////////////////////////////////////
100//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
101{
102 for(unsigned int index=0;index<fTables.size();index++) delete fTables[index];
103}
104//////////////////////////////////////////////////////////////////////////////
105void* G4Lab::PhysicsTableAccessor::cast(
106 const std::string& aClass
107) const
108//////////////////////////////////////////////////////////////////////////////
109//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
110{
111 if_Lib_SCast(Slash::Data::IVisualizer)
112 else return Lib::BaseAccessor::cast(aClass);
113}
114//////////////////////////////////////////////////////////////////////////////
115std::string G4Lab::PhysicsTableAccessor::name(
116) const
117//////////////////////////////////////////////////////////////////////////////
118//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
119{
120 return fType;
121}
122namespace G4Lab {
123 class PhysicsTableIterator : public virtual Slash::Data::IIterator {
124 public: //Slash::Data::IIterator
125 virtual Slash::Data::IAccessor::Data object() {
126 if(fIterator==fVector.end()) return 0;
127 return *fIterator;
128 }
129 virtual void next() { ++fIterator; }
130 virtual void* tag() { return 0;}
131 public:
132 PhysicsTableIterator(std::vector<G4Lab::PhysicsTable*>& aVector)
133 :fVector(aVector) {
134 fIterator = fVector.begin();
135 }
136 virtual ~PhysicsTableIterator() {}
137 private:
138 std::vector<G4Lab::PhysicsTable*>& fVector;
139 std::vector<G4Lab::PhysicsTable*>::iterator fIterator;
140 };
141}
142//////////////////////////////////////////////////////////////////////////////
143Slash::Data::IIterator* G4Lab::PhysicsTableAccessor::iterator(
144)
145//////////////////////////////////////////////////////////////////////////////
146//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
147{
148 if(!fTables.size()) {
149 if(!buildProcessTables()) return 0;
150 }
151 return new PhysicsTableIterator(fTables);
152}
153//////////////////////////////////////////////////////////////////////////////
154Slash::Core::IValue* G4Lab::PhysicsTableAccessor::findValue(
155 Slash::Data::IAccessor::Data aData
156,const std::string& aName
157,void*
158)
159//////////////////////////////////////////////////////////////////////////////
160//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
161{
162 PhysicsTable* obj = (PhysicsTable*)aData;
163 if(aName=="table") {
164 return new Lib::Value(obj->fTableName);
165 } else if(aName=="process") {
166 return new Lib::Value(obj->fProcessName);
167 } else if(aName=="particle") {
168 return new Lib::Value(obj->fParticleName);
169 } else if(aName=="number") {
170 return new Lib::Value((int)obj->fTable->size());
171 } else if(aName=="name") {
172 std::string name =
173 obj->fTableName + "." + obj->fProcessName + "." + obj->fParticleName;
174 return new Lib::Value(name);
175 } else {
176 return new Lib::Value();
177 }
178}
179//////////////////////////////////////////////////////////////////////////////
180void G4Lab::PhysicsTableAccessor::beginVisualize(
181)
182//////////////////////////////////////////////////////////////////////////////
183//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
184{
185 if(!fAIDA) {
186 Lib::Out out(fSession.printer());
187 out << "G4Lab::PhysicsTableAccessor::beginVisualize :"
188 << " AIDA not found. Can't then find an AIDA plotter."
189 << Lib::endl;
190 }
191}
192//////////////////////////////////////////////////////////////////////////////
193void G4Lab::PhysicsTableAccessor::visualize(
194 Slash::Data::IAccessor::Data aData
195,void*
196)
197//////////////////////////////////////////////////////////////////////////////
198//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
199{
200 PhysicsTable* obj = (PhysicsTable*)aData;
201
202 std::string value;
203
204 bool superpose = true;
205 if(fSession.parameterValue("modeling.superpose",value))
206 if(!Lib::smanip::tobool(value,superpose)) superpose = true;
207
208 if(!fAIDA) return;
209
210 AIDA::ITreeFactory* tf = fAIDA->createTreeFactory();
211 if(!tf) return;
212 AIDA::ITree* memoryTree = tf->create();
213 delete tf;
214 if(!memoryTree) return;
215
216 AIDA::IFunctionFactory* functionFactory =
217 fAIDA->createFunctionFactory(*memoryTree);
218 if(!functionFactory) {
219 Lib::Out out(fSession.printer());
220 out << "Can't create an function factory." << Lib::endl;
221 return;
222 }
223
224 AIDA::IPlotterFactory* plotterFactory = fAIDA->createPlotterFactory();
225 if(!plotterFactory) {
226 Lib::Out out(fSession.printer());
227 out << "Can't create a plotter factory." << Lib::endl;
228 delete functionFactory;
229 return;
230 }
231
232 AIDA::IPlotter* plotter = plotterFactory->create();
233 if(!plotter) {
234 Lib::Out out(fSession.printer());
235 out << "Can't create a plotter factory." << Lib::endl;
236 delete functionFactory;
237 delete plotterFactory;
238 return;
239 }
240
241 int number = obj->fTable->entries();
242 if(!superpose) {
243 if(number>0) plotter->createRegions(1,number,0);
244 else plotter->createRegions(1,1,0);
245 }
246
247 for (int count=0;count<number;count++){
248 G4PhysicsVector* physicsVector = (*(obj->fTable))[count];
249 size_t n = physicsVector->GetVectorLength();
250 if(n<2) {
251 Lib::Out out(fSession.printer());
252 out << "PhysicsVector with " << int(n) << " entries. Can't plot."
253 << Lib::endl;
254 continue;
255 }
256
257 std::string name;
258 Lib::smanip::printf(name,1024,"%s.%s.%s.%d",
259 obj->fTableName.c_str(),
260 obj->fProcessName.c_str(),
261 obj->fParticleName.c_str(),
262 count);
263 std::vector<double> params;
264 params.push_back(1); // Dimension of the grid (here one).
265 params.push_back(n); // Number of entries for first grid axis.
266 size_t i;
267 for (i=0;i<n;i++) params.push_back(i); // Xs.
268 for (i=0;i<n;i++) params.push_back((*physicsVector)[i]); //Values.
269 AIDA::IFunction* function =
270 functionFactory->createFunctionByName(name,"Grid1D");
271 if(!function) {
272 Lib::Out out(fSession.printer());
273 out << "Can't create \"" << name << "\" function." << Lib::endl;
274 break;
275 }
276 function->setParameters(params);
277
278 AIDA::IPlotterRegion& region = plotter->currentRegion();
279
280 region.plot(*function);
281 region.setParameter("plotter.xAxisAutomated","FALSE");
282 std::string value1;
283 Lib::smanip::printf(value1,32,"%g",0.);
284 region.setParameter("plotter.xAxisMinimum",value1);
285 std::string value2;
286 Lib::smanip::printf(value2,32,"%g",(double)(n-1));
287 region.setParameter("plotter.xAxisMaximum",value2);
288 //std::string value3;
289 //Lib::smanip::printf(value3,32,"%d",4 * n);
290 //region.setParameter("plotter.xNumberOfPoints",value3);
291 /*
292 fSession.out().println("debug : min : %s, max : %s, steps : %s\n.",
293 value1.c_str(),value2.c_str(),value3.c_str());
294 */
295 if(!superpose) plotter->next();
296
297 }
298
299 delete plotter;
300 delete functionFactory;
301 delete plotterFactory;
302
303}
304//////////////////////////////////////////////////////////////////////////////
305void G4Lab::PhysicsTableAccessor::endVisualize(
306)
307//////////////////////////////////////////////////////////////////////////////
308//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
309{
310}
311//////////////////////////////////////////////////////////////////////////////
312bool G4Lab::PhysicsTableAccessor::buildProcessTables(
313)
314//////////////////////////////////////////////////////////////////////////////
315//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
316{
317 if(fTables.size()) return true; //done.
318
319 if(!fManager.isRunning()) {
320 Lib::Out out(fSession.printer());
321 out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
322 << " It is needed to have started a run to get physics tables."
323 << Lib::endl;
324 return false;
325 }
326
327 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
328 if(!particleTable) {
329 Lib::Out out(fSession.printer());
330 out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
331 << " Can't get particle table."
332 << Lib::endl;
333 return false;
334 }
335
336 G4ParticleTable::G4PTblDicIterator* particleTableIterator =
337 particleTable->GetIterator();
338 if(!particleTableIterator) {
339 Lib::Out out(fSession.printer());
340 out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
341 << " Can't get particle table iterator."
342 << Lib::endl;
343 return false;
344 }
345
346 particleTableIterator->reset();
347 while( (*particleTableIterator)() ){
348 G4ParticleDefinition* particle = particleTableIterator->value();
349 if(!particle) {
350 Lib::Out out(fSession.printer());
351 out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
352 << " Can't get particle table particle."
353 << Lib::endl;
354 return false;
355 }
356 G4String particleName = particle->GetParticleName();
357 //Lib::Out out(fSession.printer());
358 //out << "debug : Particle : " << particleName << Lib::endl;
359
360 // Processes for this particle :
361 G4ProcessManager* processManager = particle->GetProcessManager();
362 if(!processManager) {
363 Lib::Out out(fSession.printer());
364 out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
365 << " Can't get process manager of particle "
366 << particleName << "."
367 << Lib::endl;
368 return false;
369 }
370
371 G4ProcessVector* processList = processManager->GetProcessList();
372 if(!processList) {
373 Lib::Out out(fSession.printer());
374 out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
375 << " Can't get process list for process manager of particle "
376 << particleName << "."
377 << Lib::endl;
378 return false;
379 }
380
381 int number = processList->entries();
382 for (int index=0;index<number;index++){
383 G4VProcess* process = (*processList)(index);
384 std::vector<std::string> names;
385 std::vector<G4PhysicsTable*> tables;
386 if(!findProcessTables(particle,process,names,tables)) {
387 return false;
388 }
389
390 for(unsigned int i=0;i<tables.size();i++) {
391 std::vector<std::string> words;
392 Lib::smanip::words(names[i],".",words);
393 if(words.size()!=3) {
394 Lib::Out out(fSession.printer());
395 out << "G4Lab::PhysicsTableAccessor::buildProcessTables :"
396 << " Three words separted with a dot expected in "
397 << names[i] << " for particle " << particleName << "."
398 << Lib::endl;
399 return false;
400 }
401 fTables.push_back(new PhysicsTable(words[0],
402 words[1],
403 words[2],
404 tables[i],
405 process,
406 particle));
407 }
408 }
409 }
410
411 return true;
412}
413//////////////////////////////////////////////////////////////////////////////
414bool G4Lab::PhysicsTableAccessor::findProcessTables(
415 G4ParticleDefinition* aParticle
416,G4VProcess* aProcess
417,std::vector<std::string>& aTables
418,std::vector<G4PhysicsTable*>& aPhysicsTables
419)
420//////////////////////////////////////////////////////////////////////////////
421//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
422{
423 aTables.clear();
424 aPhysicsTables.clear();
425
426 std::string tmpName;
427 if(!Lib::fmanip::tmpname(tmpName,".","G4Lab_tmp","")) {
428 Lib::Out out(fSession.printer());
429 out << "Can't build a temporary directory name." << Lib::endl;
430 return false;
431 }
432
433 if(!Lib::dirmanip::create(tmpName)) {
434 Lib::Out out(fSession.printer());
435 out << "Can't create \"" << tmpName << "\" directory." << Lib::endl;
436 return false;
437 }
438
439#ifdef WIN32 //binary mode out on Windows.
440 G4bool ascii = true;
441#else
442 G4bool ascii = false;
443 //G4bool ascii = true;
444#endif
445
446 if(!aProcess->StorePhysicsTable(aParticle,tmpName,ascii)) {
447 Lib::Out out(fSession.printer());
448 out << "Can't store physics table for particle "
449 << aParticle->GetParticleName()
450 << " and process " << aProcess->GetProcessName() << "."
451 << Lib::endl;
452 //return false;
453 } else {
454
455 // Analyse the content of the directory :
456 std::vector<std::string> files;
457 if(!Lib::dirmanip::entries(tmpName,files,false)) {
458 Lib::Out out(fSession.printer());
459 out << "Can't get files in \"" << tmpName << "\" directory."
460 << Lib::endl;
461 } else {
462 // Analyse files :
463 Lib::smanip::remove(files,".");
464 Lib::smanip::remove(files,"..");
465/*
466 if(!files.size()) {
467 Lib::Out out(fSession.printer());
468 out << "No files in \"" << tmpName << "\" directory for particle "
469 << aParticle->GetParticleName()
470 << " and process " << aProcess->GetProcessName() << "."
471 << Lib::endl;
472 //return false;
473 }
474*/
475 for(unsigned int index=0;index<files.size();index++) {
476 std::vector<std::string> words;
477 Lib::smanip::words(files[index],".",words);
478 //Lib::Out out(fSession.printer());
479 //out << "debug : Analyse \"" << files[index] << "\" ." << Lib::endl;
480 if(words.size()!=4) {
481 Lib::Out out(fSession.printer());
482 out << "Bad syntax in file name " << Lib::sout(files[index])
483 << " for particle " << aParticle->GetParticleName()
484 << " and process " << aProcess->GetProcessName() << "."
485 << Lib::endl;
486 //aTables.clear();
487 //aPhysicsTables.clear();
488 //return false;
489 continue;
490 }
491 std::string table = words[0];
492 std::string process = words[1];
493 std::string particle = words[2];
494 //Lib::smanip::replace(aTables[i],".dat","");
495 std::string filename =
496 aProcess->GetPhysicsTableFileName(aParticle,
497 tmpName,
498 table,
499 ascii);
500 G4PhysicsTable* physicsTable =
501 new G4PhysicsTable(G4Element::GetNumberOfElements());
502 G4bool status = physicsTable->RetrievePhysicsTable(filename,ascii);
503 if(!status){
504 Lib::Out out(fSession.printer());
505 out << "Unable to read \"" << filename << "\" file." << Lib::endl;
506 delete physicsTable;
507 physicsTable = 0;
508 }
509 if(physicsTable) {
510 aTables.push_back(table+"."+process+"."+particle);
511 aPhysicsTables.push_back(physicsTable);
512 }
513 }
514 }
515
516 }
517
518 // Remove the temporary directory :
519 if(!Lib::dirmanip::remove(tmpName)) {
520 Lib::Out out(fSession.printer());
521 out << "Can't remove \"" << tmpName << "\" directory." << Lib::endl;
522 aTables.clear();
523 aPhysicsTables.clear();
524 return false;
525 }
526
527 return true;
528}
Note: See TracBrowser for help on using the repository browser.