| 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