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 HistoManager.cc
0027 /// \brief Implementation of the HistoManager class
0028
0029 #include "HistoManager.hh"
0030
0031 #include "HistoMessenger.hh"
0032
0033 #include "G4PhysicalConstants.hh"
0034 #include "G4SystemOfUnits.hh"
0035 #include "G4UnitsTable.hh"
0036
0037 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0038
0039 HistoManager::HistoManager() : fHistoMessenger(0)
0040 {
0041   fileName[0] = "testem17";
0042   factoryOn = false;
0043   fNbHist = 0;
0044
0045   // histograms
0046   for (G4int k = 0; k < kMaxHisto; k++) {
0047     fHistId[k] = 0;
0048     fHistPt[k] = 0;
0049     fExist[k] = false;
0050     fUnit[k] = 1.0;
0051     fWidth[k] = 1.0;
0052     fAscii[k] = false;
0053   }
0054
0055   fHistoMessenger = new HistoMessenger(this);
0056 }
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059
0060 HistoManager::~HistoManager()
0061 {
0062   delete fHistoMessenger;
0063 }
0064
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066
0067 void HistoManager::Book()
0068 {
0069   // if no histos, do nothing
0070   if (fNbHist == 0) return;
0071
0072   // Create or get analysis manager
0073   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0074   analysisManager->SetDefaultFileType("root");
0075   analysisManager->SetVerboseLevel(0);
0076   G4String extension = analysisManager->GetFileType();
0077   fileName[1] = fileName[0] + "." + extension;
0078
0079   // Open an output file
0080   //
0081   G4bool fileOpen = analysisManager->OpenFile(fileName[0]);
0082   if (!fileOpen) {
0083     G4cout << "\n---> HistoManager::book(): cannot open " << fileName[1] << G4endl;
0084     return;
0085   }
0086
0087   // create selected histograms
0088   //
0089   analysisManager->SetFirstHistoId(1);
0090
0091   for (G4int k = 0; k < kMaxHisto; k++) {
0092     if (fExist[k]) {
0093       fHistId[k] = analysisManager->CreateH1(fLabel[k], fTitle[k], fNbins[k], fVmin[k], fVmax[k]);
0094       fHistPt[k] = analysisManager->GetH1(fHistId[k]);
0095       factoryOn = true;
0096     }
0097   }
0098
0099   if (factoryOn) G4cout << "\n----> Histogram file is opened in " << fileName[1] << G4endl;
0100 }
0101
0102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0103
0104 void HistoManager::Save()
0105 {
0106   if (factoryOn) {
0107     G4cout << "\n----> HistoManager::save " << G4endl;
0108     G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0109     analysisManager->Write();
0110     analysisManager->CloseFile();
0111     SaveAscii();  // Write fAscii file, if any
0112     G4cout << "\n----> Histograms are saved in " << fileName[1] << G4endl;
0113
0114     analysisManager->Clear();
0115     factoryOn = false;
0116   }
0117 }
0118
0119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0120
0121 void HistoManager::FillHisto(G4int ih, G4double e, G4double weight)
0122 {
0123   if (ih > kMaxHisto) {
0124     G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih
0125            << "does not fExist; e= " << e << " w= " << weight << G4endl;
0126     return;
0127   }
0128
0129   if (fHistPt[ih]) fHistPt[ih]->fill(e / fUnit[ih], weight);
0130 }
0131
0132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0133
0134 void HistoManager::SetHisto(G4int ih, G4int nbins, G4double vmin, G4double vmax,
0135                             const G4String& unit)
0136 {
0137   if (ih > kMaxHisto) {
0138     G4cout << "---> warning from HistoManager::SetHisto() : histo " << ih << "does not fExist"
0139            << G4endl;
0140     return;
0141   }
0142
0143   const G4String id[] = {"0",  "1",  "2",  "3",  "4",  "5",  "6",  "7",  "8",  "9",
0144                          "10", "11", "12", "13", "14", "15", "16", "17", "18", "19",
0145                          "20", "21", "22", "23", "24", "25", "26", "27", "28", "29"};
0146   const G4String title[] = {
0147     "dummy",  // 0
0148     "log10(Eloss/Emu) muIonization",  // 1
0149     "log10(Eloss/Emu) muPair",  // 2
0150     "log10(Eloss/Emu) muBrems",  // 3
0151     "log10(Eloss/Emu) muNuclear",  // 4
0152     "log10(Eloss/Emu) hIonization",  // 5
0153     "log10(Eloss/Emu) hPair",  // 6
0154     "log10(Eloss/Emu) hBrems",  // 7
0155     "log10(Eloss/Emu) muToMuonPair",  // 8
0156     "dummy",  // 9
0157     "dummy",  // 10
0158     "log10(Eloss/Emu) muIonization",  // 11
0159     "log10(Eloss/Emu) muPair",  // 12
0160     "log10(Eloss/Emu) muBrems",  // 13
0161     "log10(Eloss/Emu) muNuclear",  // 14
0162     "dummy",  // 15
0163     "dummy",  // 16
0164     "dummy",  // 17
0165     "log10(Eloss/Emu) muToMuonPair",  // 18
0166     "dummy",  // 19
0167     "dummy",  // 20
0168     "CS(1/mm) vs ekin muIonisation",  // 21
0169     "CS(1/mm) vs ekin muPair",  // 22
0170     "CS(1/mm) vs ekin muBrems",  // 23
0171     "CS(1/mm) vs ekin muNuclear",  // 24
0172     "dummy",  // 25
0173     "dummy",  // 26
0174     "dummy",  // 27
0175     "CS(1/mm) vs ekin muToMuonPair",  // 28
0176     "dummy"  // 29
0177   };
0178
0179   G4String titl = title[ih];
0180   fUnit[ih] = 1.;
0181
0182   if (unit != "none") {
0183     titl = title[ih] + " (" + unit + ")";
0184     fUnit[ih] = G4UnitDefinition::GetValueOf(unit);
0185   }
0186
0187   fExist[ih] = true;
0188   fLabel[ih] = "h" + id[ih];
0189   fTitle[ih] = titl;
0190   fNbins[ih] = nbins;
0191   fVmin[ih] = vmin;
0192   fVmax[ih] = vmax;
0193   fWidth[ih] = fUnit[ih] * (vmax - vmin) / nbins;
0194
0195   fNbHist++;
0196
0197   G4cout << "----> SetHisto " << ih << ": " << titl << ";  " << nbins << " bins from " << vmin
0198          << " " << unit << " to " << vmax << " " << unit << G4endl;
0199 }
0200
0201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0202
0203 void HistoManager::Normalize(G4int ih, G4double fac)
0204 {
0205   if (ih >= kMaxHisto) {
0206     G4cout << "---> warning from HistoManager::Normalize() : histo " << ih << "  fac= " << fac
0207            << G4endl;
0208     return;
0209   }
0210
0211   if (fHistPt[ih]) fHistPt[ih]->scale(fac);
0212 }
0213
0214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0215
0216 void HistoManager::PrintHisto(G4int ih)
0217 {
0218   if (ih < kMaxHisto) {
0219     fAscii[ih] = true;
0220     fAscii[0] = true;
0221   }
0222   else
0223     G4cout << "---> warning from HistoManager::PrintHisto() : histo " << ih << "does not exist"
0224            << G4endl;
0225 }
0226
0227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0228
0229 #include <fstream>
0230
0231 void HistoManager::SaveAscii()
0232 {
0233   if (!fAscii[0]) return;
0234
0235   G4String name = fileName[0] + ".ascii";
0236   std::ofstream File(name, std::ios::out);
0237   if (!File) {
0238     G4cout << "\n---> HistoManager::saveAscii(): cannot open " << name << G4endl;
0239     return;
0240   }
0241
0242   File.setf(std::ios::scientific, std::ios::floatfield);
0243
0244   // write selected histograms
0245   for (G4int ih = 0; ih < kMaxHisto; ih++) {
0246     if (fHistPt[ih] && fAscii[ih]) {
0247       File << "\n  1D histogram " << ih << ": " << fTitle[ih] << "\n \n \t     X \t\t     Y"
0248            << G4endl;
0249
0250       for (G4int iBin = 0; iBin < fNbins[ih]; iBin++) {
0251         File << "  " << iBin << "\t" << fHistPt[ih]->axis().bin_center(iBin) << "\t"
0252              << fHistPt[ih]->bin_height(iBin) << G4endl;
0253       }
0254     }
0255   }
0256 }
0257
0258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......