File indexing completed on 2025-02-23 09:21:54
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039 #include "ScoreStrandBreaks.hh"
0040
0041 #include "G4AnalysisManager.hh"
0042 #include "G4DNAChemistryManager.hh"
0043 #include "G4Event.hh"
0044 #include "G4EventManager.hh"
0045 #include "G4Scheduler.hh"
0046 #include "G4UImessenger.hh"
0047 #include "G4UnitsTable.hh"
0048
0049 #include <G4EventManager.hh>
0050 #include <G4SystemOfUnits.hh>
0051 #include <globals.hh>
0052
0053
0054
0055 ScoreStrandBreaks::ScoreStrandBreaks(G4String name, DetectorConstruction* detect, G4double* radius)
0056 : G4VPrimitiveScorer(name), G4UImessenger(), fpDetector(detect)
0057 {
0058 fpOutputFileUI = new G4UIcmdWithAString("/scorer/StrandBreak/OutputFile", this);
0059 fpOutputFileUI->SetGuidance("Set output file name");
0060
0061 fpOutputTypeUI = new G4UIcmdWithAString("/scorer/StrandBreak/OutputFormat", this);
0062 fpOutputTypeUI->SetGuidance("Set output file format: ASCII by default");
0063
0064 fpBreakEnergyUI = new G4UIcmdWithADoubleAndUnit("/scorer/StrandBreak/BreakEnergy", this);
0065 fpBreakEnergyUI->SetDefaultUnit("eV");
0066 fpBreakEnergyUI->SetGuidance("Direct SSB energy treshold");
0067
0068 fRadius = radius;
0069 fpGun = new MoleculeInserter(true);
0070
0071 }
0072
0073
0074
0075 ScoreStrandBreaks::~ScoreStrandBreaks()
0076 {
0077 delete fpOutputFileUI;
0078 delete fpBreakEnergyUI;
0079 delete fpOutputTypeUI;
0080 delete fpGun;
0081 }
0082
0083
0084
0085 void ScoreStrandBreaks::SetNewValue(G4UIcommand* command, G4String newValue)
0086 {
0087 if (command == fpOutputFileUI) {
0088 fOutputName = newValue;
0089 }
0090
0091 if (command == fpBreakEnergyUI) {
0092 fBreakEnergy = fpBreakEnergyUI->GetNewDoubleValue(newValue);
0093 }
0094
0095 if (command == fpOutputTypeUI) {
0096 fOutputType = newValue;
0097 G4StrUtil::to_lower(fOutputType);
0098 }
0099 }
0100
0101
0102
0103 G4bool ScoreStrandBreaks::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0104 {
0105 G4double edep = aStep->GetTotalEnergyDeposit();
0106 if (edep <= 0) {
0107 return FALSE;
0108 }
0109 fEnergyDeposit += edep;
0110
0111 G4StepPoint* aStepPoint = aStep->GetTrack()->GetStep()->GetPreStepPoint();
0112 G4TouchableHandle aTouchable = aStepPoint->GetTouchableHandle();
0113 G4String volumeName = aTouchable->GetVolume()->GetName();
0114 G4StrUtil::to_lower(volumeName);
0115 G4int strand = 0;
0116
0117 if (G4StrUtil::contains(volumeName, "deoxyribose")
0118 || G4StrUtil::contains(volumeName, "phosphate"))
0119 {
0120 if (G4StrUtil::contains(volumeName, "1"))
0121 strand = 1;
0122 else if (G4StrUtil::contains(volumeName, "2"))
0123 strand = 2;
0124
0125 G4int StrandID = strand;
0126 G4int BaseID = aTouchable->GetReplicaNumber(0);
0127 G4int PlasmidID = aTouchable->GetReplicaNumber(1);
0128 G4int EventID = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
0129
0130 fEnergyDepositMap[EventID][PlasmidID][BaseID][StrandID] += edep;
0131 return TRUE;
0132 }
0133
0134 if (!fDNAInserted) {
0135
0136 std::vector<std::vector<G4int>> DNADetails = fpDetector->GetDNADetails();
0137 std::vector<G4ThreeVector> DNAPositions = fpDetector->GetDNAPositions();
0138 std::vector<G4String> DNAMolecules = fpDetector->GetDNANames();
0139
0140 fpGun->Clean();
0141 for (size_t i = 0; i < DNAPositions.size(); i++) {
0142 fpGun->AddMolecule(DNAMolecules[i], DNAPositions[i], 1.0 * ps);
0143 }
0144 fpGun->DefineTracks();
0145
0146 fDNAInserted = true;
0147 }
0148 return FALSE;
0149 }
0150
0151
0152
0153 void ScoreStrandBreaks::Initialize(G4HCofThisEvent*) {}
0154
0155
0156
0157 void ScoreStrandBreaks::EndOfEvent(G4HCofThisEvent*)
0158 {
0159
0160 G4int EventID = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
0161
0162
0163 G4double volume = 4.0 / 3 * 3.14159 * (*fRadius) * (*fRadius) * (*fRadius) / 8;
0164 G4double density = 1.0 * g / cm3;
0165 G4double mass = density * volume;
0166 G4double Dose = fEnergyDeposit / mass;
0167 fDoseArray[EventID] = Dose / gray;
0168
0169
0170 for (auto EventAndElse : fEnergyDepositMap) {
0171 G4int event = EventAndElse.first;
0172 for (auto PlasmidAndElse : EventAndElse.second) {
0173 G4int plasmid = PlasmidAndElse.first;
0174 for (auto BaseAndElse : PlasmidAndElse.second) {
0175 G4int base = BaseAndElse.first;
0176 for (auto StrandAndEdep : BaseAndElse.second) {
0177 G4int strand = StrandAndEdep.first;
0178 G4double edep = StrandAndEdep.second;
0179 if (edep > fBreakEnergy) fDirectDamageMap[event].push_back({plasmid, base, strand});
0180 }
0181 }
0182 }
0183 }
0184 fEnergyDepositMap.clear();
0185
0186
0187 std::vector<std::vector<G4int>> DNADetails = fpDetector->GetDNADetails();
0188 std::vector<G4Track*> InsertedTracks = fpGun->GetInsertedTracks();
0189 std::vector<std::vector<G4int>> IndirectBreaks;
0190
0191 for (size_t i = 0; i < InsertedTracks.size(); i++) {
0192 if (InsertedTracks[i]->GetTrackStatus() == fStopAndKill)
0193 IndirectBreaks.push_back(DNADetails[i]);
0194 }
0195
0196 for (size_t i = 0; i < IndirectBreaks.size(); i++) {
0197 fIndirectDamageMap[EventID].push_back(IndirectBreaks[i]);
0198 }
0199
0200 fEnergyDeposit = 0;
0201 fDNAInserted = false;
0202 fnbOfEvents++;
0203 }
0204
0205
0206
0207 void ScoreStrandBreaks::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0208 {
0209 ScoreStrandBreaks* right =
0210 dynamic_cast<ScoreStrandBreaks*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0211
0212 if (right == 0) return;
0213 if (right == this) return;
0214
0215 DamageMap WorkerDirectDamageMap = right->fDirectDamageMap;
0216 DamageMap WorkerIndirectDamageMap = right->fIndirectDamageMap;
0217 std::map<G4int, G4double> WorkerDose = right->fDoseArray;
0218
0219 for (auto EventAndBreaks : WorkerDirectDamageMap) {
0220 G4int EventID = EventAndBreaks.first;
0221 std::vector<std::vector<G4int>> Breaks = EventAndBreaks.second;
0222 for (size_t i = 0; i < Breaks.size(); i++) {
0223 fDirectDamageMap[EventID].push_back(Breaks[i]);
0224 }
0225 }
0226
0227 for (auto EventAndBreaks : WorkerIndirectDamageMap) {
0228 G4int EventID = EventAndBreaks.first;
0229 std::vector<std::vector<G4int>> Breaks = EventAndBreaks.second;
0230 for (size_t i = 0; i < Breaks.size(); i++) {
0231 fIndirectDamageMap[EventID].push_back(Breaks[i]);
0232 }
0233 }
0234
0235 for (auto EventAndDose : WorkerDose) {
0236 G4int EventID = EventAndDose.first;
0237 G4double dose = EventAndDose.second;
0238 fDoseArray[EventID] = dose;
0239 }
0240
0241 fnbOfEvents += right->fnbOfEvents;
0242 right->Clear();
0243 }
0244
0245
0246
0247 void ScoreStrandBreaks::Clear()
0248 {
0249 if (fDirectDamageMap.size() != 0) fDirectDamageMap.clear();
0250
0251 if (fIndirectDamageMap.size() != 0) fIndirectDamageMap.clear();
0252
0253 if (fEnergyDepositMap.size() != 0) fEnergyDepositMap.clear();
0254
0255 if (fDoseArray.size() != 0) fDoseArray.clear();
0256 }
0257
0258
0259
0260 void ScoreStrandBreaks::DrawAll() {}
0261
0262
0263
0264 void ScoreStrandBreaks::PrintAll() {}
0265
0266
0267
0268 void ScoreStrandBreaks::OutputAndClear(G4double LET, G4double LET_STD)
0269 {
0270 if (G4Threading::IsWorkerThread()) return;
0271
0272 if (fOutputType != "ascii") {
0273 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0274
0275 if (fOutputType == "csv")
0276 analysisManager->SetDefaultFileType("csv");
0277
0278 else if (fOutputType == "root")
0279 analysisManager->SetDefaultFileType("root");
0280
0281 else if (fOutputType == "xml")
0282 analysisManager->SetDefaultFileType("xml");
0283
0284 WriteWithAnalysisManager(analysisManager);
0285 }
0286
0287 else
0288 ASCII(LET, LET_STD);
0289
0290 Clear();
0291 }
0292
0293
0294
0295 void ScoreStrandBreaks::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
0296 {
0297 analysisManager->OpenFile(fOutputName);
0298 int fNtupleID = analysisManager->CreateNtuple("SB", "Direct And Indirect SBs");
0299 analysisManager->CreateNtupleIColumn(fNtupleID, "EventID");
0300 analysisManager->CreateNtupleIColumn(fNtupleID, "PlasmidID");
0301 analysisManager->CreateNtupleIColumn(fNtupleID, "StrandID");
0302 analysisManager->CreateNtupleIColumn(fNtupleID, "BaseID");
0303 analysisManager->CreateNtupleIColumn(fNtupleID, "DamageID");
0304 analysisManager->CreateNtupleIColumn(fNtupleID, "BreakID");
0305 analysisManager->FinishNtuple(fNtupleID);
0306
0307 for (G4int i = 0; i < fnbOfEvents; i++) {
0308 if (fDirectDamageMap.find(i) != fDirectDamageMap.end()) {
0309 for (size_t j = 0; j < fDirectDamageMap[i].size(); j++) {
0310 std::vector<G4int> Break = fDirectDamageMap[i][j];
0311 analysisManager->FillNtupleIColumn(fNtupleID, 0, i);
0312 analysisManager->FillNtupleIColumn(fNtupleID, 1, Break[0]);
0313 analysisManager->FillNtupleIColumn(fNtupleID, 2, Break[1]);
0314 analysisManager->FillNtupleIColumn(fNtupleID, 3, Break[2]);
0315 analysisManager->FillNtupleIColumn(fNtupleID, 4, 1);
0316 analysisManager->AddNtupleRow(fNtupleID);
0317 }
0318 }
0319
0320 if (fIndirectDamageMap.find(i) != fIndirectDamageMap.end()) {
0321 for (size_t j = 0; j < fIndirectDamageMap[i].size(); j++) {
0322 std::vector<G4int> Break = fIndirectDamageMap[i][j];
0323 analysisManager->FillNtupleIColumn(fNtupleID, 0, i);
0324 analysisManager->FillNtupleIColumn(fNtupleID, 1, Break[0]);
0325 analysisManager->FillNtupleIColumn(fNtupleID, 2, Break[1]);
0326 analysisManager->FillNtupleIColumn(fNtupleID, 3, Break[2]);
0327 analysisManager->FillNtupleIColumn(fNtupleID, 4, 2);
0328 analysisManager->AddNtupleRow(fNtupleID);
0329 }
0330 }
0331 }
0332
0333 analysisManager->Write();
0334 analysisManager->CloseFile();
0335 }
0336
0337
0338
0339 void ScoreStrandBreaks::ASCII(G4double LET, G4double LET_STD)
0340 {
0341 std::ofstream dna(fOutputName + ".txt");
0342
0343 dna << "# DNA SB Map File" << G4endl;
0344 dna << "# LET = " << LET / (keV / um) << " +- " << LET_STD / (keV / um) << " keV / um" << G4endl;
0345 dna << "#" << std::setw(11) << "Event-ID" << std::setw(12) << "Plasmid-ID" << std::setw(12)
0346 << "BP-ID" << std::setw(12) << "Strand-ID" << std::setw(10) << "Break-ID" << std::setw(12)
0347 << "Dose (Gy) " << G4endl;
0348
0349 for (G4int i = 0; i < fnbOfEvents; i++) {
0350 if (fDirectDamageMap.find(i) != fDirectDamageMap.end()) {
0351 for (size_t j = 0; j < fDirectDamageMap[i].size(); j++) {
0352 std::vector<G4int> Break = fDirectDamageMap[i][j];
0353 dna << std::setw(12) << i << std::setw(12) << Break[0] << std::setw(12) << Break[1]
0354 << std::setw(12) << Break[2] << std::setw(10) << 1 << std::setw(12) << fDoseArray[i]
0355 << G4endl;
0356 }
0357 }
0358
0359 if (fIndirectDamageMap.find(i) != fIndirectDamageMap.end()) {
0360 for (size_t j = 0; j < fIndirectDamageMap[i].size(); j++) {
0361 std::vector<G4int> Break = fIndirectDamageMap[i][j];
0362 dna << std::setw(12) << i << std::setw(12) << Break[0] << std::setw(12) << Break[1]
0363 << std::setw(12) << Break[2] << std::setw(10) << 2 << std::setw(12) << fDoseArray[i]
0364 << G4endl;
0365 }
0366 }
0367 }
0368
0369 dna.close();
0370 }
0371
0372
