ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ScoringBox.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ScoringBox.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 
29 #include "G4ScoringBox.hh"
30 
31 #include "G4Box.hh"
32 #include "G4LogicalVolume.hh"
33 #include "G4VPhysicalVolume.hh"
34 #include "G4PVPlacement.hh"
35 #include "G4PVReplica.hh"
36 #include "G4PVDivision.hh"
37 #include "G4VisAttributes.hh"
38 #include "G4VVisManager.hh"
39 #include "G4VScoreColorMap.hh"
40 
41 #include "G4MultiFunctionalDetector.hh"
42 #include "G4SDParticleFilter.hh"
43 #include "G4VPrimitiveScorer.hh"
44 #include "G4Polyhedron.hh"
45 
46 #include "G4ScoringManager.hh"
47 #include "G4StatDouble.hh"
48 
49 #include "G4SystemOfUnits.hh"
50 
51 #include <map>
52 #include <fstream>
53 
54 G4ScoringBox::G4ScoringBox(G4String wName)
55  :G4VScoringMesh(wName), fSegmentDirection(-1)
56 {
57  fShape = MeshShape::box;
58  fDivisionAxisNames[0] = "X";
59  fDivisionAxisNames[1] = "Y";
60  fDivisionAxisNames[2] = "Z";
61 }
62 
63 G4ScoringBox::~G4ScoringBox()
64 {
65 }
66 
67 void G4ScoringBox::SetupGeometry(G4VPhysicalVolume * fWorldPhys) {
68 
69  if(verboseLevel > 9) G4cout << "G4ScoringBox::SetupGeometry() ..." << G4endl;
70 
71  // World
72  G4VPhysicalVolume * scoringWorld = fWorldPhys;
73  G4LogicalVolume * worldLogical = scoringWorld->GetLogicalVolume();
74 
75  // Scoring Mesh
76  if(verboseLevel > 9) G4cout << fWorldName << G4endl;
77  G4String boxName = fWorldName;
78 
79  if(verboseLevel > 9)
80  G4cout << fSize[0] << ", " << fSize[1] << ", " << fSize[2] << G4endl;
81  G4VSolid * boxSolid = new G4Box(boxName+"0", fSize[0], fSize[1], fSize[2]);
82  G4LogicalVolume * boxLogical = new G4LogicalVolume(boxSolid, 0, boxName+"_0");
83  new G4PVPlacement(fRotationMatrix, fCenterPosition,
84  boxLogical, boxName+"0", worldLogical, false, 0);
85 
86  //G4double fsegment[3][3];
87  //G4int segOrder[3];
88  //GetSegmentOrder(fSegmentDirection, fNSegment, segOrder, fsegment);
89  //EAxis axis[3] = {kXAxis, kYAxis, kZAxis};
90 
91  G4String layerName[2] = {boxName + "_1", boxName + "_2"};
92  G4VSolid * layerSolid[2];
93  G4LogicalVolume * layerLogical[2];
94 
95  //-- fisrt nested layer (replicated to x direction)
96  if(verboseLevel > 9) G4cout << "layer 1 :" << G4endl;
97  layerSolid[0] = new G4Box(layerName[0],
98  fSize[0]/fNSegment[0],
99  fSize[1],
100  fSize[2]);
101  layerLogical[0] = new G4LogicalVolume(layerSolid[0], 0, layerName[0]);
102  if(fNSegment[0] > 1) {
103  if(verboseLevel > 9)
104  G4cout << "G4ScoringBox::Construct() : Replicate to x direction" << G4endl;
105  if(G4ScoringManager::GetReplicaLevel()>0)
106  {
107  new G4PVReplica(layerName[0], layerLogical[0], boxLogical, kXAxis,
108  fNSegment[0], fSize[0]/fNSegment[0]*2.);
109  }
110  else
111  {
112  new G4PVDivision(layerName[0], layerLogical[0], boxLogical, kXAxis,
113  fNSegment[0], 0.);
114  }
115  } else if(fNSegment[0] == 1) {
116  if(verboseLevel > 9)
117  G4cout << "G4ScoringBox::Construct() : Placement" << G4endl;
118  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), layerLogical[0], layerName[0],
119  boxLogical, false, 0);
120  } else
121  G4cerr << "ERROR : G4ScoringBox::SetupGeometry() : invalid parameter ("
122  << fNSegment[0] << ") "
123  << "in placement of the first nested layer." << G4endl;
124 
125  if(verboseLevel > 9) {
126  G4cout << fSize[0]/fNSegment[0] << ", "
127  << fSize[1] << ", "
128  << fSize[2] << G4endl;
129  G4cout << layerName[0] << ": kXAxis, "
130  << fNSegment[0] << ", "
131  << 2.*fSize[0]/fNSegment[0] << G4endl;
132  }
133 
134  // second nested layer (replicated to y direction)
135  if(verboseLevel > 9) G4cout << "layer 2 :" << G4endl;
136  layerSolid[1] = new G4Box(layerName[1],
137  fSize[0]/fNSegment[0],
138  fSize[1]/fNSegment[1],
139  fSize[2]);
140  layerLogical[1] = new G4LogicalVolume(layerSolid[1], 0, layerName[1]);
141  if(fNSegment[1] > 1) {
142  if(verboseLevel > 9)
143  G4cout << "G4ScoringBox::Construct() : Replicate to y direction" << G4endl;
144  if(G4ScoringManager::GetReplicaLevel()>1)
145  {
146  new G4PVReplica(layerName[1], layerLogical[1], layerLogical[0], kYAxis,
147  fNSegment[1], fSize[1]/fNSegment[1]*2.);
148  }
149  else
150  {
151  new G4PVDivision(layerName[1], layerLogical[1], layerLogical[0], kYAxis,
152  fNSegment[1], 0.);
153  }
154  } else if(fNSegment[1] == 1) {
155  if(verboseLevel > 9)
156  G4cout << "G4ScoringBox::Construct() : Placement" << G4endl;
157  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), layerLogical[1], layerName[1],
158  layerLogical[0], false, 0);
159  } else
160  G4cerr << "ERROR : G4ScoringBox::SetupGeometry() : invalid parameter ("
161  << fNSegment[1] << ") "
162  << "in placement of the second nested layer." << G4endl;
163 
164  if(verboseLevel > 9) {
165  G4cout << fSize[0]/fNSegment[0] << ", "
166  << fSize[1]/fNSegment[1] << ", "
167  << fSize[2] << G4endl;
168  G4cout << layerName[1] << ": kYAxis, "
169  << fNSegment[1] << ", "
170  << 2.*fSize[1]/fNSegment[1] << G4endl;
171  }
172 
173  // mesh elements (replicated to z direction)
174  if(verboseLevel > 9) G4cout << "mesh elements :" << G4endl;
175  G4String elementName = boxName +"_3";
176  G4VSolid * elementSolid = new G4Box(elementName,
177  fSize[0]/fNSegment[0],
178  fSize[1]/fNSegment[1],
179  fSize[2]/fNSegment[2]);
180  fMeshElementLogical = new G4LogicalVolume(elementSolid, 0, elementName);
181  if(fNSegment[2] > 1) {
182  if(verboseLevel > 9)
183  G4cout << "G4ScoringBox::Construct() : Replicate to z direction" << G4endl;
184 
185  if(G4ScoringManager::GetReplicaLevel()>2)
186  {
187  new G4PVReplica(elementName, fMeshElementLogical, layerLogical[1], kZAxis,
188  fNSegment[2], 2.*fSize[2]/fNSegment[2]);
189  }
190  else
191  {
192  new G4PVDivision(elementName, fMeshElementLogical, layerLogical[1], kZAxis,
193  fNSegment[2], 0.);
194  }
195  } else if(fNSegment[2] == 1) {
196  if(verboseLevel > 9)
197  G4cout << "G4ScoringBox::Construct() : Placement" << G4endl;
198  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), fMeshElementLogical,
199  elementName, layerLogical[1], false, 0);
200  } else
201  G4cerr << "ERROR : G4ScoringBox::SetupGeometry() : "
202  << "invalid parameter (" << fNSegment[2] << ") "
203  << "in mesh element placement." << G4endl;
204 
205  if(verboseLevel > 9) {
206  G4cout << fSize[0]/fNSegment[0] << ", "
207  << fSize[1]/fNSegment[1] << ", "
208  << fSize[2]/fNSegment[2] << G4endl;
209  G4cout << elementName << ": kZAxis, "
210  << fNSegment[2] << ", "
211  << 2.*fSize[2]/fNSegment[2] << G4endl;
212  }
213 
214 
215  // set the sensitive detector
216  fMeshElementLogical->SetSensitiveDetector(fMFD);
217 
218 
219  // vis. attributes
220  G4VisAttributes * visatt = new G4VisAttributes(G4Colour(.5,.5,.5));
221  visatt->SetVisibility(false);
222  layerLogical[0]->SetVisAttributes(visatt);
223  layerLogical[1]->SetVisAttributes(visatt);
224  visatt->SetVisibility(true);
225  fMeshElementLogical->SetVisAttributes(visatt);
226 }
227 
228 
229 void G4ScoringBox::List() const {
230  G4cout << "G4ScoringBox : " << fWorldName << " --- Shape: Box mesh" << G4endl;
231  G4cout << " Size (x, y, z): ("
232  << fSize[0]/cm << ", "
233  << fSize[1]/cm << ", "
234  << fSize[2]/cm << ") [cm]"
235  << G4endl;
236 
237  G4VScoringMesh::List();
238 }
239 
240 void G4ScoringBox::Draw(RunScore * map,
241  G4VScoreColorMap* colorMap, G4int axflg) {
242 
243  G4VVisManager * pVisManager = G4VVisManager::GetConcreteInstance();
244  if(pVisManager) {
245 
246  // cell vectors
247  std::vector<std::vector<std::vector<double> > > cell; // cell[X][Y][Z]
248  std::vector<double> ez;
249  for(int z = 0; z < fNSegment[2]; z++) ez.push_back(0.);
250  std::vector<std::vector<double> > eyz;
251  for(int y = 0; y < fNSegment[1]; y++) eyz.push_back(ez);
252  for(int x = 0; x < fNSegment[0]; x++) cell.push_back(eyz);
253 
254  std::vector<std::vector<double> > xycell; // xycell[X][Y]
255  std::vector<double> ey;
256  for(int y = 0; y < fNSegment[1]; y++) ey.push_back(0.);
257  for(int x = 0; x < fNSegment[0]; x++) xycell.push_back(ey);
258 
259  std::vector<std::vector<double> > yzcell; // yzcell[Y][Z]
260  for(int y = 0; y < fNSegment[1]; y++) yzcell.push_back(ez);
261 
262  std::vector<std::vector<double> > xzcell; // xzcell[X][Z]
263  for(int x = 0; x < fNSegment[0]; x++) xzcell.push_back(ez);
264 
265  // projections
266  G4int q[3];
267  std::map<G4int, G4StatDouble*>::iterator itr = map->GetMap()->begin();
268  for(; itr != map->GetMap()->end(); itr++) {
269  GetXYZ(itr->first, q);
270 
271  xycell[q[0]][q[1]] += (itr->second->sum_wx())/fDrawUnitValue;
272  yzcell[q[1]][q[2]] += (itr->second->sum_wx())/fDrawUnitValue;
273  xzcell[q[0]][q[2]] += (itr->second->sum_wx())/fDrawUnitValue;
274  }
275 
276  // search max. & min. values in each slice
277  G4double xymin = DBL_MAX, yzmin = DBL_MAX, xzmin = DBL_MAX;
278  G4double xymax = 0., yzmax = 0., xzmax = 0.;
279  for(int x = 0; x < fNSegment[0]; x++) {
280  for(int y = 0; y < fNSegment[1]; y++) {
281  if(xymin > xycell[x][y]) xymin = xycell[x][y];
282  if(xymax < xycell[x][y]) xymax = xycell[x][y];
283  }
284  for(int z = 0; z < fNSegment[2]; z++) {
285  if(xzmin > xzcell[x][z]) xzmin = xzcell[x][z];
286  if(xzmax < xzcell[x][z]) xzmax = xzcell[x][z];
287  }
288  }
289  for(int y = 0; y < fNSegment[1]; y++) {
290  for(int z = 0; z < fNSegment[2]; z++) {
291  if(yzmin > yzcell[y][z]) yzmin = yzcell[y][z];
292  if(yzmax < yzcell[y][z]) yzmax = yzcell[y][z];
293  }
294  }
295 
296 
297  G4VisAttributes att;
298  att.SetForceSolid(true);
299  att.SetForceAuxEdgeVisible(true);
300  G4double thick = 0.01;
301 
302  G4Scale3D scale;
303  if(axflg/100==1) {
304  pVisManager->BeginDraw();
305 
306  // xy plane
307  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(xymin ,xymax); }
308  G4ThreeVector zhalf(0., 0., fSize[2]/fNSegment[2]-thick);
309  for(int x = 0; x < fNSegment[0]; x++) {
310  for(int y = 0; y < fNSegment[1]; y++) {
311 
312  G4ThreeVector pos(GetReplicaPosition(x, y, 0) - zhalf);
313  G4ThreeVector pos2(GetReplicaPosition(x, y, fNSegment[2]-1) + zhalf);
314  G4Transform3D trans, trans2;
315  if(fRotationMatrix) {
316  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos);
317  trans = G4Translate3D(fCenterPosition)*trans;
318  trans2 = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos2);
319  trans2 = G4Translate3D(fCenterPosition)*trans2;
320  } else {
321  trans = G4Translate3D(pos)*G4Translate3D(fCenterPosition);
322  trans2 = G4Translate3D(pos2)*G4Translate3D(fCenterPosition);
323  }
324  G4double c[4];
325  colorMap->GetMapColor(xycell[x][y], c);
326  att.SetColour(c[0], c[1], c[2]);//, c[3]);
327 
328  G4Box xyplate("xy", fSize[0]/fNSegment[0], fSize[1]/fNSegment[1],
329  thick);
330  G4Polyhedron * poly = xyplate.GetPolyhedron();
331  poly->Transform(trans);
332  poly->SetVisAttributes(&att);
333  pVisManager->Draw(*poly);
334 
335  G4Box xyplate2 = xyplate;
336  G4Polyhedron * poly2 = xyplate2.GetPolyhedron();
337  poly2->Transform(trans2);
338  poly2->SetVisAttributes(&att);
339  pVisManager->Draw(*poly2);
340 
341  /*
342  G4double nodes[][3] =
343  {{-fSize[0]/fNSegment[0], -fSize[1]/fNSegment[1], 0.},
344  { fSize[0]/fNSegment[0], -fSize[1]/fNSegment[1], 0.},
345  { fSize[0]/fNSegment[0], fSize[1]/fNSegment[1], 0.},
346  {-fSize[0]/fNSegment[0], fSize[1]/fNSegment[1], 0.}};
347  G4int facets[][4] = {{4, 3, 2, 1}};
348  G4int facets2[][4] = {{1, 2, 3, 4}};
349 
350  G4Polyhedron poly, poly2;
351  poly.createPolyhedron(4, 1, nodes, facets);
352  poly.Transform(trans);
353  poly.SetVisAttributes(att);
354  pVisManager->Draw(poly);
355 
356  poly2.createPolyhedron(4, 1, nodes, facets2);
357  poly2.Transform(trans2);
358  poly2.SetVisAttributes(att);
359  pVisManager->Draw(poly2);
360  */
361  }
362  }
363  pVisManager->EndDraw();
364  }
365  axflg = axflg%100;
366  if(axflg/10==1) {
367  pVisManager->BeginDraw();
368 
369  // yz plane
370  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(yzmin, yzmax); }
371  G4ThreeVector xhalf(fSize[0]/fNSegment[0]-thick, 0., 0.);
372  for(int y = 0; y < fNSegment[1]; y++) {
373  for(int z = 0; z < fNSegment[2]; z++) {
374 
375  G4ThreeVector pos(GetReplicaPosition(0, y, z) - xhalf);
376  G4ThreeVector pos2(GetReplicaPosition(fNSegment[0]-1, y, z) + xhalf);
377  G4Transform3D trans, trans2;
378  if(fRotationMatrix) {
379  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos);
380  trans = G4Translate3D(fCenterPosition)*trans;
381  trans2 = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos2);
382  trans2 = G4Translate3D(fCenterPosition)*trans2;
383  } else {
384  trans = G4Translate3D(pos)*G4Translate3D(fCenterPosition);
385  trans2 = G4Translate3D(pos2)*G4Translate3D(fCenterPosition);
386  }
387  G4double c[4];
388  colorMap->GetMapColor(yzcell[y][z], c);
389  att.SetColour(c[0], c[1], c[2]);//, c[3]);
390 
391  G4Box yzplate("yz", thick,//fSize[0]/fNSegment[0]*0.001,
392  fSize[1]/fNSegment[1],
393  fSize[2]/fNSegment[2]);
394  G4Polyhedron * poly = yzplate.GetPolyhedron();
395  poly->Transform(trans);
396  poly->SetVisAttributes(&att);
397  pVisManager->Draw(*poly);
398 
399  G4Box yzplate2 = yzplate;
400  G4Polyhedron * poly2 = yzplate2.GetPolyhedron();
401  poly2->Transform(trans2);
402  poly2->SetVisAttributes(&att);
403  pVisManager->Draw(*poly2);
404 
405  /*
406  G4double nodes[][3] =
407  {{0., -fSize[1]/fNSegment[1], -fSize[2]/fNSegment[2]},
408  {0., fSize[1]/fNSegment[1], -fSize[2]/fNSegment[2]},
409  {0., fSize[1]/fNSegment[1], fSize[2]/fNSegment[2]},
410  {0., -fSize[1]/fNSegment[1], fSize[2]/fNSegment[2]}};
411  G4int facets[][4] = {{4, 3, 2, 1}};
412  G4int facets2[][4] = {{1, 2, 3, 4}};
413 
414  G4Polyhedron poly, poly2;
415  poly.createPolyhedron(4, 1, nodes, facets);
416  poly.Transform(trans);
417  poly.SetVisAttributes(att);
418  pVisManager->Draw(poly);
419 
420  poly2.createPolyhedron(4, 1, nodes, facets2);
421  poly2.Transform(trans2);
422  poly2.SetVisAttributes(att);
423  pVisManager->Draw(poly2);
424  */
425  }
426  }
427  pVisManager->EndDraw();
428  }
429  axflg = axflg%10;
430  if(axflg==1) {
431  pVisManager->BeginDraw();
432 
433  // xz plane
434  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(xzmin,xzmax); }
435  G4ThreeVector yhalf(0., fSize[1]/fNSegment[1]-thick, 0.);
436  for(int x = 0; x < fNSegment[0]; x++) {
437  for(int z = 0; z < fNSegment[2]; z++) {
438 
439  G4ThreeVector pos(GetReplicaPosition(x, 0, z) - yhalf);
440  G4ThreeVector pos2(GetReplicaPosition(x, fNSegment[1]-1, z) + yhalf);
441  G4Transform3D trans, trans2;
442  if(fRotationMatrix) {
443  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos);
444  trans = G4Translate3D(fCenterPosition)*trans;
445  trans2 = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos2);
446  trans2 = G4Translate3D(fCenterPosition)*trans2;
447  } else {
448  trans = G4Translate3D(pos)*G4Translate3D(fCenterPosition);
449  trans2 = G4Translate3D(pos2)*G4Translate3D(fCenterPosition);
450  }
451  G4double c[4];
452  colorMap->GetMapColor(xzcell[x][z], c);
453  att.SetColour(c[0], c[1], c[2]);//, c[3]);
454 
455  G4Box xzplate("xz", fSize[0]/fNSegment[0], thick,//fSize[1]/fNSegment[1]*0.001,
456  fSize[2]/fNSegment[2]);
457  G4Polyhedron * poly = xzplate.GetPolyhedron();
458  poly->Transform(trans);
459  poly->SetVisAttributes(&att);
460  pVisManager->Draw(*poly);
461 
462  G4Box xzplate2 = xzplate;
463  G4Polyhedron * poly2 = xzplate2.GetPolyhedron();
464  poly2->Transform(trans2);
465  poly2->SetVisAttributes(&att);
466  pVisManager->Draw(*poly2);
467 
468 
469  /*
470  G4double nodes[][3] =
471  {{-fSize[1]/fNSegment[1], 0., -fSize[2]/fNSegment[2]},
472  { fSize[1]/fNSegment[1], 0., -fSize[2]/fNSegment[2]},
473  { fSize[1]/fNSegment[1], 0., fSize[2]/fNSegment[2]},
474  {-fSize[1]/fNSegment[1], 0., fSize[2]/fNSegment[2]}};
475  G4int facets[][4] = {{1, 2, 3, 4}};
476  G4int facets2[][4] = {{4, 3, 2, 1}};
477 
478  G4Polyhedron poly, poly2;
479  poly.createPolyhedron(4, 1, nodes, facets);
480  poly.Transform(trans);
481  poly.SetVisAttributes(att);
482  pVisManager->Draw(poly);
483 
484  poly2.createPolyhedron(4, 1, nodes, facets2);
485  poly2.Transform(trans2);
486  poly2.SetVisAttributes(att);
487  pVisManager->Draw(poly2);
488  */
489  }
490  }
491  pVisManager->EndDraw();
492  }
493  }
494  colorMap->SetPSUnit(fDrawUnit);
495  colorMap->SetPSName(fDrawPSName);
496  colorMap->DrawColorChart();
497 }
498 
499 G4ThreeVector G4ScoringBox::GetReplicaPosition(G4int x, G4int y, G4int z) {
500 
501  G4ThreeVector width(fSize[0]/fNSegment[0], fSize[1]/fNSegment[1],
502  fSize[2]/fNSegment[2]);
503  G4ThreeVector pos(-fSize[0] + 2*(x+0.5)*width.x(),
504  -fSize[1] + 2*(y+0.5)*width.y(),
505  -fSize[2] + 2*(z+0.5)*width.z());
506 
507  return pos;
508 }
509 
510 void G4ScoringBox::GetXYZ(G4int index, G4int q[3]) const {
511 
512  q[0] = index/(fNSegment[2]*fNSegment[1]);
513  q[1] = (index - q[0]*fNSegment[2]*fNSegment[1])/fNSegment[2];
514  q[2] = index - q[1]*fNSegment[2] - q[0]*fNSegment[2]*fNSegment[1];
515 
516 }
517 
518 G4int G4ScoringBox::GetIndex(G4int x, G4int y, G4int z) const {
519  return x + y*fNSegment[0] + z*fNSegment[0]*fNSegment[1];
520 }
521 
522 void G4ScoringBox::DrawColumn(RunScore * map,
523  G4VScoreColorMap* colorMap,
524  G4int idxProj, G4int idxColumn)
525 {
526  G4int iColumn[3] = {2, 0, 1};
527  if(idxColumn<0 || idxColumn>=fNSegment[iColumn[idxProj]])
528  {
529  G4cerr << "ERROR : Column number " << idxColumn
530  << " is out of scoring mesh [0," << fNSegment[iColumn[idxProj]]-1
531  << "]. Method ignored." << G4endl;
532  return;
533  }
534  G4VVisManager * pVisManager = G4VVisManager::GetConcreteInstance();
535  if(pVisManager) {
536  pVisManager->BeginDraw();
537 
538  // cell vectors
539  std::vector<std::vector<std::vector<double> > > cell; // cell[X][Y][Z]
540  std::vector<double> ez;
541  for(int z = 0; z < fNSegment[2]; z++) ez.push_back(0.);
542  std::vector<std::vector<double> > eyz;
543  for(int y = 0; y < fNSegment[1]; y++) eyz.push_back(ez);
544  for(int x = 0; x < fNSegment[0]; x++) cell.push_back(eyz);
545 
546  std::vector<std::vector<double> > xycell; // xycell[X][Y]
547  std::vector<double> ey;
548  for(int y = 0; y < fNSegment[1]; y++) ey.push_back(0.);
549  for(int x = 0; x < fNSegment[0]; x++) xycell.push_back(ey);
550 
551  std::vector<std::vector<double> > yzcell; // yzcell[Y][Z]
552  for(int y = 0; y < fNSegment[1]; y++) yzcell.push_back(ez);
553 
554  std::vector<std::vector<double> > xzcell; // xzcell[X][Z]
555  for(int x = 0; x < fNSegment[0]; x++) xzcell.push_back(ez);
556 
557  // projections
558  G4int q[3];
559  std::map<G4int, G4StatDouble*>::iterator itr = map->GetMap()->begin();
560  for(; itr != map->GetMap()->end(); itr++) {
561  GetXYZ(itr->first, q);
562 
563  if(idxProj == 0 && q[2] == idxColumn) { // xy plane
564  xycell[q[0]][q[1]] += (itr->second->sum_wx())/fDrawUnitValue;
565  }
566  if(idxProj == 1 && q[0] == idxColumn) { // yz plane
567  yzcell[q[1]][q[2]] += (itr->second->sum_wx())/fDrawUnitValue;
568  }
569  if(idxProj == 2 && q[1] == idxColumn) { // zx plane
570  xzcell[q[0]][q[2]] += (itr->second->sum_wx())/fDrawUnitValue;
571  }
572  }
573 
574  // search max. & min. values in each slice
575  G4double xymin = DBL_MAX, yzmin = DBL_MAX, xzmin = DBL_MAX;
576  G4double xymax = 0., yzmax = 0., xzmax = 0.;
577  for(int x = 0; x < fNSegment[0]; x++) {
578  for(int y = 0; y < fNSegment[1]; y++) {
579  if(xymin > xycell[x][y]) xymin = xycell[x][y];
580  if(xymax < xycell[x][y]) xymax = xycell[x][y];
581  }
582  for(int z = 0; z < fNSegment[2]; z++) {
583  if(xzmin > xzcell[x][z]) xzmin = xzcell[x][z];
584  if(xzmax < xzcell[x][z]) xzmax = xzcell[x][z];
585  }
586  }
587  for(int y = 0; y < fNSegment[1]; y++) {
588  for(int z = 0; z < fNSegment[2]; z++) {
589  if(yzmin > yzcell[y][z]) yzmin = yzcell[y][z];
590  if(yzmax < yzcell[y][z]) yzmax = yzcell[y][z];
591  }
592  }
593 
594 
595  G4VisAttributes att;
596  att.SetForceSolid(true);
597  att.SetForceAuxEdgeVisible(true);
598 
599  G4Scale3D scale;
600  // xy plane
601  if(idxProj == 0) {
602  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(xymin,xymax); }
603  for(int x = 0; x < fNSegment[0]; x++) {
604  for(int y = 0; y < fNSegment[1]; y++) {
605  G4Box xyplate("xy", fSize[0]/fNSegment[0], fSize[1]/fNSegment[1],
606  fSize[2]/fNSegment[2]);
607 
608  G4ThreeVector pos(GetReplicaPosition(x, y, idxColumn));
609  G4Transform3D trans;
610  if(fRotationMatrix) {
611  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos);
612  trans = G4Translate3D(fCenterPosition)*trans;
613  } else {
614  trans = G4Translate3D(pos)*G4Translate3D(fCenterPosition);
615  }
616  G4double c[4];
617  colorMap->GetMapColor(xycell[x][y], c);
618  att.SetColour(c[0], c[1], c[2]);
619 
620  G4Polyhedron * poly = xyplate.GetPolyhedron();
621  poly->Transform(trans);
622  poly->SetVisAttributes(att);
623  pVisManager->Draw(*poly);
624  }
625  }
626 
627  } else
628  // yz plane
629  if(idxProj == 1) {
630  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(yzmin,yzmax); }
631  for(int y = 0; y < fNSegment[1]; y++) {
632  for(int z = 0; z < fNSegment[2]; z++) {
633  G4Box yzplate("yz", fSize[0]/fNSegment[0], fSize[1]/fNSegment[1],
634  fSize[2]/fNSegment[2]);
635 
636  G4ThreeVector pos(GetReplicaPosition(idxColumn, y, z));
637  G4Transform3D trans;
638  if(fRotationMatrix) {
639  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos);
640  trans = G4Translate3D(fCenterPosition)*trans;
641  } else {
642  trans = G4Translate3D(pos)*G4Translate3D(fCenterPosition);
643  }
644  G4double c[4];
645  colorMap->GetMapColor(yzcell[y][z], c);
646  att.SetColour(c[0], c[1], c[2]);//, c[3]);
647 
648  G4Polyhedron * poly = yzplate.GetPolyhedron();
649  poly->Transform(trans);
650  poly->SetVisAttributes(att);
651  pVisManager->Draw(*poly);
652  }
653  }
654  } else
655  // xz plane
656  if(idxProj == 2) {
657  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(xzmin,xzmax);}
658  for(int x = 0; x < fNSegment[0]; x++) {
659  for(int z = 0; z < fNSegment[2]; z++) {
660  G4Box xzplate("xz", fSize[0]/fNSegment[0], fSize[1]/fNSegment[1],
661  fSize[2]/fNSegment[2]);
662 
663  G4ThreeVector pos(GetReplicaPosition(x, idxColumn, z));
664  G4Transform3D trans;
665  if(fRotationMatrix) {
666  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(pos);
667  trans = G4Translate3D(fCenterPosition)*trans;
668  } else {
669  trans = G4Translate3D(pos)*G4Translate3D(fCenterPosition);
670  }
671  G4double c[4];
672  colorMap->GetMapColor(xzcell[x][z], c);
673  att.SetColour(c[0], c[1], c[2]);//, c[3]);
674 
675  G4Polyhedron * poly = xzplate.GetPolyhedron();
676  poly->Transform(trans);
677  poly->SetVisAttributes(att);
678  pVisManager->Draw(*poly);
679  }
680  }
681  }
682  pVisManager->EndDraw();
683 
684  }
685 
686  colorMap->SetPSUnit(fDrawUnit);
687  colorMap->SetPSName(fDrawPSName);
688  colorMap->DrawColorChart();
689 }
690 
691