Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 09:22:02 

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 PhysAnalysis.cc
0028 /// \brief Implementation of the PhysAnalysis class
0029 
0030 #include "PhysAnalysis.hh"
0031 #include "G4AutoDelete.hh"
0032 #include "InformationKeeper.hh"
0033 
0034 #include <sstream>
0035 #include <fstream>
0036 
0037 
0038 G4ThreadLocal PhysAnalysis* the_analysis = nullptr;
0039 
0040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0041 
0042 PhysAnalysis* PhysAnalysis::GetAnalysis()
0043 {
0044     if (!the_analysis) {
0045         the_analysis = new PhysAnalysis();
0046         G4AutoDelete::Register(the_analysis);
0047     }
0048 
0049     return the_analysis;
0050 }
0051 
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053 
0054 void PhysAnalysis::OpenFile(const G4String& fname)
0055 {
0056     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0057     anManager->SetDefaultFileType("root");
0058     anManager->OpenFile(fname.c_str());
0059 }
0060 
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0062 
0063 void PhysAnalysis::Save()
0064 {
0065     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0066     anManager->Write();
0067 }
0068 
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0070 
0071 void PhysAnalysis::Close(G4bool reset)
0072 {
0073     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0074     anManager->CloseFile(reset);
0075 }
0076 
0077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0078 
0079 void PhysAnalysis::Book()
0080 {
0081     G4AnalysisManager* anManager = G4AnalysisManager::Instance();
0082     anManager->SetVerboseLevel(1);
0083 #ifdef G4MULTITHREADED
0084     // MT ntuple merging
0085     anManager->SetNtupleMerging(true);//for future use MT+MPI
0086 #endif
0087     anManager->SetNtupleDirectoryName("ntuple");
0088     // create ntuple
0089     if (anManager->GetFirstNtupleId() != 1) anManager->SetFirstNtupleId(1);
0090     anManager->CreateNtuple("ntuple_1","physical_stage");
0091     anManager->CreateNtupleDColumn("flagParticle");
0092     anManager->CreateNtupleDColumn("flagParentID");
0093     anManager->CreateNtupleDColumn("flagProcess");
0094     anManager->CreateNtupleDColumn("x");
0095     anManager->CreateNtupleDColumn("y");
0096     anManager->CreateNtupleDColumn("z");
0097     anManager->CreateNtupleDColumn("edep");
0098     anManager->CreateNtupleDColumn("eventNumber");
0099     anManager->CreateNtupleDColumn("volumeName");
0100     anManager->CreateNtupleDColumn("copyNumber");
0101     anManager->CreateNtupleDColumn("lastMetVoxelCopyNum");
0102     anManager->FinishNtuple(1);
0103 
0104     // For total edep
0105     anManager->CreateNtuple("ntuple_3","total_edep");
0106     anManager->CreateNtupleDColumn("eventNumber");
0107     anManager->CreateNtupleDColumn("edep");
0108     anManager->FinishNtuple(2);
0109 }
0110 
0111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0112 
0113 G4AnalysisManager* PhysAnalysis::GetAnalysisManager()
0114 {
0115     return G4AnalysisManager::Instance();
0116 }
0117 
0118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0119 
0120 void PhysAnalysis::AddInfoForChemGeo(InfoForChemGeo b)
0121 {
0122     fInfoForChemGeoVector.push_back(b);
0123 }
0124 
0125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0126 
0127 void PhysAnalysis::ClearVector()
0128 {
0129     fInfoForChemGeoVector.clear();
0130     fInfoInPhysStageVector.clear();
0131     fOutputFiles.clear();
0132 }
0133 
0134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0135 
0136 void PhysAnalysis::AddInfoInPhysStage(InfoInPhysStage b)
0137 {
0138     fInfoInPhysStageVector.push_back(b);
0139 }
0140 
0141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0142 
0143 void PhysAnalysis::UpdateChemInputDataAndFillNtuple()
0144 {
0145     //if (fInfoForChemGeoVector.size() == 0) return;
0146 
0147     // We will loop on all the "element" of the vector and create several files:
0148     // - New file for each event/voxel couple
0149     // Create a map to define all the output file
0150 
0151     for (auto const& binfo : fInfoForChemGeoVector) {
0152         if( binfo.fVolume == 161  // voxelStraight
0153             || binfo.fVolume == 162 // voxelRight
0154             || binfo.fVolume == 163 // voxelLeft
0155             || binfo.fVolume == 164 // voxelUp
0156             || binfo.fVolume == 165 // voxelDown
0157             || binfo.fVolume == 261 // voxelStraight2
0158             || binfo.fVolume == 262 // voxelRight2
0159             || binfo.fVolume == 263 // voxelLeft2
0160             || binfo.fVolume == 264 // voxelUp2
0161             || binfo.fVolume == 265) // voxelDown2
0162         {
0163             // We are in a voxel
0164 
0165             std::string voxelName("");
0166 
0167             if(binfo.fVolume==161) voxelName = "VoxelStraight";
0168             else if(binfo.fVolume==162) voxelName = "VoxelRight";
0169             else if(binfo.fVolume==163) voxelName = "VoxelLeft";
0170             else if(binfo.fVolume==164) voxelName = "VoxelUp";
0171             else if(binfo.fVolume==165) voxelName = "VoxelDown";
0172             else if(binfo.fVolume==261) voxelName = "VoxelStraight2";
0173             else if(binfo.fVolume==262) voxelName = "VoxelRight2";
0174             else if(binfo.fVolume==263) voxelName = "VoxelLeft2";
0175             else if(binfo.fVolume==264) voxelName = "VoxelUp2";
0176             else if(binfo.fVolume==265) voxelName = "VoxelDown2";
0177             
0178             // Get the event number
0179             int eventNum = int(binfo.fEventNumber);
0180             // Here we have all the information for one ntuple line
0181             // We want to know if we have to create a new output file (new eventNumber/voxel couple)
0182             // or if we already have one.
0183 
0184             // Check if the event has already been registered
0185             if(fOutputFiles.find(eventNum)==fOutputFiles.end() )
0186             {
0187                 // If not then create the event and voxel case
0188                 fOutputFiles[eventNum][binfo.fVolumeCopyNumber] = 
0189                 CreateChemInputFile(eventNum, int(binfo.fVolumeCopyNumber), voxelName);
0190             }
0191             // If the event has been registered then we need to check that the current voxel has also been
0192             // registered.
0193             else
0194             {
0195                 if(fOutputFiles[eventNum].find(binfo.fVolumeCopyNumber)==fOutputFiles[eventNum].end())
0196                 {
0197                     // We register the event-voxel couple
0198                     fOutputFiles[eventNum][binfo.fVolumeCopyNumber] = 
0199                     CreateChemInputFile(eventNum, int(binfo.fVolumeCopyNumber), voxelName);
0200                 }
0201             }
0202 
0203             // Write in the file the information needed by the chemistry simulation to build
0204             // the input water molecule or solvated electron.
0205             UpdatingChemInputFile(binfo);
0206         }
0207     }
0208 
0209     if (fInfoInPhysStageVector.size() == 0) return;
0210     for (auto const& binfo : fInfoInPhysStageVector) {
0211         // Check if the process is an ionisation
0212         // Only ionisation should trigger the removal of a DNA molecule from the chemical step
0213             if(binfo.fFlagProcess == 13
0214             || binfo.fFlagProcess == 113
0215             || binfo.fFlagProcess == 18
0216             || binfo.fFlagProcess == 21
0217             || binfo.fFlagProcess == 24
0218             || binfo.fFlagProcess == 27
0219             || binfo.fFlagProcess == 31) {
0220                 // Check the interaction happened in a dna molecule or its hydration shell
0221                 if( binfo.fVolumeName == 1 // d1
0222                     || binfo.fVolumeName == 11 // p1
0223                     || binfo.fVolumeName == 2 // d2
0224                     || binfo.fVolumeName == 22 // p2
0225                     || binfo.fVolumeName == 3 // cyto
0226                     || binfo.fVolumeName == 4 // gua
0227                     || binfo.fVolumeName == 5 // thy
0228                     || binfo.fVolumeName == 6 // ade
0229                     || binfo.fVolumeName == 7 // d1_w
0230                     || binfo.fVolumeName == 71 // p1_w
0231                     || binfo.fVolumeName == 8 // d2_w
0232                     || binfo.fVolumeName == 81 // p2_w
0233                     || binfo.fVolumeName == 9 // ade_w
0234                     || binfo.fVolumeName == 10 // gua_w
0235                     || binfo.fVolumeName == 13 // cyto_w
0236                     || binfo.fVolumeName == 12) // thy_w
0237                 {
0238                     // Retrieve the voxel copy number
0239                     double voxelCopyNumber  = binfo.fLastMetVoxelCopyNum;
0240                     // Get the event number
0241                     double eventNum = binfo.fEventNumber;
0242 
0243                     // Check if the event has already been registered
0244                     if(fOutputFiles.find(eventNum) != fOutputFiles.end() )
0245                     {
0246                         // Check if the volume number has been registered
0247                         if(fOutputFiles.at(eventNum).find(voxelCopyNumber) 
0248                             != fOutputFiles.at(eventNum).end() )
0249                         {
0250                             // If we are here then the event and volume couple has a already generated 
0251                             //file in which we should add a dna molecule to be removed
0252                             UpdatingChemInputFile(binfo);
0253                         }
0254                     }
0255                 }
0256         }
0257 
0258         // fill ntuple
0259         auto analysisManager =G4AnalysisManager::Instance();
0260         analysisManager->FillNtupleDColumn(1, 0, binfo.fFlagParticle);
0261         analysisManager->FillNtupleDColumn(1, 1, binfo.fFlagParentID);
0262         analysisManager->FillNtupleDColumn(1, 2, binfo.fFlagProcess);
0263         analysisManager->FillNtupleDColumn(1, 3, binfo.fX);
0264         analysisManager->FillNtupleDColumn(1, 4, binfo.fY);
0265         analysisManager->FillNtupleDColumn(1, 5, binfo.fZ);
0266         analysisManager->FillNtupleDColumn(1, 6, binfo.fEdep);
0267         analysisManager->FillNtupleDColumn(1, 7, binfo.fEventNumber);
0268         analysisManager->FillNtupleDColumn(1, 8, binfo.fVolumeName);
0269         analysisManager->FillNtupleDColumn(1, 9, binfo.fCopyNumber);
0270         analysisManager->FillNtupleDColumn(1, 10, binfo.fLastMetVoxelCopyNum);
0271         analysisManager->AddNtupleRow(1);
0272     }
0273 
0274 }
0275 
0276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0277 
0278 G4String PhysAnalysis::CreateChemInputFile(G4int eventNum, G4int volumeCopyNumber, const G4String &voxelName)
0279 {
0280     if (fOutputFolder.empty()) fOutputFolder = InformationKeeper::Instance()->GetChemInputFolderName();
0281     std::stringstream sstream;
0282     sstream<<"./"<<fOutputFolder<<"/event_"<<eventNum<<"_voxel_"<<volumeCopyNumber<<".dat";
0283     std::ofstream oFile;
0284     oFile.open(sstream.str().c_str() );
0285 
0286     oFile<<"_eventNum"<<"\t\t"<<eventNum<<std::endl;
0287     oFile<<"_voxelType"<<"\t\t"<<voxelName<<std::endl;
0288     oFile<<"_voxelCopyNumber"<<"\t\t"<<volumeCopyNumber<<std::endl;
0289     oFile<<std::endl;
0290 
0291     oFile
0292             <<"# Chemistry input informations\n"
0293             <<"# "<<"_input, type, state, electronicLevel, x, y, z, parentTrackID"<<std::endl;
0294     oFile<<"# type=1 -> water molecule && type=2 -> solvated electron"<<std::endl;
0295     oFile<<std::endl;
0296 
0297     oFile.close();
0298 
0299     return sstream.str().c_str();
0300 }
0301 
0302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0303 
0304 void PhysAnalysis::UpdatingChemInputFile(InfoForChemGeo b)
0305 {
0306     std::ofstream out;
0307     out.open( fOutputFiles[b.fEventNumber][b.fVolumeCopyNumber].c_str(), std::ios::app); // Open the file and go at the end
0308 
0309     if (b.fType==1)
0310     {
0311     out<<"_input"<<"\t"
0312                                             <<b.fType<<"\t\t"
0313                                             <<b.fState<<"\t\t"
0314                                         <<4-b.fElectronicLevel<<"\t\t"
0315                                         <<b.fRelX<<"\t\t"
0316                                         <<b.fRelY<<"\t\t"
0317                                         <<b.fRelZ<<"\t\t"
0318                                     <<b.fParentTrackID<<"\t\t"
0319                                 <<"\n";
0320     } 
0321     if (b.fType==2)
0322     {
0323     out<<"_input"<<"\t"
0324                                             <<b.fType<<"\t\t"
0325                                             <<b.fState<<"\t\t"
0326                                         <<b.fElectronicLevel<<"\t\t"
0327 // If we are here then the event and volume couple has a already 
0328 //generated file in which we should add a dna molecule to be removed
0329                                         <<b.fRelX<<"\t\t"
0330                                         <<b.fRelY<<"\t\t"
0331                                         <<b.fRelZ<<"\t\t"
0332                                     <<b.fParentTrackID<<"\t\t"
0333                                 <<"\n";
0334     }                              
0335     out.close();
0336 }
0337 
0338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0339 
0340 void PhysAnalysis::UpdatingChemInputFile(InfoInPhysStage b)
0341 {
0342     G4String name;
0343     G4double volumeName =  b.fVolumeName;
0344     if(volumeName == 1 || volumeName == 2 || volumeName == 7 || 
0345             volumeName == 8) name =  "Deoxyribose";
0346     else if(volumeName == 11 || volumeName == 22 || 
0347             volumeName == 71 || volumeName == 81) name =  "Phosphate";
0348     else if(volumeName == 6 || volumeName == 9) name =  "base";
0349     else if(volumeName == 4 || volumeName == 10) name =  "base";
0350     else if(volumeName == 5 || volumeName == 12) name =  "base";
0351     else if(volumeName == 3 || volumeName == 13) name =  "base";
0352     else
0353     {
0354         G4ExceptionDescription msg;
0355         msg <<"Volume number "<<volumeName<<" not registered.";
0356         G4Exception("PhysAnalysis::UpdatingChemInputFile", 
0357                     "", FatalException, msg);
0358     }
0359 
0360     G4double strand (-1);
0361 
0362     // Determine the strand
0363     if(volumeName==1
0364             || volumeName==11
0365             || volumeName==7
0366             || volumeName==71
0367             || volumeName==6 // ade
0368             || volumeName==9 // ade
0369             || volumeName==4 // gua
0370             || volumeName==10) // gua
0371         strand = 1;
0372     else if(volumeName==2
0373             || volumeName==22
0374             || volumeName==8
0375             || volumeName==81
0376             || volumeName==5 // thy
0377             || volumeName==12 // thy
0378             || volumeName==3 // cyto
0379             || volumeName==13) // cyto
0380         strand = 2;
0381     std::ofstream out;
0382     out.open(fOutputFiles[b.fEventNumber][b.fLastMetVoxelCopyNum].c_str(), std::ios::app); // Open the file and go at the end
0383     out<<"_remove"<<"\t"<<name<<"\t\t"<<b.fLastMetVoxelCopyNum<<"\t\t"<<strand<<"\t\t"<<std::endl; 
0384     out.close();
0385 }
0386 
0387 //....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!