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 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........oo 37 38 SDDData::SDDData(std::string p_name): 39 filename_(p_name) 40 { 41 } 42 43 //....oooOO0OOooo........oooOO0OOooo........oo 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_ << s 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........oo 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:: 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........oo 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........oo 196 197 std::map<unsigned int,std::map<unsigned int,st 198 { 199 if(data_.size()<=0) 200 { 201 ParseData(); 202 } 203 204 std::map<unsigned int,std::map<unsigned int, 205 206 for(auto it=data_.begin();it!=data_.end();it 207 { 208 209 Damage::DamageType pType = Damage::DamageT 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::Damag 216 Damage::DamageChromatin pChromatin = Damag 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:: 234 break; 235 case 1: 236 pChromatin = Damage::DamageChromatin:: 237 break; 238 case 2: 239 pChromatin = Damage::DamageChromatin:: 240 break; 241 case 3: 242 pChromatin = Damage::DamageChromatin:: 243 break; 244 case 4: 245 pChromatin = Damage::DamageChromatin:: 246 break; 247 default: 248 pChromatin = Damage::DamageChromatin:: 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 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::fIndirec 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, 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 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(aDamag 292 } 293 } 294 } 295 296 return fmDamage; 297 } 298 299 //....oooOO0OOooo........oooOO0OOooo........oo 300 301 void SDDData::ExtractInts(std::string& strLine 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........oo 334 335 void SDDData::ExtractDoubles(std::string& strL 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........oo 368 369 void SDDData::ReadString(std::ifstream& file,s 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........oo 384 385 void SDDData::ReadInt(std::ifstream& file,int& 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........oo 402 403 void SDDData::ReadInts(std::ifstream& file,std 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........oo 420 421 void SDDData::ReadDouble(std::ifstream& file,d 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........oo 438 439 void SDDData::ReadDoubles(std::ifstream& file, 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........oo 456 457 void SDDData::ReadBools(std::ifstream& file,st 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........oo 474 475 double SDDData::GetDose() 476 { 477 return header_.dose[1]; 478 } 479 480 //....oooOO0OOooo........oooOO0OOooo........oo 481 482 std::map<int,unsigned long long int> SDDData:: 483 { 484 sum=0; 485 std::map<int,unsigned long long int> chromap 486 for (int i=1;i<header_.chromo_size.size(); i 487 double nbp = header_.chromo_size[i]; // in 488 nbp *= 1e+6; // to bp 489 sum += nbp; 490 chromap.insert({(i-1),nbp}); 491 } 492 return chromap; 493 } 494 495 //....oooOO0OOooo........oooOO0OOooo........oo 496