dsbandrepair/src/ChemGeoImport.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 /// \file ChemGeoImport.cc
0028 /// \brief Implementation of the ChemGeoImport class
0029
0030 #include "ChemGeoImport.hh"
0031 #include "G4Filesystem.hh"
0032
0033 ChemGeoImport::ChemGeoImport()
0034 {
0035     GetVoxelDefFilePathList();
0036     fpGun = new UserMoleculeGun();
0037 }
0038
0039 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0040
0041 ChemGeoImport::~ChemGeoImport()
0042 {
0043     if(fpGun)
0044         delete fpGun;
0045 }
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 void ChemGeoImport::InsertMoleculeInWorld()
0050 {
0051     // The idea is to add all the molecules specified in the input files
0052
0053     if(fIsParsed)
0054     {
0055         // Create the molecules
0056
0057         // Loop on all the parsed molecules
0058         for(G4int i=0, ie=fMolecules.size(); i<ie; i++)
0059         {
0060             // Retrieve general molecule informations
0061             //
0062             G4String name = fMolecules[i].fName;
0063
0064             G4ThreeVector moleculePosition = fMolecules[i].fPosition;
0065             G4int copyNum = fMolecules[i].fCopyNumber;
0066
0067             G4int strand = fMolecules[i].fStrand;
0068             ChemMolecule Bmolecule(name, copyNum, moleculePosition, strand, -1, -1, -1);
0069
0070             if(name=="phosphate1")
0071             {
0072                 name="Phosphate";
0073                 Bmolecule.fName=name;
0074             }
0075             else if(name=="phosphate2")
0076             {
0077                 name="Phosphate";
0078                 Bmolecule.fName=name;
0079             }
0080             else if(name=="deoxyribose1")
0081             {
0082                 name="Deoxyribose";
0083                 Bmolecule.fName=name;
0084             }
0085             else if(name=="deoxyribose2")
0086             {
0087                 name="Deoxyribose";
0088                 Bmolecule.fName=name;
0089             }
0090             else if(name=="base_adenine")
0091             {
0092                 name="Adenine";
0093                 Bmolecule.fName="base";
0094             }
0095             else if(name=="base_guanine")
0096             {
0097                 name="Guanine";
0098                 Bmolecule.fName="base";
0099             }
0100             else if(name=="base_thymine")
0101             {
0102                 name="Thymine";
0103                 Bmolecule.fName="base";
0104             }
0105             else if(name=="base_cytosine")
0106             {
0107                 name="Cytosine";
0108                 Bmolecule.fName="base";
0109             }
0110             else if(name=="histone")
0111             {
0112                 name="Histone";
0113                 Bmolecule.fName=name;
0114             }
0115             else if(name=="solvatedElectron")
0116             {
0117                 name=G4Electron_aq::Definition()->GetName();
0118             }
0119             else if(name=="water")
0120             {
0121                 name=G4H2O::Definition()->GetName();
0122             }
0123             else
0124             {
0125                 G4String msg =
0126                 "The name "+ name+" is not specified in the listed chemical molecules";
0127                 G4Exception("ChemGeoImport::BuildGeometry", "", FatalException, msg);
0128             }
0129
0130
0131
0132             // Check if the molecule is on the "remove list"
0133             G4bool toBeRemoved = IsMoleculeInTheRemoveTable(Bmolecule);
0134             if(!toBeRemoved)
0135             {
0136                 // Molecule is not in the "remove list" and we can add it to the simulation
0137                 // Check the molecule to be added is not a water molecule (special case)
0138                 if(name != G4H2O::Definition()->GetName() )
0139                     fpGun->AddMolecule(name, moleculePosition, 1.e-12*s, copyNum, strand);
0140                 else // Water molecule case
0141                     fpGun->AddWaterMolecule(moleculePosition, fMolecules.at(i).fTrackId,
0142                         ElectronicModification(fMolecules.at(i).fState),
0143                         fMolecules.at(i).fElectronicLevel);
0144             }
0145
0146         }
0147         G4DNAChemistryManager::Instance()->SetGun(fpGun);
0148     }
0149     else
0150     {
0151         G4String msg =
0152         "ChemGeoImport::BuildGeometry: The parse method needs to be called first.";
0153         G4Exception("ChemGeoImport::BuildGeometry", "", FatalException, msg);
0154     }
0155 }
0156
0157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0158
0159 void ChemGeoImport::Reset()
0160 {
0161     // Clear the containers
0162     if(fpGun){
0163         delete fpGun;
0164         fpGun = new UserMoleculeGun();
0165     }
0166     fMolecules.clear();
0167     fIsParsed = false;
0168     fFactor = 1.;
0169 }
0170
0171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0172
0173 void ChemGeoImport::ParseFiles(const G4String& chemInputFile)
0174 {
0175     G4fs::path aP{std::string(chemInputFile)};
0176     if (G4fs::exists(aP)) {
0177         Reset();
0178         ParseChemInputFile(chemInputFile);
0179         auto geoPathFileName = GetVoxelDefFilePath(fGeoNameFromChemInput);
0180
0181         ParseGeoFile(geoPathFileName);
0182         fIsParsed = true;
0183     }
0184 }
0185
0186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0187
0188 void ChemGeoImport::ParseChemInputFile(const G4String& fileName)
0189 {
0190     // Setup the input stream
0191     std::ifstream file;
0192     file.open(fileName.c_str() );
0193
0194     if(!file.good() )
0195     {
0196         // Geant4 exception
0197         G4String msg = fileName+" could not be opened";
0198         G4Exception("ChemGeoImport::ParseChemInputFile", "", FatalException, msg);
0199     }
0200
0201     // Define the line string variable
0202     G4String line;
0203
0204     // Read the file line per line
0205     while(std::getline(file, line) )
0206     {
0207         // Check the line to determine if it is empty
0208         if(line.empty() )
0209             continue; // skip the line if it is empty
0210
0211         // Data string stream
0212         std::istringstream issLine(line);
0213
0214         // String to determine the first letter/word
0215         G4String firstItem;
0216
0217         // Put the first letter/word within the string
0218         issLine >> firstItem;
0219
0220         // Check first letter to determine if the line is data or comment
0221         if(firstItem=="#")
0222             continue; // skip the line if it is comment
0223
0224         else if(firstItem=="_input")
0225         {
0226             G4int type(-1), state(-1), electronicLevel(-1), parentTrackId(-1);
0227             G4double x, y, z;
0228             issLine >> type >> state >> electronicLevel;
0229             issLine >> x >> y >> z;
0230             issLine >> parentTrackId;
0231
0232             x *= fFactor*nm;
0233             y *= fFactor*nm;
0234             z *= fFactor*nm;
0235
0236             G4String name;
0237             if(type==1)
0238                 name="water";
0239             else if(type==2)
0240                 name="solvatedElectron";
0241             else
0242             {
0243                 G4ExceptionDescription description;
0244                 description <<  "The type " << type <<" is not recognized";
0245                 G4Exception("ChemGeoImport::ParseFile", "Fatal", FatalException, description, "");
0246             }
0247
0248             ChemMolecule molecule(name, -1, G4ThreeVector(x,y,z), -1,
0249                 state, electronicLevel, parentTrackId);
0250
0251             fMolecules.push_back(molecule);
0252         }
0253
0254         else if(firstItem=="_remove")
0255         {
0256             G4String name;
0257             issLine >> name;
0258
0259             G4int copyNumber;
0260             issLine >> copyNumber;
0261
0262             G4int strand;
0263             issLine >> strand;
0264
0265             fToBeRemovedMol.push_back(ChemMolecule(name,copyNumber,G4ThreeVector(),strand,-1,-1,-1));
0266         }
0267
0268         else if(firstItem=="_eventNum")
0269         {
0270             // Nothing
0271         }
0272
0273         else if(firstItem=="_voxelType")
0274         {
0275             issLine >> fGeoNameFromChemInput;
0276         }
0277
0278         else if(firstItem=="_voxelCopyNumber")
0279         {
0280             // Nothing
0281         }
0282
0283         else if(firstItem=="_Version")
0284         {
0285             // Nothing
0286         }
0287
0288         else
0289         {
0290             // Geant4 exception
0291             G4String msg =
0292             firstItem+" is not defined in the parser. Check the input file: "+fileName+".";
0293             G4Exception("ChemGeoImport::ParseChemInputFile", "Geo_WrongParse",FatalException, msg);
0294         }
0295     }
0296 }
0297
0298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0299
0300 void ChemGeoImport::ParseGeoFile(const G4String& fileName)
0301 {
0302     // Setup the input stream
0303     std::ifstream file(fileName.c_str());
0304
0305     // Check if the file was correctly opened
0306     if(!file.is_open() )
0307     {
0308         // Geant4 exception
0309         G4String msg = fileName+" could not be opened";
0310         G4Exception("ChemGeoImport::ParseGeoFile", "", FatalException, msg);
0311     }
0312
0313     // Define the line string variable
0314     G4String line;
0315
0316     // Read the file line per line
0317     while(std::getline(file, line) )
0318     {
0319         // Check the line to determine if it is empty
0320         if(line.empty() )
0321             continue; // skip the line if it is empty
0322
0323         // Data string stream
0324         std::istringstream issLine(line);
0325
0326         // String to determine the first letter/word
0327         G4String firstItem;
0328
0329         // Put the first letter/word within the string
0330         issLine >> firstItem;
0331
0332         // Check first letter to determine if the line is data or comment
0333         if(firstItem=="#")
0334             continue; // skip the line if it is comment
0335
0336         // Use the file
0337         else if(firstItem=="_Name")
0338         {
0339             G4String name;
0340             issLine >> name;
0341         }
0342         else if(firstItem=="_Size")
0343         {
0344             G4double size;
0345             issLine >> size;
0346             size *= fFactor*nm;
0347
0348             fSize = size;
0349         }
0350         else if(firstItem=="_Number")
0351         {
0352             // Nothing
0353         }
0354         else if(firstItem=="_Radius")
0355         {
0356             // Nothing
0357         }
0358         else if(firstItem=="_Version")
0359         {
0360             // Nothing
0361         }
0362         else if(firstItem=="_pl")
0363         {
0364             G4String name;
0365             issLine >> name;
0366
0367             G4String material;
0368             issLine >> material;
0369
0370             G4int strand;
0371             issLine >> strand;
0372
0373             G4int copyNumber;
0374             issLine >> copyNumber;
0375
0376             G4double x;
0377             issLine >> x;
0378             x *= fFactor*nm;
0379
0380             G4double y;
0381             issLine >> y;
0382             y *= fFactor*nm;
0383
0384             G4double z;
0385             issLine >> z;
0386             z *= fFactor*nm;
0387
0388             ChemMolecule molecule(name, copyNumber, G4ThreeVector(x, y, z), strand, -1, -1, -1);
0389
0390             fMolecules.push_back(molecule);
0391         }
0392
0393         else
0394         {
0395             // Geant4 exception
0396             G4String msg =
0397             firstItem+" is not defined in the parser. Check the input file: "+fileName+".";
0398             G4Exception("ChemGeoImport::ParseGeoFile", "Geo_WrongParse", FatalException, msg);
0399         }
0400     }
0401 }
0402
0403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0404
0405 G4bool ChemGeoImport::IsMoleculeInTheRemoveTable(const ChemMolecule& molecule)
0406 {
0407     if(std::find(fToBeRemovedMol.begin(),fToBeRemovedMol.end(),molecule) != fToBeRemovedMol.end())
0408         return true;
0409     else
0410         return false;
0411 }
0412
0413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0414
0415 G4String ChemGeoImport::GetVoxelDefFilePath(G4String bareName)
0416 {
0417     G4String strRes = "";
0418     for (auto const &entry : fVoxelDefFilesList) {
0419         G4fs::path voxelP{std::string(entry)};
0420         if (voxelP.stem().string() == bareName) {
0421             strRes = entry;
0422         }
0423     }
0424     return strRes;
0425 }
0426
0427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0428
0429 void ChemGeoImport::GetVoxelDefFilePathList()
0430 {
0431     G4fs::path thisP = G4fs::current_path();
0432     G4bool doesWantedFileExist = false;
0433     for (const auto &entry : G4fs::directory_iterator(thisP)){
0434         if (entry.path().filename() == "imp.info") {
0435             std::ifstream file(entry.path().c_str());
0436             if(!file.good() ){
0437                 G4String msg =
0438                 "File imp.info is broken. Check its content or try to rerun the PhysicalStage?";
0439                 G4Exception("ChemGeoImport::GetVoxelDefFilePathList()", "", FatalException, msg);
0440             }
0441             doesWantedFileExist = true;
0442             G4String line;
0443             while(std::getline(file, line) ){
0444                 std::istringstream iss(line);
0445                 G4String flag;
0446                 G4String voxelDefFile;
0447                 iss >> flag;
0448                 if ( flag == "_geovolxelpath") {
0449                     iss >> voxelDefFile;
0450                     fVoxelDefFilesList.insert(voxelDefFile);
0451                 }
0452             }
0453             file.close();
0454         }
0455     }
0456
0457     if (!doesWantedFileExist) {
0458         G4String msg = "File imp.info does not exist. Did you run the Physical Stage?";
0459         G4Exception("ChemGeoImport::GetVoxelDefFilePathList()", "", FatalException, msg);
0460     }
0461 }
0462
0463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......