source: trunk/examples/advanced/underground_physics/src/DMXEventAction.cc@ 1303

Last change on this file since 1303 was 807, checked in by garnier, 17 years ago

update
File size: 12.5 KB
RevLine 
[807]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 "G4Event.hh"
61#include "G4EventManager.hh"
62#include "G4HCofThisEvent.hh"
63#include "G4VHitsCollection.hh"
64#include "G4TrajectoryContainer.hh"
65#include "G4Trajectory.hh"
66#include "G4VVisManager.hh"
67#include "G4SDManager.hh"
68#include "G4UImanager.hh"
69#include "G4UnitsTable.hh"
70#include "G4ios.hh"
71#include <fstream>
72#include <iomanip>
73
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76DMXEventAction::DMXEventAction(DMXRunAction* DMXRun,
77 DMXPrimaryGeneratorAction* DMXGenerator)
78 : runAct(DMXRun),genAction(DMXGenerator)
79{
80
81 // create messenger
82 eventMessenger = new DMXEventActionMessenger(this);
83
84 // defaults for messenger
85 drawColsFlag = "standard";
86 drawTrksFlag = "all";
87 drawHitsFlag = 1;
88 savePmtFlag = 0;
89 saveHitsFlag = 1;
90
91 printModulo = 1;
92
93 // hits collections
94 scintillatorCollID = -1;
95 pmtCollID = -1;
96
97 energy_pri=0;
98 seeds=NULL;
99
100}
101
102
103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104DMXEventAction::~DMXEventAction() {
105
106 delete eventMessenger;
107}
108
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111void DMXEventAction::BeginOfEventAction(const G4Event* evt) {
112
113
114 // grab seeds
115 seeds = genAction->GetEventSeeds();
116
117 // grab energy of primary
118 energy_pri = genAction->GetEnergyPrimary();
119
120 event_id = evt->GetEventID();
121
122 // print this information event by event (modulo n)
123 if (event_id%printModulo == 0)
124 {
125 G4cout << "\n---> Begin of event: " << event_id << G4endl;
126 G4cout << "\n Primary Energy: " << G4BestUnit(energy_pri,"Energy")
127 << G4endl;
128 // HepRandom::showEngineStatus();
129 }
130
131
132 // get ID for scintillator hits collection
133 if (scintillatorCollID==-1) {
134 G4SDManager *SDman = G4SDManager::GetSDMpointer();
135 scintillatorCollID = SDman->GetCollectionID("scintillatorCollection");
136 }
137
138 // get ID for pmt hits collection
139 if (pmtCollID==-1) {
140 G4SDManager *SDman = G4SDManager::GetSDMpointer();
141 pmtCollID = SDman->GetCollectionID("pmtCollection");
142 }
143
144}
145
146
147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148void DMXEventAction::EndOfEventAction(const G4Event* evt) {
149
150 // check that both hits collections have been defined
151 if(scintillatorCollID<0||pmtCollID<0) return;
152
153 // address hits collections
154 DMXScintHitsCollection* SHC = NULL;
155 DMXPmtHitsCollection* PHC = NULL;
156 G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
157 if(HCE) {
158 SHC = (DMXScintHitsCollection*)(HCE->GetHC(scintillatorCollID));
159 PHC = (DMXPmtHitsCollection*)(HCE->GetHC(pmtCollID));
160 }
161
162 // event summary
163 totEnergy = 0.;
164 totEnergyGammas = 0.;
165 totEnergyNeutrons = 0.;
166 firstParticleE = 0.;
167 particleEnergy = 0.;
168 firstLXeHitTime = 0.;
169 aveTimePmtHits = 0.;
170
171 firstParticleName = "";
172 particleName = "";
173
174
175 // particle flags
176 gamma_ev = false;
177 neutron_ev = false;
178 positron_ev = false;
179 electron_ev = false;
180 proton_ev = false;
181 other_ev = false;
182 start_gamma = false;
183 start_neutron = false;
184
185
186 // scintillator hits
187 if(SHC) {
188 S_hits = SHC->entries();
189
190 for (G4int i=0; i<S_hits; i++) {
191 if(i==0) {
192 firstParticleName = (*SHC)[0]->GetParticle();
193 firstLXeHitTime = (*SHC)[0]->GetTime();
194 firstParticleE = (*SHC)[0]->GetParticleEnergy();
195 if (event_id%printModulo == 0 && S_hits > 0) {
196 G4cout << " First hit in LXe: " << firstParticleName << G4endl;
197 G4cout << " Number of hits in LXe: " << S_hits << G4endl;
198 }
199 }
200 hitEnergy = (*SHC)[i]->GetEdep();
201 totEnergy += hitEnergy;
202
203 particleName = (*SHC)[i]->GetParticle();
204 particleEnergy = (*SHC)[i]->GetParticleEnergy();
205
206 if(particleName == "gamma") {
207 gamma_ev = true;
208 start_gamma = true;
209 start_neutron = false;
210 }
211 else if(particleName == "neutron")
212 neutron_ev = true;
213 else if(particleName == "e+")
214 positron_ev = true;
215 else if(particleName == "e-")
216 electron_ev = true;
217 else if(particleName == "proton")
218 proton_ev = true;
219 else {
220 other_ev = true;
221 start_gamma = false;
222 start_neutron = true;
223 }
224
225 if(start_gamma && !start_neutron)
226 totEnergyGammas += hitEnergy;
227 if(start_neutron && !start_gamma)
228 totEnergyNeutrons += hitEnergy;
229 }
230
231 if (event_id%printModulo == 0)
232 G4cout << " Total energy in LXe: "
233 << G4BestUnit(totEnergy,"Energy") << G4endl;
234
235 }
236
237
238 // PMT hits
239 if(PHC) {
240 P_hits = PHC->entries();
241
242 // average time of PMT hits
243 for (G4int i=0; i<P_hits; i++) {
244 G4double time = ( (*PHC)[i]->GetTime() - firstLXeHitTime );
245 aveTimePmtHits += time / (G4double)P_hits;
246 //// if (event_id == 0) {
247 if(P_hits >= 0) {
248#ifdef G4ANALYSIS_USE
249 // pass first event for histogram record:
250 DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
251 analysis->HistTime(time);
252#endif
253 }
254 }
255
256 if (event_id%printModulo == 0 && P_hits > 0) {
257 G4cout << " Average light collection time: "
258 << G4BestUnit(aveTimePmtHits,"Time") << G4endl;
259 G4cout << " Number of PMT hits (photoelectron equivalent): "
260 << P_hits << G4endl;
261 //#ifdef G4ANALYSIS_USE
262 // plot histograms, interactively:
263 //G4bool plotevent=runAct->Getplotevent();
264 //if(plotevent) {
265 // DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
266 // analysis->PlotHistosInter(P_hits);
267 //}
268 //#endif
269 }
270 // write out (x,y,z) of PMT hits
271 if (savePmtFlag)
272 writePmtHitsToFile(PHC);
273 }
274
275
276 // write out event summary
277 if(saveHitsFlag)
278 writeScintHitsToFile();
279
280 // draw trajectories
281 if(drawColsFlag=="standard" && drawTrksFlag!="none")
282 drawTracks(evt);
283
284 // hits in PMT
285 if(drawHitsFlag && PHC)
286 PHC->DrawAllHits();
287
288 // print this event by event (modulo n)
289 if (event_id%printModulo == 0)
290 G4cout << "\n---> End of event: " << event_id << G4endl << G4endl;
291
292}
293
294
295
296//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
297void DMXEventAction::writeScintHitsToFile(void) {
298
299 // G4String filename="hits.out";
300 G4String filename=runAct->GetsavehitsFile();
301 std::ofstream hitsfile(filename, std::ios::app);
302 if(!event_id) {
303 std::ofstream hitsfile(filename);
304 hitsfile <<"Evt Eprim Etot LXe LXeTime PMT PMTTime Seed1 Seed2 First Flags"
305 << G4endl;
306 hitsfile <<"# MeV MeV hits ns hits ns hit"
307 << G4endl
308 << G4endl;
309 }
310
311 if(S_hits) {
312
313 if(hitsfile.is_open()) {
314
315
316 hitsfile << std::setiosflags(std::ios::fixed)
317 << std::setprecision(4)
318 << std::setiosflags(std::ios::left)
319 << std::setw(6)
320 << event_id << "\t"
321 << energy_pri/MeV << "\t"
322 << totEnergy/MeV << "\t"
323 << S_hits << "\t"
324 << std::setiosflags(std::ios::scientific)
325 << std::setprecision(2)
326 << firstLXeHitTime/nanosecond << "\t"
327 << P_hits << "\t"
328 << std::setiosflags(std::ios::fixed)
329 << std::setprecision(4)
330 << aveTimePmtHits/nanosecond << "\t"
331 << *seeds << "\t"
332 << *(seeds+1) << "\t"
333 << firstParticleName << "\t"
334 << (gamma_ev ? "gamma " : "")
335 << (neutron_ev ? "neutron " : "")
336 << (positron_ev ? "positron " : "")
337 << (electron_ev ? "electron " : "")
338 << (other_ev ? "other " : "")
339 << G4endl;
340
341 if (event_id%printModulo == 0)
342 G4cout << " Event summary in file " << filename << G4endl;
343 hitsfile.close();
344 }
345
346#ifdef G4ANALYSIS_USE
347 long seed1 = *seeds;
348 long seed2 = *(seeds+1);
349 // pass event summary to analysis manager for booking into histos and ntple
350 DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
351 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);
352#endif
353 }
354
355}
356
357
358//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359void DMXEventAction::writePmtHitsToFile(const DMXPmtHitsCollection* hits) {
360
361 // G4String filename="pmt.out";
362 G4String filename=runAct->GetsavepmtFile();
363 std::ofstream pmtfile(filename, std::ios::app);
364 G4double x; G4double y; G4double z;
365
366 if(pmtfile.is_open()) {
367 pmtfile << "Hit# X, mm Y, mm Z, mm" << G4endl;
368 pmtfile << std::setiosflags(std::ios::fixed)
369 << std::setprecision(3)
370 << std::setiosflags(std::ios::left)
371 << std::setw(6);
372 for (G4int i=0; i<P_hits; i++)
373 {
374 x = ((*hits)[i]->GetPos()).x()/mm;
375 y = ((*hits)[i]->GetPos()).y()/mm;
376 z = ((*hits)[i]->GetPos()).z()/mm;
377 pmtfile << i << "\t"
378 << x << "\t"
379 << y << "\t"
380 << z << G4endl;
381#ifdef G4ANALYSIS_USE
382 // pass pmt hit summary to analysis manager for booking in ntples
383 DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
384 analysis->analysePMTHits(event_id,i,x,y,z);
385#endif
386 }
387 if (event_id%printModulo == 0 && P_hits > 0)
388 G4cout << " " << P_hits << " PMT hits in " << filename << G4endl;
389 pmtfile.close();
390 }
391
392}
393
394
395
396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
397void DMXEventAction::drawTracks(const G4Event* evt) {
398
399 if(G4VVisManager::GetConcreteInstance()) {
400 G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers");
401 G4TrajectoryContainer* trajContainer = evt->GetTrajectoryContainer();
402 G4int n_trajectories = 0;
403
404 if(trajContainer) n_trajectories = trajContainer->entries();
405 for (G4int i=0; i<n_trajectories; i++) {
406 G4Trajectory* trj = (G4Trajectory*)(*trajContainer)[i];
407 if (drawTrksFlag == "all")
408 trj->DrawTrajectory();
409 else if ((drawTrksFlag == "charged") && (trj->GetCharge() != 0.))
410 trj->DrawTrajectory();
411 else if ((drawTrksFlag == "noscint")
412 && (trj->GetParticleName() != "opticalphoton"))
413 trj->DrawTrajectory();
414 }
415
416 // G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
417 }
418
419}
Note: See TracBrowser for help on using the repository browser.