| 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 LEMIVModel.cc
28 /// \brief Implementation of the LEMIVModel class
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........oooOO0OOooo........oooOO0OOooo.....
42
43 LEMIVModel::LEMIVModel(double pLoopLength,double pNi, double pNc,
44 double pNDSB,double pFunrej,double pTfast,double pTslow):
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, 198, 198, 190, 190, 182,
55 182, 171, 171, 159, 159, 145, 145, 138, 138, 134,
56 134, 135, 135, 133, 133, 114, 114, 107, 107, 102,
57 102, 90, 90, 83, 83, 80, 80, 59, 59, 64, 64, 47, 47, 51, 51, 156, 57}; // in Mbp
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
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/fTfast)+
69 (Fslow-fFunrej)*std::exp(-std::log(2)*pTime/fTslow)+
70 fFunrej;
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75 void LEMIVModel::ComputeAndSetDamageInput(std::vector<Damage> vecDamage)
76 {
77
78 double nidsb=0;
79 double ncdsb=0;
80 double ndsb=0;
81
82 DamageClassifier damclass = DamageClassifier();
83
84 auto sortedDamage = damclass.SortDamageByChromo(vecDamage);
85 // Check chromosome sizes:
86 if (fChromosomeBpMap.size() == 0) {// using default values
87 std::cout<<"=====> LEMIV Calculation will Using Default Choromosome sizes"<<std::endl;
88 for (int ii=0;ii<fDefaultsChromosomeSizes.size();ii++) {
89 unsigned long long int nBp = fDefaultsChromosomeSizes[ii]*1E6; // convert MBp tp Bp
90 fChromosomeBpMap.insert({ii,nBp});
91 }
92 }
93 // Loop on each chromosome
94 int i=0;
95 for(auto it=sortedDamage.begin();it!=sortedDamage.end();it++)
96 {
97 i++;
98 // Damage are now sorted by event, push all the damage in the same vector
99 std::vector<Damage> chromoDamage;
100 for(auto itt=it->second.begin();itt!=it->second.end();itt++)
101 {
102 std::move(itt->second.begin(), itt->second.end(), std::back_inserter(chromoDamage));
103 }
104
105 // sort the list of damage by ascending bp
106 std::sort(chromoDamage.begin(), chromoDamage.end(),
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) == fChromosomeBpMap.end()) {
114 std::cerr<<"**** Fatal Error *****"<<std::endl;
115 std::cerr<<"LEMIVModel::ComputeAndSetDamageInput: Cannot find size info for chrom ID "
116 <<chromID<<std::endl;
117 std::cerr<<"*************** *****"<<std::endl;
118 exit(EXIT_FAILURE);
119 }
120 for(auto startLoop=0;startLoop<fChromosomeBpMap[chromID];startLoop+=fLoopLength)
121 {
122
123 int n = GetDSBPerLoop(chromoDamage,startLoop);
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........oooOO0OOooo........oooOO0OOooo......
141
142 int LEMIVModel::GetDSBPerLoop(std::vector<Damage> vecDamage,unsigned int startLoop)
143 {
144 // Start to fill a vector with damage having bp between startLopp and startLoop+2Mbp
145 std::vector<Damage> loopDamage;
146
147 for(int i=0;i<vecDamage.size();i++)
148 {
149 if((vecDamage[i].GetCopyNb()>startLoop)&&(vecDamage[i].GetCopyNb()<startLoop+fLoopLength))
150 {
151 loopDamage.push_back(vecDamage[i]);
152 }
153 }
154
155 // Make cluster
156 DamageClassifier dam;
157 auto classifiedDamage = dam.MakeCluster(loopDamage,fBpForDSB,false);
158
159 // Return the number of DSB in this loop
160 return dam.GetNumDSB(classifiedDamage);
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164
165 void LEMIVModel::CalculateRepair(double pTMax, double pDeltaT)
166 {
167 if (pTMax <= 0.) {
168 std::cout<<"LEMIVModel::CalculateRepair() wrong value for timeMax !!!\n"
169 <<"Plese check the input macro file!!!"<<std::endl;
170 exit(0);
171 }
172 if (pDeltaT <= 0.) {
173 std::cout<<"LEMIVModel::CalculateRepair() wrong value for deltaTime !!!\n"
174 <<"Plese check the input macro file!!!"<<std::endl;
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,ComputeUnrej(time)));
182 }
183
184 }
185
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
187
188 void LEMIVModel::WriteOutput(std::string pFileName)
189 {
190 std::fstream file;
191 file.open(pFileName.c_str(), std::ios_base::out);
192 //Header part
193 file <<"#============================================= LEMIV MODEL =============================================#\n";
194 file << " LEMIV Model, CalculateRepair with:\n";
195 file << "#Number of DSBs: " << fNDSB * fDose << " DSBs.\n";
196 file << "#Number of domains with clustered DSB, Nc = " << fNc * fDose << " domains.\n";
197 file << "#Number of domains with isolated DSB, Ni = " << fNi * fDose<< " domains.\n";
198 file << "#Funrej = " << fFunrej << "\n";
199 file << "#Tfast = " << fTfast << " h-1 " << "#Tslow = " << fTslow << " h-1\n";
200 file << "#LoopLength (length of domain) = " << fLoopLength<<" bp \n";
201 file <<"#========================================================================================================#\n";
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" << fUCurve[i].second << "\n";
207 }
208
209 file.close();
210 }