gammaray_telescope/src/GammaRayTelEventAction.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 //      ------------ GammaRayTelEventAction  ------
0033 //           by  R.Giannitrapani, F.Longo & G.Santin (13 nov 2000)
0034 //
0035 // - inclusion of Digits by F.Longo & R.Giannitrapani (24 oct 2001)
0036 //
0037 // - Modification of analysis management by G.Santin (18 Nov 2001)
0038 //
0039 // ************************************************************
0040
0041 #include "GammaRayTelEventAction.hh"
0042 #include "GammaRayTelTrackerHit.hh"
0043 #include "GammaRayTelAnticoincidenceHit.hh"
0044 #include "GammaRayTelCalorimeterHit.hh"
0045 #include "GammaRayTelAnalysis.hh"
0046 #include "GammaRayTelDigi.hh"
0047 #include "GammaRayTelDigitizer.hh"
0048
0049 #include "G4DigiManager.hh"
0050 #include "G4Event.hh"
0051 #include "G4EventManager.hh"
0052 #include "G4HCofThisEvent.hh"
0053 #include "G4ios.hh"
0054 #include "G4SDManager.hh"
0055 #include "G4SystemOfUnits.hh"
0056 #include "G4UImanager.hh"
0057 #include "G4UnitsTable.hh"
0058 #include "G4VHitsCollection.hh"
0059 #include "Randomize.hh"
0060
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0062
0063 // This file is a global variable in which we store energy deposition per hit
0064 // and other relevant information
0065
0066 GammaRayTelEventAction::GammaRayTelEventAction(GammaRayTelRunAction *runAction) : theRunAction(runAction) {
0067     auto *digitizer = new GammaRayTelDigitizer("GammaRayTelDigitizer");
0068     G4DigiManager::GetDMpointer()->AddNewModule(digitizer);
0069 }
0070
0071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0072
0073 GammaRayTelEventAction::~GammaRayTelEventAction() {
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0077
0078 void GammaRayTelEventAction::BeginOfEventAction(const G4Event *event) {
0079     G4int eventIdentifer = event->GetEventID();
0080     G4cout << "Event: " << eventIdentifer << G4endl;
0081     auto *sensitiveDetectorManager = G4SDManager::GetSDMpointer();
0082
0083     if (trackerCollectionID == -1) {
0084         trackerCollectionID = sensitiveDetectorManager->GetCollectionID("TrackerCollection");
0085     }
0086     if (anticoincidenceCollectionID == -1) {
0087         anticoincidenceCollectionID = sensitiveDetectorManager->GetCollectionID("AnticoincidenceCollection");
0088     }
0089     if (calorimeterCollectionID == -1) {
0090         calorimeterCollectionID = sensitiveDetectorManager->GetCollectionID("CalorimeterCollection");
0091     }
0092 }
0093
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0095
0096 void GammaRayTelEventAction::EndOfEventAction(const G4Event *event) {
0097     G4int eventIdentifier = event->GetEventID();
0098
0099     if (theRunAction == nullptr) {
0100         G4Exception("GammaRayTelEventAction::BeginOfEventAction()", "GTR0001", FatalException, "Null pointer to Run Action: this should not be");
0101     }
0102
0103 #ifdef G4STORE_DATA
0104     auto *outputFile = theRunAction->GetOutputFile();
0105 #endif
0106
0107     auto *collection = event->GetHCofThisEvent();
0108     GammaRayTelTrackerHitsCollection *trackerCollection{nullptr}; // TKR
0109     // GammaRayTelCalorimeterHitsCollection* calorimeterCollection{nullptr}; // CAL
0110     // GammaRayTelAnticoincidenceHitsCollection* anticoincidenceCollection{nullptr}; // ACD
0111
0112     auto *fDM = G4DigiManager::GetDMpointer();
0113     GammaRayTelAnalysis *analysis = GammaRayTelAnalysis::getInstance();
0114
0115     if (collection != nullptr) {
0116         trackerCollection = (GammaRayTelTrackerHitsCollection*) (collection->GetHC(trackerCollectionID));
0117         // calorimeterCollection = (GammaRayTelCalorimeterHitsCollection*) (collection->GetHC(calorimeterCollectionID));
0118         // anticoincidenceCollection = (GammaRayTelAnticoincidenceHitsCollection*) (collection->GetHC(anticoincidenceCollectionID));
0119
0120         if (trackerCollection != nullptr) {
0121             G4int numberOfHits = trackerCollection->entries();
0122             G4cout << "Number of tracker hits in this event: " << numberOfHits << G4endl;
0123
0124             G4double depositedEnergy{0};
0125             G4int stripNumber;
0126             G4int planeNumber;
0127             G4int isXPlane;
0128
0129             // This is a cycle on all the tracker hits of this event
0130
0131             for (auto i = 0; i < numberOfHits; i++) {
0132                 // Here we put the hit data in an ASCII file for later analysis
0133                 depositedEnergy = (*trackerCollection)[i]->GetDepositedEnergy();
0134                 stripNumber = (*trackerCollection)[i]->GetStripNumber();
0135                 planeNumber = (*trackerCollection)[i]->GetSiliconPlaneNumber();
0136                 isXPlane = (*trackerCollection)[i]->GetPlaneType();
0137
0138 #ifdef G4STORE_DATA
0139                 (*outputFile) << std::setw(7) << eventIdentifier
0140                     << " " << depositedEnergy/keV
0141                     << " " << stripNumber
0142                     << " " << planeNumber
0143                     << " " << isXPlane
0144                     << " " << (*trackerCollection)[i]->GetPosition().x() / mm
0145                     << " " << (*trackerCollection)[i]->GetPosition().y() / mm
0146                     << " " << (*trackerCollection)[i]->GetPosition().z() / mm
0147                     << " " << G4endl;
0148 #else
0149                 G4cout << std::setw(7) << eventIdentifier
0150                     << " " << depositedEnergy / keV
0151                     << " " << stripNumber
0152                     << " " << planeNumber
0153                     << " " << isXPlane
0154                     << " " << (*trackerCollection)[i]->GetPosition().x() / mm
0155                     << " " << (*trackerCollection)[i]->GetPosition().y() / mm
0156                     << " " << (*trackerCollection)[i]->GetPosition().z() / mm
0157                     << " " << G4endl;
0158 #endif
0159
0160                 // Here we fill the histograms of the Analysis manager
0161                 if (isXPlane != 0) {
0162                     if (analysis->GetHisto2DMode() == "position") {
0163                         analysis->InsertPositionXZ((*trackerCollection)[i]->GetPosition().x() / mm, (*trackerCollection)[i]->GetPosition().z() / mm);
0164                     } else {
0165                         analysis->InsertPositionXZ(stripNumber, planeNumber);
0166                     }
0167
0168                     if (planeNumber == 0) {
0169                         analysis->InsertEnergy(depositedEnergy / keV);
0170                     }
0171                     analysis->InsertHits(planeNumber);
0172                 } else {
0173                     if (analysis->GetHisto2DMode() == "position") {
0174                         analysis->InsertPositionYZ((*trackerCollection)[i]->GetPosition().y() / mm, (*trackerCollection)[i]->GetPosition().z() / mm);
0175                     } else {
0176                         analysis->InsertPositionYZ(stripNumber, planeNumber);
0177                     }
0178                     if (planeNumber == 0) {
0179                         analysis->InsertEnergy(depositedEnergy / keV);
0180                     }
0181                     analysis->InsertHits(planeNumber);
0182                 }
0183                 analysis->setNtuple(depositedEnergy / keV, planeNumber,
0184                     (*trackerCollection)[i]->GetPosition().x() / mm,
0185                     (*trackerCollection)[i]->GetPosition().y() / mm,
0186                     (*trackerCollection)[i]->GetPosition().z() / mm
0187                 );
0188             }
0189             analysis->EndOfEvent(numberOfHits);
0190         }
0191
0192         auto *myDM = (GammaRayTelDigitizer*) fDM->FindDigitizerModule("GammaRayTelDigitizer");
0193         myDM->Digitize();
0194
0195 #ifdef G4STORE_DATA
0196         // The whole block is needed only when outputFile is active;
0197         // protect block to avoid compilations warnings from gcc4.6, Gunter Folger
0198
0199         auto digitsCollectionIdentifier = fDM->GetDigiCollectionID("DigitsCollection");
0200         // G4cout << "Digits collection: " << digitsCollectionIdentifier << G4endl;
0201
0202         auto *digitsCollection = (GammaRayTelDigitsCollection*) fDM->GetDigiCollection(digitsCollectionIdentifier);
0203
0204         if (digitsCollection != nullptr) {
0205             auto numberOfDigits = digitsCollection->entries();
0206             // G4cout << "Total number of digits: " << numberOfDigits << G4endl;
0207
0208             G4int stripNumber;
0209             G4int planeNumber;
0210             G4int isXPlane;
0211
0212             for (auto i = 0; i < numberOfDigits; i++) {
0213                 // Here we put the digi data in an ASCII file for later analysis
0214                 stripNumber = (*digitsCollection)[i]->GetStripNumber();
0215                 planeNumber = (*digitsCollection)[i]->GetPlaneNumber();
0216                 isXPlane = (*digitsCollection)[i]->GetPlaneType();
0217
0218                 (*outputFile) << std::setw(7)
0219                     << eventIdentifier
0220                     << " " << stripNumber
0221                     << " " << planeNumber
0222                     << " " << isXPlane
0223                     << " " << G4endl;
0224             }
0225         }
0226 #endif
0227     }
0228 }