Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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), fSegmentDirection(-1)
56{
57 fShape = MeshShape::box;
58 fDivisionAxisNames[0] = "X";
59 fDivisionAxisNames[1] = "Y";
60 fDivisionAxisNames[2] = "Z";
61}
62
63G4ScoringBox::~G4ScoringBox()
64{
65}
66
67void 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
229void 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
240void 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
499G4ThreeVector 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
510void 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
518G4int G4ScoringBox::GetIndex(G4int x, G4int y, G4int z) const {
519 return x + y*fNSegment[0] + z*fNSegment[0]*fNSegment[1];
520}
521
522void 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
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
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:468
G4ScoringBox::Draw
void Draw(RunScore *map, G4VScoreColorMap *colorMap, G4int axflg=111)
Definition: G4ScoringBox.cc:240
G4ScoringBox::List
void List() const
Definition: G4ScoringBox.cc:229
G4ScoringBox::SetupGeometry
virtual void SetupGeometry(G4VPhysicalVolume *fWorldPhys)
Definition: G4ScoringBox.cc:67
G4ScoringBox::G4ScoringBox
G4ScoringBox(G4String wName)
Definition: G4ScoringBox.cc:54
G4ScoringBox::DrawColumn
void DrawColumn(RunScore *map, G4VScoreColorMap *colorMap, G4int idxProj, G4int idxColumn)
Definition: G4ScoringBox.cc:522
G4ScoringManager::GetReplicaLevel
static G4int GetReplicaLevel()
Definition: G4ScoringManager.cc:79
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4VScoreColorMap::DrawColorChart
virtual void DrawColorChart(G4int nPoint=5)
Definition: G4VScoreColorMap.cc:49
G4VScoreColorMap::IfFloatMinMax
G4bool IfFloatMinMax() const
Definition: G4VScoreColorMap.hh:50
G4VScoreColorMap::SetMinMax
void SetMinMax(G4double minVal, G4double maxVal)
Definition: G4VScoreColorMap.hh:52
G4VScoreColorMap::SetPSUnit
void SetPSUnit(G4String &unit)
Definition: G4VScoreColorMap.hh:76
G4VScoreColorMap::GetMapColor
virtual void GetMapColor(G4double val, G4double color[4])=0
G4VScoreColorMap::SetPSName
void SetPSName(G4String &psName)
Definition: G4VScoreColorMap.hh:77
G4VScoringMesh::fRotationMatrix
G4RotationMatrix * fRotationMatrix
Definition: G4VScoringMesh.hh:190
G4VScoringMesh::List
virtual void List() const
Definition: G4VScoringMesh.cc:236
G4VScoringMesh::fDrawUnitValue
G4double fDrawUnitValue
Definition: G4VScoringMesh.hh:202
G4VScoringMesh::fDrawPSName
G4String fDrawPSName
Definition: G4VScoringMesh.hh:203
G4VScoringMesh::fMFD
G4MultiFunctionalDetector * fMFD
Definition: G4VScoringMesh.hh:194
G4VScoringMesh::fDivisionAxisNames
G4String fDivisionAxisNames[3]
Definition: G4VScoringMesh.hh:205
G4VScoringMesh::fMeshElementLogical
G4LogicalVolume * fMeshElementLogical
Definition: G4VScoringMesh.hh:207
G4VScoringMesh::fSize
G4double fSize[3]
Definition: G4VScoringMesh.hh:188
G4VScoringMesh::fCenterPosition
G4ThreeVector fCenterPosition
Definition: G4VScoringMesh.hh:189
G4VTHitsMap::GetMap
Map_t * GetMap() const
Definition: G4THitsMap.hh:151
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:209
G4VisAttributes::SetForceSolid
void SetForceSolid(G4bool=true)
Definition: G4VisAttributes.cc:180
G4Visible::SetVisAttributes
void SetVisAttributes(const G4VisAttributes *)
Definition: G4Visible.cc:79
HepGeom::Transform3D::inverse
Transform3D inverse() const
Definition: Transform3D.cc:141
HepPolyhedron::Transform
HepPolyhedron & Transform(const G4Transform3D &t)
Definition: HepPolyhedron.cc:893
kYAxis
@ kYAxis
Definition: geomdefs.hh:56
kXAxis
@ kXAxis
Definition: geomdefs.hh:55
kZAxis
@ kZAxis
Definition: geomdefs.hh:57
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62