| 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 ScanDamage.cc
28 /// \brief Implementation of the ScanDamage cl
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........oo
47
48 ScanDamage::ScanDamage(): fSkipScanningIndirec
49 {
50 fThresholdEnergy = 17.5; //eV
51 fEdepSumInNucleus = 0;
52 fProbabilityForIndirectSB = 0.40; // 40%
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 std::map<unsigned int,std::map<unsigned int,st
58 ReadCellandVoxelDefFilePaths();
59 fMergedTables.clear();
60 fDamage.clear();
61 RetrieveVoxelBp();
62 FillVoxelData();
63 ScanDamageFromPhys();
64 if (!fSkipScanningIndirectDamage) ScanDama
65 MergeDamageFromPhysChem();
66 for (const auto& [pChrom,table] : fMergedT
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<Dama
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::fBac
87 else
88 pType=Damage::DamageType::fBas
89 if(origin==0)
90 pOrigin=Damage::DamageCause::f
91 else
92 pOrigin=Damage::DamageCause::f
93 if (perChromoDamage.find(evt) == p
94 std::vector<Damage> dm{Damage(
95 pO
96 perChromoDamage.insert({evt,dm
97 } else perChromoDamage[evt].push_b
98 pO
99 }
100 fDamage.insert({pChrom,perChromoDamage
101 }
102 return fDamage;
103 }
104
105 //....oooOO0OOooo........oooOO0OOooo........oo
106
107 void ScanDamage::RetrieveVoxelBp()
108 {
109 fBpPerVoxel.clear();
110
111 for (const auto & entry : fVoxelDefFilesLi
112 std::ifstream file(entry);
113 if(!file.good() )
114 {
115 std::cerr<<"**** Fatal Error *****
116 std::cerr<<"ScanDamage::RetrieveVo
117 std::cerr<<"*************** *****"
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 "voxelBa
140 if(flag=="_Number" && charac=="vox
141 {
142 int numOfBp;
143 iss >> numOfBp;
144 if(!foundName)
145 {
146 std::cerr<<"*** Fatal Erro
147 std::cerr<<"ScanDamage::Re
148 <<"of the voxel... This is
149 std::cerr<<"******"<<std::
150 exit(EXIT_FAILURE);
151 }
152 else
153 {
154 fBpPerVoxel[voxelName] = n
155 std::cout<<voxelName<<" ha
156 foundNumOfBp = true;
157 }
158 }
159 }
160 file.close();
161 }
162
163 if (fBpPerVoxel.size() == 0) {
164 std::cerr<<"**** Fatal Error *****"<<s
165 std::cerr<<"ScanDamage::RetrieveVoxelB
166 <<" sure that file imp.info exists in
167 std::cerr<<"*************** *****"<<st
168 exit(EXIT_FAILURE);
169 }
170 }
171
172 //....oooOO0OOooo........oooOO0OOooo........oo
173
174 void ScanDamage::FillVoxelData()
175 {
176 std::ifstream file(fCellDefFilePath);
177 if(!file.good() )
178 {
179 std::cerr<<"**** Fatal Error *****"<<s
180 std::cerr<<"FillVoxelData: No file nam
181 std::cerr<<"*************** *****"<<st
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 place
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
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(chrom
212 int numBpinthisvoxel = int(fBpPerV
213 bpCount += numBpinthisvoxel;
214 if (fChromosomeBpMap.find(chromo)
215 fChromosomeBpMap.insert({chrom
216 } else {
217 fChromosomeBpMap[chromo] += nu
218 }
219 voxelCount++;
220 }
221 }
222 file.close();
223
224
225 if (fVoxels.size() == 0) {
226 std::cerr<<"**** Fatal Error *****"<<s
227 std::cerr<<"ScanDamage::FillVoxelData:
228 std::cerr<<"*************** *****"<<st
229 exit(EXIT_FAILURE);
230 }
231 std::cout<<"Num of voxels: "<< fVoxels.siz
232 std::cout<<"=====Choromosome sizes====="<<
233 std::cout<<"Chromosome ID\t Number of Bp"<
234 for (auto const& [chrom, nBp] : fChromosom
235 std::cout<<chrom<< "\t"<<nBp<<std::en
236 }
237 std::cout<<"==========================="<<
238 }
239
240 //....oooOO0OOooo........oooOO0OOooo........oo
241
242 void ScanDamage::ScanDamageFromPhys()
243 {
244 std::cout<<"===== Start Scanning Damages F
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::e
251 if (isExist) {
252 bool isFoundRootFiles = false;
253 for (const auto entry : fs::directory_
254 if (entry.path().extension() == ".
255 std::cout <<"ScanDamageFromPhy
256 AnaPhysRootFile(entry.path());
257 if (!isFoundRootFiles) isFound
258 }
259 }
260 if (!isFoundRootFiles) {
261 std::cout<<"=====>> No root files
262 }
263 if (fphysTables.size() > 0) {
264 SortPhysTableWithSelection();
265 }
266 } else {
267 std::cout<<"=====>> Cannot find folder
268 }
269
270 std::cout<<"===== End Scanning Damages Fro
271 }
272
273 //....oooOO0OOooo........oooOO0OOooo........oo
274
275 void ScanDamage::ScanDamageFromChem()
276 {
277 std::cout<<"===== Start Scanning Damages F
278 std::string fChemOutFolderName = Parameter
279 if (fChemOutFolderName == "") fChemOutFold
280 fs::path currentP{fChemOutFolderName};
281 fs::file_status s = fs::file_status{};
282 auto isExist = fs::status_known(s) ? fs::e
283 if (isExist) {
284 bool isFoundRootFiles = false;
285 for (const auto entry : fs::directory_
286 if (entry.path().extension() == ".
287 AnaChemRootFile(entry);
288 if (!isFoundRootFiles) isFound
289 }
290 }
291 if (!isFoundRootFiles) {
292 std::cout<<"=====>> No root files
293 fSkipScanningIndirectDamage = true
294 }
295 if (fchemTables.size() > 0) {
296 SortChemTableWithSelection();
297 }
298 } else {
299 std::cout<<"=====>> Cannot find folder
300 fSkipScanningIndirectDamage = true;
301 }
302
303 std::cout<<"===== End Scanning Damages Fro
304 }
305
306 //....oooOO0OOooo........oooOO0OOooo........oo
307
308 void ScanDamage::AnaPhysRootFile(const std::st
309 {
310 TFile* f = new TFile(fileName.c_str());
311 if(f->IsZombie() ){
312 // File is corrupted
313 std::cerr<<"*********** Warning ******
314 std::cerr<<"The file "<<fileName<<" se
315 std::cerr<<"We will skip it."<<std::en
316 std::cerr<<"**************************
317 return;
318 }
319 AnaPhysRootTree1(f);
320 AnaPhysRootTree2(f);
321 f->Close();
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oo
325
326 void ScanDamage::AnaChemRootFile(fs::directory
327 {
328 std::string fnameWithoutExtension = entry.
329 auto [eventNumber, voxelNumber] = GetEvent
330
331 // Voxel data retrieval
332 VoxelData& voxelData = fVoxels.at(voxelNum
333
334 int chromo = voxelData.fChromosome;
335 ullint firstBpNum = voxelData.fFirstBpCopy
336 // *******************
337 // Analyse of the ntuple to detect SB and
338 // *******************
339
340 // Load the file, the directory and the nt
341 TFile f(entry.path().c_str());
342 if(f.IsZombie() ){
343 corruptedFiles++;
344 // File is corrupted
345 std::cerr<<"*********** Warning ******
346 std::cerr<<"The file "<<entry.path().s
347 std::cerr<<"We will skip it."<<std::en
348 std::cerr<<"Number of corrupted files:
349 std::cerr<<"**************************
350 } else {
351 TDirectoryFile *d = dynamic_cast<TDire
352 TTree* chemTree = (TTree*) d->Get("ntu
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
364 chemTree->SetBranchAddress("copyNu
365 chemTree->SetBranchAddress("xp", &
366 chemTree->SetBranchAddress("yp", &
367 chemTree->SetBranchAddress("zp", &
368 chemTree->SetBranchAddress("time",
369 chemTree->SetBranchAddress("base",
370 unsigned int entryNumber = (int) c
371 for (unsigned int e=0;e<entryNumbe
372 {
373 chemTree->GetEntry(e);
374 ullint cpNumCorrected = firstB
375 std::vector<ullint> newLine{
376 (ullint)eventNumber,(ullin
377 (ullint)base,(ullint)(time
378 auto itr = fchemTables.find(ch
379 if ( itr == fchemTables.end())
380 Table tbforThisChro{newLin
381 fchemTables.insert({chromo
382 } else (itr->second).push_back
383 }
384 }
385 }//
386 }
387
388 //....oooOO0OOooo........oooOO0OOooo........oo
389
390 std::tuple<unsigned int, unsigned int> ScanDam
391 const std::string fileNameWithoutExtension
392 {
393 unsigned int evnN, volxelN;
394 auto fristPos = fileNameWithoutExtension.f
395 auto secondPos = fileNameWithoutExtension.
396 auto lastPos = fileNameWithoutExtension.fi
397 evnN = std::stoul(fileNameWithoutExtension
398 volxelN = std::stoul(fileNameWithoutExtens
399 return {evnN, volxelN};
400 }
401
402 //....oooOO0OOooo........oooOO0OOooo........oo
403
404 void ScanDamage::AnaPhysRootTree1(TFile* f)
405 {
406 TDirectoryFile* d = dynamic_cast<TDirector
407 TTree* tPhys = dynamic_cast<TTree*> (d->Ge
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"
424 tPhys->SetBranchAddress("flagParentID"
425 tPhys->SetBranchAddress("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",
431 tPhys->SetBranchAddress("volumeName",
432 tPhys->SetBranchAddress("copyNumber",
433 tPhys->SetBranchAddress("lastMetVoxelC
434
435 unsigned int entryNumber = tPhys->Get
436
437 // Loop on all the "lines" (ie entry)
438 for(unsigned int e=0; e<entryNumber; e
439 {
440 // Set all the variables to the va
441 tPhys->GetEntry(e);
442 // Check if the process is an ioni
443 // Only ionisation should trigger
444 if(flagProcess == 13 // e-_DNAIoni
445 || flagProcess == 113 // e
446 || flagProcess == 18 // pr
447 || flagProcess == 21 // hy
448 || flagProcess == 24 // al
449 || flagProcess == 27 // al
450 || flagProcess == 31 // he
451 || flagProcess == 12 // e-
452 || flagProcess == 112 // e
453 || flagProcess == 15 // e-
454 || flagProcess == 17 // pr
455 || flagProcess == 20 // hy
456 || flagProcess == 23 // al
457 || flagProcess == 26 // al
458 || flagProcess == 30 // he
459 ) {
460 // Check the interaction happe
461 if(volumeName == 1 // d1
462 || volumeName == 11 //
463 || volumeName == 2 //
464 || volumeName == 22 //
465 || volumeName == 7 //
466 || volumeName == 71 //
467 || volumeName == 8 //
468 || volumeName == 81 //
469 )
470 {
471 // *************
472
473 // Retrieve the voxel copy
474 double voxelCopyNumber =
475 if (voxelCopyNumber >= 0 &
476 // Chromosome, domain
477 const VoxelData& voxel
478 int chromo = voxelData
479 int domain = voxelData
480 ullint firstBpCN = vox
481 ullint cpNumCorrected
482
483 // Get the event numbe
484 double eventNum = even
485
486 // Determine the stran
487 double strand (-1);
488 if(volumeName==1
489 || volumeName=
490 || volumeName=
491 || volumeName=
492 || volumeName=
493 || volumeName=
494 || volumeName=
495 || volumeName=
496 strand = 1;
497 else if(volumeName==2
498 || volumeName=
499 || volumeName=
500 || volumeName=
501 || volumeName=
502 || volumeName=
503 || volumeName=
504 || volumeName=
505 strand = 2;
506 // Check if the chromo
507 std::vector<ullint> ne
508 (ullint)eventNum,(
509 (ullint)volumeName
510 // *1000000 in ede
511 // because the num
512 auto itr = fphysTables
513 if ( itr == fphysTable
514 Table tbforThisChr
515 fphysTables.insert
516 } else (itr->second).p
517 }
518 }
519 }
520 }
521 }
522 }
523
524 //....oooOO0OOooo........oooOO0OOooo........oo
525
526 void ScanDamage::AnaPhysRootTree2(TFile* f)
527 {
528 TDirectoryFile* d2 = dynamic_cast<TDirecto
529 TTree* tPhys2 = dynamic_cast<TTree*> (d2->
530
531 if( int(tPhys2->GetEntries() ) > 0)
532 {
533 double edep;
534 int eventNumber;
535
536 tPhys2->SetBranchAddress("edep", &edep
537 tPhys2->SetBranchAddress("eventNumber"
538
539 unsigned int entryNumber = int( tPhys2
540
541 // Loop on all the "lines" (ie entry)
542 for(unsigned int e=0; e<entryNumber; e
543 {
544 tPhys2->GetEntry(e);
545 fEdepSumInNucleus += edep;
546 }
547 }
548 }
549
550 //....oooOO0OOooo........oooOO0OOooo........oo
551
552 void ScanDamage::SortPhysTableWithSelection()
553 {
554 for(const auto [chrom, physTable] : fphysT
555 {
556 // Final table
557 Table physTableWithSelection;
558
559 std::map<ullint,std::map<ullint,std::m
560
561 // Loop on all the lines of the table
562 for(unsigned int line=0, eline=physTab
563 {
564 ullint eventNum = physTable[line][
565 ullint strand = physTable[line][1]
566 ullint copyNumber = physTable[line
567 ullint volumeFlag = physTable[line
568 ullint processFlagr = physTable[li
569 ullint energy = physTable[line][5]
570
571 // Cumulate the energy value
572 energyMap[eventNum][strand][copyNu
573 }
574
575 int notDuplicatedLine = 0;
576
577 // Loop on all the events
578 std::map<ullint, std::map<ullint, std:
579 std::map<ullint, std::map<ullint, std:
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,
586 std::map<ullint, std::map<ullint,
587 for(; itt!=itte;++itt)
588 {
589 ullint strand = itt->first;
590 // Loop on all the copy number
591 std::map<ullint, ullint>::iter
592 std::map<ullint, ullint>::iter
593 for(; ittt!=ittte;++ittt)
594 {
595 ullint copyNumber = ittt->
596 double currentE = double(i
597 // Energy condition(s) are
598 bool fill = false;
599 // Threshold condition
600 if(currentE < fThresholdE
601 fill=false;
602 else
603 fill=true;
604
605 if(fill)
606 {
607 // Add a line
608 physTableWithSelection
609 // Fill the line
610 physTableWithSelection
611 physTableWithSelection
612 physTableWithSelection
613 physTableWithSelection
614 physTableWithSelection
615 physTableWithSelection
616 physTableWithSelection
617 ++notDuplicatedLine;
618 }
619 }
620 }
621 }
622 // ***********************************
623 // Print the "physTableWithSelection"
624 // ***********************************
625 std::cout << "### Phys SB for chromoso
626 if (physTableWithSelection.size()>0) f
627 }
628 fphysTables.clear();//Free memory
629 }
630
631 //....oooOO0OOooo........oooOO0OOooo........oo
632
633 void ScanDamage::SortChemTableWithSelection()
634 {
635 for (const auto& [chromo,chemTable] : fche
636 {
637 // Final table
638 Table chemTableWithSelection;
639
640 int notDuplicatedLine = 0;
641
642 for(unsigned int line=0, eline=chemTab
643 {
644 ullint eventNum = chemTable[line][
645 ullint strand = chemTable[line][1]
646 ullint copyNumber = chemTable[line
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_ba
654 // Fill the line
655 chemTableWithSelection[notDupl
656 chemTableWithSelection[notDupl
657 chemTableWithSelection[notDupl
658 chemTableWithSelection[notDupl
659 chemTableWithSelection[notDupl
660 chemTableWithSelection[notDupl
661 chemTableWithSelection[notDupl
662 ++notDuplicatedLine;
663 }
664 else {
665 //double r = double(std::rand() ) /
666 double r = gRandomGen.Rndm();
667 //std::cout<<r<<std::endl;
668 if(r <= fProbabilityForIndirectSB){
669 // Add a line
670 chemTableWithSelection.pus
671 // Fill the line
672 chemTableWithSelection[not
673 chemTableWithSelection[not
674 chemTableWithSelection[not
675 chemTableWithSelection[not
676 chemTableWithSelection[not
677 chemTableWithSelection[not
678 chemTableWithSelection[not
679 ++notDuplicatedLine;
680 }
681 }
682 }
683 if (chemTableWithSelection.size() > 0)
684 }
685 fchemTables.clear();//free memory
686 }
687
688 //....oooOO0OOooo........oooOO0OOooo........oo
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(chr
702 if(itr != fphysSlectedTables.end())
703 {
704 Table physTable= itr->second;
705 for(auto const& vec : physTable) m
706 }
707 itr = fchemSlectedTables.find(chromo);
708 if(itr != fchemSlectedTables.end())
709 {
710 Table chemTable = itr->second;
711 for(auto const& vec : chemTable) m
712 }
713
714 // ***********************************
715 // Sort the merged table to put the ev
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::m
727
728 // Put all the values of the mergedTab
729 // Loop on all the lines of the merged
730 for(unsigned int line=0; line<mergedTa
731 {
732 ullint eventNum = mergedTable[line
733 ullint strand = mergedTable[line][
734 ullint copyNumber = mergedTable[li
735 ullint base = mergedTable[line][3]
736 std::vector<ullint> lineV;
737 // If more elements are presents,
738 int lineSize = mergedTable[line].s
739 if(lineSize > 4)
740 {
741 for(int i=4; i<lineSize; i++)
742 {
743 lineV.push_back(mergedTabl
744 }
745 removeDuplicateValueMap[eventN
746 }
747 else
748 {
749 lineV.push_back(0);
750 removeDuplicateValueMap[eventN
751 }
752 }
753
754 // ***********************************
755 // Create the "mergedTableWithoutDupli
756 // ***********************************
757
758 // At this point, we have created a ma
759 // AND that does not contain any dupli
760 // Indeed, doing map[2] = "hello" foll
761 // This is a cheap way to remove dupli
762
763 // The next part is dedicated to the c
764 // Only the event, strand and copynumb
765 Table mergedTableWithoutDuplicatedSB;
766 Table mergedTableWithoutDuplicatedSBan
767 int notDuplicatedLine = 0;
768 int notDuplicatedLine2= 0;
769 // Loop on all the events
770 std::map<ullint,std::map<ullint,std::m
771 iit = removeDuplicateValueMap.begi
772 std::map<ullint,std::map<ullint,std::m
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,
779 itt = iit->second.begin();
780 std::map<ullint, std::map<ullint,
781 itte = iit->second.end();
782 for(; itt!=itte;++itt)
783 {
784 ullint strand = itt->first;
785 // Loop on all the copy number
786 std::map<ullint, std::map<ulli
787 std::map<ullint, std::map<ulli
788 for(; ittt!=ittte;++ittt)
789 {
790 ullint copyNumber = ittt->
791 // Loop on all the base flag
792 std::map<ullint, std::vect
793 std::map<ullint, std::vect
794 for(; itttt!=itttte;++ittt
795 {
796 ullint base = itttt->f
797 // Fill the table
798 if (strand>0 && strand
799 {
800 // Add a line
801 mergedTableWithout
802 // Fill the line
803 mergedTableWithout
804 mergedTableWithout
805 mergedTableWithout
806 mergedTableWithout
807 mergedTableWithout
808 removeDuplicat
809 mergedTableWithout
810 removeDuplicat
811 mergedTableWithout
812 removeDuplicat
813 ++notDuplicatedLin
814 if (base==0)
815 {
816 // Add a line
817 mergedTableWit
818 // Fill the li
819 mergedTableWit
820 mergedTableWit
821 mergedTableWit
822 ++notDuplicate
823 }
824 }
825 }
826 }
827 }
828 }
829
830 // ***********************************
831 // Print the "mergedTableWithoutDuplic
832 // ***********************************
833 if (mergedTableWithoutDuplicatedSB.siz
834 //PrintTable(fMergeFolder + "/chro
835 //mergedTableWithoutDuplicatedSB,
836 //time(ns*1000000), edep, phy:0 ch
837 fMergedTables.insert({chromo,merge
838 }
839 mergedTable.clear();
840 mergedTableWithoutDuplicatedSBandbases
841 mergedTableWithoutDuplicatedSB.clear()
842 }
843 //free memory:
844 fphysSlectedTables.clear();
845 fchemSlectedTables.clear();
846 }
847
848 //....oooOO0OOooo........oooOO0OOooo........oo
849
850 void ScanDamage::ReadCellandVoxelDefFilePaths(
851 {
852 fs::path thisP = fs::current_path();
853 for (const auto entry : fs::directory_iter
854 if (entry.path().filename() == "imp.in
855 std::ifstream file(entry.path().c_
856 if(!file.good() ){
857 std::cerr<<"**** Fatal Error *
858 std::cerr<<"ScanDamage::ReadCe
859 <<entry.path()<<std::endl;
860 std::cerr<<"*************** **
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(
873 }
874 if ( flag == "_geocellpath") {
875 std::string cellpname;
876 iss >> cellpname;
877 fCellDefFilePath = cellpna
878 }
879 if ( flag == "_numberOfBasepai
880 iss >> fTotalNbBpPlacedInG
881 }
882 if ( flag == "_numberOfHistone
883 iss >> fTotalNbHistonePlac
884 }
885 if ( flag == "_nucleusVolume")
886 iss >> fNucleusVolume;
887 }
888 if ( flag == "_nucleusMassDens
889 iss >> fNucleusMassDensity
890 }
891 if ( flag == "_nucleusMass") {
892 iss >> fNucleusMass;
893 }
894 }
895 file.close();
896 }
897 }
898 }
899
900 //....oooOO0OOooo........oooOO0OOooo........oo