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 electromagnetic/TestEm1/src/HistoMan << 26 // $Id: HistoManager.cc,v 1.4 2009/07/24 13:02:02 maire Exp $ 27 /// \brief Implementation of the HistoManager << 27 // GEANT4 tag $Name: geant4-09-03-patch-02 $ 28 // << 29 // 28 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 31 33 #include "HistoManager.hh" 32 #include "HistoManager.hh" >> 33 #include "HistoMessenger.hh" 34 #include "G4UnitsTable.hh" 34 #include "G4UnitsTable.hh" 35 35 >> 36 #ifdef G4ANALYSIS_USE >> 37 #include "AIDA/AIDA.h" >> 38 #endif >> 39 36 //....oooOO0OOooo........oooOO0OOooo........oo 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 37 41 38 HistoManager::HistoManager() 42 HistoManager::HistoManager() 39 : fFileName("amsEcal") << 43 :factoryOn(false),af(0),tree(0) 40 { 44 { 41 Book(); << 45 #ifdef G4ANALYSIS_USE >> 46 // Creating the analysis factory >> 47 af = AIDA_createAnalysisFactory(); >> 48 if(!af) { >> 49 G4cout << " HistoManager::HistoManager :" >> 50 << " problem creating the AIDA analysis factory." >> 51 << G4endl; >> 52 } >> 53 #endif >> 54 >> 55 fileName[0] = "ecal"; >> 56 fileType = "root"; >> 57 fileOption = "--noErrors export=root uncompress"; >> 58 >> 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 // ntuples >> 69 for (G4int k=0; k<MaxNtupl; k++) { >> 70 ntupl[k] = 0; >> 71 existNt[k] = false; >> 72 } >> 73 >> 74 histoMessenger = new HistoMessenger(this); 42 } 75 } 43 76 44 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 78 46 HistoManager::~HistoManager() 79 HistoManager::~HistoManager() >> 80 { >> 81 delete histoMessenger; >> 82 >> 83 #ifdef G4ANALYSIS_USE >> 84 delete af; >> 85 #endif >> 86 } >> 87 >> 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 89 >> 90 void HistoManager::book() >> 91 { >> 92 #ifdef G4ANALYSIS_USE >> 93 if(!af) return; >> 94 >> 95 // Creating a tree mapped to an histogram file. >> 96 fileName[1] = fileName[0] + "." + fileType; >> 97 G4bool readOnly = false; >> 98 G4bool createNew = true; >> 99 AIDA::ITreeFactory* tf = af->createTreeFactory(); >> 100 tree = tf->create(fileName[1], fileType, readOnly, createNew, fileOption); >> 101 delete tf; >> 102 if(!tree) { >> 103 G4cout << " HistoManager::book :" >> 104 << " problem creating the AIDA tree with " >> 105 << " storeName = " << fileName[1] >> 106 << " storeType = " << fileType >> 107 << " readOnly = " << readOnly >> 108 << " createNew = " << createNew >> 109 << " options = " << fileOption >> 110 << G4endl; >> 111 return; >> 112 } >> 113 >> 114 // Creating a histogram factory, whose histograms will be handled by the tree >> 115 AIDA::IHistogramFactory* hf = af->createHistogramFactory(*tree); >> 116 >> 117 // create selected histograms >> 118 for (G4int k=1; k<MaxHisto; k++) { >> 119 if (exist[k]) { >> 120 histo[k] = hf->createHistogram1D( Label[k], Title[k], >> 121 Nbins[k], Vmin[k], Vmax[k]); >> 122 factoryOn = true; >> 123 } >> 124 } >> 125 delete hf; >> 126 >> 127 // Creating a ntuple factory, handled by the tree >> 128 AIDA::ITupleFactory* ntf = af->createTupleFactory(*tree); >> 129 >> 130 // create selected ntuples >> 131 for (G4int k=1; k<MaxNtupl; k++) { >> 132 if (existNt[k]) { >> 133 ntupl[k] = ntf->create( LabelNt[k], TitleNt[k], ColumnNt[k]); >> 134 factoryOn = true; >> 135 } >> 136 } >> 137 >> 138 delete ntf; >> 139 >> 140 if (factoryOn) >> 141 G4cout << "\n----> Histogram Tree is opened in " << fileName[1] << G4endl; >> 142 #endif >> 143 } >> 144 >> 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 146 >> 147 void HistoManager::save() >> 148 { >> 149 #ifdef G4ANALYSIS_USE >> 150 if (factoryOn) { >> 151 saveAscii(); // Write ascii file, if any >> 152 tree->commit(); // Writing the histograms to the file >> 153 tree->close(); // and closing the tree (and the file) >> 154 G4cout << "\n----> Histogram Tree is saved in " << fileName[1] << G4endl; >> 155 >> 156 delete tree; >> 157 tree = 0; >> 158 factoryOn = false; >> 159 } >> 160 #endif >> 161 } >> 162 >> 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 164 >> 165 void HistoManager::FillHisto(G4int ih, G4double xbin, G4double weight) >> 166 { >> 167 if (ih >= MaxHisto) { >> 168 G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih >> 169 << " xbin= " << xbin << " weight= " << weight << G4endl; >> 170 return; >> 171 } >> 172 #ifdef G4ANALYSIS_USE >> 173 if(exist[ih]) histo[ih]->fill(xbin/Unit[ih], weight); >> 174 #endif >> 175 } >> 176 >> 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 178 >> 179 void HistoManager::SetHisto(G4int ih, >> 180 G4int nbins, G4double valmin, G4double valmax, const G4String& unit) >> 181 { >> 182 if (ih < 1 || ih >= MaxHisto) { >> 183 G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih >> 184 << "does not exist" << G4endl; >> 185 return; >> 186 } >> 187 >> 188 const G4String id[] = { "0", "1", "2", "3", "4", "5" }; >> 189 const G4String title[] = >> 190 { "dummy", //0 >> 191 "total Evis in Ecal", //1 >> 192 "total Edep in Ecal", //2 >> 193 "Evis profile", //3 >> 194 "Edep profile", //4 >> 195 "Nb of Radiation Length" //5 >> 196 }; >> 197 >> 198 G4String titl = title[ih]; >> 199 G4double vmin = valmin, vmax = valmax; >> 200 Unit[ih] = 1.; >> 201 >> 202 if (unit != "none") { >> 203 titl = title[ih] + " (" + unit + ")"; >> 204 Unit[ih] = G4UnitDefinition::GetValueOf(unit); >> 205 vmin = valmin/Unit[ih]; vmax = valmax/Unit[ih]; >> 206 } >> 207 >> 208 exist[ih] = true; >> 209 Label[ih] = id[ih]; >> 210 Title[ih] = titl; >> 211 Nbins[ih] = nbins; >> 212 Vmin[ih] = vmin; >> 213 Vmax[ih] = vmax; >> 214 Width[ih] = (valmax-valmin)/nbins; >> 215 >> 216 G4cout << "----> SetHisto " << ih << ": " << titl << "; " >> 217 << nbins << " bins from " >> 218 << vmin << " " << unit << " to " << vmax << " " << unit << G4endl; >> 219 >> 220 } >> 221 >> 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 223 >> 224 void HistoManager::Normalize(G4int ih, G4double fac) 47 { 225 { >> 226 if (ih >= MaxHisto) { >> 227 G4cout << "---> warning from HistoManager::Normalize() : histo " << ih >> 228 << " fac= " << fac << G4endl; >> 229 return; >> 230 } >> 231 #ifdef G4ANALYSIS_USE >> 232 if(exist[ih]) histo[ih]->scale(fac); >> 233 #endif >> 234 } >> 235 >> 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 237 >> 238 void HistoManager::RemoveHisto(G4int ih) >> 239 { >> 240 if (ih >= MaxHisto) { >> 241 G4cout << "---> warning from HistoManager::RemoveHisto() : histo " << ih >> 242 << "does not exist" << G4endl; >> 243 return; >> 244 } >> 245 >> 246 histo[ih] = 0; exist[ih] = false; >> 247 } >> 248 >> 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 250 >> 251 void HistoManager::PrintHisto(G4int ih) >> 252 { >> 253 if (ih < MaxHisto) { ascii[ih] = true; ascii[0] = true; } >> 254 else >> 255 G4cout << "---> warning from HistoManager::PrintHisto() : histo " << ih >> 256 << "does not exist" << G4endl; 48 } 257 } 49 258 50 //....oooOO0OOooo........oooOO0OOooo........oo 259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 51 260 52 void HistoManager::Book() << 261 #include <fstream> >> 262 >> 263 void HistoManager::saveAscii() 53 { 264 { 54 // Create or get analysis manager << 265 #ifdef G4ANALYSIS_USE 55 // The choice of analysis technology is done << 266 56 // in HistoManager.hh << 267 if (!ascii[0]) return; 57 G4AnalysisManager* analysisManager = G4Analy << 268 58 analysisManager->SetDefaultFileType("root"); << 269 G4String name = fileName[0] + ".ascii"; 59 analysisManager->SetFileName(fFileName); << 270 std::ofstream File(name, std::ios::out); 60 analysisManager->SetVerboseLevel(1); << 271 File.setf( std::ios::scientific, std::ios::floatfield ); 61 analysisManager->SetActivation(true); // << 272 62 << 273 //write selected histograms 63 // Define histograms start values << 274 for (G4int ih=0; ih<MaxHisto; ih++) { 64 const G4int kMaxHisto = 6; << 275 if (exist[ih] && ascii[ih]) { 65 const G4String id[] = {"0", "1", "2", "3" , << 276 File << "\n 1D histogram " << ih << ": " << Title[ih] >> 277 << "\n \n \t X \t\t Y" << G4endl; >> 278 >> 279 for (G4int iBin=0; iBin<Nbins[ih]; iBin++) { >> 280 File << " " << iBin << "\t" >> 281 << 0.5*(histo[ih]->axis().binLowerEdge(iBin) + >> 282 histo[ih]->axis().binUpperEdge(iBin)) << "\t" >> 283 << histo[ih]->binHeight(iBin) >> 284 << G4endl; >> 285 } >> 286 } >> 287 } >> 288 #endif >> 289 } >> 290 >> 291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 292 >> 293 void HistoManager::SetNtuple(G4int nt) >> 294 { >> 295 if (nt < 1 || nt >= MaxNtupl) { >> 296 G4cout << "---> warning from HistoManager::SetNtuple() : Ntuple " << nt >> 297 << "does not exist" << G4endl; >> 298 return; >> 299 } >> 300 >> 301 const G4String id[] = { "100", "101" }; 66 const G4String title[] = 302 const G4String title[] = 67 { "dummy", << 303 { "dummy", //0 68 "total Etot in Ecal", //1 << 304 "Energy deposit in subModule" //1 69 "total Evis in Ecal", << 70 "Etot profile", << 71 "Evis profile", << 72 "Evis per fiber" << 73 }; 305 }; 74 << 306 75 // Default values (to be reset via /analysis << 307 G4String column[MaxNtupl]; 76 G4int nbins = 100; << 308 77 G4double vmin = 0.; << 309 column[0] = " int dum=0 "; //0 78 G4double vmax = 100.; << 310 79 << 311 column[1] = 80 // Create all histograms as inactivated << 312 "double Evis0, Evis1, Evis2, Evis3, Evis4, Evis5, Evis6, Evis7, Evis8, Evis9"; 81 // as we have not yet set nbins, vmin, vmax << 313 column[1] += 82 for (G4int k=0; k<kMaxHisto; k++) { << 314 ", Evis10, Evis11, Evis12, Evis13, Evis14, Evis15, Evis16, Evis17"; 83 G4int ih = analysisManager->CreateH1(id[k] << 315 column[1] += 84 analysisManager->SetH1Activation(ih, false << 316 ", Edep0, Edep1, Edep2, Edep3, Edep4, Edep5, Edep6, Edep7, Edep8, Edep9"; >> 317 column[1] += >> 318 ", Edep10, Edep11, Edep12, Edep13, Edep14, Edep15, Edep16, Edep17"; >> 319 >> 320 G4String titl = title[nt]; >> 321 >> 322 existNt[nt] = true; >> 323 LabelNt[nt] = id[nt]; >> 324 TitleNt[nt] = titl; >> 325 ColumnNt[nt] = column[nt]; >> 326 >> 327 G4cout << "----> SetNtuple " << nt << ": " << titl << "; " << G4endl; >> 328 >> 329 } >> 330 >> 331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 332 >> 333 void HistoManager::FillNtuple(G4int nt, G4int column, G4double value) >> 334 { >> 335 if (nt >= MaxNtupl) { >> 336 G4cout << "---> warning from HistoManager::FillNtuple() : Ntuple " << nt >> 337 << " does not exist " << column << value << G4endl; >> 338 return; 85 } 339 } >> 340 #ifdef G4ANALYSIS_USE >> 341 if(existNt[nt]) ntupl[nt]->fill(column, value); >> 342 #endif 86 } 343 } 87 344 88 //....oooOO0OOooo........oooOO0OOooo........oo 345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 346 >> 347 void HistoManager::AddRowNtuple(G4int nt) >> 348 { >> 349 if (nt >= MaxNtupl) { >> 350 G4cout << "---> warning from HistoManager::AddRowNtuple() : Ntuple " << nt >> 351 << " do not exist" << G4endl; >> 352 return; >> 353 } >> 354 #ifdef G4ANALYSIS_USE >> 355 if(existNt[nt]) ntupl[nt]->addRow(); >> 356 #endif >> 357 } >> 358 >> 359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 360 >> 361 89 362