Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/visualization/gMocren/src/G4GMocrenFileSceneHandler.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 //
 29 // Created:  Mar. 31, 2009  Akinori Kimura  
 30 //           Sep. 22, 2009  Akinori Kimura : modify and fix code to support 
 31 //                                          PrimitiveScorers and clean up
 32 //
 33 // GMocrenFile scene handler
 34 
 35 
 36 //----- header files
 37 #include <fstream>
 38 #include <cstdlib>
 39 #include <cstring>
 40 #include <sstream>
 41 #include <sstream>
 42 #include <iomanip>
 43 
 44 #include "globals.hh"
 45 #include "G4VisManager.hh"
 46 
 47 #include "G4GMocrenFile.hh"
 48 #include "G4GMocrenFileSceneHandler.hh"
 49 #include "G4GMocrenFileViewer.hh"
 50 #include "G4Point3D.hh"
 51 #include "G4VisAttributes.hh"
 52 #include "G4Scene.hh"
 53 #include "G4Transform3D.hh"
 54 #include "G4Polyhedron.hh"
 55 #include "G4Box.hh"
 56 #include "G4Cons.hh"
 57 #include "G4Polyline.hh"
 58 #include "G4Trd.hh"
 59 #include "G4Tubs.hh"
 60 #include "G4Trap.hh"
 61 #include "G4Torus.hh"
 62 #include "G4Sphere.hh"
 63 #include "G4Para.hh"
 64 #include "G4Text.hh"
 65 #include "G4Circle.hh"
 66 #include "G4Square.hh"
 67 #include "G4VPhysicalVolume.hh"
 68 #include "G4PhysicalVolumeModel.hh"
 69 #include "G4LogicalVolume.hh"
 70 #include "G4Material.hh"
 71 
 72 #include "G4VPVParameterisation.hh"
 73 #include "G4VVolumeMaterialScanner.hh"
 74 #include "G4VisTrajContext.hh"
 75 #include "G4TrajectoriesModel.hh"
 76 #include "G4VTrajectoryModel.hh"
 77 #include "G4TrajectoryDrawByCharge.hh"
 78 #include "G4HitsModel.hh"
 79 #include "G4GMocrenMessenger.hh"
 80 //#include "G4PSHitsModel.hh"
 81 #include "G4GMocrenIO.hh"
 82 #include "G4VNestedParameterisation.hh"
 83 #include "G4GMocrenTouchable.hh"
 84 #include "G4GMocrenFileCTtoDensityMap.hh"
 85 #include "G4PhantomParameterisation.hh"
 86 
 87 #include "G4ScoringManager.hh"
 88 #include "G4ScoringBox.hh"
 89 
 90 #include "G4PhysicalConstants.hh"
 91 #include "G4SystemOfUnits.hh"
 92 
 93 #include <utility>
 94 
 95 //----- constants
 96 const char  GDD_FILE_HEADER      [] = "g4_";
 97 const char  DEFAULT_GDD_FILE_NAME[] = "g4_00.gdd";
 98 
 99 const G4int FR_MAX_FILE_NUM = 100 ;
100 const G4int MAX_NUM_TRAJECTORIES = 100000;
101 
102 //-- for a debugging
103 const G4bool GFDEBUG = false;
104 const G4bool GFDEBUG_TRK = false;//true;
105 const G4bool GFDEBUG_HIT = false;//true;
106 const G4bool GFDEBUG_DIGI = false;//true;
107 const G4int GFDEBUG_DET = 0; // 0: false 
108 
109 //////////////////////
110 // static variables //
111 //////////////////////
112 
113 //----- static variables
114 G4int G4GMocrenFileSceneHandler::kSceneIdCount = 0; 
115 
116 ///////////////////////////
117 // Driver-dependent part //
118 ///////////////////////////
119 
120 
121 //----- G4GMocrenFileSceneHandler, constructor
122 G4GMocrenFileSceneHandler::G4GMocrenFileSceneHandler(G4GMocrenFile& system,
123                  G4GMocrenMessenger & messenger,
124                  const G4String& name)
125   : G4VSceneHandler(system, kSceneIdCount++, name),
126     kSystem(system),
127     kMessenger(messenger),
128     kgMocrenIO(new G4GMocrenIO()),
129     kbSetModalityVoxelSize(false),
130     kbModelingTrajectory(false),
131 //    kGddDest(0),
132     kFlagInModeling(false),
133     kFlagSaving_g4_gdd(false),
134     kFlagParameterization(0),
135     kFlagProcessedInteractiveScorer(false) {
136 
137   // g4.gdd filename and its directory
138   if(std::getenv("G4GMocrenFile_DEST_DIR") == NULL) {
139     kGddDestDir[0] = '\0';
140     //std::strcpy(kGddDestDir , "");                    // output dir
141     //std::strcpy(kGddFileName, DEFAULT_GDD_FILE_NAME); // filename
142     std::strncpy(kGddFileName, DEFAULT_GDD_FILE_NAME,
143      std::strlen(DEFAULT_GDD_FILE_NAME)+1);   // filename
144   } else {
145     const char * env = std::getenv("G4GMocrenFile_DEST_DIR");
146     G4int len = (G4int)std::strlen(env);
147     if(len > 256) {
148       G4Exception("G4GMocrenFileSceneHandler::G4GMocrenFileSceneHandler(*)",
149                   "gMocren1000", FatalException,
150                   "Invalid length of string set in G4GMocrenFile_DEST_DIR");
151     }
152     std::strncpy(kGddDestDir, env, len+1);  // output dir
153     std::strncpy(kGddFileName, DEFAULT_GDD_FILE_NAME,
154                  std::strlen(DEFAULT_GDD_FILE_NAME)+1);  // filename
155   }
156     
157   // maximum number of g4.gdd files in the dest directory
158   kMaxFileNum = FR_MAX_FILE_NUM ; // initialization
159   if ( std::getenv( "G4GMocrenFile_MAX_FILE_NUM" ) != NULL ) {  
160     char * pcFileNum = std::getenv("G4GMocrenFile_MAX_FILE_NUM");
161     char c10FileNum[10];
162     std::strncpy(c10FileNum, pcFileNum, 9);
163     c10FileNum[9] = '\0';
164     kMaxFileNum = std::atoi(c10FileNum);
165 
166   } else {
167     kMaxFileNum = FR_MAX_FILE_NUM ;
168   }
169   if( kMaxFileNum < 1 ) { kMaxFileNum = 1 ; }
170 
171   InitializeParameters();
172 
173 } 
174 
175 
176 //----- G4GMocrenFileSceneHandler, destructor
177 G4GMocrenFileSceneHandler::~G4GMocrenFileSceneHandler () 
178 {
179   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
180       G4cout << "***** ~G4GMocrenFileSceneHandler" << G4endl;
181 
182   if(kGddDest) {
183     //----- End of modeling
184     // close g4.gdd
185     GFEndModeling();
186   }
187   if(kgMocrenIO != NULL) delete kgMocrenIO;
188 
189 }
190 
191 //----- initialize all parameters
192 void G4GMocrenFileSceneHandler::InitializeParameters() {
193 
194   kbSetModalityVoxelSize = false;
195 
196   for(G4int i = 0; i < 3; i++) {
197     kModalitySize[i] = 0;
198     kNestedVolumeDimension[i] = 0;
199     kNestedVolumeDirAxis[i] = -1;
200   }
201 
202   // delete kgMocrenIO;
203 
204 }
205 
206 //-----
207 void G4GMocrenFileSceneHandler::SetGddFileName() 
208 {
209   // g4_00.gdd, g4_01.gdd, ..., g4_MAX_FILE_INDEX.gdd
210   const G4int MAX_FILE_INDEX = kMaxFileNum - 1 ;
211 
212   // dest directory (null if no environmental variables is set)
213   std::strncpy(kGddFileName, kGddDestDir, sizeof(kGddFileName)-1);
214   kGddFileName[sizeof(kGddFileName)-1] = '\0';
215 
216   // create full path name (default)
217   std::strncat ( kGddFileName, DEFAULT_GDD_FILE_NAME,
218                 sizeof(kGddFileName) - std::strlen(kGddFileName) - 1);
219 
220   // Automatic updation of file names
221   static G4int currentNumber = 0;
222   for( G4int i = currentNumber ; i < kMaxFileNum ; i++) { 
223 
224     // Message in the final execution
225     if( i == MAX_FILE_INDEX ) 
226       {
227   if (G4VisManager::GetVerbosity() >= G4VisManager::warnings) {
228     G4cout << "==========================================="   << G4endl; 
229     G4cout << "WARNING MESSAGE from GMocrenFile driver:   "   << G4endl;
230     G4cout << "  This file name is the final one in the   "   << G4endl;
231     G4cout << "  automatic updation of the output file name." << G4endl; 
232     G4cout << "  You may overwrite existing files, i.e.   "   << G4endl; 
233     G4cout << "  g4_XX.gdd."   << G4endl;
234     G4cout << "==========================================="   << G4endl; 
235   }
236       }
237 
238     // re-determine file name as G4GMocrenFile_DEST_DIR/g4_XX.gdd 
239     std::ostringstream filename;
240     filename
241     << kGddDestDir << GDD_FILE_HEADER
242     << std::setw(2) << std::setfill('0') << i << ".wrl";
243     strncpy(kGddFileName,filename.str().c_str(),sizeof(kGddFileName)-1);
244     kGddFileName[sizeof(kGddFileName)-1] = '\0';
245 
246     // check validity of the file name
247     std::ifstream fin(kGddFileName); 
248     if(GFDEBUG)
249       G4cout << "FILEOPEN: " << i << " : " << kGddFileName << fin.fail()
250        << G4endl;
251     if(!fin) { 
252       // new file 
253       fin.close();
254       currentNumber = i+1;
255       break; 
256     } else { 
257       // already exists (try next) 
258       fin.close(); 
259     } 
260 
261   } // for 
262 
263   G4cout << "======================================================================" << G4endl; 
264   G4cout << "Output file: " << kGddFileName                          << G4endl; 
265   G4cout << "Destination directory (current dir if NULL): " << kGddDestDir << G4endl; 
266   G4cout << "Maximum number of files in the destination directory: " << kMaxFileNum << G4endl; 
267   G4cout << "Note:" << G4endl; 
268   G4cout << "  * The maximum number is customizable as:           " << G4endl;
269   G4cout << "      % setenv  G4GMocrenFile_MAX_FILE_NUM  number " << G4endl;        
270   G4cout << "  * The destination directory is customizable as:" << G4endl;
271   G4cout << "      % setenv  G4GMocrenFile_DEST_DIR  dir_name/  " << G4endl;        
272   G4cout << "     ** Do not forget \"/\" at the end of the dir_name, e.g. \"./tmp/\"." << G4endl;              
273   //G4cout << "        dir_name, e.g. \"./tmp/\"."                 << G4endl;              
274   G4cout << G4endl;
275   G4cout << "Maximum number of trajectories is set to " << MAX_NUM_TRAJECTORIES << "."<< G4endl;
276   G4cout << "======================================================================" << G4endl; 
277 
278 } // G4GMocrenFileSceneHandler::SetGddFileName()
279 
280 
281 //-----
282 void  G4GMocrenFileSceneHandler::BeginSavingGdd( void )
283 {
284   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
285       G4cout << "***** BeginSavingGdd (called)" << G4endl;
286 
287   if( !IsSavingGdd() ) {
288 
289     if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations) {
290       G4cout << "*****                   (started) " ;
291       G4cout << "(open g4.gdd, ##)"  << G4endl;
292     }
293 
294     SetGddFileName() ; // result set to kGddFileName
295     kFlagSaving_g4_gdd = true; 
296 
297 
298     G4GMocrenFileCTtoDensityMap ctdens;
299     short minmax[2];
300     minmax[0] = ctdens.GetMinCT();
301     minmax[1] = ctdens.GetMaxCT();
302     kgMocrenIO->setModalityImageMinMax(minmax);
303     std::vector<G4float> map;
304     G4float dens;
305     for(G4int i = minmax[0]; i <= minmax[1]; i++) {
306       dens = ctdens.GetDensity(i);
307       map.push_back(dens);
308     }
309     kgMocrenIO->setModalityImageDensityMap(map);
310 
311     /*
312     G4String fname = "modality-map.dat";
313     std::ifstream ifile(fname);
314     if(ifile) {
315       short minmax[2];
316       ifile >> minmax[0] >> minmax[1];
317       kgMocrenIO->setModalityImageMinMax(minmax);
318       std::vector<G4float> map;
319       G4float dens;
320       for(G4int i = minmax[0]; i <= minmax[1]; i++) {
321   ifile >> dens;
322   map.push_back(dens);
323       }
324       kgMocrenIO->setModalityImageDensityMap(map);
325       
326     } else {
327       G4cout << "cann't open the file : " << fname << G4endl;
328     }
329     */
330 
331     // mesh size
332     //kMessenger.getNoVoxels(kModalitySize[0], kModalitySize[1], kModalitySize[2]);
333     //kgMocrenIO->setModalityImageSize(kModalitySize);
334     
335     // initializations
336     //kgMocrenIO->clearModalityImage();
337     kgMocrenIO->clearDoseDistAll();
338     kgMocrenIO->clearROIAll();
339     kgMocrenIO->clearTracks();
340     kgMocrenIO->clearDetector();
341     std::vector<Detector>::iterator itr = kDetectors.begin();
342     for(; itr != kDetectors.end(); itr++) {
343       itr->clear();
344     }
345     kDetectors.clear();
346     
347     kNestedHitsList.clear();
348     kNestedVolumeNames.clear();
349       
350   }
351 }
352 
353 void  G4GMocrenFileSceneHandler::EndSavingGdd  ( void ) 
354 {
355   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
356     G4cout << "***** EndSavingGdd (called)" << G4endl;
357 
358   if(IsSavingGdd()) {
359     if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
360       G4cout << "*****                 (started) (close "
361        << kGddFileName << ")" << G4endl;
362 
363     if(kGddDest) kGddDest.close();
364     kFlagSaving_g4_gdd = false; 
365 
366     std::map<Index3D, G4float>::iterator itr = kNestedModality.begin();
367     G4int xmax=0, ymax=0, zmax=0;
368     for(; itr != kNestedModality.end(); itr++) {
369       if(itr->first.x > xmax) xmax = itr->first.x;
370       if(itr->first.y > ymax) ymax = itr->first.y;
371       if(itr->first.z > zmax) zmax = itr->first.z;
372     }
373     // mesh size
374     kModalitySize[0] = xmax+1;
375     kModalitySize[1] = ymax+1;
376     kModalitySize[2] = zmax+1;
377     kgMocrenIO->setModalityImageSize(kModalitySize);
378     if(GFDEBUG) G4cout << "gMocren-file driver : modality size : "
379            << kModalitySize[0] << " x "
380            << kModalitySize[1] << " x "
381            << kModalitySize[2] << G4endl;
382 
383     G4int nxy = kModalitySize[0]*kModalitySize[1];
384     //std::map<G4int, G4float>::iterator itr;
385     for(G4int z = 0; z < kModalitySize[2]; z++) {
386       short * modality = new short[nxy];
387       for(G4int y = 0; y < kModalitySize[1]; y++) {
388   for(G4int x = 0; x < kModalitySize[0]; x++) {
389     //for(G4int x = kModalitySize[0]-1; x >= 0 ; x--) {
390     //G4int ixy = x + (kModalitySize[1]-y-1)*kModalitySize[0];
391 
392     G4int ixy = x + y*kModalitySize[0];
393     Index3D idx(x,y,z);
394     itr = kNestedModality.find(idx);
395     if(itr != kNestedModality.end()) {
396 
397       modality[ixy] = kgMocrenIO->convertDensityToHU(itr->second);
398     } else {
399       modality[ixy] = -1024;
400     }
401 
402   }
403       }
404       kgMocrenIO->setModalityImage(modality);
405     }
406 
407     //-- dose
408     size_t nhits = kNestedHitsList.size();
409     if(GFDEBUG) G4cout << "gMocren-file driver : # hits = " << nhits << G4endl;
410 
411     std::map<Index3D, G4double>::iterator hitsItr;
412     std::map<G4String, std::map<Index3D, G4double> >::iterator hitsListItr = kNestedHitsList.begin();
413 
414     for(G4int n = 0; hitsListItr != kNestedHitsList.end(); hitsListItr++, n++) {
415 
416       kgMocrenIO->newDoseDist();
417       kgMocrenIO->setDoseDistName(hitsListItr->first, n);
418       kgMocrenIO->setDoseDistSize(kModalitySize, n);
419 
420       G4double minmax[2] = {DBL_MAX, -DBL_MAX};
421       for(G4int z = 0 ; z < kModalitySize[2]; z++) {
422   G4double * values = new G4double[nxy];
423   for(G4int y = 0; y < kModalitySize[1]; y++) {
424     for(G4int x = 0; x < kModalitySize[0]; x++) {
425 
426       G4int ixy = x + y*kModalitySize[0];
427       Index3D idx(x,y,z);
428       hitsItr = hitsListItr->second.find(idx);
429       if(hitsItr != hitsListItr->second.end()) {
430 
431         values[ixy] = hitsItr->second;
432       } else {
433         values[ixy] = 0.;
434       }
435       if(values[ixy] < minmax[0]) minmax[0] = values[ixy];
436       if(values[ixy] > minmax[1]) minmax[1] = values[ixy];
437     }
438   }
439   kgMocrenIO->setDoseDist(values, n);
440       }
441       kgMocrenIO->setDoseDistMinMax(minmax, n);
442       G4double lower = 0.;
443       if(minmax[0] < 0)  lower = minmax[0];
444       G4double scale = (minmax[1]-lower)/25000.;
445       kgMocrenIO->setDoseDistScale(scale, n);
446       G4String sunit("unit?"); //temporarily
447       kgMocrenIO->setDoseDistUnit(sunit, n);
448     }
449       
450 
451     //-- draw axes
452     if(false) {//true,false
453       G4ThreeVector trans;
454       G4RotationMatrix rot;
455       trans = kVolumeTrans3D.getTranslation();
456       rot = kVolumeTrans3D.getRotation().inverse();
457       // x
458       std::vector<G4float *> tracks;
459       unsigned char colors[3];
460       G4float * trk = new G4float[6];
461       tracks.push_back(trk);
462 
463       G4ThreeVector orig(0.,0.,0), xa(2000.,0.,0.), ya(0.,2000.,0.), za(0.,0.,2000.);
464       orig -= trans;
465       orig.transform(rot);
466       xa -= trans;
467       xa.transform(rot);
468       ya -= trans;
469       ya.transform(rot);
470       za -= trans;
471       za.transform(rot);
472       for(G4int i = 0; i < 3; i++) trk[i] = orig[i];
473       for(G4int i = 0; i < 3; i++) trk[i+3] = xa[i];
474       colors[0] = 255; colors[1] = 0; colors[2] = 0;
475       kgMocrenIO->addTrack(tracks, colors);
476       // y
477       for(G4int i = 0; i < 3; i++) trk[i+3] = ya[i];
478       colors[0] = 0; colors[1] = 255; colors[2] = 0;
479       kgMocrenIO->addTrack(tracks, colors);
480       // z
481       for(G4int i = 0; i < 3; i++) trk[i+3] = za[i];
482       colors[0] = 0; colors[1] = 0; colors[2] = 255;
483       kgMocrenIO->addTrack(tracks, colors);
484     }
485 
486     //-- detector
487     ExtractDetector();
488 
489 
490     if(GFDEBUG_DET) G4cout << ">>>>>>>>>>>>>>>>>>>>>>   (";
491     std::vector<G4float> transformObjects;
492     for(G4int i = 0; i < 3; i++) {
493       // need to check!!
494       transformObjects.push_back((kVolumeSize[i]/2. - kVoxelDimension[i]/2.));
495       if(GFDEBUG_DET) G4cout << transformObjects[i] << ", ";
496     }
497     if(GFDEBUG_DET) G4cout << ")" << G4endl;
498 
499 
500     kgMocrenIO->translateTracks(transformObjects);
501     kgMocrenIO->translateDetector(transformObjects);
502 
503     // store
504     kgMocrenIO->storeData(kGddFileName);
505   } 
506 
507 }
508 
509 
510 //----- 
511 void G4GMocrenFileSceneHandler::GFBeginModeling( void )
512 {
513   G4VSceneHandler::BeginModeling();
514 
515   if( !GFIsInModeling() ) {
516 
517     if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
518       G4cout << "***** G4GMocrenFileSceneHandler::GFBeginModeling (called & started)" << G4endl;
519 
520     //----- Send saving command and heading comment
521     BeginSavingGdd();
522 
523     kFlagInModeling = true ;
524 
525     // These models are entrusted to user commands /vis/scene/add/psHits or hits
526     //GetScene()->AddEndOfEventModel(new G4PSHitsModel());
527     //GetScene()->AddEndOfRunModel(new G4PSHitsModel());
528     //scene->AddEndOfEventModel(new G4HitsModel());
529     if(GFDEBUG_HIT) {
530       G4Scene * scene = GetScene();
531       std::vector<G4Scene::Model> vmodel = scene->GetEndOfEventModelList();
532       std::vector<G4Scene::Model>::iterator itr = vmodel.begin();
533       for(; itr != vmodel.end(); itr++) {
534         G4cout << " IIIIII model name: " << itr->fpModel->GetGlobalTag() << G4endl;
535       }
536     }
537   }
538 }
539 
540 
541 //========== AddPrimitive() functions ==========//
542 
543 //----- Add polyline 
544 void G4GMocrenFileSceneHandler::AddPrimitive (const G4Polyline& polyline) 
545 {
546   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
547     G4cout << "***** AddPrimitive" << G4endl;
548 
549   if (fProcessing2D) {
550     static G4bool warned = false;
551     if (!warned) {
552       warned = true;
553       G4Exception
554   ("G4GMocrenFileSceneHandler::AddPrimitive (const G4Polyline&)",
555    "gMocren1001", JustWarning,
556    "2D polylines not implemented.  Ignored.");
557     }
558     return;
559   }
560 
561   //----- Initialize if necessary
562   GFBeginModeling();
563 
564   static G4int numTrajectories = 0;
565   if(numTrajectories >= MAX_NUM_TRAJECTORIES) return;
566 
567   // draw trajectories
568   if(kbModelingTrajectory) {
569     
570     G4TrajectoriesModel * pTrModel = dynamic_cast<G4TrajectoriesModel*>(fpModel);
571     if (!pTrModel) { 
572       G4Exception 
573   ("G4VSceneHandler::AddCompound(const G4Polyline&)",
574    "gMocren0002", FatalException, "Not a G4TrajectoriesModel.");
575     }
576 
577     G4ThreeVector trans;
578     G4RotationMatrix rot;
579     trans = kVolumeTrans3D.getTranslation();
580     rot = kVolumeTrans3D.getRotation().inverse();
581 
582     if(GFDEBUG_TRK) G4cout << "   trajectory points : " << G4endl;
583     std::vector<G4float *> trajectory;
584     if(polyline.size() < 2) return;
585     G4Polyline::const_iterator preitr = polyline.begin();
586     G4Polyline::const_iterator postitr = preitr; postitr++;
587     for(; postitr != polyline.end(); preitr++, postitr++) {
588       G4ThreeVector prePts(preitr->x(), preitr->y(), preitr->z());
589       prePts -= trans;
590       prePts.transform(rot);
591       G4ThreeVector postPts(postitr->x(), postitr->y(), postitr->z());
592       postPts -= trans;
593       postPts.transform(rot);
594       G4float * stepPts = new G4float[6];
595       stepPts[0] = prePts.x();
596       stepPts[1] = prePts.y();
597       stepPts[2] = prePts.z();
598       stepPts[3] = postPts.x();
599       stepPts[4] = postPts.y();
600       stepPts[5] = postPts.z();
601       trajectory.push_back(stepPts);
602 
603       if(GFDEBUG_TRK) {
604   G4cout << "                       ("
605          << stepPts[0] << ", "
606          << stepPts[1] << ", "
607          << stepPts[2] << ") - (" 
608          << stepPts[3] << ", "
609          << stepPts[4] << ", "
610          << stepPts[5] << ")" << G4endl;
611       }
612     }
613 
614     const G4VisAttributes * att = polyline.GetVisAttributes();
615     G4Color color = att->GetColor();
616     unsigned char trkcolor[3];
617     trkcolor[0] = (unsigned char)(color.GetRed()*255);
618     trkcolor[1] = (unsigned char)(color.GetGreen()*255);
619     trkcolor[2] = (unsigned char)(color.GetBlue()*255);
620     if(GFDEBUG_TRK) {
621       G4cout << "              color  : ["
622        << color.GetRed() << ", "
623        << color.GetGreen() << ", "
624        << color.GetBlue() << "]" << G4endl;
625     }
626 
627     kgMocrenIO->addTrack(trajectory, trkcolor);
628 
629     numTrajectories++;
630   }
631 
632 } // G4GMocrenFileSceneHandler::AddPrimitive (polyline)
633 
634 
635 //----- Add text
636 void G4GMocrenFileSceneHandler::AddPrimitive ( const G4Text& text )
637 {
638   if (fProcessing2D) {
639     static G4bool warned = false;
640     if (!warned) {
641       warned = true;
642       G4Exception
643   ("G4GMocrenFileSceneHandler::AddPrimitive (const G4Text&)",
644    "gMocren1002", JustWarning,
645    "2D text not implemented.  Ignored.");
646     }
647     return;
648   }
649 
650   // to avoid a warning in the compile process
651   G4Text dummytext = text;
652 
653   //----- 
654   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
655     G4cout << "***** AddPrimitive( G4Text )" << G4endl;
656 
657   //----- Initialize IF NECESSARY
658   GFBeginModeling();
659 
660 } // G4GMocrenFileSceneHandler::AddPrimitive ( text )
661 
662 
663 //----- Add circle
664 void G4GMocrenFileSceneHandler::AddPrimitive ( const G4Circle& mark_circle )
665 {
666   // to avoid a warning in the compile process
667   G4Circle dummycircle = mark_circle;
668 
669   if (fProcessing2D) {
670     static G4bool warned = false;
671     if (!warned) {
672       warned = true;
673       G4Exception
674   ("G4GMocrenFileSceneHandler::AddPrimitive (const G4Circle&)",
675    "gMocren1003", JustWarning,
676    "2D circles not implemented.  Ignored.");
677     }
678     return;
679   }
680 
681   //----- 
682   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
683     G4cout << "***** AddPrimitive( G4Circle )" << G4endl;
684 
685   //----- Initialize IF NECESSARY
686   GFBeginModeling();
687 
688 
689 } // G4GMocrenFileSceneHandler::AddPrimitive ( mark_circle )
690 
691 
692 //----- Add square
693 void G4GMocrenFileSceneHandler::AddPrimitive (const G4Square& mark_square )
694 {
695   // to avoid a warning in the compile process
696   G4Square dummysquare = mark_square;
697 
698   if (fProcessing2D) {
699     static G4bool warned = false;
700     if (!warned) {
701       warned = true;
702       G4Exception
703   ("G4GMocrenFileSceneHandler::AddPrimitive (const G4Square&)",
704    "gMocren1004", JustWarning,
705    "2D squares not implemented.  Ignored.");
706     }
707     return;
708   }
709 
710   //----- 
711   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
712     G4cout << "***** AddPrimitive( G4Square )" << G4endl;
713 
714   //----- Initialize if necessary
715   GFBeginModeling();
716 
717 } // G4GMocrenFileSceneHandler::AddPrimitive ( mark_square )
718 
719 
720 //----- Add polyhedron
721 void G4GMocrenFileSceneHandler::AddPrimitive ( const G4Polyhedron& polyhedron ) 
722 {
723   //----- 
724   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
725     G4cout << "***** AddPrimitive( G4Polyhedron )" << G4endl;
726 
727 
728   if (polyhedron.GetNoFacets() == 0) return;
729 
730   if (fProcessing2D) {
731     static G4bool warned = false;
732     if (!warned) {
733       warned = true;
734       G4Exception
735   ("G4GMocrenFileSceneHandler::AddPrimitive (const G4Polyhedron&)",
736    "gMocren1005", JustWarning,
737    "2D polyhedra not implemented.  Ignored.");
738     }
739     return;
740   }
741 
742   //----- Initialize if necessary
743   GFBeginModeling();
744 
745   //---------- (3) Facet block
746   for (G4int f = polyhedron.GetNoFacets(); f; f--){
747     G4bool notLastEdge = true;
748     G4int index = -1; // initialization
749     G4int edgeFlag = 1;
750     //G4int preedgeFlag = 1;
751     //G4int work[4], i = 0;
752     G4int i = 0;
753     do {
754       //preedgeFlag = edgeFlag;
755       notLastEdge = polyhedron.GetNextVertexIndex(index, edgeFlag);
756       //work[i++] = index;
757       i++;
758     }while (notLastEdge);
759     switch (i){
760     case 3:
761       //SendStrInt3(FR_FACET, work[0], work[1], work[2] );
762       break;
763     case 4:
764       //SendStrInt4(FR_FACET, work[0], work[1], work[2], work[3] );
765       break;
766     default:
767       if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
768   G4cout <<
769     "ERROR G4GMocrenFileSceneHandler::AddPrimitive(G4Polyhedron)" << G4endl;
770       G4PhysicalVolumeModel* pPVModel =
771         dynamic_cast<G4PhysicalVolumeModel*>(fpModel);
772       if (pPVModel)   
773   if(G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
774     G4cout << "Volume " << pPVModel->GetCurrentPV()->GetName() <<
775       ", Solid " << pPVModel->GetCurrentLV()->GetSolid()->GetName() <<
776       " (" << pPVModel->GetCurrentLV()->GetSolid()->GetEntityType();
777 
778       if(G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
779   G4cout <<
780     "\nG4Polyhedron facet with " << i << " edges" << G4endl;  
781     }
782   }
783 
784 } // G4GMocrenFileSceneHandler::AddPrimitive (polyhedron) 
785 
786 
787 //----- 
788 void G4GMocrenFileSceneHandler::GFEndModeling ()
789 {
790   G4VSceneHandler::EndModeling();
791 
792   //-----     
793   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
794     G4cout << "***** GFEndModeling (called)" << G4endl;
795 
796   if( GFIsInModeling() ) {
797 
798     if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations) {
799       G4cout << "***** GFEndModeling (started) " ; 
800       G4cout << "(/EndModeling, /DrawAll, /CloseDevice)" << G4endl;
801     }
802 
803     //----- End saving data to g4.gdd
804     EndSavingGdd() ;
805 
806     //------ Reset flag 
807     kFlagInModeling = false ;
808 
809   }
810 
811 }
812 
813 
814 //----- 
815 void G4GMocrenFileSceneHandler::BeginPrimitives (const G4Transform3D& objectTransformation)
816 {
817   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
818     G4cout << "***** BeginPrimitives " << G4endl;
819 
820   GFBeginModeling();
821 
822   G4VSceneHandler::BeginPrimitives (objectTransformation);
823 
824 }
825 
826 
827 //----- 
828 void G4GMocrenFileSceneHandler::EndPrimitives ()
829 {
830   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
831     G4cout << "***** EndPrimitives " << G4endl;
832 
833   G4VSceneHandler::EndPrimitives ();
834 }
835 
836 
837 //========== AddSolid() functions ==========//
838 
839 //----- Add box
840 void G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )
841 {
842   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
843     G4cout << "***** AddSolid ( box )" << G4endl;
844   
845   if(GFDEBUG_DET > 0)
846     G4cout << "G4GMocrenFileSceneHandler::AddSolid(const G4Box&)  : "
847      << box.GetName() << G4endl;
848   
849   //----- skip drawing invisible primitive
850   if( !IsVisible() ) { return ; }
851   
852   //----- Initialize if necessary
853   GFBeginModeling();
854   
855   
856   //--
857   if(GFDEBUG_DET > 1) {
858     G4cout << "-------" << G4endl;
859     G4cout << "    " << box.GetName() << G4endl;
860     G4Polyhedron * poly = box.CreatePolyhedron();
861     poly->Transform(fObjectTransformation);
862     //G4int nv = poly->GetNoVertices();
863     G4Point3D v1, v2;
864     G4int next;
865     //while(1) { // next flag isn't functional.
866     for(G4int i = 0; i < 12; i++) { // # of edges is 12.
867       poly->GetNextEdge(v1, v2, next);
868       if(next == 0) break;
869       G4cout << "    (" << v1.x() << ", "
870        << v1.y() << ", "
871        << v1.z() << ") - ("
872        << v2.x() << ", "
873        << v2.y() << ", "
874        << v2.z() << ") [" << next << "]"
875        << G4endl;
876     }
877     delete poly;
878   }
879   
880   
881   // the volume name set by /vis/gMocren/setVolumeName
882   G4String volName = kMessenger.getVolumeName();
883   
884   
885   if(kFlagParameterization != 2) {
886     G4ScoringManager * pScrMan = G4ScoringManager::GetScoringManager();
887     if(pScrMan) {
888       G4ScoringBox * pScBox = dynamic_cast<G4ScoringBox*>(pScrMan->FindMesh(volName));
889       G4bool bMesh = false;
890       if(pScBox != NULL) bMesh = true;
891       if(bMesh) kFlagParameterization = 2;
892       if(GFDEBUG_DET > 0) G4cout << "   G4ScoringManager::FindMesh() : "
893         << volName << " - " << bMesh << G4endl;
894     }
895   }
896   
897   const G4VModel* pv_model  = GetModel();
898   if (!pv_model) { return ; }
899   G4PhysicalVolumeModel* pPVModel =
900   dynamic_cast<G4PhysicalVolumeModel*>(fpModel);
901   if (!pPVModel) { return ; }
902 
903   
904   //-- debug information
905   if(GFDEBUG_DET > 0) {
906     G4Material * mat = pPVModel->GetCurrentMaterial();
907     G4String name = mat->GetName();
908     G4double dens = mat->GetDensity()/(g/cm3);
909     G4int copyNo = pPVModel->GetCurrentPV()->GetCopyNo();
910     G4int depth = pPVModel->GetCurrentDepth();
911     G4cout << "    copy no.: " << copyNo << G4endl;
912     G4cout << "    depth   : " << depth << G4endl;
913     G4cout << "    density : " << dens << " [g/cm3]" << G4endl;
914     G4cout << "    location: " << pPVModel->GetCurrentPV()->GetObjectTranslation() << G4endl;
915     G4cout << "    Multiplicity        : " << pPVModel->GetCurrentPV()->GetMultiplicity() << G4endl;
916     G4cout << "    Is replicated?      : " << pPVModel->GetCurrentPV()->IsReplicated() << G4endl;
917     G4cout << "    Is parameterised?   : " << pPVModel->GetCurrentPV()->IsParameterised() << G4endl;
918     G4cout << "    top phys. vol. name : " << pPVModel->GetTopPhysicalVolume()->GetName() << G4endl;
919   }
920   
921   //-- check the parameterised volume
922   if(box.GetName() == volName) {
923     
924     kVolumeTrans3D = fObjectTransformation;
925     // coordination system correction for gMocren
926     G4ThreeVector raxis(1., 0., 0.), dummy(0.,0.,0.);
927     G4RotationMatrix rot(raxis, pi*rad);
928     G4Transform3D trot(rot, dummy);
929     if(GFDEBUG_DET) {
930       G4ThreeVector trans1 = kVolumeTrans3D.getTranslation();
931       G4RotationMatrix rot1 = kVolumeTrans3D.getRotation().inverse();
932       G4cout << "kVolumeTrans3D: " << trans1 << G4endl << rot1 << G4endl;
933     }
934     kVolumeTrans3D = kVolumeTrans3D*trot;
935     if(GFDEBUG_DET) G4cout << " Parameterised volume : " << box.GetName() << G4endl;
936     
937     
938     
939     //
940     G4VPhysicalVolume * pv[3] = {0,0,0};
941     pv[0] = pPVModel->GetCurrentPV()->GetLogicalVolume()->GetDaughter(0);
942     if(!pv[0]) {
943       G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
944                   "gMocren0003", FatalException, "Unexpected volume.");
945     }
946     G4int dirAxis[3] = {-1,-1,-1};
947     G4int nDaughters[3] = {0,0,0};
948     
949     EAxis axis; G4int nReplicas; G4double width; G4double offset; G4bool consuming;
950     pv[0]->GetReplicationData(axis, nReplicas, width, offset, consuming);
951     nDaughters[0] = nReplicas;
952     switch(axis) {
953       case kXAxis: dirAxis[0] = 0; break;
954       case kYAxis: dirAxis[0] = 1; break;
955       case kZAxis: dirAxis[0] = 2; break;
956       default:
957         G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
958                     "gMocren0004", FatalException, "Error.");
959     }
960     kNestedVolumeNames.push_back(pv[0]->GetName());
961     if(GFDEBUG_DET)
962       G4cout << "        daughter name :  " << pv[0]->GetName()
963        << "   # : " << nDaughters[0] << G4endl;
964     
965     //
966     if(GFDEBUG_DET) {
967       if(pv[0]->GetLogicalVolume()->GetNoDaughters()) {
968         G4cout << "# of daughters : "
969          << pv[0]->GetLogicalVolume()->GetNoDaughters() << G4endl;
970       } else {
971         //G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
972         //      "gMocren0005", FatalException, "Error.");
973       }
974     }
975     
976     // check whether nested or regular parameterization
977     if(GFDEBUG_DET) G4cout << "# of daughters : "
978          << pv[0]->GetLogicalVolume()->GetNoDaughters() << G4endl;
979     if(pv[0]->GetLogicalVolume()->GetNoDaughters() == 0) {
980       kFlagParameterization = 1;
981       //G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
982       //      "gMocren0006", FatalException, "Error.");
983     }
984     
985     if(kFlagParameterization == 0) {
986       
987       pv[1] = pv[0]->GetLogicalVolume()->GetDaughter(0);
988       if(pv[1]) {
989         pv[1]->GetReplicationData(axis, nReplicas, width, offset, consuming);
990         nDaughters[1] = nReplicas;
991         switch(axis) {
992           case kXAxis: dirAxis[1] = 0; break;
993           case kYAxis: dirAxis[1] = 1; break;
994           case kZAxis: dirAxis[1] = 2; break;
995           default:
996             G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
997                         "gMocren0007", FatalException, "Error.");
998         }
999         kNestedVolumeNames.push_back(pv[1]->GetName());
1000         if(GFDEBUG_DET)
1001           G4cout << "        sub-daughter name :  " << pv[1]->GetName()
1002           << "   # : " << nDaughters[1]<< G4endl;
1003         
1004         //
1005         pv[2] = pv[1]->GetLogicalVolume()->GetDaughter(0);
1006         if(pv[2]) {
1007           nDaughters[2] = pv[2]->GetMultiplicity();
1008           kNestedVolumeNames.push_back(pv[2]->GetName());
1009           if(GFDEBUG_DET)
1010             G4cout << "        sub-sub-daughter name :  " << pv[2]->GetName()
1011             << "   # : " << nDaughters[2] << G4endl;
1012           
1013           if(nDaughters[2] > 1) {
1014             G4VNestedParameterisation * nestPara
1015             = dynamic_cast<G4VNestedParameterisation*>(pv[2]->GetParameterisation());
1016             if(nestPara == NULL)
1017               G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
1018                           "gMocren0008", FatalException, "Non-nested parameterisation");
1019             
1020             nestPara->ComputeTransformation(0, pv[2]);
1021             G4ThreeVector trans0 = pv[2]->GetObjectTranslation();
1022             nestPara->ComputeTransformation(1, pv[2]);
1023             G4ThreeVector trans1 = pv[2]->GetObjectTranslation();
1024             G4ThreeVector diff(trans0 - trans1);
1025             if(GFDEBUG_DET)
1026               G4cout << trans0 << " - " << trans1 << " - " << diff << G4endl;
1027             
1028             if(diff.x() != 0.) dirAxis[2] = 0;
1029             else if(diff.y() != 0.) dirAxis[2] = 1;
1030             else if(diff.z() != 0.) dirAxis[2] = 2;
1031             else
1032               G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
1033                           "gMocren0009", FatalException, "Unexpected nested parameterisation");
1034           }
1035         }
1036       }
1037       
1038       for(G4int i = 0; i < 3; i++) {
1039         kNestedVolumeDimension[i] = nDaughters[i];
1040         //kNestedVolumeDimension[i] = nDaughters[dirAxis[i]];
1041         kNestedVolumeDirAxis[i] = dirAxis[i];
1042       }
1043       //G4cout << "@@@@@@@@@ "
1044       //       << dirAxis[0] << ", " << dirAxis[1] << ", " << dirAxis[2] << G4endl;
1045       
1046       // get densities
1047       G4VNestedParameterisation * nestPara
1048       = dynamic_cast<G4VNestedParameterisation*>(pv[2]->GetParameterisation());
1049       if(nestPara != NULL) {
1050         G4double prexyz[3] = {0.,0.,0.}, xyz[3] = {0.,0.,0.};
1051         for(G4int n0 = 0; n0 < nDaughters[0]; n0++) {
1052           for(G4int n1 = 0; n1 < nDaughters[1]; n1++) {
1053             for(G4int n2 = 0; n2 < nDaughters[2]; n2++) {
1054               
1055               G4GMocrenTouchable * touch = new G4GMocrenTouchable(n1, n0);
1056               if(GFDEBUG_DET)
1057                 G4cout << "   retrieve volume : copy # : " << n0
1058                 << ", " << n1 << ", " << n2 << G4endl;
1059               G4Material * mat = nestPara->ComputeMaterial(pv[2], n2, touch);
1060               delete touch;
1061               G4double dens = mat->GetDensity()/(g/cm3);
1062               
1063               if(GFDEBUG_DET)
1064                 G4cout << "           density :" << dens << " [g/cm3]" << G4endl;
1065               
1066               G4Box tbox(box);
1067               nestPara->ComputeDimensions(tbox, n2, pv[2]);
1068               xyz[0] = tbox.GetXHalfLength()/mm;
1069               xyz[1] = tbox.GetYHalfLength()/mm;
1070               xyz[2] = tbox.GetZHalfLength()/mm;
1071               if(n0 != 0 || n1 != 0 || n2 != 0) {
1072                 for(G4int i = 0; i < 3; i++) {
1073                   if(xyz[i] != prexyz[i])
1074                     G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
1075                                 "gMocren0010", FatalException, "Unsupported parameterisation");
1076                 }
1077               }
1078               if(GFDEBUG_DET)
1079                 G4cout << "              size : " << tbox.GetXHalfLength()/mm << " x "
1080                 << tbox.GetYHalfLength()/mm << " x "
1081                 << tbox.GetZHalfLength()/mm << " [mm3]" << G4endl;
1082               
1083               G4int idx[3];
1084               idx[dirAxis[0]] = n0;
1085               idx[dirAxis[1]] = n1;
1086               idx[dirAxis[2]] = n2;
1087               Index3D i3d(idx[0],idx[1],idx[2]);
1088               kNestedModality[i3d] = dens;
1089               if(GFDEBUG_DET)
1090                 G4cout << " index: " << idx[0] << ", " << idx[1] << ", " << idx[2]
1091                 << "  density: " << dens << G4endl;
1092               
1093               for(G4int i = 0; i < 3; i++) prexyz[i] = xyz[i];
1094             }
1095           }
1096         }
1097         
1098         kVolumeSize.set(box.GetXHalfLength()*2/mm,
1099                         box.GetYHalfLength()*2/mm,
1100                         box.GetZHalfLength()*2/mm);
1101         // mesh size
1102         if(!kbSetModalityVoxelSize) {
1103           G4float spacing[3] = {static_cast<G4float>(2*xyz[0]),
1104             static_cast<G4float>(2*xyz[1]),
1105             static_cast<G4float>(2*xyz[2])};
1106           kgMocrenIO->setVoxelSpacing(spacing);
1107           kVoxelDimension.set(spacing[0], spacing[1], spacing[2]);
1108           kbSetModalityVoxelSize = true;
1109         }
1110         
1111       } else {
1112         if(GFDEBUG_DET)
1113           G4cout << pv[2]->GetName() << G4endl;
1114         G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
1115                     "gMocren0011", FatalException, "Non-nested parameterisation");
1116       }
1117       
1118       
1119       
1120       //-- debug
1121       if(GFDEBUG_DET > 1) {
1122         if(pPVModel->GetCurrentPV()->IsParameterised()) {
1123           G4VPVParameterisation * para = pPVModel->GetCurrentPV()->GetParameterisation();
1124           G4cout << " Is nested parameterisation? : " << para->IsNested() << G4endl;
1125           
1126           
1127           G4int npvp = pPVModel->GetDrawnPVPath().size();
1128           G4cout << "     physical volume node id : "
1129           << "size: " << npvp << ", PV name: ";
1130           for(G4int i = 0; i < npvp; i++) {
1131             G4cout << pPVModel->GetDrawnPVPath()[i].GetPhysicalVolume()->GetName()
1132             << " [param:"
1133             << pPVModel->GetDrawnPVPath()[i].GetPhysicalVolume()->IsParameterised()
1134             << ",rep:"
1135             << pPVModel->GetDrawnPVPath()[i].GetPhysicalVolume()->IsReplicated();
1136             if(pPVModel->GetDrawnPVPath()[i].GetPhysicalVolume()->GetParameterisation()) {
1137               G4cout << ",nest:"
1138               << pPVModel->GetDrawnPVPath()[i].GetPhysicalVolume()->GetParameterisation()->IsNested();
1139             }
1140             G4cout << ",copyno:"
1141             << pPVModel->GetDrawnPVPath()[i].GetPhysicalVolume()->GetCopyNo();
1142             G4cout << "] - ";
1143           }
1144           G4cout << G4endl;
1145           
1146           
1147           pPVModel->GetCurrentPV()->GetReplicationData(axis, nReplicas, width, offset, consuming);
1148           G4cout << "     # replicas : " << nReplicas << G4endl;
1149           G4double pareDims[3] = {0.,0.,0.};
1150           G4Box * pbox = dynamic_cast<G4Box *>(pPVModel->GetDrawnPVPath()[npvp-2].GetPhysicalVolume()->GetLogicalVolume()->GetSolid());
1151           if(pbox) {
1152             pareDims[0] = 2.*pbox->GetXHalfLength()*mm;
1153             pareDims[1] = 2.*pbox->GetYHalfLength()*mm;
1154             pareDims[2] = 2.*pbox->GetZHalfLength()*mm;
1155             G4cout << "     mother size ["
1156             << pPVModel->GetDrawnPVPath()[npvp-2].GetPhysicalVolume()->GetName()
1157             << "] : "
1158             << pareDims[0] << " x "
1159             << pareDims[1] << " x "
1160             << pareDims[2] << " [mm3]"
1161             << G4endl;
1162           }
1163           G4double paraDims[3];
1164           G4Box * boxP = dynamic_cast<G4Box *>(pPVModel->GetDrawnPVPath()[npvp-1].GetPhysicalVolume()->GetLogicalVolume()->GetSolid());
1165           if(boxP) {
1166             paraDims[0] = 2.*boxP->GetXHalfLength()*mm;
1167             paraDims[1] = 2.*boxP->GetYHalfLength()*mm;
1168             paraDims[2] = 2.*boxP->GetZHalfLength()*mm;
1169             G4cout << "     parameterised volume? ["
1170             << pPVModel->GetDrawnPVPath()[npvp-1].GetPhysicalVolume()->GetName()
1171             << "] : "
1172             << paraDims[0] << " x "
1173             << paraDims[1] << " x "
1174             << paraDims[2] << " [mm3]  : "
1175             << G4int(pareDims[0]/paraDims[0]) << " x "
1176             << G4int(pareDims[1]/paraDims[1]) << " x "
1177             << G4int(pareDims[2]/paraDims[2]) << G4endl;
1178           } else {
1179             G4cout << pPVModel->GetDrawnPVPath()[npvp-2].GetPhysicalVolume()->GetName()
1180             << " isn't a G4Box." << G4endl;
1181           }
1182         }
1183       }
1184       
1185       
1186     } else if(kFlagParameterization == 1) { // G4PhantomParameterisation based geom. construnction
1187       
1188       // get the dimension of the parameterized patient geometry
1189       G4PhantomParameterisation * phantomPara
1190       = dynamic_cast<G4PhantomParameterisation*>(pv[0]->GetParameterisation());
1191       if(phantomPara == NULL) {
1192         G4Exception("G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )",
1193                     "gMocren0012", FatalException, "no G4PhantomParameterisation");
1194       } else {
1195         ;
1196       }
1197       
1198       kNestedVolumeDimension[0] = (G4int)phantomPara->GetNoVoxelsX();
1199       kNestedVolumeDimension[1] = (G4int)phantomPara->GetNoVoxelsY();
1200       kNestedVolumeDimension[2] = (G4int)phantomPara->GetNoVoxelsZ();
1201       kNestedVolumeDirAxis[0] = 0;
1202       kNestedVolumeDirAxis[1] = 1;
1203       kNestedVolumeDirAxis[2] = 2;
1204       
1205       // get densities of the parameterized patient geometry
1206       G4int nX = kNestedVolumeDimension[0];
1207       G4int nXY = kNestedVolumeDimension[0]*kNestedVolumeDimension[1];
1208       
1209       for(G4int n0 = 0; n0 < kNestedVolumeDimension[0]; n0++) {
1210         for(G4int n1 = 0; n1 < kNestedVolumeDimension[1]; n1++) {
1211           for(G4int n2 = 0; n2 < kNestedVolumeDimension[2]; n2++) {
1212             
1213             G4int repNo = n0 + n1*nX + n2*nXY;
1214             G4Material * mat = phantomPara->ComputeMaterial(repNo, pv[0]);
1215             G4double dens = mat->GetDensity()/(g/cm3);
1216             
1217             
1218             G4int idx[3];
1219             idx[kNestedVolumeDirAxis[0]] = n0;
1220             idx[kNestedVolumeDirAxis[1]] = n1;
1221             idx[kNestedVolumeDirAxis[2]] = n2;
1222             Index3D i3d(idx[0],idx[1],idx[2]);
1223             kNestedModality[i3d] = dens;
1224             
1225             if(GFDEBUG_DET)
1226               G4cout << " index: " << idx[0] << ", " << idx[1] << ", " << idx[2]
1227               << "  density: " << dens << G4endl;
1228             
1229           }
1230         }
1231       }
1232       
1233       kVolumeSize.set(box.GetXHalfLength()*2/mm,
1234                       box.GetYHalfLength()*2/mm,
1235                       box.GetZHalfLength()*2/mm);
1236       
1237       // mesh size
1238       if(!kbSetModalityVoxelSize) {
1239         G4float spacing[3] = {static_cast<G4float>(2*phantomPara->GetVoxelHalfX()),
1240           static_cast<G4float>(2*phantomPara->GetVoxelHalfY()),
1241           static_cast<G4float>(2*phantomPara->GetVoxelHalfZ())};
1242         kgMocrenIO->setVoxelSpacing(spacing);
1243         kVoxelDimension.set(spacing[0], spacing[1], spacing[2]);
1244         kbSetModalityVoxelSize = true;
1245       }
1246     }
1247     
1248   } // if(box.GetName() == volName)
1249   
1250   
1251   // processing geometry construction based on the interactive PS
1252   if(!kFlagProcessedInteractiveScorer) {
1253     
1254     
1255     // get the dimension of the geometry defined in G4VScoringMesh
1256     G4ScoringManager * pScrMan = G4ScoringManager::GetScoringManager();
1257     //if(!pScrMan) return;
1258     if(pScrMan) {
1259     G4ScoringBox * scoringBox
1260     = dynamic_cast<G4ScoringBox*>(pScrMan->FindMesh(volName));
1261     //if(scoringBox == NULL) return;
1262     if(scoringBox) {
1263 
1264     
1265     
1266     G4int nVoxels[3];
1267     scoringBox->GetNumberOfSegments(nVoxels);
1268     // this order depends on the G4ScoringBox
1269     kNestedVolumeDimension[0] = nVoxels[2];
1270     kNestedVolumeDimension[1] = nVoxels[1];
1271     kNestedVolumeDimension[2] = nVoxels[0];
1272     kNestedVolumeDirAxis[0] = 2;
1273     kNestedVolumeDirAxis[1] = 1;
1274     kNestedVolumeDirAxis[2] = 0;
1275     
1276     // get densities of the parameterized patient geometry
1277     for(G4int n0 = 0; n0 < kNestedVolumeDimension[0]; n0++) {
1278       for(G4int n1 = 0; n1 < kNestedVolumeDimension[1]; n1++) {
1279         for(G4int n2 = 0; n2 < kNestedVolumeDimension[2]; n2++) {
1280           
1281           G4double dens = 0.*(g/cm3);
1282           
1283           G4int idx[3];
1284           idx[kNestedVolumeDirAxis[0]] = n0;
1285           idx[kNestedVolumeDirAxis[1]] = n1;
1286           idx[kNestedVolumeDirAxis[2]] = n2;
1287           Index3D i3d(idx[0],idx[1],idx[2]);
1288           kNestedModality[i3d] = dens;
1289           
1290         }
1291       }
1292     }
1293     
1294     G4ThreeVector boxSize = scoringBox->GetSize();
1295     if(GFDEBUG_DET > 1) {
1296       G4cout << "Interactive Scorer : size - "
1297        << boxSize.x()/cm << " x "
1298        << boxSize.y()/cm << " x "
1299        << boxSize.z()/cm << " [cm3]" << G4endl;
1300       G4cout << "Interactive Scorer : # voxels - "
1301        << nVoxels[0] << " x "
1302        << nVoxels[1] << " x "
1303        << nVoxels[2] << G4endl;
1304     }
1305     kVolumeSize.set(boxSize.x()*2,
1306                     boxSize.y()*2,
1307                     boxSize.z()*2);
1308     
1309     // mesh size
1310     if(!kbSetModalityVoxelSize) {
1311       G4float spacing[3] = {static_cast<G4float>(boxSize.x()*2/nVoxels[0]),
1312         static_cast<G4float>(boxSize.y()*2/nVoxels[1]),
1313         static_cast<G4float>(boxSize.z()*2/nVoxels[2])};
1314       
1315       kgMocrenIO->setVoxelSpacing(spacing);
1316       kVoxelDimension.set(spacing[0], spacing[1], spacing[2]);
1317       kbSetModalityVoxelSize = true;
1318       
1319     }
1320     
1321     
1322     kVolumeTrans3D = fObjectTransformation;
1323     
1324     // translation for the scoring mesh
1325     G4ThreeVector sbth = scoringBox->GetTranslation();
1326     G4Translate3D sbtranslate(sbth);
1327     kVolumeTrans3D = kVolumeTrans3D*sbtranslate;
1328     
1329     // rotation matrix for the scoring mesh 
1330     G4RotationMatrix sbrm;
1331     sbrm = scoringBox->GetRotationMatrix();
1332     if(!sbrm.isIdentity()) {
1333       G4ThreeVector sbdummy(0.,0.,0.); 
1334       G4Transform3D sbrotate(sbrm.inverse(), sbdummy);
1335       kVolumeTrans3D = kVolumeTrans3D*sbrotate;
1336     }
1337     
1338     
1339     // coordination system correction for gMocren
1340     G4ThreeVector raxisY(0., 1., 0.), dummyY(0.,0.,0.); 
1341     G4RotationMatrix rotY(raxisY, pi*rad);
1342     G4Transform3D trotY(rotY, dummyY);
1343     G4ThreeVector raxisZ(0., 0., 1.), dummyZ(0.,0.,0.); 
1344     G4RotationMatrix rotZ(raxisZ, pi*rad);
1345     G4Transform3D trotZ(rotZ, dummyZ);
1346     
1347     kVolumeTrans3D = kVolumeTrans3D*trotY*trotZ;
1348     
1349     }
1350     }
1351     //
1352     kFlagProcessedInteractiveScorer = true;
1353   }
1354   
1355   
1356   static G4VPhysicalVolume * volPV = NULL;
1357   if(pPVModel->GetCurrentPV()->GetName() == volName) {
1358     volPV = pPVModel->GetCurrentPV();
1359   }
1360   
1361   //-- add detectors
1362   G4bool bAddDet = true;
1363   if(!kMessenger.getDrawVolumeGrid()) {
1364 
1365     if(kFlagParameterization == 0) { // nested parameterisation
1366       /*
1367       G4String volDSolidName;
1368       if(volPV) {
1369         G4int nDaughter = volPV->GetLogicalVolume()->GetNoDaughters();
1370         G4VPhysicalVolume * volDPV = NULL;
1371         if(nDaughter > 0) volDPV = volPV->GetLogicalVolume()->GetDaughter(0);
1372         if(volDPV) {
1373           nDaughter = volDPV->GetLogicalVolume()->GetNoDaughters();
1374           if(nDaughter > 0)
1375             volDSolidName = volDPV->GetLogicalVolume()->GetDaughter(0)
1376                            ->GetLogicalVolume()->GetSolid()->GetName();
1377         }
1378       }
1379       */
1380       
1381       //std::cout << "Parameterization volume: " << volName << " - "
1382       //          << box.GetName() << std::endl;
1383       
1384       if(volName == box.GetName()) {
1385         bAddDet = false;
1386       }
1387       
1388       std::vector<G4String>::iterator itr = kNestedVolumeNames.begin();
1389       for(; itr != kNestedVolumeNames.end(); itr++) {
1390         if(*itr == box.GetName())  {
1391           bAddDet = false;
1392           break;
1393         }
1394       }
1395     } else if(kFlagParameterization == 1) { // phantom paramemterisation
1396       
1397       G4String volDSolidName;
1398       if(volPV) {
1399         volDSolidName = volPV->GetLogicalVolume()->GetDaughter(0)
1400                         ->GetLogicalVolume()->GetSolid()->GetName();
1401       }
1402 
1403       //std::cout << "Phantom Parameterization volume: " << volDSolidName
1404       //          << " - " << box.GetName() << std::endl;
1405       
1406       if(volDSolidName == box.GetName()) {
1407         bAddDet = false;
1408       }
1409       
1410     } else if(kFlagParameterization == 2) { // interactive primitive scorer
1411       //std::cout << "Regular Parameterization volume: " << box.GetName() << std::endl;
1412     }
1413     
1414   }
1415   if(bAddDet) AddDetector(box);
1416   
1417 
1418 } // void G4GMocrenFileSceneHandler::AddSolid( const G4Box& box )
1419 
1420 
1421 //----- Add tubes
1422 void 
1423 G4GMocrenFileSceneHandler::AddSolid( const G4Tubs& tubes )
1424 {
1425   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
1426     G4cout << "***** AddSolid ( tubes )" << G4endl;
1427 
1428   //----- skip drawing invisible primitive
1429   if( !IsVisible() ) { return ; }
1430 
1431   //----- Initialize if necessary
1432   GFBeginModeling();
1433 
1434   //
1435   AddDetector(tubes);
1436 
1437 
1438   // for a debug
1439   if(GFDEBUG_DET > 0) {
1440     G4cout << "-------" << G4endl;
1441     G4cout << "    " << tubes.GetName() << G4endl;
1442     G4Polyhedron * poly = tubes.CreatePolyhedron();
1443     G4int nv = poly->GetNoVertices();
1444     for(G4int i = 0; i < nv; i++) {
1445       G4cout << "    (" << poly->GetVertex(i).x() << ", "
1446        << poly->GetVertex(i).y() << ", "
1447        << poly->GetVertex(i).z() << ")" << G4endl;
1448     }
1449     delete poly;
1450   }
1451 
1452   const G4VModel* pv_model  = GetModel();
1453   if (!pv_model) { return ; } 
1454   G4PhysicalVolumeModel* pPVModel =
1455     dynamic_cast<G4PhysicalVolumeModel*>(fpModel);
1456   if (!pPVModel) { return ; }
1457   G4Material * mat = pPVModel->GetCurrentMaterial();
1458   G4String name = mat->GetName();
1459 
1460 } // void G4GMocrenFileSceneHandler::AddSolid( const G4Tubs& )
1461 
1462 
1463 
1464 //----- Add cons
1465 void 
1466 G4GMocrenFileSceneHandler::AddSolid( const G4Cons& cons )
1467 {
1468   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
1469     G4cout << "***** AddSolid ( cons )" << G4endl;
1470 
1471   //----- skip drawing invisible primitive
1472   if( !IsVisible() ) { return ; }
1473 
1474   //----- Initialize if necessary
1475   GFBeginModeling();
1476 
1477   //
1478   AddDetector(cons);
1479 
1480 }// G4GMocrenFileSceneHandler::AddSolid( cons )
1481 
1482 
1483 //----- Add trd
1484 void G4GMocrenFileSceneHandler::AddSolid ( const G4Trd& trd )
1485 {
1486   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
1487     G4cout << "***** AddSolid ( trd )" << G4endl;
1488 
1489 
1490   //----- skip drawing invisible primitive
1491   if( !IsVisible() ) { return ; }
1492 
1493   //----- Initialize if necessary
1494   GFBeginModeling();
1495 
1496   //
1497   AddDetector(trd);
1498 
1499 } // G4GMocrenFileSceneHandler::AddSolid ( trd )
1500 
1501 
1502 //----- Add sphere
1503 void G4GMocrenFileSceneHandler::AddSolid ( const G4Sphere& sphere )
1504 {
1505   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
1506     G4cout << "***** AddSolid ( sphere )" << G4endl;
1507 
1508   //----- skip drawing invisible primitive
1509   if( !IsVisible() ) { return ; }
1510 
1511   //----- Initialize if necessary
1512   GFBeginModeling();
1513 
1514   //
1515   AddDetector(sphere);
1516 
1517 } // G4GMocrenFileSceneHandler::AddSolid ( sphere )
1518 
1519 
1520 //----- Add para
1521 void G4GMocrenFileSceneHandler::AddSolid (const G4Para& para)
1522 {
1523   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
1524     G4cout << "***** AddSolid ( para )" << G4endl;
1525 
1526   //----- skip drawing invisible primitive
1527   if( !IsVisible() ) { return ; }
1528 
1529   //----- Initialize if necessary
1530   GFBeginModeling();
1531 
1532   //
1533   AddDetector(para);
1534 
1535 } // G4GMocrenFileSceneHandler::AddSolid ( para )
1536 
1537 
1538 //----- Add trap
1539 void G4GMocrenFileSceneHandler::AddSolid (const G4Trap& trap)
1540 {
1541   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
1542     G4cout << "***** AddSolid ( trap )" << G4endl;
1543 
1544   //----- skip drawing invisible primitive
1545   if( !IsVisible() ) { return ; }
1546 
1547   //----- Initialize if necessary
1548   GFBeginModeling();
1549 
1550   //
1551   AddDetector(trap);
1552 
1553 } // G4GMocrenFileSceneHandler::AddSolid (const G4Trap& trap)
1554 
1555 
1556 //----- Add torus
1557 void 
1558 G4GMocrenFileSceneHandler::AddSolid( const G4Torus& torus )
1559 {
1560   if(GFDEBUG || G4VisManager::GetVerbosity() >= G4VisManager::confirmations)
1561     G4cout << "***** AddSolid ( torus )" << G4endl;
1562 
1563   //----- skip drawing invisible primitive
1564   if( !IsVisible() ) { return ; }
1565 
1566   //----- Initialize if necessary
1567   GFBeginModeling();
1568 
1569   //
1570   AddDetector(torus);
1571 
1572 } // void G4GMocrenFileSceneHandler::AddSolid( const G4Torus& )
1573 
1574 
1575 
1576 //----- Add a shape which is not treated above
1577 void G4GMocrenFileSceneHandler::AddSolid ( const G4VSolid& solid  )
1578 {
1579   //----- skip drawing invisible primitive
1580   if( !IsVisible() ) { return ; }
1581 
1582   //----- Initialize if necessary
1583   GFBeginModeling();
1584 
1585   //
1586   AddDetector(solid);
1587 
1588   //----- Send a primitive
1589   G4VSceneHandler::AddSolid( solid ) ; 
1590 
1591 } //G4GMocrenFileSceneHandler::AddSolid ( const G4VSolid& ) 
1592 
1593 #include "G4TrajectoriesModel.hh"
1594 #include "G4VTrajectory.hh"
1595 #include "G4VTrajectoryPoint.hh"
1596 
1597 //----- Add a trajectory
1598 void G4GMocrenFileSceneHandler::AddCompound(const G4VTrajectory & traj) {
1599 
1600   kbModelingTrajectory = true;
1601 
1602   G4VSceneHandler::AddCompound(traj);
1603 
1604   if(GFDEBUG_TRK) {
1605     G4cout << " ::AddCompound(const G4VTrajectory&) >>>>>>>>> " << G4endl;
1606     G4TrajectoriesModel * pTrModel = dynamic_cast<G4TrajectoriesModel*>(fpModel);
1607     if (!pTrModel) { 
1608       G4Exception 
1609   ("G4VSceneHandler::AddCompound(const G4VTrajectory&)",
1610    "gMocren0013", FatalException, "Not a G4TrajectoriesModel.");
1611     } else {
1612       traj.DrawTrajectory();
1613 
1614       const G4VTrajectory * trj = pTrModel->GetCurrentTrajectory();
1615       G4cout << "------ track" << G4endl;
1616       G4cout << "    name:     " << trj->GetParticleName() << G4endl;
1617       G4cout << "    id:       " << trj->GetTrackID() << G4endl;
1618       G4cout << "    charge:   " << trj->GetCharge() << G4endl;
1619       G4cout << "    momentum: " << trj->GetInitialMomentum() << G4endl;
1620       
1621       G4int nPnt = trj->GetPointEntries();
1622       G4cout << "    point:    ";
1623       for(G4int i = 0; i < nPnt; i++) {
1624   G4cout << trj->GetPoint(i)->GetPosition() << ", ";
1625       }
1626       G4cout << G4endl;
1627     }
1628     G4cout << G4endl;
1629   }
1630 
1631   kbModelingTrajectory = false;
1632 }
1633 
1634 #include <vector>
1635 #include "G4VHit.hh"
1636 #include "G4AttValue.hh"
1637 //----- Add a hit 
1638 void G4GMocrenFileSceneHandler::AddCompound( const G4VHit & hit) {
1639   if(GFDEBUG_HIT) G4cout << " ::AddCompound(const G4VHit&) >>>>>>>>> " << G4endl;
1640 
1641   G4VSceneHandler::AddCompound(hit);
1642 
1643   /*
1644     const std::map<G4String, G4AttDef> * map = hit.GetAttDefs();
1645     if(!map) return;
1646     std::map<G4String, G4AttDef>::const_iterator itr = map->begin();
1647     for(; itr != map->end(); itr++) {
1648     G4cout << itr->first << " : " << itr->second.GetName()
1649     << " , " << itr->second.GetDesc() << G4endl;
1650     }
1651   */
1652 
1653   std::vector<G4String> hitNames = kMessenger.getHitNames();
1654   if(GFDEBUG_HIT) {
1655     std::vector<G4String>::iterator itr = hitNames.begin();
1656     for(; itr != hitNames.end(); itr++) 
1657       G4cout << "  hit name : " << *itr << G4endl;
1658   }
1659   
1660   std::vector<G4AttValue> * attval = hit.CreateAttValues();
1661   if(!attval) {G4cout << "0 empty " << G4endl;}
1662   else {
1663 
1664     G4bool bid[3] = {false, false, false};
1665     Index3D id;
1666 
1667     std::vector<G4AttValue>::iterator itr;
1668     // First, get IDs
1669     for(itr = attval->begin(); itr != attval->end(); itr++) {
1670       std::string stmp = itr->GetValue();
1671       std::istringstream sval(stmp.c_str());
1672 
1673       if(itr->GetName() == G4String("XID")) {
1674   sval >> id.x;
1675   bid[0] = true;
1676   continue;
1677       }
1678       if(itr->GetName() == G4String("YID")) {
1679   sval >> id.y;
1680   bid[1] = true;
1681   continue;
1682       }
1683       if(itr->GetName() == G4String("ZID")) {
1684   sval >> id.z;
1685   bid[2] = true;
1686   continue;
1687       }
1688     }
1689 
1690     G4int nhitname = (G4int)hitNames.size();
1691 
1692     if(bid[0] && bid[1] && bid[2]) {
1693 
1694       if(GFDEBUG_HIT)
1695   G4cout << " Hit : index(" << id.x << ", " << id.y << ", "
1696          << id.z << ")" << G4endl;
1697 
1698       // Get attributes
1699       for(itr = attval->begin(); itr != attval->end(); itr++) {
1700   for(G4int i = 0; i < nhitname; i++) {
1701     if(itr->GetName() == hitNames[i]) {
1702 
1703       std::string stmp = itr->GetValue();
1704       std::istringstream sval(stmp.c_str());
1705       G4double value;
1706       G4String unit;
1707       sval >> value >> unit;
1708 
1709       std::map<G4String, std::map<Index3D, G4double> >::iterator kNestedHitsListItr;
1710       kNestedHitsListItr = kNestedHitsList.find(hitNames[i]);
1711       if(kNestedHitsListItr != kNestedHitsList.end()) {
1712         //fTempNestedHits = &kNestedHitsListItr->second;
1713         //(*fTempNestedHits)[id] = value;
1714         kNestedHitsListItr->second[id] = value;
1715       } else {
1716         std::map<Index3D, G4double> hits;
1717         hits.insert(std::map<Index3D, G4double>::value_type(id, value));
1718         kNestedHitsList[hitNames[i]] = std::move(hits);
1719       }
1720 
1721       
1722       if(GFDEBUG_HIT)
1723         G4cout << "     : " << hitNames[i] << " -> " << value
1724          << " [" << unit << "]" << G4endl;
1725     }
1726   }
1727       }
1728     } else {
1729       G4Exception("G4GMocrenFileSceneHandler::AddCompound(const G4VHit &)",
1730       "gMocren0014", FatalException, "Error");
1731     }
1732 
1733     delete attval;
1734   }
1735 
1736 }
1737 
1738 void G4GMocrenFileSceneHandler::AddCompound( const G4VDigi & digi) {
1739   if(GFDEBUG_DIGI) G4cout << " ::AddCompound(const G4VDigi&) >>>>>>>>> " << G4endl;
1740   G4VSceneHandler::AddCompound(digi);
1741 }
1742 
1743 void G4GMocrenFileSceneHandler::AddCompound(const G4THitsMap<G4double> & hits) {
1744   if(GFDEBUG_HIT)
1745     G4cout << " ::AddCompound(const std::map<G4int, G4double*> &) >>>>>>>>> " << G4endl;
1746 
1747 
1748   std::vector<G4String> hitScorerNames = kMessenger.getHitScorerNames();
1749   G4int nhitname = (G4int)hitScorerNames.size();
1750   G4String scorername = static_cast<G4VHitsCollection>(hits).GetName();
1751 
1752   //-- --//
1753   /*
1754   std::vector<G4String> hitScorerNames = kMessenger.getHitScorerNames();
1755   if(GFDEBUG_HIT) {
1756     std::vector<G4String>::iterator itr = hitScorerNames.begin();
1757     for(; itr != hitScorerNames.end(); itr++) 
1758       G4cout << "  PS name : " << *itr << G4endl;
1759   }
1760   */
1761   
1762   {  // Scope bracket to avoid compiler messages about shadowing (JA).
1763   //for(G4int i = 0; i < nhitname; i++) {       // this selection trusts
1764     //if(scorername == hitScorerNames[i]) {   // thea command /vis/scene/add/psHits hit_name.
1765 
1766       G4int idx[3];
1767       std::map<G4int, G4double*> * map = hits.GetMap();
1768       std::map<G4int, G4double*>::const_iterator itr = map->begin();
1769       for(; itr != map->end(); itr++) {
1770   GetNestedVolumeIndex(itr->first, idx);
1771   Index3D id(idx[0], idx[1], idx[2]);
1772   
1773   std::map<G4String, std::map<Index3D, G4double> >::iterator nestedHitsListItr;
1774   nestedHitsListItr = kNestedHitsList.find(scorername);
1775   if(nestedHitsListItr != kNestedHitsList.end()) {
1776     nestedHitsListItr->second[id] = *(itr->second);
1777   } else {
1778     std::map<Index3D, G4double> hit;
1779     hit.insert(std::map<Index3D, G4double>::value_type(id, *(itr->second)));
1780     kNestedHitsList[scorername] = std::move(hit);
1781   }
1782       }
1783  
1784       //break;
1785     //}
1786   //}
1787   }
1788 
1789   if(GFDEBUG_HIT) {
1790     G4String meshname = static_cast<G4VHitsCollection>(hits).GetSDname();
1791     G4cout << "       >>>>> " << meshname << " : " << scorername  << G4endl;
1792 
1793     for(G4int i = 0; i < nhitname; i++)
1794       if(scorername == hitScorerNames[i]) 
1795   G4cout << "       !!!! Hit scorer !!!! " << scorername << G4endl;
1796 
1797     G4cout << " dimension: "
1798      << kNestedVolumeDimension[0] << " x "
1799      << kNestedVolumeDimension[1] << " x "
1800      << kNestedVolumeDimension[2] << G4endl;
1801 
1802     G4int id[3];
1803     std::map<G4int, G4double*> * map = hits.GetMap();
1804     std::map<G4int, G4double*>::const_iterator itr = map->begin();
1805     for(; itr != map->end(); itr++) {
1806       GetNestedVolumeIndex(itr->first, id);
1807       G4cout << "[" << itr->first << "] "
1808        << "("<< id[0] << "," << id[1] << "," << id[2] << ")"
1809        << *(itr->second) << ", ";
1810     }
1811     G4cout << G4endl;
1812   }
1813 }
1814 
1815 void G4GMocrenFileSceneHandler::AddCompound(const G4THitsMap<G4StatDouble> & hits) {
1816   if(GFDEBUG_HIT)
1817     G4cout << " ::AddCompound(const std::map<G4int, G4StatDouble*> &) >>>>>>>>> " << G4endl;
1818 
1819 
1820   std::vector<G4String> hitScorerNames = kMessenger.getHitScorerNames();
1821   G4int nhitname = (G4int)hitScorerNames.size();
1822   G4String scorername = static_cast<G4VHitsCollection>(hits).GetName();
1823 
1824   //-- --//
1825   /*
1826   std::vector<G4String> hitScorerNames = kMessenger.getHitScorerNames();
1827   if(GFDEBUG_HIT) {
1828     std::vector<G4String>::iterator itr = hitScorerNames.begin();
1829     for(; itr != hitScorerNames.end(); itr++) 
1830       G4cout << "  PS name : " << *itr << G4endl;
1831   }
1832   */
1833   
1834   {  // Scope bracket to avoid compiler messages about shadowing (JA).
1835   //for(G4int i = 0; i < nhitname; i++) {       // this selection trusts
1836     //if(scorername == hitScorerNames[i]) {   // thea command /vis/scene/add/psHits hit_name.
1837 
1838       G4int idx[3];
1839       std::map<G4int, G4StatDouble*> * map = hits.GetMap();
1840       std::map<G4int, G4StatDouble*>::const_iterator itr = map->begin();
1841       for(; itr != map->end(); itr++) {
1842   GetNestedVolumeIndex(itr->first, idx);
1843   Index3D id(idx[0], idx[1], idx[2]);
1844   
1845   std::map<G4String, std::map<Index3D, G4double> >::iterator nestedHitsListItr;
1846   nestedHitsListItr = kNestedHitsList.find(scorername);
1847   if(nestedHitsListItr != kNestedHitsList.end()) {
1848     nestedHitsListItr->second[id] = itr->second->sum_wx();
1849   } else {
1850     std::map<Index3D, G4double> hit;
1851     hit.insert(std::map<Index3D, G4double>::value_type(id, itr->second->sum_wx()));
1852     kNestedHitsList[scorername] = std::move(hit);
1853   }
1854       }
1855  
1856       //break;
1857     //}
1858   //}
1859   }
1860 
1861   if(GFDEBUG_HIT) {
1862     G4String meshname = static_cast<G4VHitsCollection>(hits).GetSDname();
1863     G4cout << "       >>>>> " << meshname << " : " << scorername  << G4endl;
1864 
1865     for(G4int i = 0; i < nhitname; i++)
1866       if(scorername == hitScorerNames[i]) 
1867   G4cout << "       !!!! Hit scorer !!!! " << scorername << G4endl;
1868 
1869     G4cout << " dimension: "
1870      << kNestedVolumeDimension[0] << " x "
1871      << kNestedVolumeDimension[1] << " x "
1872      << kNestedVolumeDimension[2] << G4endl;
1873 
1874     G4int id[3];
1875     std::map<G4int, G4StatDouble*> * map = hits.GetMap();
1876     std::map<G4int, G4StatDouble*>::const_iterator itr = map->begin();
1877     for(; itr != map->end(); itr++) {
1878       GetNestedVolumeIndex(itr->first, id);
1879       G4cout << "[" << itr->first << "] "
1880        << "("<< id[0] << "," << id[1] << "," << id[2] << ")"
1881        << itr->second->sum_wx() << ", ";
1882     }
1883     G4cout << G4endl;
1884   }
1885 }
1886 
1887 //----- 
1888 G4bool G4GMocrenFileSceneHandler::IsVisible()
1889 {
1890   //----- 
1891   G4bool  visibility  = true ;
1892 
1893   const G4VisAttributes* pVisAttribs =
1894     fpViewer->GetApplicableVisAttributes( fpVisAttribs );
1895 
1896   if(pVisAttribs) {
1897     visibility = pVisAttribs->IsVisible();
1898   } 
1899 
1900   return visibility ;
1901 
1902 } // G4GMocrenFileSceneHandler::IsVisible()
1903 
1904 
1905 //----- 
1906 void G4GMocrenFileSceneHandler::ClearTransientStore() 
1907 {
1908   // This is typically called after an update and before drawing hits
1909   // of the next event.  To simulate the clearing of "transients"
1910   // (hits, etc.) the detector is redrawn...
1911   if (fpViewer) {
1912     fpViewer -> SetView ();
1913     fpViewer -> ClearView ();
1914     fpViewer -> DrawView ();
1915   }
1916 }
1917 
1918 //----- 
1919 void G4GMocrenFileSceneHandler::AddDetector(const G4VSolid & solid) {
1920 
1921   Detector detector;
1922 
1923   // detector name
1924   detector.name = solid.GetName();
1925   if(GFDEBUG_DET > 1)
1926     G4cout << "0 Detector name : " << detector.name << G4endl;
1927  
1928   const G4VModel* pv_model  = GetModel();
1929   if (!pv_model) { return ; } 
1930   G4PhysicalVolumeModel* pPVModel =
1931     dynamic_cast<G4PhysicalVolumeModel*>(fpModel);
1932   if (!pPVModel) { return ; }
1933 
1934   // edge points of the detector
1935   std::vector<G4float *> dedges;
1936   G4Polyhedron * poly = solid.CreatePolyhedron();
1937   detector.polyhedron = poly;
1938   detector.transform3D = fObjectTransformation;
1939 
1940   // retrieve color
1941   unsigned char uccolor[3] = {30, 30, 30};
1942   if(pPVModel->GetCurrentLV()->GetVisAttributes()) {
1943     G4Color color = pPVModel->GetCurrentLV()->GetVisAttributes()->GetColor();
1944     uccolor[0] = (unsigned char)(color.GetRed()*255);
1945     uccolor[1] = (unsigned char)(color.GetGreen()*255);
1946     uccolor[2] = (unsigned char)(color.GetBlue()*255);
1947     //if(uccolor[0] < 2 && uccolor[1] < 2 && uccolor[2] < 2)
1948     //uccolor[0] = uccolor[1] = uccolor[2] = 30; // dark grey
1949   }
1950   for(G4int i = 0; i < 3; i++) detector.color[i] = uccolor[i];
1951   //
1952   kDetectors.push_back(detector);
1953 
1954   if(GFDEBUG_DET > 1) {
1955     G4cout << "0     color:   (" << (G4int)uccolor[0] << ", "
1956      << (G4int)uccolor[1] << ", " << (G4int)uccolor[2] << ")"
1957      << G4endl;
1958   }
1959 
1960 }
1961 
1962 //----- 
1963 void G4GMocrenFileSceneHandler::ExtractDetector() {
1964 
1965   std::vector<Detector>::iterator itr = kDetectors.begin();
1966 
1967   for(; itr != kDetectors.end(); itr++) {
1968 
1969     // detector name
1970     G4String detname = itr->name;
1971     if(GFDEBUG_DET > 1)
1972       G4cout << "Detector name : " << detname << G4endl;
1973 
1974     // edge points of the detector
1975     std::vector<G4float *> dedges;
1976     G4Polyhedron * poly = itr->polyhedron;
1977     poly->Transform(itr->transform3D);
1978     G4Transform3D invVolTrans = kVolumeTrans3D.inverse();
1979     poly->Transform(invVolTrans);
1980 
1981     G4Point3D v1, v2;
1982     G4bool bnext = true;
1983     G4int next;
1984     G4int nedges = 0;
1985     //
1986     while(bnext) {
1987       if(!(poly->GetNextEdge(v1, v2, next))) bnext =false;
1988       G4float * edge = new G4float[6];
1989       edge[0] = v1.x()/mm;
1990       edge[1] = v1.y()/mm;
1991       edge[2] = v1.z()/mm;
1992       edge[3] = v2.x()/mm;
1993       edge[4] = v2.y()/mm;
1994       edge[5] = v2.z()/mm;
1995       dedges.push_back(edge);
1996       nedges++;
1997     }
1998     //delete poly;
1999     // detector color
2000     unsigned char uccolor[3] = {itr->color[0],
2001         itr->color[1],
2002         itr->color[2]};
2003     //
2004     kgMocrenIO->addDetector(detname, dedges, uccolor);
2005     for(G4int i = 0; i < nedges; i++) { // # of edges is 12.
2006       delete [] dedges[i];
2007     }
2008     dedges.clear(); 
2009 
2010     if(GFDEBUG_DET > 1) {
2011       G4cout << "    color:   (" << (G4int)uccolor[0] << ", "
2012        << (G4int)uccolor[1] << ", " << (G4int)uccolor[2] << ")"
2013        << G4endl;
2014     }
2015   }
2016 }
2017 
2018 void G4GMocrenFileSceneHandler::GetNestedVolumeIndex(G4int _idx, G4int _idx3d[3]) {
2019   if(kNestedVolumeDimension[0] == 0 ||
2020      kNestedVolumeDimension[1] == 0 ||
2021      kNestedVolumeDimension[2] == 0) {
2022     for(G4int i = 0; i < 3; i++) _idx3d[i] = 0;
2023     return;
2024   }
2025 
2026 
2027   if(kFlagParameterization == 0) {
2028 
2029     G4int plane = kNestedVolumeDimension[2]*kNestedVolumeDimension[1];
2030     G4int line = kNestedVolumeDimension[2];
2031  
2032   /*
2033   G4int idx3d[3];
2034   idx3d[0] = _idx/plane;
2035   idx3d[1] = (_idx%plane)/line;
2036   idx3d[2] = (_idx%plane)%line;
2037   _idx3d[0] = idx3d[kNestedVolumeDirAxis[0]];
2038   _idx3d[1] = idx3d[kNestedVolumeDirAxis[1]];
2039   _idx3d[2] = idx3d[kNestedVolumeDirAxis[2]];
2040   */
2041 
2042     _idx3d[kNestedVolumeDirAxis[0]] = _idx/plane;
2043     _idx3d[kNestedVolumeDirAxis[1]] = (_idx%plane)/line;
2044     _idx3d[kNestedVolumeDirAxis[2]] = (_idx%plane)%line;
2045 
2046 
2047 
2048   /*
2049 
2050   G4cout << "G4GMocrenFileSceneHandler::GetNestedVolumeIndex : " << G4endl;
2051   G4cout << "(depi, depj, depk) : "
2052    << kNestedVolumeDirAxis[0] << ", "
2053    << kNestedVolumeDirAxis[1] << ", "
2054    << kNestedVolumeDirAxis[2] << G4endl;
2055   G4cout << "(ni, nj, nk) :"
2056    << kNestedVolumeDimension[0] << ", " 
2057    << kNestedVolumeDimension[1] << ", "
2058    << kNestedVolumeDimension[2] << " - " << G4endl;
2059 
2060   G4cout << " _idx = " << _idx << "  :  plane = "
2061    << plane << " ,   line = " << line << G4endl;
2062   G4cout << "(idx,idy,idz) + " << _idx3d[0] << ", "
2063    << _idx3d[1] << ", " << _idx3d[2] << " + " << G4endl;
2064 
2065   */
2066 
2067 
2068 
2069   } else {
2070     
2071     G4int plane = kNestedVolumeDimension[0]*kNestedVolumeDimension[1];
2072     G4int line = kNestedVolumeDimension[0];
2073     _idx3d[kNestedVolumeDirAxis[2]] = _idx/plane;
2074     _idx3d[kNestedVolumeDirAxis[1]] = (_idx%plane)/line;
2075     _idx3d[kNestedVolumeDirAxis[0]] = (_idx%plane)%line;
2076 
2077   }
2078 
2079 }
2080 
2081 
2082 //-- --//
2083 G4GMocrenFileSceneHandler::Detector::Detector()
2084   : polyhedron(0) {
2085   color[0] = color[1] = color[2] = 255;
2086 }
2087 G4GMocrenFileSceneHandler::Detector::~Detector() {
2088   if(!polyhedron) delete polyhedron;
2089 }
2090 void G4GMocrenFileSceneHandler::Detector::clear() {
2091   name.clear();
2092   if(!polyhedron) delete polyhedron;
2093   color[0] = color[1] = color[2] = 255;
2094   transform3D = G4Transform3D::Identity;
2095 }
2096 
2097 //-- --//
2098 G4GMocrenFileSceneHandler::Index3D::Index3D()
2099   : x(0), y(0), z(0) {
2100   ;
2101 }
2102 
2103 G4GMocrenFileSceneHandler::Index3D::Index3D(const Index3D & _index3D) 
2104   : x(_index3D.x), y(_index3D.y), z(_index3D.z) {
2105   //: x(_index3D.X()),
2106   //y(_index3D.Y()),
2107   //z(_index3D.Z()) {
2108   //  : x(static_cast<Index3D>(_index3D).x),
2109   //    y(static_cast<Index3D>(_index3D).y),
2110   //    z(static_cast<Index3D>(_index3D).z) {
2111   ;
2112 }
2113 
2114 G4GMocrenFileSceneHandler::Index3D::Index3D(G4int _x, G4int _y, G4int _z) 
2115   : x(_x), y(_y), z(_z) {
2116   ;
2117 }
2118 G4bool G4GMocrenFileSceneHandler::Index3D::operator < (const Index3D & _right) const {
2119   if(z < static_cast<Index3D>(_right).z) {
2120      return true;
2121   } else if(z == _right.z) {
2122     if(y < static_cast<Index3D>(_right).y) return true;
2123     else if(y == _right.y) 
2124       if(x < static_cast<Index3D>(_right).x) return true;
2125   } 
2126   return false;
2127 }
2128 G4bool G4GMocrenFileSceneHandler::Index3D::operator == (const Index3D & _right) const {
2129   if(z == _right.z && y == _right.y && x == _right.x) return true;
2130   return false;
2131 }
2132