AnaEx02/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 "G4UnitsTable.hh"
0032
0033 #include <CLHEP/Units/SystemOfUnits.h>
0034 #include <TFile.h>
0035 #include <TH1D.h>
0036 #include <TTree.h>
0037
0038 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0039
0040 HistoManager::~HistoManager()
0041 {
0042   delete fRootFile;
0043 }
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 void HistoManager::Book()
0048 {
0049   // Creating a tree container to handle histograms and ntuples.
0050   // This tree is associated to an output file.
0051   //
0052   G4String fileName = "AnaEx02.root";
0053   fRootFile = new TFile(fileName, "RECREATE");
0054   if (fRootFile == nullptr) {
0055     G4cout << " HistoManager::Book :"
0056            << " problem creating the ROOT TFile " << G4endl;
0057     return;
0058   }
0059
0060   // id = 0
0061   fHisto[0] = new TH1D("EAbs", "Edep in absorber (MeV)", 100, 0., 800 * CLHEP::MeV);
0062   // id = 1
0063   fHisto[1] = new TH1D("EGap", "Edep in gap (MeV)", 100, 0., 100 * CLHEP::MeV);
0064   // id = 2
0065   fHisto[2] = new TH1D("LAbs", "trackL in absorber (mm)", 100, 0., 1 * CLHEP::m);
0066   // id = 3
0067   fHisto[3] = new TH1D("LGap", "trackL in gap (mm)", 100, 0., 50 * CLHEP::cm);
0068
0069   for (G4int i = 0; i < kMaxHisto; ++i) {
0070     if (fHisto[i] == nullptr) {
0071       G4cout << "\n can't create histo " << i << G4endl;
0072     }
0073   }
0074
0075   // create 1st ntuple
0076   fNtuple1 = new TTree("Ntuple1", "Edep");
0077   fNtuple1->Branch("Eabs", &fEabs, "Eabs/D");
0078   fNtuple1->Branch("Egap", &fEgap, "Egap/D");
0079
0080   // create 2nd ntuple
0081   fNtuple2 = new TTree("Ntuple2", "TrackL");
0082   fNtuple2->Branch("Labs", &fLabs, "Labs/D");
0083   fNtuple2->Branch("Lgap", &fLgap, "Lgap/D");
0084
0085   G4cout << "\n----> Output file is open in " << fileName << G4endl;
0086 }
0087
0088 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0089
0090 void HistoManager::Save()
0091 {
0092   if (fRootFile == nullptr) {
0093     return;
0094   }
0095
0096   fRootFile->Write();  // Writing the histograms to the file
0097   fRootFile->Close();  // and closing the tree (and the file)
0098
0099   G4cout << "\n----> Histograms and ntuples are saved\n" << G4endl;
0100 }
0101
0102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0103
0104 void HistoManager::FillHisto(G4int ih, G4double xbin, G4double weight)
0105 {
0106   if (ih >= kMaxHisto) {
0107     G4cerr << "---> warning from HistoManager::FillHisto() : histo " << ih
0108            << " does not exist. (xbin=" << xbin << " weight=" << weight << ")" << G4endl;
0109     return;
0110   }
0111   if (fHisto[ih] != nullptr) {
0112     fHisto[ih]->Fill(xbin, weight);
0113   }
0114 }
0115
0116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0117
0118 void HistoManager::Normalize(G4int ih, G4double fac)
0119 {
0120   if (ih >= kMaxHisto) {
0121     G4cout << "---> warning from HistoManager::Normalize() : histo " << ih
0122            << " does not exist. (fac=" << fac << ")" << G4endl;
0123     return;
0124   }
0125   if (fHisto[ih] != nullptr) {
0126     fHisto[ih]->Scale(fac);
0127   }
0128 }
0129
0130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0131
0132 void HistoManager::FillNtuple(G4double energyAbs, G4double energyGap, G4double trackLAbs,
0133                               G4double trackLGap)
0134 {
0135   fEabs = energyAbs;
0136   fEgap = energyGap;
0137   fLabs = trackLAbs;
0138   fLgap = trackLGap;
0139
0140   if (fNtuple1 != nullptr) {
0141     fNtuple1->Fill();
0142   }
0143   if (fNtuple2 != nullptr) {
0144     fNtuple2->Fill();
0145   }
0146 }
0147
0148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0149
0150 void HistoManager::PrintStatistic()
0151 {
0152   G4cout << "\n ----> print histograms statistic \n" << G4endl;
0153   for (auto h1 : fHisto) {
0154     const G4String name = h1->GetName();
0155
0156     G4String unitCategory;
0157     if (name[0] == 'E') {
0158       unitCategory = "Energy";
0159     }
0160     if (name[0] == 'L') {
0161       unitCategory = "Length";
0162     }
0163
0164     G4cout << name << ": mean = " << G4BestUnit(h1->GetMean(), unitCategory)
0165            << " rms = " << G4BestUnit(h1->GetRMS(), unitCategory) << G4endl;
0166   }
0167 }
0168
0169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......