TestEm17/src/HistoManager.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file electromagnetic/TestEm17/src/HistoManager.cc
0027 /// \brief Implementation of the HistoManager class
0028 //
0029 //
0030 //
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0033
0034 #include "HistoManager.hh"
0035
0036 #include "HistoMessenger.hh"
0037
0038 #include "G4PhysicalConstants.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4UnitsTable.hh"
0041
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0043
0044 HistoManager::HistoManager() : fHistoMessenger(0)
0045 {
0046   fileName[0] = "testem17";
0047   factoryOn = false;
0048   fNbHist = 0;
0049
0050   // histograms
0051   for (G4int k = 0; k < kMaxHisto; k++) {
0052     fHistId[k] = 0;
0053     fHistPt[k] = 0;
0054     fExist[k] = false;
0055     fUnit[k] = 1.0;
0056     fWidth[k] = 1.0;
0057     fAscii[k] = false;
0058   }
0059
0060   fHistoMessenger = new HistoMessenger(this);
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 HistoManager::~HistoManager()
0066 {
0067   delete fHistoMessenger;
0068 }
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071
0072 void HistoManager::Book()
0073 {
0074   // if no histos, do nothing
0075   if (fNbHist == 0) return;
0076
0077   // Create or get analysis manager
0078   // The choice of analysis technology is done via selection of a namespace
0079   // in HistoManager.hh
0080   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0081   analysisManager->SetDefaultFileType("root");
0082   analysisManager->SetVerboseLevel(0);
0083   G4String extension = analysisManager->GetFileType();
0084   fileName[1] = fileName[0] + "." + extension;
0085
0086   // Open an output file
0087   //
0088   G4bool fileOpen = analysisManager->OpenFile(fileName[0]);
0089   if (!fileOpen) {
0090     G4cout << "\n---> HistoManager::book(): cannot open " << fileName[1] << G4endl;
0091     return;
0092   }
0093
0094   // create selected histograms
0095   //
0096   analysisManager->SetFirstHistoId(1);
0097
0098   for (G4int k = 0; k < kMaxHisto; k++) {
0099     if (fExist[k]) {
0100       fHistId[k] = analysisManager->CreateH1(fLabel[k], fTitle[k], fNbins[k], fVmin[k], fVmax[k]);
0101       fHistPt[k] = analysisManager->GetH1(fHistId[k]);
0102       factoryOn = true;
0103     }
0104   }
0105
0106   if (factoryOn) G4cout << "\n----> Histogram file is opened in " << fileName[1] << G4endl;
0107 }
0108
0109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0110
0111 void HistoManager::Save()
0112 {
0113   if (factoryOn) {
0114     G4cout << "\n----> HistoManager::save " << G4endl;
0115     G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0116     analysisManager->Write();
0117     analysisManager->CloseFile();
0118     SaveAscii();  // Write fAscii file, if any
0119     G4cout << "\n----> Histograms are saved in " << fileName[1] << G4endl;
0120
0121     analysisManager->Clear();
0122     factoryOn = false;
0123   }
0124 }
0125
0126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0127
0128 void HistoManager::FillHisto(G4int ih, G4double e, G4double weight)
0129 {
0130   if (ih > kMaxHisto) {
0131     G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih
0132            << "does not fExist; e= " << e << " w= " << weight << G4endl;
0133     return;
0134   }
0135
0136   if (fHistPt[ih]) fHistPt[ih]->fill(e / fUnit[ih], weight);
0137 }
0138
0139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0140
0141 void HistoManager::SetHisto(G4int ih, G4int nbins, G4double vmin, G4double vmax,
0142                             const G4String& unit)
0143 {
0144   if (ih > kMaxHisto) {
0145     G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih << "does not fExist"
0146            << G4endl;
0147     return;
0148   }
0149
0150   const G4String id[] = {"0",  "1",  "2",  "3",  "4",  "5",  "6",  "7",  "8",  "9",
0151                          "10", "11", "12", "13", "14", "15", "16", "17", "18", "19",
0152                          "20", "21", "22", "23", "24", "25", "26", "27", "28", "29"};
0153   const G4String title[] = {
0154     "dummy",  // 0
0155     "log10(Eloss/Emu) muIonization",  // 1
0156     "log10(Eloss/Emu) muPair",  // 2
0157     "log10(Eloss/Emu) muBrems",  // 3
0158     "log10(Eloss/Emu) muNuclear",  // 4
0159     "log10(Eloss/Emu) hIonization",  // 5
0160     "log10(Eloss/Emu) hPair",  // 6
0161     "log10(Eloss/Emu) hBrems",  // 7
0162     "log10(Eloss/Emu) muToMuonPair",  // 8
0163     "dummy",  // 9
0164     "dummy",  // 10
0165     "log10(Eloss/Emu) muIonization",  // 11
0166     "log10(Eloss/Emu) muPair",  // 12
0167     "log10(Eloss/Emu) muBrems",  // 13
0168     "log10(Eloss/Emu) muNuclear",  // 14
0169     "dummy",  // 15
0170     "dummy",  // 16
0171     "dummy",  // 17
0172     "log10(Eloss/Emu) muToMuonPair",  // 18
0173     "dummy",  // 19
0174     "dummy",  // 20
0175     "CS(1/mm) vs ekin muIonisation",  // 21
0176     "CS(1/mm) vs ekin muPair",  // 22
0177     "CS(1/mm) vs ekin muBrems",  // 23
0178     "CS(1/mm) vs ekin muNuclear",  // 24
0179     "dummy",  // 25
0180     "dummy",  // 26
0181     "dummy",  // 27
0182     "CS(1/mm) vs ekin muToMuonPair",  // 28
0183     "dummy"  // 29
0184   };
0185
0186   G4String titl = title[ih];
0187   fUnit[ih] = 1.;
0188
0189   if (unit != "none") {
0190     titl = title[ih] + " (" + unit + ")";
0191     fUnit[ih] = G4UnitDefinition::GetValueOf(unit);
0192   }
0193
0194   fExist[ih] = true;
0195   fLabel[ih] = "h" + id[ih];
0196   fTitle[ih] = titl;
0197   fNbins[ih] = nbins;
0198   fVmin[ih] = vmin;
0199   fVmax[ih] = vmax;
0200   fWidth[ih] = fUnit[ih] * (vmax - vmin) / nbins;
0201
0202   fNbHist++;
0203
0204   G4cout << "----> SetHisto " << ih << ": " << titl << ";  " << nbins << " bins from " << vmin
0205          << " " << unit << " to " << vmax << " " << unit << G4endl;
0206 }
0207
0208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0209
0210 void HistoManager::Normalize(G4int ih, G4double fac)
0211 {
0212   if (ih >= kMaxHisto) {
0213     G4cout << "---> warning from HistoManager::Normalize() : histo " << ih << "  fac= " << fac
0214            << G4endl;
0215     return;
0216   }
0217
0218   if (fHistPt[ih]) fHistPt[ih]->scale(fac);
0219 }
0220
0221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0222
0223 void HistoManager::PrintHisto(G4int ih)
0224 {
0225   if (ih < kMaxHisto) {
0226     fAscii[ih] = true;
0227     fAscii[0] = true;
0228   }
0229   else
0230     G4cout << "---> warning from HistoManager::PrintHisto() : histo " << ih << "does not exist"
0231            << G4endl;
0232 }
0233
0234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0235
0236 #include <fstream>
0237
0238 void HistoManager::SaveAscii()
0239 {
0240   if (!fAscii[0]) return;
0241
0242   G4String name = fileName[0] + ".ascii";
0243   std::ofstream File(name, std::ios::out);
0244   if (!File) {
0245     G4cout << "\n---> HistoManager::saveAscii(): cannot open " << name << G4endl;
0246     return;
0247   }
0248
0249   File.setf(std::ios::scientific, std::ios::floatfield);
0250
0251   // write selected histograms
0252   for (G4int ih = 0; ih < kMaxHisto; ih++) {
0253     if (fHistPt[ih] && fAscii[ih]) {
0254       File << "\n  1D histogram " << ih << ": " << fTitle[ih] << "\n \n \t     X \t\t     Y"
0255            << G4endl;
0256
0257       for (G4int iBin = 0; iBin < fNbins[ih]; iBin++) {
0258         File << "  " << iBin << "\t" << fHistPt[ih]->axis().bin_center(iBin) << "\t"
0259              << fHistPt[ih]->bin_height(iBin) << G4endl;
0260       }
0261     }
0262   }
0263 }
0264
0265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......