ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DMXEventAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DMXEventAction.cc
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 
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 
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 
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 }