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