Geant4  10.00.p01
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"
51 
52 // note DMXPmtHit.hh and DMXScintHit.hh are included in DMXEventAction.hh
53 
55 #include "DMXAnalysisManager.hh"
56 
57 #include "G4SystemOfUnits.hh"
58 #include "G4Event.hh"
59 #include "G4EventManager.hh"
60 #include "G4HCofThisEvent.hh"
61 #include "G4VHitsCollection.hh"
62 #include "G4TrajectoryContainer.hh"
63 #include "G4Trajectory.hh"
64 #include "G4VVisManager.hh"
65 #include "G4SDManager.hh"
66 #include "G4UImanager.hh"
67 #include "G4UnitsTable.hh"
68 #include "G4ios.hh"
69 #include <fstream>
70 #include <iomanip>
71 
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75  DMXPrimaryGeneratorAction* DMXGenerator)
76  : runAct(DMXRun),genAction(DMXGenerator)
77 {
78 
79  // create messenger
81 
82  // defaults for messenger
83  drawColsFlag = "standard";
84  drawTrksFlag = "all";
85  drawHitsFlag = 1;
86  savePmtFlag = 0;
87  saveHitsFlag = 1;
88 
89  printModulo = 1;
90 
91  // hits collections
92  scintillatorCollID = -1;
93  pmtCollID = -1;
94 
95  energy_pri=0;
96  seeds=NULL;
97 
98 }
99 
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103 
104  delete eventMessenger;
105 }
106 
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
111 
112  // grab seeds
114 
115  // grab energy of primary
117 
118  event_id = evt->GetEventID();
119 
120  // print this information event by event (modulo n)
121  if (event_id%printModulo == 0)
122  {
123  G4cout << "\n---> Begin of event: " << event_id << G4endl;
124  G4cout << "\n Primary Energy: " << G4BestUnit(energy_pri,"Energy")
125  << G4endl;
126  // HepRandom::showEngineStatus();
127  }
128 
129 
130  // get ID for scintillator hits collection
131  if (scintillatorCollID==-1) {
133  scintillatorCollID = SDman->GetCollectionID("scintillatorCollection");
134  }
135 
136  // get ID for pmt hits collection
137  if (pmtCollID==-1) {
139  pmtCollID = SDman->GetCollectionID("pmtCollection");
140  }
141 
142 }
143 
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147 
148  // check that both hits collections have been defined
149  if(scintillatorCollID<0||pmtCollID<0) return;
150 
151  G4AnalysisManager* man = G4AnalysisManager::Instance();
152 
153  // address hits collections
154  DMXScintHitsCollection* SHC = NULL;
155  DMXPmtHitsCollection* PHC = NULL;
156  G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
157  if(HCE) {
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)
227  if(start_neutron && !start_gamma)
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;
247  if(P_hits >= 0) {
248  man->FillH1(7,time);
249  }
250  }
251 
252  if (event_id%printModulo == 0 && P_hits > 0) {
253  G4cout << " Average light collection time: "
254  << G4BestUnit(aveTimePmtHits,"Time") << G4endl;
255  G4cout << " Number of PMT hits (photoelectron equivalent): "
256  << P_hits << G4endl;
257  }
258  // write out (x,y,z) of PMT hits
259  if (savePmtFlag)
260  writePmtHitsToFile(PHC);
261  }
262 
263 
264  // write out event summary
265  if(saveHitsFlag)
267 
268  // draw trajectories
269  if(drawColsFlag=="standard" && drawTrksFlag!="none")
270  drawTracks(evt);
271 
272  // hits in PMT
273  if(drawHitsFlag && PHC)
274  PHC->DrawAllHits();
275 
276  // print this event by event (modulo n)
277  if (event_id%printModulo == 0)
278  G4cout << "\n---> End of event: " << event_id << G4endl << G4endl;
279 
280 }
281 
282 
283 
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286 
287  // G4String filename="hits.out";
288  G4String filename=runAct->GetsavehitsFile();
289  std::ofstream hitsfile(filename, std::ios::app);
290  if(!event_id) {
291  std::ofstream hitsfileInit(filename);
292  hitsfileInit <<"Evt Eprim Etot LXe LXeTime PMT PMTTime Seed1 Seed2 First Flags"
293  << G4endl;
294  hitsfileInit <<"# MeV MeV hits ns hits ns hit"
295  << G4endl
296  << G4endl;
297  }
298 
299  if(S_hits) {
300 
301  if(hitsfile.is_open()) {
302 
303 
304  hitsfile << std::setiosflags(std::ios::fixed)
305  << std::setprecision(4)
306  << std::setiosflags(std::ios::left)
307  << std::setw(6)
308  << event_id << "\t"
309  << energy_pri/MeV << "\t"
310  << totEnergy/MeV << "\t"
311  << S_hits << "\t"
312  << std::setiosflags(std::ios::scientific)
313  << std::setprecision(2)
314  << firstLXeHitTime/nanosecond << "\t"
315  << P_hits << "\t"
316  << std::setiosflags(std::ios::fixed)
317  << std::setprecision(4)
318  << aveTimePmtHits/nanosecond << "\t"
319  << *seeds << "\t"
320  << *(seeds+1) << "\t"
321  << firstParticleName << "\t"
322  << (gamma_ev ? "gamma " : "")
323  << (neutron_ev ? "neutron " : "")
324  << (positron_ev ? "positron " : "")
325  << (electron_ev ? "electron " : "")
326  << (other_ev ? "other " : "")
327  << G4endl;
328 
329  if (event_id%printModulo == 0)
330  G4cout << " Event summary in file " << filename << G4endl;
331  hitsfile.close();
332  }
333 
334  G4AnalysisManager* man = G4AnalysisManager::Instance();
335  G4int firstparticleIndex = 0;
336  if(firstParticleName == "gamma") firstparticleIndex = 1;
337  else if (firstParticleName == "neutron") firstparticleIndex = 2;
338  else if(firstParticleName == "electron") firstparticleIndex = 3;
339  else if(firstParticleName == "positron") firstparticleIndex = 4;
340  else{
341  firstparticleIndex = 5;
342  man->FillH1(3,totEnergy);
343  }
344 
345  man->FillH1(4,P_hits,10); //weight
346  man->FillH1(5,P_hits);
347 
348  man->FillH1(1,energy_pri/keV);
349  man->FillH1(2,totEnergy/keV);
350  man->FillH1(6,aveTimePmtHits/ns);
351 
352  long seed1 = *seeds;
353  long seed2 = *(seeds+1);
354 
355  //Fill ntuple #2
356  man->FillNtupleDColumn(2,0,event_id);
357  man->FillNtupleDColumn(2,1,energy_pri/keV);
358  man->FillNtupleDColumn(2,2,totEnergy);
359  man->FillNtupleDColumn(2,3,S_hits);
361  man->FillNtupleDColumn(2,5,P_hits);
363  man->FillNtupleDColumn(2,7,firstparticleIndex);
365  man->FillNtupleDColumn(2,9,gamma_ev);
366  man->FillNtupleDColumn(2,10,neutron_ev);
367  man->FillNtupleDColumn(2,11,positron_ev);
368  man->FillNtupleDColumn(2,12,electron_ev);
369  man->FillNtupleDColumn(2,13,other_ev);
370  man->FillNtupleDColumn(2,14,seed1);
371  man->FillNtupleDColumn(2,15,seed2);
372  man->AddNtupleRow(2);
373 
374  }
375 
376 }
377 
378 
379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381 
382  // G4String filename="pmt.out";
383  G4String filename=runAct->GetsavepmtFile();
384  std::ofstream pmtfile(filename, std::ios::app);
385  G4double x; G4double y; G4double z;
386  G4AnalysisManager* man = G4AnalysisManager::Instance();
387 
388  if(pmtfile.is_open()) {
389  pmtfile << "Hit# X, mm Y, mm Z, mm" << G4endl;
390  pmtfile << std::setiosflags(std::ios::fixed)
391  << std::setprecision(3)
392  << std::setiosflags(std::ios::left)
393  << std::setw(6);
394  for (G4int i=0; i<P_hits; i++)
395  {
396  x = ((*hits)[i]->GetPos()).x()/mm;
397  y = ((*hits)[i]->GetPos()).y()/mm;
398  z = ((*hits)[i]->GetPos()).z()/mm;
399  pmtfile << i << "\t"
400  << x << "\t"
401  << y << "\t"
402  << z << G4endl;
403 
404  man->FillH2(1,x/mm, y/mm); // fill(x,y)
405  if (event_id == 0 ) {
406  man->FillH2(2,x,y);
407  }
408 
409  //Fill ntuple #3
410  man->FillNtupleDColumn(3,0,event_id);
411  man->FillNtupleDColumn(3,1,(G4double) i);
412  man->FillNtupleDColumn(3,2,x);
413  man->FillNtupleDColumn(3,3,y);
414  man->FillNtupleDColumn(3,4,z);
415  man->AddNtupleRow(3);
416 
417  }
418  if (event_id%printModulo == 0 && P_hits > 0)
419  G4cout << " " << P_hits << " PMT hits in " << filename << G4endl;
420  pmtfile.close();
421  }
422 
423 }
424 
425 
426 
427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
429 
431  G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers");
432  G4TrajectoryContainer* trajContainer = evt->GetTrajectoryContainer();
433  G4int n_trajectories = 0;
434 
435  if(trajContainer) n_trajectories = trajContainer->entries();
436  for (G4int i=0; i<n_trajectories; i++) {
437  G4Trajectory* trj = (G4Trajectory*)(*trajContainer)[i];
438  if (drawTrksFlag == "all")
439  trj->DrawTrajectory();
440  else if ((drawTrksFlag == "charged") && (trj->GetCharge() != 0.))
441  trj->DrawTrajectory();
442  else if ((drawTrksFlag == "noscint")
443  && (trj->GetParticleName() != "opticalphoton"))
444  trj->DrawTrajectory();
445  }
446 
447  // G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
448  }
449 
450 }
static const double MeV
Definition: G4SIunits.hh:193
G4double totEnergyNeutrons
virtual void BeginOfEventAction(const G4Event *)
G4VHitsCollection * GetHC(G4int i)
G4int GetCollectionID(G4String colName)
Definition: G4SDManager.cc:131
G4String drawTrksFlag
static G4VVisManager * GetConcreteInstance()
G4double z
Definition: TRTMaterials.hh:39
G4double GetCharge() const
G4int entries() const
virtual ~DMXEventAction()
G4double totEnergyGammas
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4double totEnergy
int G4int
Definition: G4Types.hh:78
G4TrajectoryContainer * GetTrajectoryContainer() const
Definition: G4Event.hh:178
static G4UImanager * GetUIpointer()
Definition: G4UImanager.cc:58
G4int GetEventID() const
Definition: G4Event.hh:140
void writeScintHitsToFile(void)
G4int scintillatorCollID
G4GLOB_DLL std::ostream G4cout
G4bool FillNtupleDColumn(G4int id, G4double value)
G4double firstParticleE
virtual void EndOfEventAction(const G4Event *)
G4String firstParticleName
G4String drawColsFlag
G4String GetsavehitsFile()
Definition: DMXRunAction.hh:70
void drawTracks(const G4Event *)
static const double nanosecond
Definition: G4SIunits.hh:139
G4double firstLXeHitTime
virtual void DrawTrajectory() const
G4bool FillH2(G4int id, G4double xvalue, G4double yvalue, G4double weight=1.0)
G4bool FillH1(G4int id, G4double value, G4double weight=1.0)
static G4SDManager * GetSDMpointer()
Definition: G4SDManager.cc:40
virtual void DrawAllHits()
DMXEventAction(DMXRunAction *, DMXPrimaryGeneratorAction *)
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:195
G4HCofThisEvent * GetHCofThisEvent() const
Definition: G4Event.hh:174
const long * seeds
double G4double
Definition: G4Types.hh:76
DMXRunAction * runAct
G4String GetsavepmtFile()
Definition: DMXRunAction.hh:71
DMXPrimaryGeneratorAction * genAction
G4String particleName
static const double mm
Definition: G4SIunits.hh:102
G4double energy_pri
#define ns
Definition: xmlparse.cc:597
G4String GetParticleName() const
Definition: G4Trajectory.hh:99
void writePmtHitsToFile(const DMXPmtHitsCollection *)
G4double aveTimePmtHits
G4int ApplyCommand(const char *aCommand)
Definition: G4UImanager.cc:419
DMXEventActionMessenger * eventMessenger
G4double particleEnergy