cellularPhantom/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 // -----------------------------------------------------------------------------
0027 //       MONTE CARLO SIMULATION OF REALISTIC GEOMETRY FROM MICROSCOPES IMAGES
0028 //
0029 // Authors and contributors:
0030 // P. Barberet (a), S. Incerti (a), N. H. Tran (a), L. Morelli (a,b)
0031 //
0032 // a) University of Bordeaux, CNRS, LP2i, UMR5797, Gradignan, France
0033 // b) Politecnico di Milano, Italy
0034 //
0035 // If you use this code, please cite the following publication:
0036 // P. Barberet et al.,
0037 // "Monte-Carlo dosimetry on a realistic cell monolayer
0038 // geometry exposed to alpha particles."
0039 // Ph. Barberet et al 2012 Phys. Med. Biol. 57 2189
0040 // doi: 110.1088/0031-9155/57/8/2189
0041 // -----------------------------------------------------------------------------
0042
0043 #include "CellParameterisation.hh"
0044
0045 #include "G4Material.hh"
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 CellParameterisation *CellParameterisation::gInstance = nullptr;
0050
0051 CellParameterisation::CellParameterisation
0052 (G4String fileName,
0053  G4Material *RedMat, G4Material *GreenMat, G4Material *BlueMat,
0054  G4double shiftX, G4double shiftY, G4double shiftZ
0055 )
0056 :fRedMaterial(RedMat), fGreenMaterial(GreenMat), fBlueMaterial(BlueMat),
0057  fShiftX(shiftX), fShiftY(shiftY), fShiftZ(shiftZ)
0058 {
0059   Initialize(fileName);
0060 }
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0063
0064 void CellParameterisation::Initialize(const G4String &fileName)
0065 {
0066   G4int ncols, l, mat;
0067   G4int pixelX, pixelY, pixelZ;
0068   G4double x, y, z, den1, den2, den3;
0069
0070   ncols = 0;
0071   l = 0;
0072
0073   // Read phantom
0074
0075   FILE *fMap;
0076   fMap = fopen(fileName, "r");
0077
0078   fRedMass = 0;
0079   fGreenMass = 0;
0080   fBlueMass = 0;
0081
0082   ncols = fscanf(fMap, "%d  %d  %d  %d", &fPhantomTotalPixels, &fRedTotalPixels, &fGreenTotalPixels,
0083                                          &fBlueTotalPixels);
0084   ncols = fscanf(fMap, "%lf  %lf  %lf  %s", &fSizeRealX, &fSizeRealY, &fSizeRealZ, &fRealUnit);
0085   ncols = fscanf(fMap, "%lf  %lf  %lf  %s", &fDimCellBoxX, &fDimCellBoxY, &fDimCellBoxZ, &fRealUnit);
0086
0087   fMapCell = new G4ThreeVector[fPhantomTotalPixels]; //geant4 coordinates space
0088   fMapCellPxl = new G4ThreeVector[fPhantomTotalPixels]; //voxel space
0089   fMapCellOriginal = new G4ThreeVector[fPhantomTotalPixels]; //original coordinates space
0090   fMaterial = new G4int[fPhantomTotalPixels];
0091
0092   fDimCellBoxX = fDimCellBoxX * um;
0093   fDimCellBoxY = fDimCellBoxY * um;
0094   fDimCellBoxZ = fDimCellBoxZ * um;
0095
0096   den1 = fRedMaterial->GetDensity();
0097   den2 = fGreenMaterial->GetDensity();
0098   den3 = fBlueMaterial->GetDensity();
0099
0100   fOffsetX = -fSizeRealX / 2 *um;
0101   fOffsetY = -fSizeRealY / 2 *um;
0102   fOffsetZ = -fSizeRealZ / 2 *um;
0103
0104   G4cout << G4endl;
0105   G4cout << " #########################################################################" << G4endl;
0106   G4cout << "                               Phantom placement and density              " << G4endl;
0107   G4cout << " #########################################################################" << G4endl;
0108   G4cout << G4endl;
0109   G4cout << " ==========> Phantom origin - X (um) = " << (fOffsetX + fShiftX)/um << G4endl;
0110   G4cout << " ==========> Phantom origin - Y (um) = " << (fOffsetY + fShiftY)/um << G4endl;
0111   G4cout << " ==========> Phantom origin - Z (um) = " << (fOffsetZ + fShiftZ)/um << G4endl;
0112   G4cout << G4endl;
0113   G4cout << " ==========> Red density (g/cm3)   = " << den1/(g/cm3) << G4endl;
0114   G4cout << " ==========> Green density (g/cm3) = " << den2/(g/cm3) << G4endl;
0115   G4cout << " ==========> Blue density (g/cm3)  = " << den3/(g/cm3) << G4endl;
0116   G4cout << G4endl;
0117   G4cout << " #########################################################################" << G4endl;
0118   G4cout << G4endl;
0119
0120   while (1)
0121   {
0122     ncols = fscanf(fMap, "%lf %lf %lf %d", &x, &y, &z, &mat);
0123     if (ncols < 0) break;
0124
0125     G4ThreeVector v( x*um + fOffsetX + fShiftX, // phantom shift
0126                    -(y*um + fOffsetY + fShiftY),
0127                      z*um + fOffsetZ + fShiftZ );
0128
0129     // Pixel coordinates
0130     pixelX = (x*um)/fDimCellBoxX;
0131     pixelY = (y*um)/fDimCellBoxY;
0132     pixelZ = (z*um)/fDimCellBoxZ;
0133
0134     G4ThreeVector w(pixelX, pixelY, pixelZ);
0135
0136     G4ThreeVector v_original(x*um, y*um, z*um);
0137
0138     fMapCell[l] = v;
0139     fMapCellPxl[l] = w;
0140     fMapCellOriginal[l] = v_original;
0141
0142     fMaterial[l] = mat;
0143
0144     if (mat == 1){
0145       fRedMass += den1 * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ;
0146     }
0147     else if (mat == 2){
0148       fGreenMass += den2 * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ;
0149     }
0150     else if (mat == 3){
0151       fBlueMass += den3 * fDimCellBoxX * fDimCellBoxY * fDimCellBoxZ;
0152     }
0153     l++;
0154   }
0155
0156   fclose(fMap);
0157
0158   fRedAttributes = new G4VisAttributes;
0159   fRedAttributes->SetColour(G4Colour(1, 0, 0));
0160   fRedAttributes->SetForceSolid(false);
0161
0162   fGreenAttributes = new G4VisAttributes;
0163   fGreenAttributes->SetColour(G4Colour(0, 1, 0));
0164   fGreenAttributes->SetForceSolid(false);
0165
0166   fBlueAttributes = new G4VisAttributes;
0167   fBlueAttributes->SetColour(G4Colour(0, 0, 1));
0168   fBlueAttributes->SetForceSolid(false);
0169
0170   gInstance = this;
0171 }
0172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0173
0174 CellParameterisation::~CellParameterisation()
0175 {
0176   delete[] fMapCell;
0177   delete[] fMapCellPxl;
0178   delete[] fMaterial;
0179 }
0180
0181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0182
0183 void CellParameterisation::ComputeTransformation
0184 (const G4int copyNo, G4VPhysicalVolume *physVol) const
0185 {
0186   if(fMapCell == nullptr)
0187   {
0188     G4ExceptionDescription ex;
0189     ex<< "fMapCell == nullptr ";
0190     G4Exception("CellParameterisation::ComputeTransformation",
0191                 "CellParameterisation001",
0192                 FatalException,
0193                 ex);
0194   }
0195   else
0196   {
0197     G4ThreeVector
0198       origin(fMapCell[copyNo].x(), fMapCell[copyNo].y(), fMapCell[copyNo].z());
0199
0200     physVol->SetTranslation(origin);
0201   }
0202 }
0203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0204
0205 G4Material *
0206 CellParameterisation::ComputeMaterial(const G4int copyNo,
0207                                       G4VPhysicalVolume *physVol,
0208                                       const G4VTouchable *)
0209 {
0210   if (fMaterial[copyNo] == 3) // fMaterial 3 is blue
0211   {
0212     physVol->SetName("physicalMat3");
0213     physVol->GetLogicalVolume()->SetVisAttributes(fBlueAttributes);
0214     return fBlueMaterial;
0215   }
0216   else if (fMaterial[copyNo] == 2) // fMaterial 2 is green
0217   {
0218     physVol->SetName("physicalMat2");
0219     physVol->GetLogicalVolume()->SetVisAttributes(fGreenAttributes);
0220     return fGreenMaterial;
0221   }
0222   else if (fMaterial[copyNo] == 1) // fMaterial 1 is red
0223   {
0224     physVol->SetName("physicalMat1");
0225     physVol->GetLogicalVolume()->SetVisAttributes(fRedAttributes);
0226     return fRedMaterial;
0227   }
0228
0229   return physVol->GetLogicalVolume()->GetMaterial();
0230 }