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