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