| Geant4 Cross Reference |
1 // ------------------------------------------- 1
2 // $Id: scandamages.C
3 // -------------------------------------------
4 // 13/10/2023 - Le Tuan Anh renamed and improv
5 // *******************************************
6 // To execute this macro under ROOT after your
7 // 1 - launch ROOT v6.xx (usually type 'root
8 // 2 - type '.X analysis.C' at the ROOT sess
9 // *******************************************
10
11 //....oooOO0OOooo........oooOO0OOooo........oo
12
13 struct Damage;
14 struct Cluster;
15 void CreateSDDfile(std::map<int,std::vector<Da
16 void ExtractDirectAndIndirectDamages(std::map<
17 void ClassifyDamages(std::map<int,std::vector<
18 std::map<int,vector<Cluste
19
20 //....oooOO0OOooo........oooOO0OOooo........oo
21 // play with scandamages() by changing firect/
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
32 double fProbabilityForIndirectionSelection
33 int fClusterdistance = 10; // bp, minium d
34
35 std::map<int,std::vector<Damage>> fDamages
36 std::map<int,vector<Cluster>> fClustersMap
37 long int totalDamages = 0;
38
39 ExtractDirectAndIndirectDamages(fDamagesEv
40 fEnergyThr
41 fProbabili
42
43 ClassifyDamages(fDamagesEventMap,fClusterd
44 // print results in each event for checkin
45 //ClassifyDamages(fDamagesEventMap,fCluste
46
47 CreateSDDfile(fDamagesEventMap,totalDamage
48 }
49
50 //....oooOO0OOooo........oooOO0OOooo........oo
51
52 struct Damage
53 {
54 Damage(int a, int b, int c=-1): strand(a),
55 int strand;
56 int copyNb;
57 int cause; // Unknown = -1, Direct = 0, In
58 bool operator<(const Damage& r) const {ret
59 };
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 struct Cluster // to score SSB or DSB
64 {
65 std::vector<Damage> fDamagesVector;
66 bool HasLeftStrand=false, HasRightStrand=f
67 int firstCopyNb;
68 };
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71
72 void
73 ClassifyDamages(std::map<int,std::vector<Damag
74 std::map<int,vector<Cluster>>&
75 {
76 long int dirSB = 0, indirSB = 0;
77 long int nDSB = 0, nsDSB = 0, ncDSB=0, nSS
78 std::cout<<"------------------------------
79 if (verbose>0)std::cout<<" ------> Print f
80 if (verbose>0)std::cout<<"\t\tCopyNb\tStra
81 for (auto &[eventID, dmgVector] : fDamages
82 std::sort(dmgVector.begin(),dmgVector.
83 if (verbose>0) std::cout<<"Event #"<<e
84 unsigned long fSBcountThisEvent = 0;
85 long copyNBofPrevDamage = 0;
86 long firstCopyNb;
87 std::map<int,Cluster> clustersInthisEv
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"<<
92 <<dmg.stra
93 if (clustersInthisEvent.size()==0
94 ((dmg.copyNb - copyNBofPrevDam
95 Cluster newCluster;
96 newCluster.fDamagesVector.push
97 copyNBofPrevDamage = dmg.copyN
98 firstCopyNb = dmg.copyNb;
99 if (dmg.strand == 0) newCluste
100 if (dmg.strand == 1) newCluste
101 clustersInthisEvent.insert({fi
102 } else {
103 clustersInthisEvent[firstCopyN
104 copyNBofPrevDamage = dmg.copyN
105 if (dmg.strand == 0 ) clusters
106 if (dmg.strand == 1 ) clusters
107 }
108 }
109
110 for (auto [firstCpNb,acluster] : clust
111 if (fClustersMap.find(eventID) ==
112 std::vector<Cluster> sclv{aclu
113 fClustersMap.insert({eventID,s
114 } else {
115 fClustersMap[eventID].push_bac
116 }
117 }
118 clustersInthisEvent.clear();
119 }
120 if (verbose>0)std::cout<<" \n------> Check
121 if (verbose>0)std::cout<<"\t\t#SSB\t#DSB\t
122
123 for (auto [evt,clustervector] : fClustersM
124 long int nDSBinThisEvt = 0, nSSBinThis
125 long int nsDSBinThisEvt = 0, ncDSBinTh
126 if (verbose>0) std::cout<<"Event #"<<e
127 for (auto sb : clustervector) {
128 if (sb.HasLeftStrand && sb.HasRigh
129 nDSBinThisEvt++;
130 if (sb.fDamagesVector.size()>2
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"<<nSSB
140 <<"\t"<<ncDSBi
141 }
142 std::cout<<"------------------------------
143 std::cout<<"\tSummary of results:"
144 <<"\n #Total SBs = "<<(dirSB+i
145 <<"\n #dirSB = "<<dirSB <<
146 <<"\n #SSB = "<<nSSB <<
147 <<"\n #sDSB = "<<nsDSB <<
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oo
151
152 void ExtractDirectAndIndirectDamages(std::map<
153 long int &totalDamages,double Ethresh,double C
154 {
155 TRandom rdomN;
156
157 vector<Molecule> fMolecules = molecule();
158 using Table = map<int,array<map<int,double
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<TDirector
177 //analyse Physical stage
178 TTree* tPhys = dynamic_cast<TTree*> (d->Ge
179 tPhys->SetBranchAddress("x",&xphy);
180 tPhys->SetBranchAddress("y",&yphy);
181 tPhys->SetBranchAddress("z",&zphy);
182 tPhys->SetBranchAddress("edep",&EnergyDepo
183 tPhys->SetBranchAddress("diffKin",&kinetic
184 tPhys->SetBranchAddress("volumeName",&flag
185 tPhys->SetBranchAddress("CopyNumber",©
186 tPhys->SetBranchAddress("EventID",&EventID
187
188 auto entryPhyNumber = tPhys->GetEntries();
189
190 bool addCopyN = true;
191 for(decltype(entryPhyNumber) i = 0; i < en
192 {
193 tPhys->GetEntry(i);
194 //avoid histone
195 if(flagVolume == DNAVolumeType::histon
196 {
197 continue;
198 }
199 //Determine the strand
200 int strand = -1;
201 bool isLeftStrand = DNAVolumeType::deo
202 DNAVolumeType::pho
203 DNAVolumeType::deo
204 DNAVolumeType::pho
205
206 bool isRightStrand = DNAVolumeType::de
207 DNAVolumeType::ph
208 DNAVolumeType::de
209 DNAVolumeType::ph
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 assign
223 continue;
224 }
225
226 //Determine the name
227 bool isDeoxyribode = flagVolume == DNA
228 flagVolume == DNA
229 flagVolume == DNA
230 flagVolume == DNA
231
232 bool isPhosphate = flagVolume == DNA
233 flagVolume == DNA
234 flagVolume == DNA
235 flagVolume == DNA
236
237 if(isDeoxyribode || isPhosphate)
238 {
239 physTable[EventIDp][strand][copyN]
240 }
241
242 if(physTable[EventIDp][strand][copyN]
243 {
244 continue;
245 }
246
247 if(CopyNumTable.empty())
248 {
249 CopyNumTable.insert(pair<int,int>(
250 if (fDamagesEventMap.find(EventIDp
251 std::vector<Damage> admg{Damag
252 fDamagesEventMap.insert({Event
253 } else {
254 fDamagesEventMap[EventIDp].pus
255 }
256 totalDamages++;
257 }
258 else
259 {
260 addCopyN = true;
261 }
262
263 auto itCopyNumTable = CopyNumTable.fin
264 if (itCopyNumTable != CopyNumTable.end
265 {
266 if (itCopyNumTable->second == stra
267 {
268 addCopyN = false;
269 }
270 }
271
272 if(addCopyN)
273 {
274 CopyNumTable.insert(pair<int,int>(
275 if (fDamagesEventMap.find(EventIDp
276 std::vector<Damage> admg{Damag
277 fDamagesEventMap.insert({Event
278 } else {
279 fDamagesEventMap[EventIDp].pus
280 }
281 totalDamages++;
282 }
283 }
284 //Chemistry analyse
285 TTree* tChem = dynamic_cast<TTree*> (d->Ge
286 tChem->SetBranchAddress("x",&xrad);
287 tChem->SetBranchAddress("y",&yrad);
288 tChem->SetBranchAddress("z",&zrad);
289 tChem->SetBranchAddress("EventID",&EventID
290
291 auto entryNChem = tChem->GetEntries() ;
292
293 for(decltype(entryNChem) i=0; i < entryNCh
294 {
295 tChem->GetEntry(i);
296 ThreeVector<double> DNAElement(xrad,yr
297
298 for(const auto& moleculeIt : fMolecule
299 {
300 if(moleculeIt.fPosition == DNAElem
301 {
302 int strand = -1;
303 if(moleculeIt.fName == "deoxyr
304 {
305 strand = 0;
306 }
307 else if(moleculeIt.fName == "d
308 {
309 strand = 1;
310 }
311 else
312 {
313 string msg = "Unknown DNA
314 throw std::runtime_error(m
315 }
316
317 if(rdomN.Rndm() > ChemPro)//pr
318 {
319 continue;
320 }
321
322 if (!CopyNumTable.empty())
323 {
324 addCopyN = true;
325 }
326
327 auto itCopyNumTable = CopyNumT
328
329 if (itCopyNumTable != CopyNumT
330 {
331 if (itCopyNumTable->second
332 {
333 addCopyN = false;
334 }
335 }
336 if(addCopyN)
337 {
338 if (fDamagesEventMap.find(
339 std::vector<Damage> ad
340 fDamagesEventMap.inser
341 } else {
342 fDamagesEventMap[Event
343
344 }
345 totalDamages++;
346 }
347 }
348 }
349 }
350 }
351
352 //....oooOO0OOooo........oooOO0OOooo........oo
353
354 void CreateSDDfile(std::map<int,std::vector<Da
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
362 <<"Author contact, Carmen Villagra
363 "02/04/2018, Melyn et al.,Sc. Re
364 <<"Simulation Details, DNA damages
365 "indirect effects;\n"
366 <<"Source, Monoenergetic cylinder-
367 "in 5 nm radius from left side o
368 <<" exposing chromatin fiber. Ener
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 sphe
377 " in a voxel 40 nm;\n"
378 <<"Volumes, 0,20,20,20,-20,-20,-20
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
387 <<"Time, 2.5ns;\n"
388 <<"Damage and primary count, "<<to
389 <<"Data entries, 1, 0, 1, 1, 1, 1,
390 <<"Data field explaination, Field
391 <<"Field 4:damage position (copynb
392 <<" (indirect: [1]=1), Field 6: Da
393 <<"\n"
394 <<"***EndOfHeader***;\n"
395 <<endl;
396 //For detail, please refer to Schuemann, J., e
397 //Damage (SDD) Data Format. Radiation Research
398
399 bool Primaryflag = true;
400 for (auto const [eventID,dmgvector] : fDam
401 {
402 Primaryflag = true; // update Primaryf
403 for (auto const dmg:dmgvector) {
404 //Field 1 Calassification
405 int newEvtFlag = 0; // = 2 if new
406 if (Primaryflag) {
407 newEvtFlag = 2;
408 Primaryflag = false;
409 }
410 ofile<<newEvtFlag<<", "<<eventID<<
411 //Field 3 Chromosome IDs and stran
412 ofile<<0<<", "<<0<<", "<<0<<", "<<
413 //Field 4, Chromosome position
414 ofile<<dmg.copyNb<<"; ";
415 //Field 5, Cause: Unknown = -1, Di
416 ofile<<dmg.cause<<", "<<0<<", "<<0
417 //Field 6, Damage types: only Back
418 ofile<<0<<", "<<1<<", "<<0<<"; ";
419 ofile<<"\n";
420 }
421 }
422
423 ofile<<endl;
424 }
425
426 ofile.close();
427 }
428
429 //....oooOO0OOooo........oooOO0OOooo........oo