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

Last change on this file since 1288 was 1288, checked in by garnier, 14 years ago

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