Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DMXEventAction.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // --------------------------------------------------------------
28 // GEANT 4 - Underground Dark Matter Detector Advanced Example
29 //
30 // For information related to this code contact: Alex Howard
31 // e-mail: alexander.howard@cern.ch
32 // --------------------------------------------------------------
33 // Comments
34 //
35 // Underground Advanced
36 // by A. Howard and H. Araujo
37 // (27th November 2001)
38 //
39 // History/Additions:
40 // 16 Jan 2002 Added analysis
41 //
42 //
43 // EventAction program
44 // --------------------------------------------------------------
45 
46 #include "DMXEventAction.hh"
47 
48 // pass parameters for messengers:
49 #include "DMXRunAction.hh"
50 #include "DMXPrimaryGeneratorAction.hh"
51 
52 // note DMXPmtHit.hh and DMXScintHit.hh are included in DMXEventAction.hh
53 
54 #include "DMXEventActionMessenger.hh"
55 
56 #ifdef G4ANALYSIS_USE
57 #include "DMXAnalysisManager.hh"
58 #endif
59 
60 #include "G4SystemOfUnits.hh"
61 #include "G4Event.hh"
62 #include "G4EventManager.hh"
63 #include "G4HCofThisEvent.hh"
64 #include "G4VHitsCollection.hh"
65 #include "G4TrajectoryContainer.hh"
66 #include "G4Trajectory.hh"
67 #include "G4VVisManager.hh"
68 #include "G4SDManager.hh"
69 #include "G4UImanager.hh"
70 #include "G4UnitsTable.hh"
71 #include "G4ios.hh"
72 #include <fstream>
73 #include <iomanip>
74 
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 DMXEventAction::DMXEventAction(DMXRunAction* DMXRun,
78  DMXPrimaryGeneratorAction* DMXGenerator)
79  : runAct(DMXRun),genAction(DMXGenerator)
80 {
81 
82  // create messenger
83  eventMessenger = new DMXEventActionMessenger(this);
84 
85  // defaults for messenger
86  drawColsFlag = "standard";
87  drawTrksFlag = "all";
88  drawHitsFlag = 1;
89  savePmtFlag = 0;
90  saveHitsFlag = 1;
91 
92  printModulo = 1;
93 
94  // hits collections
95  scintillatorCollID = -1;
96  pmtCollID = -1;
97 
98  energy_pri=0;
99  seeds=NULL;
100 
101 }
102 
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 DMXEventAction::~DMXEventAction() {
106 
107  delete eventMessenger;
108 }
109 
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 void DMXEventAction::BeginOfEventAction(const G4Event* evt) {
113 
114 
115  // grab seeds
116  seeds = genAction->GetEventSeeds();
117 
118  // grab energy of primary
119  energy_pri = genAction->GetEnergyPrimary();
120 
121  event_id = evt->GetEventID();
122 
123  // print this information event by event (modulo n)
124  if (event_id%printModulo == 0)
125  {
126  G4cout << "\n---> Begin of event: " << event_id << G4endl;
127  G4cout << "\n Primary Energy: " << G4BestUnit(energy_pri,"Energy")
128  << G4endl;
129  // HepRandom::showEngineStatus();
130  }
131 
132 
133  // get ID for scintillator hits collection
134  if (scintillatorCollID==-1) {
135  G4SDManager *SDman = G4SDManager::GetSDMpointer();
136  scintillatorCollID = SDman->GetCollectionID("scintillatorCollection");
137  }
138 
139  // get ID for pmt hits collection
140  if (pmtCollID==-1) {
141  G4SDManager *SDman = G4SDManager::GetSDMpointer();
142  pmtCollID = SDman->GetCollectionID("pmtCollection");
143  }
144 
145 }
146 
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149 void DMXEventAction::EndOfEventAction(const G4Event* evt) {
150 
151  // check that both hits collections have been defined
152  if(scintillatorCollID<0||pmtCollID<0) return;
153 
154  // address hits collections
155  DMXScintHitsCollection* SHC = NULL;
156  DMXPmtHitsCollection* PHC = NULL;
157  G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
158  if(HCE) {
159  SHC = (DMXScintHitsCollection*)(HCE->GetHC(scintillatorCollID));
160  PHC = (DMXPmtHitsCollection*)(HCE->GetHC(pmtCollID));
161  }
162 
163  // event summary
164  totEnergy = 0.;
165  totEnergyGammas = 0.;
166  totEnergyNeutrons = 0.;
167  firstParticleE = 0.;
168  particleEnergy = 0.;
169  firstLXeHitTime = 0.;
170  aveTimePmtHits = 0.;
171 
172  firstParticleName = "";
173  particleName = "";
174 
175 
176  // particle flags
177  gamma_ev = false;
178  neutron_ev = false;
179  positron_ev = false;
180  electron_ev = false;
181  proton_ev = false;
182  other_ev = false;
183  start_gamma = false;
184  start_neutron = false;
185 
186 
187  // scintillator hits
188  if(SHC) {
189  S_hits = SHC->entries();
190 
191  for (G4int i=0; i<S_hits; i++) {
192  if(i==0) {
193  firstParticleName = (*SHC)[0]->GetParticle();
194  firstLXeHitTime = (*SHC)[0]->GetTime();
195  firstParticleE = (*SHC)[0]->GetParticleEnergy();
196  if (event_id%printModulo == 0 && S_hits > 0) {
197  G4cout << " First hit in LXe: " << firstParticleName << G4endl;
198  G4cout << " Number of hits in LXe: " << S_hits << G4endl;
199  }
200  }
201  hitEnergy = (*SHC)[i]->GetEdep();
202  totEnergy += hitEnergy;
203 
204  particleName = (*SHC)[i]->GetParticle();
205  particleEnergy = (*SHC)[i]->GetParticleEnergy();
206 
207  if(particleName == "gamma") {
208  gamma_ev = true;
209  start_gamma = true;
210  start_neutron = false;
211  }
212  else if(particleName == "neutron")
213  neutron_ev = true;
214  else if(particleName == "e+")
215  positron_ev = true;
216  else if(particleName == "e-")
217  electron_ev = true;
218  else if(particleName == "proton")
219  proton_ev = true;
220  else {
221  other_ev = true;
222  start_gamma = false;
223  start_neutron = true;
224  }
225 
226  if(start_gamma && !start_neutron)
227  totEnergyGammas += hitEnergy;
228  if(start_neutron && !start_gamma)
229  totEnergyNeutrons += hitEnergy;
230  }
231 
232  if (event_id%printModulo == 0)
233  G4cout << " Total energy in LXe: "
234  << G4BestUnit(totEnergy,"Energy") << G4endl;
235 
236  }
237 
238 
239  // PMT hits
240  if(PHC) {
241  P_hits = PHC->entries();
242 
243  // average time of PMT hits
244  for (G4int i=0; i<P_hits; i++) {
245  G4double time = ( (*PHC)[i]->GetTime() - firstLXeHitTime );
246  aveTimePmtHits += time / (G4double)P_hits;
248  if(P_hits >= 0) {
249 #ifdef G4ANALYSIS_USE
250  // pass first event for histogram record:
251  DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
252  analysis->HistTime(time);
253 #endif
254  }
255  }
256 
257  if (event_id%printModulo == 0 && P_hits > 0) {
258  G4cout << " Average light collection time: "
259  << G4BestUnit(aveTimePmtHits,"Time") << G4endl;
260  G4cout << " Number of PMT hits (photoelectron equivalent): "
261  << P_hits << G4endl;
262  //#ifdef G4ANALYSIS_USE
263  // plot histograms, interactively:
264  //G4bool plotevent=runAct->Getplotevent();
265  //if(plotevent) {
266  // DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
267  // analysis->PlotHistosInter(P_hits);
268  //}
269  //#endif
270  }
271  // write out (x,y,z) of PMT hits
272  if (savePmtFlag)
273  writePmtHitsToFile(PHC);
274  }
275 
276 
277  // write out event summary
278  if(saveHitsFlag)
279  writeScintHitsToFile();
280 
281  // draw trajectories
282  if(drawColsFlag=="standard" && drawTrksFlag!="none")
283  drawTracks(evt);
284 
285  // hits in PMT
286  if(drawHitsFlag && PHC)
287  PHC->DrawAllHits();
288 
289  // print this event by event (modulo n)
290  if (event_id%printModulo == 0)
291  G4cout << "\n---> End of event: " << event_id << G4endl << G4endl;
292 
293 }
294 
295 
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298 void DMXEventAction::writeScintHitsToFile(void) {
299 
300  // G4String filename="hits.out";
301  G4String filename=runAct->GetsavehitsFile();
302  std::ofstream hitsfile(filename, std::ios::app);
303  if(!event_id) {
304  std::ofstream hitsfile(filename);
305  hitsfile <<"Evt Eprim Etot LXe LXeTime PMT PMTTime Seed1 Seed2 First Flags"
306  << G4endl;
307  hitsfile <<"# MeV MeV hits ns hits ns hit"
308  << G4endl
309  << G4endl;
310  }
311 
312  if(S_hits) {
313 
314  if(hitsfile.is_open()) {
315 
316 
317  hitsfile << std::setiosflags(std::ios::fixed)
318  << std::setprecision(4)
319  << std::setiosflags(std::ios::left)
320  << std::setw(6)
321  << event_id << "\t"
322  << energy_pri/MeV << "\t"
323  << totEnergy/MeV << "\t"
324  << S_hits << "\t"
325  << std::setiosflags(std::ios::scientific)
326  << std::setprecision(2)
327  << firstLXeHitTime/nanosecond << "\t"
328  << P_hits << "\t"
329  << std::setiosflags(std::ios::fixed)
330  << std::setprecision(4)
331  << aveTimePmtHits/nanosecond << "\t"
332  << *seeds << "\t"
333  << *(seeds+1) << "\t"
334  << firstParticleName << "\t"
335  << (gamma_ev ? "gamma " : "")
336  << (neutron_ev ? "neutron " : "")
337  << (positron_ev ? "positron " : "")
338  << (electron_ev ? "electron " : "")
339  << (other_ev ? "other " : "")
340  << G4endl;
341 
342  if (event_id%printModulo == 0)
343  G4cout << " Event summary in file " << filename << G4endl;
344  hitsfile.close();
345  }
346 
347 #ifdef G4ANALYSIS_USE
348  long seed1 = *seeds;
349  long seed2 = *(seeds+1);
350  // pass event summary to analysis manager for booking into histos and ntple
351  DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
352  analysis->analyseScintHits(event_id,energy_pri,totEnergy,S_hits,firstLXeHitTime,P_hits,aveTimePmtHits,firstParticleName,firstParticleE,gamma_ev,neutron_ev,positron_ev,electron_ev,other_ev,seed1,seed2);
353 #endif
354  }
355 
356 }
357 
358 
359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
360 void DMXEventAction::writePmtHitsToFile(const DMXPmtHitsCollection* hits) {
361 
362  // G4String filename="pmt.out";
363  G4String filename=runAct->GetsavepmtFile();
364  std::ofstream pmtfile(filename, std::ios::app);
365  G4double x; G4double y; G4double z;
366 
367  if(pmtfile.is_open()) {
368  pmtfile << "Hit# X, mm Y, mm Z, mm" << G4endl;
369  pmtfile << std::setiosflags(std::ios::fixed)
370  << std::setprecision(3)
371  << std::setiosflags(std::ios::left)
372  << std::setw(6);
373  for (G4int i=0; i<P_hits; i++)
374  {
375  x = ((*hits)[i]->GetPos()).x()/mm;
376  y = ((*hits)[i]->GetPos()).y()/mm;
377  z = ((*hits)[i]->GetPos()).z()/mm;
378  pmtfile << i << "\t"
379  << x << "\t"
380  << y << "\t"
381  << z << G4endl;
382 #ifdef G4ANALYSIS_USE
383  // pass pmt hit summary to analysis manager for booking in ntples
384  DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
385  analysis->analysePMTHits(event_id,i,x,y,z);
386 #endif
387  }
388  if (event_id%printModulo == 0 && P_hits > 0)
389  G4cout << " " << P_hits << " PMT hits in " << filename << G4endl;
390  pmtfile.close();
391  }
392 
393 }
394 
395 
396 
397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
398 void DMXEventAction::drawTracks(const G4Event* evt) {
399 
400  if(G4VVisManager::GetConcreteInstance()) {
401  G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers");
402  G4TrajectoryContainer* trajContainer = evt->GetTrajectoryContainer();
403  G4int n_trajectories = 0;
404 
405  if(trajContainer) n_trajectories = trajContainer->entries();
406  for (G4int i=0; i<n_trajectories; i++) {
407  G4Trajectory* trj = (G4Trajectory*)(*trajContainer)[i];
408  if (drawTrksFlag == "all")
409  trj->DrawTrajectory();
410  else if ((drawTrksFlag == "charged") && (trj->GetCharge() != 0.))
411  trj->DrawTrajectory();
412  else if ((drawTrksFlag == "noscint")
413  && (trj->GetParticleName() != "opticalphoton"))
414  trj->DrawTrajectory();
415  }
416 
417  // G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
418  }
419 
420 }
Generated on Sat May 25 2013 14:32:15 for Geant4 by   doxygen 1.8.4