Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:21:59 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// file:
0028 /// brief:
0029 
0030 #include "AnalysisManager.hh"
0031 
0032 #include "AnalysisMessenger.hh"
0033 #include "ChromosomeHit.hh"
0034 #include "DNAGeometry.hh"
0035 #include "DNAHit.hh"
0036 #include "DetectorConstruction.hh"
0037 #include "UtilityFunctions.hh"
0038 
0039 #include "G4Event.hh"
0040 #include "G4EventManager.hh"
0041 #include "G4MolecularConfiguration.hh"
0042 #include "G4Run.hh"
0043 #include "G4RunManager.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "Randomize.hh"
0046 
0047 #include <fstream>
0048 #include <utility>
0049 
0050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0051 
0052 AnalysisManager::AnalysisManager()
0053 {
0054   fAnalysisManager = G4AnalysisManager::Instance();
0055   fpAnalysisMessenger = new AnalysisMessenger(this);
0056 }
0057 
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059 
0060 AnalysisManager::~AnalysisManager()
0061 {
0062   delete fpAnalysisMessenger;
0063 }
0064 
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066 
0067 void AnalysisManager::Initialize()
0068 {
0069   const G4Run* run = G4RunManager::GetRunManager()->GetCurrentRun();
0070   if (run->GetRunID() == 0) {
0071     G4cout << "AnalysisManager::Initialize() GetRunID : " << run->GetRunID() << G4endl;
0072     fAnalysisManager->SetDefaultFileType("root");
0073     fAnalysisManager->SetVerboseLevel(1);
0074 
0075     fAnalysisManager->SetNtupleDirectoryName("tuples");
0076     fAnalysisManager->SetHistoDirectoryName("hists");
0077 
0078     fAnalysisManager->CreateH1("ssb_counts", "SSBs", 1000, 0, 1000);
0079     fAnalysisManager->CreateH1("ssb_energies_ev", "SSB Energies", 500, 0, 500);
0080     // fAnalysisManager->CreateH1("fragments", "Fragment sizes", 500, 0,
0081     // 500);//ORG
0082     fAnalysisManager->CreateH1("fragments", "Fragment sizes", 10000, 0,
0083                                10000);  // dousatsu
0084     fAnalysisManager->CreateH3("local_pos", "Strand Interaction Positions", 50, -25, 25, 50, -25,
0085                                25, 50, -25, 25);
0086     fAnalysisManager->CreateH3("global_pos", "Strand Interaction Positions", 50, -100, 100, 50,
0087                                -100, 100, 50, -100, 100);
0088     fAnalysisManager->CreateH1("e_cell_kev", "Energy Deposits in Cell", 2500, 0,
0089                                2500);  // dousatsu
0090 
0091     fAnalysisManager->CreateNtuple("primary_source", "Primary Source");
0092     fAnalysisManager->CreateNtupleSColumn("Primary");
0093     fAnalysisManager->CreateNtupleDColumn("Energy");
0094     fAnalysisManager->CreateNtupleDColumn("PosX_um");
0095     fAnalysisManager->CreateNtupleDColumn("PosY_um");
0096     fAnalysisManager->CreateNtupleDColumn("PosZ_um");
0097     fAnalysisManager->CreateNtupleDColumn("MomX");
0098     fAnalysisManager->CreateNtupleDColumn("MomY");
0099     fAnalysisManager->CreateNtupleDColumn("MomZ");
0100     fAnalysisManager->CreateNtupleDColumn("StopPosX_um");
0101     fAnalysisManager->CreateNtupleDColumn("StopPosY_um");
0102     fAnalysisManager->CreateNtupleDColumn("StopPosZ_um");
0103     fAnalysisManager->CreateNtupleDColumn("TraLen_cell_um");
0104     fAnalysisManager->CreateNtupleDColumn("TraLen_chro_um");
0105     fAnalysisManager->FinishNtuple();
0106 
0107     fAnalysisManager->CreateNtuple("chromosome_hits", "Energy Deposits in Chromosomes");
0108     fAnalysisManager->CreateNtupleSColumn("chromosome");
0109     fAnalysisManager->CreateNtupleDColumn("e_chromosome_kev");
0110     fAnalysisManager->CreateNtupleDColumn("e_dna_kev");
0111     fAnalysisManager->FinishNtuple();
0112 
0113     fAnalysisManager->CreateNtuple("classification", "Break Complexity Frequency");
0114     fAnalysisManager->CreateNtupleSColumn("Primary");
0115     fAnalysisManager->CreateNtupleDColumn("Energy");
0116     fAnalysisManager->CreateNtupleIColumn("None");
0117     fAnalysisManager->CreateNtupleIColumn("SSB");
0118     fAnalysisManager->CreateNtupleIColumn("SSBp");
0119     fAnalysisManager->CreateNtupleIColumn("2SSB");
0120     fAnalysisManager->CreateNtupleIColumn("DSB");
0121     fAnalysisManager->CreateNtupleIColumn("DSBp");
0122     fAnalysisManager->CreateNtupleIColumn("DSBpp");
0123     fAnalysisManager->FinishNtuple();
0124 
0125     fAnalysisManager->CreateNtuple("source", "Break Source Frequency");
0126     fAnalysisManager->CreateNtupleSColumn("Primary");
0127     fAnalysisManager->CreateNtupleDColumn("Energy");
0128     fAnalysisManager->CreateNtupleIColumn("None");
0129     fAnalysisManager->CreateNtupleIColumn("SSBd");
0130     fAnalysisManager->CreateNtupleIColumn("SSBi");
0131     fAnalysisManager->CreateNtupleIColumn("SSBm");
0132     fAnalysisManager->CreateNtupleIColumn("DSBd");
0133     fAnalysisManager->CreateNtupleIColumn("DSBi");
0134     fAnalysisManager->CreateNtupleIColumn("DSBm");
0135     fAnalysisManager->CreateNtupleIColumn("DSBh");
0136     fAnalysisManager->FinishNtuple();
0137 
0138     fAnalysisManager->CreateNtuple("damage", "DNA damage locations");
0139     fAnalysisManager->CreateNtupleIColumn("Event");
0140     fAnalysisManager->CreateNtupleSColumn("Primary");
0141     fAnalysisManager->CreateNtupleDColumn("Energy");
0142     fAnalysisManager->CreateNtupleSColumn("TypeClassification");
0143     fAnalysisManager->CreateNtupleSColumn("SourceClassification");
0144     fAnalysisManager->CreateNtupleDColumn("Position_x_um");
0145     fAnalysisManager->CreateNtupleDColumn("Position_y_um");
0146     fAnalysisManager->CreateNtupleDColumn("Position_z_um");
0147     fAnalysisManager->CreateNtupleDColumn("Size_nm");
0148     fAnalysisManager->CreateNtupleIColumn("FragmentLength");
0149     fAnalysisManager->CreateNtupleIColumn("BaseDamage");
0150     fAnalysisManager->CreateNtupleIColumn("StrandDamage");
0151     fAnalysisManager->CreateNtupleIColumn("DirectBreaks");
0152     fAnalysisManager->CreateNtupleIColumn("IndirectBreaks");
0153     fAnalysisManager->CreateNtupleIColumn("EaqBaseHits");
0154     fAnalysisManager->CreateNtupleIColumn("EaqStrandHits");
0155     fAnalysisManager->CreateNtupleIColumn("OHBaseHits");
0156     fAnalysisManager->CreateNtupleIColumn("OHStrandHits");
0157     fAnalysisManager->CreateNtupleIColumn("HBaseHits");
0158     fAnalysisManager->CreateNtupleIColumn("HStrandHits");
0159     fAnalysisManager->CreateNtupleDColumn("EnergyDeposited_eV");
0160     fAnalysisManager->CreateNtupleIColumn("InducedBreaks");
0161     fAnalysisManager->CreateNtupleIColumn("Chain");
0162     fAnalysisManager->CreateNtupleIColumn("Strand");  // WG
0163     fAnalysisManager->CreateNtupleDColumn("BasePair");  // dousatsu
0164     fAnalysisManager->CreateNtupleSColumn("Name");
0165     fAnalysisManager->FinishNtuple();
0166   }
0167 
0168   if (!fFileName) {
0169     fFileName = "molecular-dna";
0170   }
0171   fAnalysisManager->OpenFile(fFileName);
0172 }
0173 
0174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0175 
0176 void AnalysisManager::Close()
0177 {
0178   fAnalysisManager->Write();
0179   fAnalysisManager->CloseFile();
0180 }
0181 
0182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0183 
0184 void AnalysisManager::ProcessDNAHitsVector(const std::vector<const DNAHit*>& dnaHits)
0185 {
0186   // Check we have DNA Geometry (runs once)
0187   if (fpDNAGeometry == nullptr) {
0188     auto det = dynamic_cast<const DetectorConstruction*>(
0189       G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0190     fpDNAGeometry = det->GetDNAGeometry();
0191   }
0192 
0193   if (dnaHits.empty()) {
0194     return;
0195   }
0196   const G4Event* event = G4EventManager::GetEventManager()->GetConstCurrentEvent();
0197   // SSBs
0198   for (auto dnaHit : dnaHits) {
0199     fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("ssb_energies_ev"),
0200                              dnaHit->GetEnergy() / eV);
0201     if (fChainToSave == -1) {
0202       // save all chains
0203       fAnalysisManager->FillH3(
0204         fAnalysisManager->GetH3Id("local_pos"), dnaHit->GetLocalPosition().getX() / nm,
0205         dnaHit->GetLocalPosition().getY() / nm, dnaHit->GetLocalPosition().getZ() / nm);
0206       fAnalysisManager->FillH3(fAnalysisManager->GetH3Id("global_pos"),
0207                                dnaHit->GetPosition().getX() / nm, dnaHit->GetPosition().getY() / nm,
0208                                dnaHit->GetPosition().getZ() / nm);
0209     }
0210     else if (fChainToSave == dnaHit->GetChainIdx()) {  // what is different ?
0211       // save only specified chain
0212       fAnalysisManager->FillH3(
0213         fAnalysisManager->GetH3Id("local_pos"), dnaHit->GetLocalPosition().getX() / nm,
0214         dnaHit->GetLocalPosition().getY() / nm, dnaHit->GetLocalPosition().getZ() / nm);
0215       fAnalysisManager->FillH3(fAnalysisManager->GetH3Id("global_pos"),
0216                                dnaHit->GetPosition().getX() / nm, dnaHit->GetPosition().getY() / nm,
0217                                dnaHit->GetPosition().getZ() / nm);
0218     }
0219   }
0220   fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("ssb_counts"), dnaHits.size());
0221 
0222   // DSBs
0223   // Sorting
0224   // Make a map relating chromosome/chain to its binary tree
0225   std::map<std::pair<G4String, G4int>, BinaryTree*> treemap;
0226 
0227   // Populate the binary tree
0228   // G4cout << "Building Binary Trees for " << dnaHits.size() << " Nodes" <<
0229   // G4endl;
0230   for (auto it = dnaHits.begin(); it != dnaHits.end(); ++it) {
0231     std::pair<G4String, G4int> key = std::make_pair((*it)->GetChromosome(), (*it)->GetChainIdx());
0232     if (treemap.find(key) == treemap.end()) {
0233       treemap[key] = new BinaryTree();
0234       if (fAnalysisManager->GetVerboseLevel() >= 2) {
0235         G4cout << "Constructing binary tree. Chromosome: " << key.first << " Chain: " << key.second
0236                << G4endl;
0237       }
0238     }
0239     treemap.at(key)->Insert(*it);
0240     if (fAnalysisManager->GetVerboseLevel() >= 2) {
0241       G4cout << "Adding hit to binary tree. Chromosome: " << key.first << " Chain: " << key.second
0242              << G4endl;
0243     }
0244   }
0245 
0246   // Analysis
0247   // Create a vector to stock damage records for the event
0248   // G4cout << "Building Damage Records" << G4endl;
0249   std::vector<DamageRecord*> damageRecords;
0250   for (auto& it : treemap) {
0251     if (fAnalysisManager->GetVerboseLevel() >= 2) {
0252       G4cout << "Analysing hits for chromosome: " << it.first.first << " Chain: " << it.first.second
0253              << G4endl;
0254     }
0255     DNAHit* currentHit = it.second->First();
0256     DNAHit* nextHit = currentHit;
0257 
0258     // Runs while their are still hits
0259     while (nextHit) {
0260       // Create record for this fragment
0261       if (currentHit->GetBasePairIdx() > 9223372036854775807) {
0262         G4cout << " SEE AnalysisManager !!!" << G4endl;
0263         abort();
0264       }  // dousatsu
0265       // if (currentHit->GetBasePairIdx()>2147483647){G4cout<<" SEE
0266       // AnalysisManager !!!"<<G4endl;abort();}//ORG
0267       auto idx = (int64_t)currentHit->GetBasePairIdx();  // dousatsu
0268       int64_t startidx = (idx > 10) ? (idx - 10) : 0;  // dousatsu
0269       // G4int idx = currentHit->GetBasePairIdx();//ORG
0270       // G4int startidx = (idx > 10) ? (idx - 10) : 0;//ORG
0271       G4String name = it.first.first + "_" + std::to_string(it.first.second) + "_"
0272                       + std::to_string(startidx);  // check!!!!
0273       auto* record =
0274         new DamageRecord(name, idx, currentHit->GetPlacementIdx(), currentHit->GetChainIdx());
0275       damageRecords.push_back(record);
0276       BasePairDamageRecord* bp;
0277       int64_t gap = 0;  // dousatsu
0278 
0279       // bookend with 10 bps on either side of no damage
0280       record->AddEmptyBPDamage(idx - startidx);
0281       // continues until fragment ends
0282       while (true) {
0283         bp = new BasePairDamageRecord;
0284         bp->fStrand1Energy = currentHit->GetStrand1Energy();
0285         bp->fStrand2Energy = currentHit->GetStrand2Energy();
0286         bp->fBp1Energy = currentHit->GetBP1Energy();
0287         bp->fBp2Energy = currentHit->GetBP2Energy();
0288         bp->fBp1IndirectEvt = (currentHit->GetBase1Rad() != nullptr);
0289         bp->fBp2IndirectEvt = (currentHit->GetBase2Rad() != nullptr);
0290         bp->fStrand1IndirectEvt = (currentHit->GetStrand1Rad() != nullptr);
0291         bp->fStrand2IndirectEvt = (currentHit->GetStrand2Rad() != nullptr);
0292         bp->fbp1DirectDmg = false;  // No physical damage on BPs
0293         bp->fbp2DirectDmg = false;  // No physical damage on BPs
0294 
0295         bp->fStrand1DirectDmg =
0296           fpDNAGeometry->GetDamageModel()->IsDirectStrandBreak(bp->fStrand1Energy);
0297         bp->fStrand2DirectDmg =
0298           fpDNAGeometry->GetDamageModel()->IsDirectStrandBreak(bp->fStrand2Energy);
0299 
0300         //        if(bp->fStrand1DirectDmg) bp->fStrandIdx = 1;
0301         //        else if(bp->fStrand2DirectDmg) strandID = 2;
0302 
0303         if (bp->fBp1IndirectEvt) {
0304           bp->fBp1IndirectDmg =
0305             fpDNAGeometry->GetDamageModel()->IsIndirectBaseDamage(currentHit->GetBase1Rad());
0306           if (bp->fBp1IndirectDmg) {
0307             bp->fbp1InducedBreak =
0308               fpDNAGeometry->GetDamageModel()->IsInducedStrandBreak(currentHit->GetBase1Rad());
0309           }
0310         }
0311 
0312         if (bp->fBp2IndirectEvt) {
0313           bp->fBp2IndirectDmg =
0314             fpDNAGeometry->GetDamageModel()->IsIndirectBaseDamage(currentHit->GetBase2Rad());
0315           if (bp->fBp2IndirectDmg) {
0316             bp->fbp2InducedBreak =
0317               fpDNAGeometry->GetDamageModel()->IsInducedStrandBreak(currentHit->GetBase2Rad());
0318           }
0319         }
0320 
0321         if (bp->fStrand1IndirectEvt) {
0322           bp->fStrand1IndirectDmg =
0323             fpDNAGeometry->GetDamageModel()->IsIndirectStrandDamage(currentHit->GetStrand1Rad());
0324         }
0325         if (bp->fStrand2IndirectEvt) {
0326           bp->fStrand2IndirectDmg =
0327             fpDNAGeometry->GetDamageModel()->IsIndirectStrandDamage(currentHit->GetStrand2Rad());
0328         }
0329         // Record radical-base/strand damage.
0330         // it is possible for base/strand damage to not correspond to
0331         // direct hits due to the damage induction probabilities.
0332         // This mainly affects strands.
0333         if (bp->fBp1IndirectDmg) {
0334           record->AddBaseHit(currentHit->GetBase1Rad()->GetDefinition());
0335         }
0336         if (bp->fBp2IndirectDmg) {
0337           record->AddBaseHit(currentHit->GetBase2Rad()->GetDefinition());
0338         }
0339         if (bp->fStrand1IndirectDmg) {
0340           record->AddStrandHit(currentHit->GetStrand1Rad()->GetDefinition());
0341           //         strandID = 1;
0342         }
0343         if (bp->fStrand2IndirectDmg) {
0344           record->AddStrandHit(currentHit->GetStrand2Rad()->GetDefinition());
0345           //          strandID = 2;
0346         }
0347 
0348         record->AddBasePairDamage(bp, currentHit->GetPosition());
0349 
0350         nextHit = it.second->Next(currentHit);
0351         if (nextHit == nullptr) {
0352           break;
0353         }
0354 
0355         gap = nextHit->GetBasePairIdx() - currentHit->GetBasePairIdx();
0356         if (fFragmentGap > 0) {
0357           // case 1: continuous strand
0358           if (gap > fFragmentGap) {
0359             currentHit = nextHit;
0360             break;
0361           }
0362           else {
0363             record->AddEmptyBPDamage(gap);
0364           }
0365         }
0366         else if (fFragmentGap == 0) {
0367           // case 2: individual placements, not joined
0368           // ie. plasmids etc. Each placement is a separate strand
0369           if (currentHit->GetPlacementIdx() != nextHit->GetPlacementIdx()) {
0370             if (fpDNAGeometry->GetVerbosity() > 1) {
0371               G4cout << "Analysis passing to a new placement" << G4endl;
0372             }
0373             currentHit = nextHit;
0374             break;
0375           }
0376           else {
0377             record->AddEmptyBPDamage(gap);
0378           }
0379         }
0380         else {
0381           G4Exception("MolecularAnalaysisManager", "ERR_FRAGMENT_VALUE", FatalException,
0382                       "The value set in /analysisDNA/fragmentGap is bad");
0383         }
0384         currentHit = nextHit;
0385       }
0386       record->AddEmptyBPDamage(10);
0387       // end bookend
0388     }
0389   }
0390   // Count the number of breaks and print out the records
0391   // std::map<G4String, G4int> breakmap;
0392   std::map<complexityEnum, G4int> breakmap;
0393 
0394   breakmap[DSBplusplus] = 0;
0395   breakmap[DSBplus] = 0;
0396   breakmap[DSB] = 0;
0397   breakmap[twoSSB] = 0;
0398   breakmap[SSBplus] = 0;
0399   breakmap[SSB] = 0;
0400   breakmap[NoneComplexity] = 0;
0401   std::map<sourceEnum, G4int> sourcemap;
0402   sourcemap[SSBd] = 0;
0403   sourcemap[SSBi] = 0;
0404   sourcemap[SSBm] = 0;
0405   sourcemap[DSBh] = 0;
0406   sourcemap[DSBm] = 0;
0407   sourcemap[DSBd] = 0;
0408   sourcemap[DSBi] = 0;
0409   sourcemap[undefined] = 0;
0410 
0411   for (auto& damageRecord : damageRecords) {
0412     DamageClassification* classif = damageRecord->GetClassification(fDSBDistance);
0413     breakmap[classif->fComplexity]++;
0414     sourcemap[classif->fSource]++;
0415     fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("fragments"), damageRecord->GetSize());
0416 
0417     ///////////dousatsu--------------------------------------------
0418     fAnalysisManager->FillNtupleIColumn(4, 0, event->GetEventID());
0419     fAnalysisManager->FillNtupleSColumn(
0420       4, 1, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
0421     fAnalysisManager->FillNtupleDColumn(
0422       4, 2, event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
0423 
0424     G4String complexityString = "None";
0425     switch (classif->fComplexity) {
0426       case SSB:
0427         complexityString = "SSB";
0428         break;
0429       case SSBplus:
0430         complexityString = "SSB+";
0431         break;
0432       case twoSSB:
0433         complexityString = "2SSB";
0434         break;
0435       case DSB:
0436         complexityString = "DSB";
0437         break;
0438       case DSBplus:
0439         complexityString = "DSB+";
0440         break;
0441       case DSBplusplus:
0442         complexityString = "DSB++";
0443         break;
0444       default:
0445         complexityString = "None";
0446         break;
0447     }
0448     fAnalysisManager->FillNtupleSColumn(4, 3, complexityString);
0449 
0450     G4String sourceString = "undefined";
0451     switch (classif->fSource) {
0452       case SSBd:
0453         sourceString = "SSBd";
0454         break;
0455       case SSBi:
0456         sourceString = "SSBi";
0457         break;
0458       case SSBm:
0459         sourceString = "SSBm";
0460         break;
0461       case DSBh:
0462         sourceString = "DSBh";
0463         break;
0464       case DSBm:
0465         sourceString = "DSBm";
0466         break;
0467       case DSBd:
0468         sourceString = "DSBd";
0469         break;
0470       case DSBi:
0471         sourceString = "DSBi";
0472         break;
0473       default:
0474         sourceString = "undefined";
0475         break;
0476     }
0477 
0478     fAnalysisManager->FillNtupleSColumn(4, 4, sourceString);
0479     fAnalysisManager->FillNtupleDColumn(4, 5, damageRecord->GetMeanPosition().getX() / um);
0480     fAnalysisManager->FillNtupleDColumn(4, 6, damageRecord->GetMeanPosition().getY() / um);
0481     fAnalysisManager->FillNtupleDColumn(4, 7, damageRecord->GetMeanPosition().getZ() / um);
0482     fAnalysisManager->FillNtupleDColumn(4, 8, damageRecord->GetMeanDistance() / nm);
0483     fAnalysisManager->FillNtupleIColumn(4, 9, damageRecord->GetSize());
0484     fAnalysisManager->FillNtupleIColumn(4, 10, classif->fbaseDmg);
0485     fAnalysisManager->FillNtupleIColumn(4, 11, classif->fStrandDmg);
0486     fAnalysisManager->FillNtupleIColumn(4, 12, classif->fDirectBreaks);
0487     fAnalysisManager->FillNtupleIColumn(4, 13, classif->fIndirectBreaks);
0488     fAnalysisManager->FillNtupleIColumn(4, 14, damageRecord->GetEaqBaseHits());
0489     fAnalysisManager->FillNtupleIColumn(4, 15, damageRecord->GetEaqStrandHits());
0490     fAnalysisManager->FillNtupleIColumn(4, 16, damageRecord->GetOHBaseHits());
0491     fAnalysisManager->FillNtupleIColumn(4, 17, damageRecord->GetOHStrandHits());
0492     fAnalysisManager->FillNtupleIColumn(4, 18, damageRecord->GetHBaseHits());
0493     fAnalysisManager->FillNtupleIColumn(4, 19, damageRecord->GetHStrandHits());
0494     fAnalysisManager->FillNtupleDColumn(4, 20, damageRecord->GetEnergy() / eV);
0495     fAnalysisManager->FillNtupleIColumn(4, 21, classif->fInducedBreaks);
0496     fAnalysisManager->FillNtupleIColumn(4, 22, damageRecord->GetChainIdx());
0497     fAnalysisManager->FillNtupleIColumn(4, 23, damageRecord->GetPlacementIdx());
0498     fAnalysisManager->FillNtupleDColumn(4, 24, (G4double)damageRecord->GetStartBPIdx());
0499     // fAnalysisManager->FillNtupleIColumn(4, 23, (*it)->GetStartBPIdx());//ORG
0500     fAnalysisManager->FillNtupleSColumn(4, 25, damageRecord->GetName());
0501     fAnalysisManager->AddNtupleRow(4);
0502     delete classif;
0503     //
0504     // TODO: Use (*it)->GetPlacementIdx() to work out if the damage
0505     // is continuously joined to a past event or not.
0506     // If yes, save the fragment length.
0507     //
0508     if (fSaveStrands) {
0509       damageRecord->PrintRecord(fStrandDirectory + "/" + std::to_string(event->GetEventID()) + "_"
0510                                 + damageRecord->GetName() + ".txt");
0511     }
0512     if (fAnalysisManager->GetVerboseLevel() >= 2) damageRecord->PrintRecord("");
0513   }
0514 
0515   fAnalysisManager->FillNtupleSColumn(
0516     2, 0, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
0517   fAnalysisManager->FillNtupleDColumn(2, 1,
0518                                       event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
0519   fAnalysisManager->FillNtupleIColumn(2, 2, breakmap[NoneComplexity]);
0520   fAnalysisManager->FillNtupleIColumn(2, 3, breakmap[SSB]);
0521   fAnalysisManager->FillNtupleIColumn(2, 4, breakmap[SSBplus]);
0522   fAnalysisManager->FillNtupleIColumn(2, 5, breakmap[twoSSB]);
0523   fAnalysisManager->FillNtupleIColumn(2, 6, breakmap[DSB]);
0524   fAnalysisManager->FillNtupleIColumn(2, 7, breakmap[DSBplus]);
0525   fAnalysisManager->FillNtupleIColumn(2, 8, breakmap[DSBplusplus]);
0526   fAnalysisManager->AddNtupleRow(2);
0527 
0528   fAnalysisManager->FillNtupleSColumn(
0529     3, 0, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
0530   fAnalysisManager->FillNtupleDColumn(3, 1,
0531                                       event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
0532   fAnalysisManager->FillNtupleIColumn(3, 2, sourcemap[undefined]);
0533   fAnalysisManager->FillNtupleIColumn(3, 3, sourcemap[SSBd]);
0534   fAnalysisManager->FillNtupleIColumn(3, 4, sourcemap[SSBi]);
0535   fAnalysisManager->FillNtupleIColumn(3, 5, sourcemap[SSBm]);
0536   fAnalysisManager->FillNtupleIColumn(3, 6, sourcemap[DSBd]);
0537   fAnalysisManager->FillNtupleIColumn(3, 7, sourcemap[DSBi]);
0538   fAnalysisManager->FillNtupleIColumn(3, 8, sourcemap[DSBm]);
0539   fAnalysisManager->FillNtupleIColumn(3, 9, sourcemap[DSBh]);
0540   fAnalysisManager->AddNtupleRow(3);
0541 
0542   for (auto it : damageRecords) {
0543     delete it;
0544   }
0545 
0546   // Cleanup, delete binary trees
0547   for (const auto& it : treemap) {
0548     delete it.second;
0549   }
0550 }
0551 
0552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0553 
0554 void AnalysisManager::ProcessChromosomeHitMap(const std::map<uint32_t, ChromosomeHit*>& chromosomes)
0555 {
0556   for (auto chromosome : chromosomes) {
0557     fAnalysisManager->FillNtupleSColumn(1, 0, chromosome.second->GetName());
0558     fAnalysisManager->FillNtupleDColumn(1, 1, chromosome.second->GetChromosomeEdep() / keV);
0559     fAnalysisManager->FillNtupleDColumn(1, 2, chromosome.second->GetDNAEdep() / keV);
0560     fAnalysisManager->AddNtupleRow(1);
0561   }
0562 }
0563 
0564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0565 
0566 void AnalysisManager::ProcessPrimary(const G4ThreeVector& primstoppos, const G4double& tlcell,
0567                                      const G4double& tlchro)
0568 {
0569   const G4Event* event = G4EventManager::GetEventManager()->GetConstCurrentEvent();
0570   fAnalysisManager->FillNtupleSColumn(
0571     0, 0, event->GetPrimaryVertex()->GetPrimary()->GetParticleDefinition()->GetParticleName());
0572   fAnalysisManager->FillNtupleDColumn(0, 1,
0573                                       event->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy());
0574   fAnalysisManager->FillNtupleDColumn(0, 2, event->GetPrimaryVertex()->GetX0() / um);
0575   fAnalysisManager->FillNtupleDColumn(0, 3, event->GetPrimaryVertex()->GetY0() / um);
0576   fAnalysisManager->FillNtupleDColumn(0, 4, event->GetPrimaryVertex()->GetZ0() / um);
0577 
0578   G4ThreeVector mom = event->GetPrimaryVertex()->GetPrimary()->GetMomentumDirection();
0579   fAnalysisManager->FillNtupleDColumn(0, 5, mom.x() / mom.mag());
0580   fAnalysisManager->FillNtupleDColumn(0, 6, mom.y() / mom.mag());
0581   fAnalysisManager->FillNtupleDColumn(0, 7, mom.z() / mom.mag());
0582   fAnalysisManager->FillNtupleDColumn(0, 8, primstoppos.x() / um);
0583   fAnalysisManager->FillNtupleDColumn(0, 9, primstoppos.y() / um);
0584   fAnalysisManager->FillNtupleDColumn(0, 10, primstoppos.z() / um);
0585   fAnalysisManager->FillNtupleDColumn(0, 11, tlcell / um);
0586   fAnalysisManager->FillNtupleDColumn(0, 12, tlchro / um);
0587   fAnalysisManager->AddNtupleRow(0);
0588 }
0589 
0590 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0591 
0592 void AnalysisManager::ProcessCellEdep(const G4double& edep)
0593 {
0594   fAnalysisManager->FillH1(fAnalysisManager->GetH1Id("e_cell_kev"), edep / keV);
0595 }
0596 
0597 ////////////////////////////////////////////////////////////////////////////////
0598 ////// Binary Tree Methods
0599 ////////////////////////////////////////////////////////////////////////////////
0600 
0601 BinaryTree::BinaryTree()
0602 {
0603   fRoot = nullptr;
0604 }
0605 
0606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0607 
0608 BinaryTree::~BinaryTree()
0609 {
0610   Destroy_tree();
0611 }
0612 
0613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0614 
0615 // Makes own internal copy of hit
0616 void BinaryTree::Insert(const DNAHit* hit)
0617 {
0618   auto newHit = new DNAHit(*hit);
0619   if (fRoot == nullptr) {
0620     Node* node = new Node;
0621     node->fleft = nullptr;
0622     node->fright = nullptr;
0623     node->fparent = nullptr;
0624     node->fdata = newHit;
0625     node->fkey = newHit->GetBasePairIdx();
0626     fRoot = node;
0627   }
0628   else {
0629     Insert_(newHit, fRoot);
0630   }
0631 }
0632 
0633 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0634 
0635 DNAHit* BinaryTree::Search(int64_t index)  // dousatsu
0636 {
0637   return Search_(index, fRoot);
0638 }
0639 
0640 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0641 
0642 void BinaryTree::Destroy_tree()
0643 {
0644   Destroy_tree_(fRoot);
0645 }
0646 
0647 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0648 
0649 void BinaryTree::Destroy_tree_(Node* node)
0650 {
0651   if (node != nullptr) {
0652     Destroy_tree_(node->fleft);
0653     Destroy_tree_(node->fright);
0654     delete node->fdata;
0655     delete node;
0656   }
0657 }
0658 
0659 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0660 
0661 void BinaryTree::Insert_(DNAHit* hit, Node* node)
0662 {
0663   if (hit->GetBasePairIdx() > 9223372036854775807) {
0664     G4cout << " SEE AnalysisManager !!!" << G4endl;
0665     G4ExceptionDescription exceptionDescription;
0666     exceptionDescription << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0667     G4cout << "!!! aborting ... " << G4endl;
0668     G4cout << "This pair ID is so big " << hit->GetBasePairIdx() << G4endl;
0669     G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0670     G4Exception(
0671       "BinaryTree"
0672       "::insert()",
0673       "insert000", FatalException, exceptionDescription);
0674   }
0675   int64_t key = hit->GetBasePairIdx();  // dousatsu
0676   // G4int key = hit->GetBasePairIdx();//ORG
0677   if (key < node->fkey) {
0678     if (node->fleft != nullptr) {
0679       Insert_(hit, node->fleft);
0680     }
0681     else {
0682       Node* daughter = new Node;
0683       daughter->fparent = node;
0684       daughter->fleft = nullptr;
0685       daughter->fright = nullptr;
0686       daughter->fdata = hit;
0687       daughter->fkey = hit->GetBasePairIdx();
0688       node->fleft = daughter;
0689     }
0690   }
0691   else if (key > node->fkey) {
0692     if (node->fright != nullptr) {
0693       Insert_(hit, node->fright);
0694     }
0695     else {
0696       Node* daughter = new Node;
0697       daughter->fparent = node;
0698       daughter->fleft = nullptr;
0699       daughter->fright = nullptr;
0700       daughter->fdata = hit;
0701       daughter->fkey = hit->GetBasePairIdx();
0702       node->fright = daughter;
0703     }
0704   }
0705   else {
0706     node->fdata->AddHit(*hit);
0707   }
0708 }
0709 
0710 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0711 
0712 DNAHit* BinaryTree::Search_(int64_t key, Node* node)  // dousatsu check
0713 {
0714   if (node != nullptr) {
0715     if (key < node->fkey) {
0716       return Search_(key, node->fleft);
0717     }
0718     else if (key > node->fkey) {
0719       return Search_(key, node->fright);
0720     }
0721     else {
0722       return node->fdata;
0723     }
0724   }
0725   else {
0726     return nullptr;
0727   }
0728 }
0729 
0730 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0731 
0732 DNAHit* BinaryTree::First() const
0733 {
0734   if (fRoot == nullptr) {
0735     return nullptr;
0736   }
0737   else {
0738     return First_(fRoot);
0739   }
0740 }
0741 
0742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0743 
0744 DNAHit* BinaryTree::First_(Node* node) const
0745 {
0746   if (node->fleft == nullptr) {
0747     return node->fdata;
0748   }
0749   else {
0750     return First_(node->fleft);
0751   }
0752 }
0753 
0754 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0755 
0756 DNAHit* BinaryTree::Next(const DNAHit* hit) const
0757 {
0758   if (hit->GetBasePairIdx() > 9223372036854775807) {
0759     G4cout << " SEE AnalysisManager !!!" << G4endl;
0760     G4ExceptionDescription exceptionDescription;
0761     exceptionDescription << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0762     G4cout << "!!! aborting ... " << G4endl;
0763     G4cout << "This pair ID is so big " << hit->GetBasePairIdx() << G4endl;
0764     G4cout << "!!!!!!!!!!!!!!!!!!!!" << G4endl;
0765     G4Exception(
0766       "BinaryTree"
0767       "::insert()",
0768       "insert000", FatalException, exceptionDescription);
0769   }
0770   int64_t key = hit->GetBasePairIdx();  // dousatsu
0771   if (fRoot == nullptr) {
0772     return nullptr;
0773   }
0774   else {
0775     return Next_(key, fRoot);
0776   }
0777 }
0778 
0779 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0780 
0781 // Looks more complicated than it is
0782 // The interface to these functions with integer keys rather than node
0783 // objects makes this look more complicated than standard algorithms
0784 DNAHit* BinaryTree::Next_(int64_t key, Node* node) const  // dousatsu
0785 {
0786   if (key < node->fkey) {
0787     if (node->fleft != nullptr) {
0788       return Next_(key, node->fleft);
0789     }
0790     else  // left is NULL
0791     {
0792       return node->fdata;
0793     }
0794   }
0795   else  // (key >= node->key)
0796   {
0797     if (node->fright != nullptr) {
0798       return Next_(key, node->fright);
0799     }
0800     else  // right is NULL; solution is higher in tree
0801     {
0802       while (true) {
0803         node = node->fparent;
0804         if (node == nullptr) {
0805           return nullptr;
0806         }
0807         if (key < node->fkey) {
0808           return node->fdata;
0809         }
0810       }
0811     }
0812   }
0813 }
0814 
0815 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0816 
0817 ////////////////////////////////////////////////////////////////////////////////
0818 ////////// Damage Record
0819 ////////////////////////////////////////////////////////////////////////////////
0820 
0821 const char* DamageRecord::fDirectDamageChar = "D";
0822 const char* DamageRecord::fIndirectDamageChar = "I";
0823 const char* DamageRecord::fHitNoDamageChar = "~";
0824 const char* DamageRecord::fNotHitChar = "-";
0825 const char* DamageRecord::fBothDamageChar = "X";
0826 
0827 DamageRecord::DamageRecord(const G4String& name, int64_t startIndex, G4int place_idx,
0828                            G4int chain_idx)
0829   : fName(name), fStartIndex(startIndex), fStartPlacement(place_idx), fChainIdx(chain_idx)
0830 {}
0831 
0832 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0833 
0834 DamageRecord::~DamageRecord()
0835 {
0836   for (auto& fDamageRecord : fDamageRecords) {
0837     delete fDamageRecord;
0838   }
0839 }
0840 
0841 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0842 
0843 void DamageRecord::PrintRecord(const G4String& filename, const G4double& dsbDistance)
0844 {
0845   std::stringstream strand1;
0846   std::stringstream bp1;
0847   std::stringstream bp2;
0848   std::stringstream strand2;
0849   for (auto& fDamageRecord : fDamageRecords) {
0850     strand1 << GetChar(fDamageRecord->fStrand1DirectDmg, fDamageRecord->fStrand1IndirectDmg,
0851                        fDamageRecord->fStrand1Energy);
0852     bp1 << GetChar(fDamageRecord->fbp1DirectDmg, fDamageRecord->fBp1IndirectDmg,
0853                    fDamageRecord->fBp1Energy);
0854     bp2 << GetChar(fDamageRecord->fbp2DirectDmg, fDamageRecord->fBp2IndirectDmg,
0855                    fDamageRecord->fBp2Energy);
0856     strand2 << GetChar(fDamageRecord->fStrand2DirectDmg, fDamageRecord->fStrand2IndirectDmg,
0857                        fDamageRecord->fStrand2Energy);
0858   }
0859 
0860   // print to with no filename
0861   if (filename.empty()) {
0862     DamageClassification* classification = this->GetClassification(dsbDistance);
0863     G4cout << "Sequences starts at base pair " << fStartIndex << " in placement " << fStartPlacement
0864            << ", Chain " << fChainIdx << ". Total length: " << fDamageRecords.size()
0865            << ". Classification: " << classification->fComplexity << " " << classification->fSource
0866            << G4endl;
0867     G4cout << strand1.str() << G4endl;
0868     G4cout << bp1.str() << G4endl;
0869     G4cout << bp2.str() << G4endl;
0870     G4cout << strand2.str() << G4endl;
0871     return;
0872   }
0873   // To Text output
0874   std::fstream fs(filename, std::fstream::in | std::fstream::out | std::fstream::trunc);
0875   if (fs.fail()) {
0876     G4cout << "Could not open filestream" << G4endl;
0877   }
0878   DamageClassification* classification = this->GetClassification(dsbDistance);
0879   fs << "Sequences starts at base pair " << fStartIndex << " in placement " << fStartPlacement
0880      << ", Chain " << fChainIdx << ". Total length: " << fDamageRecords.size()
0881      << ". Classification: " << classification->fComplexity << " " << classification->fSource
0882      << std::endl;
0883   fs << strand1.str() << std::endl;
0884   fs << bp1.str() << std::endl;
0885   fs << bp2.str() << std::endl;
0886   fs << strand2.str() << std::endl;
0887   fs.close();
0888   delete classification;
0889 }
0890 
0891 const char* DamageRecord::GetChar(const G4bool& direct, const G4bool& indirect, const G4double& e)
0892 {
0893   if (direct && indirect) {
0894     return fBothDamageChar;
0895   }
0896   else if (direct) {
0897     return fDirectDamageChar;
0898   }
0899   else if (indirect) {
0900     return fIndirectDamageChar;
0901   }
0902   else if (e > 0) {
0903     return fHitNoDamageChar;
0904   }
0905   else {
0906     return fNotHitChar;
0907   }
0908 }
0909 
0910 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0911 
0912 DamageClassification* DamageRecord::GetClassification(const G4double& dsbDistance)
0913 {
0914   /* Explanation of class because we do some sneaky bit-level manipulation to
0915      track breaks along the strand.
0916 
0917      Things to note straight off
0918      1) This routine makes no promise to count SSBs and DSBs correctly
0919       once they cease being relevant to the break complexity. That is, once
0920       a DSB has occurred, the SSB count becomes irrelevant
0921      2) The previous ten base pairs are kept in the ints lastTenStrandX
0922       Each bit in the number represents one base pair precedent. The bit
0923       is set if the strand was broken.
0924      3) The countX variables are used to track the number of breaks in the
0925       preceeding 10 base pairs on each strand. This is needed to calculate
0926       DSB+ events
0927      4) The lambda function count_bytes is an algorithm to count the number
0928       of set bytes in constant time for 32-bit integers. It yields for
0929       us the count of breaks in the previous ten bases.
0930 
0931      Classification of complexity:
0932      DSB++  2DSBs in fragment
0933      DSB+   1 DSB and then an SSB within 10 base pairs
0934      DSB  Each strand hit once within <=10 base pairs
0935      2SSB   Each strand hit once with  >10 base pairs separation
0936      SSB+   One strand hit twice or more
0937      SSB  One strand hit only once
0938      None   Nothing hit
0939 
0940      Classification by Source
0941 
0942   */
0943   // lambda function to count the number of bytes set in a 32-bit uint
0944   // https://blogs.msdn.microsoft.com/jeuge/2005/06/08/bit-fiddling-3/
0945   // Note: integers beginning with 0 are entered as octal!
0946   auto count_bytes = [](uint32_t u) {
0947     uint32_t uCount = u - ((u >> 1) & 033333333333) - ((u >> 2) & 011111111111);
0948     return ((uCount + (uCount >> 3)) & 030707070707) % 63;
0949   };
0950 
0951   auto classification = new DamageClassification();
0952 
0953   G4int nDSB = 0;
0954   G4int nSSB = 0;
0955   // Source classification
0956   G4int oldnDSB = 0;
0957   G4int nDSBm = 0;
0958   G4int nDSBi = 0;
0959   G4int nDSBd = 0;
0960   G4int nDSBh = 0;
0961 
0962   // Complexity Classification
0963   G4int nDSBPlus = 0;
0964   G4bool SSB1 = false;
0965   G4bool SSB2 = false;
0966   // Counters for the last ten bp
0967   uint32_t lastTenDirectStrand1 = 0;
0968   uint32_t lastTenDirectStrand2 = 0;
0969   uint32_t lastTenIndirectStrand1 = 0;
0970   uint32_t lastTenIndirectStrand2 = 0;
0971   uint32_t lastTenStrand1 = 0;
0972   uint32_t lastTenStrand2 = 0;
0973   uint32_t lastTenTracked1 = 0;  // We remove entries from these if they have
0974   uint32_t lastTenTracked2 = 0;  // been counted in DSBs
0975   uint32_t count1 = 0;
0976   uint32_t count2 = 0;
0977   G4bool strand1Indirect = false;
0978   G4bool strand1Direct = false;
0979   G4bool strand2Indirect = false;
0980   G4bool strand2Direct = false;
0981   G4int baseDamage = 0;
0982   G4int strandDamage = 0;
0983   uint32_t truncator = std::pow(2, dsbDistance + 1) - 1;
0984   for (G4int ii = 0; ii != (G4int)fDamageRecords.size(); ii++) {
0985     lastTenDirectStrand1 = (lastTenDirectStrand1 << 1) & truncator;
0986     lastTenDirectStrand2 = (lastTenDirectStrand2 << 1) & truncator;
0987     lastTenIndirectStrand1 = (lastTenIndirectStrand1 << 1) & truncator;
0988     lastTenIndirectStrand2 = (lastTenIndirectStrand2 << 1) & truncator;
0989 
0990     lastTenStrand1 = lastTenStrand1 << 1;
0991     lastTenStrand2 = lastTenStrand2 << 1;
0992     lastTenTracked1 = lastTenTracked1 << 1;
0993     lastTenTracked2 = lastTenTracked2 << 1;
0994     // DSB distance by default is 10
0995     // 2047 = 0b11111111111, ie. truncation to 11bp
0996     lastTenStrand1 = lastTenStrand1 & truncator;
0997     lastTenStrand2 = lastTenStrand2 & truncator;
0998     lastTenTracked1 = lastTenTracked1 & truncator;
0999     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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!