microbeam/src/CellParameterisation.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software
0028 // shall cite the following Geant4-DNA collaboration publication:
0029 // Med. Phys. 37 (2010) 4692-4708
0030 // The Geant4-DNA web site is available at http://geant4-dna.org
0031 //
0032 // If you use this example, please cite the following publication:
0033 // Rad. Prot. Dos. 133 (2009) 2-11
0034
0035 #include "CellParameterisation.hh"
0036 #include "G4LogicalVolume.hh"
0037 #include "G4SystemOfUnits.hh"
0038
0039 // SINGLETON
0040 CellParameterisation * CellParameterisation::gInstance = 0;
0041
0042 CellParameterisation::CellParameterisation
0043 (G4Material * nucleus1,  G4Material * cytoplasm1,
0044  G4Material * nucleus2,  G4Material * cytoplasm2,
0045  G4Material * nucleus3,  G4Material * cytoplasm3
0046  )
0047 {
0048    fNucleusMaterial1   = nucleus1;
0049    fCytoplasmMaterial1 = cytoplasm1;
0050    fNucleusMaterial2   = nucleus2;
0051    fCytoplasmMaterial2 = cytoplasm2;
0052    fNucleusMaterial3   = nucleus3;
0053    fCytoplasmMaterial3 = cytoplasm3;
0054
0055    G4int ncols,nlines;
0056    G4int shiftX, shiftY, shiftZ;
0057    G4double x,y,z,mat,den,tmp,density;
0058    G4double denCyto1, denCyto2, denCyto3, denNucl1, denNucl2, denNucl3;
0059
0060    ncols = nlines = shiftX = shiftY = shiftZ = 0;
0061    x = y = z = mat = den = tmp = density =
0062      denCyto1 = denCyto2 = denCyto3 = denNucl1 = denNucl2 = denNucl3 = 0.0;
0063
0064    // READ PHANTOM
0065    fNucleusMass = 0;
0066    fCytoplasmMass = 0;
0067
0068    fDimCellBoxX = fDimCellBoxY = fDimCellBoxZ = micrometer;
0069
0070    FILE *fMap;
0071    fMap = fopen("phantom.dat","r");
0072
0073    while (1)
0074    {
0075       if (nlines == 0)
0076       {
0077         ncols = fscanf(fMap,"%i %i %i",&fPhantomTotalPixels,&fNucleusTotalPixels,&fCytoplasmTotalPixels);
0078     fMapCell    = new G4ThreeVector[fPhantomTotalPixels];
0079         fMaterial   = new G4double[fPhantomTotalPixels];
0080         fMass       = new G4double[fPhantomTotalPixels];
0081         fTissueType = new G4int[fPhantomTotalPixels];
0082       }
0083
0084       if (nlines == 1)
0085       {
0086         ncols = fscanf(fMap,"%lf %lf %lf",&fDimCellBoxX,&fDimCellBoxY,&fDimCellBoxZ);
0087     fDimCellBoxX=fDimCellBoxX*micrometer;
0088         fDimCellBoxY=fDimCellBoxY*micrometer;
0089         fDimCellBoxZ=fDimCellBoxZ*micrometer;
0090       }
0091
0092       // VOXEL SHIFT IN Z ASSUMED TO BE NEGATIVE
0093       if (nlines == 2) ncols = fscanf(fMap,"%i %i %i",&shiftX,&shiftY,&shiftZ);
0094
0095       if (nlines == 3) ncols = fscanf(fMap,"%lf %lf %lf",&denCyto1, &denCyto2, &denCyto3);
0096
0097       if (nlines == 4) ncols = fscanf(fMap,"%lf %lf %lf",&denNucl1, &denNucl2, &denNucl3);
0098
0099       if (nlines >  4) ncols = fscanf(fMap,"%lf %lf %lf %lf %lf %lf",&x,&y,&z,&mat,&den,&tmp);
0100
0101       if (ncols  <  0) break;
0102
0103       // VOXEL SHIFT IN ORDER TO CENTER PHANTOM
0104       G4ThreeVector v(x+shiftX,y+shiftY,z-1500/(fDimCellBoxZ/micrometer)-shiftZ);
0105       if (nlines>4)
0106       {
0107
0108       fMapCell[nlines-5]=v;
0109       fMaterial[nlines-5]=mat;
0110       fMass[nlines-5]=den;
0111
0112       // fTissueType: 1 is Cytoplasm - 2 is Nucleus
0113
0114           if( fMaterial[nlines-5] == 2 ) // fMaterial 2 is nucleus
0115           {
0116         if( fMass[nlines-5] == 1 )
0117         {
0118           fTissueType[nlines-5]=2;
0119         }
0120         if( fMass[nlines-5] == 2 )
0121         {
0122           fTissueType[nlines-5]=2;
0123         }
0124         if( fMass[nlines-5] == 3 )
0125         {
0126           fTissueType[nlines-5]=2;
0127         }
0128       }
0129
0130           else if( fMaterial[nlines-5] == 1 ) // fMaterial 1 is cytoplasm
0131       {
0132         if( fMass[nlines-5] == 1 )
0133         {
0134           fTissueType[nlines-5]=1;
0135         }
0136         if( fMass[nlines-5] == 2 )
0137         {
0138           fTissueType[nlines-5]=2;
0139         }
0140         if( fMass[nlines-5] == 3 )
0141         {
0142           fTissueType[nlines-5]=1;
0143         }
0144       }
0145
0146       //
0147
0148       if (std::abs(mat-2)<1.e-30) // NUCLEUS
0149           {
0150         if (std::abs(den-1)<1.e-30) density = denNucl1*(g/cm3);
0151         if (std::abs(den-2)<1.e-30) density = denNucl2*(g/cm3);
0152         if (std::abs(den-3)<1.e-30) density = denNucl3*(g/cm3);
0153         fNucleusMass   = fNucleusMass   + density * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ ;
0154           }
0155
0156           if (std::abs(mat-1)<1.e-30) // CYTOPLASM
0157           {
0158         if (std::abs(den-1)<1e-30) density = denCyto1*(g/cm3);
0159         if (std::abs(den-2)<1e-30) density = denCyto2*(g/cm3);
0160         if (std::abs(den-3)<1e-30) density = denCyto3*(g/cm3);
0161         fCytoplasmMass = fCytoplasmMass + density * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ ;
0162       }
0163
0164       }
0165
0166       nlines++;
0167    }
0168    fclose(fMap);
0169
0170   // NUCLEUS IN GREEN
0171
0172   fNucleusAttributes1 = new G4VisAttributes;
0173   fNucleusAttributes1->SetColour(G4Colour(0,.8,0));
0174   fNucleusAttributes1->SetForceSolid(false);
0175
0176   fNucleusAttributes2 = new G4VisAttributes;
0177   fNucleusAttributes2->SetColour(G4Colour(0,.9,0));
0178   fNucleusAttributes2->SetForceSolid(false);
0179
0180   fNucleusAttributes3 = new G4VisAttributes;
0181   fNucleusAttributes3->SetColour(G4Colour(0,1,0));
0182   fNucleusAttributes3->SetForceSolid(false);
0183
0184   // CYTOPLASM IN RED
0185
0186   fCytoplasmAttributes1 = new G4VisAttributes;
0187   fCytoplasmAttributes1->SetColour(G4Colour(1,0,0));
0188   fCytoplasmAttributes1->SetForceSolid(false);
0189
0190   fCytoplasmAttributes2 = new G4VisAttributes; // nucleoli in yellow
0191   fCytoplasmAttributes2->SetColour(G4Colour(1.,1.,0));
0192   fCytoplasmAttributes2->SetForceSolid(false);
0193
0194   fCytoplasmAttributes3 = new G4VisAttributes;
0195   fCytoplasmAttributes3->SetColour(G4Colour(1,0,0));
0196   fCytoplasmAttributes3->SetForceSolid(false);
0197
0198   //
0199   gInstance = this;
0200  }
0201
0202 CellParameterisation::~CellParameterisation()
0203 {
0204   delete[] fMapCell;
0205   delete[] fMaterial;
0206   delete[] fMass;
0207   delete[] fTissueType;
0208 }
0209
0210 void CellParameterisation::ComputeTransformation
0211 (const G4int copyNo, G4VPhysicalVolume* physVol) const
0212 {
0213   G4ThreeVector origin
0214   (
0215       fMapCell[copyNo].x()*fDimCellBoxX,
0216       fMapCell[copyNo].y()*fDimCellBoxY,
0217       fMapCell[copyNo].z()*fDimCellBoxZ
0218   );
0219
0220   physVol->SetTranslation(origin);
0221 }
0222
0223 void CellParameterisation::ComputeDimensions
0224 (G4Box&, const G4int, const G4VPhysicalVolume*) const
0225 {}
0226
0227 G4Material*
0228 CellParameterisation::ComputeMaterial(const G4int copyNo,
0229                                       G4VPhysicalVolume* physVol,
0230                                       const G4VTouchable*)
0231 {
0232     if( fMaterial[copyNo] == 2 ) // fMaterial 2 is nucleus
0233         {
0234      if( fMass[copyNo] == 1 )
0235         {
0236           physVol->GetLogicalVolume()->SetVisAttributes( fNucleusAttributes1 );
0237           return fNucleusMaterial1;
0238         }
0239      if( fMass[copyNo] == 2 )
0240         {
0241           physVol->GetLogicalVolume()->SetVisAttributes( fNucleusAttributes2 );
0242           return fNucleusMaterial2;
0243         }
0244      if( fMass[copyNo] == 3 )
0245         {
0246           physVol->GetLogicalVolume()->SetVisAttributes( fNucleusAttributes3 );
0247           return fNucleusMaterial3;
0248         }
0249     }
0250
0251     else if( fMaterial[copyNo] == 1 ) // fMaterial 1 is cytoplasm
0252     {
0253      if( fMass[copyNo] == 1 )
0254         {
0255           physVol->GetLogicalVolume()->SetVisAttributes( fCytoplasmAttributes1 );
0256           return fCytoplasmMaterial1;
0257         }
0258      if( fMass[copyNo] == 2 )
0259         {
0260           // nucleoli so taken as nucleus !
0261           physVol->GetLogicalVolume()->SetVisAttributes( fCytoplasmAttributes2 );
0262           return fCytoplasmMaterial2;
0263         }
0264      if( fMass[copyNo] == 3 )
0265         {
0266           physVol->GetLogicalVolume()->SetVisAttributes( fCytoplasmAttributes3 );
0267           return fCytoplasmMaterial3;
0268         }
0269     }
0270
0271     return physVol->GetLogicalVolume()->GetMaterial();
0272 }