gammaray_telescope/src/GammaRayTelAnalysis.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 //
0027 // ------------------------------------------------------------
0028 //      GEANT 4 class implementation file
0029 //      CERN Geneva Switzerland
0030 //
0031 //
0032 //      ------------ GammaRayAnalysisManager  ------
0033 //           by R.Giannitrapani, F.Longo & G.Santin (03 dic 2000)
0034 //
0035 // 03.04.2013 F.Longo/L.Pandola
0036 // - migrated to G4tools
0037 //
0038 // 29.05.2003 F.Longo
0039 // - Anaphe 5.0.5 compliant
0040 //
0041 // 18.06.2002 R.Giannitrapani, F.Longo & G.Santin
0042 // - new release for Anaphe 4.0.3
0043 //
0044 // 07.12.2001 A.Pfeiffer
0045 // - integrated Guy's addition of the ntuple
0046 //
0047 // 06.12.2001 A.Pfeiffer
0048 // - updating to new design (singleton)
0049 //
0050 // 22.11.2001 G.Barrand
0051 // - Adaptation to AIDA
0052 //
0053 // ************************************************************
0054
0055 #include <fstream>
0056 #include <iomanip>
0057
0058 #include "GammaRayTelAnalysis.hh"
0059 #include "GammaRayTelAnalysisMessenger.hh"
0060 #include "GammaRayTelDetectorConstruction.hh"
0061
0062 #include "G4RunManager.hh"
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0065
0066 GammaRayTelAnalysis *GammaRayTelAnalysis::instance{nullptr};
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0069
0070 GammaRayTelAnalysis::GammaRayTelAnalysis() : detector(nullptr), histo2DMode("strip") {
0071     detector = static_cast<const GammaRayTelDetectorConstruction*>(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0072
0073     // Define the messenger and the analysis system
0074     analysisMessenger = new GammaRayTelAnalysisMessenger(this);
0075     histogramFileName = "gammaraytel";
0076 }
0077
0078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0079
0080 GammaRayTelAnalysis::~GammaRayTelAnalysis() {
0081     Finish();
0082 }
0083
0084 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0085
0086 void GammaRayTelAnalysis::Init() {
0087 }
0088
0089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0090
0091 void GammaRayTelAnalysis::Finish() {
0092     delete analysisMessenger;
0093     analysisMessenger = nullptr;
0094 }
0095
0096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0097
0098 auto GammaRayTelAnalysis::getInstance() -> GammaRayTelAnalysis* {
0099     if (instance == nullptr) {
0100         instance = new GammaRayTelAnalysis();
0101     }
0102     return instance;
0103 }
0104
0105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0106
0107 // This function fill the 2d histogram of the XZ positions
0108 void GammaRayTelAnalysis::InsertPositionXZ(G4double x, G4double z) {
0109     auto *manager = G4AnalysisManager::Instance();
0110     manager->FillH2(1, x, z);
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0114
0115 // This function fill the 2d histogram of the YZ positions
0116 void GammaRayTelAnalysis::InsertPositionYZ(G4double y, G4double z) {
0117     auto *manager = G4AnalysisManager::Instance();
0118     manager->FillH2(2, y, z);
0119 }
0120
0121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0122
0123 // This function fill the 1d histogram of the energy released in the last Si plane
0124 void GammaRayTelAnalysis::InsertEnergy(G4double energy) {
0125     auto *manager = G4AnalysisManager::Instance();
0126     manager->FillH1(1, energy);
0127 }
0128
0129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0130
0131 // This function fill the 1d histogram of the hits distribution along the TKR planes
0132 void GammaRayTelAnalysis::InsertHits(G4int planeNumber) {
0133     auto *manager = G4AnalysisManager::Instance();
0134     manager->FillH1(2, planeNumber);
0135 }
0136
0137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0138
0139 void GammaRayTelAnalysis::setNtuple(G4double energy, G4int planeNumber, G4double x, G4double y, G4double z) {
0140     auto *manager = G4AnalysisManager::Instance();
0141     manager->FillNtupleDColumn(0, energy);
0142     manager->FillNtupleDColumn(1, planeNumber);
0143     manager->FillNtupleDColumn(2, x);
0144     manager->FillNtupleDColumn(3, y);
0145     manager->FillNtupleDColumn(4, z);
0146     manager->AddNtupleRow();
0147 }
0148
0149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0150
0151 /*
0152  This member reset the histograms and it is called at the begin of each run;
0153  here we put the inizialization so that the histograms have
0154  always the right dimensions depending from the detector geometry
0155  */
0156 void GammaRayTelAnalysis::BeginOfRun() {
0157     auto *manager = G4AnalysisManager::Instance();
0158     manager->SetDefaultFileType("root");
0159
0160     // Open an output file
0161
0162     G4cout << "Opening output file " << histogramFileName << " ... ";
0163     manager->OpenFile(histogramFileName);
0164     manager->SetFirstHistoId(1);
0165     G4cout << " done" << G4endl;
0166
0167     auto Nplane = detector->GetNbOfTKRLayers();
0168     auto numberOfStrips = detector->GetNbOfTKRStrips();
0169     auto numberOfTiles = detector->GetNbOfTKRTiles();
0170     auto sizeXY = detector->GetTKRSizeXY();
0171     auto sizeZ = detector->GetTKRSizeZ();
0172     auto N = numberOfStrips * numberOfTiles;
0173
0174     // Book 1D histograms
0175     //-------------------
0176
0177     constexpr auto NUMBER_OF_BINS{100};
0178     constexpr auto LOWER_BOUND{50};
0179     constexpr auto UPPER_BOUND{200};
0180
0181     // 1D histogram that store the energy deposition of the particle in the last (number 0) TKR X-plane
0182     manager->CreateH1("1", "Deposited energy in the last X plane (keV)", NUMBER_OF_BINS, LOWER_BOUND, UPPER_BOUND);
0183
0184     // 1D histogram that store the hits distribution along the TKR X-planes
0185     manager->CreateH1("2", "Hits distribution in TKR X planes", Nplane, 0, Nplane - 1);
0186
0187     // Book 2D histograms
0188     //-------------------
0189
0190     // 2D histogram that store the position (mm) of the hits (XZ projection)
0191
0192     if (histo2DMode == "strip") {
0193         manager->CreateH2("1", "Tracker hits XZ (strip, plane)", N, 0, N - 1, 2 * Nplane, 0, Nplane - 1);
0194     } else {
0195         manager->CreateH2("1", "Tracker hits XZ (x, z) in mm", G4int(sizeXY / 5), -sizeXY / 2, sizeXY / 2, G4int(sizeZ / 5), -sizeZ / 2, sizeZ / 2);
0196     }
0197
0198     // 2D histogram that store the position (mm) of the hits (YZ projection)
0199     if (histo2DMode == "strip") {
0200         manager->CreateH2("2", "Tracker hits YZ (strip, plane)", N, 0, N - 1, 2 * Nplane, 0, Nplane - 1);
0201     } else {
0202         manager->CreateH2("2", "Tracker hits YZ (x, z) in mm", G4int(sizeXY / 5), -sizeXY / 2, sizeXY / 2, G4int(sizeZ / 5), -sizeZ / 2, sizeZ / 2);
0203     }
0204
0205     // Book n-tuple (energy, plane, x, y, z)
0206     //------------------------------------------
0207     manager->CreateNtuple("1", "Track n-tuple");
0208     manager->CreateNtupleDColumn("energy");
0209     manager->CreateNtupleDColumn("plane");
0210     manager->CreateNtupleDColumn("x");
0211     manager->CreateNtupleDColumn("y");
0212     manager->CreateNtupleDColumn("z");
0213     manager->FinishNtuple();
0214 }
0215
0216 /*
0217  This member is called at the end of each run
0218  */
0219 void GammaRayTelAnalysis::EndOfRun() {
0220     // Save histograms
0221     auto *manager = G4AnalysisManager::Instance();
0222     manager->Write();
0223     manager->CloseFile();
0224 }
0225
0226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0227
0228 // This member is called at the end of every event
0229 void GammaRayTelAnalysis::EndOfEvent(G4int /* flag */) {
0230 }