Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4GMocrenFileSceneHandler.cc 110513 2018-05-28 07:37:38Z gcosmo $
28 //
29 //
30 // Created: Mar. 31, 2009 Akinori Kimura
31 // Sep. 22, 2009 Akinori Kimura : modify and fix code to support
32 // PrimitiveScorers and clean up
33 //
34 // GMocrenFile scene handler
35 
36 
37 //----- header files
38 #include <fstream>
39 #include <cstdlib>
40 #include <cstring>
41 #include <sstream>
42 #include <sstream>
43 #include <iomanip>
44 
45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4VisManager.hh"
49 
50 #include "G4GMocrenFile.hh"
52 #include "G4GMocrenFileViewer.hh"
53 #include "G4Point3D.hh"
54 #include "G4VisAttributes.hh"
55 #include "G4Scene.hh"
56 #include "G4Transform3D.hh"
57 #include "G4Polyhedron.hh"
58 #include "G4Box.hh"
59 #include "G4Cons.hh"
60 #include "G4Polyline.hh"
61 #include "G4Trd.hh"
62 #include "G4Tubs.hh"
63 #include "G4Trap.hh"
64 #include "G4Torus.hh"
65 #include "G4Sphere.hh"
66 #include "G4Para.hh"
67 #include "G4Text.hh"
68 #include "G4Circle.hh"
69 #include "G4Square.hh"
70 #include "G4VPhysicalVolume.hh"
71 #include "G4PhysicalVolumeModel.hh"
72 #include "G4LogicalVolume.hh"
73 #include "G4Material.hh"
74 
75 #include "G4VPVParameterisation.hh"
77 #include "G4VisTrajContext.hh"
78 #include "G4TrajectoriesModel.hh"
79 #include "G4VTrajectoryModel.hh"
81 #include "G4HitsModel.hh"
82 #include "G4GMocrenMessenger.hh"
83 //#include "G4PSHitsModel.hh"
84 #include "G4GMocrenIO.hh"
86 #include "G4GMocrenTouchable.hh"
90 
91 #include "G4ScoringManager.hh"
92 #include "G4ScoringBox.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
114 
116 // Driver-dependent part //
118 
119 
120 //----- G4GMocrenFileSceneHandler, constructor
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 = 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 {
167  }
168  if( kMaxFileNum < 1 ) { kMaxFileNum = 1 ; }
169 
171 
172 }
173 
174 
175 //----- G4GMocrenFileSceneHandler, destructor
177 {
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
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 //-----
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  {
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
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 //-----
282 {
284  G4cout << "***** BeginSavingGdd (called)" << G4endl;
285 
286  if( !IsSavingGdd() ) {
287 
289  G4cout << "***** (started) " ;
290  G4cout << "(open g4.gdd, ##)" << G4endl;
291  }
292 
293  SetGddFileName() ; // result set to kGddFileName
294  kFlagSaving_g4_gdd = true;
295 
296 
298  short minmax[2];
299  minmax[0] = ctdens.GetMinCT();
300  minmax[1] = ctdens.GetMaxCT();
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  }
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();
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 
353 {
355  G4cout << "***** EndSavingGdd (called)" << G4endl;
356 
357  if(IsSavingGdd()) {
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;
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 
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();
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
504  }
505 
506 }
507 
508 
509 //-----
511 {
513 
514  if( !GFIsInModeling() ) {
515 
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
544 {
546  G4cout << "***** AddPrimitive" << G4endl;
547 
548  if (fProcessing2D) {
549  static G4bool warned = false;
550  if (!warned) {
551  warned = true;
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
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();
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
636 {
637  if (fProcessing2D) {
638  static G4bool warned = false;
639  if (!warned) {
640  warned = true;
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  //-----
654  G4cout << "***** AddPrimitive( G4Text )" << G4endl;
655 
656  //----- Initialize IF NECESSARY
657  GFBeginModeling();
658 
659 } // G4GMocrenFileSceneHandler::AddPrimitive ( text )
660 
661 
662 //----- Add 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;
673  ("G4GMocrenFileSceneHandler::AddPrimitive (const G4Circle&)",
674  "gMocren1003", JustWarning,
675  "2D circles not implemented. Ignored.");
676  }
677  return;
678  }
679 
680  //-----
682  G4cout << "***** AddPrimitive( G4Circle )" << G4endl;
683 
684  //----- Initialize IF NECESSARY
685  GFBeginModeling();
686 
687 
688 } // G4GMocrenFileSceneHandler::AddPrimitive ( mark_circle )
689 
690 
691 //----- Add 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;
702  ("G4GMocrenFileSceneHandler::AddPrimitive (const G4Square&)",
703  "gMocren1004", JustWarning,
704  "2D squares not implemented. Ignored.");
705  }
706  return;
707  }
708 
709  //-----
711  G4cout << "***** AddPrimitive( G4Square )" << G4endl;
712 
713  //----- Initialize if necessary
714  GFBeginModeling();
715 
716 } // G4GMocrenFileSceneHandler::AddPrimitive ( mark_square )
717 
718 
719 //----- Add polyhedron
721 {
722  //-----
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;
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:
767  G4cout <<
768  "ERROR G4GMocrenFileSceneHandler::AddPrimitive(G4Polyhedron)" << G4endl;
769  G4PhysicalVolumeModel* pPVModel =
770  dynamic_cast<G4PhysicalVolumeModel*>(fpModel);
771  if (pPVModel)
773  G4cout << "Volume " << pPVModel->GetCurrentPV()->GetName() <<
774  ", Solid " << pPVModel->GetCurrentLV()->GetSolid()->GetName() <<
775  " (" << pPVModel->GetCurrentLV()->GetSolid()->GetEntityType();
776 
778  G4cout <<
779  "\nG4Polyhedron facet with " << i << " edges" << G4endl;
780  }
781  }
782 
783 } // G4GMocrenFileSceneHandler::AddPrimitive (polyhedron)
784 
785 
786 //-----
788 {
790 
791  //-----
793  G4cout << "***** GFEndModeling (called)" << G4endl;
794 
795  if( GFIsInModeling() ) {
796 
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 //-----
815 {
817  G4cout << "***** BeginPrimitives " << G4endl;
818 
819  GFBeginModeling();
820 
821  G4VSceneHandler::BeginPrimitives (objectTransformation);
822 
823 }
824 
825 
826 //-----
828 {
830  G4cout << "***** EndPrimitives " << G4endl;
831 
833 }
834 
835 
836 //========== AddSolid() functions ==========//
837 
838 //----- Add box
840 {
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();
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) {
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 
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) {
931  G4cout << "kVolumeTrans3D: " << trans1 << G4endl << rot1 << G4endl;
932  }
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) {
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];
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
1252 
1253 
1254  // get the dimension of the geometry defined in G4VScoringMesh
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 
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  //
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
1423 {
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
1466 {
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
1484 {
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
1503 {
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
1521 {
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
1539 {
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
1558 {
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
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
1598 
1599  kbModelingTrajectory = true;
1600 
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
1638  if(GFDEBUG_HIT) G4cout << " ::AddCompound(const G4VHit&) >>>>>>>>> " << G4endl;
1639 
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 
1738  if(GFDEBUG_DIGI) G4cout << " ::AddCompound(const G4VDigi&) >>>>>>>>> " << G4endl;
1740 }
1741 
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 
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 //-----
1888 {
1889  //-----
1890  G4bool visibility = true ;
1891 
1892  const G4VisAttributes* pVisAttribs =
1894 
1895  if(pVisAttribs) {
1896  visibility = pVisAttribs->IsVisible();
1897  }
1898 
1899  return visibility ;
1900 
1901 } // G4GMocrenFileSceneHandler::IsVisible()
1902 
1903 
1904 //-----
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 //-----
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;
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 //-----
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 
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 
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 
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 //-- --//
2083  : polyhedron(0) {
2084  color[0] = color[1] = color[2] = 255;
2085 }
2087  if(!polyhedron) delete polyhedron;
2088 }
2090  name.clear();
2091  if(!polyhedron) delete polyhedron;
2092  color[0] = color[1] = color[2] = 255;
2093  transform3D = G4Transform3D::Identity;
2094 }
2095 
2096 //-- --//
2098  : x(0), y(0), z(0) {
2099  ;
2100 }
2101 
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 
2114  : x(_x), y(_y), z(_z) {
2115  ;
2116 }
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 }
2128  if(z == _right.z && y == _right.y && x == _right.x) return true;
2129  return false;
2130 }
G4ThreeVector GetSize() const
Float_t x
Definition: compare.C:6
const XML_Char int len
Definition: expat.h:262
const G4String & GetName() const
Hep3Vector & transform(const HepRotation &)
Definition: ThreeVectorR.cc:24
const G4bool GFDEBUG
G4Point3D GetVertex(G4int index) const
void set(double x, double y, double z)
CLHEP::HepRotation getRotation() const
const XML_Char * name
Definition: expat.h:151
G4double GetYHalfLength() const
const G4bool GFDEBUG_TRK
virtual const G4ThreeVector GetPosition() const =0
bool isIdentity() const
Definition: Rotation.cc:172
CLHEP::Hep3Vector getTranslation() const
const G4Color & GetColor() const
void clearDoseDistAll()
G4Polyhedron * CreatePolyhedron() const
Definition: G4Box.cc:557
size_t GetNoVoxelZ() const
G4Transform3D fObjectTransformation
system("rm -rf microbeam.root")
G4LogicalVolume * GetLogicalVolume() const
virtual G4bool IsParameterised() const =0
virtual void BeginPrimitives(const G4Transform3D &objectTransformation)
virtual G4Material * ComputeMaterial(const G4int repNo, G4VPhysicalVolume *currentVol, const G4VTouchable *parentTouch=0)
void setModalityImageMinMax(short _minmax[2])
void setDoseDistMinMax(short _minmax[2], int _num=0)
static constexpr double mm
Definition: G4SIunits.hh:115
Definition: G4Tubs.hh:85
const G4VisAttributes * GetApplicableVisAttributes(const G4VisAttributes *) const
short convertDensityToHU(float &_dens)
virtual G4VPVParameterisation * GetParameterisation() const =0
const G4VisAttributes * GetVisAttributes() const
#define G4endl
Definition: G4ios.hh:61
Transform3D inverse() const
Definition: Transform3D.cc:142
G4VViewer * fpViewer
Float_t y
Definition: compare.C:6
Definition: G4VHit.hh:48
virtual G4ThreeVector GetInitialMomentum() const =0
float G4float
Definition: G4Types.hh:77
Double_t z
G4int GetNoDaughters() const
G4double GetVoxelHalfX() const
const G4VTrajectory * GetCurrentTrajectory() const
size_t GetNoVoxelX() const
virtual std::vector< G4AttValue > * CreateAttValues() const
Definition: G4VHit.hh:67
double z() const
void setModalityImageSize(int _size[3])
G4bool IsVisible() const
G4double GetVoxelHalfZ() const
void setModalityImage(short *_image)
void setDoseDistScale(double &_scale, int _num=0)
virtual std::vector< G4String > getHitScorerNames()
void clearDetector()
Definition: G4GMocrenIO.hh:451
void setDoseDist(double *_image, int _num=0)
void clearTracks()
Definition: G4GMocrenIO.hh:439
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
virtual G4int GetMultiplicity() const
const G4int MAX_NUM_TRAJECTORIES
const G4int GFDEBUG_DET
G4bool GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
fin
Definition: AddMC.C:2
const G4String & GetName() const
Definition: G4Material.hh:179
virtual G4String getVolumeName()
bool storeData(char *_filename)
Definition: G4GMocrenIO.cc:458
static constexpr double g
Definition: G4SIunits.hh:183
void AddDetector(const G4VSolid &solid)
virtual void AddCompound(const G4VTrajectory &)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
Double_t scale
const char DEFAULT_GDD_FILE_NAME[]
virtual G4double GetCharge() const =0
#define width
G4bool GetNextEdge(G4Point3D &p1, G4Point3D &p2, G4int &edgeFlag) const
void hits()
Definition: readHits.C:15
void translateTracks(std::vector< float > &_translateo)
G4VPhysicalVolume * GetTopPhysicalVolume() const
Definition: G4Text.hh:73
const std::vector< Model > & GetEndOfEventModelList() const
const G4bool GFDEBUG_HIT
const XML_Char int const XML_Char * value
Definition: expat.h:331
Definition: G4Cons.hh:83
void GetNumberOfSegments(G4int nSegment[3])
virtual G4bool IsNested() const
virtual int GetPointEntries() const =0
virtual void BeginModeling()
HepPolyhedron & Transform(const G4Transform3D &t)
virtual G4String GetParticleName() const =0
Definition: G4Box.hh:64
const G4bool GFDEBUG_DIGI
static Verbosity GetVerbosity()
G4bool operator==(const Index3D &_right) const
virtual G4int GetCopyNo() const =0
static constexpr double rad
Definition: G4SIunits.hh:149
std::vector< Detector > kDetectors
void translateDetector(std::vector< float > &_translate)
Float_t mat
const G4VisAttributes * GetVisAttributes() const
G4double GetRed() const
Definition: G4Colour.hh:151
void AddCompound(const G4VTrajectory &traj)
G4double GetXHalfLength() const
virtual void EndModeling()
std::vector< G4String > kNestedVolumeNames
G4GMocrenFileSceneHandler(G4GMocrenFile &system, G4GMocrenMessenger &messenger, const G4String &name="")
virtual std::vector< G4String > getHitNames()
G4LogicalVolume * GetCurrentLV() const
static DLL_API const Transform3D Identity
Definition: Transform3D.h:197
const G4VisAttributes * fpVisAttribs
virtual G4Material * ComputeMaterial(G4VPhysicalVolume *currentVol, const G4int repNo, const G4VTouchable *parentTouch=0)=0
G4Material * GetCurrentMaterial() const
virtual void ComputeTransformation(const G4int no, G4VPhysicalVolume *currentPV) const =0
virtual void AddSolid(const G4Box &)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
void setVoxelSpacing(float _spacing[3])
G4int GetNoVertices() const
G4ThreeVector GetTranslation() const
G4bool operator<(const Index3D &_right) const
void addDetector(std::string &_name, std::vector< float * > &_det, unsigned char _color[3])
int G4int
Definition: G4Types.hh:78
G4VPhysicalVolume * GetCurrentPV() const
G4Polyhedron * CreatePolyhedron() const
Definition: G4Tubs.cc:1768
EAxis
Definition: geomdefs.hh:54
G4String GetName() const
virtual void DrawTrajectory() const
virtual void BeginPrimitives(const G4Transform3D &objectTransformation)
Definition: G4Trd.hh:72
void newDoseDist()
void setDoseDistName(std::string _name, int _num=0)
G4VModel * GetModel() const
const std::vector< G4PhysicalVolumeNodeID > & GetDrawnPVPath() const
virtual void GetReplicationData(EAxis &axis, G4int &nReplicas, G4double &width, G4double &offset, G4bool &consuming) const =0
virtual G4int GetTrackID() const =0
const int FR_MAX_FILE_NUM
G4RotationMatrix GetRotationMatrix() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
G4double GetVoxelHalfY() const
double x() const
virtual G4bool IsReplicated() const =0
G4VSolid * GetSolid() const
void clearROIAll()
std::map< G4String, std::map< Index3D, G4double > > kNestedHitsList
void setDoseDistUnit(std::string &_unit, int _num=0)
void setDoseDistSize(int _size[3], int _num=0)
virtual G4Polyhedron * CreatePolyhedron() const
Definition: G4VSolid.cc:660
G4double GetBlue() const
Definition: G4Colour.hh:153
void AddPrimitive(const G4Polyline &line)
Char_t n[5]
G4double GetGreen() const
Definition: G4Colour.hh:152
void setModalityImageDensityMap(std::vector< float > &_map)
static constexpr double pi
Definition: G4SIunits.hh:75
G4double GetDensity(G4int &_ct) const
double y() const
G4int GetNoFacets() const
static G4ScoringManager * GetScoringManager()
G4VScoringMesh * FindMesh(G4VHitsCollection *map)
virtual G4GeometryType GetEntityType() const =0
G4ThreeVector GetObjectTranslation() const
virtual G4bool getDrawVolumeGrid()
#define DBL_MAX
Definition: templates.hh:83
Definition: G4Para.hh:86
static constexpr double cm3
Definition: G4SIunits.hh:121
void addTrack(float *_tracks)
G4double GetZHalfLength() const
virtual G4VTrajectoryPoint * GetPoint(G4int i) const =0
HepRotation inverse() const
const char GDD_FILE_HEADER[]
virtual void EndPrimitives()
G4Scene * GetScene() const
Map_t * GetMap() const
Definition: G4THitsMap.hh:117
const G4String & GetName() const
G4VPhysicalVolume * GetDaughter(const G4int i) const
void GetNestedVolumeIndex(G4int, G4int[3])
std::map< Index3D, float > kNestedModality
G4double GetDensity() const
Definition: G4Material.hh:181
size_t GetNoVoxelY() const