ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CellParameterisation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CellParameterisation.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 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 //
32 // If you use this example, please cite the following publication:
33 // Rad. Prot. Dos. 133 (2009) 2-11
34 
35 #include "CellParameterisation.hh"
36 #include "G4LogicalVolume.hh"
37 #include "G4SystemOfUnits.hh"
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40 
41 // SINGLETON
42 CellParameterisation * CellParameterisation::gInstance = 0;
43 
44 
45 CellParameterisation::CellParameterisation
46 (G4Material * nucleus1, G4Material * cytoplasm1,
47  G4Material * nucleus2, G4Material * cytoplasm2,
48  G4Material * nucleus3, G4Material * cytoplasm3
49  )
50 {
51  fNucleusMaterial1 = nucleus1;
52  fCytoplasmMaterial1 = cytoplasm1;
53  fNucleusMaterial2 = nucleus2;
54  fCytoplasmMaterial2 = cytoplasm2;
55  fNucleusMaterial3 = nucleus3;
56  fCytoplasmMaterial3 = cytoplasm3;
57 
58  G4int ncols,nlines;
59  G4int shiftX, shiftY, shiftZ;
60  G4double x,y,z,mat,den,tmp,density;
61  G4double denCyto1, denCyto2, denCyto3, denNucl1, denNucl2, denNucl3;
62 
63  density=0;
64  denCyto1=0;
65  denCyto2=0;
66  denCyto3=0;
67  denNucl1=0;
68  denNucl2=0;
69  denNucl3=0;
70 
71  ncols=0; nlines=0;
72 
73  // READ PHANTOM
74 
75  fNucleusMass = 0;
76  fCytoplasmMass = 0;
77 
78  FILE *fMap;
79  fMap = fopen("phantom.dat","r");
80 
81  while (1)
82  {
83  if (nlines == 0)
84  {
85  ncols = fscanf(fMap,"%i %i %i",&fPhantomTotalPixels,&fNucleusTotalPixels,&fCytoplasmTotalPixels);
86  fMapCell = new G4ThreeVector[fPhantomTotalPixels];
87  fMaterial = new G4double[fPhantomTotalPixels];
88  fMass = new G4double[fPhantomTotalPixels];
89  fTissueType = new G4int[fPhantomTotalPixels];
90  }
91 
92  if (nlines == 1)
93  {
94  ncols = fscanf(fMap,"%lf %lf %lf",&fDimCellBoxX,&fDimCellBoxY,&fDimCellBoxZ);
95 
96  fDimCellBoxX=fDimCellBoxX*micrometer;
97  fDimCellBoxY=fDimCellBoxY*micrometer;
98  fDimCellBoxZ=fDimCellBoxZ*micrometer;
99  }
100 
101  if (nlines == 2) ncols = fscanf(fMap,"%i %i %i",&shiftX,&shiftY,&shiftZ); // VOXEL SHIFT IN Z ASSUMED TO BE NEGATIVE
102 
103  if (nlines == 3) ncols = fscanf(fMap,"%lf %lf %lf",&denCyto1, &denCyto2, &denCyto3);
104 
105  if (nlines == 4) ncols = fscanf(fMap,"%lf %lf %lf",&denNucl1, &denNucl2, &denNucl3);
106 
107  if (nlines > 4) ncols = fscanf(fMap,"%lf %lf %lf %lf %lf %lf",&x,&y,&z,&mat,&den,&tmp);
108 
109  if (ncols < 0) break;
110 
111  G4ThreeVector v(x+shiftX,y+shiftY,z-1500/(fDimCellBoxZ/micrometer)-shiftZ); // VOXEL SHIFT IN ORDER TO CENTER PHANTOM
112 
113  if (nlines>4)
114  {
115 
116  fMapCell[nlines-5]=v;
117  fMaterial[nlines-5]=mat;
118  fMass[nlines-5]=den;
119 
120  // fTissueType: 1 is Cytoplasm - 2 is Nucleus
121 
122  if( fMaterial[nlines-5] == 2 ) // fMaterial 2 is nucleus
123  {
124  if( fMass[nlines-5] == 1 )
125  {
126  fTissueType[nlines-5]=2;
127  }
128  if( fMass[nlines-5] == 2 )
129  {
130  fTissueType[nlines-5]=2;
131  }
132  if( fMass[nlines-5] == 3 )
133  {
134  fTissueType[nlines-5]=2;
135  }
136  }
137 
138  else if( fMaterial[nlines-5] == 1 ) // fMaterial 1 is cytoplasm
139  {
140  if( fMass[nlines-5] == 1 )
141  {
142  fTissueType[nlines-5]=1;
143  }
144  if( fMass[nlines-5] == 2 )
145  {
146  fTissueType[nlines-5]=2;
147  }
148  if( fMass[nlines-5] == 3 )
149  {
150  fTissueType[nlines-5]=1;
151  }
152  }
153 
154  //
155 
156  if (std::abs(mat-2)<1.e-30) // NUCLEUS
157  {
158  if (std::abs(den-1)<1.e-30) density = denNucl1*(g/cm3);
159  if (std::abs(den-2)<1.e-30) density = denNucl2*(g/cm3);
160  if (std::abs(den-3)<1.e-30) density = denNucl3*(g/cm3);
161  fNucleusMass = fNucleusMass + density * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ ;
162  }
163 
164  if (std::abs(mat-1)<1.e-30) // CYTOPLASM
165  {
166  if (std::abs(den-1)<1e-30) density = denCyto1*(g/cm3);
167  if (std::abs(den-2)<1e-30) density = denCyto2*(g/cm3);
168  if (std::abs(den-3)<1e-30) density = denCyto3*(g/cm3);
169  fCytoplasmMass = fCytoplasmMass + density * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ ;
170  }
171 
172  }
173 
174  nlines++;
175  }
176  fclose(fMap);
177 
178  // NUCLEUS IN GREEN
179 
180  fNucleusAttributes1 = new G4VisAttributes;
181  fNucleusAttributes1->SetColour(G4Colour(0,.8,0));
182  fNucleusAttributes1->SetForceSolid(false);
183 
184  fNucleusAttributes2 = new G4VisAttributes;
185  fNucleusAttributes2->SetColour(G4Colour(0,.9,0));
186  fNucleusAttributes2->SetForceSolid(false);
187 
188  fNucleusAttributes3 = new G4VisAttributes;
189  fNucleusAttributes3->SetColour(G4Colour(0,1,0));
190  fNucleusAttributes3->SetForceSolid(false);
191 
192  // CYTOPLASM IN RED
193 
194  fCytoplasmAttributes1 = new G4VisAttributes;
195  fCytoplasmAttributes1->SetColour(G4Colour(1,0,0));
196  fCytoplasmAttributes1->SetForceSolid(false);
197 
198  fCytoplasmAttributes2 = new G4VisAttributes; // nucleoli in yellow
199  fCytoplasmAttributes2->SetColour(G4Colour(1.,1.,0));
200  fCytoplasmAttributes2->SetForceSolid(false);
201 
202  fCytoplasmAttributes3 = new G4VisAttributes;
203  fCytoplasmAttributes3->SetColour(G4Colour(1,0,0));
204  fCytoplasmAttributes3->SetForceSolid(false);
205 
206  //
207 
208  gInstance = this;
209 
210  }
211 
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
213 
214 CellParameterisation::~CellParameterisation()
215 {
216  delete[] fMapCell;
217  delete[] fMaterial;
218  delete[] fMass;
219  delete[] fTissueType;
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
223 
224 void CellParameterisation::ComputeTransformation
225 (const G4int copyNo, G4VPhysicalVolume* physVol) const
226 {
227  G4ThreeVector origin
228  (
229  fMapCell[copyNo].x()*fDimCellBoxX,
230  fMapCell[copyNo].y()*fDimCellBoxY,
231  fMapCell[copyNo].z()*fDimCellBoxZ
232  );
233 
234  physVol->SetTranslation(origin);
235 }
236 
237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
238 
239 void CellParameterisation::ComputeDimensions
240 (G4Box&, const G4int, const G4VPhysicalVolume*) const
241 {}
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244 
245 G4Material*
246 CellParameterisation::ComputeMaterial(const G4int copyNo,
247  G4VPhysicalVolume* physVol,
248  const G4VTouchable*)
249 {
250  if( fMaterial[copyNo] == 2 ) // fMaterial 2 is nucleus
251  {
252  if( fMass[copyNo] == 1 )
253  {
254  physVol->GetLogicalVolume()->SetVisAttributes( fNucleusAttributes1 );
255  return fNucleusMaterial1;
256  }
257  if( fMass[copyNo] == 2 )
258  {
259  physVol->GetLogicalVolume()->SetVisAttributes( fNucleusAttributes2 );
260  return fNucleusMaterial2;
261  }
262  if( fMass[copyNo] == 3 )
263  {
264  physVol->GetLogicalVolume()->SetVisAttributes( fNucleusAttributes3 );
265  return fNucleusMaterial3;
266  }
267  }
268 
269  else if( fMaterial[copyNo] == 1 ) // fMaterial 1 is cytoplasm
270  {
271  if( fMass[copyNo] == 1 )
272  {
273  physVol->GetLogicalVolume()->SetVisAttributes( fCytoplasmAttributes1 );
274  return fCytoplasmMaterial1;
275  }
276  if( fMass[copyNo] == 2 )
277  {
278  // nucleoli so taken as nucleus !
279  physVol->GetLogicalVolume()->SetVisAttributes( fCytoplasmAttributes2 );
280  return fCytoplasmMaterial2;
281  }
282  if( fMass[copyNo] == 3 )
283  {
284  physVol->GetLogicalVolume()->SetVisAttributes( fCytoplasmAttributes3 );
285  return fCytoplasmMaterial3;
286  }
287  }
288 
289  return physVol->GetLogicalVolume()->GetMaterial();
290 }