Geant4 Cross Reference |
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 // << 27 /// \file optical/LXe/src/LXeEventAction.cc 26 /// \file optical/LXe/src/LXeEventAction.cc 28 /// \brief Implementation of the LXeEventActio 27 /// \brief Implementation of the LXeEventAction class 29 // 28 // 30 // 29 // 31 #include "LXeEventAction.hh" 30 #include "LXeEventAction.hh" 32 << 33 #include "LXeDetectorConstruction.hh" << 34 #include "LXeHistoManager.hh" << 35 #include "LXePMTHit.hh" << 36 #include "LXeRun.hh" << 37 #include "LXeScintHit.hh" 31 #include "LXeScintHit.hh" >> 32 #include "LXePMTHit.hh" >> 33 #include "LXeUserEventInformation.hh" 38 #include "LXeTrajectory.hh" 34 #include "LXeTrajectory.hh" >> 35 #include "LXeRecorderBase.hh" 39 36 40 #include "G4Event.hh" << 41 #include "G4EventManager.hh" 37 #include "G4EventManager.hh" 42 #include "G4RunManager.hh" << 43 #include "G4SDManager.hh" 38 #include "G4SDManager.hh" 44 #include "G4SystemOfUnits.hh" << 39 #include "G4RunManager.hh" 45 #include "G4Trajectory.hh" << 40 #include "G4Event.hh" >> 41 #include "G4EventManager.hh" 46 #include "G4TrajectoryContainer.hh" 42 #include "G4TrajectoryContainer.hh" 47 #include "G4UImanager.hh" << 43 #include "G4Trajectory.hh" 48 #include "G4VVisManager.hh" 44 #include "G4VVisManager.hh" 49 #include "G4ios.hh" 45 #include "G4ios.hh" >> 46 #include "G4UImanager.hh" >> 47 #include "G4SystemOfUnits.hh" 50 48 51 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 52 50 53 LXeEventAction::LXeEventAction(const LXeDetect << 51 LXeEventAction::LXeEventAction(LXeRecorderBase* r) >> 52 : fRecorder(r),fSaveThreshold(0),fScintCollID(-1),fPMTCollID(-1),fVerbose(0), >> 53 fPMTThreshold(1),fForcedrawphotons(false),fForcenophotons(false) 54 { 54 { 55 fEventMessenger = new LXeEventMessenger(this 55 fEventMessenger = new LXeEventMessenger(this); 56 } 56 } 57 << 57 58 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 59 59 60 LXeEventAction::~LXeEventAction() << 60 LXeEventAction::~LXeEventAction(){} 61 { << 62 delete fEventMessenger; << 63 } << 64 61 65 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 66 63 67 void LXeEventAction::BeginOfEventAction(const << 64 void LXeEventAction::BeginOfEventAction(const G4Event* anEvent){ 68 { << 65 69 fHitCount = 0; << 66 //New event, add the user information object 70 fPhotonCount_Scint = 0; << 67 G4EventManager:: 71 fPhotonCount_Ceren = 0; << 68 GetEventManager()->SetUserInformation(new LXeUserEventInformation); 72 fAbsorptionCount = 0; << 73 fBoundaryAbsorptionCount = 0; << 74 fTotE = 0.0; << 75 << 76 fConvPosSet = false; << 77 fEdepMax = 0.0; << 78 << 79 fPMTsAboveThreshold = 0; << 80 69 81 G4SDManager* SDman = G4SDManager::GetSDMpoin 70 G4SDManager* SDman = G4SDManager::GetSDMpointer(); 82 if (fScintCollID < 0) fScintCollID = SDman-> << 71 if(fScintCollID<0) 83 if (fPMTCollID < 0) fPMTCollID = SDman->GetC << 72 fScintCollID=SDman->GetCollectionID("scintCollection"); 84 } << 73 if(fPMTCollID<0) >> 74 fPMTCollID=SDman->GetCollectionID("pmtHitCollection"); 85 75 >> 76 if(fRecorder)fRecorder->RecordBeginOfEvent(anEvent); >> 77 } >> 78 86 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 87 80 88 void LXeEventAction::EndOfEventAction(const G4 << 81 void LXeEventAction::EndOfEventAction(const G4Event* anEvent){ 89 { << 90 G4TrajectoryContainer* trajectoryContainer = << 91 82 >> 83 LXeUserEventInformation* eventInformation >> 84 =(LXeUserEventInformation*)anEvent->GetUserInformation(); >> 85 >> 86 G4TrajectoryContainer* trajectoryContainer=anEvent->GetTrajectoryContainer(); >> 87 92 G4int n_trajectories = 0; 88 G4int n_trajectories = 0; 93 if (trajectoryContainer) n_trajectories = tr 89 if (trajectoryContainer) n_trajectories = trajectoryContainer->entries(); 94 90 95 // extract the trajectories and draw them 91 // extract the trajectories and draw them 96 if (G4VVisManager::GetConcreteInstance()) { << 92 if (G4VVisManager::GetConcreteInstance()){ 97 for (G4int i = 0; i < n_trajectories; ++i) << 93 for (G4int i=0; i<n_trajectories; i++){ 98 auto trj = (LXeTrajectory*)((*(anEvent-> << 94 LXeTrajectory* trj = (LXeTrajectory*) 99 if (trj->GetParticleName() == "opticalph << 95 ((*(anEvent->GetTrajectoryContainer()))[i]); >> 96 if(trj->GetParticleName()=="opticalphoton"){ 100 trj->SetForceDrawTrajectory(fForcedraw 97 trj->SetForceDrawTrajectory(fForcedrawphotons); 101 trj->SetForceNoDrawTrajectory(fForceno 98 trj->SetForceNoDrawTrajectory(fForcenophotons); 102 } 99 } 103 trj->DrawTrajectory(); << 100 trj->DrawTrajectory(50); 104 } 101 } 105 } 102 } 106 << 103 107 LXeScintHitsCollection* scintHC = nullptr; << 104 LXeScintHitsCollection* scintHC = 0; 108 LXePMTHitsCollection* pmtHC = nullptr; << 105 LXePMTHitsCollection* pmtHC = 0; 109 G4HCofThisEvent* hitsCE = anEvent->GetHCofTh 106 G4HCofThisEvent* hitsCE = anEvent->GetHCofThisEvent(); 110 << 107 111 // Get the hit collections << 108 //Get the hit collections 112 if (hitsCE) { << 109 if(hitsCE){ 113 if (fScintCollID >= 0) { << 110 if(fScintCollID>=0)scintHC = (LXeScintHitsCollection*)(hitsCE->GetHC(fScintCollID)); 114 scintHC = (LXeScintHitsCollection*)(hits << 111 if(fPMTCollID>=0)pmtHC = (LXePMTHitsCollection*)(hitsCE->GetHC(fPMTCollID)); 115 } << 116 if (fPMTCollID >= 0) { << 117 pmtHC = (LXePMTHitsCollection*)(hitsCE-> << 118 } << 119 } 112 } 120 113 121 // Hits in scintillator << 114 //Hits in scintillator 122 if (scintHC) { << 115 if(scintHC){ 123 size_t n_hit = scintHC->entries(); << 116 int n_hit = scintHC->entries(); 124 G4ThreeVector eWeightPos(0.); << 117 G4ThreeVector eWeightPos(0.); 125 G4double edep; 118 G4double edep; 126 G4double edepMax = 0; << 119 G4double edepMax=0; 127 120 128 for (size_t i = 0; i < n_hit; ++i) { // g << 121 for(int i=0;i<n_hit;i++){ //gather info on hits in scintillator 129 edep = (*scintHC)[i]->GetEdep(); << 122 edep=(*scintHC)[i]->GetEdep(); 130 fTotE += edep; << 123 eventInformation->IncEDep(edep); //sum up the edep 131 eWeightPos += (*scintHC)[i]->GetPos() * << 124 eWeightPos += (*scintHC)[i]->GetPos()*edep;//calculate energy weighted pos 132 if (edep > edepMax) { << 125 if(edep>edepMax){ 133 edepMax = edep; // store max energy d << 126 edepMax=edep;//store max energy deposit 134 G4ThreeVector posMax = (*scintHC)[i]-> << 127 G4ThreeVector posMax=(*scintHC)[i]->GetPos(); 135 fPosMax = posMax; << 128 eventInformation->SetPosMax(posMax,edep); 136 fEdepMax = edep; << 129 } 137 } << 130 } 138 } << 131 if(eventInformation->GetEDep()==0.){ 139 << 132 if(fVerbose>0)G4cout<<"No hits in the scintillator this event."<<G4endl; 140 G4AnalysisManager::Instance()->FillH1(7, f << 133 } 141 << 134 else{ 142 if (fTotE == 0.) { << 135 //Finish calculation of energy weighted position 143 if (fVerbose > 0) G4cout << "No hits in << 136 eWeightPos/=eventInformation->GetEDep(); 144 } << 137 eventInformation->SetEWeightPos(eWeightPos); 145 else { << 138 if(fVerbose>0){ 146 // Finish calculation of energy weighted << 139 G4cout << "\tEnergy weighted position of hits in LXe : " 147 eWeightPos /= fTotE; << 140 << eWeightPos/mm << G4endl; 148 fEWeightPos = eWeightPos; << 141 } 149 if (fVerbose > 0) { << 142 } 150 G4cout << "\tEnergy weighted position << 143 if(fVerbose>0){ >> 144 G4cout << "\tTotal energy deposition in scintillator : " >> 145 << eventInformation->GetEDep() / keV << " (keV)" << G4endl; >> 146 } >> 147 } >> 148 >> 149 if(pmtHC){ >> 150 G4ThreeVector reconPos(0.,0.,0.); >> 151 G4int pmts=pmtHC->entries(); >> 152 //Gather info from all PMTs >> 153 for(G4int i=0;i<pmts;i++){ >> 154 eventInformation->IncHitCount((*pmtHC)[i]->GetPhotonCount()); >> 155 reconPos+=(*pmtHC)[i]->GetPMTPos()*(*pmtHC)[i]->GetPhotonCount(); >> 156 if((*pmtHC)[i]->GetPhotonCount()>=fPMTThreshold){ >> 157 eventInformation->IncPMTSAboveThreshold(); 151 } 158 } 152 } << 159 else{//wasnt above the threshold, turn it back off 153 if (fVerbose > 0) { << 154 G4cout << "\tTotal energy deposition in << 155 } << 156 } << 157 << 158 if (pmtHC) { << 159 G4ThreeVector reconPos(0., 0., 0.); << 160 size_t pmts = pmtHC->entries(); << 161 // Gather info from all PMTs << 162 for (size_t i = 0; i < pmts; ++i) { << 163 fHitCount += (*pmtHC)[i]->GetPhotonCount << 164 reconPos += (*pmtHC)[i]->GetPMTPos() * ( << 165 if ((*pmtHC)[i]->GetPhotonCount() >= fPM << 166 ++fPMTsAboveThreshold; << 167 } << 168 else { // wasn't above the threshold, t << 169 (*pmtHC)[i]->SetDrawit(false); 160 (*pmtHC)[i]->SetDrawit(false); 170 } 161 } 171 } 162 } 172 << 163 173 G4AnalysisManager::Instance()->FillH1(1, f << 164 if(eventInformation->GetHitCount()>0){//dont bother unless there were hits 174 G4AnalysisManager::Instance()->FillH1(2, f << 165 reconPos/=eventInformation->GetHitCount(); 175 << 166 if(fVerbose>0){ 176 if (fHitCount > 0) { // don't bother unle << 167 G4cout << "\tReconstructed position of hits in LXe : " 177 reconPos /= fHitCount; << 168 << reconPos/mm << G4endl; 178 if (fVerbose > 0) { << 179 G4cout << "\tReconstructed position of << 180 } 169 } 181 fReconPos = reconPos; << 170 eventInformation->SetReconPos(reconPos); 182 } 171 } 183 pmtHC->DrawAllHits(); 172 pmtHC->DrawAllHits(); 184 } 173 } 185 174 186 G4AnalysisManager::Instance()->FillH1(3, fPh << 175 if(fVerbose>0){ 187 G4AnalysisManager::Instance()->FillH1(4, fPh << 176 //End of event output. later to be controlled by a verbose level 188 G4AnalysisManager::Instance()->FillH1(5, fAb << 177 G4cout << "\tNumber of photons that hit PMTs in this event : " 189 G4AnalysisManager::Instance()->FillH1(6, fBo << 178 << eventInformation->GetHitCount() << G4endl; 190 << 179 G4cout << "\tNumber of PMTs above threshold("<<fPMTThreshold<<") : " 191 if (fVerbose > 0) { << 180 << eventInformation->GetPMTSAboveThreshold() << G4endl; 192 // End of event output. later to be contro << 181 G4cout << "\tNumber of photons produced by scintillation in this event : " 193 G4cout << "\tNumber of photons that hit PM << 182 << eventInformation->GetPhotonCount_Scint() << G4endl; 194 G4cout << "\tNumber of PMTs above threshol << 183 G4cout << "\tNumber of photons produced by cerenkov in this event : " 195 << G4endl; << 184 << eventInformation->GetPhotonCount_Ceren() << G4endl; 196 G4cout << "\tNumber of photons produced by << 185 G4cout << "\tNumber of photons absorbed (OpAbsorption) in this event : " 197 << G4endl; << 186 << eventInformation->GetAbsorptionCount() << G4endl; 198 G4cout << "\tNumber of photons produced by << 199 << G4endl; << 200 G4cout << "\tNumber of photons absorbed (O << 201 << G4endl; << 202 G4cout << "\tNumber of photons absorbed at 187 G4cout << "\tNumber of photons absorbed at boundaries (OpBoundary) in " 203 << "this event : " << fBoundaryAbso << 188 << "this event : " << eventInformation->GetBoundaryAbsorptionCount() 204 G4cout << "Unaccounted for photons in this << 189 << G4endl; 205 << (fPhotonCount_Scint + fPhotonCou << 190 G4cout << "Unacounted for photons in this event : " 206 - fBoundaryAbsorptionCount) << 191 << (eventInformation->GetPhotonCount_Scint() + >> 192 eventInformation->GetPhotonCount_Ceren() - >> 193 eventInformation->GetAbsorptionCount() - >> 194 eventInformation->GetHitCount() - >> 195 eventInformation->GetBoundaryAbsorptionCount()) 207 << G4endl; 196 << G4endl; 208 } 197 } 209 << 198 //If we have set the flag to save 'special' events, save here 210 // update the run statistics << 199 if(fSaveThreshold&&eventInformation->GetPhotonCount() <= fSaveThreshold) 211 auto run = static_cast<LXeRun*>(G4RunManager << 212 << 213 run->IncHitCount(fHitCount); << 214 run->IncPhotonCount_Scint(fPhotonCount_Scint << 215 run->IncPhotonCount_Ceren(fPhotonCount_Ceren << 216 run->IncEDep(fTotE); << 217 run->IncAbsorption(fAbsorptionCount); << 218 run->IncBoundaryAbsorption(fBoundaryAbsorpti << 219 run->IncHitsAboveThreshold(fPMTsAboveThresho << 220 << 221 // If we have set the flag to save 'special' << 222 if (fPhotonCount_Scint + fPhotonCount_Ceren << 223 G4RunManager::GetRunManager()->rndmSaveThi 200 G4RunManager::GetRunManager()->rndmSaveThisEvent(); 224 } << 201 >> 202 if(fRecorder)fRecorder->RecordEndOfEvent(anEvent); 225 } 203 } 226 204 227 //....oooOO0OOooo........oooOO0OOooo........oo 205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 206 >> 207 void LXeEventAction::SetSaveThreshold(G4int save){ >> 208 /*Sets the save threshold for the random number seed. If the number of photons >> 209 generated in an event is lower than this, then save the seed for this event >> 210 in a file called run###evt###.rndm >> 211 */ >> 212 fSaveThreshold=save; >> 213 G4RunManager::GetRunManager()->SetRandomNumberStore(true); >> 214 G4RunManager::GetRunManager()->SetRandomNumberStoreDir("random/"); >> 215 // G4UImanager::GetUIpointer()->ApplyCommand("/random/setSavingFlag 1"); >> 216 } 228 217