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