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:SteppingAction.cc 28 /// brief: 29 #include "SteppingAction.hh" 30 31 #include "ChromosomeMapper.hh" 32 #include "DNAGeometry.hh" 33 #include "DNAHit.hh" 34 #include "DetectorConstruction.hh" 35 #include "EventAction.hh" 36 #include "OctreeNode.hh" 37 #include "PlacementVolumeInfo.hh" 38 #include "UtilityFunctions.hh" 39 40 #include "G4AffineTransform.hh" 41 #include "G4LogicalVolume.hh" 42 #include "G4RunManager.hh" 43 #include "G4StepPoint.hh" 44 #include "G4ThreeVector.hh" 45 #include "G4VPhysicalVolume.hh" 46 #include "G4VTouchable.hh" 47 #include "Randomize.hh" 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 50 51 SteppingAction::SteppingAction(EventAction* ev 52 { 53 const auto* det = dynamic_cast<const Detecto 54 G4RunManager::GetRunManager()->GetUserDete 55 56 fDNAGeometry = det->GetDNAGeometry(); 57 } 58 59 //....oooOO0OOooo........oooOO0OOooo........oo 60 61 void SteppingAction::UserSteppingAction(const 62 { 63 G4double edep = aStep->GetTotalEnergyDeposit 64 if (edep == 0) { 65 return; 66 } 67 68 // Work out if we are in the Logical volume 69 const G4VTouchable* touch = aStep->GetPreSte 70 71 // remember the chem volume is the "real" vo 72 G4LogicalVolume* dnaVol = fDNAGeometry->GetD 73 G4int depth = touch->GetHistoryDepth(); 74 G4bool isInDNARegion = false; 75 76 if (depth >= 0) { 77 fEventAction->AddCellEdep(edep); // dousa 78 } 79 80 while ((depth >= 0) && !isInDNARegion) { 81 if (touch->GetVolume(depth) != nullptr) { 82 isInDNARegion = isInDNARegion || (touch- 83 } 84 else { 85 G4cout << "Null pointer volume in mother 86 << " looking for volume at depth 87 assert(depth < 0); 88 } 89 depth--; 90 } 91 92 if (isInDNARegion) { 93 ProcessDNARegionHit(aStep); 94 } 95 96 G4ThreeVector prepos = aStep->GetPreStepPoin 97 if (fDNAGeometry->GetChromosomeMapper()->Get 98 DoChromosomeRegionHit(prepos, edep, false) 99 } 100 101 if (aStep->GetTrack()->GetTrackID() == 1) { 102 fEventAction->SetPrimStopPos(prepos); 103 G4double tl = aStep->GetStepLength(); 104 fEventAction->AddTrackLengthCell(tl); 105 if (fDNAGeometry->GetChromosomeMapper()->G 106 fEventAction->AddTrackLengthChro(tl); 107 } 108 } 109 } 110 111 //....oooOO0OOooo........oooOO0OOooo........oo 112 113 void SteppingAction::ProcessDNARegionHit(const 114 { 115 G4double edep = aStep->GetTotalEnergyDeposit 116 const G4VTouchable* touch = aStep->GetPreSte 117 G4ThreeVector pos = 118 aStep->GetPreStepPoint()->GetPosition() 119 + G4UniformRand() 120 * (aStep->GetPostStepPoint()->GetPosit 121 G4VPhysicalVolume* pv = aStep->GetPreStepPoi 122 G4LogicalVolume* lv = pv->GetLogicalVolume() 123 124 std::vector<G4VPhysicalVolume*> molvec; 125 OctreeNode* parentNode; 126 G4AffineTransform transform = touch->GetHist 127 parentNode = fDNAGeometry->GetTopOctreeNode( 128 if (parentNode == nullptr) { // See if pare 129 parentNode = fDNAGeometry->GetTopOctreeNod 130 if (parentNode != nullptr) { // If mother 131 // These methods are unsafe if the pv is 132 // need to be in this if statement. 133 // lv = pv->GetMotherLogical(); 134 pv = touch->GetVolume(1); 135 transform = touch->GetHistory()->GetTran 136 } 137 } 138 139 G4ThreeVector localPos = transform.Transform 140 141 if (parentNode != nullptr) { // Get Molecul 142 molvec = parentNode->SearchOctree(localPos 143 } 144 145 // debug output for finding logical volumes. 146 if (fDNAGeometry->GetVerbosity() > 1) { 147 if (parentNode != nullptr) { 148 G4String lvname = pv->GetLogicalVolume() 149 G4cout << "Found octree for logical volu 150 } 151 else { 152 // Don't bother doing error printouts if 153 // is at the DNA mother volume level. 154 // Remember the Chemistry Volume is the 155 if (pv->GetLogicalVolume() != fDNAGeomet 156 G4String lvname = pv->GetLogicalVolume 157 G4String motherlvname; 158 if (lv != nullptr) { 159 motherlvname = lv->GetName(); 160 } 161 G4cout << "Could not find octree for l 162 << "Particle in LV: " << lvname 163 } 164 } 165 } 166 167 if (!molvec.empty()) { 168 // NOTE: This loop can be optimised to avo 169 G4VPhysicalVolume* closest = molvec[0]; 170 G4double dist = (closest->GetTranslation() 171 if (molvec.size() > 1) { 172 G4double newdist; 173 for (auto it = molvec.begin() + 1; it != 174 newdist = ((*it)->GetTranslation() - l 175 if (newdist < dist) { 176 dist = newdist; 177 closest = (*it); 178 } 179 } 180 } 181 G4bool isDNAHit = (dist < fDNAGeometry->Ge 182 if (isDNAHit) { 183 DoChromosomeDNAHit(pos, localPos, edep, 184 DoChromosomeRegionHit(pos, edep, true); 185 } 186 } 187 } 188 189 //....oooOO0OOooo........oooOO0OOooo........oo 190 191 void SteppingAction::DoChromosomeRegionHit(con 192 con 193 { 194 // TODO: chromosome should probably be based 195 // parent volume, not of current pos 196 G4String chromosome = fDNAGeometry->GetChrom 197 if (!(chromosome.empty())) { 198 fEventAction->AddChromosomeEdep(chromosome 199 } 200 } 201 202 //....oooOO0OOooo........oooOO0OOooo........oo 203 204 void SteppingAction::DoChromosomeDNAHit(const 205 const 206 const 207 { 208 G4int placementIdx = std::stoi(utility::Spli 209 std::vector<G4String> pvNameVec = utility::S 210 molecule mol = utility::GetMoleculeEnum(pvNa 211 G4int chainIdx = std::stoi(pvNameVec.at(1)); 212 G4int strandIdx = std::stoi(pvNameVec.at(2)) 213 int64_t baseIdx = std::stoll(pvNameVec.at(3) 214 215 // Convert to local chain and base index 216 chainIdx = fDNAGeometry->GetGlobalChain(plac 217 baseIdx += fDNAGeometry->GetStartIdx(placeme 218 219 if (fDNAGeometry->GetStrandsFlipped(placemen 220 strandIdx = (strandIdx + 1) % 2; 221 } 222 G4String chromosome = fDNAGeometry->GetChrom 223 const DNAHit* dnaHit = new DNAHit(mol, place 224 edep, dist 225 226 fEventAction->AddDNAHit(dnaHit); 227 } 228 229 //....oooOO0OOooo........oooOO0OOooo........oo 230