Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ScoringCylinder.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id$
28 //
29 
30 #include "G4ScoringCylinder.hh"
31 
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4VPhysicalVolume.hh"
35 #include "G4Tubs.hh"
36 #include "G4LogicalVolume.hh"
37 #include "G4VPhysicalVolume.hh"
38 #include "G4PVPlacement.hh"
39 #include "G4PVReplica.hh"
40 #include "G4PVDivision.hh"
41 #include "G4VisAttributes.hh"
42 #include "G4VVisManager.hh"
43 #include "G4VScoreColorMap.hh"
44 
45 #include "G4SDManager.hh"
47 #include "G4SDParticleFilter.hh"
48 #include "G4VPrimitiveScorer.hh"
49 #include "G4PSEnergyDeposit.hh"
50 #include "G4PSTrackLength.hh"
51 #include "G4PSNofStep.hh"
52 #include "G4ScoringManager.hh"
53 
54 
56  :G4VScoringMesh(wName), fMeshElementLogical(0)
57 {
59 
60  fDivisionAxisNames[0] = "Z";
61  fDivisionAxisNames[1] = "PHI";
62  fDivisionAxisNames[2] = "R";
63 }
64 
66 {;}
67 
69 {
70  if(fConstructed) {
71 
72  if(verboseLevel > 0)
73  G4cout << fWorldPhys->GetName() << " --- All quantities are reset." << G4endl;
74  ResetScore();
75 
76  } else {
77  fConstructed = true;
78  SetupGeometry(fWorldPhys);
79  }
80 }
81 
82 
83 
84 void G4ScoringCylinder::SetupGeometry(G4VPhysicalVolume * fWorldPhys) {
85 
86  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::SetupGeometry() ..." << G4endl;
87 
88  // World
89  G4VPhysicalVolume * scoringWorld = fWorldPhys;
90  G4LogicalVolume * worldLogical = scoringWorld->GetLogicalVolume();
91 
92  // Scoring Mesh
93  if(verboseLevel > 9) G4cout << fWorldName << G4endl;
94  G4String tubsName = fWorldName+"_mesh";
95 
96  if(verboseLevel > 9) G4cout << "R max., Dz =: " << fSize[0] << ", " << fSize[1] << G4endl;
97  G4VSolid * tubsSolid = new G4Tubs(tubsName+"0", // name
98  0., // R min
99  fSize[0], // R max
100  fSize[1], // Dz
101  0., // starting phi
102  twopi*rad); // segment phi
103  G4LogicalVolume * tubsLogical = new G4LogicalVolume(tubsSolid, 0, tubsName);
105  tubsLogical, tubsName+"0", worldLogical, false, 0);
106 
107  if(verboseLevel > 9) G4cout << " # of segments : r, phi, z =: "
108  << fNSegment[IR] << ", " << fNSegment[IPHI] << ", " << fNSegment[IZ] << G4endl;
109 
110  G4String layerName[2] = {tubsName + "1", tubsName + "2"};
111  G4VSolid * layerSolid[2];
112  G4LogicalVolume * layerLogical[2];
113 
114  //-- fisrt nested layer (replicated along z direction)
115  if(verboseLevel > 9) G4cout << "layer 1 :" << G4endl;
116  layerSolid[0] = new G4Tubs(layerName[0], // name
117  0., // inner radius
118  fSize[0], // outer radius
119  fSize[1]/fNSegment[IZ], // half len. in z
120  0., // starting phi angle
121  twopi*rad); // delta angle of the segment
122  layerLogical[0] = new G4LogicalVolume(layerSolid[0], 0, layerName[0]);
123  if(fNSegment[IZ] > 1) {
124  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replicate along z direction" << G4endl;
126  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replica" << G4endl;
127  new G4PVReplica(layerName[0], layerLogical[0], tubsLogical, kZAxis, fNSegment[IZ], 2.*fSize[1]/fNSegment[IZ]);
128  } else {
129  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Division" << G4endl;
130  new G4PVDivision(layerName[0], layerLogical[0], tubsLogical, kZAxis, fNSegment[IZ], 0.);
131  }
132  } else if(fNSegment[IZ] == 1) {
133  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Placement" << G4endl;
134  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), layerLogical[0], layerName[0], tubsLogical, false, 0);
135  } else {
136  G4cerr << "G4ScoringCylinder::SetupGeometry() : invalid parameter ("
137  << fNSegment[IZ] << ") "
138  << "in placement of the first nested layer." << G4endl;
139  }
140 
141  // second nested layer (replicated along phi direction)
142  if(verboseLevel > 9) G4cout << "layer 2 :" << G4endl;
143  layerSolid[1] = new G4Tubs(layerName[1],
144  0.,
145  fSize[0],
146  fSize[1]/fNSegment[IZ],
147  0.,
148  twopi*rad/fNSegment[IPHI]);
149  layerLogical[1] = new G4LogicalVolume(layerSolid[1], 0, layerName[1]);
150  if(fNSegment[IPHI] > 1) {
151  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replicate along phi direction" << G4endl;
153  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replica" << G4endl;
154  new G4PVReplica(layerName[1], layerLogical[1], layerLogical[0], kPhi,
155  fNSegment[IPHI], twopi*rad/fNSegment[IPHI]);
156  } else {
157  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Division" << G4endl;
158  new G4PVDivision(layerName[1], layerLogical[1], layerLogical[0], kPhi, fNSegment[IPHI], 0.);
159  }
160  } else if(fNSegment[IPHI] == 1) {
161  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Placement" << G4endl;
162  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), layerLogical[1], layerName[1], layerLogical[0], false, 0);
163  } else
164  G4cerr << "ERROR : G4ScoringCylinder::SetupGeometry() : invalid parameter ("
165  << fNSegment[IPHI] << ") "
166  << "in placement of the second nested layer." << G4endl;
167 
168  // mesh elements
169  if(verboseLevel > 9) G4cout << "mesh elements :" << G4endl;
170  G4String elementName = tubsName +"3";
171  G4VSolid * elementSolid = new G4Tubs(elementName,
172  0.,
173  fSize[0]/fNSegment[IR],
174  fSize[1]/fNSegment[IZ],
175  0.,
176  twopi*rad/fNSegment[IPHI]);
177  fMeshElementLogical = new G4LogicalVolume(elementSolid, 0, elementName);
178  if(fNSegment[IR] > 1) {
179 
180  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replicate along r direction" << G4endl;
181 
183  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Replica" << G4endl;
184  new G4PVReplica(elementName, fMeshElementLogical, layerLogical[1], kRho,
185  fNSegment[IR], fSize[0]/fNSegment[IR]);
186  } else {
187  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Division" << G4endl;
188  new G4PVDivision(elementName, fMeshElementLogical, layerLogical[1], kRho, fNSegment[IR], 0.);
189  }
190  } else if(fNSegment[IR] == 1) {
191  if(verboseLevel > 9) G4cout << "G4ScoringCylinder::Construct() : Placement" << G4endl;
192  new G4PVPlacement(0, G4ThreeVector(0.,0.,0.), fMeshElementLogical, elementName, layerLogical[1], false, 0);
193  } else {
194  G4cerr << "G4ScoringCylinder::SetupGeometry() : "
195  << "invalid parameter (" << fNSegment[IR] << ") "
196  << "in mesh element placement." << G4endl;
197  }
198 
199  // set the sensitive detector
200  fMeshElementLogical->SetSensitiveDetector(fMFD);
201 
202 
203  // vis. attributes
204  G4VisAttributes * visatt = new G4VisAttributes(G4Colour(.5,.5,.5));
205  visatt->SetVisibility(true);
206  layerLogical[0]->SetVisAttributes(visatt);
207  layerLogical[1]->SetVisAttributes(visatt);
208  visatt = new G4VisAttributes(G4Colour(.5,.5,.5,0.01));
209  //visatt->SetForceSolid(true);
210  fMeshElementLogical->SetVisAttributes(visatt);
211 }
212 
214  G4cout << "G4ScoringCylinder : " << fWorldName << " --- Shape: Cylindrical mesh" << G4endl;
215 
216  G4cout << " Size (R, Dz): ("
217  << fSize[0]/cm << ", "
218  << fSize[1]/cm << ") [cm]"
219  << G4endl;
220 
222 }
223 
224 
225 void G4ScoringCylinder::Draw(std::map<G4int, G4double*> * map, G4VScoreColorMap* colorMap, G4int axflg) {
226 
228  if(pVisManager) {
229 
230  // cell vectors
231  std::vector<double> ephi;
232  for(int phi = 0; phi < fNSegment[IPHI]; phi++) ephi.push_back(0.);
233  //-
234  std::vector<std::vector<double> > zphicell; // zphicell[Z][PHI]
235  for(int z = 0; z < fNSegment[IZ]; z++) zphicell.push_back(ephi);
236  //-
237  std::vector<std::vector<double> > rphicell; // rphicell[R][PHI]
238  for(int r = 0; r < fNSegment[IR]; r++) rphicell.push_back(ephi);
239 
240  // projections
241  G4int q[3];
242  std::map<G4int, G4double*>::iterator itr = map->begin();
243  for(; itr != map->end(); itr++) {
244  if(itr->first < 0) {
245  G4cout << itr->first << G4endl;
246  continue;
247  }
248  GetRZPhi(itr->first, q);
249 
250  zphicell[q[IZ]][q[IPHI]] += *(itr->second)/fDrawUnitValue;
251  rphicell[q[IR]][q[IPHI]] += *(itr->second)/fDrawUnitValue;
252  }
253 
254  // search min./max. values
255  G4double zphimin = DBL_MAX, rphimin = DBL_MAX;
256  G4double zphimax = 0., rphimax = 0.;
257  for(int iphi = 0; iphi < fNSegment[IPHI]; iphi++) {
258  for(int iz = 0; iz < fNSegment[IZ]; iz++) {
259  if(zphimin > zphicell[iz][iphi]) zphimin = zphicell[iz][iphi];
260  if(zphimax < zphicell[iz][iphi]) zphimax = zphicell[iz][iphi];
261  }
262  for(int ir = 0; ir < fNSegment[IR]; ir++) {
263  if(rphimin > rphicell[ir][iphi]) rphimin = rphicell[ir][iphi];
264  if(rphimax < rphicell[ir][iphi]) rphimax = rphicell[ir][iphi];
265  }
266  }
267 
268  G4VisAttributes att;
269  att.SetForceSolid(true);
270  att.SetForceAuxEdgeVisible(true);
271 
272 
274  if(axflg/100==1) {
275  // rz plane
276  }
277  axflg = axflg%100;
278  if(axflg/10==1) {
279  // z-phi plane
280  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(zphimin, zphimax); }
281 
282  G4double zhalf = fSize[1]/fNSegment[IZ];
283  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
284  for(int z = 0; z < fNSegment[IZ]; z++) {
285  //-
286  G4double angle = twopi/fNSegment[IPHI]*phi;
287  G4double dphi = twopi/fNSegment[IPHI];
288  G4Tubs cylinder("z-phi", // name
289  fSize[0]*0.99, fSize[0], // inner radius, outer radius
290  zhalf, // half length in z
291  angle, dphi*0.99999); // starting phi angle, delta angle
292  //-
293  G4ThreeVector zpos(0., 0., -fSize[1] + fSize[1]/fNSegment[IZ]*(1 + 2.*z));
294  G4Transform3D trans;
295  if(fRotationMatrix) {
296  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zpos);
297  trans = G4Translate3D(fCenterPosition)*trans;
298  } else {
300  }
301  G4double c[4];
302  colorMap->GetMapColor(zphicell[z][phi], c);
303  att.SetColour(c[0], c[1], c[2]);//, c[3]);
304  //-
305  pVisManager->Draw(cylinder, att, trans);
306  }
307  }
308  }
309  axflg = axflg%10;
310  if(axflg==1) {
311  // r-phi plane
312  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(rphimin, rphimax); }
313 
314  G4double rsize = fSize[0]/fNSegment[IR];
315  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
316  for(int r = 0; r < fNSegment[IR]; r++) {
317 
318  G4double rs[2] = {rsize*r, rsize*(r+1)};
319  G4double angle = twopi/fNSegment[IPHI]*phi;
320  G4double dphi = twopi/fNSegment[IPHI];
321  G4Tubs cylinder("z-phi", rs[0], rs[1], 0.001,
322  angle, dphi*0.99999);
323  /*
324  G4cout << ">>>> "
325  << rs[0] << " - " << rs[1] << " : "
326  << angle << " - " << angle + dphi
327  << G4endl;
328  */
329 
330  G4ThreeVector zposn(0., 0., -fSize[1]);
331  G4ThreeVector zposp(0., 0., fSize[1]);
332  G4Transform3D transn, transp;
333  if(fRotationMatrix) {
334  transn = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zposn);
335  transn = G4Translate3D(fCenterPosition)*transn;
336  transp = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zposp);
337  transp = G4Translate3D(fCenterPosition)*transp;
338  } else {
341  }
342  G4double c[4];
343  colorMap->GetMapColor(rphicell[r][phi], c);
344  att.SetColour(c[0], c[1], c[2]);//, c[3]);
345  /*
346  G4cout << " " << c[0] << ", "
347  << c[1] << ", " << c[2] << G4endl;
348  */
349  pVisManager->Draw(cylinder, att, transn);
350  pVisManager->Draw(cylinder, att, transp);
351  }
352  }
353  }
354 
355  colorMap->SetPSUnit(fDrawUnit);
356  colorMap->SetPSName(fDrawPSName);
357  colorMap->DrawColorChart();
358 
359  }
360 }
361 
362 void G4ScoringCylinder::DrawColumn(std::map<G4int, G4double*> * map, G4VScoreColorMap* colorMap,
363  G4int idxProj, G4int idxColumn)
364 {
365  G4int projAxis = 0;
366  switch(idxProj) {
367  case 0:
368  projAxis = IR;
369  break;
370  case 1:
371  projAxis = IZ;
372  break;
373  case 2:
374  projAxis = IPHI;
375  break;
376  }
377 
378  if(idxColumn<0 || idxColumn>=fNSegment[projAxis])
379  {
380  G4cerr << "Warning : Column number " << idxColumn << " is out of scoring mesh [0," << fNSegment[projAxis]-1 <<
381  "]. Method ignored." << G4endl;
382  return;
383  }
385  if(pVisManager) {
386 
387  // cell vectors
388  std::vector<std::vector<std::vector<double> > > cell; // cell[R][Z][PHI]
389  std::vector<double> ephi;
390  for(int phi = 0; phi < fNSegment[IPHI]; phi++) ephi.push_back(0.);
391  std::vector<std::vector<double> > ezphi;
392  for(int z = 0; z < fNSegment[IZ]; z++) ezphi.push_back(ephi);
393  for(int r = 0; r < fNSegment[IR]; r++) cell.push_back(ezphi);
394 
395  std::vector<std::vector<double> > rzcell; // rzcell[R][Z]
396  std::vector<double> ez;
397  for(int z = 0; z < fNSegment[IZ]; z++) ez.push_back(0.);
398  for(int r = 0; r < fNSegment[IR]; r++) rzcell.push_back(ez);
399 
400  std::vector<std::vector<double> > zphicell; // zphicell[Z][PHI]
401  for(int z = 0; z < fNSegment[IZ]; z++) zphicell.push_back(ephi);
402 
403  std::vector<std::vector<double> > rphicell; // rphicell[R][PHI]
404  for(int r = 0; r < fNSegment[IR]; r++) rphicell.push_back(ephi);
405 
406  // projections
407  G4int q[3];
408  std::map<G4int, G4double*>::iterator itr = map->begin();
409  for(; itr != map->end(); itr++) {
410  if(itr->first < 0) {
411  G4cout << itr->first << G4endl;
412  continue;
413  }
414  GetRZPhi(itr->first, q);
415 
416  if(projAxis == IR && q[IR] == idxColumn) { // zphi plane
417  zphicell[q[IZ]][q[IPHI]] += *(itr->second)/fDrawUnitValue;
418  }
419  if(projAxis == IZ && q[IZ] == idxColumn) { // rphi plane
420  rphicell[q[IR]][q[IPHI]] += *(itr->second)/fDrawUnitValue;
421  }
422  if(projAxis == IPHI && q[IPHI] == idxColumn) { // rz plane
423  rzcell[q[IR]][q[IZ]] += *(itr->second)/fDrawUnitValue;
424  }
425  }
426 
427  // search min./max. values
428  G4double rzmin = DBL_MAX, zphimin = DBL_MAX, rphimin = DBL_MAX;
429  G4double rzmax = 0., zphimax = 0., rphimax = 0.;
430  for(int r = 0; r < fNSegment[IR]; r++) {
431  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
432  if(rphimin > rphicell[r][phi]) rphimin = rphicell[r][phi];
433  if(rphimax < rphicell[r][phi]) rphimax = rphicell[r][phi];
434  }
435  for(int z = 0; z < fNSegment[IZ]; z++) {
436  if(rzmin > rzcell[r][z]) rzmin = rzcell[r][z];
437  if(rzmax < rzcell[r][z]) rzmax = rzcell[r][z];
438  }
439  }
440  for(int z = 0; z < fNSegment[IZ]; z++) {
441  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
442  if(zphimin > zphicell[z][phi]) zphimin = zphicell[z][phi];
443  if(zphimax < zphicell[z][phi]) zphimax = zphicell[z][phi];
444  }
445  }
446 
447 
448  G4VisAttributes att;
449  att.SetForceSolid(true);
450  att.SetForceAuxEdgeVisible(true);
451 
452 
454  // z-phi plane
455  if(projAxis == IR) {
456  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(zphimin,zphimax); }
457 
458  G4double zhalf = fSize[1]/fNSegment[IZ];
459  G4double rsize[2] = {fSize[0]/fNSegment[IR]*idxColumn,
460  fSize[0]/fNSegment[IR]*(idxColumn+1)};
461  for(int phi = 0; phi < fNSegment[IPHI]; phi++) {
462  for(int z = 0; z < fNSegment[IZ]; z++) {
463 
464  G4double angle = twopi/fNSegment[IPHI]*phi*radian;
465  G4double dphi = twopi/fNSegment[IPHI]*radian;
466  G4Tubs cylinder("z-phi", rsize[0], rsize[1], zhalf,
467  angle, dphi*0.99999);
468 
469  G4ThreeVector zpos(0., 0., -fSize[1] + fSize[1]/fNSegment[IZ]*(1 + 2.*z));
470  G4Transform3D trans;
471  if(fRotationMatrix) {
472  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zpos);
473  trans = G4Translate3D(fCenterPosition)*trans;
474  } else {
476  }
477  G4double c[4];
478  colorMap->GetMapColor(zphicell[z][phi], c);
479  att.SetColour(c[0], c[1], c[2]);//, c[3]);
480  pVisManager->Draw(cylinder, att, trans);
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 {
506  }
507  G4double c[4];
508  colorMap->GetMapColor(rphicell[r][phi], c);
509  att.SetColour(c[0], c[1], c[2]);//, c[3]);
510  pVisManager->Draw(cylinder, att, trans);
511  }
512  }
513 
514  // r-z plane
515  } else if(projAxis == IPHI) {
516  if(colorMap->IfFloatMinMax()) { colorMap->SetMinMax(rzmin,rzmax); }
517 
518  G4double rsize = fSize[0]/fNSegment[IR];
519  G4double zhalf = fSize[1]/fNSegment[IZ];
520  G4double angle = twopi/fNSegment[IPHI]*idxColumn*radian;
521  G4double dphi = twopi/fNSegment[IPHI]*radian;
522  for(int z = 0; z < fNSegment[IZ]; z++) {
523  for(int r = 0; r < fNSegment[IR]; r++) {
524 
525  G4double rs[2] = {rsize*r, rsize*(r+1)};
526  G4Tubs cylinder("z-phi", rs[0], rs[1], zhalf,
527  angle, dphi);
528 
529  G4ThreeVector zpos(0., 0.,
530  -fSize[1]+fSize[1]/fNSegment[IZ]*(2.*z+1));
531  G4Transform3D trans;
532  if(fRotationMatrix) {
533  trans = G4Rotate3D(*fRotationMatrix).inverse()*G4Translate3D(zpos);
534  trans = G4Translate3D(fCenterPosition)*trans;
535  } else {
537  }
538  G4double c[4];
539  colorMap->GetMapColor(rzcell[r][z], c);
540  att.SetColour(c[0], c[1], c[2]);//, c[3]);
541  pVisManager->Draw(cylinder, att, trans);
542  }
543  }
544  }
545  }
546 
547  colorMap->SetPSUnit(fDrawUnit);
548  colorMap->SetPSName(fDrawPSName);
549  colorMap->DrawColorChart();
550 
551 }
552 
554  // index = k + j * k-size + i * jk-plane-size
555 
556  // nested : z -> phi -> r
557  G4int i = IZ;
558  G4int j = IPHI;
559  G4int k = IR;
560  G4int jk = fNSegment[j]*fNSegment[k];
561  q[i] = index/jk;
562  q[j] = (index - q[i]*jk)/fNSegment[k];
563  q[k] = index - q[j]*fNSegment[k] - q[i]*jk;
564 }