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