Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/HGCal_testbeam/src/EventAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/HGCal_testbeam/src/EventAction.cc (Version 11.3.0) and /examples/advanced/HGCal_testbeam/src/EventAction.cc (Version 10.7.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 #include "EventAction.hh"                          26 #include "EventAction.hh"
 27                                                    27 
 28 #include "SiPMHit.hh"                              28 #include "SiPMHit.hh"
 29 #include "SiliconPixelHit.hh"                      29 #include "SiliconPixelHit.hh"
 30                                                    30 
 31 #include "G4Event.hh"                              31 #include "G4Event.hh"
 32 #include "G4SDManager.hh"                          32 #include "G4SDManager.hh"
 33 #include "G4GenericMessenger.hh"                   33 #include "G4GenericMessenger.hh"
 34 #include "G4AnalysisManager.hh"                <<  34 #include "g4root.hh"
 35                                                    35 
 36 //....oooOO0OOooo........oooOO0OOooo........oo     36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 37                                                    37 
 38 EventAction::EventAction() : G4UserEventAction     38 EventAction::EventAction() : G4UserEventAction() { DefineCommands(); }
 39                                                    39 
 40 //....oooOO0OOooo........oooOO0OOooo........oo     40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 41                                                    41 
 42 EventAction::~EventAction() {}                     42 EventAction::~EventAction() {}
 43                                                    43 
 44 //....oooOO0OOooo........oooOO0OOooo........oo     44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 45                                                    45 
 46 void EventAction::BeginOfEventAction(const G4E     46 void EventAction::BeginOfEventAction(const G4Event *) {
 47   fPrimariesPDG.clear();                           47   fPrimariesPDG.clear();
 48   fPrimariesEnergy.clear();                        48   fPrimariesEnergy.clear();
 49   fPrimariesX.clear();                             49   fPrimariesX.clear();
 50   fPrimariesY.clear();                             50   fPrimariesY.clear();
 51   fPrimariesZ.clear();                             51   fPrimariesZ.clear();
 52                                                    52 
 53   fSiHitsID.clear();                               53   fSiHitsID.clear();
 54   fSiHitsX.clear();                                54   fSiHitsX.clear();
 55   fSiHitsY.clear();                                55   fSiHitsY.clear();
 56   fSiHitsZ.clear();                                56   fSiHitsZ.clear();
 57   fSiHitsEdep.clear();                             57   fSiHitsEdep.clear();
 58   fSiHitsEdepNonIonising.clear();                  58   fSiHitsEdepNonIonising.clear();
 59   fSiHitsTOA.clear();                              59   fSiHitsTOA.clear();
 60   fSiHitsTOAlast.clear();                          60   fSiHitsTOAlast.clear();
 61   fSiHitsType.clear();                             61   fSiHitsType.clear();
 62                                                    62 
 63   fSiPMhitsID.clear();                             63   fSiPMhitsID.clear();
 64   fSiPMhitsX.clear();                              64   fSiPMhitsX.clear();
 65   fSiPMhitsY.clear();                              65   fSiPMhitsY.clear();
 66   fSiPMhitsZ.clear();                              66   fSiPMhitsZ.clear();
 67   fSiPMhitsEdep.clear();                           67   fSiPMhitsEdep.clear();
 68   fSiPMhitsEdepNonIonising.clear();                68   fSiPMhitsEdepNonIonising.clear();
 69   fSiPMhitsTOA.clear();                            69   fSiPMhitsTOA.clear();
 70   fSiPMhitsType.clear();                           70   fSiPMhitsType.clear();
 71 }                                                  71 }
 72                                                    72 
 73 //....oooOO0OOooo........oooOO0OOooo........oo     73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 74                                                    74 
 75 void EventAction::EndOfEventAction(const G4Eve     75 void EventAction::EndOfEventAction(const G4Event *event) {
 76   // sanity check                                  76   // sanity check
 77   if (event->GetNumberOfPrimaryVertex() == 0)      77   if (event->GetNumberOfPrimaryVertex() == 0)
 78     return;                                        78     return;
 79                                                    79 
 80   auto analysisManager = G4AnalysisManager::In     80   auto analysisManager = G4AnalysisManager::Instance();
 81   analysisManager->FillNtupleIColumn(0, event-     81   analysisManager->FillNtupleIColumn(0, event->GetEventID());
 82   // fill for all primary input particles          82   // fill for all primary input particles
 83   for (G4int iVertex = 0; iVertex < event->Get     83   for (G4int iVertex = 0; iVertex < event->GetNumberOfPrimaryVertex();
 84        iVertex++) {                                84        iVertex++) {
 85     auto vertex = event->GetPrimaryVertex(iVer     85     auto vertex = event->GetPrimaryVertex(iVertex);
 86     fPrimariesX.push_back(vertex->GetX0() / CL     86     fPrimariesX.push_back(vertex->GetX0() / CLHEP::cm);
 87     fPrimariesY.push_back(vertex->GetY0() / CL     87     fPrimariesY.push_back(vertex->GetY0() / CLHEP::cm);
 88     fPrimariesZ.push_back(vertex->GetZ0() / CL     88     fPrimariesZ.push_back(vertex->GetZ0() / CLHEP::cm);
 89     for (G4int iParticle = 0; iParticle < vert     89     for (G4int iParticle = 0; iParticle < vertex->GetNumberOfParticle();
 90          iParticle++) {                            90          iParticle++) {
 91       auto particle = vertex->GetPrimary(iPart     91       auto particle = vertex->GetPrimary(iParticle);
 92       fPrimariesPDG.push_back(particle->GetPDG     92       fPrimariesPDG.push_back(particle->GetPDGcode());
 93       fPrimariesEnergy.push_back(particle->Get     93       fPrimariesEnergy.push_back(particle->GetTotalEnergy() / CLHEP::GeV);
 94     }                                              94     }
 95   }                                                95   }
 96                                                    96 
 97   auto hce = event->GetHCofThisEvent();            97   auto hce = event->GetHCofThisEvent();
 98   auto sdManager = G4SDManager::GetSDMpointer(     98   auto sdManager = G4SDManager::GetSDMpointer();
 99   G4int collId;                                    99   G4int collId;
100                                                   100 
101   // HGCAL EE + FH                                101   // HGCAL EE + FH
102   collId = sdManager->GetCollectionID("Silicon    102   collId = sdManager->GetCollectionID("SiliconPixelHitCollection");
103   auto hc = hce->GetHC(collId);                   103   auto hc = hce->GetHC(collId);
104   if (!hc)                                        104   if (!hc)
105     return;                                       105     return;
106   double esumHGCAL = 0;                           106   double esumHGCAL = 0;
107   double cogzHGCAL = 0;                           107   double cogzHGCAL = 0;
108   int NhitsHGCAL = 0;                             108   int NhitsHGCAL = 0;
109   for (unsigned int i = 0; i < hc->GetSize();     109   for (unsigned int i = 0; i < hc->GetSize(); ++i) {
110     auto hit = static_cast<SiliconPixelHit *>(    110     auto hit = static_cast<SiliconPixelHit *>(hc->GetHit(i));
111     hit->Digitise(fHitTimeCut / CLHEP::ns, fTo    111     hit->Digitise(fHitTimeCut / CLHEP::ns, fToaThreshold / CLHEP::keV);
112                                                   112 
113     if (hit->isValidHit()) {                      113     if (hit->isValidHit()) {
114       fSiHitsID.push_back(hit->ID());             114       fSiHitsID.push_back(hit->ID());
115       fSiHitsX.push_back(hit->GetX());            115       fSiHitsX.push_back(hit->GetX());
116       fSiHitsY.push_back(hit->GetY());            116       fSiHitsY.push_back(hit->GetY());
117       fSiHitsZ.push_back(hit->GetZ());            117       fSiHitsZ.push_back(hit->GetZ());
118       fSiHitsEdep.push_back(hit->GetEdep());      118       fSiHitsEdep.push_back(hit->GetEdep());
119       fSiHitsEdepNonIonising.push_back(hit->Ge    119       fSiHitsEdepNonIonising.push_back(hit->GetEdepNonIonizing());
120       fSiHitsTOA.push_back(hit->GetTOA());        120       fSiHitsTOA.push_back(hit->GetTOA());
121       fSiHitsTOAlast.push_back(hit->GetLastTOA    121       fSiHitsTOAlast.push_back(hit->GetLastTOA());
122       fSiHitsType.push_back(0);                   122       fSiHitsType.push_back(0);
123                                                   123 
124       NhitsHGCAL++;                               124       NhitsHGCAL++;
125       esumHGCAL += hit->GetEdep() * CLHEP::keV    125       esumHGCAL += hit->GetEdep() * CLHEP::keV / CLHEP::MeV;
126       cogzHGCAL += hit->GetZ() * hit->GetEdep(    126       cogzHGCAL += hit->GetZ() * hit->GetEdep();
127     }                                             127     }
128   }                                               128   }
129   if (esumHGCAL > 0)                              129   if (esumHGCAL > 0)
130     cogzHGCAL /= esumHGCAL;                       130     cogzHGCAL /= esumHGCAL;
131                                                   131 
132   analysisManager->FillNtupleDColumn(23, esumH    132   analysisManager->FillNtupleDColumn(23, esumHGCAL / CLHEP::GeV);
133   analysisManager->FillNtupleDColumn(24, cogzH    133   analysisManager->FillNtupleDColumn(24, cogzHGCAL);
134   analysisManager->FillNtupleIColumn(25, Nhits    134   analysisManager->FillNtupleIColumn(25, NhitsHGCAL);
135                                                   135 
136   // AHCAL                                        136   // AHCAL
137   collId = sdManager->GetCollectionID("SiPMHit    137   collId = sdManager->GetCollectionID("SiPMHitCollection");
138   hc = hce->GetHC(collId);                        138   hc = hce->GetHC(collId);
139   if (!hc)                                        139   if (!hc)
140     return;                                       140     return;
141   double esumAHCAL = 0;                           141   double esumAHCAL = 0;
142   double cogzAHCAL = 0;                           142   double cogzAHCAL = 0;
143   int NhitsAHCAL = 0;                             143   int NhitsAHCAL = 0;
144   for (unsigned int i = 0; i < hc->GetSize();     144   for (unsigned int i = 0; i < hc->GetSize(); ++i) {
145     auto hit = static_cast<SiPMHit *>(hc->GetH    145     auto hit = static_cast<SiPMHit *>(hc->GetHit(i));
146     hit->Digitise(-1, 0);                         146     hit->Digitise(-1, 0);
147                                                   147 
148     if (hit->isValidHit()) {                      148     if (hit->isValidHit()) {
149       fSiPMhitsID.push_back(hit->ID());           149       fSiPMhitsID.push_back(hit->ID());
150       fSiPMhitsX.push_back(hit->GetX());          150       fSiPMhitsX.push_back(hit->GetX());
151       fSiPMhitsY.push_back(hit->GetY());          151       fSiPMhitsY.push_back(hit->GetY());
152       fSiPMhitsZ.push_back(hit->GetZ());          152       fSiPMhitsZ.push_back(hit->GetZ());
153       fSiPMhitsEdep.push_back(hit->GetEdep());    153       fSiPMhitsEdep.push_back(hit->GetEdep());
154       fSiPMhitsEdepNonIonising.push_back(hit->    154       fSiPMhitsEdepNonIonising.push_back(hit->GetEdepNonIonizing());
155       fSiPMhitsTOA.push_back(hit->GetTOA());      155       fSiPMhitsTOA.push_back(hit->GetTOA());
156       fSiPMhitsType.push_back(1);                 156       fSiPMhitsType.push_back(1);
157                                                   157 
158       NhitsAHCAL++;                               158       NhitsAHCAL++;
159       esumAHCAL += hit->GetEdep() * CLHEP::keV    159       esumAHCAL += hit->GetEdep() * CLHEP::keV / CLHEP::MeV;
160       cogzAHCAL += hit->GetZ() * hit->GetEdep(    160       cogzAHCAL += hit->GetZ() * hit->GetEdep();
161     }                                             161     }
162   }                                               162   }
163   if (esumAHCAL > 0)                              163   if (esumAHCAL > 0)
164     cogzAHCAL /= esumAHCAL;                       164     cogzAHCAL /= esumAHCAL;
165                                                   165 
166   analysisManager->FillNtupleDColumn(26, esumA    166   analysisManager->FillNtupleDColumn(26, esumAHCAL / CLHEP::GeV);
167   analysisManager->FillNtupleDColumn(27, cogzA    167   analysisManager->FillNtupleDColumn(27, cogzAHCAL);
168   analysisManager->FillNtupleIColumn(28, Nhits    168   analysisManager->FillNtupleIColumn(28, NhitsAHCAL);
169                                                   169 
170   analysisManager->AddNtupleRow();                170   analysisManager->AddNtupleRow();
171 }                                                 171 }
172                                                   172 
173 //....oooOO0OOooo........oooOO0OOooo........oo    173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174                                                   174 
175 void EventAction::DefineCommands() {              175 void EventAction::DefineCommands() {
176                                                   176 
177   fMessenger = new G4GenericMessenger(this, "/    177   fMessenger = new G4GenericMessenger(this, "/HGCalTestbeam/hits/",
178                                       "Primary    178                                       "Primary generator control");
179                                                   179 
180   // time cut command                             180   // time cut command
181   auto &timeCutCmd = fMessenger->DeclareProper    181   auto &timeCutCmd = fMessenger->DeclarePropertyWithUnit(
182       "timeCut", "ns", fHitTimeCut,               182       "timeCut", "ns", fHitTimeCut,
183       "Size of time window for hit digitalisat    183       "Size of time window for hit digitalisation");
184   timeCutCmd.SetParameterName("timeCut", true)    184   timeCutCmd.SetParameterName("timeCut", true);
185   timeCutCmd.SetRange("timeCut>=-1");             185   timeCutCmd.SetRange("timeCut>=-1");
186   timeCutCmd.SetDefaultValue("-1");               186   timeCutCmd.SetDefaultValue("-1");
187                                                   187 
188   // toa threshold command                        188   // toa threshold command
189   auto &toaThresholdCmd = fMessenger->DeclareP    189   auto &toaThresholdCmd = fMessenger->DeclarePropertyWithUnit(
190       "toaThreshold", "keV", fToaThreshold, "T    190       "toaThreshold", "keV", fToaThreshold, "Threshold for TOA activation");
191   toaThresholdCmd.SetParameterName("toaThresho    191   toaThresholdCmd.SetParameterName("toaThreshold", true);
192   toaThresholdCmd.SetRange("toaThreshold>=0");    192   toaThresholdCmd.SetRange("toaThreshold>=0");
193   toaThresholdCmd.SetDefaultValue("0");           193   toaThresholdCmd.SetDefaultValue("0");
194 }                                                 194 }
195                                                   195 
196 //....oooOO0OOooo........oooOO0OOooo........oo    196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
197                                                   197