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