Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/optical/LXe/src/LXeEventAction.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/extended/optical/LXe/src/LXeEventAction.cc (Version 11.3.0) and /examples/extended/optical/LXe/src/LXeEventAction.cc (Version 9.6.p1)


  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