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 AnalysisHandler.cc 28 /// \brief Implementation of the AnalysisHandl 29 30 #include "AnalysisHandler.hh" 31 #include "ScanDamage.hh" 32 #include "DamageClassifier.hh" 33 #include "ParametersParser.hh" 34 #include "SDDData.hh" 35 36 #include <cmath> 37 //....oooOO0OOooo........oooOO0OOooo........oo 38 39 AnalysisHandler::AnalysisHandler(): pTLKDoseMa 40 { 41 fScanDamage = std::make_unique<ScanDamage> 42 fTLKModel = std::make_unique<TLKModel>(); 43 fLEMIVModel = std::make_unique<LEMIVModel> 44 fBelovModel = std::make_unique<BelovModel> 45 fBpForDSB = 10; 46 47 if (ParametersParser::Instance()->GetThres 48 auto e = std::stod(ParametersParser::I 49 if (e > 0) fScanDamage->SetThresholdEn 50 } 51 if (ParametersParser::Instance()->GetProba 52 auto p = std::stod(ParametersParser::I 53 if (p > 0 && p <= 100) fScanDamage->Se 54 } 55 } 56 57 //....oooOO0OOooo........oooOO0OOooo........oo 58 59 void AnalysisHandler::SetThresholdEnergy(doubl 60 { 61 fScanDamage->SetThresholdEnergy(e); 62 } 63 64 //....oooOO0OOooo........oooOO0OOooo........oo 65 66 void AnalysisHandler::GetAllDamageAndScanSB() 67 { 68 std::map<unsigned int,std::map<unsigned in 69 if (ParametersParser::Instance()->WannaLoa 70 SDDData sdddata(ParametersParser::Inst 71 dmMap = sdddata.GetAllDamage(); 72 fDose = sdddata.GetDose(); 73 fChromosomeBpMap = sdddata.GetChromoso 74 } else { 75 if (ParametersParser::Instance()->Wann 76 fScanDamage->SkipScanningIndirectD 77 } 78 dmMap = fScanDamage->ExtractDamage(); 79 fEdepInNucleus = fScanDamage->GetEdepS 80 double nuclesumass = fScanDamage->GetN 81 double eVtoJ = 1.60E-19; 82 fDose = fEdepInNucleus*eVtoJ/nuclesuma 83 fNBp = fScanDamage->GetTotalNbBpPlaced 84 fChromosomeBpMap = fScanDamage->GetChr 85 } 86 DamageClassifier damClass; 87 std::map<int,int> ndsbMap, ncdsbMap, nssbM 88 std::map<int,int> ndirsbMap, ndsbdirMap, n 89 for (const auto& [chromo,evtDm] : dmMap) { 90 for (const auto& [evt, dmV] : evtDm) { 91 fAllDamage.insert(fAllDamage.end() 92 std::vector<Damage> tmpV{dmV}; 93 auto classifiedDamage = damClass.M 94 95 for (auto dm : dmV) { 96 if (dm.GetDamageType() == Dama 97 if ( (nsbMap.find(evt) == 98 nsbMap.insert({evt,1}) 99 } else { 100 nsbMap[evt] ++; 101 } 102 if (dm.GetCause() == Damag 103 if ( (ndirsbMap.find(e 104 ndirsbMap.insert({ 105 } else { 106 ndirsbMap[evt] ++; 107 } 108 } 109 } 110 } 111 112 int NumDSBForThisCluster = damClas 113 if ( (ndsbMap.find(evt) == ndsbMap 114 if (NumDSBForThisCluster > 0) 115 } else { 116 ndsbMap[evt] += NumDSBForThisC 117 } 118 119 int NumcDSBForThisCluster = damCla 120 if ( (ncdsbMap.find(evt) == ncdsbM 121 if (NumcDSBForThisCluster > 0) 122 } else { 123 ncdsbMap[evt] += NumcDSBForThi 124 } 125 126 int NumSSBForThisCluster = damClas 127 if ( (nssbMap.find(evt) == nssbMap 128 if (NumSSBForThisCluster > 0) 129 } else { 130 nssbMap[evt] += NumSSBForThisC 131 } 132 133 int NumDSBdirForThisCluster = damC 134 if ( (ndsbdirMap.find(evt) == ndsb 135 if (NumDSBdirForThisCluster > 136 } else { 137 ndsbdirMap[evt] += NumDSBdirFo 138 } 139 140 int NumDSBInForThisCluster = damCl 141 if ( (ndsbInMap.find(evt) == ndsbI 142 if (NumDSBInForThisCluster > 0 143 } else { 144 ndsbInMap[evt] += NumDSBInForT 145 } 146 147 int NumDSBdirInForThisCluster = da 148 if ( (ndsbdirIMap.find(evt) == nds 149 if (NumDSBdirInForThisCluster 150 } else { 151 ndsbdirIMap[evt] += NumDSBdirI 152 } 153 } 154 } 155 156 // DSB and its error: 157 float xxtotal=0, rms, xtotal = 0; 158 if (ndsbMap.size() >0) { 159 for (auto const& [evt,numdsb] : ndsbMa 160 xtotal += (float)numdsb; 161 xxtotal += float(numdsb*numdsb); 162 } 163 if (ndsbMap.size() == 1) { 164 // try to estimate error using poi 165 rms = std::sqrt(xtotal); 166 } else rms = std::sqrt(std::fabs(xxtot 167 fNDSBandError.first = xtotal; 168 fNDSBandError.second = rms; 169 } 170 171 172 // cDSB and its error: 173 if (ncdsbMap.size() > 0) { 174 xxtotal=0, rms = 0, xtotal = 0; 175 for (auto const& [evt,numcdsb] : ncdsb 176 xtotal += (float)numcdsb; 177 xxtotal += float(numcdsb*numcdsb); 178 } 179 if (ncdsbMap.size() == 1) { 180 // try to estimate error using poi 181 rms = std::sqrt(xtotal); 182 } else rms = std::sqrt(std::fabs(xxtot 183 fNcDSBandError.first = xtotal; 184 fNcDSBandError.second = rms; 185 } 186 // sDSB and its error, using error propaga 187 if (fNDSBandError.first > 0) { 188 fNsDSBandError.first = fNDSBandError.f 189 if (fNcDSBandError.first > 0) { 190 fNsDSBandError.second = (fNsDSBand 191 (fNDSBandError.second/fNDSBand 192 (fNcDSBandError.second/fNcDSBa 193 } 194 else fNsDSBandError.second = (fNsDSBan 195 } 196 197 // DSBdir and its error: 198 if (ndsbdirMap.size() > 0) { 199 xxtotal=0, rms = 0, xtotal = 0; 200 for (auto const& [evt,numdsbdir] : nds 201 xtotal += (float)numdsbdir; 202 xxtotal += float(numdsbdir*numdsbd 203 } 204 if (ndsbdirMap.size() == 1) { 205 // try to estimate error using poi 206 rms = std::sqrt(xtotal); 207 } else rms = std::sqrt(std::fabs(xxtot 208 fNDSBdirandError.first = xtotal; 209 fNDSBdirandError.second = rms; 210 } 211 212 // DSBIn and its error: 213 if (ndsbInMap.size() > 0) { 214 xxtotal=0, rms = 0, xtotal = 0; 215 for (auto const& [evt,numdsbIn] : ndsb 216 xtotal += (float)numdsbIn; 217 xxtotal += float(numdsbIn*numdsbIn 218 } 219 if (ndsbInMap.size() == 1) { 220 // try to estimate error using poi 221 rms = std::sqrt(xtotal); 222 } else rms = std::sqrt(std::fabs(xxtot 223 fNDSBIndandError.first = xtotal; 224 fNDSBIndandError.second = rms; 225 } 226 227 // DSBdirIn and its error: 228 if (ndsbdirIMap.size() > 0) { 229 xxtotal=0, rms = 0, xtotal = 0; 230 for (auto const& [evt,numdsbdirIn] : n 231 xtotal += (float)numdsbdirIn; 232 xxtotal += float(numdsbdirIn*numds 233 } 234 if (ndsbdirIMap.size() == 1) { 235 // try to estimate error using poi 236 rms = std::sqrt(xtotal); 237 } else rms = std::sqrt(std::fabs(xxtot 238 fNDSBdirIandError.first = xtotal; 239 fNDSBdirIandError.second = rms; 240 } 241 242 243 // SSB and its error: 244 if (nssbMap.size() > 0) { 245 xxtotal=0, rms = 0, xtotal = 0; 246 for (auto const& [evt,numssb] : nssbMa 247 xtotal += (float)numssb; 248 xxtotal += float(numssb*numssb); 249 } 250 if (nssbMap.size() == 1) { 251 // try to estimate error using poi 252 rms = std::sqrt(xtotal); 253 } else rms = std::sqrt(std::fabs(xxtot 254 fNSSBandError.first = xtotal; 255 fNSSBandError.second = rms; 256 } 257 258 // SB and its error: 259 if (nsbMap.size() > 0) { 260 xxtotal=0, rms = 0, xtotal = 0; 261 for (auto const& [evt,numsb] : nsbMap) 262 xtotal += (float)numsb; 263 xxtotal += float(numsb*numsb); 264 } 265 if (nsbMap.size() == 1) { 266 // try to estimate error using poi 267 rms = std::sqrt(xtotal); 268 } else rms = std::sqrt(std::fabs(xxtot 269 fNSBandError.first = xtotal; 270 fNSBandError.second = rms; 271 } 272 273 // direct SB and its error: 274 if (ndirsbMap.size() > 0) { 275 xxtotal=0, rms = 0, xtotal = 0; 276 for (auto const& [evt,numdirsb] : ndir 277 xtotal += (float)numdirsb; 278 xxtotal += float(numdirsb*numdirsb 279 } 280 if (ndirsbMap.size() == 1) { 281 // try to estimate error using poi 282 rms = std::sqrt(xtotal); 283 } else rms = std::sqrt(std::fabs(xxtot 284 fNdirSBandError.first = xtotal; 285 fNdirSBandError.second = rms; 286 } 287 288 // indirect SB its error, using error prop 289 if (fNSBandError.first > 0) { 290 fNindirSBandError.first = fNSBandError 291 if (fNdirSBandError.first > 0) { 292 fNindirSBandError.second = (fNindi 293 (fNSBandError.second/fNSBandEr 294 (fNdirSBandError.second/fNdirS 295 } 296 else fNindirSBandError.second = (fNind 297 } 298 299 // clear Maps 300 ndsbMap.clear(); 301 ncdsbMap.clear(); 302 nssbMap.clear(); 303 nsbMap.clear(); 304 ndirsbMap.clear(); 305 ndsbdirMap.clear(); 306 ndsbdirIMap.clear(); 307 ndsbInMap.clear(); 308 dmMap.clear(); 309 310 fIsSBScanned = true; 311 } 312 313 //....oooOO0OOooo........oooOO0OOooo........oo 314 315 void AnalysisHandler::GiveMeSBs() 316 { 317 if (!fIsSBScanned) GetAllDamageAndScanSB() 318 std::string outName = ParametersParser::In 319 std::fstream file; 320 file.open(outName.c_str(), std::ios_base::ou 321 double norm = 1.0; 322 std::string normunit = ""; 323 if (ParametersParser::Instance()->GetUnitT 324 norm = 1.0/fDose; 325 normunit = "[SB/Gy]"; 326 } else { 327 double BbToGb = 1e-9; // convert Bb to 328 norm = 1.0/(fDose*fNBp*BbToGb); 329 normunit = "[SB/Gy/Gbp]"; 330 } 331 file <<"#=========================== Stran 332 file <<"Name of Cell Nucleus: "<<Parameter 333 if (ParametersParser::Instance()->WannaLoa 334 std::string outstrtmp = "No info from 335 file <<"Volume of Cell Nucleus: "<<out 336 file <<"Mass Density of Cell Nucleus: 337 file <<"Mass of Cell Nucleus: "<<outst 338 file <<"Energy deposited in Cell Nucle 339 file <<"Dose delivered in Cell Nucleus 340 file <<"Minimum Distance between two c 341 file <<"Number of basepairs in Cell Nu 342 file <<"Threshold Energy for direct da 343 file <<"Propability for indirect damag 344 } else { 345 file <<"Volume of Cell Nucleus: "<<fSc 346 file <<"Mass Density of Cell Nucleus: 347 file <<"Mass of Cell Nucleus: "<<fScan 348 file <<"Energy deposited in Cell Nucle 349 file <<"Dose delivered in Cell Nucleus 350 file <<"Minimum Distance between two c 351 file <<"Number of basepairs in Cell Nu 352 file <<"Threshold Energy for direct da 353 if (fScanDamage->SkippedScanningIndire 354 355 else file <<"Propability for indirect 356 <<fScanDamage->GetProbabilityF 357 } 358 359 file <<"#================================= 360 file << "\n"; 361 file <<"#Un-normalized results:\n"; 362 file << "TotalSB [SB] \t" << fNSBandErr 363 file << "DirSB [SB] \t" << fNdirSBandEr 364 file << "IndirSB [SB] \t" << fNindirSBand 365 file << "SSB [SB] \t" << fNSSBandEr 366 file << "DSB [SB] \t" << fNDSBandErro 367 file << "cDSB [SB] \t" << fNcDSBandErr 368 file << "sDSB [SB] \t" << fNsDSBandE 369 file << "DSBdir [SB] \t" << fNDSBdiran 370 file << "DSBind [SB] \t" << fNDSBIndan 371 file << "DSBdirIn [SB] \t" << fNDSBdirIa 372 file << "\n"; 373 file <<"#Normalized results:\n"; 374 file << "TotalSB " + normunit +" \t" < 375 file << "DirSB " + normunit +" \t" << 376 file << "IndirSB " + normunit +" \t" << 377 file << "SSB " + normunit +" \t" < 378 file << "DSB " + normunit +" \t" << 379 file << "cDSB " + normunit +" \t" << 380 file << "sDSB " + normunit +" \t" < 381 file << "DSBdir " + normunit +" \t" < 382 file << "DSBind " + normunit +" \t" < 383 file << "DSBdirIn " + normunit +" \t" < 384 file <<"#================================= 385 file << "where: \n"; 386 file << "-----> TotalSB: Total strand-brea 387 file << "-----> DirSB: Direct strand-break 388 file << "-----> IndirSB: Indirect strand-b 389 file << "-----> SSB: Single strand-breaks\ 390 file << "-----> DSB: Double strand-breaks\ 391 file << "-----> cDSB: Complex DSB\n"; 392 file << "-----> sDSB: Simple DSB\n"; 393 file << "-----> DSBdir: DSB that contains 394 file << "-----> DSBdind: DSB that contains 395 file << "-----> DSBdirIn: DSB that contain 396 file <<"#============================== En 397 file.close(); 398 } 399 400 //....oooOO0OOooo........oooOO0OOooo........oo 401 402 void AnalysisHandler::ApplyDNAModel(std::strin 403 { 404 if (!fIsSBScanned) GetAllDamageAndScanSB() 405 406 if (dnaModel == "TLK") { 407 std::cout << "Invoking TLK Model" << s 408 //fTLKModel->SetDose(fDose); 409 SetParametersForTLKModel(); 410 //fTLKModel->ComputeAndSetDamageInput( 411 fTLKModel->SetSingleDSBYield(fNsDSBand 412 fTLKModel->SetComplexDSBYield(fNcDSBan 413 if (ParametersParser::Instance()->GetT 414 auto val = std::stod(ParametersPar 415 if (val != pTLKDoseMax) pTLKDoseMa 416 } 417 if (ParametersParser::Instance()->GetT 418 auto val = std::stod(ParametersPar 419 if (val != pTLKDeltaDose) pTLKDelt 420 } 421 fTLKModel->CalculateRepair(pTLKDoseMax 422 std::string outname = "TLK_"+Parameter 423 fTLKModel->WriteOutput(outname); 424 } 425 426 if (dnaModel == "LEMIV") { 427 std::cout << "Invoking LEMIV Model" << 428 fLEMIVModel->SetChromosomeBpSizesMap(f 429 fLEMIVModel->SetDose(fDose); 430 SetParametersForLEMIVModel(); 431 fLEMIVModel->ComputeAndSetDamageInput( 432 if (ParametersParser::Instance()->GetL 433 auto val = std::stod(ParametersPar 434 if (val != pLEMIVtimeMax) pLEMIVti 435 } 436 if (ParametersParser::Instance()->GetL 437 auto val = std::stod(ParametersPar 438 if (val != pLEMIVdeltaTime) pLEMIV 439 } 440 fLEMIVModel->CalculateRepair(pLEMIVtim 441 std::string outname = "LEMIV_"+Paramet 442 fLEMIVModel->WriteOutput(outname); 443 } 444 445 if (dnaModel == "BELOV") { 446 std::cout << "Invoking Belov's Model" 447 fBelovModel->SetDSBandComDSBandDose(fN 448 if (ParametersParser::Instance()->GetB 449 auto Nirrep = std::stod(Parameters 450 fBelovModel->SetNirrep(Nirrep); 451 } 452 double Dz = 1.0; 453 if (ParametersParser::Instance()->GetB 454 Dz = std::stod(ParametersParser::I 455 } 456 fBelovModel->CalculateRepair(Dz); 457 std::string outname = "BELOV_"+Paramet 458 fBelovModel->WriteOutput(outname); 459 } 460 } 461 462 //....oooOO0OOooo........oooOO0OOooo........oo 463 464 void AnalysisHandler::CreateSDD(std::string fi 465 { 466 std::string str_tmp; 467 std::fstream ofile; 468 ofile.open(filename.c_str(), std::ios_base:: 469 if (ofile.is_open()) { 470 // Create header 471 ofile 472 <<"SDD version, SDDv1.0;\n" 473 <<"Software, dsbandrepair;\n" 474 <<"Author contact, Le Tuan Anh - a 475 <<"Simulation Details, DNA damages 476 <<"Source, ;\n" 477 <<"Source type, ;\n" 478 <<"Incident particles, "<<0<<";\n" 479 <<"Mean particle energy ("<<Parame 480 <<ParametersParser::Instance()->Ge 481 <<"Energy distribution, , ;\n" 482 <<"Particle fraction, 0;\n" 483 <<"Dose or fluence, 1, "<<fDose<<" 484 <<"Dose rate, 0;\n" 485 <<"Irradiation target, ;\n" 486 <<"Volumes, 0;\n"; 487 ofile<<"Chromosome sizes, "<<fChromoso 488 for (auto const& [chroID, nBps] :fChro 489 float nMBps = nBps*1E-6;// convert 490 ofile<<", "<<nMBps; 491 } 492 ofile<<";\n"; 493 ofile<<"DNA Density, 0;\n" 494 <<"Cell Cycle Phase, 0;\n" 495 <<"DNA Structure, 0;\n" 496 <<"In vitro / in vivo, ;\n" 497 <<"Proliferation status, ;\n" 498 <<"Microenvironment, 0, 0;\n" 499 <<"Damage definition, 0;\n" 500 <<"Time, 0;\n" 501 <<"Damage and primary count, "+std 502 <<"Data entries, 1, 0, 1, 1, 1, 1, 503 <<"Data field explaination, Field 504 <<"type [1]-ChromosomeID [3]-stran 505 <<"Cause (direct: [0]=0) (indirec 506 <<"\n" 507 <<"***EndOfHeader***;\n" 508 <<"\n"; 509 510 // Data Section 511 int prevEvt = -1; 512 for (auto &damage : fAllDamage) { 513 //Field 1 Calassification 514 int newEvtFlag = 0; // = 2 if new 515 if (prevEvt != damage.GetEvt()) { 516 newEvtFlag = 2; 517 prevEvt = damage.GetEvt(); 518 } 519 ofile<<newEvtFlag<<", "<<damage.Ge 520 //Field 3 Chromosome IDs 521 ofile<<damage.GetDamageChromatin() 522 //Field 4, Chromosome position 523 ofile<<damage.GetCopyNb()<<"; "; 524 //Field 5, Cause: Unknown = -1, Di 525 ofile<<damage.GetCause()<<", "<<0< 526 //Field 6, Damage types: 527 int firstval = 0, secval = 0; 528 if (damage.GetDamageType() == Dama 529 if (damage.GetDamageType() == Dama 530 ofile<<firstval<<", "<<secval<<", 531 ofile<<"\n"; 532 } 533 } 534 ofile.close(); 535 } 536 537 //....oooOO0OOooo........oooOO0OOooo........oo 538 539 void AnalysisHandler::SetBpForDSB(unsigned int 540 { 541 if (pVal == fBpForDSB) return; 542 fBpForDSB = pVal; 543 fTLKModel->SetBpForDSB(fBpForDSB); 544 fLEMIVModel->SetBpForDSB(fBpForDSB); 545 fBelovModel->SetBpForDSB(fBpForDSB); 546 } 547 548 //....oooOO0OOooo........oooOO0OOooo........oo 549 550 void AnalysisHandler::SetParametersForTLKModel 551 { 552 if (ParametersParser::Instance()->GetTLKLa 553 pLambda1 = std::stod(ParametersParser: 554 if (ParametersParser::Instance()->GetTLKLa 555 pLambda2 = std::stod(ParametersParser: 556 if (ParametersParser::Instance()->GetTLKBe 557 pBeta1 = std::stod(ParametersParser::I 558 if (ParametersParser::Instance()->GetTLKBe 559 pBeta2 = std::stod(ParametersParser::I 560 if (ParametersParser::Instance()->GetTLKEt 561 pEta = std::stod(ParametersParser::Ins 562 if (pBeta1 != fTLKModel->GetBeta1()) fTLKM 563 if (pBeta2 != fTLKModel->GetBeta2()) fTLKM 564 if (pLambda1 != fTLKModel->GetLambda1()) f 565 if (pLambda2 != fTLKModel->GetLambda2()) f 566 if (pEta != fTLKModel->GetEta()) fTLKModel 567 } 568 569 //....oooOO0OOooo........oooOO0OOooo........oo 570 571 void AnalysisHandler::SetParametersForLEMIVMod 572 { 573 if (ParametersParser::Instance()->GetEMIVL 574 pLoopLength = std::stod(ParametersPars 575 if (ParametersParser::Instance()->GetEMIVF 576 pFunrej = std::stod(ParametersParser:: 577 if (ParametersParser::Instance()->GetEMIVT 578 pTfast = std::stod(ParametersParser::I 579 if (ParametersParser::Instance()->GetEMIVT 580 pTslow = std::stod(ParametersParser::I 581 582 if (pLoopLength != fLEMIVModel->GetLoopLen 583 if (pFunrej != fLEMIVModel->GetFunrej()) f 584 if (pTfast != fLEMIVModel->GetTfast()) fLE 585 if (pTslow != fLEMIVModel->GetTslow()) fLE 586 } 587 588 //....oooOO0OOooo........oooOO0OOooo........oo