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