Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // << 23 // $Id: HistoManager.cc,v 1.10 2005/06/01 10:48:40 maire Exp $ >> 24 // GEANT4 tag $Name: geant4-07-01-patch-01 $ >> 25 // 27 //....oooOO0OOooo........oooOO0OOooo........oo 26 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 28 //....oooOO0OOooo........oooOO0OOooo........oo 27 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 29 28 30 #include "HistoManager.hh" 29 #include "HistoManager.hh" 31 << 30 #include "HistoMessenger.hh" 32 #include "DetectorConstruction.hh" << 33 << 34 #include "G4UnitsTable.hh" 31 #include "G4UnitsTable.hh" 35 32 >> 33 #ifdef G4ANALYSIS_USE >> 34 #include "AIDA/AIDA.h" >> 35 #endif >> 36 36 //....oooOO0OOooo........oooOO0OOooo........oo 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 37 38 38 HistoManager::HistoManager() 39 HistoManager::HistoManager() >> 40 :af(0),tree(0),factoryOn(false) 39 { 41 { 40 Book(); << 42 #ifdef G4ANALYSIS_USE >> 43 // Creating the analysis factory >> 44 af = AIDA_createAnalysisFactory(); >> 45 if(!af) { >> 46 G4cout << "TestEm1::HistoManager::HistoManager :" >> 47 << " problem creating the AIDA analysis factory." >> 48 << G4endl; >> 49 } >> 50 #endif >> 51 >> 52 fileName[0] = "testem3"; >> 53 fileType = "hbook"; >> 54 fileOption = "--noErrors uncompress"; >> 55 // histograms >> 56 for (G4int k=0; k<MaxHisto; k++) { >> 57 histo[k] = 0; >> 58 exist[k] = false; >> 59 Unit[k] = 1.0; >> 60 Width[k] = 1.0; >> 61 } >> 62 histoMessenger = new HistoMessenger(this); >> 63 } >> 64 >> 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 66 >> 67 HistoManager::~HistoManager() >> 68 { >> 69 delete histoMessenger; >> 70 >> 71 #ifdef G4ANALYSIS_USE >> 72 delete af; >> 73 #endif 41 } 74 } 42 75 43 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 44 77 45 void HistoManager::Book() << 78 void HistoManager::book() >> 79 { >> 80 #ifdef G4ANALYSIS_USE >> 81 if(!af) return; >> 82 >> 83 // Creating a tree mapped to an hbook file. >> 84 fileName[1] = fileName[0] + "." + fileType; >> 85 G4bool readOnly = false; >> 86 G4bool createNew = true; >> 87 AIDA::ITreeFactory* tf = af->createTreeFactory(); >> 88 tree = tf->create(fileName[1], fileType, readOnly, createNew, fileOption); >> 89 delete tf; >> 90 if(!tree) { >> 91 G4cout << "TestEm1::HistoManager::book :" >> 92 << " problem creating the AIDA tree with " >> 93 << " storeName = " << fileName[1] >> 94 << " storeType = " << fileType >> 95 << " readOnly = " << readOnly >> 96 << " createNew = " << createNew >> 97 << " options = " << fileOption >> 98 << G4endl; >> 99 return; >> 100 } >> 101 >> 102 // Creating a histogram factory, whose histograms will be handled by the tree >> 103 AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree); >> 104 >> 105 // create selected histograms >> 106 for (G4int k=1; k<MaxHisto; k++) { >> 107 if (exist[k]) { >> 108 histo[k] = hf->createHistogram1D( Label[k], Title[k], >> 109 Nbins[k], Vmin[k], Vmax[k]); >> 110 factoryOn = true; >> 111 } >> 112 } >> 113 delete hf; >> 114 if (factoryOn) >> 115 G4cout << "\n----> Histogram Tree is opened in " << fileName[1] << G4endl; >> 116 #endif >> 117 } >> 118 >> 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 120 >> 121 void HistoManager::save() >> 122 { >> 123 #ifdef G4ANALYSIS_USE >> 124 if (factoryOn) { >> 125 tree->commit(); // Writing the histograms to the file >> 126 tree->close(); // and closing the tree (and the file) >> 127 G4cout << "\n----> Histogram Tree is saved in " << fileName[1] << G4endl; >> 128 >> 129 delete tree; >> 130 tree = 0; >> 131 factoryOn = false; >> 132 } >> 133 #endif >> 134 } >> 135 >> 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 137 >> 138 void HistoManager::FillHisto(G4int ih, G4double xbin, G4double weight) 46 { 139 { 47 // Create or get analysis manager << 140 if (ih >= MaxHisto) { 48 // The choice of analysis technology is done << 141 G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih 49 // in HistoManager.hh << 142 << G4endl; 50 G4AnalysisManager* analysisManager = G4Analy << 143 return; 51 analysisManager->SetDefaultFileType("root"); << 144 } 52 analysisManager->SetFileName(fFileName); << 145 #ifdef G4ANALYSIS_USE 53 analysisManager->SetVerboseLevel(1); << 146 if(exist[ih]) histo[ih]->fill(xbin/Unit[ih], weight); 54 analysisManager->SetActivation(true); // en << 147 #endif 55 << 148 } 56 // Define histograms start values << 149 57 << 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 58 const G4String id[] = {"0", "1", "2", "3" << 151 59 "9", "10", "11", "12 << 152 void HistoManager::SetHisto(G4int ih, 60 "18", "19", "20", "21 << 153 G4int nbins, G4double valmin, G4double valmax, const G4String& unit) >> 154 { >> 155 if (ih < 1 || ih >= MaxHisto) { >> 156 G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih >> 157 << "does not exist" << G4endl; >> 158 return; >> 159 } >> 160 >> 161 // histo 1 : energy deposit in absorber 1 >> 162 // histo 2 : energy deposit in absorber 2 >> 163 // ...etc........... >> 164 // MaxAbsor = 10 (-1) >> 165 // >> 166 // histo 11 : longitudinal profile of energy deposit in absorber 1 (MeV) >> 167 // histo 12 : longitudinal profile of energy deposit in absorber 2 (MeV) >> 168 // ...etc........... >> 169 // >> 170 // histo 21 : forward energy flow (MeV) >> 171 // histo 22 : backward energy flow (MeV) >> 172 >> 173 const G4String id[] = { "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", >> 174 "10","11","12","13","14","15","16","17","18","19", >> 175 "20","21","22"}; >> 176 61 G4String title; 177 G4String title; >> 178 G4double vmin = valmin, vmax = valmax; >> 179 >> 180 if (ih < MaxAbsor) { >> 181 title = "Edep in absorber " + id[ih] + " (" + unit + ")"; >> 182 Unit[ih] = G4UnitDefinition::GetValueOf(unit); >> 183 vmin = valmin/Unit[ih]; vmax = valmax/Unit[ih]; >> 184 } else if (ih > MaxAbsor && ih < 2*MaxAbsor) { >> 185 title = "longit. profile of Edep (MeV/event) in absorber " >> 186 + id[ih-MaxAbsor]; >> 187 } else if (ih == 2*MaxAbsor+1) { >> 188 title = "Forward energy flow (MeV/event)"; >> 189 } else if (ih == 2*MaxAbsor+2) { >> 190 title = "Backward energy flow (MeV/event)"; >> 191 } else return; >> 192 >> 193 exist[ih] = true; >> 194 Label[ih] = id[ih]; >> 195 Title[ih] = title; >> 196 Nbins[ih] = nbins; >> 197 Vmin[ih] = vmin; >> 198 Vmax[ih] = vmax; >> 199 Width[ih] = (valmax-valmin)/nbins; >> 200 >> 201 G4cout << "----> SetHisto " << ih << ": " << title << "; " >> 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 << 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...... 62 222 63 // Default values (to be reset via /analysis << 223 void HistoManager::RemoveHisto(G4int ih) 64 G4int nbins = 100; << 224 { 65 G4double vmin = 0.; << 225 if (ih >= MaxHisto) { 66 G4double vmax = 100.; << 226 G4cout << "---> warning from HistoManager::RemoveHisto() : histo " << ih 67 << 227 << "does not exist" << G4endl; 68 // Create all histograms as inactivated << 228 return; 69 // as we have not yet set nbins, vmin, vmax << 70 for (G4int k = 0; k < kMaxHisto; k++) { << 71 if (k < kMaxAbsor) title = "Edep in absorb << 72 if (k == kMaxAbsor) title = "total energy << 73 if (k > kMaxAbsor) title = "Edep longit. p << 74 if (k == 2 * kMaxAbsor + 1) title = "energ << 75 if (k == 2 * kMaxAbsor + 2) title = "later << 76 if (k == 2 * kMaxAbsor + 3) title = "total << 77 if (k == 2 * kMaxAbsor + 4) title = "total << 78 G4int ih = analysisManager->CreateH1(id[k] << 79 analysisManager->SetH1Activation(ih, false << 80 } 229 } >> 230 >> 231 histo[ih] = 0; exist[ih] = false; 81 } 232 } >> 233 >> 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 235 >> 236 82 237