ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ScoringCylinder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ScoringCylinder.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 "G4ScoringCylinder.hh"
30 
31 #include "G4VPhysicalVolume.hh"
32 #include "G4Tubs.hh"
33 #include "G4LogicalVolume.hh"
34 #include "G4VPhysicalVolume.hh"
35 #include "G4PVPlacement.hh"
36 #include "G4PVReplica.hh"
37 #include "G4PVDivision.hh"
38 #include "G4VisAttributes.hh"
39 #include "G4VVisManager.hh"
40 #include "G4VScoreColorMap.hh"
41 
42 #include "G4SDManager.hh"
43 #include "G4MultiFunctionalDetector.hh"
44 #include "G4SDParticleFilter.hh"
45 #include "G4VPrimitiveScorer.hh"
46 #include "G4PSEnergyDeposit.hh"
47 #include "G4PSTrackLength.hh"
48 #include "G4PSNofStep.hh"
49 #include "G4ScoringManager.hh"
50 #include "G4StatDouble.hh"
51 
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 
55 G4ScoringCylinder::G4ScoringCylinder(G4String wName)
56  :G4VScoringMesh(wName)
57 {
58  fShape = MeshShape::cylinder;
59 
60  fDivisionAxisNames[0] = "Z";
61  fDivisionAxisNames[1] = "PHI";
62  fDivisionAxisNames[2] = "R";
63 }
64 
65 G4ScoringCylinder::~G4ScoringCylinder()
66 {;}
67 
68 void G4ScoringCylinder::SetupGeometry(G4VPhysicalVolume * fWorldPhys) {
69 
70  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::SetupGeometry() ..." << G4endl;
71 
72  // World
73  G4VPhysicalVolume * scoringWorld = fWorldPhys;
74  G4LogicalVolume * worldLogical = scoringWorld->GetLogicalVolume();
75 
76  // Scoring Mesh
77  if(verboseLevel > 9) G4cout << fWorldName << G4endl;
78  G4String tubsName = fWorldName+"_mesh";
79 
80  if(verboseLevel > 9) G4cout << "R max., Dz =: " << fSize[0] << ", " << fSize[1] << G4endl;
81  G4VSolid * tubsSolid = new G4Tubs(tubsName+"0", // name
82  0., // R min
83  fSize[0], // R max
84  fSize[1], // Dz
85  0., // starting phi
86  twopi*rad); // segment phi
87  G4LogicalVolume * tubsLogical = new G4LogicalVolume(tubsSolid, 0, tubsName);
88  new G4PVPlacement(fRotationMatrix, fCenterPosition,
89  tubsLogical, tubsName+"0", worldLogical, false, 0);
90 
91  if(verboseLevel > 9) G4cout << " # of segments : r, phi, z =: "
92  << fNSegment[IR] << ", " << fNSegment[IPHI] << ", " << fNSegment[IZ] << G4endl;
93 
94  G4String layerName[2] = {tubsName + "1", tubsName + "2"};
95  G4VSolid * layerSolid[2];
96  G4LogicalVolume * layerLogical[2];
97 
98  //-- fisrt nested layer (replicated along z direction)
99  if(verboseLevel > 9) G4cout << "layer 1 :" << G4endl;
100  layerSolid[0] = new G4Tubs(layerName[0], // name
101  0., // inner radius
102  fSize[0], // outer radius
103  fSize[1]/fNSegment[IZ], // half len. in z
104  0., // starting phi angle
105  twopi*rad); // delta angle of the segment
106  layerLogical[0] = new G4LogicalVolume(layerSolid[0], 0, layerName[0]);
107  if(fNSegment[IZ] > 1) {
108  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replicate along z direction" << G4endl;
109  if(G4ScoringManager::GetReplicaLevel()>0) {
110  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replica" << G4endl;
111  new G4PVReplica(layerName[0], layerLogical[0], tubsLogical, kZAxis, fNSegment[IZ], 2.*fSize[1]/fNSegment[IZ]);
112  } else {
113  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Division" << G4endl;
114  new G4PVDivision(layerName[0], layerLogical[0], tubsLogical, kZAxis, fNSegment[IZ], 0.);
115  }
116  } else if(fNSegment[IZ] == 1) {
117  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Placement" << G4endl;
118  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), layerLogical[0], layerName[0], tubsLogical, false, 0);
119  } else {
120  G4cerr << "G4ScoringCylinder::SetupGeometry() : invalid parameter ("
121  << fNSegment[IZ] << ") "
122  << "in placement of the first nested layer." << G4endl;
123  }
124 
125  // second nested layer (replicated along phi direction)
126  if(verboseLevel > 9) G4cout << "layer 2 :" << G4endl;
127  layerSolid[1] = new G4Tubs(layerName[1],
128  0.,
129  fSize[0],
130  fSize[1]/fNSegment[IZ],
131  0.,
132  twopi*rad/fNSegment[IPHI]);
133  layerLogical[1] = new G4LogicalVolume(layerSolid[1], 0, layerName[1]);
134  if(fNSegment[IPHI] > 1) {
135  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replicate along phi direction" << G4endl;
136  if(G4ScoringManager::GetReplicaLevel()>1) {
137  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replica" << G4endl;
138  new G4PVReplica(layerName[1], layerLogical[1], layerLogical[0], kPhi,
139  fNSegment[IPHI], twopi*rad/fNSegment[IPHI]);
140  } else {
141  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Division" << G4endl;
142  new G4PVDivision(layerName[1], layerLogical[1], layerLogical[0], kPhi, fNSegment[IPHI], 0.);
143  }
144  } else if(fNSegment[IPHI] == 1) {
145  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Placement" << G4endl;
146  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), layerLogical[1], layerName[1], layerLogical[0], false, 0);
147  } else
148  G4cerr << "ERROR : G4ScoringCylinder::SetupGeometry() : invalid parameter ("
149  << fNSegment[IPHI] << ") "
150  << "in placement of the second nested layer." << G4endl;
151 
152  // mesh elements
153  if(verboseLevel > 9) G4cout << "mesh elements :" << G4endl;
154  G4String elementName = tubsName +"3";
155  G4VSolid * elementSolid = new G4Tubs(elementName,
156  0.,
157  fSize[0]/fNSegment[IR],
158  fSize[1]/fNSegment[IZ],
159  0.,
160  twopi*rad/fNSegment[IPHI]);
161  fMeshElementLogical = new G4LogicalVolume(elementSolid, 0, elementName);
162  if(fNSegment[IR] > 1) {
163 
164  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replicate along r direction" << G4endl;
165 
166  if(G4ScoringManager::GetReplicaLevel()>2) {
167  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replica" << G4endl;
168  new G4PVReplica(elementName, fMeshElementLogical, layerLogical[1], kRho,
169  fNSegment[IR], fSize[0]/fNSegment[IR]);
170  } else {
171  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Division" << G4endl;
172  new G4PVDivision(elementName, fMeshElementLogical, layerLogical[1], kRho, fNSegment[IR], 0.);
173  }
174  } else if(fNSegment[IR] == 1) {
175  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Placement" << G4endl;
176  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), fMeshElementLogical, elementName, layerLogical[1], false, 0);
177  } else {
178  G4cerr << "G4ScoringCylinder::SetupGeometry() : "
179  << "invalid parameter (" << fNSegment[IR] << ") "
180  << "in mesh element placement." << G4endl;
181  }
182 
183  // set the sensitive detector
184  fMeshElementLogical->SetSensitiveDetector(fMFD);
185 
186 
187  // vis. attributes
188  G4VisAttributes * visatt = new G4VisAttributes(G4Colour(.5,.5,.5));
189  visatt->SetVisibility(true);
190  layerLogical[0]->SetVisAttributes(visatt);
191  layerLogical[1]->SetVisAttributes(visatt);
192  visatt = new G4VisAttributes(G4Colour(.5,.5,.5,0.01));
193  //visatt->SetForceSolid(true);
194  fMeshElementLogical->SetVisAttributes(visatt);
195 }
196 
197 void G4ScoringCylinder::List() const {
198  G4cout << "G4ScoringCylinder : " << fWorldName << " --- Shape: Cylindrical mesh" << G4endl;
199 
200  G4cout << " Size (R, Dz): ("
201  << fSize[0]/cm << ", "
202  << fSize[1]/cm << ") [cm]"
203  << G4endl;
204 
205  G4VScoringMesh::List();
206 }
207 
208 
209 void G4ScoringCylinder::Draw(RunScore * map,
210  G4VScoreColorMap* colorMap, G4int axflg) {
211 
212  G4VVisManager * pVisManager = G4VVisManager::GetConcreteInstance();
213  if(pVisManager) {
214 
215  // cell vectors
216  std::vector<double> ephi;
217  for(int phi = 0; phi < fNSegment[IPHI]; phi++) ephi.push_back(0.);
218  //-
219  std::vector<std::vector<double> > zphicell; // zphicell[Z][PHI]
220  for(int z = 0; z < fNSegment[IZ]; z++) zphicell.push_back(ephi);
221  //-
222  std::vector<std::vector<double> > rphicell; // rphicell[R][PHI]
223  for(int r = 0; r < fNSegment[IR]; r++) rphicell.push_back(ephi);
224 
225  // projections
226  G4int q[3];
227  std::map<G4int, G4StatDouble*>::iterator itr = map->GetMap()->begin();
228  for(; itr != map->GetMap()->end(); itr++) {
229  if(itr->first < 0) {
230  G4cout << itr->first << G4endl;
231  continue;
232  }
233  GetRZPhi(itr->first, q);
234 
235  zphicell[q[IZ]][q[IPHI]] += (itr->second->sum_wx())/fDrawUnitValue;
236  rphicell[q[IR]][q[IPHI]] += (itr->second->sum_wx())/fDrawUnitValue;
237  }
238 
239  // search min./max. values
240  G4double zphimin = DBL_MAX, rphimin = DBL_MAX;
241  G4double zphimax = 0., rphimax = 0.;
242  for(int iphi = 0; iphi < fNSegment[IPHI]; iphi++) {
243  for(int iz = 0; iz < fNSegment[IZ]; iz++) {
244  if(zphimin > zphicell[iz][iphi]) zphimin = zphicell[iz][iphi];
245  if(zphimax < zphicell[iz][iphi]) zphimax = zphicell[iz][iphi];
246  }
247  for(int ir = 0; ir < fNSegment[IR]; ir++) {
248  if(rphimin > rphicell[ir][iphi]) rphimin = rphicell[ir][iphi];
249  if(rphimax < rphicell[ir][iphi]) rphimax = rphicell[ir][iphi];
250  }
251  }
252 
253  G4VisAttributes att;
254  att.SetForceSolid(true);
255  att.SetForceAuxEdgeVisible(true);
256 
257 
258  G4Scale3D scale;
259  if(axflg/100==1) {
260  // rz plane
261  }
262  axflg = axflg%100;
263  if(axflg/10==1) {
264  pVisManager->BeginDraw();
265 
266  // z-phi plane
267  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(zphimin, zphimax); }
268 
269  G4double zhalf = fSize[1]/fNSegment[IZ];
270  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
271  for(int z = 0; z < fNSegment[IZ]; z++) {
272  //-
273  G4double angle = twopi/fNSegment[IPHI]*phi;
274  G4double dphi = twopi/fNSegment[IPHI];
275  G4Tubs cylinder("z-phi", // name
276  fSize[0]*0.99, fSize[0], // inner radius, outer radius
277  zhalf, // half length in z
278  angle, dphi*0.99999); // starting phi angle, delta angle
279  //-
280  G4ThreeVector zpos(0., 0., -fSize[1] + fSize[1]/fNSegment[IZ]*(1 + 2.*z));
281  G4Transform3D trans;
282  if(fRotationMatrix) {
283  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zpos);
284  trans = G4Translate3D(fCenterPosition)*trans;
285  } else {
286  trans = G4Translate3D(zpos)*G4Translate3D(fCenterPosition);
287  }
288  G4double c[4];
289  colorMap->GetMapColor(zphicell[z][phi], c);
290  att.SetColour(c[0], c[1], c[2]);//, c[3]);
291  //-
292  G4Polyhedron * poly = cylinder.GetPolyhedron();
293  poly->Transform(trans);
294  poly->SetVisAttributes(att);
295  pVisManager->Draw(*poly);
296  }
297  }
298  pVisManager->EndDraw();
299  }
300  axflg = axflg%10;
301  if(axflg==1) {
302  pVisManager->BeginDraw();
303 
304  // r-phi plane
305  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(rphimin, rphimax); }
306 
307  G4double rsize = fSize[0]/fNSegment[IR];
308  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
309  for(int r = 0; r < fNSegment[IR]; r++) {
310 
311  G4double rs[2] = {rsize*r, rsize*(r+1)};
312  G4double angle = twopi/fNSegment[IPHI]*phi;
313  G4double dphi = twopi/fNSegment[IPHI];
314  G4Tubs cylindern("z-phi", rs[0], rs[1], 0.001,
315  angle, dphi*0.99999);
316  G4Tubs cylinderp = cylindern;
317 
318  G4ThreeVector zposn(0., 0., -fSize[1]);
319  G4ThreeVector zposp(0., 0., fSize[1]);
320  G4Transform3D transn, transp;
321  if(fRotationMatrix) {
322  transn = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zposn);
323  transn = G4Translate3D(fCenterPosition)*transn;
324  transp = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zposp);
325  transp = G4Translate3D(fCenterPosition)*transp;
326  } else {
327  transn = G4Translate3D(zposn)*G4Translate3D(fCenterPosition);
328  transp = G4Translate3D(zposp)*G4Translate3D(fCenterPosition);
329  }
330  G4double c[4];
331  colorMap->GetMapColor(rphicell[r][phi], c);
332  att.SetColour(c[0], c[1], c[2]);//, c[3]);
333 
334  G4Polyhedron * polyn = cylindern.GetPolyhedron();
335  polyn->Transform(transn);
336  polyn->SetVisAttributes(att);
337  pVisManager->Draw(*polyn);
338 
339  G4Polyhedron * polyp = cylinderp.GetPolyhedron();
340  polyp->Transform(transp);
341  polyp->SetVisAttributes(att);
342  pVisManager->Draw(*polyp);
343  }
344  }
345 
346  pVisManager->EndDraw();
347  }
348 
349  colorMap->SetPSUnit(fDrawUnit);
350  colorMap->SetPSName(fDrawPSName);
351  colorMap->DrawColorChart();
352 
353  }
354 }
355 
356 void G4ScoringCylinder::DrawColumn(RunScore * map, G4VScoreColorMap* colorMap,
357  G4int idxProj, G4int idxColumn)
358 {
359  G4int projAxis = 0;
360  switch(idxProj) {
361  case 0:
362  projAxis = IR;
363  break;
364  case 1:
365  projAxis = IZ;
366  break;
367  case 2:
368  projAxis = IPHI;
369  break;
370  }
371 
372  if(idxColumn<0 || idxColumn>=fNSegment[projAxis])
373  {
374  G4cerr << "Warning : Column number " << idxColumn << " is out of scoring mesh [0," << fNSegment[projAxis]-1 <<
375  "]. Method ignored." << G4endl;
376  return;
377  }
378  G4VVisManager * pVisManager = G4VVisManager::GetConcreteInstance();
379  if(pVisManager) {
380 
381  // cell vectors
382  std::vector<std::vector<std::vector<double> > > cell; // cell[R][Z][PHI]
383  std::vector<double> ephi;
384  for(int phi = 0; phi < fNSegment[IPHI]; phi++) ephi.push_back(0.);
385  std::vector<std::vector<double> > ezphi;
386  for(int z = 0; z < fNSegment[IZ]; z++) ezphi.push_back(ephi);
387  for(int r = 0; r < fNSegment[IR]; r++) cell.push_back(ezphi);
388 
389  std::vector<std::vector<double> > rzcell; // rzcell[R][Z]
390  std::vector<double> ez;
391  for(int z = 0; z < fNSegment[IZ]; z++) ez.push_back(0.);
392  for(int r = 0; r < fNSegment[IR]; r++) rzcell.push_back(ez);
393 
394  std::vector<std::vector<double> > zphicell; // zphicell[Z][PHI]
395  for(int z = 0; z < fNSegment[IZ]; z++) zphicell.push_back(ephi);
396 
397  std::vector<std::vector<double> > rphicell; // rphicell[R][PHI]
398  for(int r = 0; r < fNSegment[IR]; r++) rphicell.push_back(ephi);
399 
400  // projections
401  G4int q[3];
402  std::map<G4int,G4StatDouble*>::iterator itr = map->GetMap()->begin();
403  for(; itr != map->GetMap()->end(); itr++) {
404  if(itr->first < 0) {
405  G4cout << itr->first << G4endl;
406  continue;
407  }
408  GetRZPhi(itr->first, q);
409 
410  if(projAxis == IR && q[IR] == idxColumn) { // zphi plane
411  zphicell[q[IZ]][q[IPHI]] += (itr->second->sum_wx())/fDrawUnitValue;
412  }
413  if(projAxis == IZ && q[IZ] == idxColumn) { // rphi plane
414  rphicell[q[IR]][q[IPHI]] += (itr->second->sum_wx())/fDrawUnitValue;
415  }
416  if(projAxis == IPHI && q[IPHI] == idxColumn) { // rz plane
417  rzcell[q[IR]][q[IZ]] += (itr->second->sum_wx())/fDrawUnitValue;
418  }
419  }
420 
421  // search min./max. values
422  G4double rzmin = DBL_MAX, zphimin = DBL_MAX, rphimin = DBL_MAX;
423  G4double rzmax = 0., zphimax = 0., rphimax = 0.;
424  for(int r = 0; r < fNSegment[IR]; r++) {
425  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
426  if(rphimin > rphicell[r][phi]) rphimin = rphicell[r][phi];
427  if(rphimax < rphicell[r][phi]) rphimax = rphicell[r][phi];
428  }
429  for(int z = 0; z < fNSegment[IZ]; z++) {
430  if(rzmin > rzcell[r][z]) rzmin = rzcell[r][z];
431  if(rzmax < rzcell[r][z]) rzmax = rzcell[r][z];
432  }
433  }
434  for(int z = 0; z < fNSegment[IZ]; z++) {
435  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
436  if(zphimin > zphicell[z][phi]) zphimin = zphicell[z][phi];
437  if(zphimax < zphicell[z][phi]) zphimax = zphicell[z][phi];
438  }
439  }
440 
441 
442  G4VisAttributes att;
443  att.SetForceSolid(true);
444  att.SetForceAuxEdgeVisible(true);
445 
446  pVisManager->BeginDraw();
447 
448  G4Scale3D scale;
449  // z-phi plane
450  if(projAxis == IR) {
451 
452  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(zphimin,zphimax); }
453 
454  G4double zhalf = fSize[1]/fNSegment[IZ];
455  G4double rsize[2] = {fSize[0]/fNSegment[IR]*idxColumn,
456  fSize[0]/fNSegment[IR]*(idxColumn+1)};
457  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
458  for(int z = 0; z < fNSegment[IZ]; z++) {
459 
460  G4double angle = twopi/fNSegment[IPHI]*phi*radian;
461  G4double dphi = twopi/fNSegment[IPHI]*radian;
462  G4Tubs cylinder("z-phi", rsize[0], rsize[1], zhalf,
463  angle, dphi*0.99999);
464 
465  G4ThreeVector zpos(0., 0., -fSize[1] + fSize[1]/fNSegment[IZ]*(1 + 2.*z));
466  G4Transform3D trans;
467  if(fRotationMatrix) {
468  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zpos);
469  trans = G4Translate3D(fCenterPosition)*trans;
470  } else {
471  trans = G4Translate3D(zpos)*G4Translate3D(fCenterPosition);
472  }
473  G4double c[4];
474  colorMap->GetMapColor(zphicell[z][phi], c);
475  att.SetColour(c[0], c[1], c[2]);//, c[3]);
476 
477  G4Polyhedron * poly = cylinder.GetPolyhedron();
478  poly->Transform(trans);
479  poly->SetVisAttributes(att);
480  pVisManager->Draw(*poly);
481  }
482  }
483 
484  // r-phi plane
485  } else if(projAxis == IZ) {
486  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(rphimin,rphimax); }
487 
488  G4double rsize = fSize[0]/fNSegment[IR];
489  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
490  for(int r = 0; r < fNSegment[IR]; r++) {
491 
492  G4double rs[2] = {rsize*r, rsize*(r+1)};
493  G4double angle = twopi/fNSegment[IPHI]*phi*radian;
494  G4double dz = fSize[1]/fNSegment[IZ];
495  G4double dphi = twopi/fNSegment[IPHI]*radian;
496  G4Tubs cylinder("r-phi", rs[0], rs[1], dz,
497  angle, dphi*0.99999);
498  G4ThreeVector zpos(0., 0.,
499  -fSize[1]+fSize[1]/fNSegment[IZ]*(idxColumn*2+1));
500  G4Transform3D trans;
501  if(fRotationMatrix) {
502  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zpos);
503  trans = G4Translate3D(fCenterPosition)*trans;
504  } else {
505  trans = G4Translate3D(zpos)*G4Translate3D(fCenterPosition);
506  }
507  G4double c[4];
508  colorMap->GetMapColor(rphicell[r][phi], c);
509  att.SetColour(c[0], c[1], c[2]);//, c[3]);
510 
511  G4Polyhedron * poly = cylinder.GetPolyhedron();
512  poly->Transform(trans);
513  poly->SetVisAttributes(att);
514  pVisManager->Draw(*poly);
515  }
516  }
517 
518  // r-z plane
519  } else if(projAxis == IPHI) {
520  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(rzmin,rzmax); }
521 
522  G4double rsize = fSize[0]/fNSegment[IR];
523  G4double zhalf = fSize[1]/fNSegment[IZ];
524  G4double angle = twopi/fNSegment[IPHI]*idxColumn*radian;
525  G4double dphi = twopi/fNSegment[IPHI]*radian;
526  for(int z = 0; z < fNSegment[IZ]; z++) {
527  for(int r = 0; r < fNSegment[IR]; r++) {
528 
529  G4double rs[2] = {rsize*r, rsize*(r+1)};
530  G4Tubs cylinder("z-phi", rs[0], rs[1], zhalf,
531  angle, dphi);
532 
533  G4ThreeVector zpos(0., 0.,
534  -fSize[1]+fSize[1]/fNSegment[IZ]*(2.*z+1));
535  G4Transform3D trans;
536  if(fRotationMatrix) {
537  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zpos);
538  trans = G4Translate3D(fCenterPosition)*trans;
539  } else {
540  trans = G4Translate3D(zpos)*G4Translate3D(fCenterPosition);
541  }
542  G4double c[4];
543  colorMap->GetMapColor(rzcell[r][z], c);
544  att.SetColour(c[0], c[1], c[2]);//, c[3]);
545 
546  G4Polyhedron * poly = cylinder.GetPolyhedron();
547  poly->Transform(trans);
548  poly->SetVisAttributes(att);
549  pVisManager->Draw(*poly);
550  }
551  }
552  }
553  pVisManager->EndDraw();
554  }
555 
556  colorMap->SetPSUnit(fDrawUnit);
557  colorMap->SetPSName(fDrawPSName);
558  colorMap->DrawColorChart();
559 
560 }
561 
562 void G4ScoringCylinder::GetRZPhi(G4int index, G4int q[3]) const {
563  // index = k + j * k-size + i * jk-plane-size
564 
565  // nested : z -> phi -> r
566  G4int i = IZ;
567  G4int j = IPHI;
568  G4int k = IR;
569  G4int jk = fNSegment[j]*fNSegment[k];
570  q[i] = index/jk;
571  q[j] = (index - q[i]*jk)/fNSegment[k];
572  q[k] = index - q[j]*fNSegment[k] - q[i]*jk;
573 }