| 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 LEMIVModel.cc
28 /// \brief Implementation of the LEMIVModel cl
29
30 #include "LEMIVModel.hh"
31
32 #include "ClassifiedDamage.hh"
33 #include "Damage.hh"
34 #include "DamageClassifier.hh"
35
36 #include <cmath>
37 #include <iostream>
38 #include <fstream>
39 #include <algorithm>
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42
43 LEMIVModel::LEMIVModel(double pLoopLength,doub
44 double pNDSB,double pFunrej,double pTfast,doub
45 fLoopLength(pLoopLength),
46 fNi(pNi),
47 fNc(pNc),
48 fNDSB(pNDSB),
49 fFunrej(pFunrej),
50 fTfast(pTfast),
51 fTslow(pTslow)
52 {
53 fBpForDSB = 10;
54 fDefaultsChromosomeSizes={250, 250, 242, 242
55 182, 171, 171, 159, 159, 145, 145, 1
56 134, 135, 135, 133, 133, 114, 114, 1
57 102, 90, 90, 83, 83, 80, 80, 59, 59,
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 double LEMIVModel::ComputeUnrej(double pTime)
63 {
64 double lambdac = (fNDSB-fNi)/fNc;
65 double Ffast = fNi/fNDSB;
66 double Fslow = fNc*lambdac/fNDSB;
67
68 return Ffast*std::exp(-std::log(2)*pTime/fTf
69 (Fslow-fFunrej)*std::exp(-std::log(2)*pTime/
70 fFunrej;
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oo
74
75 void LEMIVModel::ComputeAndSetDamageInput(std:
76 {
77
78 double nidsb=0;
79 double ncdsb=0;
80 double ndsb=0;
81
82 DamageClassifier damclass = DamageClassifier
83
84 auto sortedDamage = damclass.SortDamageByChr
85 // Check chromosome sizes:
86 if (fChromosomeBpMap.size() == 0) {// using
87 std::cout<<"=====> LEMIV Calculation will
88 for (int ii=0;ii<fDefaultsChromosomeSizes.
89 unsigned long long int nBp = fDefaultsCh
90 fChromosomeBpMap.insert({ii,nBp});
91 }
92 }
93 // Loop on each chromosome
94 int i=0;
95 for(auto it=sortedDamage.begin();it!=sortedD
96 {
97 i++;
98 // Damage are now sorted by event, push al
99 std::vector<Damage> chromoDamage;
100 for(auto itt=it->second.begin();itt!=it->s
101 {
102 std::move(itt->second.begin(), itt->seco
103 }
104
105 // sort the list of damage by ascending bp
106 std::sort(chromoDamage.begin(), chromoDama
107 [](const Damage& a, const Damage& b) {
108 return a.GetCopyNb() < b.GetCopyNb();
109 });
110
111 // for each loop inside the chromosome
112 auto chromID = it->first;
113 if (fChromosomeBpMap.find(chromID) == fChr
114 std::cerr<<"**** Fatal Error *****"<<std
115 std::cerr<<"LEMIVModel::ComputeAndSetDam
116 <<chromID<<std::endl;
117 std::cerr<<"*************** *****"<<std:
118 exit(EXIT_FAILURE);
119 }
120 for(auto startLoop=0;startLoop<fChromosome
121 {
122
123 int n = GetDSBPerLoop(chromoDamage,start
124
125 ndsb+=n;
126
127 if(n==1)
128 nidsb+=1.0;
129 if(n>=2)
130 ncdsb+=1.0;
131 }
132 }
133
134 // Set in the model input parameters
135 SetNumDSB(ndsb/fDose);
136 SetNumDomainIsolated(nidsb/fDose);
137 SetNumDomainClustered(ncdsb/fDose);
138 }
139
140 //....oooOO0OOooo........oooOO0OOooo........oo
141
142 int LEMIVModel::GetDSBPerLoop(std::vector<Dama
143 {
144 // Start to fill a vector with damage having
145 std::vector<Damage> loopDamage;
146
147 for(int i=0;i<vecDamage.size();i++)
148 {
149 if((vecDamage[i].GetCopyNb()>startLoop)&&(
150 {
151 loopDamage.push_back(vecDamage[i]);
152 }
153 }
154
155 // Make cluster
156 DamageClassifier dam;
157 auto classifiedDamage = dam.MakeCluster(loop
158
159 // Return the number of DSB in this loop
160 return dam.GetNumDSB(classifiedDamage);
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oo
164
165 void LEMIVModel::CalculateRepair(double pTMax,
166 {
167 if (pTMax <= 0.) {
168 std::cout<<"LEMIVModel::CalculateRepair()
169 <<"Plese check the input macro file!!!
170 exit(0);
171 }
172 if (pDeltaT <= 0.) {
173 std::cout<<"LEMIVModel::CalculateRepair()
174 <<"Plese check the input macro file!!!
175 exit(0);
176 }
177 fUCurve.clear();
178
179 for(double time=0.;time<=pTMax;time+=pDeltaT
180 {
181 fUCurve.push_back(std::make_pair(time,Comp
182 }
183
184 }
185
186 //....oooOO0OOooo........oooOO0OOooo........oo
187
188 void LEMIVModel::WriteOutput(std::string pFile
189 {
190 std::fstream file;
191 file.open(pFileName.c_str(), std::ios_base::
192 //Header part
193 file <<"#===================================
194 file << " L
195 file << "#Number of DSBs: " << fNDSB * fDose
196 file << "#Number of domains with clustered D
197 file << "#Number of domains with isolated DS
198 file << "#Funrej = " << fFunrej << "\n";
199 file << "#Tfast = " << fTfast << " h-1
200 file << "#LoopLength (length of domain) = "
201 file <<"#===================================
202 file << "Time (h)\tU\n";
203 //End Header part
204 for(int i=0;i<fUCurve.size();i++)
205 {
206 file << fUCurve[i].first << "\t" << fUCurv
207 }
208
209 file.close();
210 }