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 Define histograms and ntuples 29 30 #include "Analysis.hh" 31 32 #include "G4AutoDelete.hh" 33 #include "G4SystemOfUnits.hh" 34 35 // Note: ntuple merging is supported only with Root format 36 #ifndef G4MULTITHREADED 37 # include "G4RootMpiAnalysisManager.hh" 38 using G4AnalysisManager = G4RootMpiAnalysisManager; 39 #else 40 # include "G4RootAnalysisManager.hh" 41 using G4AnalysisManager = G4RootAnalysisManager; 42 #endif 43 44 G4ThreadLocal G4int Analysis::fincidentFlag = false; 45 G4ThreadLocal Analysis* the_analysis = 0; 46 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 48 Analysis* Analysis::GetAnalysis() 49 { 50 if (!the_analysis) { 51 the_analysis = new Analysis(); 52 G4AutoDelete::Register(the_analysis); 53 } 54 return the_analysis; 55 } 56 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 58 Analysis::Analysis() 59 : fUseNtuple(true), 60 fMergeNtuple(true), 61 fincident_x_hist(0), 62 fincident_map(0), 63 fdose_hist(0), 64 fdose_map(0), 65 fdose_prof(0), 66 fdose_map_prof(0), 67 fdose_map3d(0) 68 {} 69 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 71 void Analysis::Book() 72 { 73 G4cout << "Analysis::Book start, fUseNtuple: " << fUseNtuple << G4endl; 74 75 G4AnalysisManager* mgr = G4AnalysisManager::Instance(); 76 mgr->SetVerboseLevel(1); 77 #ifdef G4MULTITHREADED 78 // MT ntuple merging 79 mgr->SetNtupleMerging(fMergeNtuple); 80 #endif 81 82 fincident_x_hist = mgr->CreateH1("incident_x", "Incident X", 100, -5 * cm, 5 * cm, "cm"); 83 fincident_map = mgr->CreateH2("incident_map", "Incident Map", 50, -5 * cm, 5 * cm, 50, -5 * cm, 84 5 * cm, "cm", "cm"); 85 fdose_hist = mgr->CreateH1("dose", "Dose distribution", 500, 0, 50 * cm, "cm"); 86 fdose_map = mgr->CreateH2("dose_map", "Dose distribution", 500, 0, 50 * cm, 200, -10 * cm, 87 10 * cm, "cm", "cm"); 88 fdose_map3d = mgr->CreateH3("dose_map_3d", "Dose distribution", 30, 0, 50 * cm, 20, -10 * cm, 89 10 * cm, 20, -10 * cm, 10 * cm, "cm", "cm", "cm"); 90 fdose_prof = 91 mgr->CreateP1("dose_prof", "Dose distribution", 300, 0, 30 * cm, 0, 100 * MeV, "cm", "MeV"); 92 fdose_map_prof = mgr->CreateP2("dose_map_prof", "Dose distribution", 300, 0, 30 * cm, 80, -4 * cm, 93 4 * cm, 0, 100 * MeV, "cm", "cm", "MeV"); 94 95 if (fUseNtuple) { 96 mgr->CreateNtuple("Dose", "Dose distribution"); // ntuple Id = 0 97 mgr->CreateNtupleDColumn("pz"); 98 mgr->CreateNtupleDColumn("px"); 99 mgr->CreateNtupleDColumn("dose"); 100 mgr->FinishNtuple(); 101 } 102 103 G4cout << "Analysis::Book finished " << G4endl; 104 } 105 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 107 Analysis::~Analysis() {} 108 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 110 void Analysis::OpenFile(const G4String& fname) 111 { 112 G4AnalysisManager* mgr = G4AnalysisManager::Instance(); 113 mgr->OpenFile(fname.c_str()); 114 } 115 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 117 void Analysis::Save() 118 { 119 G4AnalysisManager* mgr = G4AnalysisManager::Instance(); 120 mgr->Write(); 121 } 122 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 124 void Analysis::Close(G4bool reset) 125 { 126 G4AnalysisManager* mgr = G4AnalysisManager::Instance(); 127 mgr->CloseFile(reset); 128 } 129 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 131 void Analysis::FillIncident(const G4ThreeVector& p) 132 { 133 if (!fincidentFlag) { 134 G4AnalysisManager* mgr = G4AnalysisManager::Instance(); 135 mgr->FillH2(fincident_map, p.x(), p.y()); 136 mgr->FillH1(fincident_x_hist, p.x()); 137 fincidentFlag = true; 138 } 139 } 140 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 142 void Analysis::FillDose(const G4ThreeVector& p, G4double dedx) 143 { 144 const G4double dxy = 10. * mm; 145 if (std::abs(p.y()) < dxy) { 146 G4AnalysisManager* mgr = G4AnalysisManager::Instance(); 147 const G4double Z0 = 25. * cm; 148 149 mgr->FillH2(fdose_map, p.z() + Z0, p.x(), dedx / GeV); 150 mgr->FillP2(fdose_map_prof, p.z() + Z0, p.x(), dedx); 151 mgr->FillH3(fdose_map3d, p.z() + Z0, p.x(), p.y(), dedx / GeV); 152 if (std::abs(p.x()) < dxy) { 153 mgr->FillH1(fdose_hist, p.z() + Z0, dedx / GeV); 154 mgr->FillP1(fdose_prof, p.z() + Z0, dedx); 155 } 156 157 if (fUseNtuple) { 158 // the same fill frequency as H2 "dose_map" 159 mgr->FillNtupleDColumn(0, p.z() + Z0); 160 mgr->FillNtupleDColumn(1, p.x()); 161 mgr->FillNtupleDColumn(2, dedx / GeV); 162 mgr->AddNtupleRow(); 163 } 164 } 165 } 166