| Geant4 Cross Reference |
1 // ------------------------------------------- 1 // -------------------------------------------------------------------
2 // $Id: scandamages.C 2 // $Id: scandamages.C
3 // ------------------------------------------- 3 // -------------------------------------------------------------------
4 // 13/10/2023 - Le Tuan Anh renamed and improv 4 // 13/10/2023 - Le Tuan Anh renamed and improved
5 // ******************************************* 5 // *********************************************************************
6 // To execute this macro under ROOT after your 6 // To execute this macro under ROOT after your simulation ended,
7 // 1 - launch ROOT v6.xx (usually type 'root 7 // 1 - launch ROOT v6.xx (usually type 'root' at your machine's prompt)
8 // 2 - type '.X analysis.C' at the ROOT sess 8 // 2 - type '.X analysis.C' at the ROOT session prompt
9 // ******************************************* 9 // *********************************************************************
10 10
11 //....oooOO0OOooo........oooOO0OOooo........oo 11 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
12 12
13 struct Damage; 13 struct Damage;
14 struct Cluster; 14 struct Cluster;
15 void CreateSDDfile(std::map<int,std::vector<Da 15 void CreateSDDfile(std::map<int,std::vector<Damage>>,long int);
16 void ExtractDirectAndIndirectDamages(std::map< 16 void ExtractDirectAndIndirectDamages(std::map<int,std::vector<Damage>>&,long int&,double,double);
17 void ClassifyDamages(std::map<int,std::vector< 17 void ClassifyDamages(std::map<int,std::vector<Damage>>&,int,
18 std::map<int,vector<Cluste 18 std::map<int,vector<Cluster>>&,int verbose =0);
19 19
20 //....oooOO0OOooo........oooOO0OOooo........oo 20 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
21 // play with scandamages() by changing firect/ 21 // play with scandamages() by changing firect/indirect criterion selection or Clusterdistance
22 22
23 void scandamages() 23 void scandamages()
24 { 24 {
25 ifstream fs("output_t0.root"); 25 ifstream fs("output_t0.root");
26 if ( fs.good() ) { 26 if ( fs.good() ) {
27 system ("rm -rf output.root"); 27 system ("rm -rf output.root");
28 system ("hadd output.root output_*.root" 28 system ("hadd output.root output_*.root");
29 } 29 }
30 30
31 double fEnergyThresholdForDirectSelection 31 double fEnergyThresholdForDirectSelection = 17.5; // eV
32 double fProbabilityForIndirectionSelection 32 double fProbabilityForIndirectionSelection = 0.40; // 40 %
33 int fClusterdistance = 10; // bp, minium d 33 int fClusterdistance = 10; // bp, minium distance that separate two clusters
34 34
35 std::map<int,std::vector<Damage>> fDamages 35 std::map<int,std::vector<Damage>> fDamagesEventMap;
36 std::map<int,vector<Cluster>> fClustersMap 36 std::map<int,vector<Cluster>> fClustersMap;
37 long int totalDamages = 0; 37 long int totalDamages = 0;
38 38
39 ExtractDirectAndIndirectDamages(fDamagesEv 39 ExtractDirectAndIndirectDamages(fDamagesEventMap,totalDamages,
40 fEnergyThr 40 fEnergyThresholdForDirectSelection,
41 fProbabili 41 fProbabilityForIndirectionSelection);
42 42
43 ClassifyDamages(fDamagesEventMap,fClusterd 43 ClassifyDamages(fDamagesEventMap,fClusterdistance,fClustersMap);
44 // print results in each event for checkin 44 // print results in each event for checking:
45 //ClassifyDamages(fDamagesEventMap,fCluste 45 //ClassifyDamages(fDamagesEventMap,fClusterdistance,fClustersMap,1);
46 46
47 CreateSDDfile(fDamagesEventMap,totalDamage 47 CreateSDDfile(fDamagesEventMap,totalDamages);
48 } 48 }
49 49
50 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 51
52 struct Damage 52 struct Damage
53 { 53 {
54 Damage(int a, int b, int c=-1): strand(a), 54 Damage(int a, int b, int c=-1): strand(a),copyNb(b), cause(c){}
55 int strand; 55 int strand;
56 int copyNb; 56 int copyNb;
57 int cause; // Unknown = -1, Direct = 0, In 57 int cause; // Unknown = -1, Direct = 0, Indirect = 1
58 bool operator<(const Damage& r) const {ret 58 bool operator<(const Damage& r) const {return copyNb<r.copyNb;} // for sorting
59 }; 59 };
60 60
61 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 62
63 struct Cluster // to score SSB or DSB 63 struct Cluster // to score SSB or DSB
64 { 64 {
65 std::vector<Damage> fDamagesVector; 65 std::vector<Damage> fDamagesVector;
66 bool HasLeftStrand=false, HasRightStrand=f 66 bool HasLeftStrand=false, HasRightStrand=false;
67 int firstCopyNb; 67 int firstCopyNb;
68 }; 68 };
69 69
70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 71
72 void 72 void
73 ClassifyDamages(std::map<int,std::vector<Damag 73 ClassifyDamages(std::map<int,std::vector<Damage>> &fDamages,int fClusterditance,
74 std::map<int,vector<Cluster>>& 74 std::map<int,vector<Cluster>>& fClustersMap,int verbose)
75 { 75 {
76 long int dirSB = 0, indirSB = 0; 76 long int dirSB = 0, indirSB = 0;
77 long int nDSB = 0, nsDSB = 0, ncDSB=0, nSS 77 long int nDSB = 0, nsDSB = 0, ncDSB=0, nSSB = 0;
78 std::cout<<"------------------------------ 78 std::cout<<"-----------------------------------------------"<<std::endl;
79 if (verbose>0)std::cout<<" ------> Print f 79 if (verbose>0)std::cout<<" ------> Print few events for checking <------"<<std::endl;
80 if (verbose>0)std::cout<<"\t\tCopyNb\tStra 80 if (verbose>0)std::cout<<"\t\tCopyNb\tStrand\tCause"<<std::endl;
81 for (auto &[eventID, dmgVector] : fDamages 81 for (auto &[eventID, dmgVector] : fDamages) {
82 std::sort(dmgVector.begin(),dmgVector. 82 std::sort(dmgVector.begin(),dmgVector.end());
83 if (verbose>0) std::cout<<"Event #"<<e 83 if (verbose>0) std::cout<<"Event #"<<eventID<<": "<<std::endl;
84 unsigned long fSBcountThisEvent = 0; 84 unsigned long fSBcountThisEvent = 0;
85 long copyNBofPrevDamage = 0; 85 long copyNBofPrevDamage = 0;
86 long firstCopyNb; 86 long firstCopyNb;
87 std::map<int,Cluster> clustersInthisEv 87 std::map<int,Cluster> clustersInthisEvent;
88 for (auto dmg : dmgVector) { 88 for (auto dmg : dmgVector) {
89 if (dmg.cause == 0) dirSB++; 89 if (dmg.cause == 0) dirSB++;
90 if (dmg.cause == 1) indirSB++; 90 if (dmg.cause == 1) indirSB++;
91 if (verbose>0) std::cout<<"\t\t"<< 91 if (verbose>0) std::cout<<"\t\t"<<dmg.copyNb<<"\t"
92 <<dmg.stra 92 <<dmg.strand<<"\t"<<dmg.cause<<"\t"<<std::endl;
93 if (clustersInthisEvent.size()==0 93 if (clustersInthisEvent.size()==0 ||
94 ((dmg.copyNb - copyNBofPrevDam 94 ((dmg.copyNb - copyNBofPrevDamage) >= fClusterditance)) {
95 Cluster newCluster; 95 Cluster newCluster;
96 newCluster.fDamagesVector.push 96 newCluster.fDamagesVector.push_back(dmg);
97 copyNBofPrevDamage = dmg.copyN 97 copyNBofPrevDamage = dmg.copyNb;
98 firstCopyNb = dmg.copyNb; 98 firstCopyNb = dmg.copyNb;
99 if (dmg.strand == 0) newCluste 99 if (dmg.strand == 0) newCluster.HasLeftStrand = true;
100 if (dmg.strand == 1) newCluste 100 if (dmg.strand == 1) newCluster.HasRightStrand = true;
101 clustersInthisEvent.insert({fi 101 clustersInthisEvent.insert({firstCopyNb,newCluster});
102 } else { 102 } else {
103 clustersInthisEvent[firstCopyN 103 clustersInthisEvent[firstCopyNb].fDamagesVector.push_back(dmg);
104 copyNBofPrevDamage = dmg.copyN 104 copyNBofPrevDamage = dmg.copyNb;
105 if (dmg.strand == 0 ) clusters 105 if (dmg.strand == 0 ) clustersInthisEvent[firstCopyNb].HasLeftStrand = true;
106 if (dmg.strand == 1 ) clusters 106 if (dmg.strand == 1 ) clustersInthisEvent[firstCopyNb].HasRightStrand = true;
107 } 107 }
108 } 108 }
109 109
110 for (auto [firstCpNb,acluster] : clust 110 for (auto [firstCpNb,acluster] : clustersInthisEvent) {
111 if (fClustersMap.find(eventID) == 111 if (fClustersMap.find(eventID) == fClustersMap.end()) {
112 std::vector<Cluster> sclv{aclu 112 std::vector<Cluster> sclv{acluster};
113 fClustersMap.insert({eventID,s 113 fClustersMap.insert({eventID,sclv});
114 } else { 114 } else {
115 fClustersMap[eventID].push_bac 115 fClustersMap[eventID].push_back(acluster);
116 } 116 }
117 } 117 }
118 clustersInthisEvent.clear(); 118 clustersInthisEvent.clear();
119 } 119 }
120 if (verbose>0)std::cout<<" \n------> Check 120 if (verbose>0)std::cout<<" \n------> Checking SSB and DSB:"<<std::endl;
121 if (verbose>0)std::cout<<"\t\t#SSB\t#DSB\t 121 if (verbose>0)std::cout<<"\t\t#SSB\t#DSB\t#cDSB\t#sDSB"<<std::endl;
122 122
123 for (auto [evt,clustervector] : fClustersM 123 for (auto [evt,clustervector] : fClustersMap) {
124 long int nDSBinThisEvt = 0, nSSBinThis 124 long int nDSBinThisEvt = 0, nSSBinThisEvt = 0;
125 long int nsDSBinThisEvt = 0, ncDSBinTh 125 long int nsDSBinThisEvt = 0, ncDSBinThisEvt = 0;
126 if (verbose>0) std::cout<<"Event #"<<e 126 if (verbose>0) std::cout<<"Event #"<<evt<<": "<<std::endl; // maximum 5 event for cheking
127 for (auto sb : clustervector) { 127 for (auto sb : clustervector) {
128 if (sb.HasLeftStrand && sb.HasRigh 128 if (sb.HasLeftStrand && sb.HasRightStrand) {
129 nDSBinThisEvt++; 129 nDSBinThisEvt++;
130 if (sb.fDamagesVector.size()>2 130 if (sb.fDamagesVector.size()>2) ncDSBinThisEvt++;
131 else nsDSBinThisEvt++; 131 else nsDSBinThisEvt++;
132 } 132 }
133 else nSSBinThisEvt++; 133 else nSSBinThisEvt++;
134 } 134 }
135 nSSB += nSSBinThisEvt; 135 nSSB += nSSBinThisEvt;
136 nDSB += nDSBinThisEvt; 136 nDSB += nDSBinThisEvt;
137 ncDSB += ncDSBinThisEvt; 137 ncDSB += ncDSBinThisEvt;
138 nsDSB += nsDSBinThisEvt; 138 nsDSB += nsDSBinThisEvt;
139 if (verbose>0) std::cout<<"\t\t"<<nSSB 139 if (verbose>0) std::cout<<"\t\t"<<nSSBinThisEvt<<"\t"<<nDSBinThisEvt
140 <<"\t"<<ncDSBi 140 <<"\t"<<ncDSBinThisEvt<<"\t"<<nsDSBinThisEvt<<std::endl;
141 } 141 }
142 std::cout<<"------------------------------ 142 std::cout<<"-----------------------------------------------"<<std::endl;
143 std::cout<<"\tSummary of results:" 143 std::cout<<"\tSummary of results:"
144 <<"\n #Total SBs = "<<(dirSB+i 144 <<"\n #Total SBs = "<<(dirSB+indirSB)
145 <<"\n #dirSB = "<<dirSB << 145 <<"\n #dirSB = "<<dirSB <<" \t;\t #indirSB = "<<indirSB
146 <<"\n #SSB = "<<nSSB << 146 <<"\n #SSB = "<<nSSB <<" \t;\t #DSB = "<<nDSB
147 <<"\n #sDSB = "<<nsDSB << 147 <<"\n #sDSB = "<<nsDSB <<" \t;\t #cDSB = "<<ncDSB<<std::endl;
148 } 148 }
149 149
150 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 151
152 void ExtractDirectAndIndirectDamages(std::map< 152 void ExtractDirectAndIndirectDamages(std::map<int,std::vector<Damage>> &fDamagesEventMap,
153 long int &totalDamages,double Ethresh,double C 153 long int &totalDamages,double Ethresh,double ChemPro)
154 { 154 {
155 TRandom rdomN; 155 TRandom rdomN;
156 156
157 vector<Molecule> fMolecules = molecule(); 157 vector<Molecule> fMolecules = molecule();
158 using Table = map<int,array<map<int,double 158 using Table = map<int,array<map<int,double>, 2>>;
159 double xphy, yphy, zphy; 159 double xphy, yphy, zphy;
160 double xChe, yChe, zChe; 160 double xChe, yChe, zChe;
161 double xrad, yrad, zrad; 161 double xrad, yrad, zrad;
162 162
163 map<int,int> CopyNumTable; 163 map<int,int> CopyNumTable;
164 Table physTable; 164 Table physTable;
165 165
166 double EnergyDeposit; 166 double EnergyDeposit;
167 double kineticEnergyDifference; 167 double kineticEnergyDifference;
168 int flagVolume; 168 int flagVolume;
169 double copyN; 169 double copyN;
170 int EventIDp; 170 int EventIDp;
171 int EventIDc; 171 int EventIDc;
172 172
173 TFile f("output.root"); 173 TFile f("output.root");
174 174
175 // Load the file, the directory and the ntuple 175 // Load the file, the directory and the ntuple
176 TDirectoryFile* d = dynamic_cast<TDirector 176 TDirectoryFile* d = dynamic_cast<TDirectoryFile*>(f.Get("ntuple") );
177 //analyse Physical stage 177 //analyse Physical stage
178 TTree* tPhys = dynamic_cast<TTree*> (d->Ge 178 TTree* tPhys = dynamic_cast<TTree*> (d->Get("ntuple_1") );
179 tPhys->SetBranchAddress("x",&xphy); 179 tPhys->SetBranchAddress("x",&xphy);
180 tPhys->SetBranchAddress("y",&yphy); 180 tPhys->SetBranchAddress("y",&yphy);
181 tPhys->SetBranchAddress("z",&zphy); 181 tPhys->SetBranchAddress("z",&zphy);
182 tPhys->SetBranchAddress("edep",&EnergyDepo 182 tPhys->SetBranchAddress("edep",&EnergyDeposit);
183 tPhys->SetBranchAddress("diffKin",&kinetic 183 tPhys->SetBranchAddress("diffKin",&kineticEnergyDifference);
184 tPhys->SetBranchAddress("volumeName",&flag 184 tPhys->SetBranchAddress("volumeName",&flagVolume);
185 tPhys->SetBranchAddress("CopyNumber",© 185 tPhys->SetBranchAddress("CopyNumber",©N);
186 tPhys->SetBranchAddress("EventID",&EventID 186 tPhys->SetBranchAddress("EventID",&EventIDp);
187 187
188 auto entryPhyNumber = tPhys->GetEntries(); 188 auto entryPhyNumber = tPhys->GetEntries();
189 189
190 bool addCopyN = true; 190 bool addCopyN = true;
191 for(decltype(entryPhyNumber) i = 0; i < en 191 for(decltype(entryPhyNumber) i = 0; i < entryPhyNumber; ++i)
192 { 192 {
193 tPhys->GetEntry(i); 193 tPhys->GetEntry(i);
194 //avoid histone 194 //avoid histone
195 if(flagVolume == DNAVolumeType::histon 195 if(flagVolume == DNAVolumeType::histone)
196 { 196 {
197 continue; 197 continue;
198 } 198 }
199 //Determine the strand 199 //Determine the strand
200 int strand = -1; 200 int strand = -1;
201 bool isLeftStrand = DNAVolumeType::deo 201 bool isLeftStrand = DNAVolumeType::deoxyribose1 == flagVolume ||
202 DNAVolumeType::pho 202 DNAVolumeType::phosphate1 == flagVolume ||
203 DNAVolumeType::deo 203 DNAVolumeType::deoxyribose1_water == flagVolume ||
204 DNAVolumeType::pho 204 DNAVolumeType::phosphate1_water == flagVolume;
205 205
206 bool isRightStrand = DNAVolumeType::de 206 bool isRightStrand = DNAVolumeType::deoxyribose2 == flagVolume ||
207 DNAVolumeType::ph 207 DNAVolumeType::phosphate2 == flagVolume ||
208 DNAVolumeType::de 208 DNAVolumeType::deoxyribose2_water== flagVolume ||
209 DNAVolumeType::ph 209 DNAVolumeType::phosphate2_water == flagVolume;
210 210
211 211
212 if( isLeftStrand ) 212 if( isLeftStrand )
213 { 213 {
214 strand = 0; 214 strand = 0;
215 } 215 }
216 else if( isRightStrand ) 216 else if( isRightStrand )
217 { 217 {
218 strand = 1; 218 strand = 1;
219 } 219 }
220 else 220 else
221 { 221 {
222 //in case molecules are not assign 222 //in case molecules are not assigned "strand"
223 continue; 223 continue;
224 } 224 }
225 225
226 //Determine the name 226 //Determine the name
227 bool isDeoxyribode = flagVolume == DNA 227 bool isDeoxyribode = flagVolume == DNAVolumeType::deoxyribose1 ||
228 flagVolume == DNA 228 flagVolume == DNAVolumeType::deoxyribose2 ||
229 flagVolume == DNA 229 flagVolume == DNAVolumeType::deoxyribose1_water ||
230 flagVolume == DNA 230 flagVolume == DNAVolumeType::deoxyribose2_water;
231 231
232 bool isPhosphate = flagVolume == DNA 232 bool isPhosphate = flagVolume == DNAVolumeType::phosphate1 ||
233 flagVolume == DNA 233 flagVolume == DNAVolumeType::phosphate2 ||
234 flagVolume == DNA 234 flagVolume == DNAVolumeType::phosphate1_water||
235 flagVolume == DNA 235 flagVolume == DNAVolumeType::phosphate2_water;
236 236
237 if(isDeoxyribode || isPhosphate) 237 if(isDeoxyribode || isPhosphate)
238 { 238 {
239 physTable[EventIDp][strand][copyN] 239 physTable[EventIDp][strand][copyN] += EnergyDeposit;
240 } 240 }
241 241
242 if(physTable[EventIDp][strand][copyN] 242 if(physTable[EventIDp][strand][copyN] < Ethresh)
243 { 243 {
244 continue; 244 continue;
245 } 245 }
246 246
247 if(CopyNumTable.empty()) 247 if(CopyNumTable.empty())
248 { 248 {
249 CopyNumTable.insert(pair<int,int>( 249 CopyNumTable.insert(pair<int,int>(copyN,strand));
250 if (fDamagesEventMap.find(EventIDp 250 if (fDamagesEventMap.find(EventIDp) == fDamagesEventMap.end()) {
251 std::vector<Damage> admg{Damag 251 std::vector<Damage> admg{Damage(strand,copyN,0)};
252 fDamagesEventMap.insert({Event 252 fDamagesEventMap.insert({EventIDp,admg});
253 } else { 253 } else {
254 fDamagesEventMap[EventIDp].pus 254 fDamagesEventMap[EventIDp].push_back(Damage(strand,copyN,0));
255 } 255 }
256 totalDamages++; 256 totalDamages++;
257 } 257 }
258 else 258 else
259 { 259 {
260 addCopyN = true; 260 addCopyN = true;
261 } 261 }
262 262
263 auto itCopyNumTable = CopyNumTable.fin 263 auto itCopyNumTable = CopyNumTable.find(copyN);
264 if (itCopyNumTable != CopyNumTable.end 264 if (itCopyNumTable != CopyNumTable.end())
265 { 265 {
266 if (itCopyNumTable->second == stra 266 if (itCopyNumTable->second == strand)
267 { 267 {
268 addCopyN = false; 268 addCopyN = false;
269 } 269 }
270 } 270 }
271 271
272 if(addCopyN) 272 if(addCopyN)
273 { 273 {
274 CopyNumTable.insert(pair<int,int>( 274 CopyNumTable.insert(pair<int,int>(copyN,strand));
275 if (fDamagesEventMap.find(EventIDp 275 if (fDamagesEventMap.find(EventIDp) == fDamagesEventMap.end()) {
276 std::vector<Damage> admg{Damag 276 std::vector<Damage> admg{Damage(strand,copyN,0)};
277 fDamagesEventMap.insert({Event 277 fDamagesEventMap.insert({EventIDp,admg});
278 } else { 278 } else {
279 fDamagesEventMap[EventIDp].pus 279 fDamagesEventMap[EventIDp].push_back(Damage(strand,copyN,0));
280 } 280 }
281 totalDamages++; 281 totalDamages++;
282 } 282 }
283 } 283 }
284 //Chemistry analyse 284 //Chemistry analyse
285 TTree* tChem = dynamic_cast<TTree*> (d->Ge 285 TTree* tChem = dynamic_cast<TTree*> (d->Get("ntuple_2") );
286 tChem->SetBranchAddress("x",&xrad); 286 tChem->SetBranchAddress("x",&xrad);
287 tChem->SetBranchAddress("y",&yrad); 287 tChem->SetBranchAddress("y",&yrad);
288 tChem->SetBranchAddress("z",&zrad); 288 tChem->SetBranchAddress("z",&zrad);
289 tChem->SetBranchAddress("EventID",&EventID 289 tChem->SetBranchAddress("EventID",&EventIDc);
290 290
291 auto entryNChem = tChem->GetEntries() ; 291 auto entryNChem = tChem->GetEntries() ;
292 292
293 for(decltype(entryNChem) i=0; i < entryNCh 293 for(decltype(entryNChem) i=0; i < entryNChem; ++i)
294 { 294 {
295 tChem->GetEntry(i); 295 tChem->GetEntry(i);
296 ThreeVector<double> DNAElement(xrad,yr 296 ThreeVector<double> DNAElement(xrad,yrad,zrad);
297 297
298 for(const auto& moleculeIt : fMolecule 298 for(const auto& moleculeIt : fMolecules)
299 { 299 {
300 if(moleculeIt.fPosition == DNAElem 300 if(moleculeIt.fPosition == DNAElement)
301 { 301 {
302 int strand = -1; 302 int strand = -1;
303 if(moleculeIt.fName == "deoxyr 303 if(moleculeIt.fName == "deoxyribose1")
304 { 304 {
305 strand = 0; 305 strand = 0;
306 } 306 }
307 else if(moleculeIt.fName == "d 307 else if(moleculeIt.fName == "deoxyribose2")
308 { 308 {
309 strand = 1; 309 strand = 1;
310 } 310 }
311 else 311 else
312 { 312 {
313 string msg = "Unknown DNA 313 string msg = "Unknown DNA component";
314 throw std::runtime_error(m 314 throw std::runtime_error(msg);
315 } 315 }
316 316
317 if(rdomN.Rndm() > ChemPro)//pr 317 if(rdomN.Rndm() > ChemPro)//probability of reaction make damages
318 { 318 {
319 continue; 319 continue;
320 } 320 }
321 321
322 if (!CopyNumTable.empty()) 322 if (!CopyNumTable.empty())
323 { 323 {
324 addCopyN = true; 324 addCopyN = true;
325 } 325 }
326 326
327 auto itCopyNumTable = CopyNumT 327 auto itCopyNumTable = CopyNumTable.find(moleculeIt.fCopyNumber);
328 328
329 if (itCopyNumTable != CopyNumT 329 if (itCopyNumTable != CopyNumTable.end())
330 { 330 {
331 if (itCopyNumTable->second 331 if (itCopyNumTable->second == strand)
332 { 332 {
333 addCopyN = false; 333 addCopyN = false;
334 } 334 }
335 } 335 }
336 if(addCopyN) 336 if(addCopyN)
337 { 337 {
338 if (fDamagesEventMap.find( 338 if (fDamagesEventMap.find(EventIDc) == fDamagesEventMap.end()) {
339 std::vector<Damage> ad 339 std::vector<Damage> admg{Damage(strand,moleculeIt.fCopyNumber,1)};
340 fDamagesEventMap.inser 340 fDamagesEventMap.insert({EventIDc,admg});
341 } else { 341 } else {
342 fDamagesEventMap[Event 342 fDamagesEventMap[EventIDc].push_back(
343 343 Damage(strand,moleculeIt.fCopyNumber,1));
344 } 344 }
345 totalDamages++; 345 totalDamages++;
346 } 346 }
347 } 347 }
348 } 348 }
349 } 349 }
350 } 350 }
351 351
352 //....oooOO0OOooo........oooOO0OOooo........oo 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353 353
354 void CreateSDDfile(std::map<int,std::vector<Da 354 void CreateSDDfile(std::map<int,std::vector<Damage>> fDamagesEventMap,long int totalDamages)
355 { 355 {
356 ofstream ofile("SDDFormat.txt"); 356 ofstream ofile("SDDFormat.txt");
357 if(ofile.is_open()) 357 if(ofile.is_open())
358 { 358 {
359 ofile<<'\n' 359 ofile<<'\n'
360 <<"SDD version, SDDv1.0;\n" 360 <<"SDD version, SDDv1.0;\n"
361 <<"Software, chromatin fiber Model 361 <<"Software, chromatin fiber Model;\n"
362 <<"Author contact, Carmen Villagra 362 <<"Author contact, Carmen Villagrasa,carmen.villagrasa@irsn.fr, "
363 "02/04/2018, Melyn et al.,Sc. Re 363 "02/04/2018, Melyn et al.,Sc. Rep. 7 (2017)11923;\n"
364 <<"Simulation Details, DNA damages 364 <<"Simulation Details, DNA damages from direct and "
365 "indirect effects;\n" 365 "indirect effects;\n"
366 <<"Source, Monoenergetic cylinder- 366 <<"Source, Monoenergetic cylinder-parallel proton beam uniformly "
367 "in 5 nm radius from left side o 367 "in 5 nm radius from left side of the cube"
368 <<" exposing chromatin fiber. Ener 368 <<" exposing chromatin fiber. Energy: 5 keV;\n"
369 <<"Source type, 1;\n" 369 <<"Source type, 1;\n"
370 <<"Incident particles, 2212;\n" 370 <<"Incident particles, 2212;\n"
371 <<"Mean particle energy, 5 keV;\n" 371 <<"Mean particle energy, 5 keV;\n"
372 <<"Energy distribution, M, 0;\n" 372 <<"Energy distribution, M, 0;\n"
373 <<"Particle fraction, 1.0;\n" 373 <<"Particle fraction, 1.0;\n"
374 <<"Dose or fluence, 0, 0;\n" 374 <<"Dose or fluence, 0, 0;\n"
375 <<"Dose rate, 0.0;\n" 375 <<"Dose rate, 0.0;\n"
376 <<"Irradiation target, Simple sphe 376 <<"Irradiation target, Simple spherical chromatin fiber model"
377 " in a voxel 40 nm;\n" 377 " in a voxel 40 nm;\n"
378 <<"Volumes, 0,20,20,20,-20,-20,-20 378 <<"Volumes, 0,20,20,20,-20,-20,-20,2,15,15,20,-15,-15,20;\n"
379 <<"Chromosome sizes, ;\n" 379 <<"Chromosome sizes, ;\n"
380 <<"DNA Density, ;" 380 <<"DNA Density, ;"
381 <<"Cell Cycle Phase, G0/G1;\n" 381 <<"Cell Cycle Phase, G0/G1;\n"
382 <<"DNA Structure, 4, 1;\n" 382 <<"DNA Structure, 4, 1;\n"
383 <<"In vitro / in vivo, 0;\n" 383 <<"In vitro / in vivo, 0;\n"
384 <<"Proliferation status, 1;\n" 384 <<"Proliferation status, 1;\n"
385 <<"Microenvironment, 27, 0.0;\n" 385 <<"Microenvironment, 27, 0.0;\n"
386 <<"Damage definition, 1, 1, 10, -1 386 <<"Damage definition, 1, 1, 10, -1, 17.5;\n"
387 <<"Time, 2.5ns;\n" 387 <<"Time, 2.5ns;\n"
388 <<"Damage and primary count, "<<to 388 <<"Damage and primary count, "<<totalDamages<<", 0;\n"
389 <<"Data entries, 1, 0, 1, 1, 1, 1, 389 <<"Data entries, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0;\n"
390 <<"Data field explaination, Field 390 <<"Data field explaination, Field 1: [2]-eventID, Field 3: [4]-strand, "
391 <<"Field 4:damage position (copynb 391 <<"Field 4:damage position (copynb), Field 5: Cause (direct: [1]=0) "
392 <<" (indirect: [1]=1), Field 6: Da 392 <<" (indirect: [1]=1), Field 6: Damage types (Base:[1]>0) (Backbone: [2]>0);\n"
393 <<"\n" 393 <<"\n"
394 <<"***EndOfHeader***;\n" 394 <<"***EndOfHeader***;\n"
395 <<endl; 395 <<endl;
396 //For detail, please refer to Schuemann, J., e 396 //For detail, please refer to Schuemann, J., et al. (2019). A New Standard DNA
397 //Damage (SDD) Data Format. Radiation Research 397 //Damage (SDD) Data Format. Radiation Research, 191(1), 76.
398 398
399 bool Primaryflag = true; 399 bool Primaryflag = true;
400 for (auto const [eventID,dmgvector] : fDam 400 for (auto const [eventID,dmgvector] : fDamagesEventMap)
401 { 401 {
402 Primaryflag = true; // update Primaryf 402 Primaryflag = true; // update Primaryflag for new primary
403 for (auto const dmg:dmgvector) { 403 for (auto const dmg:dmgvector) {
404 //Field 1 Calassification 404 //Field 1 Calassification
405 int newEvtFlag = 0; // = 2 if new 405 int newEvtFlag = 0; // = 2 if new event;
406 if (Primaryflag) { 406 if (Primaryflag) {
407 newEvtFlag = 2; 407 newEvtFlag = 2;
408 Primaryflag = false; 408 Primaryflag = false;
409 } 409 }
410 ofile<<newEvtFlag<<", "<<eventID<< 410 ofile<<newEvtFlag<<", "<<eventID<<"; ";
411 //Field 3 Chromosome IDs and stran 411 //Field 3 Chromosome IDs and strand
412 ofile<<0<<", "<<0<<", "<<0<<", "<< 412 ofile<<0<<", "<<0<<", "<<0<<", "<<dmg.strand<<"; ";
413 //Field 4, Chromosome position 413 //Field 4, Chromosome position
414 ofile<<dmg.copyNb<<"; "; 414 ofile<<dmg.copyNb<<"; ";
415 //Field 5, Cause: Unknown = -1, Di 415 //Field 5, Cause: Unknown = -1, Direct = 0, Indirect = 1
416 ofile<<dmg.cause<<", "<<0<<", "<<0 416 ofile<<dmg.cause<<", "<<0<<", "<<0<<"; ";
417 //Field 6, Damage types: only Back 417 //Field 6, Damage types: only Backbone is considered here
418 ofile<<0<<", "<<1<<", "<<0<<"; "; 418 ofile<<0<<", "<<1<<", "<<0<<"; ";
419 ofile<<"\n"; 419 ofile<<"\n";
420 } 420 }
421 } 421 }
422 422
423 ofile<<endl; 423 ofile<<endl;
424 } 424 }
425 425
426 ofile.close(); 426 ofile.close();
427 } 427 }
428 428
429 //....oooOO0OOooo........oooOO0OOooo........oo 429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....