Geant4  10.01.p02
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 
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 "G4RunManager.hh"
69 #include "G4Threading.hh"
70 #include <fstream>
71 #include <iomanip>
72 
73 #include "G4SystemOfUnits.hh"
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77  : runAct(0),genAction(0),hitsfile(0),pmtfile(0)
78 {
79 
80  // create messenger
82 
83  // defaults for messenger
84  drawColsFlag = "standard";
85  drawTrksFlag = "all";
86  drawHitsFlag = 1;
87  savePmtFlag = 0;
88  saveHitsFlag = 1;
89 
90  printModulo = 1;
91 
92  // hits collections
93  scintillatorCollID = -1;
94  pmtCollID = -1;
95 
96  energy_pri=0;
97  seeds=NULL;
98 
99 }
100 
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 
105  if (hitsfile)
106  {
107  hitsfile->close();
108  delete hitsfile;
109  }
110  if (pmtfile)
111  {
112  pmtfile->close();
113  delete pmtfile;
114  }
115  delete eventMessenger;
116 }
117 
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121 {
122 
123  //thread-local run action
124  if (!runAct)
125  runAct =
126  dynamic_cast<const DMXRunAction*>
128 
129  if (!genAction)
130  genAction = dynamic_cast<const DMXPrimaryGeneratorAction*>
132 
133 
134  // grab seeds
136 
137  // grab energy of primary
139 
140  event_id = evt->GetEventID();
141 
142  // print this information event by event (modulo n)
143  if (event_id%printModulo == 0)
144  {
145  G4cout << "\n---> Begin of event: " << event_id << G4endl;
146  G4cout << "\n Primary Energy: " << G4BestUnit(energy_pri,"Energy")
147  << G4endl;
148  // HepRandom::showEngineStatus();
149  }
150 
151 
152  // get ID for scintillator hits collection
153  if (scintillatorCollID==-1) {
155  scintillatorCollID = SDman->GetCollectionID("scintillatorCollection");
156  }
157 
158  // get ID for pmt hits collection
159  if (pmtCollID==-1) {
161  pmtCollID = SDman->GetCollectionID("pmtCollection");
162  }
163 
164 }
165 
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169 
170  // check that both hits collections have been defined
171  if(scintillatorCollID<0||pmtCollID<0) return;
172 
173  G4AnalysisManager* man = G4AnalysisManager::Instance();
174 
175  // address hits collections
176  DMXScintHitsCollection* SHC = NULL;
177  DMXPmtHitsCollection* PHC = NULL;
178  G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
179  if(HCE) {
181  PHC = (DMXPmtHitsCollection*)(HCE->GetHC(pmtCollID));
182  }
183 
184  // event summary
185  totEnergy = 0.;
186  totEnergyGammas = 0.;
187  totEnergyNeutrons = 0.;
188  firstParticleE = 0.;
189  particleEnergy = 0.;
190  firstLXeHitTime = 0.;
191  aveTimePmtHits = 0.;
192 
193  firstParticleName = "";
194  particleName = "";
195 
196 
197  // particle flags
198  gamma_ev = false;
199  neutron_ev = false;
200  positron_ev = false;
201  electron_ev = false;
202  proton_ev = false;
203  other_ev = false;
204  start_gamma = false;
205  start_neutron = false;
206 
207 
208  // scintillator hits
209  if(SHC) {
210  S_hits = SHC->entries();
211 
212  for (G4int i=0; i<S_hits; i++) {
213  if(i==0) {
214  firstParticleName = (*SHC)[0]->GetParticle();
215  firstLXeHitTime = (*SHC)[0]->GetTime();
216  firstParticleE = (*SHC)[0]->GetParticleEnergy();
217  if (event_id%printModulo == 0 && S_hits > 0) {
218  G4cout << " First hit in LXe: " << firstParticleName << G4endl;
219  G4cout << " Number of hits in LXe: " << S_hits << G4endl;
220  }
221  }
222  hitEnergy = (*SHC)[i]->GetEdep();
223  totEnergy += hitEnergy;
224 
225  particleName = (*SHC)[i]->GetParticle();
226  particleEnergy = (*SHC)[i]->GetParticleEnergy();
227 
228  if(particleName == "gamma") {
229  gamma_ev = true;
230  start_gamma = true;
231  start_neutron = false;
232  }
233  else if(particleName == "neutron")
234  neutron_ev = true;
235  else if(particleName == "e+")
236  positron_ev = true;
237  else if(particleName == "e-")
238  electron_ev = true;
239  else if(particleName == "proton")
240  proton_ev = true;
241  else {
242  other_ev = true;
243  start_gamma = false;
244  start_neutron = true;
245  }
246 
247  if(start_gamma && !start_neutron)
249  if(start_neutron && !start_gamma)
251  }
252 
253  if (event_id%printModulo == 0)
254  G4cout << " Total energy in LXe: "
255  << G4BestUnit(totEnergy,"Energy") << G4endl;
256 
257  }
258 
259 
260  // PMT hits
261  if(PHC) {
262  P_hits = PHC->entries();
263 
264  // average time of PMT hits
265  for (G4int i=0; i<P_hits; i++) {
266  G4double time = ( (*PHC)[i]->GetTime() - firstLXeHitTime );
267  aveTimePmtHits += time / (G4double)P_hits;
269  if(P_hits >= 0) {
270  man->FillH1(7,time);
271  }
272  }
273 
274  if (event_id%printModulo == 0 && P_hits > 0) {
275  G4cout << " Average light collection time: "
276  << G4BestUnit(aveTimePmtHits,"Time") << G4endl;
277  G4cout << " Number of PMT hits (photoelectron equivalent): "
278  << P_hits << G4endl;
279  }
280  // write out (x,y,z) of PMT hits
281  if (savePmtFlag)
282  writePmtHitsToFile(PHC);
283  }
284 
285 
286  // write out event summary
287  if(saveHitsFlag)
289 
290  // draw trajectories
291  if(drawColsFlag=="standard" && drawTrksFlag!="none")
292  drawTracks(evt);
293 
294  // hits in PMT
295  if(drawHitsFlag && PHC)
296  PHC->DrawAllHits();
297 
298  // print this event by event (modulo n)
299  if (event_id%printModulo == 0)
300  G4cout << "\n---> End of event: " << event_id << G4endl << G4endl;
301 
302 }
303 
304 
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
308 {
309 
310  G4String filename="hits.out";
311  if (runAct)
312  filename=runAct->GetsavehitsFile();
313 
314 
315 
316  //First time it is inkoved
317  if (!hitsfile)
318  {
319  //check that we are in a worker: returns -1 in a master and -2 in sequential
320  //one file per thread is produced ending with ".N", with N= thread number
321  if (G4Threading::G4GetThreadId() >= 0)
322  {
323  std::stringstream sss;
324  sss << filename.c_str() << "." << G4Threading::G4GetThreadId();
325  filename = sss.str();
326  //G4cout << "Filename is: " << filename << G4endl;
327  }
328 
329  hitsfile = new std::ofstream;
330  hitsfile->open(filename);
331  (*hitsfile) <<"Evt Eprim Etot LXe LXeTime PMT PMTTime Seed1 Seed2 First Flags"
332  << G4endl;
333  (*hitsfile) <<"# MeV MeV hits ns hits ns hit"
334  << G4endl
335  << G4endl;
336  }
337 
338  if(S_hits) {
339 
340  if(hitsfile->is_open()) {
341 
342 
343  (*hitsfile) << std::setiosflags(std::ios::fixed)
344  << std::setprecision(4)
345  << std::setiosflags(std::ios::left)
346  << std::setw(6)
347  << event_id << "\t"
348  << energy_pri/MeV << "\t"
349  << totEnergy/MeV << "\t"
350  << S_hits << "\t"
351  << std::setiosflags(std::ios::scientific)
352  << std::setprecision(2)
353  << firstLXeHitTime/nanosecond << "\t"
354  << P_hits << "\t"
355  << std::setiosflags(std::ios::fixed)
356  << std::setprecision(4)
357  << aveTimePmtHits/nanosecond << "\t"
358  << *seeds << "\t"
359  << *(seeds+1) << "\t"
360  << firstParticleName << "\t"
361  << (gamma_ev ? "gamma " : "")
362  << (neutron_ev ? "neutron " : "")
363  << (positron_ev ? "positron " : "")
364  << (electron_ev ? "electron " : "")
365  << (other_ev ? "other " : "")
366  << G4endl;
367 
368  if (event_id%printModulo == 0)
369  G4cout << " Event summary in file " << filename << G4endl;
370  }
371 
372  G4AnalysisManager* man = G4AnalysisManager::Instance();
373  G4int firstparticleIndex = 0;
374  if(firstParticleName == "gamma") firstparticleIndex = 1;
375  else if (firstParticleName == "neutron") firstparticleIndex = 2;
376  else if(firstParticleName == "electron") firstparticleIndex = 3;
377  else if(firstParticleName == "positron") firstparticleIndex = 4;
378  else{
379  firstparticleIndex = 5;
380  man->FillH1(3,totEnergy);
381  }
382 
383  man->FillH1(4,P_hits,10); //weight
384  man->FillH1(5,P_hits);
385 
386  man->FillH1(1,energy_pri/keV);
387  man->FillH1(2,totEnergy/keV);
388  man->FillH1(6,aveTimePmtHits/ns);
389 
390  long seed1 = *seeds;
391  long seed2 = *(seeds+1);
392 
393  //Fill ntuple #2
394  man->FillNtupleDColumn(2,0,event_id);
395  man->FillNtupleDColumn(2,1,energy_pri/keV);
396  man->FillNtupleDColumn(2,2,totEnergy);
397  man->FillNtupleDColumn(2,3,S_hits);
399  man->FillNtupleDColumn(2,5,P_hits);
401  man->FillNtupleDColumn(2,7,firstparticleIndex);
403  man->FillNtupleDColumn(2,9,gamma_ev);
404  man->FillNtupleDColumn(2,10,neutron_ev);
405  man->FillNtupleDColumn(2,11,positron_ev);
406  man->FillNtupleDColumn(2,12,electron_ev);
407  man->FillNtupleDColumn(2,13,other_ev);
408  man->FillNtupleDColumn(2,14,seed1);
409  man->FillNtupleDColumn(2,15,seed2);
410  man->AddNtupleRow(2);
411 
412  }
413 
414 }
415 
416 
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
419 {
420  G4String filename="pmt.out";
421  if (runAct)
422  filename=runAct->GetsavepmtFile();
423 
424 
425  //first time it is invoked: create it
426  if (!pmtfile)
427  {
428  //check that we are in a worker: returns -1 in a master and -2 in sequential
429  //one file per thread is produced ending with ".N", with N= thread number
430  if (G4Threading::G4GetThreadId() >= 0)
431  {
432  std::stringstream sss;
433  sss << filename.c_str() << "." << G4Threading::G4GetThreadId();
434  filename = sss.str();
435  //G4cout << "Filename is: " << filename << G4endl;
436  }
437  pmtfile = new std::ofstream;
438  pmtfile->open(filename);
439  }
440 
441 
442  G4double x; G4double y; G4double z;
443  G4AnalysisManager* man = G4AnalysisManager::Instance();
444 
445  if(pmtfile->is_open()) {
446  (*pmtfile) << "Hit# X, mm Y, mm Z, mm" << G4endl;
447  (*pmtfile) << std::setiosflags(std::ios::fixed)
448  << std::setprecision(3)
449  << std::setiosflags(std::ios::left)
450  << std::setw(6);
451  for (G4int i=0; i<P_hits; i++)
452  {
453  x = ((*hits)[i]->GetPos()).x()/mm;
454  y = ((*hits)[i]->GetPos()).y()/mm;
455  z = ((*hits)[i]->GetPos()).z()/mm;
456  (*pmtfile) << i << "\t"
457  << x << "\t"
458  << y << "\t"
459  << z << G4endl;
460 
461  man->FillH2(1,x/mm, y/mm); // fill(x,y)
462  if (event_id == 0 ) {
463  man->FillH2(2,x,y);
464  }
465 
466  //Fill ntuple #3
467  man->FillNtupleDColumn(3,0,event_id);
468  man->FillNtupleDColumn(3,1,(G4double) i);
469  man->FillNtupleDColumn(3,2,x);
470  man->FillNtupleDColumn(3,3,y);
471  man->FillNtupleDColumn(3,4,z);
472  man->AddNtupleRow(3);
473 
474  }
475  if (event_id%printModulo == 0 && P_hits > 0)
476  G4cout << " " << P_hits << " PMT hits in " << filename << G4endl;
477  }
478 
479 }
480 
481 
482 
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
485 
487  G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers");
488  G4TrajectoryContainer* trajContainer = evt->GetTrajectoryContainer();
489  G4int n_trajectories = 0;
490 
491  if(trajContainer) n_trajectories = trajContainer->entries();
492  for (G4int i=0; i<n_trajectories; i++) {
493  G4Trajectory* trj = (G4Trajectory*)(*trajContainer)[i];
494  if (drawTrksFlag == "all")
495  trj->DrawTrajectory();
496  else if ((drawTrksFlag == "charged") && (trj->GetCharge() != 0.))
497  trj->DrawTrajectory();
498  else if ((drawTrksFlag == "noscint")
499  && (trj->GetParticleName() != "opticalphoton"))
500  trj->DrawTrajectory();
501  }
502 
503  // G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");
504  }
505 
506 }
const G4VUserPrimaryGeneratorAction * GetUserPrimaryGeneratorAction() const
static const double MeV
Definition: G4SIunits.hh:193
G4double totEnergyNeutrons
G4String GetsavehitsFile() const
Definition: DMXRunAction.hh:69
void writeScintHitsToFile()
virtual void BeginOfEventAction(const G4Event *)
G4VHitsCollection * GetHC(G4int i)
G4int GetCollectionID(G4String colName)
Definition: G4SDManager.cc:131
G4String GetsavepmtFile() const
Definition: DMXRunAction.hh:70
G4String drawTrksFlag
static G4VVisManager * GetConcreteInstance()
G4double z
Definition: TRTMaterials.hh:39
G4double GetCharge() const
const DMXPrimaryGeneratorAction * genAction
G4int entries() const
virtual ~DMXEventAction()
G4double totEnergyGammas
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
const DMXRunAction * runAct
G4double totEnergy
int G4int
Definition: G4Types.hh:78
G4TrajectoryContainer * GetTrajectoryContainer() const
Definition: G4Event.hh:189
static G4UImanager * GetUIpointer()
Definition: G4UImanager.cc:58
G4int GetEventID() const
Definition: G4Event.hh:151
G4int scintillatorCollID
G4GLOB_DLL std::ostream G4cout
G4bool FillNtupleDColumn(G4int id, G4double value)
G4double firstParticleE
virtual void EndOfEventAction(const G4Event *)
G4String firstParticleName
G4String drawColsFlag
void drawTracks(const G4Event *)
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
static const double nanosecond
Definition: G4SIunits.hh:139
static const char sss[MAX_N_PAR+2]
Definition: Evaluator.cc:64
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
std::ofstream * pmtfile
virtual void DrawAllHits()
#define G4endl
Definition: G4ios.hh:61
std::ofstream * hitsfile
static const double keV
Definition: G4SIunits.hh:195
G4HCofThisEvent * GetHCofThisEvent() const
Definition: G4Event.hh:185
const G4UserRunAction * GetUserRunAction() const
const long * seeds
double G4double
Definition: G4Types.hh:76
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 G4GetThreadId()
Definition: G4Threading.cc:127
G4int ApplyCommand(const char *aCommand)
Definition: G4UImanager.cc:432
DMXEventActionMessenger * eventMessenger
G4double particleEnergy