| 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 ScanDamage.cc
28 /// \brief Implementation of the ScanDamage class
29
30 #include "ScanDamage.hh"
31 #include "ParametersParser.hh"
32
33 #include "TSystemDirectory.h"
34 #include "TFile.h"
35 #include "TTree.h"
36 #include "TRandom.h"
37
38
39 #include <iostream>
40 #include <fstream>
41 #include <sstream>
42 #include <iostream>
43 #include <filesystem>
44 namespace fs = std::filesystem;
45 TRandom gRandomGen;
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
47
48 ScanDamage::ScanDamage(): fSkipScanningIndirectDamage(false)
49 {
50 fThresholdEnergy = 17.5; //eV
51 fEdepSumInNucleus = 0;
52 fProbabilityForIndirectSB = 0.40; // 40%
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
56
57 std::map<unsigned int,std::map<unsigned int,std::vector<Damage> > > ScanDamage::ExtractDamage(){
58 ReadCellandVoxelDefFilePaths();
59 fMergedTables.clear();
60 fDamage.clear();
61 RetrieveVoxelBp();
62 FillVoxelData();
63 ScanDamageFromPhys();
64 if (!fSkipScanningIndirectDamage) ScanDamageFromChem();
65 MergeDamageFromPhysChem();
66 for (const auto& [pChrom,table] : fMergedTables) {
67 unsigned int evt;
68 unsigned int strand;
69 ullint cpyNb;
70 unsigned int isBase;
71 double time;
72 double edep;
73 double origin;
74 Damage::DamageType pType;
75 Damage::DamageCause pOrigin;
76 std::map<unsigned int,std::vector<Damage> > perChromoDamage;
77 for (const auto& v : table) {
78 evt = v.at(0);
79 strand = v.at(1);
80 cpyNb = v.at(2);
81 isBase = v.at(3);
82 time = v.at(4);
83 edep = v.at(5);
84 origin = v.at(6);
85 if(isBase==0)
86 pType=Damage::DamageType::fBackbone;
87 else
88 pType=Damage::DamageType::fBase;
89 if(origin==0)
90 pOrigin=Damage::DamageCause::fDirect;
91 else
92 pOrigin=Damage::DamageCause::fIndirect;
93 if (perChromoDamage.find(evt) == perChromoDamage.end()) {
94 std::vector<Damage> dm{Damage(pType,pChrom,evt,strand,cpyNb,Position(0,0,0),
95 pOrigin,Damage::DamageChromatin::fUnspecified)};
96 perChromoDamage.insert({evt,dm});
97 } else perChromoDamage[evt].push_back(Damage(pType,pChrom,evt,strand,cpyNb,Position(0,0,0),
98 pOrigin,Damage::DamageChromatin::fUnspecified));
99 }
100 fDamage.insert({pChrom,perChromoDamage});
101 }
102 return fDamage;
103 }
104
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
106
107 void ScanDamage::RetrieveVoxelBp()
108 {
109 fBpPerVoxel.clear();
110
111 for (const auto & entry : fVoxelDefFilesList) {
112 std::ifstream file(entry);
113 if(!file.good() )
114 {
115 std::cerr<<"**** Fatal Error *****"<<std::endl;
116 std::cerr<<"ScanDamage::RetrieveVoxelBp: No file named "<<entry<<std::endl;
117 std::cerr<<"*************** *****"<<std::endl;
118 exit(EXIT_FAILURE);
119 }
120 std::string voxelName = "noName";
121 std::string line;
122 bool foundName = false;
123 bool foundNumOfBp = false;
124 while(std::getline(file, line)
125 && !foundNumOfBp)
126 {
127 std::istringstream iss(line);
128 std::string flag;
129 iss >> flag;
130 std::string charac;
131 iss >> charac;
132 // Look for the name of the voxel
133 if(flag=="_Name")
134 {
135 voxelName = charac;
136 foundName = true;
137 }
138 // Look for the flag "_Number"
139 // And the characteristic "voxelBasePair"
140 if(flag=="_Number" && charac=="voxelBasePair")
141 {
142 int numOfBp;
143 iss >> numOfBp;
144 if(!foundName)
145 {
146 std::cerr<<"*** Fatal Error ***"<<std::endl;
147 std::cerr<<"ScanDamage::RetrieveVoxelBp: The number of bp was found before the name "
148 <<"of the voxel... This is an unexpected case."<<std::endl;
149 std::cerr<<"******"<<std::endl;
150 exit(EXIT_FAILURE);
151 }
152 else
153 {
154 fBpPerVoxel[voxelName] = numOfBp;
155 std::cout<<voxelName<<" has "<<numOfBp<<" bp"<<std::endl;
156 foundNumOfBp = true;
157 }
158 }
159 }
160 file.close();
161 }
162
163 if (fBpPerVoxel.size() == 0) {
164 std::cerr<<"**** Fatal Error *****"<<std::endl;
165 std::cerr<<"ScanDamage::RetrieveVoxelBp: No Bp found in voxel definition files. \n Or make"
166 <<" sure that file imp.info exists in working directory!"<<std::endl;
167 std::cerr<<"*************** *****"<<std::endl;
168 exit(EXIT_FAILURE);
169 }
170 }
171
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
173
174 void ScanDamage::FillVoxelData()
175 {
176 std::ifstream file(fCellDefFilePath);
177 if(!file.good() )
178 {
179 std::cerr<<"**** Fatal Error *****"<<std::endl;
180 std::cerr<<"FillVoxelData: No file named "<<fCellDefFilePath<<std::endl;
181 std::cerr<<"*************** *****"<<std::endl;
182 exit(EXIT_FAILURE);
183 }
184 ullint bpCount = 0;
185 unsigned int voxelCount = 0;
186 int chromo_previous = 0;
187 // Read the file line by line
188 std::string line;
189 while(std::getline(file, line) )
190 {
191 std::istringstream iss(line);
192 std::string flag;
193 iss >> flag;
194 // If the flag correspond to the placement of a voxel
195 if(flag == "_pl")
196 {
197 std::string voxelName;
198 iss >> voxelName;
199 int chromo;
200 iss >> chromo;
201 int domain;
202 iss >> domain;
203 // If we change of chromosome then reset the number of bp.
204 // Each chromosome starts at 0 bp.
205 if(chromo != chromo_previous)
206 {
207 bpCount = 0;
208 chromo_previous = chromo;
209 }
210 // Fill the data structure
211 fVoxels.push_back( VoxelData(chromo, domain, bpCount) );
212 int numBpinthisvoxel = int(fBpPerVoxel[voxelName]);
213 bpCount += numBpinthisvoxel;
214 if (fChromosomeBpMap.find(chromo) == fChromosomeBpMap.end()) {
215 fChromosomeBpMap.insert({chromo,numBpinthisvoxel});
216 } else {
217 fChromosomeBpMap[chromo] += numBpinthisvoxel;
218 }
219 voxelCount++;
220 }
221 }
222 file.close();
223
224
225 if (fVoxels.size() == 0) {
226 std::cerr<<"**** Fatal Error *****"<<std::endl;
227 std::cerr<<"ScanDamage::FillVoxelData: NofVoxels info found in files "<<fCellDefFilePath<<std::endl;
228 std::cerr<<"*************** *****"<<std::endl;
229 exit(EXIT_FAILURE);
230 }
231 std::cout<<"Num of voxels: "<< fVoxels.size()<<" placed in Cell Nucleus."<<std::endl;
232 std::cout<<"=====Choromosome sizes====="<<std::endl;
233 std::cout<<"Chromosome ID\t Number of Bp"<<std::endl;
234 for (auto const& [chrom, nBp] : fChromosomeBpMap) {
235 std::cout<<chrom<< "\t"<<nBp<<std::endl;
236 }
237 std::cout<<"==========================="<<std::endl;
238 }
239
240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
241
242 void ScanDamage::ScanDamageFromPhys()
243 {
244 std::cout<<"===== Start Scanning Damages From Phys =====\n";
245 fEdepSumInNucleus = 0;
246 fphysTables.clear();
247 fphysSlectedTables.clear();
248 fs::path currentP{"phys_output"};
249 fs::file_status s = fs::file_status{};
250 auto isExist = fs::status_known(s) ? fs::exists(s) : fs::exists(currentP);
251 if (isExist) {
252 bool isFoundRootFiles = false;
253 for (const auto entry : fs::directory_iterator(currentP)) {
254 if (entry.path().extension() == ".root") {
255 std::cout <<"ScanDamageFromPhys(): Processing file: "<< entry.path().filename()<< std::endl;
256 AnaPhysRootFile(entry.path());
257 if (!isFoundRootFiles) isFoundRootFiles=true;
258 }
259 }
260 if (!isFoundRootFiles) {
261 std::cout<<"=====>> No root files found in folder \"phys_ouput\"!!! Skip Scanning Damages From Phys =====\n";
262 }
263 if (fphysTables.size() > 0) {
264 SortPhysTableWithSelection();
265 }
266 } else {
267 std::cout<<"=====>> Cannot find folder \"phys_ouput\"!!! Skip Scanning Damages From Phys =====\n";
268 }
269
270 std::cout<<"===== End Scanning Damages From Phys =====\n";
271 }
272
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
274
275 void ScanDamage::ScanDamageFromChem()
276 {
277 std::cout<<"===== Start Scanning Damages From Chem =====\n";
278 std::string fChemOutFolderName = ParametersParser::Instance()->GetChemOutFolderName();
279 if (fChemOutFolderName == "") fChemOutFolderName = "chem_output";
280 fs::path currentP{fChemOutFolderName};
281 fs::file_status s = fs::file_status{};
282 auto isExist = fs::status_known(s) ? fs::exists(s) : fs::exists(currentP);
283 if (isExist) {
284 bool isFoundRootFiles = false;
285 for (const auto entry : fs::directory_iterator(currentP)) {
286 if (entry.path().extension() == ".root") {
287 AnaChemRootFile(entry);
288 if (!isFoundRootFiles) isFoundRootFiles=true;
289 }
290 }
291 if (!isFoundRootFiles) {
292 std::cout<<"=====>> No root files found in folder \""<<fChemOutFolderName<<"\"!!! Skip Scanning Damages From Chem =====\n";
293 fSkipScanningIndirectDamage = true;
294 }
295 if (fchemTables.size() > 0) {
296 SortChemTableWithSelection();
297 }
298 } else {
299 std::cout<<"=====>> Cannot find folder \""<<fChemOutFolderName<<"\"!!! Skip Scanning Damages From Chem =====\n";
300 fSkipScanningIndirectDamage = true;
301 }
302
303 std::cout<<"===== End Scanning Damages From Chem =====\n";
304 }
305
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
307
308 void ScanDamage::AnaPhysRootFile(const std::string fileName)
309 {
310 TFile* f = new TFile(fileName.c_str());
311 if(f->IsZombie() ){
312 // File is corrupted
313 std::cerr<<"*********** Warning *************"<<std::endl;
314 std::cerr<<"The file "<<fileName<<" seems to be corrupted..."<<std::endl;
315 std::cerr<<"We will skip it."<<std::endl;
316 std::cerr<<"**********************************"<<std::endl;
317 return;
318 }
319 AnaPhysRootTree1(f);
320 AnaPhysRootTree2(f);
321 f->Close();
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
325
326 void ScanDamage::AnaChemRootFile(fs::directory_entry entry)
327 {
328 std::string fnameWithoutExtension = entry.path().stem().string();
329 auto [eventNumber, voxelNumber] = GetEventNberAndVoxelNberFromChemRoot(fnameWithoutExtension);
330
331 // Voxel data retrieval
332 VoxelData& voxelData = fVoxels.at(voxelNumber);
333
334 int chromo = voxelData.fChromosome;
335 ullint firstBpNum = voxelData.fFirstBpCopyNum;
336 // *******************
337 // Analyse of the ntuple to detect SB and associate them with a bpNumCorrected
338 // *******************
339
340 // Load the file, the directory and the ntuple
341 TFile f(entry.path().c_str());
342 if(f.IsZombie() ){
343 corruptedFiles++;
344 // File is corrupted
345 std::cerr<<"*********** Warning *************"<<std::endl;
346 std::cerr<<"The file "<<entry.path().string()<<" seems to be corrupted..."<<std::endl;
347 std::cerr<<"We will skip it."<<std::endl;
348 std::cerr<<"Number of corrupted files: "<< corruptedFiles<<std::endl;
349 std::cerr<<"**********************************"<<std::endl;
350 } else {
351 TDirectoryFile *d = dynamic_cast<TDirectoryFile*> (f.Get("ntuple") );
352 TTree* chemTree = (TTree*) d->Get("ntuple_2");
353
354 if( (int) chemTree->GetEntries() >0)
355 {
356 int strand;
357 int copyNumber;
358 double xp;
359 double yp;
360 double zp;
361 double time;
362 int base;
363 chemTree->SetBranchAddress("strand", &strand);
364 chemTree->SetBranchAddress("copyNumber", ©Number);
365 chemTree->SetBranchAddress("xp", &xp);
366 chemTree->SetBranchAddress("yp", &yp);
367 chemTree->SetBranchAddress("zp", &zp);
368 chemTree->SetBranchAddress("time", &time);
369 chemTree->SetBranchAddress("base", &base);
370 unsigned int entryNumber = (int) chemTree->GetEntries();
371 for (unsigned int e=0;e<entryNumber;e++)
372 {
373 chemTree->GetEntry(e);
374 ullint cpNumCorrected = firstBpNum+int(copyNumber);
375 std::vector<ullint> newLine{
376 (ullint)eventNumber,(ullint)strand,cpNumCorrected,
377 (ullint)base,(ullint)(time*1000000)};
378 auto itr = fchemTables.find(chromo);
379 if ( itr == fchemTables.end()) {
380 Table tbforThisChro{newLine};
381 fchemTables.insert({chromo,tbforThisChro});
382 } else (itr->second).push_back(newLine);
383 }
384 }
385 }//
386 }
387
388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
389
390 std::tuple<unsigned int, unsigned int> ScanDamage::GetEventNberAndVoxelNberFromChemRoot(
391 const std::string fileNameWithoutExtension)
392 {
393 unsigned int evnN, volxelN;
394 auto fristPos = fileNameWithoutExtension.find_first_of("_");
395 auto secondPos = fileNameWithoutExtension.substr(fristPos+1).find_first_of("_");
396 auto lastPos = fileNameWithoutExtension.find_last_of("_");
397 evnN = std::stoul(fileNameWithoutExtension.substr(fristPos+1,secondPos));
398 volxelN = std::stoul(fileNameWithoutExtension.substr(lastPos+1));
399 return {evnN, volxelN};
400 }
401
402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
403
404 void ScanDamage::AnaPhysRootTree1(TFile* f)
405 {
406 TDirectoryFile* d = dynamic_cast<TDirectoryFile*> (f->Get("ntuple") );
407 TTree* tPhys = dynamic_cast<TTree*> (d->Get("ntuple_1") );
408
409 if( tPhys->GetEntries() > 0)
410 {
411 int flagParticle;
412 int flagParentID;
413 int flagProcess;
414 double x;
415 double y;
416 double z;
417 double edep;
418 int eventNumber;
419 int volumeName;
420 int copyNumber;
421 int lastMetVoxelCopyNum;
422
423 tPhys->SetBranchAddress("flagParticle", &flagParticle);
424 tPhys->SetBranchAddress("flagParentID", &flagParentID);
425 tPhys->SetBranchAddress("flagProcess", &flagProcess);
426 tPhys->SetBranchAddress("x", &x);
427 tPhys->SetBranchAddress("y", &y);
428 tPhys->SetBranchAddress("z", &z);
429 tPhys->SetBranchAddress("edep", &edep);
430 tPhys->SetBranchAddress("eventNumber", &eventNumber);
431 tPhys->SetBranchAddress("volumeName", &volumeName);
432 tPhys->SetBranchAddress("copyNumber", ©Number);
433 tPhys->SetBranchAddress("lastMetVoxelCopyNum", &lastMetVoxelCopyNum);
434
435 unsigned int entryNumber = tPhys->GetEntries() ;
436
437 // Loop on all the "lines" (ie entry) of the ntuple
438 for(unsigned int e=0; e<entryNumber; e++)
439 {
440 // Set all the variables to the values corresponding to the entry number
441 tPhys->GetEntry(e);
442 // Check if the process is an ionisation
443 // Only ionisation should trigger the removal of a DNA molecule from the chemical step
444 if(flagProcess == 13 // e-_DNAIonisation
445 || flagProcess == 113 // e-_DNAPTBIonisation
446 || flagProcess == 18 // proton_DNAIonisation
447 || flagProcess == 21 // hydrogen_DNAIonisation
448 || flagProcess == 24 // alpha_DNAIonisation
449 || flagProcess == 27 // alpha+_DNAIonisation
450 || flagProcess == 31 // helium_DNAIonisation
451 || flagProcess == 12 // e-_DNAExcitation
452 || flagProcess == 112 // e-_DNAPTBExcitation
453 || flagProcess == 15 // e-_DNAVibExcitation
454 || flagProcess == 17 // proton_DNAExcitation
455 || flagProcess == 20 // hydrogen_DNAExcitation
456 || flagProcess == 23 // alpha_DNAExcitation
457 || flagProcess == 26 // alpha+_DNAExcitation
458 || flagProcess == 30 // helium_DNAExcitation
459 ) {
460 // Check the interaction happened in a dna molecule or its hydration shell
461 if(volumeName == 1 // d1
462 || volumeName == 11 // p1
463 || volumeName == 2 // d2
464 || volumeName == 22 // p2
465 || volumeName == 7 // d1_w
466 || volumeName == 71 // p1_w
467 || volumeName == 8 // d2_w
468 || volumeName == 81 // p2_w
469 )
470 {
471 // *************
472
473 // Retrieve the voxel copy number
474 double voxelCopyNumber = lastMetVoxelCopyNum;
475 if (voxelCopyNumber >= 0 && voxelCopyNumber<fVoxels.size()){
476 // Chromosome, domain and firstNucleotideNum
477 const VoxelData& voxelData = fVoxels.at(size_t(voxelCopyNumber) );
478 int chromo = voxelData.fChromosome;
479 int domain = voxelData.fDomain;
480 ullint firstBpCN = voxelData.fFirstBpCopyNum;
481 ullint cpNumCorrected = firstBpCN+int(copyNumber);
482
483 // Get the event number
484 double eventNum = eventNumber;
485
486 // Determine the strand
487 double strand (-1);
488 if(volumeName==1
489 || volumeName==11
490 || volumeName==7
491 || volumeName==71
492 || volumeName==6 // ade
493 || volumeName==9 // ade
494 || volumeName==4 // gua
495 || volumeName==10) // gua
496 strand = 1;
497 else if(volumeName==2
498 || volumeName==22
499 || volumeName==8
500 || volumeName==81
501 || volumeName==5 // thy
502 || volumeName==12 // thy
503 || volumeName==3 // cyto
504 || volumeName==13) // cyto
505 strand = 2;
506 // Check if the chromo has already been registered
507 std::vector<ullint> newLine{
508 (ullint)eventNum,(ullint)strand,cpNumCorrected,
509 (ullint)volumeName,(ullint)flagProcess,(ullint)(edep*1000000)};
510 // *1000000 in edep is taken back after. This is done
511 // because the number is "unsigned long long int" instead of "double".
512 auto itr = fphysTables.find(chromo);
513 if ( itr == fphysTables.end()) {
514 Table tbforThisChro{newLine};
515 fphysTables.insert({chromo,tbforThisChro});
516 } else (itr->second).push_back(newLine);
517 }
518 }
519 }
520 }
521 }
522 }
523
524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
525
526 void ScanDamage::AnaPhysRootTree2(TFile* f)
527 {
528 TDirectoryFile* d2 = dynamic_cast<TDirectoryFile*> (f->Get("ntuple") );
529 TTree* tPhys2 = dynamic_cast<TTree*> (d2->Get("ntuple_3") );
530
531 if( int(tPhys2->GetEntries() ) > 0)
532 {
533 double edep;
534 int eventNumber;
535
536 tPhys2->SetBranchAddress("edep", &edep);
537 tPhys2->SetBranchAddress("eventNumber", &eventNumber);
538
539 unsigned int entryNumber = int( tPhys2->GetEntries() );
540
541 // Loop on all the "lines" (ie entry) of the ntuple
542 for(unsigned int e=0; e<entryNumber; e++)
543 {
544 tPhys2->GetEntry(e);
545 fEdepSumInNucleus += edep;
546 }
547 }
548 }
549
550 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
551
552 void ScanDamage::SortPhysTableWithSelection()
553 {
554 for(const auto [chrom, physTable] : fphysTables)
555 {
556 // Final table
557 Table physTableWithSelection;
558
559 std::map<ullint,std::map<ullint,std::map<ullint, ullint > > > energyMap;
560
561 // Loop on all the lines of the table
562 for(unsigned int line=0, eline=physTable.size(); line<eline; ++line)
563 {
564 ullint eventNum = physTable[line][0];
565 ullint strand = physTable[line][1];
566 ullint copyNumber = physTable[line][2];
567 ullint volumeFlag = physTable[line][3];
568 ullint processFlagr = physTable[line][4];
569 ullint energy = physTable[line][5];
570
571 // Cumulate the energy value
572 energyMap[eventNum][strand][copyNumber] += energy;
573 }
574
575 int notDuplicatedLine = 0;
576
577 // Loop on all the events
578 std::map<ullint, std::map<ullint, std::map<ullint, ullint> > >::iterator iit = energyMap.begin();
579 std::map<ullint, std::map<ullint, std::map<ullint, ullint> > >::iterator iite = energyMap.end();
580 for(; iit!=iite;++iit)
581 {
582 ullint eventNum = iit->first;
583
584 // Loop on all the strands
585 std::map<ullint, std::map<ullint, ullint> >::iterator itt = iit->second.begin();
586 std::map<ullint, std::map<ullint, ullint> >::iterator itte = iit->second.end();
587 for(; itt!=itte;++itt)
588 {
589 ullint strand = itt->first;
590 // Loop on all the copy numbers
591 std::map<ullint, ullint>::iterator ittt = itt->second.begin();
592 std::map<ullint, ullint>::iterator ittte = itt->second.end();
593 for(; ittt!=ittte;++ittt)
594 {
595 ullint copyNumber = ittt->first;
596 double currentE = double(ittt->second) / 1000000; // eV
597 // Energy condition(s) are set here
598 bool fill = false;
599 // Threshold condition
600 if(currentE < fThresholdEnergy)
601 fill=false;
602 else
603 fill=true;
604
605 if(fill)
606 {
607 // Add a line
608 physTableWithSelection.push_back(std::vector<ullint>());
609 // Fill the line
610 physTableWithSelection[notDuplicatedLine].push_back(eventNum);
611 physTableWithSelection[notDuplicatedLine].push_back(strand);
612 physTableWithSelection[notDuplicatedLine].push_back(copyNumber);
613 physTableWithSelection[notDuplicatedLine].push_back(0);
614 physTableWithSelection[notDuplicatedLine].push_back(0);
615 physTableWithSelection[notDuplicatedLine].push_back((ullint)(currentE*1000000) );
616 physTableWithSelection[notDuplicatedLine].push_back(0);
617 ++notDuplicatedLine;
618 }
619 }
620 }
621 }
622 // *******************************************
623 // Print the "physTableWithSelection" table for the current chromosome
624 // *******************************************
625 std::cout << "### Phys SB for chromosome "<<chrom<<" : " << physTableWithSelection.size() << " ###" << std::endl;
626 if (physTableWithSelection.size()>0) fphysSlectedTables.insert({chrom,physTableWithSelection});
627 }
628 fphysTables.clear();//Free memory
629 }
630
631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
632
633 void ScanDamage::SortChemTableWithSelection()
634 {
635 for (const auto& [chromo,chemTable] : fchemTables)
636 {
637 // Final table
638 Table chemTableWithSelection;
639
640 int notDuplicatedLine = 0;
641
642 for(unsigned int line=0, eline=chemTable.size(); line<eline; ++line)
643 {
644 ullint eventNum = chemTable[line][0];
645 ullint strand = chemTable[line][1];
646 ullint copyNumber = chemTable[line][2];
647 ullint base = chemTable[line][3];
648 ullint time = chemTable[line][4];
649
650 // Random number between 0 and <1
651 if (base==1) {
652 // Add a line
653 chemTableWithSelection.push_back(std::vector<ullint>());
654 // Fill the line
655 chemTableWithSelection[notDuplicatedLine].push_back(eventNum);
656 chemTableWithSelection[notDuplicatedLine].push_back(strand);
657 chemTableWithSelection[notDuplicatedLine].push_back(copyNumber);
658 chemTableWithSelection[notDuplicatedLine].push_back(base);
659 chemTableWithSelection[notDuplicatedLine].push_back(time);
660 chemTableWithSelection[notDuplicatedLine].push_back(0);
661 chemTableWithSelection[notDuplicatedLine].push_back(1);
662 ++notDuplicatedLine;
663 }
664 else {
665 //double r = double(std::rand() ) / RAND_MAX;
666 double r = gRandomGen.Rndm();
667 //std::cout<<r<<std::endl;
668 if(r <= fProbabilityForIndirectSB){
669 // Add a line
670 chemTableWithSelection.push_back(std::vector<ullint>());
671 // Fill the line
672 chemTableWithSelection[notDuplicatedLine].push_back(eventNum);
673 chemTableWithSelection[notDuplicatedLine].push_back(strand);
674 chemTableWithSelection[notDuplicatedLine].push_back(copyNumber);
675 chemTableWithSelection[notDuplicatedLine].push_back(base);
676 chemTableWithSelection[notDuplicatedLine].push_back(time);
677 chemTableWithSelection[notDuplicatedLine].push_back(0);
678 chemTableWithSelection[notDuplicatedLine].push_back(1);
679 ++notDuplicatedLine;
680 }
681 }
682 }
683 if (chemTableWithSelection.size() > 0) fchemSlectedTables.insert({chromo,chemTableWithSelection});
684 }
685 fchemTables.clear();//free memory
686 }
687
688 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
689
690 void ScanDamage::MergeDamageFromPhysChem()
691 {
692 for(int chromo=0; chromo<46; chromo++)
693 {
694 // *****************************
695 // Add one table after the other
696 // *****************************
697
698 // MergedTable to be built
699 Table mergedTable;
700
701 auto itr = fphysSlectedTables.find(chromo);
702 if(itr != fphysSlectedTables.end())
703 {
704 Table physTable= itr->second;
705 for(auto const& vec : physTable) mergedTable.push_back(vec);
706 }
707 itr = fchemSlectedTables.find(chromo);
708 if(itr != fchemSlectedTables.end())
709 {
710 Table chemTable = itr->second;
711 for(auto const& vec : chemTable) mergedTable.push_back(vec);
712 }
713
714 // *******************************************
715 // Sort the merged table to put the event in the correct order
716 // *******************************************
717 Trier tri;
718 Table::iterator it = mergedTable.begin();
719 Table::iterator ite = mergedTable.end();
720 std::sort(it, ite, tri);
721
722 // *******************************************
723 // Delete duplicate SB
724 // *******************************************
725
726 std::map<ullint,std::map<ullint,std::map<ullint,std::map<ullint,std::vector<ullint>>>>> removeDuplicateValueMap;
727
728 // Put all the values of the mergedTable in the map created just above.
729 // Loop on all the lines of the mergedTable
730 for(unsigned int line=0; line<mergedTable.size(); ++line)
731 {
732 ullint eventNum = mergedTable[line][0];
733 ullint strand = mergedTable[line][1];
734 ullint copyNumber = mergedTable[line][2];
735 ullint base = mergedTable[line][3];
736 std::vector<ullint> lineV;
737 // If more elements are presents, add them here
738 int lineSize = mergedTable[line].size();
739 if(lineSize > 4)
740 {
741 for(int i=4; i<lineSize; i++)
742 {
743 lineV.push_back(mergedTable[line][i]);
744 }
745 removeDuplicateValueMap[eventNum][strand][copyNumber][base]=lineV;
746 }
747 else
748 {
749 lineV.push_back(0);
750 removeDuplicateValueMap[eventNum][strand][copyNumber][base]=lineV;
751 }
752 }
753
754 // *******************************************
755 // Create the "mergedTableWithoutDuplicatedSB" table
756 // *******************************************
757
758 // At this point, we have created a map named "removeDuplicateValueMap" that organized all the mergedTable content
759 // AND that does not contain any duplicate because of the override caracteristic of a map.
760 // Indeed, doing map[2] = "hello" followed by map[2] = "bye" will put map[2] value to "bye" because it overrided the first "hello".
761 // This is a cheap way to remove duplicates by overriding them.
762
763 // The next part is dedicated to the creation of the "mergedTableWithoutDuplicatedSB" table from the "removeDuplicateValueMap" map.
764 // Only the event, strand and copynumber and will be put in this final table.
765 Table mergedTableWithoutDuplicatedSB;
766 Table mergedTableWithoutDuplicatedSBandbases;
767 int notDuplicatedLine = 0;
768 int notDuplicatedLine2= 0;
769 // Loop on all the events
770 std::map<ullint,std::map<ullint,std::map<ullint,std::map<ullint,std::vector<ullint>>>>>::iterator
771 iit = removeDuplicateValueMap.begin();
772 std::map<ullint,std::map<ullint,std::map<ullint,std::map<ullint,std::vector<ullint>>>>>::iterator
773 iite = removeDuplicateValueMap.end();
774 for(; iit!=iite;++iit)
775 {
776 ullint eventNum = iit->first;
777 // Loop on all the strands
778 std::map<ullint, std::map<ullint, std::map<ullint, std::vector<ullint > > > >::iterator
779 itt = iit->second.begin();
780 std::map<ullint, std::map<ullint, std::map<ullint, std::vector<ullint > > > >::iterator
781 itte = iit->second.end();
782 for(; itt!=itte;++itt)
783 {
784 ullint strand = itt->first;
785 // Loop on all the copy numbers
786 std::map<ullint, std::map<ullint, std::vector<ullint > > >::iterator ittt = itt->second.begin();
787 std::map<ullint, std::map<ullint, std::vector<ullint > > >::iterator ittte = itt->second.end();
788 for(; ittt!=ittte;++ittt)
789 {
790 ullint copyNumber = ittt->first;
791 // Loop on all the base flag
792 std::map<ullint, std::vector<ullint > >::iterator itttt = ittt->second.begin();
793 std::map<ullint, std::vector<ullint > >::iterator itttte = ittt->second.end();
794 for(; itttt!=itttte;++itttt)
795 {
796 ullint base = itttt->first;
797 // Fill the table
798 if (strand>0 && strand<3)
799 {
800 // Add a line
801 mergedTableWithoutDuplicatedSB.push_back(std::vector<ullint>());
802 // Fill the line
803 mergedTableWithoutDuplicatedSB[notDuplicatedLine].push_back(eventNum);
804 mergedTableWithoutDuplicatedSB[notDuplicatedLine].push_back(strand);
805 mergedTableWithoutDuplicatedSB[notDuplicatedLine].push_back(copyNumber);
806 mergedTableWithoutDuplicatedSB[notDuplicatedLine].push_back(base);
807 mergedTableWithoutDuplicatedSB[notDuplicatedLine].push_back(
808 removeDuplicateValueMap[eventNum][strand][copyNumber][base][0]);
809 mergedTableWithoutDuplicatedSB[notDuplicatedLine].push_back(
810 removeDuplicateValueMap[eventNum][strand][copyNumber][base][1]);
811 mergedTableWithoutDuplicatedSB[notDuplicatedLine].push_back(
812 removeDuplicateValueMap[eventNum][strand][copyNumber][base][2]);
813 ++notDuplicatedLine;
814 if (base==0)
815 {
816 // Add a line
817 mergedTableWithoutDuplicatedSBandbases.push_back(std::vector<ullint>());
818 // Fill the line
819 mergedTableWithoutDuplicatedSBandbases[notDuplicatedLine2].push_back(eventNum);
820 mergedTableWithoutDuplicatedSBandbases[notDuplicatedLine2].push_back(strand);
821 mergedTableWithoutDuplicatedSBandbases[notDuplicatedLine2].push_back(copyNumber);
822 ++notDuplicatedLine2;
823 }
824 }
825 }
826 }
827 }
828 }
829
830 // *******************************************
831 // Print the "mergedTableWithoutDuplicatedSB" table for the current chromosome
832 // *******************************************
833 if (mergedTableWithoutDuplicatedSB.size() > 0) {
834 //PrintTable(fMergeFolder + "/chromo_"+std::to_string(chromo)+".dat",
835 //mergedTableWithoutDuplicatedSB, "eventNum, strand, copyNumber, isbase,
836 //time(ns*1000000), edep, phy:0 chem:1");
837 fMergedTables.insert({chromo,mergedTableWithoutDuplicatedSB});
838 }
839 mergedTable.clear();
840 mergedTableWithoutDuplicatedSBandbases.clear();
841 mergedTableWithoutDuplicatedSB.clear();
842 }
843 //free memory:
844 fphysSlectedTables.clear();
845 fchemSlectedTables.clear();
846 }
847
848 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
849
850 void ScanDamage::ReadCellandVoxelDefFilePaths()
851 {
852 fs::path thisP = fs::current_path();
853 for (const auto entry : fs::directory_iterator(thisP)){
854 if (entry.path().filename() == "imp.info") {
855 std::ifstream file(entry.path().c_str());
856 if(!file.good() ){
857 std::cerr<<"**** Fatal Error *****"<<std::endl;
858 std::cerr<<"ScanDamage::ReadCellandVoxelDefFilePaths(): File corupted: "
859 <<entry.path()<<std::endl;
860 std::cerr<<"*************** *****"<<std::endl;
861 exit(EXIT_FAILURE);
862 }
863
864 std::string line;
865 while(std::getline(file, line) ){
866 std::istringstream iss(line);
867 std::string flag;
868 iss >> flag;
869 if ( flag == "_geovolxelpath") {
870 std::string voxname;
871 iss >> voxname;
872 fVoxelDefFilesList.insert(voxname);
873 }
874 if ( flag == "_geocellpath") {
875 std::string cellpname;
876 iss >> cellpname;
877 fCellDefFilePath = cellpname;
878 }
879 if ( flag == "_numberOfBasepairs") {
880 iss >> fTotalNbBpPlacedInGeo;
881 }
882 if ( flag == "_numberOfHistones") {
883 iss >> fTotalNbHistonePlacedInGeo;
884 }
885 if ( flag == "_nucleusVolume") {
886 iss >> fNucleusVolume;
887 }
888 if ( flag == "_nucleusMassDensity") {
889 iss >> fNucleusMassDensity;
890 }
891 if ( flag == "_nucleusMass") {
892 iss >> fNucleusMass;
893 }
894 }
895 file.close();
896 }
897 }
898 }
899
900 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....