Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/underground_physics/src/DMXEventAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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