source: trunk/source/visualization/gMocren/src/G4GMocrenFileSceneHandler.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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