| 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 SDDData.cc
28 /// \brief Implementation of the SDDData class
29
30 #include "SDDData.hh"
31 #include <iostream>
32 #include <fstream>
33 #include <sstream>
34 #include <map>
35
36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
37
38 SDDData::SDDData(std::string p_name):
39 filename_(p_name)
40 {
41 }
42
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44
45 SDDData::SDDHeader SDDData::ReadHeader()
46 {
47
48 std::ifstream file(filename_.c_str());
49
50 if(!file.is_open())
51 {
52 std::cout << "No file: " << filename_ << std::endl;
53 return header_;
54 }
55
56 //1-SDD version
57 ReadString(file,header_.sdd_version);
58 //2-Software
59 ReadString(file,header_.software);
60 //3-Author
61 ReadString(file,header_.author);
62 //4-Simulation details
63 ReadString(file,header_.sim_details);
64
65 //5- Source
66 ReadString(file,header_.src_details);
67 //6-Source type
68 ReadInt(file,header_.src_type);
69 //7-Incident particles
70 ReadInts(file,header_.src_pdg);
71 //8-Mean Particle energy
72 ReadDoubles(file,header_.src_energy);
73 //9-Energy distribution
74 ReadString(file,header_.energy_dist);
75 //10-Particle fraction
76 ReadDoubles(file,header_.part_fraction);
77 //11-Dose or fluence
78 ReadDoubles(file,header_.dose);
79 //12-Dose rate
80 ReadDouble(file,header_.dose_rate);
81
82 //13-Irradiation target
83 ReadString(file,header_.target);
84 //14-Volumes
85 ReadDoubles(file,header_.volumes);
86 //15-Chromosome sizes
87 ReadDoubles(file,header_.chromo_size);
88 //16-DNA density
89 ReadDouble(file,header_.dna_density);
90
91 //17-Cell cycle phase
92 ReadDoubles(file,header_.cell_cycle);
93
94 //18-DNA strcuture
95 ReadInts(file,header_.dna_struct);
96
97 //19- in vitro/in vivo
98 ReadInt(file,header_.vitro_vivo);
99
100 //20-Proliferation status
101 ReadString(file,header_.proliferation);
102
103 //21-Microenvironment
104 ReadDoubles(file,header_.microenv);
105
106 //22-Damage definition
107 ReadDoubles(file,header_.damage_def);
108
109 //23-Time
110 ReadDouble(file,header_.time);
111
112 //24-Damage and primary count
113 ReadInts(file,header_.damage_prim_count);
114
115 //25-Data entries
116 ReadBools(file,header_.entries);
117
118 //26-Additional information
119 ReadString(file,header_.info);
120
121 file.close();
122
123 return header_;
124 }
125
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
127
128 void SDDData::ParseData()
129 {
130 ReadHeader();
131
132 data_.clear();
133
134 std::ifstream file(filename_.c_str());
135
136 std::string line;
137
138 // Pass the header
139 while(std::getline(file,line))
140 {
141 if(line.find("EndOfHeader")!=std::string::npos)
142 break;
143 }
144
145 // Start to read the data
146
147 while(std::getline(file,line))
148 {
149 if(!line.empty())
150 ParseLineData(line);
151 }
152
153 file.close();
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157
158 void SDDData::ParseLineData(std::string& line)
159 {
160
161 SDDDamage dmg;
162
163 if(header_.entries[0])
164 ExtractInts(line,2,dmg.classification);
165 if(header_.entries[1])
166 ExtractDoubles(line,3,dmg.coordinates);
167 if(header_.entries[2])
168 ExtractInts(line,4,dmg.chromo_ID);
169 if(header_.entries[3])
170 ExtractDoubles(line,1,dmg.chromo_position);
171 if(header_.entries[4])
172 ExtractInts(line,3,dmg.cause);
173 if(header_.entries[5])
174 ExtractInts(line,3,dmg.types);
175 if(header_.entries[6])
176 ExtractInts(line,3,dmg.break_spec);
177 if(header_.entries[7])
178 ExtractInts(line,3,dmg.dna_seq);
179 if(header_.entries[8])
180 ExtractDoubles(line,1,dmg.lesion_time);
181 if(header_.entries[9])
182 ExtractInts(line,1,dmg.particles);
183 if(header_.entries[10])
184 ExtractDoubles(line,1,dmg.energy);
185 if(header_.entries[11])
186 ExtractDoubles(line,3,dmg.translation);
187 if(header_.entries[12])
188 ExtractDoubles(line,1,dmg.direction);
189 if(header_.entries[13])
190 ExtractDoubles(line,1,dmg.particle_time);
191
192 data_.push_back(dmg);
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196
197 std::map<unsigned int,std::map<unsigned int,std::vector<Damage> > > SDDData::GetAllDamage()
198 {
199 if(data_.size()<=0)
200 {
201 ParseData();
202 }
203
204 std::map<unsigned int,std::map<unsigned int,std::vector<Damage> > > fmDamage;
205
206 for(auto it=data_.begin();it!=data_.end();it++)
207 {
208
209 Damage::DamageType pType = Damage::DamageType::fOther;
210 unsigned int pChromo = -1;
211 unsigned int pEvt = -1;
212 unsigned int pStrand = -1;
213 unsigned long int pCopyNb = -1;
214 Position pPos(0,0,0);
215 Damage::DamageCause pCause = Damage::DamageCause::fUnknown;
216 Damage::DamageChromatin pChromatin = Damage::DamageChromatin::fUnspecified;
217
218 if(header_.entries[0])
219 {
220 pEvt = it->classification[1];
221 }
222 if(header_.entries[1])
223 {
224 pPos.setX(it->coordinates[0]);
225 pPos.setY(it->coordinates[1]);
226 pPos.setZ(it->coordinates[2]);
227 }
228 if(header_.entries[2])
229 {
230 switch(it->chromo_ID[0])
231 {
232 case 0:
233 pChromatin = Damage::DamageChromatin::fUnspecified;
234 break;
235 case 1:
236 pChromatin = Damage::DamageChromatin::fHeterochromatin;
237 break;
238 case 2:
239 pChromatin = Damage::DamageChromatin::fEuchromatin;
240 break;
241 case 3:
242 pChromatin = Damage::DamageChromatin::fFreeDNA;
243 break;
244 case 4:
245 pChromatin = Damage::DamageChromatin::fOtherDNA;
246 break;
247 default:
248 pChromatin = Damage::DamageChromatin::fUnspecified;
249 }
250 pChromo = it->chromo_ID[1];
251 pStrand = it->chromo_ID[3];
252 }
253 if(header_.entries[3])
254 {
255 pCopyNb = (unsigned long int)(it->chromo_position[0]);
256 }
257 if(header_.entries[4])
258 {
259 switch(it->cause[0])
260 {
261 case 0:
262 pCause = Damage::DamageCause::fDirect;
263 break;
264 case 1:
265 pCause = Damage::DamageCause::fIndirect;
266 break;
267 default:
268 pCause = Damage::DamageCause::fUnknown;
269
270 }
271 }
272 if(header_.entries[5])
273 {
274 if(it->types[0]>0)
275 pType = Damage::DamageType::fBase;
276 if(it->types[1]>0)
277 pType = Damage::DamageType::fBackbone;
278 }
279 Damage aDamage(pType,pChromo,pEvt,pStrand,pCopyNb,pPos,pCause,pChromatin);
280 auto chroPos = fmDamage.find(pChromo);
281 if (chroPos == fmDamage.end()) {
282 std::vector<Damage> dmv{aDamage};
283 std::map<unsigned int,std::vector<Damage> > evtDamages{{pEvt,dmv}};
284 fmDamage.insert({pChromo,evtDamages});
285 } else {
286 auto evtPos = chroPos->second.find(pEvt);
287 if (evtPos == chroPos->second.end()) {
288 std::vector<Damage> dmv{aDamage};
289 chroPos->second.insert({pEvt,dmv});
290 } else {
291 chroPos->second[pEvt].push_back(aDamage);
292 }
293 }
294 }
295
296 return fmDamage;
297 }
298
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
300
301 void SDDData::ExtractInts(std::string& strLine,int numInt,std::vector<int>& field)
302 {
303 std::string delimiter = ";";
304 std::string token;
305
306 size_t pos = strLine.find(delimiter);
307
308 if(pos!=std::string::npos)
309 {
310 token = strLine.substr(0,pos);
311 strLine.erase(0,pos+delimiter.length());
312 }
313 else
314 {
315 token = strLine;
316 strLine = "";
317 }
318
319 std::string value;
320 delimiter = ",";
321
322 for(int i=1;i<numInt;i++)
323 {
324 pos = token.find(delimiter);
325 value = token.substr(0,pos);
326 field.push_back(std::atoi(value.c_str()));
327 token.erase(0,pos+delimiter.length());
328 }
329
330 field.push_back(std::atoi(token.c_str()));
331 }
332
333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
334
335 void SDDData::ExtractDoubles(std::string& strLine,int numInt,std::vector<double>& field)
336 {
337 std::string delimiter = ";";
338 std::string token;
339
340 size_t pos = strLine.find(delimiter);
341
342 if(pos!=std::string::npos)
343 {
344 token = strLine.substr(0,pos);
345 strLine.erase(0,pos+delimiter.length());
346 }
347 else
348 {
349 token = strLine;
350 strLine = "";
351 }
352
353 std::string value;
354 delimiter = ",";
355
356 for(int i=1;i<numInt;i++)
357 {
358 pos = token.find(delimiter);
359 value = token.substr(0,pos);
360 field.push_back(std::atoi(value.c_str()));
361 token.erase(0,pos+delimiter.length());
362 }
363
364 field.push_back(std::stod(token.c_str()));
365 }
366
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
368
369 void SDDData::ReadString(std::ifstream& file,std::string& field)
370 {
371 std::string line;
372 std::string token;
373
374 std::getline(file,line);
375 std::istringstream ss(line);
376
377 std::getline(ss,token,',');
378 std::getline(ss,token,',');
379
380 field=token;
381 }
382
383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
384
385 void SDDData::ReadInt(std::ifstream& file,int& field)
386 {
387 std::string line;
388 std::string token;
389
390 std::getline(file,line);
391 line = line.substr(0, line.size()-1);
392
393 std::istringstream ss(line);
394
395 std::getline(ss,token,',');
396 std::getline(ss,token,',');
397
398 field=std::atoi(token.c_str());
399 }
400
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402
403 void SDDData::ReadInts(std::ifstream& file,std::vector<int>& field)
404 {
405 std::string line;
406 std::string token;
407
408 std::getline(file,line);
409 line = line.substr(0, line.size()-1);
410
411 std::istringstream ss(line);
412
413 std::getline(ss,token,',');
414
415 while(std::getline(ss,token,','))
416 field.push_back(std::atoi(token.c_str()));
417 }
418
419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
420
421 void SDDData::ReadDouble(std::ifstream& file,double& field)
422 {
423 std::string line;
424 std::string token;
425
426 std::getline(file,line);
427 line = line.substr(0, line.size()-1);
428
429 std::istringstream ss(line);
430
431 std::getline(ss,token,',');
432 std::getline(ss,token,',');
433
434 field=std::stod(token.c_str());
435 }
436
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
438
439 void SDDData::ReadDoubles(std::ifstream& file,std::vector<double>& field)
440 {
441 std::string line;
442 std::string token;
443
444 std::getline(file,line);
445 line = line.substr(0, line.size()-1);
446
447 std::istringstream ss(line);
448
449 std::getline(ss,token,',');
450
451 while(std::getline(ss,token,','))
452 field.push_back(std::stod(token.c_str()));
453 }
454
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456
457 void SDDData::ReadBools(std::ifstream& file,std::vector<bool>& field)
458 {
459 std::string line;
460 std::string token;
461
462 std::getline(file,line);
463 line = line.substr(0, line.size()-1);
464
465 std::istringstream ss(line);
466
467 std::getline(ss,token,',');
468
469 while(std::getline(ss,token,','))
470 field.push_back(std::stoi(token.c_str()));
471 }
472
473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474
475 double SDDData::GetDose()
476 {
477 return header_.dose[1];
478 }
479
480 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
481
482 std::map<int,unsigned long long int> SDDData::GetChromosomeBpSizesMap(double &sum)
483 {
484 sum=0;
485 std::map<int,unsigned long long int> chromap;
486 for (int i=1;i<header_.chromo_size.size(); i++) { // index 0 of header_.chromo_size is number of chromosomes
487 double nbp = header_.chromo_size[i]; // in Mbp
488 nbp *= 1e+6; // to bp
489 sum += nbp;
490 chromap.insert({(i-1),nbp});
491 }
492 return chromap;
493 }
494
495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
496