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