Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 /// \file HistoManager.cc << 27 /// \brief Implementation of the HistoManager << 28 // 26 // 29 // << 27 // $Id: HistoManager.cc,v 1.1 2010-09-16 16:26:13 gcosmo Exp $ 30 //....oooOO0OOooo........oooOO0OOooo........oo << 28 // GEANT4 tag $Name: geant4-09-04-patch-01 $ >> 29 // 31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 32 33 #include "HistoManager.hh" 33 #include "HistoManager.hh" 34 << 34 #include "HistoMessenger.hh" 35 #include "G4UnitsTable.hh" 35 #include "G4UnitsTable.hh" 36 36 >> 37 #ifdef G4ANALYSIS_USE >> 38 #include "AIDA/AIDA.h" >> 39 #endif >> 40 37 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 38 42 39 HistoManager::HistoManager() 43 HistoManager::HistoManager() >> 44 :af(0),tree(0),factoryOn(false) >> 45 { >> 46 #ifdef G4ANALYSIS_USE >> 47 // Creating the analysis factory >> 48 af = AIDA_createAnalysisFactory(); >> 49 if(!af) { >> 50 G4cout << " HistoManager::HistoManager() :" >> 51 << " problem creating the AIDA analysis factory." >> 52 << G4endl; >> 53 } >> 54 #endif >> 55 >> 56 fileName[0] = "rdecay1"; >> 57 fileType = "root"; >> 58 fileOption = ""; >> 59 // histograms >> 60 for (G4int k=0; k<MaxHisto; k++) { >> 61 histo[k] = 0; >> 62 exist[k] = false; >> 63 Unit[k] = 1.0; >> 64 Width[k] = 1.0; >> 65 ascii[k] = false; >> 66 } >> 67 >> 68 histoMessenger = new HistoMessenger(this); >> 69 } >> 70 >> 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 72 >> 73 HistoManager::~HistoManager() >> 74 { >> 75 delete histoMessenger; >> 76 >> 77 #ifdef G4ANALYSIS_USE >> 78 delete af; >> 79 #endif >> 80 } >> 81 >> 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 83 >> 84 void HistoManager::book() >> 85 { >> 86 #ifdef G4ANALYSIS_USE >> 87 if(!af) return; >> 88 >> 89 // Creating a tree mapped to an hbook file. >> 90 fileName[1] = fileName[0] + "." + fileType; >> 91 G4bool readOnly = false; >> 92 G4bool createNew = true; >> 93 AIDA::ITreeFactory* tf = af->createTreeFactory(); >> 94 tree = tf->create(fileName[1], fileType, readOnly, createNew, fileOption); >> 95 delete tf; >> 96 if(!tree) { >> 97 G4cout << "HistoManager::book() :" >> 98 << " problem creating the AIDA tree with " >> 99 << " storeName = " << fileName[1] >> 100 << " storeType = " << fileType >> 101 << " readOnly = " << readOnly >> 102 << " createNew = " << createNew >> 103 << " options = " << fileOption >> 104 << G4endl; >> 105 return; >> 106 } >> 107 >> 108 // Creating a histogram factory, whose histograms will be handled by the tree >> 109 AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree); >> 110 >> 111 // create selected histograms >> 112 for (G4int k=0; k<MaxHisto; k++) { >> 113 if (exist[k]) { >> 114 histo[k] = hf->createHistogram1D( Label[k], Title[k], >> 115 Nbins[k], Vmin[k], Vmax[k]); >> 116 factoryOn = true; >> 117 } >> 118 } >> 119 delete hf; >> 120 if (factoryOn) >> 121 G4cout << "\n----> Histogram Tree is opened in " << fileName[1] << G4endl; >> 122 #endif >> 123 } >> 124 >> 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 126 >> 127 void HistoManager::save() 40 { 128 { 41 Book(); << 129 #ifdef G4ANALYSIS_USE >> 130 if (factoryOn) { >> 131 saveAscii(); // Write ascii file, if any >> 132 tree->commit(); // Writing the histograms to the file >> 133 tree->close(); // and closing the tree (and the file) >> 134 G4cout << "\n----> Histogram Tree is saved in " << fileName[1] << G4endl; >> 135 >> 136 delete tree; >> 137 tree = 0; >> 138 factoryOn = false; >> 139 } >> 140 #endif 42 } 141 } 43 142 44 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 144 46 void HistoManager::Book() << 145 void HistoManager::FillHisto(G4int ih, G4double e, G4double weight) 47 { 146 { 48 // Create or get analysis manager << 147 if (ih > MaxHisto) { 49 G4AnalysisManager* analysisManager = G4Analy << 148 G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih 50 analysisManager->SetDefaultFileType("root"); << 149 << "does not exist; e= " << e << " w= " << weight << G4endl; 51 analysisManager->SetFileName(fFileName); << 150 return; 52 analysisManager->SetVerboseLevel(1); << 53 analysisManager->SetActivation(true); // en << 54 << 55 // Define histograms start values << 56 const G4int kMaxHisto = 10; << 57 const G4String id[] = {"0", "1", "2", "3", " << 58 const G4String title[] = { << 59 "dummy", // 0 << 60 "energy spectrum (%): e+ e-", // 1 << 61 "energy spectrum (%): nu_e anti_nu_e", // << 62 "energy spectrum (%): gamma", // 3 << 63 "energy spectrum (%): alpha", // 4 << 64 "energy spectrum (%): ions", // 5 << 65 "total kinetic energy per single decay (Q) << 66 "momentum balance", // 7 << 67 "total time of life of decay chain", // 8 << 68 "total visible energy in decay chain" // << 69 }; << 70 << 71 // Default values (to be reset via /analysis << 72 G4int nbins = 100; << 73 G4double vmin = 0.; << 74 G4double vmax = 100.; << 75 << 76 // Create all histograms as inactivated << 77 // as we have not yet set nbins, vmin, vmax << 78 for (G4int k = 0; k < kMaxHisto; k++) { << 79 G4int ih = analysisManager->CreateH1(id[k] << 80 analysisManager->SetH1Activation(ih, false << 81 } 151 } >> 152 #ifdef G4ANALYSIS_USE >> 153 if(exist[ih]) histo[ih]->fill(e/Unit[ih], weight); >> 154 #endif 82 } 155 } 83 156 84 //....oooOO0OOooo........oooOO0OOooo........oo 157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 158 >> 159 void HistoManager::SetHisto(G4int ih, >> 160 G4int nbins, G4double valmin, G4double valmax, const G4String& unit) >> 161 { >> 162 if (ih > MaxHisto) { >> 163 G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih >> 164 << "does not exist" << G4endl; >> 165 return; >> 166 } >> 167 >> 168 const G4String id[] = { "0", "1", "2", "3", "4", "5", "6", "7" , "8" , "9" }; >> 169 const G4String title[] = >> 170 { "dummy", //0 >> 171 "energy spectrum: e+ e-", //1 >> 172 "energy spectrum: nu_e anti_nu_e", //2 >> 173 "energy spectrum: gamma", //3 >> 174 "energy spectrum: alpha", //4 >> 175 "energy spectrum: ions", //5 >> 176 "energy balance", //6 >> 177 "momentum balance", //7 >> 178 "total time of life of decay chain",//8 >> 179 "activity in Bq/g" //9 >> 180 }; >> 181 >> 182 >> 183 G4String titl = title[ih]; >> 184 G4double vmin = valmin, vmax = valmax; >> 185 Unit[ih] = 1.; >> 186 >> 187 if (unit != "none") { >> 188 titl = title[ih] + " (" + unit + ")"; >> 189 Unit[ih] = G4UnitDefinition::GetValueOf(unit); >> 190 vmin = valmin/Unit[ih]; vmax = valmax/Unit[ih]; >> 191 } >> 192 >> 193 exist[ih] = true; >> 194 Label[ih] = id[ih]; >> 195 Title[ih] = titl; >> 196 Nbins[ih] = nbins; >> 197 Vmin[ih] = vmin; >> 198 Vmax[ih] = vmax; >> 199 Width[ih] = (valmax-valmin)/nbins; >> 200 >> 201 G4cout << "----> SetHisto " << ih << ": " << titl << "; " >> 202 << nbins << " bins from " >> 203 << vmin << " " << unit << " to " << vmax << " " << unit << G4endl; >> 204 >> 205 } >> 206 >> 207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 208 >> 209 void HistoManager::Normalize(G4int ih, G4double fac) >> 210 { >> 211 if (ih >= MaxHisto) { >> 212 G4cout << "---> warning from HistoManager::Normalize() : histo " << ih >> 213 << " fac= " << fac << G4endl; >> 214 return; >> 215 } >> 216 #ifdef G4ANALYSIS_USE >> 217 if(exist[ih]) histo[ih]->scale(fac); >> 218 #endif >> 219 } >> 220 >> 221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 222 >> 223 void HistoManager::RemoveHisto(G4int ih) >> 224 { >> 225 if (ih > MaxHisto) { >> 226 G4cout << "---> warning from HistoManager::RemoveHisto() : histo " << ih >> 227 << "does not exist" << G4endl; >> 228 return; >> 229 } >> 230 >> 231 histo[ih] = 0; exist[ih] = false; >> 232 } >> 233 >> 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 235 >> 236 void HistoManager::PrintHisto(G4int ih) >> 237 { >> 238 if (ih < MaxHisto) { ascii[ih] = true; ascii[0] = true; } >> 239 else >> 240 G4cout << "---> warning from HistoManager::PrintHisto() : histo " << ih >> 241 << "does not exist" << G4endl; >> 242 } >> 243 >> 244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 245 >> 246 #include <fstream> >> 247 >> 248 void HistoManager::saveAscii() >> 249 { >> 250 #ifdef G4ANALYSIS_USE >> 251 >> 252 if (!ascii[0]) return; >> 253 >> 254 G4String name = fileName[0] + ".ascii"; >> 255 std::ofstream File(name, std::ios::out); >> 256 File.setf( std::ios::scientific, std::ios::floatfield ); >> 257 >> 258 //write selected histograms >> 259 for (G4int ih=0; ih<MaxHisto; ih++) { >> 260 if (exist[ih] && ascii[ih]) { >> 261 File << "\n 1D histogram " << ih << ": " << Title[ih] >> 262 << "\n \n \t X \t\t Y" << G4endl; >> 263 >> 264 for (G4int iBin=0; iBin<Nbins[ih]; iBin++) { >> 265 File << " " << iBin << "\t" >> 266 << 0.5*(histo[ih]->axis().binLowerEdge(iBin) + >> 267 histo[ih]->axis().binUpperEdge(iBin)) << "\t" >> 268 << histo[ih]->binHeight(iBin) >> 269 << G4endl; >> 270 } >> 271 } >> 272 } >> 273 #endif >> 274 } >> 275 >> 276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 277 85 278