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