| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file AnalysisHandler.cc
28 /// \brief Implementation of the AnalysisHandler class
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........oooOO0OOooo........oooOO0OOooo.....
38
39 AnalysisHandler::AnalysisHandler(): pTLKDoseMax(6.0), pTLKDeltaDose(0.25), fNBp (-1)
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()->GetThresholdE() != "") {
48 auto e = std::stod(ParametersParser::Instance()->GetThresholdE());
49 if (e > 0) fScanDamage->SetThresholdEnergy(e);
50 }
51 if (ParametersParser::Instance()->GetProbabilityForIndirectSB() != "") {
52 auto p = std::stod(ParametersParser::Instance()->GetProbabilityForIndirectSB());
53 if (p > 0 && p <= 100) fScanDamage->SetProbabilityForIndirectSBSelection(p/100.);
54 }
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
58
59 void AnalysisHandler::SetThresholdEnergy(double e)
60 {
61 fScanDamage->SetThresholdEnergy(e);
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
65
66 void AnalysisHandler::GetAllDamageAndScanSB()
67 {
68 std::map<unsigned int,std::map<unsigned int,std::vector<Damage> > > dmMap;
69 if (ParametersParser::Instance()->WannaLoadDamagesFromSDD()) {
70 SDDData sdddata(ParametersParser::Instance()->GetSDDFileName());
71 dmMap = sdddata.GetAllDamage();
72 fDose = sdddata.GetDose();
73 fChromosomeBpMap = sdddata.GetChromosomeBpSizesMap(fNBp);
74 } else {
75 if (ParametersParser::Instance()->WannaSkipScanningIndirectDamage()) {
76 fScanDamage->SkipScanningIndirectDamage();
77 }
78 dmMap = fScanDamage->ExtractDamage();
79 fEdepInNucleus = fScanDamage->GetEdepSumInNucleus();//eV
80 double nuclesumass = fScanDamage->GetNucleusMass(); // kg
81 double eVtoJ = 1.60E-19;
82 fDose = fEdepInNucleus*eVtoJ/nuclesumass;
83 fNBp = fScanDamage->GetTotalNbBpPlacedInGeo();
84 fChromosomeBpMap = fScanDamage->GetChromosomeBpSizesMap();
85 }
86 DamageClassifier damClass;
87 std::map<int,int> ndsbMap, ncdsbMap, nssbMap, nsbMap;
88 std::map<int,int> ndirsbMap, ndsbdirMap, ndsbdirIMap, ndsbInMap;
89 for (const auto& [chromo,evtDm] : dmMap) {
90 for (const auto& [evt, dmV] : evtDm) {
91 fAllDamage.insert(fAllDamage.end(),dmV.begin(),dmV.end()); // to write SDD file and for LEM-IV
92 std::vector<Damage> tmpV{dmV};
93 auto classifiedDamage = damClass.MakeCluster(tmpV,fBpForDSB,false);
94
95 for (auto dm : dmV) {
96 if (dm.GetDamageType() == Damage::Damage::fBackbone ) {
97 if ( (nsbMap.find(evt) == nsbMap.end()) ) {
98 nsbMap.insert({evt,1});
99 } else {
100 nsbMap[evt] ++;
101 }
102 if (dm.GetCause() == Damage::Damage::fDirect) {
103 if ( (ndirsbMap.find(evt) == ndirsbMap.end()) ) {
104 ndirsbMap.insert({evt,1});
105 } else {
106 ndirsbMap[evt] ++;
107 }
108 }
109 }
110 }
111
112 int NumDSBForThisCluster = damClass.GetNumDSB(classifiedDamage);
113 if ( (ndsbMap.find(evt) == ndsbMap.end()) ) {
114 if (NumDSBForThisCluster > 0) ndsbMap.insert({evt,NumDSBForThisCluster});
115 } else {
116 ndsbMap[evt] += NumDSBForThisCluster;
117 }
118
119 int NumcDSBForThisCluster = damClass.GetNumComplexDSB(classifiedDamage);
120 if ( (ncdsbMap.find(evt) == ncdsbMap.end()) ) {
121 if (NumcDSBForThisCluster > 0) ncdsbMap.insert({evt,NumcDSBForThisCluster});
122 } else {
123 ncdsbMap[evt] += NumcDSBForThisCluster;
124 }
125
126 int NumSSBForThisCluster = damClass.GetNumSSB(classifiedDamage);
127 if ( (nssbMap.find(evt) == nssbMap.end()) ) {
128 if (NumSSBForThisCluster > 0) nssbMap.insert({evt,NumSSBForThisCluster});
129 } else {
130 nssbMap[evt] += NumSSBForThisCluster;
131 }
132
133 int NumDSBdirForThisCluster = damClass.GetNumDSBwithDirectDamage(classifiedDamage);
134 if ( (ndsbdirMap.find(evt) == ndsbdirMap.end()) ) {
135 if (NumDSBdirForThisCluster > 0) ndsbdirMap.insert({evt,NumDSBdirForThisCluster});
136 } else {
137 ndsbdirMap[evt] += NumDSBdirForThisCluster;
138 }
139
140 int NumDSBInForThisCluster = damClass.GetNumDSBwithIndirectDamage(classifiedDamage);
141 if ( (ndsbInMap.find(evt) == ndsbInMap.end()) ) {
142 if (NumDSBInForThisCluster > 0) ndsbInMap.insert({evt,NumDSBInForThisCluster});
143 } else {
144 ndsbInMap[evt] += NumDSBInForThisCluster;
145 }
146
147 int NumDSBdirInForThisCluster = damClass.GetNumDSBwithBothDirectIndirectDamage(classifiedDamage);
148 if ( (ndsbdirIMap.find(evt) == ndsbdirIMap.end()) ) {
149 if (NumDSBdirInForThisCluster > 0) ndsbdirIMap.insert({evt,NumDSBdirInForThisCluster});
150 } else {
151 ndsbdirIMap[evt] += NumDSBdirInForThisCluster;
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] : ndsbMap) {
160 xtotal += (float)numdsb;
161 xxtotal += float(numdsb*numdsb);
162 }
163 if (ndsbMap.size() == 1) {
164 // try to estimate error using poisson distribution
165 rms = std::sqrt(xtotal);
166 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbMap.size()));
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] : ncdsbMap) {
176 xtotal += (float)numcdsb;
177 xxtotal += float(numcdsb*numcdsb);
178 }
179 if (ncdsbMap.size() == 1) {
180 // try to estimate error using poisson distribution
181 rms = std::sqrt(xtotal);
182 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ncdsbMap.size()));
183 fNcDSBandError.first = xtotal;
184 fNcDSBandError.second = rms;
185 }
186 // sDSB and its error, using error propagation method:
187 if (fNDSBandError.first > 0) {
188 fNsDSBandError.first = fNDSBandError.first - fNcDSBandError.first;
189 if (fNcDSBandError.first > 0) {
190 fNsDSBandError.second = (fNsDSBandError.first)*std::sqrt(
191 (fNDSBandError.second/fNDSBandError.first)*(fNDSBandError.second/fNDSBandError.first) +
192 (fNcDSBandError.second/fNcDSBandError.first)*(fNcDSBandError.second/fNcDSBandError.first));
193 }
194 else fNsDSBandError.second = (fNsDSBandError.first)*(fNDSBandError.second/fNDSBandError.first);
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] : ndsbdirMap) {
201 xtotal += (float)numdsbdir;
202 xxtotal += float(numdsbdir*numdsbdir);
203 }
204 if (ndsbdirMap.size() == 1) {
205 // try to estimate error using poisson distribution
206 rms = std::sqrt(xtotal);
207 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbdirMap.size()));
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] : ndsbInMap) {
216 xtotal += (float)numdsbIn;
217 xxtotal += float(numdsbIn*numdsbIn);
218 }
219 if (ndsbInMap.size() == 1) {
220 // try to estimate error using poisson distribution
221 rms = std::sqrt(xtotal);
222 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbInMap.size()));
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] : ndsbdirIMap) {
231 xtotal += (float)numdsbdirIn;
232 xxtotal += float(numdsbdirIn*numdsbdirIn);
233 }
234 if (ndsbdirIMap.size() == 1) {
235 // try to estimate error using poisson distribution
236 rms = std::sqrt(xtotal);
237 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbdirIMap.size()));
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] : nssbMap) {
247 xtotal += (float)numssb;
248 xxtotal += float(numssb*numssb);
249 }
250 if (nssbMap.size() == 1) {
251 // try to estimate error using poisson distribution
252 rms = std::sqrt(xtotal);
253 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(nssbMap.size()));
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 poisson distribution
267 rms = std::sqrt(xtotal);
268 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(nsbMap.size()));
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] : ndirsbMap) {
277 xtotal += (float)numdirsb;
278 xxtotal += float(numdirsb*numdirsb);
279 }
280 if (ndirsbMap.size() == 1) {
281 // try to estimate error using poisson distribution
282 rms = std::sqrt(xtotal);
283 } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndirsbMap.size()));
284 fNdirSBandError.first = xtotal;
285 fNdirSBandError.second = rms;
286 }
287
288 // indirect SB its error, using error propagation method:
289 if (fNSBandError.first > 0) {
290 fNindirSBandError.first = fNSBandError.first - fNdirSBandError.first;
291 if (fNdirSBandError.first > 0) {
292 fNindirSBandError.second = (fNindirSBandError.first)*std::sqrt(
293 (fNSBandError.second/fNSBandError.first)*(fNSBandError.second/fNSBandError.first) +
294 (fNdirSBandError.second/fNdirSBandError.first)*(fNdirSBandError.second/fNdirSBandError.first));
295 }
296 else fNindirSBandError.second = (fNindirSBandError.first)*(fNSBandError.second/fNSBandError.first);
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........oooOO0OOooo........oooOO0OOooo.....
314
315 void AnalysisHandler::GiveMeSBs()
316 {
317 if (!fIsSBScanned) GetAllDamageAndScanSB();
318 std::string outName = ParametersParser::Instance()->GetOutputName();
319 std::fstream file;
320 file.open(outName.c_str(), std::ios_base::out);
321 double norm = 1.0;
322 std::string normunit = "";
323 if (ParametersParser::Instance()->GetUnitTypeOfNormalization() == 2) {
324 norm = 1.0/fDose;
325 normunit = "[SB/Gy]";
326 } else {
327 double BbToGb = 1e-9; // convert Bb to Gb
328 norm = 1.0/(fDose*fNBp*BbToGb);
329 normunit = "[SB/Gy/Gbp]";
330 }
331 file <<"#=========================== Strand Breaks ============================#\n";
332 file <<"Name of Cell Nucleus: "<<ParametersParser::Instance()->GetCellNucleusName()<<"\n";
333 if (ParametersParser::Instance()->WannaLoadDamagesFromSDD()) {
334 std::string outstrtmp = "No info from SDD file!!!\n";
335 file <<"Volume of Cell Nucleus: "<<outstrtmp;
336 file <<"Mass Density of Cell Nucleus: "<<outstrtmp;
337 file <<"Mass of Cell Nucleus: "<<outstrtmp;
338 file <<"Energy deposited in Cell Nucleus: "<<outstrtmp;
339 file <<"Dose delivered in Cell Nucleus: "<<fDose <<" (Gy)\n";
340 file <<"Minimum Distance between two clusters: "<<fBpForDSB <<" (bp)\n";
341 file <<"Number of basepairs in Cell Nucleus: "<<fNBp <<" (bp)\n";
342 file <<"Threshold Energy for direct damage selection: "<<outstrtmp;
343 file <<"Propability for indirect damage selection: "<<outstrtmp;
344 } else {
345 file <<"Volume of Cell Nucleus: "<<fScanDamage->GetNucleusVolume()<<" (m3)\n";
346 file <<"Mass Density of Cell Nucleus: "<<fScanDamage->GetNucleusMassDensity()<<" (kg/m3)\n";
347 file <<"Mass of Cell Nucleus: "<<fScanDamage->GetNucleusMass()<<" (kg)\n";
348 file <<"Energy deposited in Cell Nucleus: "<<fEdepInNucleus <<" (eV)\n";
349 file <<"Dose delivered in Cell Nucleus: "<<fDose <<" (Gy)\n";
350 file <<"Minimum Distance between two clusters: "<<fBpForDSB <<" (bp)\n";
351 file <<"Number of basepairs in Cell Nucleus: "<<fNBp <<" (bp)\n";
352 file <<"Threshold Energy for direct damage selection: "<<fScanDamage->GetThresholdEnergy() <<" (eV)\n";
353 if (fScanDamage->SkippedScanningIndirectDamage()) file <<"Propability for indirect SB selection: "
354 <<" Skipped the indirect analysis\n";
355 else file <<"Propability for indirect damage selection: "
356 <<fScanDamage->GetProbabilityForIndirectSBSelection()*100.<<" (%)\n";
357 }
358
359 file <<"#======================================================================#\n";
360 file << "\n";
361 file <<"#Un-normalized results:\n";
362 file << "TotalSB [SB] \t" << fNSBandError.first <<"\t+/-\t"<<fNSBandError.second<< "\n";
363 file << "DirSB [SB] \t" << fNdirSBandError.first <<"\t+/-\t"<<fNdirSBandError.second<< "\n";
364 file << "IndirSB [SB] \t" << fNindirSBandError.first <<"\t+/-\t"<<fNindirSBandError.second<< "\n";
365 file << "SSB [SB] \t" << fNSSBandError.first <<"\t+/-\t"<<fNSSBandError.second<< "\n";
366 file << "DSB [SB] \t" << fNDSBandError.first <<"\t+/-\t"<<fNDSBandError.second<< "\n";
367 file << "cDSB [SB] \t" << fNcDSBandError.first <<"\t+/-\t"<<fNcDSBandError.second<< "\n";
368 file << "sDSB [SB] \t" << fNsDSBandError.first <<"\t+/-\t"<<fNsDSBandError.second<< "\n";
369 file << "DSBdir [SB] \t" << fNDSBdirandError.first <<"\t+/-\t"<<fNDSBdirandError.second<< "\n";
370 file << "DSBind [SB] \t" << fNDSBIndandError.first <<"\t+/-\t"<<fNDSBIndandError.second<< "\n";
371 file << "DSBdirIn [SB] \t" << fNDSBdirIandError.first <<"\t+/-\t"<<fNDSBdirIandError.second<< "\n";
372 file << "\n";
373 file <<"#Normalized results:\n";
374 file << "TotalSB " + normunit +" \t" << fNSBandError.first * norm <<"\t+/-\t"<<fNSBandError.second * norm<< "\n";
375 file << "DirSB " + normunit +" \t" << fNdirSBandError.first * norm<<"\t+/-\t"<<fNdirSBandError.second * norm<< "\n";
376 file << "IndirSB " + normunit +" \t" << fNindirSBandError.first * norm<<"\t+/-\t"<<fNindirSBandError.second * norm<< "\n";
377 file << "SSB " + normunit +" \t" << fNSSBandError.first * norm<<"\t+/-\t"<<fNSSBandError.second * norm<< "\n";
378 file << "DSB " + normunit +" \t" << fNDSBandError.first * norm<<"\t+/-\t"<<fNDSBandError.second * norm<< "\n";
379 file << "cDSB " + normunit +" \t" << fNcDSBandError.first * norm<<"\t+/-\t"<<fNcDSBandError.second * norm<< "\n";
380 file << "sDSB " + normunit +" \t" << fNsDSBandError.first * norm<<"\t+/-\t"<<fNsDSBandError.second * norm<< "\n";
381 file << "DSBdir " + normunit +" \t" << fNDSBdirandError.first * norm<<"\t+/-\t"<<fNDSBdirandError.second * norm<< "\n";
382 file << "DSBind " + normunit +" \t" << fNDSBIndandError.first * norm<<"\t+/-\t"<<fNDSBIndandError.second * norm<< "\n";
383 file << "DSBdirIn " + normunit +" \t" << fNDSBdirIandError.first * norm<<"\t+/-\t"<<fNDSBdirIandError.second * norm<< "\n";
384 file <<"#======================================================================#\n";
385 file << "where: \n";
386 file << "-----> TotalSB: Total strand-breaks\n";
387 file << "-----> DirSB: Direct strand-breaks\n";
388 file << "-----> IndirSB: Indirect strand-breaks\n";
389 file << "-----> SSB: Single strand-breaks\n";
390 file << "-----> DSB: Double strand-breaks\n";
391 file << "-----> cDSB: Complex DSB\n";
392 file << "-----> sDSB: Simple DSB\n";
393 file << "-----> DSBdir: DSB that contains at least one direct SB\n";
394 file << "-----> DSBdind: DSB that contains at least one indirect SB\n";
395 file << "-----> DSBdirIn: DSB that contains at both direct and indirect SB\n";
396 file <<"#============================== End ===================================#\n";
397 file.close();
398 }
399
400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
401
402 void AnalysisHandler::ApplyDNAModel(std::string dnaModel)
403 {
404 if (!fIsSBScanned) GetAllDamageAndScanSB();
405
406 if (dnaModel == "TLK") {
407 std::cout << "Invoking TLK Model" << std::endl;
408 //fTLKModel->SetDose(fDose);
409 SetParametersForTLKModel();
410 //fTLKModel->ComputeAndSetDamageInput(fAllDamage);
411 fTLKModel->SetSingleDSBYield(fNsDSBandError.first/fDose);
412 fTLKModel->SetComplexDSBYield(fNcDSBandError.first/fDose);
413 if (ParametersParser::Instance()->GetTLKdoseMax() != "") {
414 auto val = std::stod(ParametersParser::Instance()->GetTLKdoseMax());
415 if (val != pTLKDoseMax) pTLKDoseMax = val;
416 }
417 if (ParametersParser::Instance()->GetTLKdeltaDose() != "") {
418 auto val = std::stod(ParametersParser::Instance()->GetTLKdeltaDose());
419 if (val != pTLKDeltaDose) pTLKDeltaDose = val;
420 }
421 fTLKModel->CalculateRepair(pTLKDoseMax,pTLKDeltaDose);
422 std::string outname = "TLK_"+ParametersParser::Instance()->GetOutputName();
423 fTLKModel->WriteOutput(outname);
424 }
425
426 if (dnaModel == "LEMIV") {
427 std::cout << "Invoking LEMIV Model" << std::endl;
428 fLEMIVModel->SetChromosomeBpSizesMap(fChromosomeBpMap);
429 fLEMIVModel->SetDose(fDose);
430 SetParametersForLEMIVModel();
431 fLEMIVModel->ComputeAndSetDamageInput(fAllDamage);
432 if (ParametersParser::Instance()->GetLEMtimeMax() != "") {
433 auto val = std::stod(ParametersParser::Instance()->GetLEMtimeMax());
434 if (val != pLEMIVtimeMax) pLEMIVtimeMax = val;
435 }
436 if (ParametersParser::Instance()->GetLEMdeltaTime() != "") {
437 auto val = std::stod(ParametersParser::Instance()->GetLEMdeltaTime());
438 if (val != pLEMIVdeltaTime) pLEMIVdeltaTime = val;
439 }
440 fLEMIVModel->CalculateRepair(pLEMIVtimeMax,pLEMIVdeltaTime);
441 std::string outname = "LEMIV_"+ParametersParser::Instance()->GetOutputName();
442 fLEMIVModel->WriteOutput(outname);
443 }
444
445 if (dnaModel == "BELOV") {
446 std::cout << "Invoking Belov's Model" << std::endl;
447 fBelovModel->SetDSBandComDSBandDose(fNDSBandError.first,fNcDSBandError.first,fDose);
448 if (ParametersParser::Instance()->GetBELOVNirrep() != "") {
449 auto Nirrep = std::stod(ParametersParser::Instance()->GetBELOVNirrep());
450 fBelovModel->SetNirrep(Nirrep);
451 }
452 double Dz = 1.0;
453 if (ParametersParser::Instance()->GetBELOVDz() != "") {
454 Dz = std::stod(ParametersParser::Instance()->GetBELOVDz());
455 }
456 fBelovModel->CalculateRepair(Dz);
457 std::string outname = "BELOV_"+ParametersParser::Instance()->GetOutputName();
458 fBelovModel->WriteOutput(outname);
459 }
460 }
461
462 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
463
464 void AnalysisHandler::CreateSDD(std::string filename)
465 {
466 std::string str_tmp;
467 std::fstream ofile;
468 ofile.open(filename.c_str(), std::ios_base::out);
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 - anh.letuan@irsn.fr, , ;\n"
475 <<"Simulation Details, DNA damages from direct and indirect effects;\n"
476 <<"Source, ;\n"
477 <<"Source type, ;\n"
478 <<"Incident particles, "<<0<<";\n"
479 <<"Mean particle energy ("<<ParametersParser::Instance()->GetEnergyUnit()<<"), "
480 <<ParametersParser::Instance()->GetParticleEnergy()<<";\n"
481 <<"Energy distribution, , ;\n"
482 <<"Particle fraction, 0;\n"
483 <<"Dose or fluence, 1, "<<fDose<<";\n"
484 <<"Dose rate, 0;\n"
485 <<"Irradiation target, ;\n"
486 <<"Volumes, 0;\n";
487 ofile<<"Chromosome sizes, "<<fChromosomeBpMap.size();
488 for (auto const& [chroID, nBps] :fChromosomeBpMap) {
489 float nMBps = nBps*1E-6;// convert from Bp to MBp
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::to_string(fAllDamage.size())+", 0;\n"
502 <<"Data entries, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0;\n"
503 <<"Data field explaination, Field 1: [1]-eventID, Field 3: [0]-Chromatin "
504 <<"type [1]-ChromosomeID [3]-strand, Field 4:chrom position (copynb), Field 5: "
505 <<"Cause (direct: [0]=0) (indirect: [0]=1), Field 6: Damage types (Base:[0]>0) (Backbone: [1]>0);\n"
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 event;
515 if (prevEvt != damage.GetEvt()) {
516 newEvtFlag = 2;
517 prevEvt = damage.GetEvt();
518 }
519 ofile<<newEvtFlag<<", "<<damage.GetEvt()<<"; ";
520 //Field 3 Chromosome IDs
521 ofile<<damage.GetDamageChromatin()<<", "<<damage.GetChromo()<<", "<<0<<", "<<damage.GetStrand()<<"; ";
522 //Field 4, Chromosome position
523 ofile<<damage.GetCopyNb()<<"; ";
524 //Field 5, Cause: Unknown = -1, Direct = 0, Indirect = 1
525 ofile<<damage.GetCause()<<", "<<0<<", "<<0<<"; ";
526 //Field 6, Damage types:
527 int firstval = 0, secval = 0;
528 if (damage.GetDamageType() == Damage::DamageType::fBase) firstval = 1;
529 if (damage.GetDamageType() == Damage::DamageType::fBackbone) secval = 1;
530 ofile<<firstval<<", "<<secval<<", "<<0<<"; ";
531 ofile<<"\n";
532 }
533 }
534 ofile.close();
535 }
536
537 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
538
539 void AnalysisHandler::SetBpForDSB(unsigned int pVal)
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........oooOO0OOooo........oooOO0OOooo.....
549
550 void AnalysisHandler::SetParametersForTLKModel(double pLambda1,double pLambda2, double pBeta1, double pBeta2,double pEta)
551 {
552 if (ParametersParser::Instance()->GetTLKLambda1() != "")
553 pLambda1 = std::stod(ParametersParser::Instance()->GetTLKLambda1());
554 if (ParametersParser::Instance()->GetTLKLambda2() != "")
555 pLambda2 = std::stod(ParametersParser::Instance()->GetTLKLambda2());
556 if (ParametersParser::Instance()->GetTLKBeta1() != "")
557 pBeta1 = std::stod(ParametersParser::Instance()->GetTLKBeta1());
558 if (ParametersParser::Instance()->GetTLKBeta2() != "")
559 pBeta2 = std::stod(ParametersParser::Instance()->GetTLKBeta2());
560 if (ParametersParser::Instance()->GetTLKEta() != "")
561 pEta = std::stod(ParametersParser::Instance()->GetTLKEta());
562 if (pBeta1 != fTLKModel->GetBeta1()) fTLKModel->SetBeta1(pBeta1);
563 if (pBeta2 != fTLKModel->GetBeta2()) fTLKModel->SetBeta2(pBeta2);
564 if (pLambda1 != fTLKModel->GetLambda1()) fTLKModel->SetLambda1(pLambda1);
565 if (pLambda2 != fTLKModel->GetLambda2()) fTLKModel->SetLambda2(pLambda2);
566 if (pEta != fTLKModel->GetEta()) fTLKModel->SetEta(pEta);
567 }
568
569 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
570
571 void AnalysisHandler::SetParametersForLEMIVModel(double pLoopLength,double pFunrej,double pTfast,double pTslow)
572 {
573 if (ParametersParser::Instance()->GetEMIVLoopLength() != "")
574 pLoopLength = std::stod(ParametersParser::Instance()->GetEMIVLoopLength());
575 if (ParametersParser::Instance()->GetEMIVFunrej() != "")
576 pFunrej = std::stod(ParametersParser::Instance()->GetEMIVFunrej());
577 if (ParametersParser::Instance()->GetEMIVTFast() != "")
578 pTfast = std::stod(ParametersParser::Instance()->GetEMIVTFast());
579 if (ParametersParser::Instance()->GetEMIVTSlow() != "")
580 pTslow = std::stod(ParametersParser::Instance()->GetEMIVTSlow());
581
582 if (pLoopLength != fLEMIVModel->GetLoopLength()) fLEMIVModel->SetLoopLength(pLoopLength);
583 if (pFunrej != fLEMIVModel->GetFunrej()) fLEMIVModel->SetFunrej(pFunrej);
584 if (pTfast != fLEMIVModel->GetTfast()) fLEMIVModel->SetTfast(pTfast);
585 if (pTslow != fLEMIVModel->GetTslow()) fLEMIVModel->SetTslow(pTslow);
586 }
587
588 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....