| 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