Geant4-11
G4ScoringBox.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//
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
54G4ScoringBox::G4ScoringBox(G4String wName)
55 : G4VScoringMesh(wName)
56 , fSegmentDirection(-1)
57{
58 fShape = MeshShape::box;
59 fDivisionAxisNames[0] = "X";
60 fDivisionAxisNames[1] = "Y";
61 fDivisionAxisNames[2] = "Z";
62}
63
64G4ScoringBox::~G4ScoringBox() {}
65
66void G4ScoringBox::SetupGeometry(G4VPhysicalVolume* fWorldPhys)
67{
68 if(verboseLevel > 9)
69 G4cout << "G4ScoringBox::SetupGeometry() ..." << G4endl;
70
71 // World
72 G4VPhysicalVolume* scoringWorld = fWorldPhys;
73 G4LogicalVolume* worldLogical = scoringWorld->GetLogicalVolume();
74
75 // Scoring Mesh
76 if(verboseLevel > 9)
77 G4cout << fWorldName << G4endl;
78 G4String boxName = fWorldName;
79
80 if(verboseLevel > 9)
81 G4cout << fSize[0] << ", " << fSize[1] << ", " << fSize[2] << G4endl;
82 G4VSolid* boxSolid = new G4Box(boxName + "0", fSize[0], fSize[1], fSize[2]);
83 G4LogicalVolume* boxLogical =
84 new G4LogicalVolume(boxSolid, 0, boxName + "_0");
85 new G4PVPlacement(fRotationMatrix, fCenterPosition, boxLogical, boxName + "0",
86 worldLogical, false, 0);
87
88 // G4double fsegment[3][3];
89 // G4int segOrder[3];
90 // GetSegmentOrder(fSegmentDirection, fNSegment, segOrder, fsegment);
91 // EAxis axis[3] = {kXAxis, kYAxis, kZAxis};
92
93 G4String layerName[2] = { boxName + "_1", boxName + "_2" };
94 G4VSolid* layerSolid[2];
95 G4LogicalVolume* layerLogical[2];
96
97 //-- fisrt nested layer (replicated to x direction)
98 if(verboseLevel > 9)
99 G4cout << "layer 1 :" << G4endl;
100 layerSolid[0] =
101 new G4Box(layerName[0], fSize[0] / fNSegment[0], fSize[1], fSize[2]);
102 layerLogical[0] = new G4LogicalVolume(layerSolid[0], 0, layerName[0]);
103 if(fNSegment[0] > 1)
104 {
105 if(verboseLevel > 9)
106 G4cout << "G4ScoringBox::Construct() : Replicate to x direction"
107 << G4endl;
108 if(G4ScoringManager::GetReplicaLevel() > 0)
109 {
110 new G4PVReplica(layerName[0], layerLogical[0], boxLogical, kXAxis,
111 fNSegment[0], fSize[0] / fNSegment[0] * 2.);
112 }
113 else
114 {
115 new G4PVDivision(layerName[0], layerLogical[0], boxLogical, kXAxis,
116 fNSegment[0], 0.);
117 }
118 }
119 else if(fNSegment[0] == 1)
120 {
121 if(verboseLevel > 9)
122 G4cout << "G4ScoringBox::Construct() : Placement" << G4endl;
123 new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), layerLogical[0],
124 layerName[0], boxLogical, false, 0);
125 }
126 else
127 G4cerr << "ERROR : G4ScoringBox::SetupGeometry() : invalid parameter ("
128 << fNSegment[0] << ") "
129 << "in placement of the first nested layer." << G4endl;
130
131 if(verboseLevel > 9)
132 {
133 G4cout << fSize[0] / fNSegment[0] << ", " << fSize[1] << ", " << fSize[2]
134 << G4endl;
135 G4cout << layerName[0] << ": kXAxis, " << fNSegment[0] << ", "
136 << 2. * fSize[0] / fNSegment[0] << G4endl;
137 }
138
139 // second nested layer (replicated to y direction)
140 if(verboseLevel > 9)
141 G4cout << "layer 2 :" << G4endl;
142 layerSolid[1] = new G4Box(layerName[1], fSize[0] / fNSegment[0],
143 fSize[1] / fNSegment[1], fSize[2]);
144 layerLogical[1] = new G4LogicalVolume(layerSolid[1], 0, layerName[1]);
145 if(fNSegment[1] > 1)
146 {
147 if(verboseLevel > 9)
148 G4cout << "G4ScoringBox::Construct() : Replicate to y direction"
149 << G4endl;
150 if(G4ScoringManager::GetReplicaLevel() > 1)
151 {
152 new G4PVReplica(layerName[1], layerLogical[1], layerLogical[0], kYAxis,
153 fNSegment[1], fSize[1] / fNSegment[1] * 2.);
154 }
155 else
156 {
157 new G4PVDivision(layerName[1], layerLogical[1], layerLogical[0], kYAxis,
158 fNSegment[1], 0.);
159 }
160 }
161 else if(fNSegment[1] == 1)
162 {
163 if(verboseLevel > 9)
164 G4cout << "G4ScoringBox::Construct() : Placement" << G4endl;
165 new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), layerLogical[1],
166 layerName[1], layerLogical[0], false, 0);
167 }
168 else
169 G4cerr << "ERROR : G4ScoringBox::SetupGeometry() : invalid parameter ("
170 << fNSegment[1] << ") "
171 << "in placement of the second nested layer." << G4endl;
172
173 if(verboseLevel > 9)
174 {
175 G4cout << fSize[0] / fNSegment[0] << ", " << fSize[1] / fNSegment[1] << ", "
176 << fSize[2] << G4endl;
177 G4cout << layerName[1] << ": kYAxis, " << fNSegment[1] << ", "
178 << 2. * fSize[1] / fNSegment[1] << G4endl;
179 }
180
181 // mesh elements (replicated to z direction)
182 if(verboseLevel > 9)
183 G4cout << "mesh elements :" << G4endl;
184 G4String elementName = boxName + "_3";
185 G4VSolid* elementSolid =
186 new G4Box(elementName, fSize[0] / fNSegment[0], fSize[1] / fNSegment[1],
187 fSize[2] / fNSegment[2]);
188 fMeshElementLogical = new G4LogicalVolume(elementSolid, 0, elementName);
189 if(fNSegment[2] > 1)
190 {
191 if(verboseLevel > 9)
192 G4cout << "G4ScoringBox::Construct() : Replicate to z direction"
193 << G4endl;
194
195 if(G4ScoringManager::GetReplicaLevel() > 2)
196 {
197 new G4PVReplica(elementName, fMeshElementLogical, layerLogical[1], kZAxis,
198 fNSegment[2], 2. * fSize[2] / fNSegment[2]);
199 }
200 else
201 {
202 new G4PVDivision(elementName, fMeshElementLogical, layerLogical[1],
203 kZAxis, fNSegment[2], 0.);
204 }
205 }
206 else if(fNSegment[2] == 1)
207 {
208 if(verboseLevel > 9)
209 G4cout << "G4ScoringBox::Construct() : Placement" << G4endl;
210 new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), fMeshElementLogical,
211 elementName, layerLogical[1], false, 0);
212 }
213 else
214 G4cerr << "ERROR : G4ScoringBox::SetupGeometry() : "
215 << "invalid parameter (" << fNSegment[2] << ") "
216 << "in mesh element placement." << G4endl;
217
218 if(verboseLevel > 9)
219 {
220 G4cout << fSize[0] / fNSegment[0] << ", " << fSize[1] / fNSegment[1] << ", "
221 << fSize[2] / fNSegment[2] << G4endl;
222 G4cout << elementName << ": kZAxis, " << fNSegment[2] << ", "
223 << 2. * fSize[2] / fNSegment[2] << G4endl;
224 }
225
226 // set the sensitive detector
227 fMeshElementLogical->SetSensitiveDetector(fMFD);
228
229 // vis. attributes
230 G4VisAttributes* visatt = new G4VisAttributes(G4Colour(.5, .5, .5));
231 visatt->SetVisibility(false);
232 layerLogical[0]->SetVisAttributes(visatt);
233 layerLogical[1]->SetVisAttributes(visatt);
234 visatt->SetVisibility(true);
235 fMeshElementLogical->SetVisAttributes(visatt);
236}
237
238void G4ScoringBox::List() const
239{
240 G4cout << "G4ScoringBox : " << fWorldName << " --- Shape: Box mesh" << G4endl;
241 G4cout << " Size (x, y, z): (" << fSize[0] / cm << ", " << fSize[1] / cm
242 << ", " << fSize[2] / cm << ") [cm]" << G4endl;
243
244 G4VScoringMesh::List();
245}
246
247void G4ScoringBox::Draw(RunScore* map, G4VScoreColorMap* colorMap, G4int axflg)
248{
249 G4VVisManager* pVisManager = G4VVisManager::GetConcreteInstance();
250 if(pVisManager)
251 {
252 // cell vectors
253 std::vector<std::vector<std::vector<double>>> cell; // cell[X][Y][Z]
254 std::vector<double> ez;
255 for(int z = 0; z < fNSegment[2]; z++)
256 ez.push_back(0.);
257 std::vector<std::vector<double>> eyz;
258 for(int y = 0; y < fNSegment[1]; y++)
259 eyz.push_back(ez);
260 for(int x = 0; x < fNSegment[0]; x++)
261 cell.push_back(eyz);
262
263 std::vector<std::vector<double>> xycell; // xycell[X][Y]
264 std::vector<double> ey;
265 for(int y = 0; y < fNSegment[1]; y++)
266 ey.push_back(0.);
267 for(int x = 0; x < fNSegment[0]; x++)
268 xycell.push_back(ey);
269
270 std::vector<std::vector<double>> yzcell; // yzcell[Y][Z]
271 for(int y = 0; y < fNSegment[1]; y++)
272 yzcell.push_back(ez);
273
274 std::vector<std::vector<double>> xzcell; // xzcell[X][Z]
275 for(int x = 0; x < fNSegment[0]; x++)
276 xzcell.push_back(ez);
277
278 // projections
279 G4int q[3];
280 std::map<G4int, G4StatDouble*>::iterator itr = map->GetMap()->begin();
281 for(; itr != map->GetMap()->end(); itr++)
282 {
283 GetXYZ(itr->first, q);
284
285 xycell[q[0]][q[1]] += (itr->second->sum_wx()) / fDrawUnitValue;
286 yzcell[q[1]][q[2]] += (itr->second->sum_wx()) / fDrawUnitValue;
287 xzcell[q[0]][q[2]] += (itr->second->sum_wx()) / fDrawUnitValue;
288 }
289
290 // search max. & min. values in each slice
291 G4double xymin = DBL_MAX, yzmin = DBL_MAX, xzmin = DBL_MAX;
292 G4double xymax = 0., yzmax = 0., xzmax = 0.;
293 for(int x = 0; x < fNSegment[0]; x++)
294 {
295 for(int y = 0; y < fNSegment[1]; y++)
296 {
297 if(xymin > xycell[x][y])
298 xymin = xycell[x][y];
299 if(xymax < xycell[x][y])
300 xymax = xycell[x][y];
301 }
302 for(int z = 0; z < fNSegment[2]; z++)
303 {
304 if(xzmin > xzcell[x][z])
305 xzmin = xzcell[x][z];
306 if(xzmax < xzcell[x][z])
307 xzmax = xzcell[x][z];
308 }
309 }
310 for(int y = 0; y < fNSegment[1]; y++)
311 {
312 for(int z = 0; z < fNSegment[2]; z++)
313 {
314 if(yzmin > yzcell[y][z])
315 yzmin = yzcell[y][z];
316 if(yzmax < yzcell[y][z])
317 yzmax = yzcell[y][z];
318 }
319 }
320
321 G4VisAttributes att;
322 att.SetForceSolid(true);
323 att.SetForceAuxEdgeVisible(true);
324 G4double thick = 0.01;
325
326 G4Scale3D scale;
327 if(axflg / 100 == 1)
328 {
329 pVisManager->BeginDraw();
330
331 // xy plane
332 if(colorMap->IfFloatMinMax())
333 {
334 colorMap->SetMinMax(xymin, xymax);
335 }
336 G4ThreeVector zhalf(0., 0., fSize[2] / fNSegment[2] - thick);
337 for(int x = 0; x < fNSegment[0]; x++)
338 {
339 for(int y = 0; y < fNSegment[1]; y++)
340 {
341 G4ThreeVector pos(GetReplicaPosition(x, y, 0) - zhalf);
342 G4ThreeVector pos2(GetReplicaPosition(x, y, fNSegment[2] - 1) +
343 zhalf);
344 G4Transform3D trans, trans2;
345 if(fRotationMatrix)
346 {
347 trans = G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos);
348 trans = G4Translate3D(fCenterPosition) * trans;
349 trans2 =
350 G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos2);
351 trans2 = G4Translate3D(fCenterPosition) * trans2;
352 }
353 else
354 {
355 trans = G4Translate3D(pos) * G4Translate3D(fCenterPosition);
356 trans2 = G4Translate3D(pos2) * G4Translate3D(fCenterPosition);
357 }
358 G4double c[4];
359 colorMap->GetMapColor(xycell[x][y], c);
360 att.SetColour(c[0], c[1], c[2]); //, c[3]);
361
362 G4Box xyplate("xy", fSize[0] / fNSegment[0], fSize[1] / fNSegment[1],
363 thick);
364 G4Polyhedron* poly = xyplate.GetPolyhedron();
365 poly->Transform(trans);
366 poly->SetVisAttributes(&att);
367 pVisManager->Draw(*poly);
368
369 G4Box xyplate2 = xyplate;
370 G4Polyhedron* poly2 = xyplate2.GetPolyhedron();
371 poly2->Transform(trans2);
372 poly2->SetVisAttributes(&att);
373 pVisManager->Draw(*poly2);
374
375 /*
376 G4double nodes[][3] =
377 {{-fSize[0]/fNSegment[0], -fSize[1]/fNSegment[1], 0.},
378 { fSize[0]/fNSegment[0], -fSize[1]/fNSegment[1], 0.},
379 { fSize[0]/fNSegment[0], fSize[1]/fNSegment[1], 0.},
380 {-fSize[0]/fNSegment[0], fSize[1]/fNSegment[1], 0.}};
381 G4int facets[][4] = {{4, 3, 2, 1}};
382 G4int facets2[][4] = {{1, 2, 3, 4}};
383
384 G4Polyhedron poly, poly2;
385 poly.createPolyhedron(4, 1, nodes, facets);
386 poly.Transform(trans);
387 poly.SetVisAttributes(att);
388 pVisManager->Draw(poly);
389
390 poly2.createPolyhedron(4, 1, nodes, facets2);
391 poly2.Transform(trans2);
392 poly2.SetVisAttributes(att);
393 pVisManager->Draw(poly2);
394 */
395 }
396 }
397 pVisManager->EndDraw();
398 }
399 axflg = axflg % 100;
400 if(axflg / 10 == 1)
401 {
402 pVisManager->BeginDraw();
403
404 // yz plane
405 if(colorMap->IfFloatMinMax())
406 {
407 colorMap->SetMinMax(yzmin, yzmax);
408 }
409 G4ThreeVector xhalf(fSize[0] / fNSegment[0] - thick, 0., 0.);
410 for(int y = 0; y < fNSegment[1]; y++)
411 {
412 for(int z = 0; z < fNSegment[2]; z++)
413 {
414 G4ThreeVector pos(GetReplicaPosition(0, y, z) - xhalf);
415 G4ThreeVector pos2(GetReplicaPosition(fNSegment[0] - 1, y, z) +
416 xhalf);
417 G4Transform3D trans, trans2;
418 if(fRotationMatrix)
419 {
420 trans = G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos);
421 trans = G4Translate3D(fCenterPosition) * trans;
422 trans2 =
423 G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos2);
424 trans2 = G4Translate3D(fCenterPosition) * trans2;
425 }
426 else
427 {
428 trans = G4Translate3D(pos) * G4Translate3D(fCenterPosition);
429 trans2 = G4Translate3D(pos2) * G4Translate3D(fCenterPosition);
430 }
431 G4double c[4];
432 colorMap->GetMapColor(yzcell[y][z], c);
433 att.SetColour(c[0], c[1], c[2]); //, c[3]);
434
435 G4Box yzplate("yz", thick, // fSize[0]/fNSegment[0]*0.001,
436 fSize[1] / fNSegment[1], fSize[2] / fNSegment[2]);
437 G4Polyhedron* poly = yzplate.GetPolyhedron();
438 poly->Transform(trans);
439 poly->SetVisAttributes(&att);
440 pVisManager->Draw(*poly);
441
442 G4Box yzplate2 = yzplate;
443 G4Polyhedron* poly2 = yzplate2.GetPolyhedron();
444 poly2->Transform(trans2);
445 poly2->SetVisAttributes(&att);
446 pVisManager->Draw(*poly2);
447
448 /*
449 G4double nodes[][3] =
450 {{0., -fSize[1]/fNSegment[1], -fSize[2]/fNSegment[2]},
451 {0., fSize[1]/fNSegment[1], -fSize[2]/fNSegment[2]},
452 {0., fSize[1]/fNSegment[1], fSize[2]/fNSegment[2]},
453 {0., -fSize[1]/fNSegment[1], fSize[2]/fNSegment[2]}};
454 G4int facets[][4] = {{4, 3, 2, 1}};
455 G4int facets2[][4] = {{1, 2, 3, 4}};
456
457 G4Polyhedron poly, poly2;
458 poly.createPolyhedron(4, 1, nodes, facets);
459 poly.Transform(trans);
460 poly.SetVisAttributes(att);
461 pVisManager->Draw(poly);
462
463 poly2.createPolyhedron(4, 1, nodes, facets2);
464 poly2.Transform(trans2);
465 poly2.SetVisAttributes(att);
466 pVisManager->Draw(poly2);
467 */
468 }
469 }
470 pVisManager->EndDraw();
471 }
472 axflg = axflg % 10;
473 if(axflg == 1)
474 {
475 pVisManager->BeginDraw();
476
477 // xz plane
478 if(colorMap->IfFloatMinMax())
479 {
480 colorMap->SetMinMax(xzmin, xzmax);
481 }
482 G4ThreeVector yhalf(0., fSize[1] / fNSegment[1] - thick, 0.);
483 for(int x = 0; x < fNSegment[0]; x++)
484 {
485 for(int z = 0; z < fNSegment[2]; z++)
486 {
487 G4ThreeVector pos(GetReplicaPosition(x, 0, z) - yhalf);
488 G4ThreeVector pos2(GetReplicaPosition(x, fNSegment[1] - 1, z) +
489 yhalf);
490 G4Transform3D trans, trans2;
491 if(fRotationMatrix)
492 {
493 trans = G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos);
494 trans = G4Translate3D(fCenterPosition) * trans;
495 trans2 =
496 G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos2);
497 trans2 = G4Translate3D(fCenterPosition) * trans2;
498 }
499 else
500 {
501 trans = G4Translate3D(pos) * G4Translate3D(fCenterPosition);
502 trans2 = G4Translate3D(pos2) * G4Translate3D(fCenterPosition);
503 }
504 G4double c[4];
505 colorMap->GetMapColor(xzcell[x][z], c);
506 att.SetColour(c[0], c[1], c[2]); //, c[3]);
507
508 G4Box xzplate("xz", fSize[0] / fNSegment[0],
509 thick, // fSize[1]/fNSegment[1]*0.001,
510 fSize[2] / fNSegment[2]);
511 G4Polyhedron* poly = xzplate.GetPolyhedron();
512 poly->Transform(trans);
513 poly->SetVisAttributes(&att);
514 pVisManager->Draw(*poly);
515
516 G4Box xzplate2 = xzplate;
517 G4Polyhedron* poly2 = xzplate2.GetPolyhedron();
518 poly2->Transform(trans2);
519 poly2->SetVisAttributes(&att);
520 pVisManager->Draw(*poly2);
521
522 /*
523 G4double nodes[][3] =
524 {{-fSize[1]/fNSegment[1], 0., -fSize[2]/fNSegment[2]},
525 { fSize[1]/fNSegment[1], 0., -fSize[2]/fNSegment[2]},
526 { fSize[1]/fNSegment[1], 0., fSize[2]/fNSegment[2]},
527 {-fSize[1]/fNSegment[1], 0., fSize[2]/fNSegment[2]}};
528 G4int facets[][4] = {{1, 2, 3, 4}};
529 G4int facets2[][4] = {{4, 3, 2, 1}};
530
531 G4Polyhedron poly, poly2;
532 poly.createPolyhedron(4, 1, nodes, facets);
533 poly.Transform(trans);
534 poly.SetVisAttributes(att);
535 pVisManager->Draw(poly);
536
537 poly2.createPolyhedron(4, 1, nodes, facets2);
538 poly2.Transform(trans2);
539 poly2.SetVisAttributes(att);
540 pVisManager->Draw(poly2);
541 */
542 }
543 }
544 pVisManager->EndDraw();
545 }
546 }
547 colorMap->SetPSUnit(fDrawUnit);
548 colorMap->SetPSName(fDrawPSName);
549 colorMap->DrawColorChart();
550}
551
552G4ThreeVector G4ScoringBox::GetReplicaPosition(G4int x, G4int y, G4int z)
553{
554 G4ThreeVector width(fSize[0] / fNSegment[0], fSize[1] / fNSegment[1],
555 fSize[2] / fNSegment[2]);
556 G4ThreeVector pos(-fSize[0] + 2 * (x + 0.5) * width.x(),
557 -fSize[1] + 2 * (y + 0.5) * width.y(),
558 -fSize[2] + 2 * (z + 0.5) * width.z());
559
560 return pos;
561}
562
563void G4ScoringBox::GetXYZ(G4int index, G4int q[3]) const
564{
565 q[0] = index / (fNSegment[2] * fNSegment[1]);
566 q[1] = (index - q[0] * fNSegment[2] * fNSegment[1]) / fNSegment[2];
567 q[2] = index - q[1] * fNSegment[2] - q[0] * fNSegment[2] * fNSegment[1];
568}
569
570G4int G4ScoringBox::GetIndex(G4int x, G4int y, G4int z) const
571{
572 return x + y * fNSegment[0] + z * fNSegment[0] * fNSegment[1];
573}
574
575void G4ScoringBox::DrawColumn(RunScore* map, G4VScoreColorMap* colorMap,
576 G4int idxProj, G4int idxColumn)
577{
578 G4int iColumn[3] = { 2, 0, 1 };
579 if(idxColumn < 0 || idxColumn >= fNSegment[iColumn[idxProj]])
580 {
581 G4cerr << "ERROR : Column number " << idxColumn
582 << " is out of scoring mesh [0," << fNSegment[iColumn[idxProj]] - 1
583 << "]. Method ignored." << G4endl;
584 return;
585 }
586 G4VVisManager* pVisManager = G4VVisManager::GetConcreteInstance();
587 if(pVisManager)
588 {
589 pVisManager->BeginDraw();
590
591 // cell vectors
592 std::vector<std::vector<std::vector<double>>> cell; // cell[X][Y][Z]
593 std::vector<double> ez;
594 for(int z = 0; z < fNSegment[2]; z++)
595 ez.push_back(0.);
596 std::vector<std::vector<double>> eyz;
597 for(int y = 0; y < fNSegment[1]; y++)
598 eyz.push_back(ez);
599 for(int x = 0; x < fNSegment[0]; x++)
600 cell.push_back(eyz);
601
602 std::vector<std::vector<double>> xycell; // xycell[X][Y]
603 std::vector<double> ey;
604 for(int y = 0; y < fNSegment[1]; y++)
605 ey.push_back(0.);
606 for(int x = 0; x < fNSegment[0]; x++)
607 xycell.push_back(ey);
608
609 std::vector<std::vector<double>> yzcell; // yzcell[Y][Z]
610 for(int y = 0; y < fNSegment[1]; y++)
611 yzcell.push_back(ez);
612
613 std::vector<std::vector<double>> xzcell; // xzcell[X][Z]
614 for(int x = 0; x < fNSegment[0]; x++)
615 xzcell.push_back(ez);
616
617 // projections
618 G4int q[3];
619 std::map<G4int, G4StatDouble*>::iterator itr = map->GetMap()->begin();
620 for(; itr != map->GetMap()->end(); itr++)
621 {
622 GetXYZ(itr->first, q);
623
624 if(idxProj == 0 && q[2] == idxColumn)
625 { // xy plane
626 xycell[q[0]][q[1]] += (itr->second->sum_wx()) / fDrawUnitValue;
627 }
628 if(idxProj == 1 && q[0] == idxColumn)
629 { // yz plane
630 yzcell[q[1]][q[2]] += (itr->second->sum_wx()) / fDrawUnitValue;
631 }
632 if(idxProj == 2 && q[1] == idxColumn)
633 { // zx plane
634 xzcell[q[0]][q[2]] += (itr->second->sum_wx()) / fDrawUnitValue;
635 }
636 }
637
638 // search max. & min. values in each slice
639 G4double xymin = DBL_MAX, yzmin = DBL_MAX, xzmin = DBL_MAX;
640 G4double xymax = 0., yzmax = 0., xzmax = 0.;
641 for(int x = 0; x < fNSegment[0]; x++)
642 {
643 for(int y = 0; y < fNSegment[1]; y++)
644 {
645 if(xymin > xycell[x][y])
646 xymin = xycell[x][y];
647 if(xymax < xycell[x][y])
648 xymax = xycell[x][y];
649 }
650 for(int z = 0; z < fNSegment[2]; z++)
651 {
652 if(xzmin > xzcell[x][z])
653 xzmin = xzcell[x][z];
654 if(xzmax < xzcell[x][z])
655 xzmax = xzcell[x][z];
656 }
657 }
658 for(int y = 0; y < fNSegment[1]; y++)
659 {
660 for(int z = 0; z < fNSegment[2]; z++)
661 {
662 if(yzmin > yzcell[y][z])
663 yzmin = yzcell[y][z];
664 if(yzmax < yzcell[y][z])
665 yzmax = yzcell[y][z];
666 }
667 }
668
669 G4VisAttributes att;
670 att.SetForceSolid(true);
671 att.SetForceAuxEdgeVisible(true);
672
673 G4Scale3D scale;
674 // xy plane
675 if(idxProj == 0)
676 {
677 if(colorMap->IfFloatMinMax())
678 {
679 colorMap->SetMinMax(xymin, xymax);
680 }
681 for(int x = 0; x < fNSegment[0]; x++)
682 {
683 for(int y = 0; y < fNSegment[1]; y++)
684 {
685 G4Box xyplate("xy", fSize[0] / fNSegment[0], fSize[1] / fNSegment[1],
686 fSize[2] / fNSegment[2]);
687
688 G4ThreeVector pos(GetReplicaPosition(x, y, idxColumn));
689 G4Transform3D trans;
690 if(fRotationMatrix)
691 {
692 trans = G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos);
693 trans = G4Translate3D(fCenterPosition) * trans;
694 }
695 else
696 {
697 trans = G4Translate3D(pos) * G4Translate3D(fCenterPosition);
698 }
699 G4double c[4];
700 colorMap->GetMapColor(xycell[x][y], c);
701 att.SetColour(c[0], c[1], c[2]);
702
703 G4Polyhedron* poly = xyplate.GetPolyhedron();
704 poly->Transform(trans);
705 poly->SetVisAttributes(att);
706 pVisManager->Draw(*poly);
707 }
708 }
709 }
710 else
711 // yz plane
712 if(idxProj == 1)
713 {
714 if(colorMap->IfFloatMinMax())
715 {
716 colorMap->SetMinMax(yzmin, yzmax);
717 }
718 for(int y = 0; y < fNSegment[1]; y++)
719 {
720 for(int z = 0; z < fNSegment[2]; z++)
721 {
722 G4Box yzplate("yz", fSize[0] / fNSegment[0], fSize[1] / fNSegment[1],
723 fSize[2] / fNSegment[2]);
724
725 G4ThreeVector pos(GetReplicaPosition(idxColumn, y, z));
726 G4Transform3D trans;
727 if(fRotationMatrix)
728 {
729 trans = G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos);
730 trans = G4Translate3D(fCenterPosition) * trans;
731 }
732 else
733 {
734 trans = G4Translate3D(pos) * G4Translate3D(fCenterPosition);
735 }
736 G4double c[4];
737 colorMap->GetMapColor(yzcell[y][z], c);
738 att.SetColour(c[0], c[1], c[2]); //, c[3]);
739
740 G4Polyhedron* poly = yzplate.GetPolyhedron();
741 poly->Transform(trans);
742 poly->SetVisAttributes(att);
743 pVisManager->Draw(*poly);
744 }
745 }
746 }
747 else
748 // xz plane
749 if(idxProj == 2)
750 {
751 if(colorMap->IfFloatMinMax())
752 {
753 colorMap->SetMinMax(xzmin, xzmax);
754 }
755 for(int x = 0; x < fNSegment[0]; x++)
756 {
757 for(int z = 0; z < fNSegment[2]; z++)
758 {
759 G4Box xzplate("xz", fSize[0] / fNSegment[0], fSize[1] / fNSegment[1],
760 fSize[2] / fNSegment[2]);
761
762 G4ThreeVector pos(GetReplicaPosition(x, idxColumn, z));
763 G4Transform3D trans;
764 if(fRotationMatrix)
765 {
766 trans = G4Rotate3D(*fRotationMatrix).inverse() * G4Translate3D(pos);
767 trans = G4Translate3D(fCenterPosition) * trans;
768 }
769 else
770 {
771 trans = G4Translate3D(pos) * G4Translate3D(fCenterPosition);
772 }
773 G4double c[4];
774 colorMap->GetMapColor(xzcell[x][z], c);
775 att.SetColour(c[0], c[1], c[2]); //, c[3]);
776
777 G4Polyhedron* poly = xzplate.GetPolyhedron();
778 poly->Transform(trans);
779 poly->SetVisAttributes(att);
780 pVisManager->Draw(*poly);
781 }
782 }
783 }
784 pVisManager->EndDraw();
785 }
786
787 colorMap->SetPSUnit(fDrawUnit);
788 colorMap->SetPSName(fDrawPSName);
789 colorMap->DrawColorChart();
790}
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66
cm
static constexpr double cm
Definition: G4SIunits.hh:99
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4Translate3D
HepGeom::Translate3D G4Translate3D
Definition: G4Transform3D.hh:41
G4Rotate3D
HepGeom::Rotate3D G4Rotate3D
Definition: G4Transform3D.hh:36
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4cerr
G4GLOB_DLL std::ostream G4cerr
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4Box
Definition: G4Box.hh:56
G4CSGSolid::GetPolyhedron
virtual G4Polyhedron * GetPolyhedron() const
Definition: G4CSGSolid.cc:129
G4LogicalVolume::SetVisAttributes
void SetVisAttributes(const G4VisAttributes *pVA)
G4LogicalVolume::SetSensitiveDetector
void SetSensitiveDetector(G4VSensitiveDetector *pSDetector)
Definition: G4LogicalVolume.cc:478
G4ScoringBox::GetXYZ
void GetXYZ(G4int index, G4int q[3]) const
Definition: G4ScoringBox.cc:563
G4ScoringBox::Draw
void Draw(RunScore *map, G4VScoreColorMap *colorMap, G4int axflg=111)
Definition: G4ScoringBox.cc:247
G4ScoringBox::GetReplicaPosition
G4ThreeVector GetReplicaPosition(G4int x, G4int y, G4int z)
Definition: G4ScoringBox.cc:552
G4ScoringBox::List
void List() const
Definition: G4ScoringBox.cc:238
G4ScoringBox::SetupGeometry
virtual void SetupGeometry(G4VPhysicalVolume *fWorldPhys)
Definition: G4ScoringBox.cc:66
G4ScoringBox::G4ScoringBox
G4ScoringBox(G4String wName)
Definition: G4ScoringBox.cc:54
G4ScoringBox::DrawColumn
void DrawColumn(RunScore *map, G4VScoreColorMap *colorMap, G4int idxProj, G4int idxColumn)
Definition: G4ScoringBox.cc:575
G4ScoringBox::GetIndex
G4int GetIndex(G4int x, G4int y, G4int z) const
Definition: G4ScoringBox.cc:570
G4ScoringManager::GetReplicaLevel
static G4int GetReplicaLevel()
Definition: G4ScoringManager.cc:84
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4VScoreColorMap::DrawColorChart
virtual void DrawColorChart(G4int nPoint=5)
Definition: G4VScoreColorMap.cc:54
G4VScoreColorMap::IfFloatMinMax
G4bool IfFloatMinMax() const
Definition: G4VScoreColorMap.hh:48
G4VScoreColorMap::SetMinMax
void SetMinMax(G4double minVal, G4double maxVal)
Definition: G4VScoreColorMap.hh:49
G4VScoreColorMap::SetPSUnit
void SetPSUnit(G4String &unit)
Definition: G4VScoreColorMap.hh:75
G4VScoreColorMap::GetMapColor
virtual void GetMapColor(G4double val, G4double color[4])=0
G4VScoreColorMap::SetPSName
void SetPSName(G4String &psName)
Definition: G4VScoreColorMap.hh:76
G4VScoringMesh::fRotationMatrix
G4RotationMatrix * fRotationMatrix
Definition: G4VScoringMesh.hh:208
G4VScoringMesh::List
virtual void List() const
Definition: G4VScoringMesh.cc:308
G4VScoringMesh::fDrawUnitValue
G4double fDrawUnitValue
Definition: G4VScoringMesh.hh:220
G4VScoringMesh::fDrawPSName
G4String fDrawPSName
Definition: G4VScoringMesh.hh:221
G4VScoringMesh::fMFD
G4MultiFunctionalDetector * fMFD
Definition: G4VScoringMesh.hh:212
G4VScoringMesh::fDivisionAxisNames
G4String fDivisionAxisNames[3]
Definition: G4VScoringMesh.hh:223
G4VScoringMesh::fMeshElementLogical
G4LogicalVolume * fMeshElementLogical
Definition: G4VScoringMesh.hh:225
G4VScoringMesh::fSize
G4double fSize[3]
Definition: G4VScoringMesh.hh:205
G4VScoringMesh::fCenterPosition
G4ThreeVector fCenterPosition
Definition: G4VScoringMesh.hh:207
G4VVisManager::EndDraw
virtual void EndDraw()=0
G4VVisManager::GetConcreteInstance
static G4VVisManager * GetConcreteInstance()
Definition: G4VVisManager.cc:41
G4VVisManager::Draw
virtual void Draw(const G4Circle &, const G4Transform3D &objectTransformation=G4Transform3D())=0
G4VVisManager::BeginDraw
virtual void BeginDraw(const G4Transform3D &objectTransformation=G4Transform3D())=0
G4VisAttributes::SetColour
void SetColour(const G4Colour &)
G4VisAttributes::SetVisibility
void SetVisibility(G4bool=true)
G4VisAttributes::SetForceAuxEdgeVisible
void SetForceAuxEdgeVisible(G4bool=true)
Definition: G4VisAttributes.cc:205
G4VisAttributes::SetForceSolid
void SetForceSolid(G4bool=true)
Definition: G4VisAttributes.cc:175
G4Visible::SetVisAttributes
void SetVisAttributes(const G4VisAttributes *)
Definition: G4Visible.cc:96
HepGeom::Transform3D::inverse
Transform3D inverse() const
Definition: Transform3D.cc:141
HepPolyhedron::Transform
HepPolyhedron & Transform(const G4Transform3D &t)
Definition: HepPolyhedron.cc:1212
kYAxis
@ kYAxis
Definition: geomdefs.hh:56
kXAxis
@ kXAxis
Definition: geomdefs.hh:55
kZAxis
@ kZAxis
Definition: geomdefs.hh:57
anonymous_namespace{G4QuasiElRatios.cc}::map
const G4double map
Definition: G4QuasiElRatios.cc:81
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62
Generated by doxygen 1.9.3