| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file Analysis.cc
28 /// \brief Implementation of the Analysis class
29 /// \file Analysis.cc
30 /// \brief Implementation of the Analysis class
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........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
60
61 void Analysis::OpenFile(const G4String outFolder)
62 {
63 G4AnalysisManager* anManager = G4AnalysisManager::Instance();
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+fFileName;
75 }
76
77 // open output file
78 G4bool fileOpen = anManager->OpenFile(fullFileName.c_str());
79 if (!fileOpen) {
80 G4cout << "\n---> HistoManager::book(): cannot open " << fFileName
81 << G4endl;
82 return;
83 }
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87
88 void Analysis::Save()
89 {
90 G4AnalysisManager* anManager = G4AnalysisManager::Instance();
91 anManager->Write();
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95
96 void Analysis::Close(G4bool reset)
97 {
98 G4AnalysisManager* anManager = G4AnalysisManager::Instance();
99 anManager->CloseFile(reset);
100 }
101
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103
104 void Analysis::Book()
105 {
106 G4AnalysisManager* anManager = G4AnalysisManager::Instance();
107 anManager->SetVerboseLevel(1);
108 if (anManager->GetFirstNtupleId() != 1) anManager->SetFirstNtupleId(1);
109
110 if (gRunMode == RunningMode::Phys) {
111 #ifdef G4MULTITHREADED
112 // MT ntuple merging
113 anManager->SetNtupleMerging(true);//for future use MT+MPI
114 #endif
115 anManager->SetNtupleDirectoryName("ntuple");
116 // create ntuple
117 anManager->CreateNtuple("ntuple_1","physical_stage");
118 anManager->CreateNtupleIColumn("flagParticle");
119 anManager->CreateNtupleIColumn("flagParentID");
120 anManager->CreateNtupleIColumn("flagProcess");
121 anManager->CreateNtupleDColumn("x");
122 anManager->CreateNtupleDColumn("y");
123 anManager->CreateNtupleDColumn("z");
124 anManager->CreateNtupleDColumn("edep");
125 anManager->CreateNtupleIColumn("eventNumber");
126 anManager->CreateNtupleIColumn("volumeName");
127 anManager->CreateNtupleIColumn("copyNumber");
128 anManager->CreateNtupleIColumn("lastMetVoxelCopyNum");
129 anManager->FinishNtuple(1);
130
131 // For total edep
132 anManager->CreateNtuple("ntuple_3","total_edep");
133 anManager->CreateNtupleIColumn("eventNumber");
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() != directoryName) {
142 anManager->SetNtupleDirectoryName(directoryName);
143 }
144
145 // DBScan
146
147 anManager->CreateNtuple("ntuple_2","DB_chemical_stage");
148 anManager->CreateNtupleIColumn(1,"strand");
149 anManager->CreateNtupleIColumn(1,"copyNumber");
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........oooOO0OOooo........oooOO0OOooo......
160
161 G4AnalysisManager* Analysis::GetAnalysisManager()
162 {
163 return G4AnalysisManager::Instance();
164 }
165
166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167
168 void Analysis::AddInfoForChemGeo(InfoForChemGeo b)
169 {
170 fInfoForChemGeoVector.push_back(b);
171 }
172
173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174
175 void Analysis::ClearVector()
176 {
177 fInfoForChemGeoVector.clear();
178 fInfoInPhysStageVector.clear();
179 fOutputFiles.clear();
180 }
181
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183
184 void Analysis::AddInfoInPhysStage(InfoInPhysStage b)
185 {
186 fInfoInPhysStageVector.push_back(b);
187 }
188
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190
191 void Analysis::UpdateChemInputDataAndFillNtuple()
192 {
193 //if (fInfoForChemGeoVector.size() == 0) return;
194
195 // We will loop on all the "element" of the vector and create several files:
196 // - New file for each event/voxel couple
197 // Create a map to define all the output file
198
199 for (auto const& binfo : fInfoForChemGeoVector) {
200 if( binfo.fVolume == 161 // voxelStraight
201 || binfo.fVolume == 162 // voxelRight
202 || binfo.fVolume == 163 // voxelLeft
203 || binfo.fVolume == 164 // voxelUp
204 || binfo.fVolume == 165 // voxelDown
205 || binfo.fVolume == 261 // voxelStraight2
206 || binfo.fVolume == 262 // voxelRight2
207 || binfo.fVolume == 263 // voxelLeft2
208 || binfo.fVolume == 264 // voxelUp2
209 || binfo.fVolume == 265) // voxelDown2
210 {
211 // We are in a voxel
212
213 std::string voxelName("");
214
215 if(binfo.fVolume==161) voxelName = "VoxelStraight";
216 else if(binfo.fVolume==162) voxelName = "VoxelRight";
217 else if(binfo.fVolume==163) voxelName = "VoxelLeft";
218 else if(binfo.fVolume==164) voxelName = "VoxelUp";
219 else if(binfo.fVolume==165) voxelName = "VoxelDown";
220 else if(binfo.fVolume==261) voxelName = "VoxelStraight2";
221 else if(binfo.fVolume==262) voxelName = "VoxelRight2";
222 else if(binfo.fVolume==263) voxelName = "VoxelLeft2";
223 else if(binfo.fVolume==264) voxelName = "VoxelUp2";
224 else if(binfo.fVolume==265) voxelName = "VoxelDown2";
225
226 // Get the event number
227 int eventNum = int(binfo.fEventNumber);
228 // Here we have all the information for one ntuple line
229 // We want to know if we have to create a new output file (new eventNumber/voxel couple)
230 // or if we already have one.
231
232 // Check if the event has already been registered
233 if(fOutputFiles.find(eventNum)==fOutputFiles.end() )
234 {
235 // If not then create the event and voxel case
236 fOutputFiles[eventNum][binfo.fVolumeCopyNumber] =
237 CreateChemInputFile(eventNum, int(binfo.fVolumeCopyNumber), voxelName);
238 }
239 // If the event has been registered then we need to check that the current voxel has also been
240 // registered.
241 else
242 {
243 if(fOutputFiles[eventNum].find(binfo.fVolumeCopyNumber)==fOutputFiles[eventNum].end())
244 {
245 // We register the event-voxel couple
246 fOutputFiles[eventNum][binfo.fVolumeCopyNumber] =
247 CreateChemInputFile(eventNum, int(binfo.fVolumeCopyNumber), voxelName);
248 }
249 }
250
251 // Write in the file the information needed by the chemistry simulation to build
252 // the input water molecule or solvated electron.
253 UpdatingChemInputFile(binfo);
254 }
255 }
256
257 if (fInfoInPhysStageVector.size() == 0) return;
258 for (auto const& binfo : fInfoInPhysStageVector) {
259 // Check if the process is an ionisation
260 // Only ionisation should trigger the removal of a DNA molecule from the chemical step
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 happened in a dna molecule or its hydration shell
269 if( binfo.fVolumeName == 1 // d1
270 || binfo.fVolumeName == 11 // p1
271 || binfo.fVolumeName == 2 // d2
272 || binfo.fVolumeName == 22 // p2
273 || binfo.fVolumeName == 3 // cyto
274 || binfo.fVolumeName == 4 // gua
275 || binfo.fVolumeName == 5 // thy
276 || binfo.fVolumeName == 6 // ade
277 || binfo.fVolumeName == 7 // d1_w
278 || binfo.fVolumeName == 71 // p1_w
279 || binfo.fVolumeName == 8 // d2_w
280 || binfo.fVolumeName == 81 // p2_w
281 || binfo.fVolumeName == 9 // ade_w
282 || binfo.fVolumeName == 10 // gua_w
283 || binfo.fVolumeName == 13 // cyto_w
284 || binfo.fVolumeName == 12) // thy_w
285 {
286 // Retrieve the voxel copy number
287 double voxelCopyNumber = binfo.fLastMetVoxelCopyNum;
288 // Get the event number
289 double eventNum = binfo.fEventNumber;
290
291 // Check if the event has already been registered
292 if(fOutputFiles.find(eventNum) != fOutputFiles.end() )
293 {
294 // Check if the volume number has been registered
295 if(fOutputFiles.at(eventNum).find(voxelCopyNumber)
296 != fOutputFiles.at(eventNum).end() )
297 {
298 // If we are here then the event and volume couple has a already generated
299 //file in which we should add a dna molecule to be removed
300 UpdatingChemInputFile(binfo);
301 }
302 }
303 }
304 }
305
306 // fill ntuple
307 auto analysisManager =G4AnalysisManager::Instance();
308 analysisManager->FillNtupleIColumn(1, 0, binfo.fFlagParticle);
309 analysisManager->FillNtupleIColumn(1, 1, binfo.fFlagParentID);
310 analysisManager->FillNtupleIColumn(1, 2, binfo.fFlagProcess);
311 analysisManager->FillNtupleDColumn(1, 3, binfo.fX);
312 analysisManager->FillNtupleDColumn(1, 4, binfo.fY);
313 analysisManager->FillNtupleDColumn(1, 5, binfo.fZ);
314 analysisManager->FillNtupleDColumn(1, 6, binfo.fEdep);
315 analysisManager->FillNtupleIColumn(1, 7, binfo.fEventNumber);
316 analysisManager->FillNtupleIColumn(1, 8, binfo.fVolumeName);
317 analysisManager->FillNtupleIColumn(1, 9, binfo.fCopyNumber);
318 analysisManager->FillNtupleIColumn(1, 10, binfo.fLastMetVoxelCopyNum);
319 analysisManager->AddNtupleRow(1);
320 }
321
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
325
326 G4String Analysis::CreateChemInputFile(G4int eventNum, G4int volumeCopyNumber, const G4String &voxelName)
327 {
328 std::stringstream sstream;
329 sstream<<"./"<<fChemInputFolderName<<"/event_"<<eventNum<<"_voxel_"<<volumeCopyNumber<<".dat";
330 std::ofstream oFile;
331 oFile.open(sstream.str().c_str() );
332
333 oFile<<"_eventNum"<<"\t\t"<<eventNum<<std::endl;
334 oFile<<"_voxelType"<<"\t\t"<<voxelName<<std::endl;
335 oFile<<"_voxelCopyNumber"<<"\t\t"<<volumeCopyNumber<<std::endl;
336 oFile<<std::endl;
337
338 oFile
339 <<"# Chemistry input informations\n"
340 <<"# "<<"_input, type, state, electronicLevel, x, y, z, parentTrackID"<<std::endl;
341 oFile<<"# type=1 -> water molecule && type=2 -> solvated electron"<<std::endl;
342 oFile<<std::endl;
343
344 oFile.close();
345
346 return sstream.str().c_str();
347 }
348
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
350
351 void Analysis::UpdatingChemInputFile(InfoForChemGeo b)
352 {
353 std::ofstream out;
354 out.open( fOutputFiles[b.fEventNumber][b.fVolumeCopyNumber].c_str(), std::ios::app); // Open the file and go at the end
355
356 if (b.fType==1)
357 {
358 out<<"_input"<<"\t"
359 <<b.fType<<"\t\t"
360 <<b.fState<<"\t\t"
361 <<4-b.fElectronicLevel<<"\t\t"
362 <<b.fRelX<<"\t\t"
363 <<b.fRelY<<"\t\t"
364 <<b.fRelZ<<"\t\t"
365 <<b.fParentTrackID<<"\t\t"
366 <<"\n";
367 }
368 if (b.fType==2)
369 {
370 out<<"_input"<<"\t"
371 <<b.fType<<"\t\t"
372 <<b.fState<<"\t\t"
373 <<b.fElectronicLevel<<"\t\t"
374 // If we are here then the event and volume couple has a already
375 //generated file in which we should add a dna molecule to be removed
376 <<b.fRelX<<"\t\t"
377 <<b.fRelY<<"\t\t"
378 <<b.fRelZ<<"\t\t"
379 <<b.fParentTrackID<<"\t\t"
380 <<"\n";
381 }
382 out.close();
383 }
384
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
386
387 void Analysis::UpdatingChemInputFile(InfoInPhysStage b)
388 {
389 G4String name;
390 G4double volumeName = b.fVolumeName;
391 if(volumeName == 1 || volumeName == 2 || volumeName == 7 ||
392 volumeName == 8) name = G4Deoxyribose::Definition()->GetName();
393 else if(volumeName == 11 || volumeName == 22 ||
394 volumeName == 71 || volumeName == 81) name = G4Phosphate::Definition()->GetName();
395 else if(volumeName == 6 || volumeName == 9) name = G4Adenine::Definition()->GetName();
396 else if(volumeName == 4 || volumeName == 10) name = G4Guanine::Definition()->GetName();
397 else if(volumeName == 5 || volumeName == 12) name = G4Thymine::Definition()->GetName();
398 else if(volumeName == 3 || volumeName == 13) name = G4Cytosine::Definition()->GetName();
399 else
400 {
401 G4ExceptionDescription msg;
402 msg <<"Volume number "<<volumeName<<" not registered.";
403 G4Exception("Analysis::UpdatingChemInputFile",
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.fLastMetVoxelCopyNum].c_str(), std::ios::app); // Open the file and go at the end
430 out<<"_remove"<<"\t"<<name<<"\t\t"<<b.fCopyNumber<<"\t\t"<<strand<<"\t\t"<<std::endl;
431 out.close();
432 }
433
434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
435
436 void Analysis::WritePhysGeo()
437 {
438 //Create a file containing geopath and Physlist in "FS build directory" for chemStage and anlysis
439 std::string fname = "imp.info";
440 std::ofstream fout(fname);
441 fout<<"====> Auto_generated file. Do not delete me !!!\n";
442 fout<<"====> The file conveys some information for chem_geo and Analysis modules!!!\n";
443 fout<<"_geocellpath "<<fCellDefFilePath<<"\n";
444 for (const auto & entry : fVoxelDefFilesList) {
445 fout<<"_geovolxelpath "<<std::string(entry)<<"\n";
446 }
447 fout<<"_numberOfBasepairs "<<fTotalNbBpPlacedInGeo<<"\n";
448 fout<<"_numberOfHistones "<<fTotalNbHistonePlacedInGeo<<"\n";
449 fout<<"_nucleusVolume "<<fNucleusVolume<<"\n"; // m3
450 fout<<"_nucleusMassDensity "<<fNucleusMassDensity<<"\n"; // kg/m3
451 fout<<"_nucleusMass "<<fNucleusVolume*fNucleusMassDensity; //kg
452 fout.close();
453 }
454
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456
457 void Analysis::DefineCommands()
458 {
459 fMessenger = std::make_unique<G4GenericMessenger>(this,
460 "/dsbandrepair/output/",
461 "cmd controld");
462 auto & fChemOutFolderCmd = fMessenger->DeclareProperty ("folderForChemOut",
463 fChemOutFolderName);
464 fChemOutFolderCmd.SetParameterName("outFolderForChem",true);
465 fChemOutFolderCmd.SetDefaultValue("chem_output");
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
470
471 void Analysis::CheckAndCreateNewFolderInChemStage()
472 {
473 G4fs::file_status myFs = G4fs::file_status{};
474 auto folderPath = G4fs::path{fChemOutFolderName.c_str()};
475 auto isExist = G4fs::status_known(myFs) ? G4fs::exists(myFs) : G4fs::exists(folderPath);
476 if (isExist && !G4fs::is_empty(folderPath)) {
477 G4ExceptionDescription msg;
478 msg <<"==>> Chem output folder "<<fChemOutFolderName
479 <<" is already existing and not empty!!!\n";
480 msg<<"==>> Please delete or rename it, or use other name for "
481 <<"Chem output folder in macro file!!!\n";
482 G4Exception("Analysis::CheckAndCreateNewFolderInChemStage()","",FatalException,msg);
483 }
484 G4fs::create_directory(folderPath);
485 }
486
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488
489 void Analysis::CheckAndCreateNewFolderInPhysStage()
490 {
491 G4fs::file_status myFs = G4fs::file_status{};
492 const G4fs::path folderPath{fPhysOutFolderName.c_str()};
493 auto isExist = G4fs::status_known(myFs) ? G4fs::exists(myFs) : G4fs::exists(folderPath);
494 if (isExist && !G4fs::is_empty(folderPath)) {
495 G4fs::remove_all(folderPath);
496 }
497 G4fs::create_directory(folderPath);
498
499 const G4fs::path folderPathPC{fChemInputFolderName.c_str()};
500 isExist = G4fs::status_known(myFs) ? G4fs::exists(myFs) : G4fs::exists(folderPathPC);
501 if (isExist && !G4fs::is_empty(folderPathPC)) {
502 G4fs::remove_all(folderPathPC);
503 }
504 G4fs::create_directory(folderPathPC);
505 }
506
507 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....