Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /geant4/examples/advanced/dna/dsbandrepair/src/Analysis.cc was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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 Analysis.cc
0028 /// \brief Implementation of the Analysis class
0029 /// \file Analysis.cc
0030 /// \brief Implementation of the Analysis class
0031 
0032 #include "Analysis.hh"
0033 #include "DetectorConstruction.hh"
0034 #include "G4AutoDelete.hh"
0035 #include "G4Filesystem.hh"
0036 #include "G4DNAMolecule.hh"
0037 
0038 
0039 #include <sstream>
0040 #include <fstream>
0041 #ifdef USE_MPI
0042 #include "G4MPImanager.hh"
0043 #endif
0044 
0045 G4ThreadLocal Analysis* the_analysis = nullptr;
0046 
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048 
0049 Analysis* Analysis::GetAnalysis()
0050 {
0051     if (!the_analysis) {
0052         the_analysis = new Analysis();
0053         G4AutoDelete::Register(the_analysis);
0054     }
0055 
0056     return the_analysis;
0057 }
0058 
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0060 
0061 void Analysis::OpenFile(const G4String outFolder)
0062 {
0063     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0064     anManager->SetDefaultFileType("root");
0065     G4String fullFileName = fFileName;
0066     if (gRunMode == RunningMode::Chem) {
0067         G4String slash = "";
0068 #if defined(_WIN32) || defined(WIN32)
0069         slash= "\\";
0070 #else
0071         G4String slashu = "/";
0072         slash = slashu;
0073 #endif
0074         fullFileName = outFolder+slash+fFileName;
0075     }
0076 
0077     // open output file
0078     G4bool fileOpen = anManager->OpenFile(fullFileName.c_str());
0079     if (!fileOpen) {
0080         G4cout << "\n---> HistoManager::book(): cannot open " << fFileName
0081             << G4endl;
0082         return;
0083     } 
0084 }
0085 
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087 
0088 void Analysis::Save()
0089 {
0090     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0091     anManager->Write();
0092 }
0093 
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095 
0096 void Analysis::Close(G4bool reset)
0097 {
0098     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0099     anManager->CloseFile(reset);
0100 }
0101 
0102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0103 
0104 void Analysis::Book()
0105 {   
0106     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0107     anManager->SetVerboseLevel(1);
0108     if (anManager->GetFirstNtupleId() != 1) anManager->SetFirstNtupleId(1);
0109     
0110     if (gRunMode == RunningMode::Phys) {     
0111 #ifdef G4MULTITHREADED
0112         // MT ntuple merging
0113         anManager->SetNtupleMerging(true);//for future use MT+MPI
0114 #endif
0115         anManager->SetNtupleDirectoryName("ntuple");
0116         // create ntuple
0117         anManager->CreateNtuple("ntuple_1","physical_stage");
0118         anManager->CreateNtupleIColumn("flagParticle");
0119         anManager->CreateNtupleIColumn("flagParentID");
0120         anManager->CreateNtupleIColumn("flagProcess");
0121         anManager->CreateNtupleDColumn("x");
0122         anManager->CreateNtupleDColumn("y");
0123         anManager->CreateNtupleDColumn("z");
0124         anManager->CreateNtupleDColumn("edep");
0125         anManager->CreateNtupleIColumn("eventNumber");
0126         anManager->CreateNtupleIColumn("volumeName");
0127         anManager->CreateNtupleIColumn("copyNumber");
0128         anManager->CreateNtupleIColumn("lastMetVoxelCopyNum");
0129         anManager->FinishNtuple(1);
0130 
0131         // For total edep
0132         anManager->CreateNtuple("ntuple_3","total_edep");
0133         anManager->CreateNtupleIColumn("eventNumber");
0134         anManager->CreateNtupleDColumn("edep");
0135         anManager->FinishNtuple(2);
0136     }
0137 
0138     if (gRunMode == RunningMode::Chem) {
0139         // Create directories
0140         const G4String directoryName = "ntuple";
0141         if (anManager->GetNtupleDirectoryName() != directoryName) {
0142             anManager->SetNtupleDirectoryName(directoryName);
0143         }
0144         
0145         // DBScan
0146 
0147         anManager->CreateNtuple("ntuple_2","DB_chemical_stage");
0148         anManager->CreateNtupleIColumn(1,"strand");
0149         anManager->CreateNtupleIColumn(1,"copyNumber");
0150         anManager->CreateNtupleDColumn(1,"xp");
0151         anManager->CreateNtupleDColumn(1,"yp");
0152         anManager->CreateNtupleDColumn(1,"zp");
0153         anManager->CreateNtupleDColumn(1,"time");
0154         anManager->CreateNtupleIColumn(1,"base");
0155         anManager->FinishNtuple(1);
0156     }
0157 }
0158 
0159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0160 
0161 G4AnalysisManager* Analysis::GetAnalysisManager()
0162 {
0163     return G4AnalysisManager::Instance();
0164 }
0165 
0166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0167 
0168 void Analysis::AddInfoForChemGeo(InfoForChemGeo b)
0169 {
0170     fInfoForChemGeoVector.push_back(b);
0171 }
0172 
0173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0174 
0175 void Analysis::ClearVector()
0176 {
0177     fInfoForChemGeoVector.clear();
0178     fInfoInPhysStageVector.clear();
0179     fOutputFiles.clear();
0180 }
0181 
0182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0183 
0184 void Analysis::AddInfoInPhysStage(InfoInPhysStage b)
0185 {
0186     fInfoInPhysStageVector.push_back(b);
0187 }
0188 
0189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0190 
0191 void Analysis::UpdateChemInputDataAndFillNtuple()
0192 {
0193     //if (fInfoForChemGeoVector.size() == 0) return;
0194 
0195     // We will loop on all the "element" of the vector and create several files:
0196     // - New file for each event/voxel couple
0197     // Create a map to define all the output file
0198 
0199     for (auto const& binfo : fInfoForChemGeoVector) {
0200         if( binfo.fVolume == 161  // voxelStraight
0201             || binfo.fVolume == 162 // voxelRight
0202             || binfo.fVolume == 163 // voxelLeft
0203             || binfo.fVolume == 164 // voxelUp
0204             || binfo.fVolume == 165 // voxelDown
0205             || binfo.fVolume == 261 // voxelStraight2
0206             || binfo.fVolume == 262 // voxelRight2
0207             || binfo.fVolume == 263 // voxelLeft2
0208             || binfo.fVolume == 264 // voxelUp2
0209             || binfo.fVolume == 265) // voxelDown2
0210         {
0211             // We are in a voxel
0212 
0213             std::string voxelName("");
0214 
0215             if(binfo.fVolume==161) voxelName = "VoxelStraight";
0216             else if(binfo.fVolume==162) voxelName = "VoxelRight";
0217             else if(binfo.fVolume==163) voxelName = "VoxelLeft";
0218             else if(binfo.fVolume==164) voxelName = "VoxelUp";
0219             else if(binfo.fVolume==165) voxelName = "VoxelDown";
0220             else if(binfo.fVolume==261) voxelName = "VoxelStraight2";
0221             else if(binfo.fVolume==262) voxelName = "VoxelRight2";
0222             else if(binfo.fVolume==263) voxelName = "VoxelLeft2";
0223             else if(binfo.fVolume==264) voxelName = "VoxelUp2";
0224             else if(binfo.fVolume==265) voxelName = "VoxelDown2";
0225             
0226             // Get the event number
0227             int eventNum = int(binfo.fEventNumber);
0228             // Here we have all the information for one ntuple line
0229             // We want to know if we have to create a new output file (new eventNumber/voxel couple)
0230             // or if we already have one.
0231 
0232             // Check if the event has already been registered
0233             if(fOutputFiles.find(eventNum)==fOutputFiles.end() )
0234             {
0235                 // If not then create the event and voxel case
0236                 fOutputFiles[eventNum][binfo.fVolumeCopyNumber] = 
0237                 CreateChemInputFile(eventNum, int(binfo.fVolumeCopyNumber), voxelName);
0238             }
0239             // If the event has been registered then we need to check that the current voxel has also been
0240             // registered.
0241             else
0242             {
0243                 if(fOutputFiles[eventNum].find(binfo.fVolumeCopyNumber)==fOutputFiles[eventNum].end())
0244                 {
0245                     // We register the event-voxel couple
0246                     fOutputFiles[eventNum][binfo.fVolumeCopyNumber] = 
0247                     CreateChemInputFile(eventNum, int(binfo.fVolumeCopyNumber), voxelName);
0248                 }
0249             }
0250 
0251             // Write in the file the information needed by the chemistry simulation to build
0252             // the input water molecule or solvated electron.
0253             UpdatingChemInputFile(binfo);
0254         }
0255     }
0256 
0257     if (fInfoInPhysStageVector.size() == 0) return;
0258     for (auto const& binfo : fInfoInPhysStageVector) {
0259         // Check if the process is an ionisation
0260         // Only ionisation should trigger the removal of a DNA molecule from the chemical step
0261             if(binfo.fFlagProcess == 13
0262             || binfo.fFlagProcess == 113
0263             || binfo.fFlagProcess == 18
0264             || binfo.fFlagProcess == 21
0265             || binfo.fFlagProcess == 24
0266             || binfo.fFlagProcess == 27
0267             || binfo.fFlagProcess == 31) {
0268                 // Check the interaction happened in a dna molecule or its hydration shell
0269                 if( binfo.fVolumeName == 1 // d1
0270                     || binfo.fVolumeName == 11 // p1
0271                     || binfo.fVolumeName == 2 // d2
0272                     || binfo.fVolumeName == 22 // p2
0273                     || binfo.fVolumeName == 3 // cyto
0274                     || binfo.fVolumeName == 4 // gua
0275                     || binfo.fVolumeName == 5 // thy
0276                     || binfo.fVolumeName == 6 // ade
0277                     || binfo.fVolumeName == 7 // d1_w
0278                     || binfo.fVolumeName == 71 // p1_w
0279                     || binfo.fVolumeName == 8 // d2_w
0280                     || binfo.fVolumeName == 81 // p2_w
0281                     || binfo.fVolumeName == 9 // ade_w
0282                     || binfo.fVolumeName == 10 // gua_w
0283                     || binfo.fVolumeName == 13 // cyto_w
0284                     || binfo.fVolumeName == 12) // thy_w
0285                 {
0286                     // Retrieve the voxel copy number
0287                     double voxelCopyNumber  = binfo.fLastMetVoxelCopyNum;
0288                     // Get the event number
0289                     double eventNum = binfo.fEventNumber;
0290 
0291                     // Check if the event has already been registered
0292                     if(fOutputFiles.find(eventNum) != fOutputFiles.end() )
0293                     {
0294                         // Check if the volume number has been registered
0295                         if(fOutputFiles.at(eventNum).find(voxelCopyNumber) 
0296                             != fOutputFiles.at(eventNum).end() )
0297                         {
0298                             // If we are here then the event and volume couple has a already generated 
0299                             //file in which we should add a dna molecule to be removed
0300                             UpdatingChemInputFile(binfo);
0301                         }
0302                     }
0303                 }
0304         }
0305 
0306         // fill ntuple
0307         auto analysisManager =G4AnalysisManager::Instance();
0308         analysisManager->FillNtupleIColumn(1, 0, binfo.fFlagParticle);
0309         analysisManager->FillNtupleIColumn(1, 1, binfo.fFlagParentID);
0310         analysisManager->FillNtupleIColumn(1, 2, binfo.fFlagProcess);
0311         analysisManager->FillNtupleDColumn(1, 3, binfo.fX);
0312         analysisManager->FillNtupleDColumn(1, 4, binfo.fY);
0313         analysisManager->FillNtupleDColumn(1, 5, binfo.fZ);
0314         analysisManager->FillNtupleDColumn(1, 6, binfo.fEdep);
0315         analysisManager->FillNtupleIColumn(1, 7, binfo.fEventNumber);
0316         analysisManager->FillNtupleIColumn(1, 8, binfo.fVolumeName);
0317         analysisManager->FillNtupleIColumn(1, 9, binfo.fCopyNumber);
0318         analysisManager->FillNtupleIColumn(1, 10, binfo.fLastMetVoxelCopyNum);
0319         analysisManager->AddNtupleRow(1);
0320     }
0321 
0322 }
0323 
0324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0325 
0326 G4String Analysis::CreateChemInputFile(G4int eventNum, G4int volumeCopyNumber, const G4String &voxelName)
0327 {
0328     std::stringstream sstream;
0329     sstream<<"./"<<fChemInputFolderName<<"/event_"<<eventNum<<"_voxel_"<<volumeCopyNumber<<".dat";
0330     std::ofstream oFile;
0331     oFile.open(sstream.str().c_str() );
0332 
0333     oFile<<"_eventNum"<<"\t\t"<<eventNum<<std::endl;
0334     oFile<<"_voxelType"<<"\t\t"<<voxelName<<std::endl;
0335     oFile<<"_voxelCopyNumber"<<"\t\t"<<volumeCopyNumber<<std::endl;
0336     oFile<<std::endl;
0337 
0338     oFile
0339             <<"# Chemistry input informations\n"
0340             <<"# "<<"_input, type, state, electronicLevel, x, y, z, parentTrackID"<<std::endl;
0341     oFile<<"# type=1 -> water molecule && type=2 -> solvated electron"<<std::endl;
0342     oFile<<std::endl;
0343 
0344     oFile.close();
0345 
0346     return sstream.str().c_str();
0347 }
0348 
0349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0350 
0351 void Analysis::UpdatingChemInputFile(InfoForChemGeo b)
0352 {
0353     std::ofstream out;
0354     out.open( fOutputFiles[b.fEventNumber][b.fVolumeCopyNumber].c_str(), std::ios::app); // Open the file and go at the end
0355 
0356     if (b.fType==1)
0357     {
0358     out<<"_input"<<"\t"
0359                                             <<b.fType<<"\t\t"
0360                                             <<b.fState<<"\t\t"
0361                                         <<4-b.fElectronicLevel<<"\t\t"
0362                                         <<b.fRelX<<"\t\t"
0363                                         <<b.fRelY<<"\t\t"
0364                                         <<b.fRelZ<<"\t\t"
0365                                     <<b.fParentTrackID<<"\t\t"
0366                                 <<"\n";
0367     } 
0368     if (b.fType==2)
0369     {
0370     out<<"_input"<<"\t"
0371                                             <<b.fType<<"\t\t"
0372                                             <<b.fState<<"\t\t"
0373                                         <<b.fElectronicLevel<<"\t\t"
0374 // If we are here then the event and volume couple has a already 
0375 //generated file in which we should add a dna molecule to be removed
0376                                         <<b.fRelX<<"\t\t"
0377                                         <<b.fRelY<<"\t\t"
0378                                         <<b.fRelZ<<"\t\t"
0379                                     <<b.fParentTrackID<<"\t\t"
0380                                 <<"\n";
0381     }                              
0382     out.close();
0383 }
0384 
0385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0386 
0387 void Analysis::UpdatingChemInputFile(InfoInPhysStage b)
0388 {
0389     G4String name;
0390     G4double volumeName =  b.fVolumeName;
0391     if(volumeName == 1 || volumeName == 2 || volumeName == 7 || 
0392             volumeName == 8) name =  G4Deoxyribose::Definition()->GetName();
0393     else if(volumeName == 11 || volumeName == 22 || 
0394             volumeName == 71 || volumeName == 81) name =  G4Phosphate::Definition()->GetName();
0395     else if(volumeName == 6 || volumeName == 9) name =  G4Adenine::Definition()->GetName();
0396     else if(volumeName == 4 || volumeName == 10) name =  G4Guanine::Definition()->GetName();
0397     else if(volumeName == 5 || volumeName == 12) name =  G4Thymine::Definition()->GetName();
0398     else if(volumeName == 3 || volumeName == 13) name =  G4Cytosine::Definition()->GetName();
0399     else
0400     {
0401         G4ExceptionDescription msg;
0402         msg <<"Volume number "<<volumeName<<" not registered.";
0403         G4Exception("Analysis::UpdatingChemInputFile", 
0404                     "", FatalException, msg);
0405     }
0406 
0407     G4double strand (-1);
0408 
0409     // Determine the strand
0410     if(volumeName==1
0411             || volumeName==11
0412             || volumeName==7
0413             || volumeName==71
0414             || volumeName==6 // ade
0415             || volumeName==9 // ade
0416             || volumeName==4 // gua
0417             || volumeName==10) // gua
0418         strand = 1;
0419     else if(volumeName==2
0420             || volumeName==22
0421             || volumeName==8
0422             || volumeName==81
0423             || volumeName==5 // thy
0424             || volumeName==12 // thy
0425             || volumeName==3 // cyto
0426             || volumeName==13) // cyto
0427         strand = 2;
0428     std::ofstream out;
0429     out.open(fOutputFiles[b.fEventNumber][b.fLastMetVoxelCopyNum].c_str(), std::ios::app); // Open the file and go at the end
0430     out<<"_remove"<<"\t"<<name<<"\t\t"<<b.fCopyNumber<<"\t\t"<<strand<<"\t\t"<<std::endl; 
0431     out.close();
0432 }
0433 
0434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0435 
0436 void Analysis::WritePhysGeo()
0437 {
0438     //Create a file containing geopath and Physlist in "FS build directory" for chemStage and anlysis
0439     std::string fname = "imp.info";
0440     std::ofstream fout(fname);
0441     fout<<"====> Auto_generated file. Do not delete me !!!\n";
0442     fout<<"====> The file conveys some information for chem_geo and Analysis modules!!!\n";
0443     fout<<"_geocellpath "<<fCellDefFilePath<<"\n";
0444     for (const auto & entry : fVoxelDefFilesList) {
0445         fout<<"_geovolxelpath "<<std::string(entry)<<"\n";
0446     }
0447     fout<<"_numberOfBasepairs "<<fTotalNbBpPlacedInGeo<<"\n";
0448     fout<<"_numberOfHistones "<<fTotalNbHistonePlacedInGeo<<"\n";
0449     fout<<"_nucleusVolume "<<fNucleusVolume<<"\n"; // m3
0450     fout<<"_nucleusMassDensity "<<fNucleusMassDensity<<"\n"; // kg/m3
0451     fout<<"_nucleusMass "<<fNucleusVolume*fNucleusMassDensity; //kg
0452     fout.close();
0453 }
0454 
0455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0456 
0457 void Analysis::DefineCommands()
0458 {
0459     fMessenger = std::make_unique<G4GenericMessenger>(this,
0460                              "/dsbandrepair/output/",
0461                              "cmd controld");
0462     auto & fChemOutFolderCmd = fMessenger->DeclareProperty ("folderForChemOut",
0463                                                  fChemOutFolderName);
0464     fChemOutFolderCmd.SetParameterName("outFolderForChem",true);
0465     fChemOutFolderCmd.SetDefaultValue("chem_output");
0466 }
0467 
0468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0470 
0471 void Analysis::CheckAndCreateNewFolderInChemStage()
0472 {
0473     G4fs::file_status myFs = G4fs::file_status{};
0474     auto folderPath = G4fs::path{fChemOutFolderName.c_str()};
0475     auto isExist = G4fs::status_known(myFs) ? G4fs::exists(myFs) : G4fs::exists(folderPath);
0476     if (isExist && !G4fs::is_empty(folderPath)) {
0477         G4ExceptionDescription msg;
0478         msg <<"==>> Chem output folder "<<fChemOutFolderName
0479             <<" is already existing and not empty!!!\n";
0480         msg<<"==>> Please delete or rename it, or use other name for  "
0481         <<"Chem output folder in macro file!!!\n";
0482         G4Exception("Analysis::CheckAndCreateNewFolderInChemStage()","",FatalException,msg);
0483     }
0484     G4fs::create_directory(folderPath);
0485 }
0486 
0487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0488 
0489 void Analysis::CheckAndCreateNewFolderInPhysStage()
0490 {
0491     G4fs::file_status myFs = G4fs::file_status{};
0492     const G4fs::path folderPath{fPhysOutFolderName.c_str()};
0493     auto isExist = G4fs::status_known(myFs) ? G4fs::exists(myFs) : G4fs::exists(folderPath);
0494     if (isExist && !G4fs::is_empty(folderPath)) {
0495         G4fs::remove_all(folderPath);
0496     }
0497     G4fs::create_directory(folderPath);
0498 
0499     const G4fs::path folderPathPC{fChemInputFolderName.c_str()};
0500     isExist = G4fs::status_known(myFs) ? G4fs::exists(myFs) : G4fs::exists(folderPathPC);
0501     if (isExist && !G4fs::is_empty(folderPathPC)) {
0502         G4fs::remove_all(folderPathPC);
0503     }
0504     G4fs::create_directory(folderPathPC);
0505 }
0506 
0507 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!