| 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 Analysis.cc
28 /// \brief Implementation of the Analysis clas
29 /// \file Analysis.cc
30 /// \brief Implementation of the Analysis clas
31
32 #include "Analysis.hh"
33 #include "DetectorConstruction.hh"
34 #include "G4AutoDelete.hh"
35 #include "G4Filesystem.hh"
36 #include "G4DNAMolecule.hh"
37
38
39 #include <sstream>
40 #include <fstream>
41 #ifdef USE_MPI
42 #include "G4MPImanager.hh"
43 #endif
44
45 G4ThreadLocal Analysis* the_analysis = nullptr
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 Analysis* Analysis::GetAnalysis()
50 {
51 if (!the_analysis) {
52 the_analysis = new Analysis();
53 G4AutoDelete::Register(the_analysis);
54 }
55
56 return the_analysis;
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 void Analysis::OpenFile(const G4String outFold
62 {
63 G4AnalysisManager* anManager = G4AnalysisM
64 anManager->SetDefaultFileType("root");
65 G4String fullFileName = fFileName;
66 if (gRunMode == RunningMode::Chem) {
67 G4String slash = "";
68 #if defined(_WIN32) || defined(WIN32)
69 slash= "\\";
70 #else
71 G4String slashu = "/";
72 slash = slashu;
73 #endif
74 fullFileName = outFolder+slash+fFileNa
75 }
76
77 // open output file
78 G4bool fileOpen = anManager->OpenFile(full
79 if (!fileOpen) {
80 G4cout << "\n---> HistoManager::book()
81 << G4endl;
82 return;
83 }
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oo
87
88 void Analysis::Save()
89 {
90 G4AnalysisManager* anManager = G4AnalysisM
91 anManager->Write();
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oo
95
96 void Analysis::Close(G4bool reset)
97 {
98 G4AnalysisManager* anManager = G4AnalysisM
99 anManager->CloseFile(reset);
100 }
101
102 //....oooOO0OOooo........oooOO0OOooo........oo
103
104 void Analysis::Book()
105 {
106 G4AnalysisManager* anManager = G4AnalysisM
107 anManager->SetVerboseLevel(1);
108 if (anManager->GetFirstNtupleId() != 1) an
109
110 if (gRunMode == RunningMode::Phys) {
111 #ifdef G4MULTITHREADED
112 // MT ntuple merging
113 anManager->SetNtupleMerging(true);//fo
114 #endif
115 anManager->SetNtupleDirectoryName("ntu
116 // create ntuple
117 anManager->CreateNtuple("ntuple_1","ph
118 anManager->CreateNtupleIColumn("flagPa
119 anManager->CreateNtupleIColumn("flagPa
120 anManager->CreateNtupleIColumn("flagPr
121 anManager->CreateNtupleDColumn("x");
122 anManager->CreateNtupleDColumn("y");
123 anManager->CreateNtupleDColumn("z");
124 anManager->CreateNtupleDColumn("edep")
125 anManager->CreateNtupleIColumn("eventN
126 anManager->CreateNtupleIColumn("volume
127 anManager->CreateNtupleIColumn("copyNu
128 anManager->CreateNtupleIColumn("lastMe
129 anManager->FinishNtuple(1);
130
131 // For total edep
132 anManager->CreateNtuple("ntuple_3","to
133 anManager->CreateNtupleIColumn("eventN
134 anManager->CreateNtupleDColumn("edep")
135 anManager->FinishNtuple(2);
136 }
137
138 if (gRunMode == RunningMode::Chem) {
139 // Create directories
140 const G4String directoryName = "ntuple
141 if (anManager->GetNtupleDirectoryName(
142 anManager->SetNtupleDirectoryName(
143 }
144
145 // DBScan
146
147 anManager->CreateNtuple("ntuple_2","DB
148 anManager->CreateNtupleIColumn(1,"stra
149 anManager->CreateNtupleIColumn(1,"copy
150 anManager->CreateNtupleDColumn(1,"xp")
151 anManager->CreateNtupleDColumn(1,"yp")
152 anManager->CreateNtupleDColumn(1,"zp")
153 anManager->CreateNtupleDColumn(1,"time
154 anManager->CreateNtupleIColumn(1,"base
155 anManager->FinishNtuple(1);
156 }
157 }
158
159 //....oooOO0OOooo........oooOO0OOooo........oo
160
161 G4AnalysisManager* Analysis::GetAnalysisManage
162 {
163 return G4AnalysisManager::Instance();
164 }
165
166 //....oooOO0OOooo........oooOO0OOooo........oo
167
168 void Analysis::AddInfoForChemGeo(InfoForChemGe
169 {
170 fInfoForChemGeoVector.push_back(b);
171 }
172
173 //....oooOO0OOooo........oooOO0OOooo........oo
174
175 void Analysis::ClearVector()
176 {
177 fInfoForChemGeoVector.clear();
178 fInfoInPhysStageVector.clear();
179 fOutputFiles.clear();
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oo
183
184 void Analysis::AddInfoInPhysStage(InfoInPhysSt
185 {
186 fInfoInPhysStageVector.push_back(b);
187 }
188
189 //....oooOO0OOooo........oooOO0OOooo........oo
190
191 void Analysis::UpdateChemInputDataAndFillNtupl
192 {
193 //if (fInfoForChemGeoVector.size() == 0) r
194
195 // We will loop on all the "element" of th
196 // - New file for each event/voxel couple
197 // Create a map to define all the output f
198
199 for (auto const& binfo : fInfoForChemGeoVe
200 if( binfo.fVolume == 161 // voxelStra
201 || binfo.fVolume == 162 // voxelRi
202 || binfo.fVolume == 163 // voxelLe
203 || binfo.fVolume == 164 // voxelUp
204 || binfo.fVolume == 165 // voxelDo
205 || binfo.fVolume == 261 // voxelSt
206 || binfo.fVolume == 262 // voxelRi
207 || binfo.fVolume == 263 // voxelLe
208 || binfo.fVolume == 264 // voxelUp
209 || binfo.fVolume == 265) // voxelD
210 {
211 // We are in a voxel
212
213 std::string voxelName("");
214
215 if(binfo.fVolume==161) voxelName =
216 else if(binfo.fVolume==162) voxelN
217 else if(binfo.fVolume==163) voxelN
218 else if(binfo.fVolume==164) voxelN
219 else if(binfo.fVolume==165) voxelN
220 else if(binfo.fVolume==261) voxelN
221 else if(binfo.fVolume==262) voxelN
222 else if(binfo.fVolume==263) voxelN
223 else if(binfo.fVolume==264) voxelN
224 else if(binfo.fVolume==265) voxelN
225
226 // Get the event number
227 int eventNum = int(binfo.fEventNum
228 // Here we have all the informatio
229 // We want to know if we have to c
230 // or if we already have one.
231
232 // Check if the event has already
233 if(fOutputFiles.find(eventNum)==fO
234 {
235 // If not then create the even
236 fOutputFiles[eventNum][binfo.f
237 CreateChemInputFile(eventNum,
238 }
239 // If the event has been registere
240 // registered.
241 else
242 {
243 if(fOutputFiles[eventNum].find
244 {
245 // We register the event-v
246 fOutputFiles[eventNum][bin
247 CreateChemInputFile(eventN
248 }
249 }
250
251 // Write in the file the informati
252 // the input water molecule or sol
253 UpdatingChemInputFile(binfo);
254 }
255 }
256
257 if (fInfoInPhysStageVector.size() == 0) re
258 for (auto const& binfo : fInfoInPhysStageV
259 // Check if the process is an ionisati
260 // Only ionisation should trigger the
261 if(binfo.fFlagProcess == 13
262 || binfo.fFlagProcess == 113
263 || binfo.fFlagProcess == 18
264 || binfo.fFlagProcess == 21
265 || binfo.fFlagProcess == 24
266 || binfo.fFlagProcess == 27
267 || binfo.fFlagProcess == 31) {
268 // Check the interaction happe
269 if( binfo.fVolumeName == 1 //
270 || binfo.fVolumeName == 11
271 || binfo.fVolumeName == 2
272 || binfo.fVolumeName == 22
273 || binfo.fVolumeName == 3
274 || binfo.fVolumeName == 4
275 || binfo.fVolumeName == 5
276 || binfo.fVolumeName == 6
277 || binfo.fVolumeName == 7
278 || binfo.fVolumeName == 71
279 || binfo.fVolumeName == 8
280 || binfo.fVolumeName == 81
281 || binfo.fVolumeName == 9
282 || binfo.fVolumeName == 10
283 || binfo.fVolumeName == 13
284 || binfo.fVolumeName == 12
285 {
286 // Retrieve the voxel copy
287 double voxelCopyNumber =
288 // Get the event number
289 double eventNum = binfo.fE
290
291 // Check if the event has
292 if(fOutputFiles.find(event
293 {
294 // Check if the volume
295 if(fOutputFiles.at(eve
296 != fOutputFiles.at
297 {
298 // If we are here
299 //file in which we
300 UpdatingChemInputF
301 }
302 }
303 }
304 }
305
306 // fill ntuple
307 auto analysisManager =G4AnalysisManage
308 analysisManager->FillNtupleIColumn(1,
309 analysisManager->FillNtupleIColumn(1,
310 analysisManager->FillNtupleIColumn(1,
311 analysisManager->FillNtupleDColumn(1,
312 analysisManager->FillNtupleDColumn(1,
313 analysisManager->FillNtupleDColumn(1,
314 analysisManager->FillNtupleDColumn(1,
315 analysisManager->FillNtupleIColumn(1,
316 analysisManager->FillNtupleIColumn(1,
317 analysisManager->FillNtupleIColumn(1,
318 analysisManager->FillNtupleIColumn(1,
319 analysisManager->AddNtupleRow(1);
320 }
321
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oo
325
326 G4String Analysis::CreateChemInputFile(G4int e
327 {
328 std::stringstream sstream;
329 sstream<<"./"<<fChemInputFolderName<<"/eve
330 std::ofstream oFile;
331 oFile.open(sstream.str().c_str() );
332
333 oFile<<"_eventNum"<<"\t\t"<<eventNum<<std:
334 oFile<<"_voxelType"<<"\t\t"<<voxelName<<st
335 oFile<<"_voxelCopyNumber"<<"\t\t"<<volumeC
336 oFile<<std::endl;
337
338 oFile
339 <<"# Chemistry input informations\
340 <<"# "<<"_input, type, state, elec
341 oFile<<"# type=1 -> water molecule && type
342 oFile<<std::endl;
343
344 oFile.close();
345
346 return sstream.str().c_str();
347 }
348
349 //....oooOO0OOooo........oooOO0OOooo........oo
350
351 void Analysis::UpdatingChemInputFile(InfoForCh
352 {
353 std::ofstream out;
354 out.open( fOutputFiles[b.fEventNumber][b.f
355
356 if (b.fType==1)
357 {
358 out<<"_input"<<"\t"
359 <<
360 <<
361 <<4-b.
362 <<b.fR
363 <<b.fR
364 <<b.fR
365 <<b.fParen
366 <<"\n";
367 }
368 if (b.fType==2)
369 {
370 out<<"_input"<<"\t"
371 <<
372 <<
373 <<b.fE
374 // If we are here then the event and volume co
375 //generated file in which we should add a dna
376 <<b.fR
377 <<b.fR
378 <<b.fR
379 <<b.fParen
380 <<"\n";
381 }
382 out.close();
383 }
384
385 //....oooOO0OOooo........oooOO0OOooo........oo
386
387 void Analysis::UpdatingChemInputFile(InfoInPhy
388 {
389 G4String name;
390 G4double volumeName = b.fVolumeName;
391 if(volumeName == 1 || volumeName == 2 || v
392 volumeName == 8) name = G4Deoxyri
393 else if(volumeName == 11 || volumeName ==
394 volumeName == 71 || volumeName ==
395 else if(volumeName == 6 || volumeName == 9
396 else if(volumeName == 4 || volumeName == 1
397 else if(volumeName == 5 || volumeName == 1
398 else if(volumeName == 3 || volumeName == 1
399 else
400 {
401 G4ExceptionDescription msg;
402 msg <<"Volume number "<<volumeName<<"
403 G4Exception("Analysis::UpdatingChemInp
404 "", FatalException, msg);
405 }
406
407 G4double strand (-1);
408
409 // Determine the strand
410 if(volumeName==1
411 || volumeName==11
412 || volumeName==7
413 || volumeName==71
414 || volumeName==6 // ade
415 || volumeName==9 // ade
416 || volumeName==4 // gua
417 || volumeName==10) // gua
418 strand = 1;
419 else if(volumeName==2
420 || volumeName==22
421 || volumeName==8
422 || volumeName==81
423 || volumeName==5 // thy
424 || volumeName==12 // thy
425 || volumeName==3 // cyto
426 || volumeName==13) // cyto
427 strand = 2;
428 std::ofstream out;
429 out.open(fOutputFiles[b.fEventNumber][b.fL
430 out<<"_remove"<<"\t"<<name<<"\t\t"<<b.fCop
431 out.close();
432 }
433
434 //....oooOO0OOooo........oooOO0OOooo........oo
435
436 void Analysis::WritePhysGeo()
437 {
438 //Create a file containing geopath and Phy
439 std::string fname = "imp.info";
440 std::ofstream fout(fname);
441 fout<<"====> Auto_generated file. Do not d
442 fout<<"====> The file conveys some informa
443 fout<<"_geocellpath "<<fCellDefFilePath<<"
444 for (const auto & entry : fVoxelDefFilesLi
445 fout<<"_geovolxelpath "<<std::string(e
446 }
447 fout<<"_numberOfBasepairs "<<fTotalNbBpPla
448 fout<<"_numberOfHistones "<<fTotalNbHiston
449 fout<<"_nucleusVolume "<<fNucleusVolume<<"
450 fout<<"_nucleusMassDensity "<<fNucleusMass
451 fout<<"_nucleusMass "<<fNucleusVolume*fNuc
452 fout.close();
453 }
454
455 //....oooOO0OOooo........oooOO0OOooo........oo
456
457 void Analysis::DefineCommands()
458 {
459 fMessenger = std::make_unique<G4GenericMes
460 "/dsbandrepair/ou
461 "cmd controld");
462 auto & fChemOutFolderCmd = fMessenger->Dec
463
464 fChemOutFolderCmd.SetParameterName("outFol
465 fChemOutFolderCmd.SetDefaultValue("chem_ou
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oo
469 //....oooOO0OOooo........oooOO0OOooo........oo
470
471 void Analysis::CheckAndCreateNewFolderInChemSt
472 {
473 G4fs::file_status myFs = G4fs::file_status
474 auto folderPath = G4fs::path{fChemOutFolde
475 auto isExist = G4fs::status_known(myFs) ?
476 if (isExist && !G4fs::is_empty(folderPath)
477 G4ExceptionDescription msg;
478 msg <<"==>> Chem output folder "<<fChe
479 <<" is already existing and not em
480 msg<<"==>> Please delete or rename it,
481 <<"Chem output folder in macro file!!!
482 G4Exception("Analysis::CheckAndCreateN
483 }
484 G4fs::create_directory(folderPath);
485 }
486
487 //....oooOO0OOooo........oooOO0OOooo........oo
488
489 void Analysis::CheckAndCreateNewFolderInPhysSt
490 {
491 G4fs::file_status myFs = G4fs::file_status
492 const G4fs::path folderPath{fPhysOutFolder
493 auto isExist = G4fs::status_known(myFs) ?
494 if (isExist && !G4fs::is_empty(folderPath)
495 G4fs::remove_all(folderPath);
496 }
497 G4fs::create_directory(folderPath);
498
499 const G4fs::path folderPathPC{fChemInputFo
500 isExist = G4fs::status_known(myFs) ? G4fs:
501 if (isExist && !G4fs::is_empty(folderPathP
502 G4fs::remove_all(folderPathPC);
503 }
504 G4fs::create_directory(folderPathPC);
505 }
506
507 //....oooOO0OOooo........oooOO0OOooo........oo