| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// file:
28 /// brief:
29
30 #include "AnalysisManager.hh"
31
32 #include "AnalysisMessenger.hh"
33 #include "ChromosomeHit.hh"
34 #include "DNAGeometry.hh"
35 #include "DNAHit.hh"
36 #include "DetectorConstruction.hh"
37 #include "UtilityFunctions.hh"
38
39 #include "G4Event.hh"
40 #include "G4EventManager.hh"
41 #include "G4MolecularConfiguration.hh"
42 #include "G4Run.hh"
43 #include "G4RunManager.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "Randomize.hh"
46
47 #include <fstream>
48 #include <utility>
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52 AnalysisManager::AnalysisManager()
53 {
54 fAnalysisManager = G4AnalysisManager::Instance();
55 fpAnalysisMessenger = new AnalysisMessenger(this);
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59
60 AnalysisManager::~AnalysisManager()
61 {
62 delete fpAnalysisMessenger;
63 }
64
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66
67 void AnalysisManager::Initialize()
68 {
69 const G4Run* run = G4RunManager::GetRunManager()->GetCurrentRun();
70 if (run->GetRunID() == 0) {
71 G4cout << "AnalysisManager::Initialize() GetRunID : " << run->GetRunID() << G4endl;
72 fAnalysisManager->SetDefaultFileType("root");
73 fAnalysisManager->SetVerboseLevel(1);
74
75 fAnalysisManager->SetNtupleDirectoryName("tuples");
76 fAnalysisManager->SetHistoDirectoryName("hists");
77
78 fAnalysisManager->CreateH1("ssb_counts", "SSBs", 1000, 0, 1000);
79 fAnalysisManager->CreateH1("ssb_energies_ev", "SSB Energies", 500, 0, 500);
80 // fAnalysisManager->CreateH1("fragments", "Fragment sizes", 500, 0,
81 // 500);//ORG
82 fAnalysisManager->CreateH1("fragments", "Fragment sizes", 10000, 0,
83 10000); // dousatsu
84 fAnalysisManager->CreateH3("local_pos", "Strand Interaction Positions", 50, -25, 25, 50, -25,
85 25, 50, -25, 25);
86 fAnalysisManager->CreateH3("global_pos", "Strand Interaction Positions", 50, -100, 100, 50,
87 -100, 100, 50, -100, 100);
88 fAnalysisManager->CreateH1("e_cell_kev", "Energy Deposits in Cell", 2500, 0,
89 2500); // dousatsu
90
91 fAnalysisManager->CreateNtuple("primary_source", "Primary Source");
92 fAnalysisManager->CreateNtupleSColumn("Primary");
93 fAnalysisManager->CreateNtupleDColumn("Energy");
94 fAnalysisManager->CreateNtupleDColumn("PosX_um");
95 fAnalysisManager->CreateNtupleDColumn("PosY_um");
96 fAnalysisManager->CreateNtupleDColumn("PosZ_um");
97 fAnalysisManager->CreateNtupleDColumn("MomX");
98 fAnalysisManager->CreateNtupleDColumn("MomY");
99 fAnalysisManager->CreateNtupleDColumn("MomZ");
100 fAnalysisManager->CreateNtupleDColumn("StopPosX_um");
101 fAnalysisManager->CreateNtupleDColumn("StopPosY_um");
102 fAnalysisManager->CreateNtupleDColumn("StopPosZ_um");
103 fAnalysisManager->CreateNtupleDColumn("TraLen_cell_um");
104 fAnalysisManager->CreateNtupleDColumn("TraLen_chro_um");
105 fAnalysisManager->FinishNtuple();
106
107 fAnalysisManager->CreateNtuple("chromosome_hits", "Energy Deposits in Chromosomes");
108 fAnalysisManager->CreateNtupleSColumn("chromosome");
109 fAnalysisManager->CreateNtupleDColumn("e_chromosome_kev");
110 fAnalysisManager->CreateNtupleDColumn("e_dna_kev");
111 fAnalysisManager->FinishNtuple();
112
113 fAnalysisManager->CreateNtuple("classification", "Break Complexity Frequency");
114 fAnalysisManager->CreateNtupleSColumn("Primary");
115 fAnalysisManager->CreateNtupleDColumn("Energy");
116 fAnalysisManager->CreateNtupleIColumn("None");
117 fAnalysisManager->CreateNtupleIColumn("SSB");
118 fAnalysisManager->CreateNtupleIColumn("SSBp");
119 fAnalysisManager->CreateNtupleIColumn("2SSB");
120 fAnalysisManager->CreateNtupleIColumn("DSB");
121 fAnalysisManager->CreateNtupleIColumn("DSBp");
122 fAnalysisManager->CreateNtupleIColumn("DSBpp");
123 fAnalysisManager->FinishNtuple();
124
125 fAnalysisManager->CreateNtuple("source", "Break Source Frequency");
126 fAnalysisManager->CreateNtupleSColumn("Primary");
127 fAnalysisManager->CreateNtupleDColumn("Energy");
128 fAnalysisManager->CreateNtupleIColumn("None");
129 fAnalysisManager->CreateNtupleIColumn("SSBd");
130 fAnalysisManager->CreateNtupleIColumn("SSBi");
131 fAnalysisManager->CreateNtupleIColumn("SSBm");
132 fAnalysisManager->CreateNtupleIColumn("DSBd");
133 fAnalysisManager->CreateNtupleIColumn("DSBi");
134 fAnalysisManager->CreateNtupleIColumn("DSBm");
135 fAnalysisManager->CreateNtupleIColumn("DSBh");
136 fAnalysisManager->FinishNtuple();
137
138 fAnalysisManager->CreateNtuple("damage", "DNA damage locations");
139 fAnalysisManager->CreateNtupleIColumn("Event");
140 fAnalysisManager->CreateNtupleSColumn("Primary");
141 fAnalysisManager->CreateNtupleDColumn("Energy");
142 fAnalysisManager->CreateNtupleSColumn("TypeClassification");
143 fAnalysisManager->CreateNtupleSColumn("SourceClassification");
144 fAnalysisManager->CreateNtupleDColumn("Position_x_um");
145 fAnalysisManager->CreateNtupleDColumn("Position_y_um");
146 fAnalysisManager->CreateNtupleDColumn("Position_z_um");
147 fAnalysisManager->CreateNtupleDColumn("Size_nm");
148 fAnalysisManager->CreateNtupleIColumn("FragmentLength");
149 fAnalysisManager->CreateNtupleIColumn("BaseDamage");
150 fAnalysisManager->CreateNtupleIColumn("StrandDamage");
151 fAnalysisManager->CreateNtupleIColumn("DirectBreaks");
152 fAnalysisManager->CreateNtupleIColumn("IndirectBreaks");
153 fAnalysisManager->CreateNtupleIColumn("EaqBaseHits");
154 fAnalysisManager->CreateNtupleIColumn("EaqStrandHits");
155 fAnalysisManager->CreateNtupleIColumn("OHBaseHits");
156 fAnalysisManager->CreateNtupleIColumn("OHStrandHits");
157 fAnalysisManager->CreateNtupleIColumn("HBaseHits");
158 fAnalysisManager->CreateNtupleIColumn("HStrandHits");
159 fAnalysisManager->CreateNtupleDColumn("EnergyDeposited_eV");
160 fAnalysisManager->CreateNtupleIColumn("InducedBreaks");
161 fAnalysisManager->CreateNtupleIColumn("Chain");
162 fAnalysisManager->CreateNtupleIColumn("Strand"); // WG
163 fAnalysisManager->CreateNtupleDColumn("BasePair"); // dousatsu
164 fAnalysisManager->CreateNtupleSColumn("Name");
165 fAnalysisManager->FinishNtuple();
166 }
167
168 if (!fFileName) {
169 fFileName = "molecular-dna";
170 }
171 fAnalysisManager->OpenFile(fFileName);
172 }
173
174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
175
176 void AnalysisManager::Close()
177 {
178 fAnalysisManager->Write();
179 fAnalysisManager->CloseFile();
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183
184 void AnalysisManager::ProcessDNAHitsVector(const std::vector<const DNAHit*>& dnaHits)
185 {
186 // Check we have DNA Geometry (runs once)
187 if (fpDNAGeometry == nullptr) {
188 auto det = dynamic_cast<const DetectorConstruction*>(
189 G4RunManager::GetRunManager()->GetUserDetectorConstruction());
190 fpDNAGeometry = det->GetDNAGeometry();
191 }
192
193 if (dnaHits.empty()) {
194 return;
195 }
196 const G4Event* event = G4EventManager::GetEventManager()->GetConstCurrentEvent();
197 // SSBs
198 for (auto dnaHit : dnaHits) {
199 fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("ssb_energies_ev"),
200 dnaHit->GetEnergy() / eV);
201 if (fChainToSave == -1) {
202 // save all chains
203 fAnalysisManager->FillH3(
204 fAnalysisManager->GetH3Id("local_pos"), dnaHit->GetLocalPosition().getX() / nm,
205 dnaHit->GetLocalPosition().getY() / nm, dnaHit->GetLocalPosition().getZ() / nm);
206 fAnalysisManager->FillH3(fAnalysisManager->GetH3Id("global_pos"),
207 dnaHit->GetPosition().getX() / nm, dnaHit->GetPosition().getY() / nm,
208 dnaHit->GetPosition().getZ() / nm);
209 }
210 else if (fChainToSave == dnaHit->GetChainIdx()) { // what is different ?
211 // save only specified chain
212 fAnalysisManager->FillH3(
213 fAnalysisManager->GetH3Id("local_pos"), dnaHit->GetLocalPosition().getX() / nm,
214 dnaHit->GetLocalPosition().getY() / nm, dnaHit->GetLocalPosition().getZ() / nm);
215 fAnalysisManager->FillH3(fAnalysisManager->GetH3Id("global_pos"),
216 dnaHit->GetPosition().getX() / nm, dnaHit->GetPosition().getY() / nm,
217 dnaHit->GetPosition().getZ() / nm);
218 }
219 }
220 fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("ssb_counts"), dnaHits.size());
221
222 // DSBs
223 // Sorting
224 // Make a map relating chromosome/chain to its binary tree
225 std::map<std::pair<G4String, G4int>, BinaryTree*> treemap;
226
227 // Populate the binary tree
228 // G4cout << "Building Binary Trees for " << dnaHits.size() << " Nodes" <<
229 // G4endl;
230 for (auto it = dnaHits.begin(); it != dnaHits.end(); ++it) {
231 std::pair<G4String, G4int> key = std::make_pair((*it)->GetChromosome(), (*it)->GetChainIdx());
232 if (treemap.find(key) == treemap.end()) {
233 treemap[key] = new BinaryTree();
234 if (fAnalysisManager->GetVerboseLevel() >= 2) {
235 G4cout << "Constructing binary tree. Chromosome: " << key.first << " Chain: " << key.second
236 << G4endl;
237 }
238 }
239 treemap.at(key)->Insert(*it);
240 if (fAnalysisManager->GetVerboseLevel() >= 2) {
241 G4cout << "Adding hit to binary tree. Chromosome: " << key.first << " Chain: " << key.second
242 << G4endl;
243 }
244 }
245
246 // Analysis
247 // Create a vector to stock damage records for the event
248 // G4cout << "Building Damage Records" << G4endl;
249 std::vector<DamageRecord*> damageRecords;
250 for (auto& it : treemap) {
251 if (fAnalysisManager->GetVerboseLevel() >= 2) {
252 G4cout << "Analysing hits for chromosome: " << it.first.first << " Chain: " << it.first.second
253 << G4endl;
254 }
255 DNAHit* currentHit = it.second->First();
256 DNAHit* nextHit = currentHit;
257
258 // Runs while their are still hits
259 while (nextHit) {
260 // Create record for this fragment
261 if (currentHit->GetBasePairIdx() > 9223372036854775807) {
262 G4cout << " SEE AnalysisManager !!!" << G4endl;
263 abort();
264 } // dousatsu
265 // if (currentHit->GetBasePairIdx()>2147483647){G4cout<<" SEE
266 // AnalysisManager !!!"<<G4endl;abort();}//ORG
267 auto idx = (int64_t)currentHit->GetBasePairIdx(); // dousatsu
268 int64_t startidx = (idx > 10) ? (idx - 10) : 0; // dousatsu
269 // G4int idx = currentHit->GetBasePairIdx();//ORG
270 // G4int startidx = (idx > 10) ? (idx - 10) : 0;//ORG
271 G4String name = it.first.first + "_" + std::to_string(it.first.second) + "_"
272 + std::to_string(startidx); // check!!!!
273 auto* record =
274 new DamageRecord(name, idx, currentHit->GetPlacementIdx(), currentHit->GetChainIdx());
275 damageRecords.push_back(record);
276 BasePairDamageRecord* bp;
277 int64_t gap = 0; // dousatsu
278
279 // bookend with 10 bps on either side of no damage
280 record->AddEmptyBPDamage(idx - startidx);
281 // continues until fragment ends
282 while (true) {
283 bp = new BasePairDamageRecord;
284 bp->fStrand1Energy = currentHit->GetStrand1Energy();
285 bp->fStrand2Energy = currentHit->GetStrand2Energy();
286 bp->fBp1Energy = currentHit->GetBP1Energy();
287 bp->fBp2Energy = currentHit->GetBP2Energy();
288 bp->fBp1IndirectEvt = (currentHit->GetBase1Rad() != nullptr);
289 bp->fBp2IndirectEvt = (currentHit->GetBase2Rad() != nullptr);
290 bp->fStrand1IndirectEvt = (currentHit->GetStrand1Rad() != nullptr);
291 bp->fStrand2IndirectEvt = (currentHit->GetStrand2Rad() != nullptr);
292 bp->fbp1DirectDmg = false; // No physical damage on BPs
293 bp->fbp2DirectDmg = false; // No physical damage on BPs
294
295 bp->fStrand1DirectDmg =
296 fpDNAGeometry->GetDamageModel()->IsDirectStrandBreak(bp->fStrand1Energy);
297 bp->fStrand2DirectDmg =
298 fpDNAGeometry->GetDamageModel()->IsDirectStrandBreak(bp->fStrand2Energy);
299
300 // if(bp->fStrand1DirectDmg) bp->fStrandIdx = 1;
301 // else if(bp->fStrand2DirectDmg) strandID = 2;
302
303 if (bp->fBp1IndirectEvt) {
304 bp->fBp1IndirectDmg =
305 fpDNAGeometry->GetDamageModel()->IsIndirectBaseDamage(currentHit->GetBase1Rad());
306 if (bp->fBp1IndirectDmg) {
307 bp->fbp1InducedBreak =
308 fpDNAGeometry->GetDamageModel()->IsInducedStrandBreak(currentHit->GetBase1Rad());
309 }
310 }
311
312 if (bp->fBp2IndirectEvt) {
313 bp->fBp2IndirectDmg =
314 fpDNAGeometry->GetDamageModel()->IsIndirectBaseDamage(currentHit->GetBase2Rad());
315 if (bp->fBp2IndirectDmg) {
316 bp->fbp2InducedBreak =
317 fpDNAGeometry->GetDamageModel()->IsInducedStrandBreak(currentHit->GetBase2Rad());
318 }
319 }
320
321 if (bp->fStrand1IndirectEvt) {
322 bp->fStrand1IndirectDmg =
323 fpDNAGeometry->GetDamageModel()->IsIndirectStrandDamage(currentHit->GetStrand1Rad());
324 }
325 if (bp->fStrand2IndirectEvt) {
326 bp->fStrand2IndirectDmg =
327 fpDNAGeometry->GetDamageModel()->IsIndirectStrandDamage(currentHit->GetStrand2Rad());
328 }
329 // Record radical-base/strand damage.
330 // it is possible for base/strand damage to not correspond to
331 // direct hits due to the damage induction probabilities.
332 // This mainly affects strands.
333 if (bp->fBp1IndirectDmg) {
334 record->AddBaseHit(currentHit->GetBase1Rad()->GetDefinition());
335 }
336 if (bp->fBp2IndirectDmg) {
337 record->AddBaseHit(currentHit->GetBase2Rad()->GetDefinition());
338 }
339 if (bp->fStrand1IndirectDmg) {
340 record->AddStrandHit(currentHit->GetStrand1Rad()->GetDefinition());
341 // strandID = 1;
342 }
343 if (bp->fStrand2IndirectDmg) {
344 record->AddStrandHit(currentHit->GetStrand2Rad()->GetDefinition());
345 // strandID = 2;
346 }
347
348 record->AddBasePairDamage(bp, currentHit->GetPosition());
349
350 nextHit = it.second->Next(currentHit);
351 if (nextHit == nullptr) {
352 break;
353 }
354
355 gap = nextHit->GetBasePairIdx() - currentHit->GetBasePairIdx();
356 if (fFragmentGap > 0) {
357 // case 1: continuous strand
358 if (gap > fFragmentGap) {
359 currentHit = nextHit;
360 break;
361 }
362 else {
363 record->AddEmptyBPDamage(gap);
364 }
365 }
366 else if (fFragmentGap == 0) {
367 // case 2: individual placements, not joined
368 // ie. plasmids etc. Each placement is a separate strand
369 if (currentHit->GetPlacementIdx() != nextHit->GetPlacementIdx()) {
370 if (fpDNAGeometry->GetVerbosity() > 1) {
371 G4cout << "Analysis passing to a new placement" << G4endl;
372 }
373 currentHit = nextHit;
374 break;
375 }
376 else {
377 record->AddEmptyBPDamage(gap);
378 }
379 }
380 else {
381 G4Exception("MolecularAnalaysisManager", "ERR_FRAGMENT_VALUE", FatalException,
382 "The value set in /analysisDNA/fragmentGap is bad");
383 }
384 currentHit = nextHit;
385 }
386 record->AddEmptyBPDamage(10);
387 // end bookend
388 }
389 }
390 // Count the number of breaks and print out the records
391 // std::map<G4String, G4int> breakmap;
392 std::map<complexityEnum, G4int> breakmap;
393
394 breakmap[DSBplusplus] = 0;
395 breakmap[DSBplus] = 0;
396 breakmap[DSB] = 0;
397 breakmap[twoSSB] = 0;
398 breakmap[SSBplus] = 0;
399 breakmap[SSB] = 0;
400 breakmap[NoneComplexity] = 0;
401 std::map<sourceEnum, G4int> sourcemap;
402 sourcemap[SSBd] = 0;
403 sourcemap[SSBi] = 0;
404 sourcemap[SSBm] = 0;
405 sourcemap[DSBh] = 0;
406 sourcemap[DSBm] = 0;
407 sourcemap[DSBd] = 0;
408 sourcemap[DSBi] = 0;
409 sourcemap[undefined] = 0;
410
411 for (auto& damageRecord : damageRecords) {
412 DamageClassification* classif = damageRecord->GetClassification(fDSBDistance);
413 breakmap[classif->fComplexity]++;
414 sourcemap[classif->fSource]++;
415 fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("fragments"), damageRecord->GetSize());
416
417 ///////////dousatsu--------------------------------------------
418 fAnalysisManager->FillNtupleIColumn(4, 0, event->GetEventID());
419 fAnalysisManager->FillNtupleSColumn(
420 4, 1, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
421 fAnalysisManager->FillNtupleDColumn(
422 4, 2, event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
423
424 G4String complexityString = "None";
425 switch (classif->fComplexity) {
426 case SSB:
427 complexityString = "SSB";
428 break;
429 case SSBplus:
430 complexityString = "SSB+";
431 break;
432 case twoSSB:
433 complexityString = "2SSB";
434 break;
435 case DSB:
436 complexityString = "DSB";
437 break;
438 case DSBplus:
439 complexityString = "DSB+";
440 break;
441 case DSBplusplus:
442 complexityString = "DSB++";
443 break;
444 default:
445 complexityString = "None";
446 break;
447 }
448 fAnalysisManager->FillNtupleSColumn(4, 3, complexityString);
449
450 G4String sourceString = "undefined";
451 switch (classif->fSource) {
452 case SSBd:
453 sourceString = "SSBd";
454 break;
455 case SSBi:
456 sourceString = "SSBi";
457 break;
458 case SSBm:
459 sourceString = "SSBm";
460 break;
461 case DSBh:
462 sourceString = "DSBh";
463 break;
464 case DSBm:
465 sourceString = "DSBm";
466 break;
467 case DSBd:
468 sourceString = "DSBd";
469 break;
470 case DSBi:
471 sourceString = "DSBi";
472 break;
473 default:
474 sourceString = "undefined";
475 break;
476 }
477
478 fAnalysisManager->FillNtupleSColumn(4, 4, sourceString);
479 fAnalysisManager->FillNtupleDColumn(4, 5, damageRecord->GetMeanPosition().getX() / um);
480 fAnalysisManager->FillNtupleDColumn(4, 6, damageRecord->GetMeanPosition().getY() / um);
481 fAnalysisManager->FillNtupleDColumn(4, 7, damageRecord->GetMeanPosition().getZ() / um);
482 fAnalysisManager->FillNtupleDColumn(4, 8, damageRecord->GetMeanDistance() / nm);
483 fAnalysisManager->FillNtupleIColumn(4, 9, damageRecord->GetSize());
484 fAnalysisManager->FillNtupleIColumn(4, 10, classif->fbaseDmg);
485 fAnalysisManager->FillNtupleIColumn(4, 11, classif->fStrandDmg);
486 fAnalysisManager->FillNtupleIColumn(4, 12, classif->fDirectBreaks);
487 fAnalysisManager->FillNtupleIColumn(4, 13, classif->fIndirectBreaks);
488 fAnalysisManager->FillNtupleIColumn(4, 14, damageRecord->GetEaqBaseHits());
489 fAnalysisManager->FillNtupleIColumn(4, 15, damageRecord->GetEaqStrandHits());
490 fAnalysisManager->FillNtupleIColumn(4, 16, damageRecord->GetOHBaseHits());
491 fAnalysisManager->FillNtupleIColumn(4, 17, damageRecord->GetOHStrandHits());
492 fAnalysisManager->FillNtupleIColumn(4, 18, damageRecord->GetHBaseHits());
493 fAnalysisManager->FillNtupleIColumn(4, 19, damageRecord->GetHStrandHits());
494 fAnalysisManager->FillNtupleDColumn(4, 20, damageRecord->GetEnergy() / eV);
495 fAnalysisManager->FillNtupleIColumn(4, 21, classif->fInducedBreaks);
496 fAnalysisManager->FillNtupleIColumn(4, 22, damageRecord->GetChainIdx());
497 fAnalysisManager->FillNtupleIColumn(4, 23, damageRecord->GetPlacementIdx());
498 fAnalysisManager->FillNtupleDColumn(4, 24, (G4double)damageRecord->GetStartBPIdx());
499 // fAnalysisManager->FillNtupleIColumn(4, 23, (*it)->GetStartBPIdx());//ORG
500 fAnalysisManager->FillNtupleSColumn(4, 25, damageRecord->GetName());
501 fAnalysisManager->AddNtupleRow(4);
502 delete classif;
503 //
504 // TODO: Use (*it)->GetPlacementIdx() to work out if the damage
505 // is continuously joined to a past event or not.
506 // If yes, save the fragment length.
507 //
508 if (fSaveStrands) {
509 damageRecord->PrintRecord(fStrandDirectory + "/" + std::to_string(event->GetEventID()) + "_"
510 + damageRecord->GetName() + ".txt");
511 }
512 if (fAnalysisManager->GetVerboseLevel() >= 2) damageRecord->PrintRecord("");
513 }
514
515 fAnalysisManager->FillNtupleSColumn(
516 2, 0, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
517 fAnalysisManager->FillNtupleDColumn(2, 1,
518 event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
519 fAnalysisManager->FillNtupleIColumn(2, 2, breakmap[NoneComplexity]);
520 fAnalysisManager->FillNtupleIColumn(2, 3, breakmap[SSB]);
521 fAnalysisManager->FillNtupleIColumn(2, 4, breakmap[SSBplus]);
522 fAnalysisManager->FillNtupleIColumn(2, 5, breakmap[twoSSB]);
523 fAnalysisManager->FillNtupleIColumn(2, 6, breakmap[DSB]);
524 fAnalysisManager->FillNtupleIColumn(2, 7, breakmap[DSBplus]);
525 fAnalysisManager->FillNtupleIColumn(2, 8, breakmap[DSBplusplus]);
526 fAnalysisManager->AddNtupleRow(2);
527
528 fAnalysisManager->FillNtupleSColumn(
529 3, 0, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
530 fAnalysisManager->FillNtupleDColumn(3, 1,
531 event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
532 fAnalysisManager->FillNtupleIColumn(3, 2, sourcemap[undefined]);
533 fAnalysisManager->FillNtupleIColumn(3, 3, sourcemap[SSBd]);
534 fAnalysisManager->FillNtupleIColumn(3, 4, sourcemap[SSBi]);
535 fAnalysisManager->FillNtupleIColumn(3, 5, sourcemap[SSBm]);
536 fAnalysisManager->FillNtupleIColumn(3, 6, sourcemap[DSBd]);
537 fAnalysisManager->FillNtupleIColumn(3, 7, sourcemap[DSBi]);
538 fAnalysisManager->FillNtupleIColumn(3, 8, sourcemap[DSBm]);
539 fAnalysisManager->FillNtupleIColumn(3, 9, sourcemap[DSBh]);
540 fAnalysisManager->AddNtupleRow(3);
541
542 for (auto it : damageRecords) {
543 delete it;
544 }
545
546 // Cleanup, delete binary trees
547 for (const auto& it : treemap) {
548 delete it.second;
549 }
550 }
551
552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553
554 void AnalysisManager::ProcessChromosomeHitMap(const std::map<uint32_t, ChromosomeHit*>& chromosomes)
555 {
556 for (auto chromosome : chromosomes) {
557 fAnalysisManager->FillNtupleSColumn(1, 0, chromosome.second->GetName());
558 fAnalysisManager->FillNtupleDColumn(1, 1, chromosome.second->GetChromosomeEdep() / keV);
559 fAnalysisManager->FillNtupleDColumn(1, 2, chromosome.second->GetDNAEdep() / keV);
560 fAnalysisManager->AddNtupleRow(1);
561 }
562 }
563
564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
565
566 void AnalysisManager::ProcessPrimary(const G4ThreeVector& primstoppos, const G4double& tlcell,
567 const G4double& tlchro)
568 {
569 const G4Event* event = G4EventManager::GetEventManager()->GetConstCurrentEvent();
570 fAnalysisManager->FillNtupleSColumn(
571 0, 0, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
572 fAnalysisManager->FillNtupleDColumn(0, 1,
573 event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
574 fAnalysisManager->FillNtupleDColumn(0, 2, event->GetPrimaryVertex()->GetX0() / um);
575 fAnalysisManager->FillNtupleDColumn(0, 3, event->GetPrimaryVertex()->GetY0() / um);
576 fAnalysisManager->FillNtupleDColumn(0, 4, event->GetPrimaryVertex()->GetZ0() / um);
577
578 G4ThreeVector mom = event->GetPrimaryVertex()->GetPrimary()->GetMomentumDirection();
579 fAnalysisManager->FillNtupleDColumn(0, 5, mom.x() / mom.mag());
580 fAnalysisManager->FillNtupleDColumn(0, 6, mom.y() / mom.mag());
581 fAnalysisManager->FillNtupleDColumn(0, 7, mom.z() / mom.mag());
582 fAnalysisManager->FillNtupleDColumn(0, 8, primstoppos.x() / um);
583 fAnalysisManager->FillNtupleDColumn(0, 9, primstoppos.y() / um);
584 fAnalysisManager->FillNtupleDColumn(0, 10, primstoppos.z() / um);
585 fAnalysisManager->FillNtupleDColumn(0, 11, tlcell / um);
586 fAnalysisManager->FillNtupleDColumn(0, 12, tlchro / um);
587 fAnalysisManager->AddNtupleRow(0);
588 }
589
590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
591
592 void AnalysisManager::ProcessCellEdep(const G4double& edep)
593 {
594 fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("e_cell_kev"), edep / keV);
595 }
596
597 ////////////////////////////////////////////////////////////////////////////////
598 ////// Binary Tree Methods
599 ////////////////////////////////////////////////////////////////////////////////
600
601 BinaryTree::BinaryTree()
602 {
603 fRoot = nullptr;
604 }
605
606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
607
608 BinaryTree::~BinaryTree()
609 {
610 Destroy_tree();
611 }
612
613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
614
615 // Makes own internal copy of hit
616 void BinaryTree::Insert(const DNAHit* hit)
617 {
618 auto newHit = new DNAHit(*hit);
619 if (fRoot == nullptr) {
620 Node* node = new Node;
621 node->fleft = nullptr;
622 node->fright = nullptr;
623 node->fparent = nullptr;
624 node->fdata = newHit;
625 node->fkey = newHit->GetBasePairIdx();
626 fRoot = node;
627 }
628 else {
629 Insert_(newHit, fRoot);
630 }
631 }
632
633 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
634
635 DNAHit* BinaryTree::Search(int64_t index) // dousatsu
636 {
637 return Search_(index, fRoot);
638 }
639
640 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
641
642 void BinaryTree::Destroy_tree()
643 {
644 Destroy_tree_(fRoot);
645 }
646
647 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
648
649 void BinaryTree::Destroy_tree_(Node* node)
650 {
651 if (node != nullptr) {
652 Destroy_tree_(node->fleft);
653 Destroy_tree_(node->fright);
654 delete node->fdata;
655 delete node;
656 }
657 }
658
659 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
660
661 void BinaryTree::Insert_(DNAHit* hit, Node* node)
662 {
663 if (hit->GetBasePairIdx() > 9223372036854775807) {
664 G4cout << " SEE AnalysisManager !!!" << G4endl;
665 G4ExceptionDescription exceptionDescription;
666 exceptionDescription << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
667 G4cout << "!!! aborting ... " << G4endl;
668 G4cout << "This pair ID is so big " << hit->GetBasePairIdx() << G4endl;
669 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
670 G4Exception(
671 "BinaryTree"
672 "::insert()",
673 "insert000", FatalException, exceptionDescription);
674 }
675 int64_t key = hit->GetBasePairIdx(); // dousatsu
676 // G4int key = hit->GetBasePairIdx();//ORG
677 if (key < node->fkey) {
678 if (node->fleft != nullptr) {
679 Insert_(hit, node->fleft);
680 }
681 else {
682 Node* daughter = new Node;
683 daughter->fparent = node;
684 daughter->fleft = nullptr;
685 daughter->fright = nullptr;
686 daughter->fdata = hit;
687 daughter->fkey = hit->GetBasePairIdx();
688 node->fleft = daughter;
689 }
690 }
691 else if (key > node->fkey) {
692 if (node->fright != nullptr) {
693 Insert_(hit, node->fright);
694 }
695 else {
696 Node* daughter = new Node;
697 daughter->fparent = node;
698 daughter->fleft = nullptr;
699 daughter->fright = nullptr;
700 daughter->fdata = hit;
701 daughter->fkey = hit->GetBasePairIdx();
702 node->fright = daughter;
703 }
704 }
705 else {
706 node->fdata->AddHit(*hit);
707 }
708 }
709
710 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
711
712 DNAHit* BinaryTree::Search_(int64_t key, Node* node) // dousatsu check
713 {
714 if (node != nullptr) {
715 if (key < node->fkey) {
716 return Search_(key, node->fleft);
717 }
718 else if (key > node->fkey) {
719 return Search_(key, node->fright);
720 }
721 else {
722 return node->fdata;
723 }
724 }
725 else {
726 return nullptr;
727 }
728 }
729
730 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
731
732 DNAHit* BinaryTree::First() const
733 {
734 if (fRoot == nullptr) {
735 return nullptr;
736 }
737 else {
738 return First_(fRoot);
739 }
740 }
741
742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
743
744 DNAHit* BinaryTree::First_(Node* node) const
745 {
746 if (node->fleft == nullptr) {
747 return node->fdata;
748 }
749 else {
750 return First_(node->fleft);
751 }
752 }
753
754 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
755
756 DNAHit* BinaryTree::Next(const DNAHit* hit) const
757 {
758 if (hit->GetBasePairIdx() > 9223372036854775807) {
759 G4cout << " SEE AnalysisManager !!!" << G4endl;
760 G4ExceptionDescription exceptionDescription;
761 exceptionDescription << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
762 G4cout << "!!! aborting ... " << G4endl;
763 G4cout << "This pair ID is so big " << hit->GetBasePairIdx() << G4endl;
764 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
765 G4Exception(
766 "BinaryTree"
767 "::insert()",
768 "insert000", FatalException, exceptionDescription);
769 }
770 int64_t key = hit->GetBasePairIdx(); // dousatsu
771 if (fRoot == nullptr) {
772 return nullptr;
773 }
774 else {
775 return Next_(key, fRoot);
776 }
777 }
778
779 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
780
781 // Looks more complicated than it is
782 // The interface to these functions with integer keys rather than node
783 // objects makes this look more complicated than standard algorithms
784 DNAHit* BinaryTree::Next_(int64_t key, Node* node) const // dousatsu
785 {
786 if (key < node->fkey) {
787 if (node->fleft != nullptr) {
788 return Next_(key, node->fleft);
789 }
790 else // left is NULL
791 {
792 return node->fdata;
793 }
794 }
795 else // (key >= node->key)
796 {
797 if (node->fright != nullptr) {
798 return Next_(key, node->fright);
799 }
800 else // right is NULL; solution is higher in tree
801 {
802 while (true) {
803 node = node->fparent;
804 if (node == nullptr) {
805 return nullptr;
806 }
807 if (key < node->fkey) {
808 return node->fdata;
809 }
810 }
811 }
812 }
813 }
814
815 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
816
817 ////////////////////////////////////////////////////////////////////////////////
818 ////////// Damage Record
819 ////////////////////////////////////////////////////////////////////////////////
820
821 const char* DamageRecord::fDirectDamageChar = "D";
822 const char* DamageRecord::fIndirectDamageChar = "I";
823 const char* DamageRecord::fHitNoDamageChar = "~";
824 const char* DamageRecord::fNotHitChar = "-";
825 const char* DamageRecord::fBothDamageChar = "X";
826
827 DamageRecord::DamageRecord(const G4String& name, int64_t startIndex, G4int place_idx,
828 G4int chain_idx)
829 : fName(name), fStartIndex(startIndex), fStartPlacement(place_idx), fChainIdx(chain_idx)
830 {}
831
832 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
833
834 DamageRecord::~DamageRecord()
835 {
836 for (auto& fDamageRecord : fDamageRecords) {
837 delete fDamageRecord;
838 }
839 }
840
841 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
842
843 void DamageRecord::PrintRecord(const G4String& filename, const G4double& dsbDistance)
844 {
845 std::stringstream strand1;
846 std::stringstream bp1;
847 std::stringstream bp2;
848 std::stringstream strand2;
849 for (auto& fDamageRecord : fDamageRecords) {
850 strand1 << GetChar(fDamageRecord->fStrand1DirectDmg, fDamageRecord->fStrand1IndirectDmg,
851 fDamageRecord->fStrand1Energy);
852 bp1 << GetChar(fDamageRecord->fbp1DirectDmg, fDamageRecord->fBp1IndirectDmg,
853 fDamageRecord->fBp1Energy);
854 bp2 << GetChar(fDamageRecord->fbp2DirectDmg, fDamageRecord->fBp2IndirectDmg,
855 fDamageRecord->fBp2Energy);
856 strand2 << GetChar(fDamageRecord->fStrand2DirectDmg, fDamageRecord->fStrand2IndirectDmg,
857 fDamageRecord->fStrand2Energy);
858 }
859
860 // print to with no filename
861 if (filename.empty()) {
862 DamageClassification* classification = this->GetClassification(dsbDistance);
863 G4cout << "Sequences starts at base pair " << fStartIndex << " in placement " << fStartPlacement
864 << ", Chain " << fChainIdx << ". Total length: " << fDamageRecords.size()
865 << ". Classification: " << classification->fComplexity << " " << classification->fSource
866 << G4endl;
867 G4cout << strand1.str() << G4endl;
868 G4cout << bp1.str() << G4endl;
869 G4cout << bp2.str() << G4endl;
870 G4cout << strand2.str() << G4endl;
871 return;
872 }
873 // To Text output
874 std::fstream fs(filename, std::fstream::in | std::fstream::out | std::fstream::trunc);
875 if (fs.fail()) {
876 G4cout << "Could not open filestream" << G4endl;
877 }
878 DamageClassification* classification = this->GetClassification(dsbDistance);
879 fs << "Sequences starts at base pair " << fStartIndex << " in placement " << fStartPlacement
880 << ", Chain " << fChainIdx << ". Total length: " << fDamageRecords.size()
881 << ". Classification: " << classification->fComplexity << " " << classification->fSource
882 << std::endl;
883 fs << strand1.str() << std::endl;
884 fs << bp1.str() << std::endl;
885 fs << bp2.str() << std::endl;
886 fs << strand2.str() << std::endl;
887 fs.close();
888 delete classification;
889 }
890
891 const char* DamageRecord::GetChar(const G4bool& direct, const G4bool& indirect, const G4double& e)
892 {
893 if (direct && indirect) {
894 return fBothDamageChar;
895 }
896 else if (direct) {
897 return fDirectDamageChar;
898 }
899 else if (indirect) {
900 return fIndirectDamageChar;
901 }
902 else if (e > 0) {
903 return fHitNoDamageChar;
904 }
905 else {
906 return fNotHitChar;
907 }
908 }
909
910 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
911
912 DamageClassification* DamageRecord::GetClassification(const G4double& dsbDistance)
913 {
914 /* Explanation of class because we do some sneaky bit-level manipulation to
915 track breaks along the strand.
916
917 Things to note straight off
918 1) This routine makes no promise to count SSBs and DSBs correctly
919 once they cease being relevant to the break complexity. That is, once
920 a DSB has occurred, the SSB count becomes irrelevant
921 2) The previous ten base pairs are kept in the ints lastTenStrandX
922 Each bit in the number represents one base pair precedent. The bit
923 is set if the strand was broken.
924 3) The countX variables are used to track the number of breaks in the
925 preceeding 10 base pairs on each strand. This is needed to calculate
926 DSB+ events
927 4) The lambda function count_bytes is an algorithm to count the number
928 of set bytes in constant time for 32-bit integers. It yields for
929 us the count of breaks in the previous ten bases.
930
931 Classification of complexity:
932 DSB++ 2DSBs in fragment
933 DSB+ 1 DSB and then an SSB within 10 base pairs
934 DSB Each strand hit once within <=10 base pairs
935 2SSB Each strand hit once with >10 base pairs separation
936 SSB+ One strand hit twice or more
937 SSB One strand hit only once
938 None Nothing hit
939
940 Classification by Source
941
942 */
943 // lambda function to count the number of bytes set in a 32-bit uint
944 // https://blogs.msdn.microsoft.com/jeuge/2005/06/08/bit-fiddling-3/
945 // Note: integers beginning with 0 are entered as octal!
946 auto count_bytes = [](uint32_t u) {
947 uint32_t uCount = u - ((u >> 1) & 033333333333) - ((u >> 2) & 011111111111);
948 return ((uCount + (uCount >> 3)) & 030707070707) % 63;
949 };
950
951 auto classification = new DamageClassification();
952
953 G4int nDSB = 0;
954 G4int nSSB = 0;
955 // Source classification
956 G4int oldnDSB = 0;
957 G4int nDSBm = 0;
958 G4int nDSBi = 0;
959 G4int nDSBd = 0;
960 G4int nDSBh = 0;
961
962 // Complexity Classification
963 G4int nDSBPlus = 0;
964 G4bool SSB1 = false;
965 G4bool SSB2 = false;
966 // Counters for the last ten bp
967 uint32_t lastTenDirectStrand1 = 0;
968 uint32_t lastTenDirectStrand2 = 0;
969 uint32_t lastTenIndirectStrand1 = 0;
970 uint32_t lastTenIndirectStrand2 = 0;
971 uint32_t lastTenStrand1 = 0;
972 uint32_t lastTenStrand2 = 0;
973 uint32_t lastTenTracked1 = 0; // We remove entries from these if they have
974 uint32_t lastTenTracked2 = 0; // been counted in DSBs
975 uint32_t count1 = 0;
976 uint32_t count2 = 0;
977 G4bool strand1Indirect = false;
978 G4bool strand1Direct = false;
979 G4bool strand2Indirect = false;
980 G4bool strand2Direct = false;
981 G4int baseDamage = 0;
982 G4int strandDamage = 0;
983 uint32_t truncator = std::pow(2, dsbDistance + 1) - 1;
984 for (G4int ii = 0; ii != (G4int)fDamageRecords.size(); ii++) {
985 lastTenDirectStrand1 = (lastTenDirectStrand1 << 1) & truncator;
986 lastTenDirectStrand2 = (lastTenDirectStrand2 << 1) & truncator;
987 lastTenIndirectStrand1 = (lastTenIndirectStrand1 << 1) & truncator;
988 lastTenIndirectStrand2 = (lastTenIndirectStrand2 << 1) & truncator;
989
990 lastTenStrand1 = lastTenStrand1 << 1;
991 lastTenStrand2 = lastTenStrand2 << 1;
992 lastTenTracked1 = lastTenTracked1 << 1;
993 lastTenTracked2 = lastTenTracked2 << 1;
994 // DSB distance by default is 10
995 // 2047 = 0b11111111111, ie. truncation to 11bp
996 lastTenStrand1 = lastTenStrand1 & truncator;
997 lastTenStrand2 = lastTenStrand2 & truncator;
998 lastTenTracked1 = lastTenTracked1 & truncator;
999 lastTenTracked2 = lastTenTracked2 & truncator;
1000 // keep counters to ten binary places
1001 count1 = count_bytes(lastTenStrand1);
1002 count2 = count_bytes(lastTenStrand2);
1003
1004 // Nightmare of if statements
1005 if (fDamageRecords.at(ii)->fStrand1DirectDmg) {
1006 strand1Direct = true;
1007 classification->fDirectBreaks++;
1008 }
1009 if (fDamageRecords.at(ii)->fStrand1IndirectDmg) {
1010 strand1Indirect = true;
1011 classification->fIndirectBreaks++;
1012 }
1013 if (fDamageRecords.at(ii)->fStrand2DirectDmg) {
1014 strand2Direct = true;
1015 classification->fDirectBreaks++;
1016 }
1017 if (fDamageRecords.at(ii)->fStrand2IndirectDmg) {
1018 strand2Indirect = true;
1019 classification->fIndirectBreaks++;
1020 }
1021 if (fDamageRecords.at(ii)->fbp1DirectDmg) {
1022 // classification->fDirectBreaks++; //???
1023 }
1024 if (fDamageRecords.at(ii)->fBp1IndirectDmg) {
1025 if (fDamageRecords.at(ii)->fbp1InducedBreak) {
1026 // Counts as a hit only if it breaks a strand
1027 classification->fIndirectBreaks++;
1028 strand1Indirect = true;
1029 classification->fInducedBreaks++;
1030 }
1031 }
1032 if (fDamageRecords.at(ii)->fbp2DirectDmg) {
1033 // classification->fDirectBreaks++; //???
1034 }
1035 if (fDamageRecords.at(ii)->fBp2IndirectDmg) {
1036 if (fDamageRecords.at(ii)->fbp2InducedBreak) {
1037 // Counts as a hit only if it breaks a strand
1038 classification->fIndirectBreaks++;
1039 strand2Indirect = true;
1040 classification->fInducedBreaks++;
1041 }
1042 }
1043
1044 G4bool s1IndirectBreak =
1045 fDamageRecords.at(ii)->fStrand1IndirectDmg || fDamageRecords.at(ii)->fbp1InducedBreak;
1046 G4bool s2IndirectBreak =
1047 fDamageRecords.at(ii)->fStrand2IndirectDmg || fDamageRecords.at(ii)->fbp2InducedBreak;
1048
1049 // strand 1 hit
1050 G4bool strand1hit = fDamageRecords.at(ii)->fStrand1DirectDmg || s1IndirectBreak;
1051 G4bool strand2hit = fDamageRecords.at(ii)->fStrand2DirectDmg || s2IndirectBreak;
1052 G4bool bp1hit = fDamageRecords.at(ii)->fbp1DirectDmg || fDamageRecords.at(ii)->fBp1IndirectDmg;
1053 G4bool bp2hit = fDamageRecords.at(ii)->fbp2DirectDmg || fDamageRecords.at(ii)->fBp2IndirectDmg;
1054 strandDamage += ((G4int)strand1hit + (G4int)strand2hit);
1055 baseDamage += ((G4int)bp1hit + (G4int)bp2hit);
1056
1057 // Update the damage type counters
1058 // strand 1
1059 if (s1IndirectBreak) {
1060 lastTenIndirectStrand1++;
1061 }
1062 if (fDamageRecords.at(ii)->fStrand1DirectDmg) {
1063 lastTenDirectStrand1++;
1064 }
1065 // strand 2
1066 if (s2IndirectBreak) {
1067 lastTenIndirectStrand2++;
1068 }
1069 if (fDamageRecords.at(ii)->fStrand2DirectDmg) lastTenDirectStrand2++;
1070
1071 oldnDSB = nDSB + nDSBPlus;
1072 if (strand1hit && strand2hit) {
1073 nDSB++;
1074 if (count1 || count2) {
1075 nDSBPlus++;
1076 }
1077 lastTenStrand1++;
1078 lastTenStrand2++;
1079 }
1080 else if (strand1hit) {
1081 if (lastTenTracked2) {
1082 nDSB++;
1083 lastTenTracked2 = (1 << (utility::Fls(lastTenTracked2) - 1)) ^ lastTenTracked2;
1084 if ((count2 > 1) || (count1)) {
1085 nDSBPlus++;
1086 }
1087 }
1088 else if (count1 && count2) {
1089 nDSBPlus++;
1090 }
1091 else {
1092 nSSB++;
1093 }
1094 SSB1 = true;
1095 lastTenTracked1++;
1096 lastTenStrand1++;
1097 }
1098 else if (strand2hit) {
1099 if (lastTenTracked1) {
1100 nDSB++;
1101 lastTenTracked1 = (1 << (utility::Fls(lastTenTracked1) - 1)) ^ lastTenTracked1;
1102 if ((count1 > 1) || (count2)) {
1103 nDSBPlus++;
1104 }
1105 }
1106 else if (count1 && count2) {
1107 nDSBPlus++;
1108 }
1109 else {
1110 nSSB++;
1111 }
1112 SSB2 = true;
1113 lastTenStrand2++;
1114 lastTenTracked2++;
1115 }
1116 if (oldnDSB != (nDSB + nDSBPlus)) {
1117 // we have had a DSB, time to classify its source.
1118 // DSB
1119 if (lastTenDirectStrand1 && lastTenDirectStrand2) {
1120 nDSBd = true;
1121 if (lastTenIndirectStrand1 || lastTenIndirectStrand2) nDSBm = true;
1122 }
1123 if (lastTenIndirectStrand1 && lastTenIndirectStrand2) {
1124 nDSBi = true;
1125 if (lastTenDirectStrand1 || lastTenDirectStrand2) nDSBh = true;
1126 if (lastTenDirectStrand1 && lastTenDirectStrand2) nDSBm = true;
1127 }
1128 if ((lastTenDirectStrand1 && lastTenIndirectStrand2)
1129 || (lastTenIndirectStrand1 && lastTenDirectStrand2))
1130 {
1131 nDSBh = true;
1132 }
1133 }
1134 }
1135
1136 // Return based on order of severity
1137
1138 // G4String complexity = "None";
1139
1140 complexityEnum complexity = NoneComplexity;
1141 if (nDSB > 1) {
1142 complexity = DSBplusplus;
1143 }
1144 else if (nDSBPlus != 0) {
1145 complexity = DSBplus;
1146 }
1147 else if (nDSB != 0) {
1148 complexity = DSB;
1149 }
1150 else if ((nSSB != 0) && SSB1 && SSB2) {
1151 complexity = twoSSB;
1152 }
1153 else if (nSSB > 1) {
1154 complexity = SSBplus;
1155 }
1156 else if (nSSB != 0) {
1157 complexity = SSB;
1158 }
1159 else {
1160 complexity = NoneComplexity;
1161 }
1162 G4bool direct = (strand1Direct || strand2Direct);
1163 G4bool indirect = (strand1Indirect || strand2Indirect);
1164 sourceEnum source = undefined;
1165 if (nDSB == 0) {
1166 if (direct && indirect) {
1167 source = SSBm;
1168 }
1169 else if (direct) {
1170 source = SSBd;
1171 }
1172 else if (indirect) {
1173 source = SSBi;
1174 }
1175 else {
1176 source = undefined;
1177 }
1178 }
1179 else {
1180 if (nDSBm != 0) {
1181 source = DSBm;
1182 }
1183 else if ((nDSBd != 0) && (nDSBi != 0)) {
1184 source = DSBm;
1185 }
1186 else if ((nDSBd != 0) && (nDSBh != 0)) {
1187 source = DSBm;
1188 }
1189 else if (nDSBd != 0) {
1190 source = DSBd;
1191 }
1192 else if (nDSBh != 0) {
1193 source = DSBh; // Catches DSBi && DSBh
1194 }
1195 else if (nDSBi != 0) {
1196 source = DSBi;
1197 }
1198 else {
1199 G4ExceptionDescription errmsg;
1200 this->PrintRecord("", dsbDistance);
1201 errmsg << "Error in source classification routine" << G4endl;
1202 errmsg << "I think this is " << complexity << "/"
1203 << "???" << G4endl;
1204 G4Exception("DamageRecord::GetClassification", "ERR_BAD_CLASSIFICATION", JustWarning, errmsg);
1205 source = undefined;
1206 }
1207 }
1208
1209 if (((complexity == NoneComplexity) && (source != undefined))
1210 || ((complexity != NoneComplexity) && (source == undefined)))
1211 {
1212 G4ExceptionDescription errmsg;
1213 errmsg << "You can't have a complexity without a source etc." << G4endl;
1214 errmsg << "I think this is " << complexity << "/" << source << G4endl;
1215 this->PrintRecord("", dsbDistance);
1216 G4Exception("DamageRecord::GetClassification", "ERR_BAD_CLASSIFICATION", JustWarning, errmsg);
1217 }
1218
1219 classification->fComplexity = complexity;
1220 classification->fSource = source;
1221 classification->fbaseDmg = baseDamage;
1222 classification->fStrandDmg = strandDamage;
1223
1224 return classification;
1225 }
1226
1227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1228
1229 // add a contrived test record to the damage record
1230 // params: s1 -> strand 1 damage code
1231 // b1 -> base 1 damage code
1232 // b2 -> base 2 damage code
1233 // s2 -> strand 2 damage code
1234 // codes:
1235 // param < 0: Nothing; param == 0 indirect damage;
1236 // param == 1 direct hit, no damage; param > 1: Direct hit physical damage
1237 void DamageRecord::AddTestDamage(G4int s1, G4int b1, G4int b2, G4int s2)
1238 {
1239 auto* bp = new BasePairDamageRecord;
1240 if (s1 == 0) {
1241 bp->fStrand1IndirectDmg = true;
1242 }
1243 if (s1 > 0) {
1244 bp->fStrand1Energy = 10 * eV;
1245 }
1246 if (s1 > 1) {
1247 bp->fStrand1DirectDmg = true;
1248 }
1249
1250 if (b1 == 0) {
1251 bp->fBp1IndirectDmg = true;
1252 }
1253 if (b1 > 0) {
1254 bp->fBp1Energy = 10 * eV;
1255 }
1256 if (b1 > 1) {
1257 bp->fbp1DirectDmg = true;
1258 }
1259
1260 if (s2 == 0) {
1261 bp->fStrand2IndirectDmg = true;
1262 }
1263 if (s2 > 0) {
1264 bp->fStrand2Energy = 10 * eV;
1265 }
1266 if (s2 > 1) {
1267 bp->fStrand2DirectDmg = true;
1268 }
1269
1270 if (b2 == 0) {
1271 bp->fBp2IndirectDmg = true;
1272 }
1273 if (b2 > 0) {
1274 bp->fBp2Energy = 10 * eV;
1275 }
1276 if (b2 > 1) {
1277 bp->fbp2DirectDmg = true;
1278 }
1279
1280 this->AddBasePairDamage(bp, G4ThreeVector(1, 1, G4UniformRand()));
1281 }
1282
1283 // Unit test for classification
1284 void AnalysisManager::TestClassification()
1285 {
1286 G4cout << G4endl;
1287 G4cout << "=================================" << G4endl;
1288 G4cout << "Classification Unit Test Begins" << G4endl;
1289
1290 // None
1291 auto* theRecord = new DamageRecord("Test", 0, 0, 0);
1292 theRecord->AddEmptyBPDamage(10);
1293 theRecord->AddTestDamage(1, -1, -1, 1);
1294 theRecord->AddEmptyBPDamage(5);
1295 theRecord->AddTestDamage(-1, -1, -1, 1);
1296 theRecord->AddEmptyBPDamage(5);
1297 theRecord->AddTestDamage(-1, -1, -1, 1);
1298 theRecord->AddEmptyBPDamage(10);
1299 DamageClassification* classification = theRecord->GetClassification();
1300 G4cout << "Test None" << G4endl;
1301 theRecord->PrintRecord("", 10);
1302 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1303 << G4endl;
1304 assert(classification->fComplexity == NoneComplexity);
1305 assert(classification->fSource == sourceEnum::undefined);
1306 delete classification;
1307 delete theRecord;
1308 G4cout << "---------------------------------" << G4endl;
1309
1310 // None, base damage
1311 theRecord = new DamageRecord("Test", 0, 0, 0);
1312 theRecord->AddEmptyBPDamage(10);
1313 theRecord->AddTestDamage(1, -1, 0, 1);
1314 theRecord->AddEmptyBPDamage(5);
1315 theRecord->AddTestDamage(-1, -1, -1, 1);
1316 theRecord->AddEmptyBPDamage(5);
1317 theRecord->AddTestDamage(-1, 0, -1, -1);
1318 theRecord->AddEmptyBPDamage(10);
1319 classification = theRecord->GetClassification();
1320 G4cout << "Test None" << G4endl;
1321 theRecord->PrintRecord("", 10);
1322 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1323 << G4endl;
1324 G4cout << "indirect Breaks = " << classification->fIndirectBreaks << G4endl;
1325 G4cout << "base Damage = " << classification->fbaseDmg << G4endl;
1326 assert(classification->fComplexity == NoneComplexity);
1327 assert(classification->fSource == sourceEnum::undefined);
1328 assert(classification->fIndirectBreaks == 0);
1329 assert(classification->fDirectBreaks == 0);
1330 assert(classification->fbaseDmg == 2);
1331 assert(classification->fStrandDmg == 0);
1332 delete classification;
1333 delete theRecord;
1334 G4cout << "---------------------------------" << G4endl;
1335
1336 // SSB
1337 theRecord = new DamageRecord("Test", 0, 0, 0);
1338 theRecord->AddEmptyBPDamage(10);
1339 theRecord->AddTestDamage(2, -1, -1, -1);
1340 theRecord->AddEmptyBPDamage(10);
1341 classification = theRecord->GetClassification();
1342 G4cout << "Test SSBd" << G4endl;
1343 theRecord->PrintRecord("", 10);
1344 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1345 << G4endl;
1346 assert(classification->fComplexity == SSB);
1347 assert(classification->fSource == SSBd);
1348 delete classification;
1349 delete theRecord;
1350 G4cout << "---------------------------------" << G4endl;
1351
1352 // SSBi
1353 theRecord = new DamageRecord("Test", 0, 0, 0);
1354 theRecord->AddEmptyBPDamage(10);
1355 theRecord->AddTestDamage(0, 0, -1, -1);
1356 theRecord->AddEmptyBPDamage(10);
1357 classification = theRecord->GetClassification();
1358 G4cout << "Test SSBi" << G4endl;
1359 theRecord->PrintRecord("", 10);
1360 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1361 << G4endl;
1362 assert(classification->fComplexity == SSB);
1363 assert(classification->fSource == SSBi);
1364 delete classification;
1365 delete theRecord;
1366 G4cout << "---------------------------------" << G4endl;
1367
1368 // SSBm / SSB+
1369 theRecord = new DamageRecord("Test", 0, 0, 0);
1370 theRecord->AddEmptyBPDamage(10);
1371 theRecord->AddTestDamage(0, -1, 0, -1);
1372 theRecord->AddEmptyBPDamage(6);
1373 theRecord->AddTestDamage(2, -1, 0, -1);
1374 theRecord->AddEmptyBPDamage(10);
1375 classification = theRecord->GetClassification();
1376 G4cout << "Test SSBm / SSB+" << G4endl;
1377 theRecord->PrintRecord("", 10);
1378 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1379 << G4endl;
1380 assert(classification->fComplexity == SSBplus);
1381 assert(classification->fSource == SSBm);
1382 delete classification;
1383 delete theRecord;
1384 G4cout << "---------------------------------" << G4endl;
1385
1386 // SSBm / 2SSB
1387 theRecord = new DamageRecord("Test", 0, 0, 0);
1388 theRecord->AddEmptyBPDamage(10);
1389 theRecord->AddTestDamage(2, -1, -1, -1);
1390 theRecord->AddEmptyBPDamage(16);
1391 theRecord->AddTestDamage(-1, -1, -1, 0);
1392 theRecord->AddEmptyBPDamage(10);
1393 classification = theRecord->GetClassification();
1394 G4cout << "Test SSBd / 2SSB" << G4endl;
1395 theRecord->PrintRecord("", 10);
1396 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1397 << G4endl;
1398 assert(classification->fComplexity == twoSSB);
1399 assert(classification->fSource == SSBm);
1400 delete classification;
1401 delete theRecord;
1402 G4cout << "---------------------------------" << G4endl;
1403
1404 // SSBi / 2SSB
1405 theRecord = new DamageRecord("Test", 0, 0, 0);
1406 theRecord->AddEmptyBPDamage(10);
1407 theRecord->AddTestDamage(0, -1, -1, -1);
1408 theRecord->AddEmptyBPDamage(16);
1409 theRecord->AddTestDamage(-1, -1, -1, 0);
1410 theRecord->AddEmptyBPDamage(16);
1411 theRecord->AddTestDamage(0, -1, -1, -1);
1412 theRecord->AddEmptyBPDamage(10);
1413 classification = theRecord->GetClassification();
1414 G4cout << "Test SSBi / 2SSB" << G4endl;
1415 theRecord->PrintRecord("", 10);
1416 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1417 << G4endl;
1418 assert(classification->fComplexity == twoSSB);
1419 assert(classification->fSource == SSBi);
1420 delete classification;
1421 delete theRecord;
1422 G4cout << "---------------------------------" << G4endl;
1423
1424 // SSBd / SSB+
1425 theRecord = new DamageRecord("Test", 0, 0, 0);
1426 theRecord->AddEmptyBPDamage(10);
1427 theRecord->AddTestDamage(2, -1, -1, -1);
1428 theRecord->AddEmptyBPDamage(6);
1429 theRecord->AddTestDamage(2, -1, 0, -1);
1430 theRecord->AddEmptyBPDamage(10);
1431 classification = theRecord->GetClassification();
1432 G4cout << "Test SSBd / SSB+" << G4endl;
1433 theRecord->PrintRecord("", 10);
1434 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1435 << G4endl;
1436 assert(classification->fComplexity == SSBplus);
1437 assert(classification->fSource == SSBd);
1438 delete classification;
1439 delete theRecord;
1440 G4cout << "---------------------------------" << G4endl;
1441
1442 // SSBm / SSB2
1443 theRecord = new DamageRecord("Test", 0, 0, 0);
1444 theRecord->AddEmptyBPDamage(10);
1445 theRecord->AddTestDamage(0, -1, -1, -1);
1446 theRecord->AddEmptyBPDamage(16);
1447 theRecord->AddTestDamage(-1, -1, -1, 2);
1448 theRecord->AddEmptyBPDamage(10);
1449 classification = theRecord->GetClassification();
1450 G4cout << "Test SSBm / 2SSB" << G4endl;
1451 theRecord->PrintRecord("", 10);
1452 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1453 << G4endl;
1454 assert(classification->fComplexity == twoSSB);
1455 assert(classification->fSource == SSBm);
1456 delete classification;
1457 delete theRecord;
1458 G4cout << "---------------------------------" << G4endl;
1459
1460 // SSBd / SSB2
1461 theRecord = new DamageRecord("Test", 0, 0, 0);
1462 theRecord->AddEmptyBPDamage(10);
1463 theRecord->AddTestDamage(2, -1, -1, -1);
1464 theRecord->AddEmptyBPDamage(16);
1465 theRecord->AddTestDamage(-1, -1, -1, 2);
1466 theRecord->AddEmptyBPDamage(10);
1467 classification = theRecord->GetClassification();
1468 G4cout << "Test SSBd / 2SSB" << G4endl;
1469 theRecord->PrintRecord("", 10);
1470 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1471 << G4endl;
1472 assert(classification->fComplexity == twoSSB);
1473 assert(classification->fSource == SSBd);
1474 delete classification;
1475 delete theRecord;
1476 G4cout << "---------------------------------" << G4endl;
1477
1478 // SSBd / SSB2
1479 theRecord = new DamageRecord("Test", 0, 0, 0);
1480 theRecord->AddEmptyBPDamage(10);
1481 theRecord->AddTestDamage(2, -1, -1, -1);
1482 theRecord->AddEmptyBPDamage(16);
1483 theRecord->AddTestDamage(-1, -1, -1, 2);
1484 theRecord->AddEmptyBPDamage(10);
1485 classification = theRecord->GetClassification();
1486 G4cout << "Test SSBd / 2SSB" << G4endl;
1487 theRecord->PrintRecord("", 10);
1488 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1489 << G4endl;
1490 assert(classification->fComplexity == twoSSB);
1491 assert(classification->fSource == SSBd);
1492 delete classification;
1493 delete theRecord;
1494 G4cout << "---------------------------------" << G4endl;
1495
1496 // DSBd / DSB
1497 theRecord = new DamageRecord("Test", 0, 0, 0);
1498 theRecord->AddEmptyBPDamage(10);
1499 theRecord->AddTestDamage(2, -1, -1, -1);
1500 theRecord->AddEmptyBPDamage(6);
1501 theRecord->AddTestDamage(-1, -1, -1, 2);
1502 theRecord->AddEmptyBPDamage(10);
1503 classification = theRecord->GetClassification();
1504 G4cout << "Test DSBd/DSB" << G4endl;
1505 theRecord->PrintRecord("", 10);
1506 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1507 << G4endl;
1508 assert(classification->fComplexity == DSB);
1509 assert(classification->fSource == DSBd);
1510 delete classification;
1511 delete theRecord;
1512 G4cout << "---------------------------------" << G4endl;
1513
1514 G4cout << G4endl;
1515
1516 // Induced Breaks DSB
1517 theRecord = new DamageRecord("Test", 0, 0, 0);
1518 theRecord->AddEmptyBPDamage(10);
1519
1520 auto* bp = new BasePairDamageRecord();
1521 bp->fBp1IndirectDmg = true;
1522 bp->fbp1InducedBreak = true;
1523 theRecord->AddBasePairDamage(bp, G4ThreeVector(0, 0, 0));
1524
1525 bp = new BasePairDamageRecord();
1526 bp->fBp2IndirectDmg = true;
1527 bp->fbp2InducedBreak = true;
1528 theRecord->AddBasePairDamage(bp, G4ThreeVector(0, 0, 0));
1529
1530 theRecord->AddEmptyBPDamage(10);
1531 classification = theRecord->GetClassification();
1532 G4cout << "Test DSBi/DSB induced" << G4endl;
1533 theRecord->PrintRecord("", 10);
1534 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1535 << G4endl;
1536 assert(classification->fComplexity == DSB);
1537 assert(classification->fSource == DSBi);
1538 delete classification;
1539 delete theRecord;
1540 G4cout << "---------------------------------" << G4endl;
1541
1542 // Induced Breaks SSB
1543 theRecord = new DamageRecord("Test", 0, 0, 0);
1544 theRecord->AddEmptyBPDamage(10);
1545
1546 bp = new BasePairDamageRecord();
1547 bp->fBp2IndirectDmg = true;
1548 bp->fbp2InducedBreak = true;
1549 theRecord->AddBasePairDamage(bp, G4ThreeVector(0, 0, 0));
1550
1551 theRecord->AddEmptyBPDamage(10);
1552 classification = theRecord->GetClassification();
1553 G4cout << "Test SSBi/SSB induced" << G4endl;
1554 theRecord->PrintRecord("", 10);
1555 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1556 << G4endl;
1557 assert(classification->fComplexity == SSB);
1558 assert(classification->fSource == SSBi);
1559 delete classification;
1560 delete theRecord;
1561 G4cout << "---------------------------------" << G4endl;
1562
1563 // DSBi / DSB
1564 theRecord = new DamageRecord("Test", 0, 0, 0);
1565 theRecord->AddEmptyBPDamage(10);
1566 theRecord->AddTestDamage(0, -1, -1, -1);
1567 theRecord->AddEmptyBPDamage(9);
1568 theRecord->AddTestDamage(-1, -1, -1, 0);
1569 theRecord->AddEmptyBPDamage(10);
1570 classification = theRecord->GetClassification();
1571 G4cout << "Test DSBi/DSB" << G4endl;
1572 theRecord->PrintRecord("", 10);
1573 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1574 << G4endl;
1575 assert(classification->fComplexity == DSB);
1576 assert(classification->fSource == DSBi);
1577 delete classification;
1578 delete theRecord;
1579 G4cout << "---------------------------------" << G4endl;
1580
1581 // DSBd / DSB
1582 theRecord = new DamageRecord("Test", 0, 0, 0);
1583 theRecord->AddEmptyBPDamage(10);
1584 theRecord->AddTestDamage(2, -1, -1, 2);
1585 theRecord->AddEmptyBPDamage(10);
1586 classification = theRecord->GetClassification();
1587 G4cout << "Test DSBd/DSB" << G4endl;
1588 theRecord->PrintRecord("", 10);
1589 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1590 << G4endl;
1591 assert(classification->fComplexity == DSB);
1592 assert(classification->fSource == DSBd);
1593 delete classification;
1594 delete theRecord;
1595 G4cout << "---------------------------------" << G4endl;
1596
1597 // DSBi/ DSB
1598 theRecord = new DamageRecord("Test", 0, 0, 0);
1599 theRecord->AddEmptyBPDamage(10);
1600 theRecord->AddTestDamage(0, -1, -1, 0);
1601 theRecord->AddEmptyBPDamage(10);
1602 classification = theRecord->GetClassification();
1603 G4cout << "Test DSBd/DSB" << G4endl;
1604 theRecord->PrintRecord("", 10);
1605 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1606 << G4endl;
1607 assert(classification->fComplexity == DSB);
1608 assert(classification->fSource == DSBi);
1609 delete classification;
1610 delete theRecord;
1611 G4cout << "---------------------------------" << G4endl;
1612
1613 // DSBh / DSB
1614 theRecord = new DamageRecord("Test", 0, 0, 0);
1615 theRecord->AddEmptyBPDamage(10);
1616 theRecord->AddTestDamage(0, -1, -1, 2);
1617 theRecord->AddEmptyBPDamage(10);
1618 classification = theRecord->GetClassification();
1619 G4cout << "Test DSBd/DSB" << G4endl;
1620 theRecord->PrintRecord("", 10);
1621 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1622 << G4endl;
1623 assert(classification->fComplexity == DSB);
1624 assert(classification->fSource == DSBh);
1625 delete classification;
1626 delete theRecord;
1627 G4cout << "---------------------------------" << G4endl;
1628
1629 // DSBh / DSB+
1630 theRecord = new DamageRecord("Test", 0, 0, 0);
1631 theRecord->AddEmptyBPDamage(10);
1632 theRecord->AddTestDamage(0, -1, -1, 2);
1633 theRecord->AddTestDamage(-1, -1, -1, 2);
1634 theRecord->AddEmptyBPDamage(10);
1635 classification = theRecord->GetClassification();
1636 G4cout << "Test DSBd/DSB+" << G4endl;
1637 theRecord->PrintRecord("", 10);
1638 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1639 << G4endl;
1640 assert(classification->fComplexity == DSBplus);
1641 assert(classification->fSource == DSBh);
1642 delete classification;
1643 delete theRecord;
1644 G4cout << "---------------------------------" << G4endl;
1645
1646 // DSBm / DSB+
1647 theRecord = new DamageRecord("Test", 0, 0, 0);
1648 theRecord->AddEmptyBPDamage(10);
1649 theRecord->AddTestDamage(0, -1, -1, 2);
1650 theRecord->AddTestDamage(2, -1, -1, -1);
1651 theRecord->AddEmptyBPDamage(10);
1652 classification = theRecord->GetClassification();
1653 G4cout << "Test DSBd/DSB+" << G4endl;
1654 theRecord->PrintRecord("", 10);
1655 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1656 << G4endl;
1657 assert(classification->fComplexity == DSBplus);
1658 assert(classification->fSource == DSBm);
1659 delete classification;
1660 delete theRecord;
1661 G4cout << "---------------------------------" << G4endl;
1662
1663 // DSBd / DSB+
1664 theRecord = new DamageRecord("Test", 0, 0, 0);
1665 theRecord->AddEmptyBPDamage(10);
1666 theRecord->AddTestDamage(2, -1, -1, -1);
1667 theRecord->AddEmptyBPDamage(6);
1668 theRecord->AddTestDamage(-1, -1, -1, 2);
1669 theRecord->AddTestDamage(-1, -1, -1, 2);
1670 theRecord->AddEmptyBPDamage(10);
1671 classification = theRecord->GetClassification();
1672 G4cout << "Test DSBd/DSB+" << G4endl;
1673 theRecord->PrintRecord("", 10);
1674 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1675 << G4endl;
1676 assert(classification->fComplexity == DSBplus);
1677 assert(classification->fSource == DSBd);
1678 delete classification;
1679 delete theRecord;
1680 G4cout << "---------------------------------" << G4endl;
1681
1682 // DSBd / DSB+
1683 theRecord = new DamageRecord("Test", 0, 0, 0);
1684 theRecord->AddEmptyBPDamage(10);
1685 theRecord->AddTestDamage(2, -1, -1, -1);
1686 theRecord->AddEmptyBPDamage(6);
1687 theRecord->AddTestDamage(-1, -1, -1, 2);
1688 theRecord->AddTestDamage(-1, -1, -1, 2);
1689 theRecord->AddEmptyBPDamage(6);
1690 theRecord->AddTestDamage(-1, -1, -1, 2);
1691 theRecord->AddEmptyBPDamage(10);
1692 classification = theRecord->GetClassification();
1693 G4cout << "Test DSBd/DSB+" << G4endl;
1694 theRecord->PrintRecord("", 10);
1695 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1696 << G4endl;
1697 assert(classification->fComplexity == DSBplus);
1698 assert(classification->fSource == DSBd);
1699 delete classification;
1700 delete theRecord;
1701 G4cout << "---------------------------------" << G4endl;
1702
1703 // DSBi / DSB+
1704 theRecord = new DamageRecord("Test", 0, 0, 0);
1705 theRecord->AddEmptyBPDamage(10);
1706 theRecord->AddTestDamage(0, -1, -1, -1);
1707 theRecord->AddEmptyBPDamage(6);
1708 theRecord->AddTestDamage(-1, -1, -1, 0);
1709 theRecord->AddTestDamage(-1, -1, -1, 0);
1710 theRecord->AddEmptyBPDamage(6);
1711 theRecord->AddTestDamage(-1, -1, -1, 2);
1712 theRecord->AddEmptyBPDamage(10);
1713 classification = theRecord->GetClassification();
1714 G4cout << "Test DSBd/DSB+" << G4endl;
1715 theRecord->PrintRecord("", 10);
1716 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1717 << G4endl;
1718 assert(classification->fComplexity == DSBplus);
1719 assert(classification->fSource == DSBi);
1720 delete classification;
1721 delete theRecord;
1722 G4cout << "---------------------------------" << G4endl;
1723
1724 // DSBh / DSB++
1725 theRecord = new DamageRecord("Test", 0, 0, 0);
1726 theRecord->AddEmptyBPDamage(10);
1727 theRecord->AddTestDamage(0, -1, -1, -1);
1728 theRecord->AddEmptyBPDamage(6);
1729 theRecord->AddTestDamage(-1, -1, -1, 0);
1730 theRecord->AddTestDamage(-1, -1, -1, 0);
1731 theRecord->AddEmptyBPDamage(6);
1732 theRecord->AddTestDamage(2, -1, -1, -1);
1733 theRecord->AddEmptyBPDamage(10);
1734 classification = theRecord->GetClassification();
1735 G4cout << "Test DSBd/DSB++" << G4endl;
1736 theRecord->PrintRecord("", 10);
1737 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1738 << G4endl;
1739 assert(classification->fComplexity == DSBplusplus);
1740 assert(classification->fSource == DSBh);
1741 delete classification;
1742 delete theRecord;
1743 G4cout << "---------------------------------" << G4endl;
1744
1745 // DSBd / DSB
1746 theRecord = new DamageRecord("Test", 0, 0, 0);
1747 theRecord->AddEmptyBPDamage(10);
1748 theRecord->AddTestDamage(2, -1, -1, -1);
1749 theRecord->AddEmptyBPDamage(6);
1750 theRecord->AddTestDamage(-1, -1, -1, 2);
1751 theRecord->AddEmptyBPDamage(6);
1752 theRecord->AddTestDamage(-1, -1, -1, 2);
1753 theRecord->AddEmptyBPDamage(10);
1754 classification = theRecord->GetClassification();
1755 G4cout << "Test DSBd/DSB" << G4endl;
1756 theRecord->PrintRecord("", 10);
1757 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1758 << G4endl;
1759 assert(classification->fComplexity == DSB);
1760 assert(classification->fSource == DSBd);
1761 delete classification;
1762 delete theRecord;
1763 G4cout << "---------------------------------" << G4endl;
1764
1765 // DSBd / DSB++
1766 theRecord = new DamageRecord("Test", 0, 0, 0);
1767 theRecord->AddEmptyBPDamage(10);
1768 theRecord->AddTestDamage(2, -1, -1, -1);
1769 theRecord->AddEmptyBPDamage(6);
1770 theRecord->AddTestDamage(-1, -1, -1, 2);
1771 theRecord->AddTestDamage(-1, -1, -1, 2);
1772 theRecord->AddEmptyBPDamage(6);
1773 theRecord->AddTestDamage(2, -1, -1, -1);
1774 theRecord->AddEmptyBPDamage(10);
1775 classification = theRecord->GetClassification();
1776 G4cout << "Test DSBd/DSB++" << G4endl;
1777 theRecord->PrintRecord("", 10);
1778 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1779 << G4endl;
1780 assert(classification->fComplexity == DSBplusplus);
1781 assert(classification->fSource == DSBd);
1782 delete classification;
1783 delete theRecord;
1784 G4cout << "---------------------------------" << G4endl;
1785
1786 // DSBh / DSB++
1787 theRecord = new DamageRecord("Test", 0, 0, 0);
1788 theRecord->AddEmptyBPDamage(10);
1789 theRecord->AddTestDamage(2, -1, -1, -1);
1790 theRecord->AddEmptyBPDamage(6);
1791 theRecord->AddTestDamage(-1, -1, -1, 0);
1792 theRecord->AddTestDamage(-1, -1, -1, 0);
1793 theRecord->AddEmptyBPDamage(6);
1794 theRecord->AddTestDamage(2, -1, -1, -1);
1795 theRecord->AddEmptyBPDamage(10);
1796 classification = theRecord->GetClassification();
1797 G4cout << "Test DSBh/DSB++" << G4endl;
1798 theRecord->PrintRecord("", 10);
1799 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1800 << G4endl;
1801 assert(classification->fComplexity == DSBplusplus);
1802 assert(classification->fSource == DSBh);
1803 delete classification;
1804 delete theRecord;
1805 G4cout << "---------------------------------" << G4endl;
1806
1807 // DSBm / DSB++
1808 theRecord = new DamageRecord("Test", 0, 0, 0);
1809 theRecord->AddEmptyBPDamage(10);
1810 theRecord->AddTestDamage(2, -1, -1, -1);
1811 theRecord->AddEmptyBPDamage(6);
1812 theRecord->AddTestDamage(-1, -1, -1, 0);
1813 theRecord->AddTestDamage(-1, -1, -1, 2);
1814 theRecord->AddEmptyBPDamage(6);
1815 theRecord->AddTestDamage(0, -1, -1, -1);
1816 theRecord->AddEmptyBPDamage(10);
1817 classification = theRecord->GetClassification();
1818 G4cout << "Test DSBm/DSB++" << G4endl;
1819 theRecord->PrintRecord("", 10);
1820 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1821 << G4endl;
1822 assert(classification->fComplexity == DSBplusplus);
1823 assert(classification->fSource == DSBm);
1824 delete classification;
1825 delete theRecord;
1826 G4cout << "---------------------------------" << G4endl;
1827
1828 // DSBi / DSB+
1829 theRecord = new DamageRecord("Test", 0, 0, 0);
1830 theRecord->AddEmptyBPDamage(10);
1831 theRecord->AddTestDamage(0, -1, -1, -1);
1832 theRecord->AddEmptyBPDamage(6);
1833 theRecord->AddTestDamage(0, -1, -1, -1);
1834 theRecord->AddTestDamage(-1, -1, -1, 0);
1835 theRecord->AddEmptyBPDamage(6);
1836 theRecord->AddTestDamage(-1, -1, -1, 0);
1837 theRecord->AddEmptyBPDamage(10);
1838 classification = theRecord->GetClassification();
1839 G4cout << "Test DSBi/DSB++" << G4endl;
1840 theRecord->PrintRecord("", 10);
1841 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1842 << G4endl;
1843 assert(classification->fComplexity == DSBplusplus);
1844 assert(classification->fSource == DSBi);
1845 delete classification;
1846 delete theRecord;
1847 G4cout << "---------------------------------" << G4endl;
1848
1849 // DSBi / DSB++
1850 theRecord = new DamageRecord("Test", 0, 0, 0);
1851 theRecord->AddEmptyBPDamage(10);
1852 theRecord->AddTestDamage(0, -1, -1, -1);
1853 theRecord->AddEmptyBPDamage(6);
1854 theRecord->AddTestDamage(-1, -1, -1, 0);
1855 theRecord->AddTestDamage(0, -1, -1, -1);
1856 theRecord->AddEmptyBPDamage(6);
1857 theRecord->AddTestDamage(-1, -1, -1, 0);
1858 theRecord->AddEmptyBPDamage(10);
1859 classification = theRecord->GetClassification();
1860 G4cout << "Test DSBi/DSB++" << G4endl;
1861 theRecord->PrintRecord("", 10);
1862 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1863 << G4endl;
1864 assert(classification->fComplexity == DSBplusplus);
1865 assert(classification->fSource == DSBi);
1866 delete classification;
1867 delete theRecord;
1868 G4cout << "---------------------------------" << G4endl;
1869
1870 // DSBd / DSB++
1871 theRecord = new DamageRecord("Test", 0, 0, 0);
1872 theRecord->AddEmptyBPDamage(10);
1873 theRecord->AddTestDamage(2, -1, -1, -1);
1874 theRecord->AddEmptyBPDamage(6);
1875 theRecord->AddTestDamage(-1, -1, -1, 2);
1876 theRecord->AddEmptyBPDamage(15);
1877 theRecord->AddTestDamage(-1, -1, -1, 2);
1878 theRecord->AddEmptyBPDamage(6);
1879 theRecord->AddTestDamage(2, -1, -1, -1);
1880 theRecord->AddEmptyBPDamage(10);
1881 classification = theRecord->GetClassification();
1882 G4cout << "Test DSBd/DSB++" << G4endl;
1883 theRecord->PrintRecord("", 10);
1884 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1885 << G4endl;
1886 assert(classification->fComplexity == DSBplusplus);
1887 assert(classification->fSource == DSBd);
1888 delete classification;
1889 delete theRecord;
1890 G4cout << "---------------------------------" << G4endl;
1891
1892 // DSBi / DSB++
1893 theRecord = new DamageRecord("Test", 0, 0, 0);
1894 theRecord->AddEmptyBPDamage(10);
1895 theRecord->AddTestDamage(0, -1, -1, -1);
1896 theRecord->AddEmptyBPDamage(6);
1897 theRecord->AddTestDamage(-1, -1, -1, 0);
1898 theRecord->AddEmptyBPDamage(12);
1899 theRecord->AddTestDamage(-1, -1, -1, 2);
1900 theRecord->AddEmptyBPDamage(12);
1901 theRecord->AddTestDamage(-1, -1, -1, 0);
1902 theRecord->AddEmptyBPDamage(6);
1903 theRecord->AddTestDamage(0, -1, -1, -1);
1904 theRecord->AddEmptyBPDamage(10);
1905 classification = theRecord->GetClassification();
1906 G4cout << "Test DSBi/DSB++" << G4endl;
1907 theRecord->PrintRecord("", 10);
1908 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1909 << G4endl;
1910 assert(classification->fComplexity == DSBplusplus);
1911 assert(classification->fSource == DSBi);
1912 delete classification;
1913 delete theRecord;
1914 G4cout << "---------------------------------" << G4endl;
1915
1916 // DSBh / DSB++
1917 theRecord = new DamageRecord("Test", 0, 0, 0);
1918 theRecord->AddEmptyBPDamage(10);
1919 theRecord->AddTestDamage(0, -1, -1, -1);
1920 theRecord->AddTestDamage(0, -1, -1, 0);
1921 theRecord->AddTestDamage(0, -1, -1, 0);
1922 theRecord->AddEmptyBPDamage(6);
1923 theRecord->AddTestDamage(-1, -1, -1, 2);
1924 theRecord->AddEmptyBPDamage(12);
1925 theRecord->AddTestDamage(-1, -1, -1, 2);
1926 theRecord->AddEmptyBPDamage(12);
1927 theRecord->AddTestDamage(-1, -1, -1, 0);
1928 theRecord->AddEmptyBPDamage(6);
1929 theRecord->AddTestDamage(0, -1, -1, -1);
1930 theRecord->AddEmptyBPDamage(10);
1931 classification = theRecord->GetClassification();
1932 G4cout << "Test DSBh/DSB++" << G4endl;
1933 theRecord->PrintRecord("", 10);
1934 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1935 << G4endl;
1936 assert(classification->fComplexity == DSBplusplus);
1937 assert(classification->fSource == DSBh);
1938 delete classification;
1939 delete theRecord;
1940 G4cout << "---------------------------------" << G4endl;
1941
1942 // DSBm / DSB++
1943 theRecord = new DamageRecord("Test", 0, 0, 0);
1944 theRecord->AddEmptyBPDamage(10);
1945 theRecord->AddTestDamage(0, -1, -1, -1);
1946 theRecord->AddTestDamage(2, -1, -1, 0);
1947 theRecord->AddTestDamage(0, -1, -1, 0);
1948 theRecord->AddEmptyBPDamage(6);
1949 theRecord->AddTestDamage(-1, -1, -1, 2);
1950 theRecord->AddEmptyBPDamage(12);
1951 theRecord->AddTestDamage(-1, -1, -1, 2);
1952 theRecord->AddEmptyBPDamage(12);
1953 theRecord->AddTestDamage(-1, -1, -1, 0);
1954 theRecord->AddEmptyBPDamage(6);
1955 theRecord->AddTestDamage(0, -1, -1, -1);
1956 theRecord->AddEmptyBPDamage(10);
1957 classification = theRecord->GetClassification();
1958 G4cout << "Test DSBm/DSB++" << G4endl;
1959 theRecord->PrintRecord("", 10);
1960 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1961 << G4endl;
1962 assert(classification->fComplexity == DSBplusplus);
1963 assert(classification->fSource == DSBm);
1964 delete classification;
1965 delete theRecord;
1966 G4cout << "---------------------------------" << G4endl;
1967
1968 // DSBi / DSB++
1969 theRecord = new DamageRecord("Test", 0, 0, 0);
1970 theRecord->AddEmptyBPDamage(10);
1971 theRecord->AddTestDamage(-1, 0, -1, -1);
1972 theRecord->AddTestDamage(-1, -1, -1, 0);
1973 theRecord->AddTestDamage(-1, -1, -1, -1);
1974 theRecord->AddTestDamage(-1, -1, -1, 0);
1975 theRecord->AddTestDamage(-1, -1, -1, -1);
1976 theRecord->AddTestDamage(-1, -1, -1, -1);
1977 theRecord->AddTestDamage(0, -1, -1, -1);
1978 theRecord->AddEmptyBPDamage(10);
1979 classification = theRecord->GetClassification();
1980 G4cout << "Test DSBi/DSB+" << G4endl;
1981 theRecord->PrintRecord("", 10);
1982 G4cout << "Classification: " << classification->fSource << "/" << classification->fComplexity
1983 << G4endl;
1984 assert(classification->fComplexity == DSBplus);
1985 assert(classification->fSource == DSBi);
1986 delete classification;
1987 delete theRecord;
1988 G4cout << "---------------------------------" << G4endl;
1989
1990 G4cout << "Classification Unit Test Complete" << G4endl;
1991 G4cout << "=================================" << G4endl;
1992
1993 G4cout << G4endl;
1994 }
1995
1996 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1997
1998 void DamageRecord::AddEmptyBPDamage(int64_t ii) // dousatsu
1999 {
2000 auto basePairNumber = ii;
2001 while (basePairNumber > 0) {
2002 fDamageRecords.push_back(new BasePairDamageRecord);
2003 basePairNumber--;
2004 }
2005 }
2006
2007 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2008
2009 void DamageRecord::AddStrandHit(const G4MoleculeDefinition* mol)
2010 {
2011 if (mol == G4OH::Definition()) {
2012 fOHStrand++;
2013 }
2014 else if (mol == G4Electron_aq::Definition()) {
2015 fEaqStrand++;
2016 }
2017 else if (mol == G4Hydrogen::Definition()) {
2018 fHStrand++;
2019 }
2020 }
2021
2022 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2023
2024 void DamageRecord::AddBaseHit(const G4MoleculeDefinition* mol)
2025 {
2026 // G4cout << "Hit by " << mol->GetParticleName() << G4endl;
2027 if (mol == fOH) {
2028 fOHBase++;
2029 }
2030 else if (mol == fe_aq) {
2031 fEaqBase++;
2032 }
2033 else if (mol == fH) {
2034 fHBase++;
2035 }
2036 }
2037
2038 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2039
2040 G4ThreeVector DamageRecord::GetMeanPosition() const
2041 {
2042 G4double x(0.);
2043 G4double y(0.);
2044 G4double z(0.);
2045 for (const auto& fPosition : fPositions) {
2046 x += fPosition.getX();
2047 y += fPosition.getY();
2048 z += fPosition.getZ();
2049 }
2050
2051 if (fPositions.empty()) {
2052 G4ExceptionDescription errmsg;
2053 errmsg << "fPositions is emply ";
2054 G4Exception("DamageRecord", "ERR_INVALID_PROB", FatalException, errmsg);
2055 return G4ThreeVector(0., 0., 0.);
2056 }
2057 else {
2058 G4double factor = 1.0 / fPositions.size();
2059 return G4ThreeVector(x * factor, y * factor, z * factor);
2060 }
2061 }
2062
2063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2064
2065 G4double DamageRecord::GetMeanDistance() const
2066 {
2067 G4double d = 0.;
2068 for (auto it = fPositions.begin(); it != fPositions.end(); ++it) {
2069 for (auto jt = it + 1; jt != fPositions.end(); ++jt) {
2070 d += ((*jt) - (*it)).mag();
2071 }
2072 }
2073 G4int number = fPositions.size();
2074 return d / number;
2075 }
2076
2077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2078
2079 G4double DamageRecord::GetEnergy() const
2080 {
2081 G4double en = 0;
2082 for (auto bp : fDamageRecords) {
2083 en = en + bp->fBp1Energy + bp->fBp2Energy + bp->fStrand1Energy + bp->fStrand2Energy;
2084 }
2085 return en;
2086 }