| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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........oo
51
52 AnalysisManager::AnalysisManager()
53 {
54 fAnalysisManager = G4AnalysisManager::Instan
55 fpAnalysisMessenger = new AnalysisMessenger(
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oo
59
60 AnalysisManager::~AnalysisManager()
61 {
62 delete fpAnalysisMessenger;
63 }
64
65 //....oooOO0OOooo........oooOO0OOooo........oo
66
67 void AnalysisManager::Initialize()
68 {
69 const G4Run* run = G4RunManager::GetRunManag
70 if (run->GetRunID() == 0) {
71 G4cout << "AnalysisManager::Initialize() G
72 fAnalysisManager->SetDefaultFileType("root
73 fAnalysisManager->SetVerboseLevel(1);
74
75 fAnalysisManager->SetNtupleDirectoryName("
76 fAnalysisManager->SetHistoDirectoryName("h
77
78 fAnalysisManager->CreateH1("ssb_counts", "
79 fAnalysisManager->CreateH1("ssb_energies_e
80 // fAnalysisManager->CreateH1("fragments",
81 // 500);//ORG
82 fAnalysisManager->CreateH1("fragments", "F
83 10000); // dou
84 fAnalysisManager->CreateH3("local_pos", "S
85 25, 50, -25, 25
86 fAnalysisManager->CreateH3("global_pos", "
87 -100, 100, 50,
88 fAnalysisManager->CreateH1("e_cell_kev", "
89 2500); // dous
90
91 fAnalysisManager->CreateNtuple("primary_so
92 fAnalysisManager->CreateNtupleSColumn("Pri
93 fAnalysisManager->CreateNtupleDColumn("Ene
94 fAnalysisManager->CreateNtupleDColumn("Pos
95 fAnalysisManager->CreateNtupleDColumn("Pos
96 fAnalysisManager->CreateNtupleDColumn("Pos
97 fAnalysisManager->CreateNtupleDColumn("Mom
98 fAnalysisManager->CreateNtupleDColumn("Mom
99 fAnalysisManager->CreateNtupleDColumn("Mom
100 fAnalysisManager->CreateNtupleDColumn("Sto
101 fAnalysisManager->CreateNtupleDColumn("Sto
102 fAnalysisManager->CreateNtupleDColumn("Sto
103 fAnalysisManager->CreateNtupleDColumn("Tra
104 fAnalysisManager->CreateNtupleDColumn("Tra
105 fAnalysisManager->FinishNtuple();
106
107 fAnalysisManager->CreateNtuple("chromosome
108 fAnalysisManager->CreateNtupleSColumn("chr
109 fAnalysisManager->CreateNtupleDColumn("e_c
110 fAnalysisManager->CreateNtupleDColumn("e_d
111 fAnalysisManager->FinishNtuple();
112
113 fAnalysisManager->CreateNtuple("classifica
114 fAnalysisManager->CreateNtupleSColumn("Pri
115 fAnalysisManager->CreateNtupleDColumn("Ene
116 fAnalysisManager->CreateNtupleIColumn("Non
117 fAnalysisManager->CreateNtupleIColumn("SSB
118 fAnalysisManager->CreateNtupleIColumn("SSB
119 fAnalysisManager->CreateNtupleIColumn("2SS
120 fAnalysisManager->CreateNtupleIColumn("DSB
121 fAnalysisManager->CreateNtupleIColumn("DSB
122 fAnalysisManager->CreateNtupleIColumn("DSB
123 fAnalysisManager->FinishNtuple();
124
125 fAnalysisManager->CreateNtuple("source", "
126 fAnalysisManager->CreateNtupleSColumn("Pri
127 fAnalysisManager->CreateNtupleDColumn("Ene
128 fAnalysisManager->CreateNtupleIColumn("Non
129 fAnalysisManager->CreateNtupleIColumn("SSB
130 fAnalysisManager->CreateNtupleIColumn("SSB
131 fAnalysisManager->CreateNtupleIColumn("SSB
132 fAnalysisManager->CreateNtupleIColumn("DSB
133 fAnalysisManager->CreateNtupleIColumn("DSB
134 fAnalysisManager->CreateNtupleIColumn("DSB
135 fAnalysisManager->CreateNtupleIColumn("DSB
136 fAnalysisManager->FinishNtuple();
137
138 fAnalysisManager->CreateNtuple("damage", "
139 fAnalysisManager->CreateNtupleIColumn("Eve
140 fAnalysisManager->CreateNtupleSColumn("Pri
141 fAnalysisManager->CreateNtupleDColumn("Ene
142 fAnalysisManager->CreateNtupleSColumn("Typ
143 fAnalysisManager->CreateNtupleSColumn("Sou
144 fAnalysisManager->CreateNtupleDColumn("Pos
145 fAnalysisManager->CreateNtupleDColumn("Pos
146 fAnalysisManager->CreateNtupleDColumn("Pos
147 fAnalysisManager->CreateNtupleDColumn("Siz
148 fAnalysisManager->CreateNtupleIColumn("Fra
149 fAnalysisManager->CreateNtupleIColumn("Bas
150 fAnalysisManager->CreateNtupleIColumn("Str
151 fAnalysisManager->CreateNtupleIColumn("Dir
152 fAnalysisManager->CreateNtupleIColumn("Ind
153 fAnalysisManager->CreateNtupleIColumn("Eaq
154 fAnalysisManager->CreateNtupleIColumn("Eaq
155 fAnalysisManager->CreateNtupleIColumn("OHB
156 fAnalysisManager->CreateNtupleIColumn("OHS
157 fAnalysisManager->CreateNtupleIColumn("HBa
158 fAnalysisManager->CreateNtupleIColumn("HSt
159 fAnalysisManager->CreateNtupleDColumn("Ene
160 fAnalysisManager->CreateNtupleIColumn("Ind
161 fAnalysisManager->CreateNtupleIColumn("Cha
162 fAnalysisManager->CreateNtupleIColumn("Str
163 fAnalysisManager->CreateNtupleDColumn("Bas
164 fAnalysisManager->CreateNtupleSColumn("Nam
165 fAnalysisManager->FinishNtuple();
166 }
167
168 if (!fFileName) {
169 fFileName = "molecular-dna";
170 }
171 fAnalysisManager->OpenFile(fFileName);
172 }
173
174 //....oooOO0OOooo........oooOO0OOooo........oo
175
176 void AnalysisManager::Close()
177 {
178 fAnalysisManager->Write();
179 fAnalysisManager->CloseFile();
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oo
183
184 void AnalysisManager::ProcessDNAHitsVector(con
185 {
186 // Check we have DNA Geometry (runs once)
187 if (fpDNAGeometry == nullptr) {
188 auto det = dynamic_cast<const DetectorCons
189 G4RunManager::GetRunManager()->GetUserDe
190 fpDNAGeometry = det->GetDNAGeometry();
191 }
192
193 if (dnaHits.empty()) {
194 return;
195 }
196 const G4Event* event = G4EventManager::GetEv
197 // SSBs
198 for (auto dnaHit : dnaHits) {
199 fAnalysisManager->FillH1(fAnalysisManager-
200 dnaHit->GetEnergy
201 if (fChainToSave == -1) {
202 // save all chains
203 fAnalysisManager->FillH3(
204 fAnalysisManager->GetH3Id("local_pos")
205 dnaHit->GetLocalPosition().getY() / nm
206 fAnalysisManager->FillH3(fAnalysisManage
207 dnaHit->GetPosi
208 dnaHit->GetPosi
209 }
210 else if (fChainToSave == dnaHit->GetChainI
211 // save only specified chain
212 fAnalysisManager->FillH3(
213 fAnalysisManager->GetH3Id("local_pos")
214 dnaHit->GetLocalPosition().getY() / nm
215 fAnalysisManager->FillH3(fAnalysisManage
216 dnaHit->GetPosi
217 dnaHit->GetPosi
218 }
219 }
220 fAnalysisManager->FillH1(fAnalysisManager->G
221
222 // DSBs
223 // Sorting
224 // Make a map relating chromosome/chain to i
225 std::map<std::pair<G4String, G4int>, BinaryT
226
227 // Populate the binary tree
228 // G4cout << "Building Binary Trees for " <<
229 // G4endl;
230 for (auto it = dnaHits.begin(); it != dnaHit
231 std::pair<G4String, G4int> key = std::make
232 if (treemap.find(key) == treemap.end()) {
233 treemap[key] = new BinaryTree();
234 if (fAnalysisManager->GetVerboseLevel()
235 G4cout << "Constructing binary tree. C
236 << G4endl;
237 }
238 }
239 treemap.at(key)->Insert(*it);
240 if (fAnalysisManager->GetVerboseLevel() >=
241 G4cout << "Adding hit to binary tree. Ch
242 << G4endl;
243 }
244 }
245
246 // Analysis
247 // Create a vector to stock damage records f
248 // G4cout << "Building Damage Records" << G4
249 std::vector<DamageRecord*> damageRecords;
250 for (auto& it : treemap) {
251 if (fAnalysisManager->GetVerboseLevel() >=
252 G4cout << "Analysing hits for chromosome
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() > 92233
262 G4cout << " SEE AnalysisManager !!!" <
263 abort();
264 } // dousatsu
265 // if (currentHit->GetBasePairIdx()>2147
266 // AnalysisManager !!!"<<G4endl;abort();
267 auto idx = (int64_t)currentHit->GetBaseP
268 int64_t startidx = (idx > 10) ? (idx - 1
269 // G4int idx = currentHit->GetBasePairId
270 // G4int startidx = (idx > 10) ? (idx -
271 G4String name = it.first.first + "_" + s
272 + std::to_string(startid
273 auto* record =
274 new DamageRecord(name, idx, currentHit
275 damageRecords.push_back(record);
276 BasePairDamageRecord* bp;
277 int64_t gap = 0; // dousatsu
278
279 // bookend with 10 bps on either side of
280 record->AddEmptyBPDamage(idx - startidx)
281 // continues until fragment ends
282 while (true) {
283 bp = new BasePairDamageRecord;
284 bp->fStrand1Energy = currentHit->GetSt
285 bp->fStrand2Energy = currentHit->GetSt
286 bp->fBp1Energy = currentHit->GetBP1Ene
287 bp->fBp2Energy = currentHit->GetBP2Ene
288 bp->fBp1IndirectEvt = (currentHit->Get
289 bp->fBp2IndirectEvt = (currentHit->Get
290 bp->fStrand1IndirectEvt = (currentHit-
291 bp->fStrand2IndirectEvt = (currentHit-
292 bp->fbp1DirectDmg = false; // No phys
293 bp->fbp2DirectDmg = false; // No phys
294
295 bp->fStrand1DirectDmg =
296 fpDNAGeometry->GetDamageModel()->IsD
297 bp->fStrand2DirectDmg =
298 fpDNAGeometry->GetDamageModel()->IsD
299
300 // if(bp->fStrand1DirectDmg) bp
301 // else if(bp->fStrand2DirectDm
302
303 if (bp->fBp1IndirectEvt) {
304 bp->fBp1IndirectDmg =
305 fpDNAGeometry->GetDamageModel()->I
306 if (bp->fBp1IndirectDmg) {
307 bp->fbp1InducedBreak =
308 fpDNAGeometry->GetDamageModel()-
309 }
310 }
311
312 if (bp->fBp2IndirectEvt) {
313 bp->fBp2IndirectDmg =
314 fpDNAGeometry->GetDamageModel()->I
315 if (bp->fBp2IndirectDmg) {
316 bp->fbp2InducedBreak =
317 fpDNAGeometry->GetDamageModel()-
318 }
319 }
320
321 if (bp->fStrand1IndirectEvt) {
322 bp->fStrand1IndirectDmg =
323 fpDNAGeometry->GetDamageModel()->I
324 }
325 if (bp->fStrand2IndirectEvt) {
326 bp->fStrand2IndirectDmg =
327 fpDNAGeometry->GetDamageModel()->I
328 }
329 // Record radical-base/strand damage.
330 // it is possible for base/strand dama
331 // direct hits due to the damage induc
332 // This mainly affects strands.
333 if (bp->fBp1IndirectDmg) {
334 record->AddBaseHit(currentHit->GetBa
335 }
336 if (bp->fBp2IndirectDmg) {
337 record->AddBaseHit(currentHit->GetBa
338 }
339 if (bp->fStrand1IndirectDmg) {
340 record->AddStrandHit(currentHit->Get
341 // strandID = 1;
342 }
343 if (bp->fStrand2IndirectDmg) {
344 record->AddStrandHit(currentHit->Get
345 // strandID = 2;
346 }
347
348 record->AddBasePairDamage(bp, currentH
349
350 nextHit = it.second->Next(currentHit);
351 if (nextHit == nullptr) {
352 break;
353 }
354
355 gap = nextHit->GetBasePairIdx() - curr
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, no
368 // ie. plasmids etc. Each placement
369 if (currentHit->GetPlacementIdx() !=
370 if (fpDNAGeometry->GetVerbosity()
371 G4cout << "Analysis passing to a
372 }
373 currentHit = nextHit;
374 break;
375 }
376 else {
377 record->AddEmptyBPDamage(gap);
378 }
379 }
380 else {
381 G4Exception("MolecularAnalaysisManag
382 "The value set in /analy
383 }
384 currentHit = nextHit;
385 }
386 record->AddEmptyBPDamage(10);
387 // end bookend
388 }
389 }
390 // Count the number of breaks and print out
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 = damageReco
413 breakmap[classif->fComplexity]++;
414 sourcemap[classif->fSource]++;
415 fAnalysisManager->FillH1(fAnalysisManager-
416
417 ///////////dousatsu-----------------------
418 fAnalysisManager->FillNtupleIColumn(4, 0,
419 fAnalysisManager->FillNtupleSColumn(
420 4, 1, event->GetPrimaryVertex()->GetPrim
421 fAnalysisManager->FillNtupleDColumn(
422 4, 2, event->GetPrimaryVertex()->GetPrim
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,
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,
479 fAnalysisManager->FillNtupleDColumn(4, 5,
480 fAnalysisManager->FillNtupleDColumn(4, 6,
481 fAnalysisManager->FillNtupleDColumn(4, 7,
482 fAnalysisManager->FillNtupleDColumn(4, 8,
483 fAnalysisManager->FillNtupleIColumn(4, 9,
484 fAnalysisManager->FillNtupleIColumn(4, 10,
485 fAnalysisManager->FillNtupleIColumn(4, 11,
486 fAnalysisManager->FillNtupleIColumn(4, 12,
487 fAnalysisManager->FillNtupleIColumn(4, 13,
488 fAnalysisManager->FillNtupleIColumn(4, 14,
489 fAnalysisManager->FillNtupleIColumn(4, 15,
490 fAnalysisManager->FillNtupleIColumn(4, 16,
491 fAnalysisManager->FillNtupleIColumn(4, 17,
492 fAnalysisManager->FillNtupleIColumn(4, 18,
493 fAnalysisManager->FillNtupleIColumn(4, 19,
494 fAnalysisManager->FillNtupleDColumn(4, 20,
495 fAnalysisManager->FillNtupleIColumn(4, 21,
496 fAnalysisManager->FillNtupleIColumn(4, 22,
497 fAnalysisManager->FillNtupleIColumn(4, 23,
498 fAnalysisManager->FillNtupleDColumn(4, 24,
499 // fAnalysisManager->FillNtupleIColumn(4,
500 fAnalysisManager->FillNtupleSColumn(4, 25,
501 fAnalysisManager->AddNtupleRow(4);
502 delete classif;
503 //
504 // TODO: Use (*it)->GetPlacementIdx() to w
505 // is continuously joined to a past event
506 // If yes, save the fragment length.
507 //
508 if (fSaveStrands) {
509 damageRecord->PrintRecord(fStrandDirecto
510 + damageRecord
511 }
512 if (fAnalysisManager->GetVerboseLevel() >=
513 }
514
515 fAnalysisManager->FillNtupleSColumn(
516 2, 0, event->GetPrimaryVertex()->GetPrimar
517 fAnalysisManager->FillNtupleDColumn(2, 1,
518 event->G
519 fAnalysisManager->FillNtupleIColumn(2, 2, br
520 fAnalysisManager->FillNtupleIColumn(2, 3, br
521 fAnalysisManager->FillNtupleIColumn(2, 4, br
522 fAnalysisManager->FillNtupleIColumn(2, 5, br
523 fAnalysisManager->FillNtupleIColumn(2, 6, br
524 fAnalysisManager->FillNtupleIColumn(2, 7, br
525 fAnalysisManager->FillNtupleIColumn(2, 8, br
526 fAnalysisManager->AddNtupleRow(2);
527
528 fAnalysisManager->FillNtupleSColumn(
529 3, 0, event->GetPrimaryVertex()->GetPrimar
530 fAnalysisManager->FillNtupleDColumn(3, 1,
531 event->G
532 fAnalysisManager->FillNtupleIColumn(3, 2, so
533 fAnalysisManager->FillNtupleIColumn(3, 3, so
534 fAnalysisManager->FillNtupleIColumn(3, 4, so
535 fAnalysisManager->FillNtupleIColumn(3, 5, so
536 fAnalysisManager->FillNtupleIColumn(3, 6, so
537 fAnalysisManager->FillNtupleIColumn(3, 7, so
538 fAnalysisManager->FillNtupleIColumn(3, 8, so
539 fAnalysisManager->FillNtupleIColumn(3, 9, so
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........oo
553
554 void AnalysisManager::ProcessChromosomeHitMap(
555 {
556 for (auto chromosome : chromosomes) {
557 fAnalysisManager->FillNtupleSColumn(1, 0,
558 fAnalysisManager->FillNtupleDColumn(1, 1,
559 fAnalysisManager->FillNtupleDColumn(1, 2,
560 fAnalysisManager->AddNtupleRow(1);
561 }
562 }
563
564 //....oooOO0OOooo........oooOO0OOooo........oo
565
566 void AnalysisManager::ProcessPrimary(const G4T
567 const G4d
568 {
569 const G4Event* event = G4EventManager::GetEv
570 fAnalysisManager->FillNtupleSColumn(
571 0, 0, event->GetPrimaryVertex()->GetPrimar
572 fAnalysisManager->FillNtupleDColumn(0, 1,
573 event->G
574 fAnalysisManager->FillNtupleDColumn(0, 2, ev
575 fAnalysisManager->FillNtupleDColumn(0, 3, ev
576 fAnalysisManager->FillNtupleDColumn(0, 4, ev
577
578 G4ThreeVector mom = event->GetPrimaryVertex(
579 fAnalysisManager->FillNtupleDColumn(0, 5, mo
580 fAnalysisManager->FillNtupleDColumn(0, 6, mo
581 fAnalysisManager->FillNtupleDColumn(0, 7, mo
582 fAnalysisManager->FillNtupleDColumn(0, 8, pr
583 fAnalysisManager->FillNtupleDColumn(0, 9, pr
584 fAnalysisManager->FillNtupleDColumn(0, 10, p
585 fAnalysisManager->FillNtupleDColumn(0, 11, t
586 fAnalysisManager->FillNtupleDColumn(0, 12, t
587 fAnalysisManager->AddNtupleRow(0);
588 }
589
590 //....oooOO0OOooo........oooOO0OOooo........oo
591
592 void AnalysisManager::ProcessCellEdep(const G4
593 {
594 fAnalysisManager->FillH1(fAnalysisManager->G
595 }
596
597 //////////////////////////////////////////////
598 ////// Binary Tree Methods
599 //////////////////////////////////////////////
600
601 BinaryTree::BinaryTree()
602 {
603 fRoot = nullptr;
604 }
605
606 //....oooOO0OOooo........oooOO0OOooo........oo
607
608 BinaryTree::~BinaryTree()
609 {
610 Destroy_tree();
611 }
612
613 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
634
635 DNAHit* BinaryTree::Search(int64_t index) //
636 {
637 return Search_(index, fRoot);
638 }
639
640 //....oooOO0OOooo........oooOO0OOooo........oo
641
642 void BinaryTree::Destroy_tree()
643 {
644 Destroy_tree_(fRoot);
645 }
646
647 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
660
661 void BinaryTree::Insert_(DNAHit* hit, Node* no
662 {
663 if (hit->GetBasePairIdx() > 9223372036854775
664 G4cout << " SEE AnalysisManager !!!" << G4
665 G4ExceptionDescription exceptionDescriptio
666 exceptionDescription << "!!!!!!!!!!!!!!!!!
667 G4cout << "!!! aborting ... " << G4endl;
668 G4cout << "This pair ID is so big " << hit
669 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl
670 G4Exception(
671 "BinaryTree"
672 "::insert()",
673 "insert000", FatalException, exceptionDe
674 }
675 int64_t key = hit->GetBasePairIdx(); // dou
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........oo
711
712 DNAHit* BinaryTree::Search_(int64_t key, Node*
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........oo
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........oo
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........oo
755
756 DNAHit* BinaryTree::Next(const DNAHit* hit) co
757 {
758 if (hit->GetBasePairIdx() > 9223372036854775
759 G4cout << " SEE AnalysisManager !!!" << G4
760 G4ExceptionDescription exceptionDescriptio
761 exceptionDescription << "!!!!!!!!!!!!!!!!!
762 G4cout << "!!! aborting ... " << G4endl;
763 G4cout << "This pair ID is so big " << hit
764 G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl
765 G4Exception(
766 "BinaryTree"
767 "::insert()",
768 "insert000", FatalException, exceptionDe
769 }
770 int64_t key = hit->GetBasePairIdx(); // dou
771 if (fRoot == nullptr) {
772 return nullptr;
773 }
774 else {
775 return Next_(key, fRoot);
776 }
777 }
778
779 //....oooOO0OOooo........oooOO0OOooo........oo
780
781 // Looks more complicated than it is
782 // The interface to these functions with integ
783 // objects makes this look more complicated th
784 DNAHit* BinaryTree::Next_(int64_t key, Node* n
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
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........oo
816
817 //////////////////////////////////////////////
818 ////////// Damage Record
819 //////////////////////////////////////////////
820
821 const char* DamageRecord::fDirectDamageChar =
822 const char* DamageRecord::fIndirectDamageChar
823 const char* DamageRecord::fHitNoDamageChar = "
824 const char* DamageRecord::fNotHitChar = "-";
825 const char* DamageRecord::fBothDamageChar = "X
826
827 DamageRecord::DamageRecord(const G4String& nam
828 G4int chain_idx)
829 : fName(name), fStartIndex(startIndex), fSta
830 {}
831
832 //....oooOO0OOooo........oooOO0OOooo........oo
833
834 DamageRecord::~DamageRecord()
835 {
836 for (auto& fDamageRecord : fDamageRecords) {
837 delete fDamageRecord;
838 }
839 }
840
841 //....oooOO0OOooo........oooOO0OOooo........oo
842
843 void DamageRecord::PrintRecord(const G4String&
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->fStrand1
851 fDamageRecord->fStrand1
852 bp1 << GetChar(fDamageRecord->fbp1DirectDm
853 fDamageRecord->fBp1Energy);
854 bp2 << GetChar(fDamageRecord->fbp2DirectDm
855 fDamageRecord->fBp2Energy);
856 strand2 << GetChar(fDamageRecord->fStrand2
857 fDamageRecord->fStrand2
858 }
859
860 // print to with no filename
861 if (filename.empty()) {
862 DamageClassification* classification = thi
863 G4cout << "Sequences starts at base pair "
864 << ", Chain " << fChainIdx << ". To
865 << ". Classification: " << classifi
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 |
875 if (fs.fail()) {
876 G4cout << "Could not open filestream" << G
877 }
878 DamageClassification* classification = this-
879 fs << "Sequences starts at base pair " << fS
880 << ", Chain " << fChainIdx << ". Total le
881 << ". Classification: " << classification
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
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........oo
911
912 DamageClassification* DamageRecord::GetClassif
913 {
914 /* Explanation of class because we do some s
915 track breaks along the strand.
916
917 Things to note straight off
918 1) This routine makes no promise to count
919 once they cease being relevant to the br
920 a DSB has occurred, the SSB count become
921 2) The previous ten base pairs are kept i
922 Each bit in the number represents one ba
923 is set if the strand was broken.
924 3) The countX variables are used to track
925 preceeding 10 base pairs on each strand.
926 DSB+ events
927 4) The lambda function count_bytes is an
928 of set bytes in constant time for 32-bit
929 us the count of breaks in the previous t
930
931 Classification of complexity:
932 DSB++ 2DSBs in fragment
933 DSB+ 1 DSB and then an SSB within 10 ba
934 DSB Each strand hit once within <=10 bas
935 2SSB Each strand hit once with >10 bas
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 by
944 // https://blogs.msdn.microsoft.com/jeuge/20
945 // Note: integers beginning with 0 are enter
946 auto count_bytes = [](uint32_t u) {
947 uint32_t uCount = u - ((u >> 1) & 03333333
948 return ((uCount + (uCount >> 3)) & 0307070
949 };
950
951 auto classification = new DamageClassificati
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
974 uint32_t lastTenTracked2 = 0; // been count
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
984 for (G4int ii = 0; ii != (G4int)fDamageRecor
985 lastTenDirectStrand1 = (lastTenDirectStran
986 lastTenDirectStrand2 = (lastTenDirectStran
987 lastTenIndirectStrand1 = (lastTenIndirectS
988 lastTenIndirectStrand2 = (lastTenIndirectS
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
996 lastTenStrand1 = lastTenStrand1 & truncato
997 lastTenStrand2 = lastTenStrand2 & truncato
998 lastTenTracked1 = lastTenTracked1 & trunca
999 lastTenTracked2 = lastTenTracked2 & trunca
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)->fStrand1Direct
1006 strand1Direct = true;
1007 classification->fDirectBreaks++;
1008 }
1009 if (fDamageRecords.at(ii)->fStrand1Indire
1010 strand1Indirect = true;
1011 classification->fIndirectBreaks++;
1012 }
1013 if (fDamageRecords.at(ii)->fStrand2Direct
1014 strand2Direct = true;
1015 classification->fDirectBreaks++;
1016 }
1017 if (fDamageRecords.at(ii)->fStrand2Indire
1018 strand2Indirect = true;
1019 classification->fIndirectBreaks++;
1020 }
1021 if (fDamageRecords.at(ii)->fbp1DirectDmg)
1022 // classification->fDirectBreaks++; //?
1023 }
1024 if (fDamageRecords.at(ii)->fBp1IndirectDm
1025 if (fDamageRecords.at(ii)->fbp1InducedB
1026 // Counts as a hit only if it breaks
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)->fBp2IndirectDm
1036 if (fDamageRecords.at(ii)->fbp2InducedB
1037 // Counts as a hit only if it breaks
1038 classification->fIndirectBreaks++;
1039 strand2Indirect = true;
1040 classification->fInducedBreaks++;
1041 }
1042 }
1043
1044 G4bool s1IndirectBreak =
1045 fDamageRecords.at(ii)->fStrand1Indirect
1046 G4bool s2IndirectBreak =
1047 fDamageRecords.at(ii)->fStrand2Indirect
1048
1049 // strand 1 hit
1050 G4bool strand1hit = fDamageRecords.at(ii)
1051 G4bool strand2hit = fDamageRecords.at(ii)
1052 G4bool bp1hit = fDamageRecords.at(ii)->fb
1053 G4bool bp2hit = fDamageRecords.at(ii)->fb
1054 strandDamage += ((G4int)strand1hit + (G4i
1055 baseDamage += ((G4int)bp1hit + (G4int)bp2
1056
1057 // Update the damage type counters
1058 // strand 1
1059 if (s1IndirectBreak) {
1060 lastTenIndirectStrand1++;
1061 }
1062 if (fDamageRecords.at(ii)->fStrand1Direct
1063 lastTenDirectStrand1++;
1064 }
1065 // strand 2
1066 if (s2IndirectBreak) {
1067 lastTenIndirectStrand2++;
1068 }
1069 if (fDamageRecords.at(ii)->fStrand2Direct
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
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
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
1118 // DSB
1119 if (lastTenDirectStrand1 && lastTenDire
1120 nDSBd = true;
1121 if (lastTenIndirectStrand1 || lastTen
1122 }
1123 if (lastTenIndirectStrand1 && lastTenIn
1124 nDSBi = true;
1125 if (lastTenDirectStrand1 || lastTenDi
1126 if (lastTenDirectStrand1 && lastTenDi
1127 }
1128 if ((lastTenDirectStrand1 && lastTenInd
1129 || (lastTenIndirectStrand1 && lastT
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 || strand2Di
1163 G4bool indirect = (strand1Indirect || stran
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 classificati
1202 errmsg << "I think this is " << complex
1203 << "???" << G4endl;
1204 G4Exception("DamageRecord::GetClassific
1205 source = undefined;
1206 }
1207 }
1208
1209 if (((complexity == NoneComplexity) && (sou
1210 || ((complexity != NoneComplexity) && (
1211 {
1212 G4ExceptionDescription errmsg;
1213 errmsg << "You can't have a complexity wi
1214 errmsg << "I think this is " << complexit
1215 this->PrintRecord("", dsbDistance);
1216 G4Exception("DamageRecord::GetClassificat
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........o
1228
1229 // add a contrived test record to the damage
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 da
1236 // param == 1 direct hit, no damage; param >
1237 void DamageRecord::AddTestDamage(G4int s1, G4
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
1281 }
1282
1283 // Unit test for classification
1284 void AnalysisManager::TestClassification()
1285 {
1286 G4cout << G4endl;
1287 G4cout << "================================
1288 G4cout << "Classification Unit Test Begins"
1289
1290 // None
1291 auto* theRecord = new DamageRecord("Test",
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 = theR
1300 G4cout << "Test None" << G4endl;
1301 theRecord->PrintRecord("", 10);
1302 G4cout << "Classification: " << classificat
1303 << G4endl;
1304 assert(classification->fComplexity == NoneC
1305 assert(classification->fSource == sourceEnu
1306 delete classification;
1307 delete theRecord;
1308 G4cout << "--------------------------------
1309
1310 // None, base damage
1311 theRecord = new DamageRecord("Test", 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->GetClassificati
1320 G4cout << "Test None" << G4endl;
1321 theRecord->PrintRecord("", 10);
1322 G4cout << "Classification: " << classificat
1323 << G4endl;
1324 G4cout << "indirect Breaks = " << classific
1325 G4cout << "base Damage = " << classificatio
1326 assert(classification->fComplexity == NoneC
1327 assert(classification->fSource == sourceEnu
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 << "--------------------------------
1335
1336 // SSB
1337 theRecord = new DamageRecord("Test", 0, 0,
1338 theRecord->AddEmptyBPDamage(10);
1339 theRecord->AddTestDamage(2, -1, -1, -1);
1340 theRecord->AddEmptyBPDamage(10);
1341 classification = theRecord->GetClassificati
1342 G4cout << "Test SSBd" << G4endl;
1343 theRecord->PrintRecord("", 10);
1344 G4cout << "Classification: " << classificat
1345 << G4endl;
1346 assert(classification->fComplexity == SSB);
1347 assert(classification->fSource == SSBd);
1348 delete classification;
1349 delete theRecord;
1350 G4cout << "--------------------------------
1351
1352 // SSBi
1353 theRecord = new DamageRecord("Test", 0, 0,
1354 theRecord->AddEmptyBPDamage(10);
1355 theRecord->AddTestDamage(0, 0, -1, -1);
1356 theRecord->AddEmptyBPDamage(10);
1357 classification = theRecord->GetClassificati
1358 G4cout << "Test SSBi" << G4endl;
1359 theRecord->PrintRecord("", 10);
1360 G4cout << "Classification: " << classificat
1361 << G4endl;
1362 assert(classification->fComplexity == SSB);
1363 assert(classification->fSource == SSBi);
1364 delete classification;
1365 delete theRecord;
1366 G4cout << "--------------------------------
1367
1368 // SSBm / SSB+
1369 theRecord = new DamageRecord("Test", 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->GetClassificati
1376 G4cout << "Test SSBm / SSB+" << G4endl;
1377 theRecord->PrintRecord("", 10);
1378 G4cout << "Classification: " << classificat
1379 << G4endl;
1380 assert(classification->fComplexity == SSBpl
1381 assert(classification->fSource == SSBm);
1382 delete classification;
1383 delete theRecord;
1384 G4cout << "--------------------------------
1385
1386 // SSBm / 2SSB
1387 theRecord = new DamageRecord("Test", 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->GetClassificati
1394 G4cout << "Test SSBd / 2SSB" << G4endl;
1395 theRecord->PrintRecord("", 10);
1396 G4cout << "Classification: " << classificat
1397 << G4endl;
1398 assert(classification->fComplexity == twoSS
1399 assert(classification->fSource == SSBm);
1400 delete classification;
1401 delete theRecord;
1402 G4cout << "--------------------------------
1403
1404 // SSBi / 2SSB
1405 theRecord = new DamageRecord("Test", 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->GetClassificati
1414 G4cout << "Test SSBi / 2SSB" << G4endl;
1415 theRecord->PrintRecord("", 10);
1416 G4cout << "Classification: " << classificat
1417 << G4endl;
1418 assert(classification->fComplexity == twoSS
1419 assert(classification->fSource == SSBi);
1420 delete classification;
1421 delete theRecord;
1422 G4cout << "--------------------------------
1423
1424 // SSBd / SSB+
1425 theRecord = new DamageRecord("Test", 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->GetClassificati
1432 G4cout << "Test SSBd / SSB+" << G4endl;
1433 theRecord->PrintRecord("", 10);
1434 G4cout << "Classification: " << classificat
1435 << G4endl;
1436 assert(classification->fComplexity == SSBpl
1437 assert(classification->fSource == SSBd);
1438 delete classification;
1439 delete theRecord;
1440 G4cout << "--------------------------------
1441
1442 // SSBm / SSB2
1443 theRecord = new DamageRecord("Test", 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->GetClassificati
1450 G4cout << "Test SSBm / 2SSB" << G4endl;
1451 theRecord->PrintRecord("", 10);
1452 G4cout << "Classification: " << classificat
1453 << G4endl;
1454 assert(classification->fComplexity == twoSS
1455 assert(classification->fSource == SSBm);
1456 delete classification;
1457 delete theRecord;
1458 G4cout << "--------------------------------
1459
1460 // SSBd / SSB2
1461 theRecord = new DamageRecord("Test", 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->GetClassificati
1468 G4cout << "Test SSBd / 2SSB" << G4endl;
1469 theRecord->PrintRecord("", 10);
1470 G4cout << "Classification: " << classificat
1471 << G4endl;
1472 assert(classification->fComplexity == twoSS
1473 assert(classification->fSource == SSBd);
1474 delete classification;
1475 delete theRecord;
1476 G4cout << "--------------------------------
1477
1478 // SSBd / SSB2
1479 theRecord = new DamageRecord("Test", 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->GetClassificati
1486 G4cout << "Test SSBd / 2SSB" << G4endl;
1487 theRecord->PrintRecord("", 10);
1488 G4cout << "Classification: " << classificat
1489 << G4endl;
1490 assert(classification->fComplexity == twoSS
1491 assert(classification->fSource == SSBd);
1492 delete classification;
1493 delete theRecord;
1494 G4cout << "--------------------------------
1495
1496 // DSBd / DSB
1497 theRecord = new DamageRecord("Test", 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->GetClassificati
1504 G4cout << "Test DSBd/DSB" << G4endl;
1505 theRecord->PrintRecord("", 10);
1506 G4cout << "Classification: " << classificat
1507 << G4endl;
1508 assert(classification->fComplexity == DSB);
1509 assert(classification->fSource == DSBd);
1510 delete classification;
1511 delete theRecord;
1512 G4cout << "--------------------------------
1513
1514 G4cout << G4endl;
1515
1516 // Induced Breaks DSB
1517 theRecord = new DamageRecord("Test", 0, 0,
1518 theRecord->AddEmptyBPDamage(10);
1519
1520 auto* bp = new BasePairDamageRecord();
1521 bp->fBp1IndirectDmg = true;
1522 bp->fbp1InducedBreak = true;
1523 theRecord->AddBasePairDamage(bp, G4ThreeVec
1524
1525 bp = new BasePairDamageRecord();
1526 bp->fBp2IndirectDmg = true;
1527 bp->fbp2InducedBreak = true;
1528 theRecord->AddBasePairDamage(bp, G4ThreeVec
1529
1530 theRecord->AddEmptyBPDamage(10);
1531 classification = theRecord->GetClassificati
1532 G4cout << "Test DSBi/DSB induced" << G4endl
1533 theRecord->PrintRecord("", 10);
1534 G4cout << "Classification: " << classificat
1535 << G4endl;
1536 assert(classification->fComplexity == DSB);
1537 assert(classification->fSource == DSBi);
1538 delete classification;
1539 delete theRecord;
1540 G4cout << "--------------------------------
1541
1542 // Induced Breaks SSB
1543 theRecord = new DamageRecord("Test", 0, 0,
1544 theRecord->AddEmptyBPDamage(10);
1545
1546 bp = new BasePairDamageRecord();
1547 bp->fBp2IndirectDmg = true;
1548 bp->fbp2InducedBreak = true;
1549 theRecord->AddBasePairDamage(bp, G4ThreeVec
1550
1551 theRecord->AddEmptyBPDamage(10);
1552 classification = theRecord->GetClassificati
1553 G4cout << "Test SSBi/SSB induced" << G4endl
1554 theRecord->PrintRecord("", 10);
1555 G4cout << "Classification: " << classificat
1556 << G4endl;
1557 assert(classification->fComplexity == SSB);
1558 assert(classification->fSource == SSBi);
1559 delete classification;
1560 delete theRecord;
1561 G4cout << "--------------------------------
1562
1563 // DSBi / DSB
1564 theRecord = new DamageRecord("Test", 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->GetClassificati
1571 G4cout << "Test DSBi/DSB" << G4endl;
1572 theRecord->PrintRecord("", 10);
1573 G4cout << "Classification: " << classificat
1574 << G4endl;
1575 assert(classification->fComplexity == DSB);
1576 assert(classification->fSource == DSBi);
1577 delete classification;
1578 delete theRecord;
1579 G4cout << "--------------------------------
1580
1581 // DSBd / DSB
1582 theRecord = new DamageRecord("Test", 0, 0,
1583 theRecord->AddEmptyBPDamage(10);
1584 theRecord->AddTestDamage(2, -1, -1, 2);
1585 theRecord->AddEmptyBPDamage(10);
1586 classification = theRecord->GetClassificati
1587 G4cout << "Test DSBd/DSB" << G4endl;
1588 theRecord->PrintRecord("", 10);
1589 G4cout << "Classification: " << classificat
1590 << G4endl;
1591 assert(classification->fComplexity == DSB);
1592 assert(classification->fSource == DSBd);
1593 delete classification;
1594 delete theRecord;
1595 G4cout << "--------------------------------
1596
1597 // DSBi/ DSB
1598 theRecord = new DamageRecord("Test", 0, 0,
1599 theRecord->AddEmptyBPDamage(10);
1600 theRecord->AddTestDamage(0, -1, -1, 0);
1601 theRecord->AddEmptyBPDamage(10);
1602 classification = theRecord->GetClassificati
1603 G4cout << "Test DSBd/DSB" << G4endl;
1604 theRecord->PrintRecord("", 10);
1605 G4cout << "Classification: " << classificat
1606 << G4endl;
1607 assert(classification->fComplexity == DSB);
1608 assert(classification->fSource == DSBi);
1609 delete classification;
1610 delete theRecord;
1611 G4cout << "--------------------------------
1612
1613 // DSBh / DSB
1614 theRecord = new DamageRecord("Test", 0, 0,
1615 theRecord->AddEmptyBPDamage(10);
1616 theRecord->AddTestDamage(0, -1, -1, 2);
1617 theRecord->AddEmptyBPDamage(10);
1618 classification = theRecord->GetClassificati
1619 G4cout << "Test DSBd/DSB" << G4endl;
1620 theRecord->PrintRecord("", 10);
1621 G4cout << "Classification: " << classificat
1622 << G4endl;
1623 assert(classification->fComplexity == DSB);
1624 assert(classification->fSource == DSBh);
1625 delete classification;
1626 delete theRecord;
1627 G4cout << "--------------------------------
1628
1629 // DSBh / DSB+
1630 theRecord = new DamageRecord("Test", 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->GetClassificati
1636 G4cout << "Test DSBd/DSB+" << G4endl;
1637 theRecord->PrintRecord("", 10);
1638 G4cout << "Classification: " << classificat
1639 << G4endl;
1640 assert(classification->fComplexity == DSBpl
1641 assert(classification->fSource == DSBh);
1642 delete classification;
1643 delete theRecord;
1644 G4cout << "--------------------------------
1645
1646 // DSBm / DSB+
1647 theRecord = new DamageRecord("Test", 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->GetClassificati
1653 G4cout << "Test DSBd/DSB+" << G4endl;
1654 theRecord->PrintRecord("", 10);
1655 G4cout << "Classification: " << classificat
1656 << G4endl;
1657 assert(classification->fComplexity == DSBpl
1658 assert(classification->fSource == DSBm);
1659 delete classification;
1660 delete theRecord;
1661 G4cout << "--------------------------------
1662
1663 // DSBd / DSB+
1664 theRecord = new DamageRecord("Test", 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->GetClassificati
1672 G4cout << "Test DSBd/DSB+" << G4endl;
1673 theRecord->PrintRecord("", 10);
1674 G4cout << "Classification: " << classificat
1675 << G4endl;
1676 assert(classification->fComplexity == DSBpl
1677 assert(classification->fSource == DSBd);
1678 delete classification;
1679 delete theRecord;
1680 G4cout << "--------------------------------
1681
1682 // DSBd / DSB+
1683 theRecord = new DamageRecord("Test", 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->GetClassificati
1693 G4cout << "Test DSBd/DSB+" << G4endl;
1694 theRecord->PrintRecord("", 10);
1695 G4cout << "Classification: " << classificat
1696 << G4endl;
1697 assert(classification->fComplexity == DSBpl
1698 assert(classification->fSource == DSBd);
1699 delete classification;
1700 delete theRecord;
1701 G4cout << "--------------------------------
1702
1703 // DSBi / DSB+
1704 theRecord = new DamageRecord("Test", 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->GetClassificati
1714 G4cout << "Test DSBd/DSB+" << G4endl;
1715 theRecord->PrintRecord("", 10);
1716 G4cout << "Classification: " << classificat
1717 << G4endl;
1718 assert(classification->fComplexity == DSBpl
1719 assert(classification->fSource == DSBi);
1720 delete classification;
1721 delete theRecord;
1722 G4cout << "--------------------------------
1723
1724 // DSBh / DSB++
1725 theRecord = new DamageRecord("Test", 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->GetClassificati
1735 G4cout << "Test DSBd/DSB++" << G4endl;
1736 theRecord->PrintRecord("", 10);
1737 G4cout << "Classification: " << classificat
1738 << G4endl;
1739 assert(classification->fComplexity == DSBpl
1740 assert(classification->fSource == DSBh);
1741 delete classification;
1742 delete theRecord;
1743 G4cout << "--------------------------------
1744
1745 // DSBd / DSB
1746 theRecord = new DamageRecord("Test", 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->GetClassificati
1755 G4cout << "Test DSBd/DSB" << G4endl;
1756 theRecord->PrintRecord("", 10);
1757 G4cout << "Classification: " << classificat
1758 << G4endl;
1759 assert(classification->fComplexity == DSB);
1760 assert(classification->fSource == DSBd);
1761 delete classification;
1762 delete theRecord;
1763 G4cout << "--------------------------------
1764
1765 // DSBd / DSB++
1766 theRecord = new DamageRecord("Test", 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->GetClassificati
1776 G4cout << "Test DSBd/DSB++" << G4endl;
1777 theRecord->PrintRecord("", 10);
1778 G4cout << "Classification: " << classificat
1779 << G4endl;
1780 assert(classification->fComplexity == DSBpl
1781 assert(classification->fSource == DSBd);
1782 delete classification;
1783 delete theRecord;
1784 G4cout << "--------------------------------
1785
1786 // DSBh / DSB++
1787 theRecord = new DamageRecord("Test", 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->GetClassificati
1797 G4cout << "Test DSBh/DSB++" << G4endl;
1798 theRecord->PrintRecord("", 10);
1799 G4cout << "Classification: " << classificat
1800 << G4endl;
1801 assert(classification->fComplexity == DSBpl
1802 assert(classification->fSource == DSBh);
1803 delete classification;
1804 delete theRecord;
1805 G4cout << "--------------------------------
1806
1807 // DSBm / DSB++
1808 theRecord = new DamageRecord("Test", 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->GetClassificati
1818 G4cout << "Test DSBm/DSB++" << G4endl;
1819 theRecord->PrintRecord("", 10);
1820 G4cout << "Classification: " << classificat
1821 << G4endl;
1822 assert(classification->fComplexity == DSBpl
1823 assert(classification->fSource == DSBm);
1824 delete classification;
1825 delete theRecord;
1826 G4cout << "--------------------------------
1827
1828 // DSBi / DSB+
1829 theRecord = new DamageRecord("Test", 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->GetClassificati
1839 G4cout << "Test DSBi/DSB++" << G4endl;
1840 theRecord->PrintRecord("", 10);
1841 G4cout << "Classification: " << classificat
1842 << G4endl;
1843 assert(classification->fComplexity == DSBpl
1844 assert(classification->fSource == DSBi);
1845 delete classification;
1846 delete theRecord;
1847 G4cout << "--------------------------------
1848
1849 // DSBi / DSB++
1850 theRecord = new DamageRecord("Test", 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->GetClassificati
1860 G4cout << "Test DSBi/DSB++" << G4endl;
1861 theRecord->PrintRecord("", 10);
1862 G4cout << "Classification: " << classificat
1863 << G4endl;
1864 assert(classification->fComplexity == DSBpl
1865 assert(classification->fSource == DSBi);
1866 delete classification;
1867 delete theRecord;
1868 G4cout << "--------------------------------
1869
1870 // DSBd / DSB++
1871 theRecord = new DamageRecord("Test", 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->GetClassificati
1882 G4cout << "Test DSBd/DSB++" << G4endl;
1883 theRecord->PrintRecord("", 10);
1884 G4cout << "Classification: " << classificat
1885 << G4endl;
1886 assert(classification->fComplexity == DSBpl
1887 assert(classification->fSource == DSBd);
1888 delete classification;
1889 delete theRecord;
1890 G4cout << "--------------------------------
1891
1892 // DSBi / DSB++
1893 theRecord = new DamageRecord("Test", 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->GetClassificati
1906 G4cout << "Test DSBi/DSB++" << G4endl;
1907 theRecord->PrintRecord("", 10);
1908 G4cout << "Classification: " << classificat
1909 << G4endl;
1910 assert(classification->fComplexity == DSBpl
1911 assert(classification->fSource == DSBi);
1912 delete classification;
1913 delete theRecord;
1914 G4cout << "--------------------------------
1915
1916 // DSBh / DSB++
1917 theRecord = new DamageRecord("Test", 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->GetClassificati
1932 G4cout << "Test DSBh/DSB++" << G4endl;
1933 theRecord->PrintRecord("", 10);
1934 G4cout << "Classification: " << classificat
1935 << G4endl;
1936 assert(classification->fComplexity == DSBpl
1937 assert(classification->fSource == DSBh);
1938 delete classification;
1939 delete theRecord;
1940 G4cout << "--------------------------------
1941
1942 // DSBm / DSB++
1943 theRecord = new DamageRecord("Test", 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->GetClassificati
1958 G4cout << "Test DSBm/DSB++" << G4endl;
1959 theRecord->PrintRecord("", 10);
1960 G4cout << "Classification: " << classificat
1961 << G4endl;
1962 assert(classification->fComplexity == DSBpl
1963 assert(classification->fSource == DSBm);
1964 delete classification;
1965 delete theRecord;
1966 G4cout << "--------------------------------
1967
1968 // DSBi / DSB++
1969 theRecord = new DamageRecord("Test", 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->GetClassificati
1980 G4cout << "Test DSBi/DSB+" << G4endl;
1981 theRecord->PrintRecord("", 10);
1982 G4cout << "Classification: " << classificat
1983 << G4endl;
1984 assert(classification->fComplexity == DSBpl
1985 assert(classification->fSource == DSBi);
1986 delete classification;
1987 delete theRecord;
1988 G4cout << "--------------------------------
1989
1990 G4cout << "Classification Unit Test Complet
1991 G4cout << "================================
1992
1993 G4cout << G4endl;
1994 }
1995
1996 //....oooOO0OOooo........oooOO0OOooo........o
1997
1998 void DamageRecord::AddEmptyBPDamage(int64_t i
1999 {
2000 auto basePairNumber = ii;
2001 while (basePairNumber > 0) {
2002 fDamageRecords.push_back(new BasePairDama
2003 basePairNumber--;
2004 }
2005 }
2006
2007 //....oooOO0OOooo........oooOO0OOooo........o
2008
2009 void DamageRecord::AddStrandHit(const G4Molec
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........o
2023
2024 void DamageRecord::AddBaseHit(const G4Molecul
2025 {
2026 // G4cout << "Hit by " << mol->GetParticleN
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........o
2039
2040 G4ThreeVector DamageRecord::GetMeanPosition()
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_
2055 return G4ThreeVector(0., 0., 0.);
2056 }
2057 else {
2058 G4double factor = 1.0 / fPositions.size()
2059 return G4ThreeVector(x * factor, y * fact
2060 }
2061 }
2062
2063 //....oooOO0OOooo........oooOO0OOooo........o
2064
2065 G4double DamageRecord::GetMeanDistance() cons
2066 {
2067 G4double d = 0.;
2068 for (auto it = fPositions.begin(); it != fP
2069 for (auto jt = it + 1; jt != fPositions.e
2070 d += ((*jt) - (*it)).mag();
2071 }
2072 }
2073 G4int number = fPositions.size();
2074 return d / number;
2075 }
2076
2077 //....oooOO0OOooo........oooOO0OOooo........o
2078
2079 G4double DamageRecord::GetEnergy() const
2080 {
2081 G4double en = 0;
2082 for (auto bp : fDamageRecords) {
2083 en = en + bp->fBp1Energy + bp->fBp2Energy
2084 }
2085 return en;
2086 }