ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot.C
1 // -------------------------------------------------------------------
2 // $Id: plot.C
3 // -------------------------------------------------------------------
4 //
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 {
12  gROOT->Reset();
13  gStyle->SetPalette(1);
14  gROOT->SetStyle("Plain");
15  TRandom rdomN;
16 
17  vector<Molecule> fMolecules = molecule();
18  using Table = map<int,array<map<int,double>, 2>>;
19  using ResultTable = map<int,array<vector<int>, 2>>;
20 
21  system ("rm -rf output.root");
22  system ("hadd output.root output_*.root");
23 
24  double xphy;
25  double yphy;
26  double zphy;
27 
28  double xChe;
29  double yChe;
30  double zChe;
31 
32  double xrad;
33  double yrad;
34  double zrad;
35 
36  map<int,int> CopyNumTable;
37  Table physTable;
38  ResultTable resultPhysTable;
39  ResultTable resultChemTable;
40  ResultTable MergedResultTable;
41 
42  double EnergyDeposit;
43  double kineticEnergyDifference;
44  int flagVolume;
45  double copyN;
46  int EventIDp;
47  int EventIDc;
48 
49  TFile f("output.root");
50 
51 // Load the file, the directory and the ntuple
52  TDirectoryFile* d = dynamic_cast<TDirectoryFile*>(f.Get("ntuple") );
53  TTree* tPhys = dynamic_cast<TTree*> (d->Get("ntuple_1") );
54  tPhys->SetBranchAddress("x",&xphy);
55  tPhys->SetBranchAddress("y",&yphy);
56  tPhys->SetBranchAddress("z",&zphy);
57  tPhys->SetBranchAddress("edep",&EnergyDeposit);
58  tPhys->SetBranchAddress("diffKin",&kineticEnergyDifference);
59  tPhys->SetBranchAddress("volumeName",&flagVolume);
60  tPhys->SetBranchAddress("CopyNumber",&copyN);
61  tPhys->SetBranchAddress("EventID",&EventIDp);
62 
63  auto entryPhyNumber = tPhys->GetEntries();
64 
65  bool addCopyN = true;
66  for(decltype(entryPhyNumber) i = 0; i < entryPhyNumber; ++i)
67  {
68  tPhys->GetEntry(i);
69 //avoid histone
70  if(flagVolume == DNAVolumeType::histone)
71  {
72  continue;
73  }
74 //Determine the strand
75  int strand = -1;
76  bool isLeftStrand = DNAVolumeType::deoxyribose1 == flagVolume ||
77  DNAVolumeType::phosphate1 == flagVolume ||
78  DNAVolumeType::deoxyribose1_water == flagVolume ||
79  DNAVolumeType::phosphate1_water == flagVolume;
80 
81  bool isRightStrand = DNAVolumeType::deoxyribose2 == flagVolume ||
82  DNAVolumeType::phosphate2 == flagVolume ||
83  DNAVolumeType::deoxyribose2_water== flagVolume ||
84  DNAVolumeType::phosphate2_water == flagVolume;
85 
86 
87  if( isLeftStrand )
88  {
89  strand = 0;
90  }
91  else if( isRightStrand )
92  {
93  strand = 1;
94  }
95  else
96  {
97  //in case molecules are not assigned "strand"
98  continue;
99  }
100 
101 //Determine the name
102  bool isDeoxyribode = flagVolume == DNAVolumeType::deoxyribose1 ||
103  flagVolume == DNAVolumeType::deoxyribose2 ||
104  flagVolume == DNAVolumeType::deoxyribose1_water ||
105  flagVolume == DNAVolumeType::deoxyribose2_water;
106 
107  bool isPhosphate = flagVolume == DNAVolumeType::phosphate1 ||
108  flagVolume == DNAVolumeType::phosphate2 ||
109  flagVolume == DNAVolumeType::phosphate1_water||
110  flagVolume == DNAVolumeType::phosphate2_water;
111 
112  if(isDeoxyribode || isPhosphate)
113  {
114  physTable[EventIDp][strand][copyN] += EnergyDeposit;
115  }
116 
117  if(physTable[EventIDp][strand][copyN] < 17.5)
118  {
119  continue;
120  }
121 
122  if(CopyNumTable.empty())
123  {
124  CopyNumTable.insert(pair<int,int>(copyN,strand));
125  resultPhysTable[EventIDp][strand].push_back(copyN);
126  }
127  else
128  {
129  addCopyN = true;
130  }
131 
132  auto itCopyNumTable = CopyNumTable.find(copyN);
133  if (itCopyNumTable != CopyNumTable.end())
134  {
135  if (itCopyNumTable->second == strand)
136  {
137  addCopyN = false;
138  }
139  }
140 
141  if(addCopyN)
142  {
143  CopyNumTable.insert(pair<int,int>(copyN,strand));
144  resultPhysTable[EventIDp][strand].push_back(copyN);
145  }
146  }
147 //Chemistry analyse
148  TTree* tChem = dynamic_cast<TTree*> (d->Get("ntuple_2") );
149  tChem->SetBranchAddress("x",&xrad);
150  tChem->SetBranchAddress("y",&yrad);
151  tChem->SetBranchAddress("z",&zrad);
152  tChem->SetBranchAddress("EventID",&EventIDc);
153 
154  auto entryNChem = tChem->GetEntries() ;
155 
156  for(decltype(entryNChem) i=0; i < entryNChem; ++i)
157  {
158  tChem->GetEntry(i);
159  ThreeVector<double> DNAElement(xrad,yrad,zrad);
160 
161  for(const auto& moleculeIt : fMolecules)
162  {
163  if(moleculeIt.fPosition == DNAElement)
164  {
165  int strand = -1;
166  if(moleculeIt.fName == "deoxyribose1")
167  {
168  strand = 0;
169  }
170  else if(moleculeIt.fName == "deoxyribose2")
171  {
172  strand = 1;
173  }
174  else
175  {
176  string msg = "Unknown DNA component";
177  throw std::runtime_error(msg);
178  }
179 
180  if(rdomN.Rndm() > 0.42)//42% of reaction make damages
181  {
182  continue;
183  }
184 
185  if (!CopyNumTable.empty())
186  {
187  addCopyN = true;
188  }
189 
190  auto itCopyNumTable = CopyNumTable.find(moleculeIt.fCopyNumber);
191 
192  if (itCopyNumTable != CopyNumTable.end())
193  {
194  if (itCopyNumTable->second == strand)
195  {
196  addCopyN = false;
197  }
198  }
199  if(addCopyN)
200  {
201  resultChemTable[EventIDc][strand].push_back(
202  moleculeIt.fCopyNumber);
203  }
204  }
205  }
206  }
207 //SDD Format
208  ofstream ofile("SDDFormat.txt");
209  if(ofile.is_open())
210  {
211  ofile<<'\n'
212  <<"SDD version, SDDv1.0;\n"
213  <<"Software, chromatin fiber Model;\n"
214  <<"Author contact, Carmen Villagrasa,carmen.villagrasa@irsn.fr, "
215  "02/04/2018, Melyn et al.,Sc. Rep. 7 (2017)11923;\n"
216  <<"Simulation Details, DNA damages from direct and "
217  "indirect effects;\n"
218  <<"Source, Monoenergetic cylinder-parallel proton beam uniformly "
219  "in 5 nm radius from left side of the cube"
220  <<" exposing chromatin fiber. Energy: 5 keV;\n"
221  <<"Source type, 1;\n"
222  <<"Incident particles, 200;\n"
223  <<"Mean particle energy, 5 keV;\n"
224  <<"Energy distribution, M, 0;\n"
225  <<"Particle fraction, 1.0;\n"
226  <<"Dose or fluence, 0, 0;\n"
227  <<"Dose rate, 0.0;\n"
228  <<"Irradiation target, Simple spherical chromatin fiber model"
229  " in a voxel 40 nm;\n"
230  <<"Volumes, 0,20,20,20,-20,-20,-20,2,15,15,20,-15,-15,20;\n"
231  <<"Chromosome sizes, ;\n"
232  <<"DNA Density, ;"
233  <<"Cell Cycle Phase, G1/G2;\n"
234  <<"DNA Structure, 4, 1;\n"
235  <<"In vitro / in vivo, 0;\n"
236  <<"Proliferation status, 1;\n"
237  <<"Microenvironment, 27, 0.0;\n"
238  <<"Damage definition, 1, 1, 10, -1, 0, 17.5;\n"
239  <<"Time, 2.5ns;\n"
240  <<"Damage and primary count, 638, 200;\n"
241  <<"Data entries, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0;\n"
242  <<"Data was generated in minimal output format and used as"
243  " an example. Please modify the primary particle number"
244  "for this file;\n"
245  <<"\n"
246  <<"***EndOfHeader***;\n"
247  <<endl;
248 //For detail, please refer to Schuemann, J., et al. (2019). A New Standard DNA
249 //Damage (SDD) Data Format. Radiation Research, 191(1), 76.
250 
251  bool Primaryflag = true;
252  for (int i = 0; i < 200; i++)//for 200 primary particles
253  {
254  for(int j = 0; j < 2; j++)//for strand
255  {
256  const std::string positionAndCopyN = "; 0, 0, 0; 3, 0, 0, ";
257  const std::string SSBType = " 0, 1, 0";
258  if (!resultPhysTable[i][j].empty())
259  {
260  for(int e = 0; e < resultPhysTable[i][j].size(); ++e)
261  {
262  if(Primaryflag)
263  {
264  ofile<<"2, "<<i<<positionAndCopyN
265  <<resultPhysTable[i][j][e]
266  <<"; 0;"<<SSBType<<"\n" ;//physics
267  Primaryflag = false;
268  }else
269  {
270  ofile<<"0, "<<i<<positionAndCopyN
271  <<resultPhysTable[i][j][e]
272  <<"; 0;"<<SSBType<<"\n" ;
273  }
274  }
275  }
276  if(resultChemTable[i][j].empty())
277  {
278  continue;
279  }
280 
281  for(int ee = 0; ee < resultChemTable[i][j].size(); ++ee)
282  {
283  if(Primaryflag)
284  {
285  ofile<<"2, "<<i<<positionAndCopyN
286  <<resultChemTable[i][j][ee]
287  <<"; 1;"<<SSBType<<"\n" ;//chemistry
288  Primaryflag = false;
289  }else
290  {
291  ofile<<"0, "<<i<<positionAndCopyN
292  <<resultChemTable[i][j][ee]
293  <<"; 1;"<<SSBType<<"\n" ;
294  }
295  }
296  }
297  }
298 
299  ofile<<endl;
300  }
301 
302  ofile.close();
303 }