brachytherapy/src/BrachyUserScoreWriter.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 /*
0028 // Code developed by:
0029 // S.Guatelli, susanna@uow.edu.au
0030 //
0031 Original code from geant4/examples/extended/runAndEvent/RE03, by M. Asai
0032 */
0033
0034 #include "BrachyUserScoreWriter.hh"
0035 #include "G4AnalysisManager.hh"
0036 #include "G4MultiFunctionalDetector.hh"
0037 #include "G4SDParticleFilter.hh"
0038 #include "G4VPrimitiveScorer.hh"
0039 #include "G4VScoringMesh.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include <map>
0042 #include <fstream>
0043 // The default output is
0044 // voxelX, voxelY, voxelZ, edep
0045 // The BrachyUserScoreWriter allows to change the format of the output file.
0046 // in the specific case:
0047 // xx (mm)  yy(mm) zz(mm) edep(keV)
0048 // The same information is stored in a ntuple, in the
0049 // brachytherapy.root file
0050
0051 BrachyUserScoreWriter::BrachyUserScoreWriter():
0052 G4VScoreWriter()
0053 {
0054 }
0055
0056 BrachyUserScoreWriter::~BrachyUserScoreWriter()
0057 {;}
0058
0059 void BrachyUserScoreWriter::DumpQuantityToFile(const G4String & psName,
0060                                                const G4String & fileName,
0061                                                const G4String & option)
0062 {
0063 using MeshScoreMap = G4VScoringMesh::MeshScoreMap;
0064
0065 if(verboseLevel > 0)
0066   {G4cout << "BrachyUserScorer-defined DumpQuantityToFile() method is invoked."
0067   << G4endl;
0068   }
0069
0070 // change the option string into lowercase to the case-insensitive.
0071 G4String opt = option;
0072 std::transform(opt.begin(), opt.end(), opt.begin(), (int (*)(int))(tolower));
0073
0074 // confirm the option
0075 if(opt.size() == 0) opt = "csv";
0076
0077 // open the file
0078 std::ofstream ofile(fileName);
0079
0080 if(!ofile)
0081 {
0082    G4cerr << "ERROR : DumpToFile : File open error -> " << fileName << G4endl;
0083    return;
0084 }
0085   ofile << "# mesh name: " << fScoringMesh->GetWorldName() << G4endl;
0086
0087 // retrieve the map
0088 MeshScoreMap fSMap = fScoringMesh -> GetScoreMap();
0089
0090 auto msMapItr = fSMap.find(psName);
0091
0092 if(msMapItr == fSMap.end())
0093   {
0094    G4cerr << "ERROR : DumpToFile : Unknown quantity, \""<< psName
0095    << "\"." << G4endl;
0096    return;
0097   }
0098
0099 auto score = msMapItr-> second-> GetMap();
0100
0101 ofile << "# primitive scorer name: " << msMapItr -> first << G4endl;
0102 //
0103 // Write quantity in the ASCII output file and in brachytherapy.root
0104 //
0105 ofile << std::setprecision(16); // for double value with 8 bytes
0106
0107 auto analysisManager = G4AnalysisManager::Instance();
0108
0109 G4bool fileOpen = analysisManager -> OpenFile("brachytherapy.root");
0110  if (! fileOpen) {
0111     G4cerr << "\n---> The ROOT output file has not been opened "
0112            << analysisManager->GetFileName() << G4endl;
0113   }
0114
0115 G4cout << "Using " << analysisManager -> GetType() << G4endl;
0116 analysisManager -> SetVerboseLevel(1);
0117 analysisManager -> SetActivation(true);
0118
0119 // Create histograms
0120 G4int histo2= analysisManager-> CreateH2("h20","edep2Dxy", 801, -100.125, 100.125, 801, -100.125, 100.125);
0121
0122 // Histo 0 with the energy spectrum will not be saved
0123 // in brachytherapy.root
0124 analysisManager->SetH1Activation(0, false);
0125 analysisManager->SetH2Activation(histo2, true);
0126
0127 for(int x = 0; x < fNMeshSegments[0]; x++) {
0128    for(int y = 0; y < fNMeshSegments[1]; y++) {
0129      for(int z = 0; z < fNMeshSegments[2]; z++){
0130         G4int numberOfVoxel_x = fNMeshSegments[0];
0131         G4int numberOfVoxel_y = fNMeshSegments[1];
0132         G4int numberOfVoxel_z =fNMeshSegments[2];
0133         // If the voxel width is changed in the macro file,
0134         // the voxel width variable must be updated
0135         G4double voxelWidth = 0.25 *CLHEP::mm;
0136         //
0137         G4double xx = ( - numberOfVoxel_x + 1+ 2*x )* voxelWidth/2;
0138         G4double yy = ( - numberOfVoxel_y + 1+ 2*y )* voxelWidth/2;
0139         G4double zz = ( - numberOfVoxel_z + 1+ 2*z )* voxelWidth/2;
0140         G4int idx = GetIndex(x, y, z);
0141         std::map<G4int, G4StatDouble*>::iterator value = score -> find(idx);
0142
0143        if (value != score -> end())
0144         {
0145          // Print in the ASCII output file the information
0146
0147          ofile << xx << "  " << yy << "  " << zz <<"  "
0148                <<(value->second->sum_wx())/keV << G4endl;
0149
0150         // Save the same information in the ROOT output file
0151
0152     if(zz> -0.125 *CLHEP::mm && zz < 0.125/mm)
0153          analysisManager->FillH2(histo2, xx, yy, (value->second->sum_wx())/keV);
0154 }}}}
0155
0156 ofile << std::setprecision(6);
0157
0158 // Close the output ASCII file
0159 ofile.close();
0160
0161 // Close the output ROOT file
0162 analysisManager -> Write();
0163 analysisManager -> CloseFile();
0164 }