Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /geant4/examples/advanced/underground_physics/src/DMXEventAction.cc was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // --------------------------------------------------------------
0028 //   GEANT 4 - Underground Dark Matter Detector Advanced Example
0029 //
0030 //      For information related to this code contact: Alex Howard
0031 //      e-mail: alexander.howard@cern.ch
0032 // --------------------------------------------------------------
0033 // Comments
0034 //
0035 //                  Underground Advanced
0036 //               by A. Howard and H. Araujo 
0037 //                    (27th November 2001)
0038 //
0039 // History/Additions:
0040 // 16 Jan 2002  Added analysis
0041 //
0042 //
0043 // EventAction program
0044 // --------------------------------------------------------------
0045 
0046 #include "DMXEventAction.hh"
0047 
0048 // pass parameters for messengers:
0049 #include "DMXRunAction.hh"
0050 #include "DMXPrimaryGeneratorAction.hh"
0051 
0052 // note DMXPmtHit.hh and DMXScintHit.hh are included in DMXEventAction.hh
0053 
0054 #include "DMXEventActionMessenger.hh"
0055 
0056 #include "G4AnalysisManager.hh"
0057 #include "G4Event.hh"
0058 #include "G4EventManager.hh"
0059 #include "G4HCofThisEvent.hh"
0060 #include "G4VHitsCollection.hh"
0061 #include "G4TrajectoryContainer.hh"
0062 #include "G4Trajectory.hh"
0063 #include "G4VVisManager.hh"
0064 #include "G4SDManager.hh"
0065 #include "G4UImanager.hh"
0066 #include "G4UnitsTable.hh"
0067 #include "G4RunManager.hh"
0068 #include "G4Threading.hh"
0069 #include <fstream>
0070 #include <iomanip>
0071 
0072 #include "G4SystemOfUnits.hh"
0073 
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0075 DMXEventAction::DMXEventAction() 
0076   : runAct(0),genAction(0),hitsfile(0),pmtfile(0)
0077 {
0078 
0079   // create messenger
0080   eventMessenger = new DMXEventActionMessenger(this);
0081 
0082   // defaults for messenger
0083   drawColsFlag = "standard";
0084   drawTrksFlag = "all";
0085   drawHitsFlag = 1;
0086   savePmtFlag  = 0;
0087   saveHitsFlag = 1;
0088   
0089   printModulo = 1;
0090 
0091   // hits collections
0092   scintillatorCollID = -1;
0093   pmtCollID = -1;
0094 
0095   energy_pri=0;
0096   seeds=NULL;
0097 
0098 }
0099 
0100 
0101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0102 DMXEventAction::~DMXEventAction() {
0103 
0104   if (hitsfile)
0105     {
0106       hitsfile->close();
0107       delete hitsfile;
0108     }
0109   if (pmtfile)
0110     {
0111       pmtfile->close();
0112       delete pmtfile;
0113     }
0114   delete eventMessenger;
0115 }
0116 
0117 
0118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0119 void DMXEventAction::BeginOfEventAction(const G4Event* evt) 
0120 {
0121 
0122   //thread-local run action
0123   if (!runAct) 
0124     runAct = 
0125       dynamic_cast<const DMXRunAction*>
0126       (G4RunManager::GetRunManager()->GetUserRunAction());
0127   
0128   if (!genAction)
0129     genAction = dynamic_cast<const DMXPrimaryGeneratorAction*>
0130       (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
0131 
0132 
0133   // grab seeds
0134   seeds = genAction->GetEventSeeds();
0135 
0136   // grab energy of primary
0137   energy_pri = genAction->GetEnergyPrimary();
0138 
0139   event_id = evt->GetEventID();
0140  
0141   // print this information event by event (modulo n)   
0142   if (event_id%printModulo == 0)
0143     {
0144       G4cout << "\n---> Begin of event: " << event_id << G4endl;
0145       G4cout << "\n     Primary Energy: " << G4BestUnit(energy_pri,"Energy") 
0146          << G4endl;
0147       //      HepRandom::showEngineStatus(); 
0148     }
0149 
0150 
0151   // get ID for scintillator hits collection
0152   if (scintillatorCollID==-1) {
0153     G4SDManager *SDman = G4SDManager::GetSDMpointer();
0154     scintillatorCollID = SDman->GetCollectionID("scintillatorCollection");
0155   }
0156 
0157   // get ID for pmt hits collection
0158   if (pmtCollID==-1) {
0159     G4SDManager *SDman = G4SDManager::GetSDMpointer();
0160     pmtCollID = SDman->GetCollectionID("pmtCollection");
0161   }
0162 
0163 }
0164 
0165 
0166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0167 void DMXEventAction::EndOfEventAction(const G4Event* evt) {
0168 
0169   // check that both hits collections have been defined
0170   if(scintillatorCollID<0||pmtCollID<0) return;
0171 
0172   G4AnalysisManager* man = G4AnalysisManager::Instance();
0173 
0174   // address hits collections
0175   DMXScintHitsCollection* SHC = NULL;
0176   DMXPmtHitsCollection*   PHC = NULL;
0177   G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
0178   if(HCE) {
0179     SHC = (DMXScintHitsCollection*)(HCE->GetHC(scintillatorCollID));
0180     PHC = (DMXPmtHitsCollection*)(HCE->GetHC(pmtCollID));
0181   }
0182 
0183   // event summary
0184   totEnergy         = 0.;
0185   totEnergyGammas   = 0.;
0186   totEnergyNeutrons = 0.;
0187   firstParticleE    = 0.;
0188   particleEnergy    = 0.;
0189   firstLXeHitTime   = 0.;
0190   aveTimePmtHits    = 0.;
0191 
0192   firstParticleName = "";
0193   particleName      = "";
0194 
0195 
0196   // particle flags
0197   gamma_ev          = false;
0198   neutron_ev        = false;
0199   positron_ev       = false;
0200   electron_ev       = false;
0201   proton_ev         = false;
0202   other_ev          = false;
0203   start_gamma       = false;
0204   start_neutron     = false;
0205 
0206 
0207   // scintillator hits
0208   if(SHC) {
0209     S_hits = SHC->entries();
0210     
0211     for (G4int i=0; i<S_hits; i++) {
0212       if(i==0) {
0213     firstParticleName = (*SHC)[0]->GetParticle();
0214     firstLXeHitTime   = (*SHC)[0]->GetTime();
0215     firstParticleE = (*SHC)[0]->GetParticleEnergy();
0216     if (event_id%printModulo == 0 && S_hits > 0) {
0217       G4cout << "     First hit in LXe: " << firstParticleName << G4endl;
0218       G4cout << "     Number of hits in LXe: " << S_hits << G4endl; 
0219     }
0220       }
0221       hitEnergy         = (*SHC)[i]->GetEdep();
0222       totEnergy        += hitEnergy;
0223       
0224       particleName      = (*SHC)[i]->GetParticle();
0225       particleEnergy    = (*SHC)[i]->GetParticleEnergy();
0226 
0227       if(particleName == "gamma") {
0228     gamma_ev = true;
0229     start_gamma = true;
0230     start_neutron = false;
0231       }
0232       else if(particleName == "neutron") 
0233     neutron_ev = true;
0234       else if(particleName == "e+") 
0235     positron_ev = true;
0236       else if(particleName == "e-") 
0237     electron_ev = true;
0238       else if(particleName == "proton") 
0239     proton_ev = true;
0240       else {
0241     other_ev = true;
0242     start_gamma = false;
0243     start_neutron = true;
0244       }
0245 
0246       if(start_gamma && !start_neutron) 
0247     totEnergyGammas += hitEnergy;
0248       if(start_neutron && !start_gamma) 
0249     totEnergyNeutrons += hitEnergy;
0250     }
0251     
0252     if (event_id%printModulo == 0)
0253       G4cout << "     Total energy in LXe: " 
0254          << G4BestUnit(totEnergy,"Energy") << G4endl;  
0255     
0256   }
0257   
0258   
0259   // PMT hits
0260   if(PHC) {
0261     P_hits = PHC->entries();
0262     
0263     // average time of PMT hits
0264     for (G4int i=0; i<P_hits; i++) {
0265       G4double time = ( (*PHC)[i]->GetTime() - firstLXeHitTime );
0266       aveTimePmtHits += time / (G4double)P_hits;
0267       ////      if (event_id == 0) {
0268       if(P_hits >= 0) {
0269     man->FillH1(7,time);
0270       }
0271     }
0272   
0273     if (event_id%printModulo == 0 && P_hits > 0) {
0274       G4cout << "     Average light collection time: "
0275          << G4BestUnit(aveTimePmtHits,"Time") << G4endl;
0276       G4cout << "     Number of PMT hits (photoelectron equivalent): " 
0277          << P_hits << G4endl;     
0278     }
0279     // write out (x,y,z) of PMT hits
0280     if (savePmtFlag)
0281       writePmtHitsToFile(PHC);
0282   }
0283   
0284 
0285   // write out event summary
0286   if(saveHitsFlag) 
0287     writeScintHitsToFile();
0288   
0289   // draw trajectories
0290   if(drawColsFlag=="standard" && drawTrksFlag!="none")
0291     drawTracks(evt);
0292 
0293   // hits in PMT
0294   if(drawHitsFlag && PHC) 
0295     PHC->DrawAllHits();
0296 
0297   // print this event by event (modulo n)   
0298   if (event_id%printModulo == 0) 
0299     G4cout << "\n---> End of event: " << event_id << G4endl << G4endl;  
0300 
0301 }
0302 
0303 
0304 
0305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0306 void DMXEventAction::writeScintHitsToFile() 
0307 {
0308 
0309   G4String filename="hits.out";
0310   if (runAct)
0311     filename=runAct->GetsavehitsFile();
0312 
0313  
0314   
0315   //First time it is inkoved
0316   if (!hitsfile)
0317     {
0318       //check that we are in a worker: returns -1 in a master and -2 in sequential
0319       //one file per thread is produced ending with ".N", with N= thread number
0320       if (G4Threading::G4GetThreadId() >= 0)
0321     {
0322       std::stringstream sss;
0323       sss << filename.c_str() << "." << G4Threading::G4GetThreadId();    
0324       filename = sss.str();
0325       //G4cout << "Filename is: " << filename << G4endl;
0326     }
0327       
0328       hitsfile = new std::ofstream;
0329       hitsfile->open(filename);
0330       (*hitsfile) <<"Evt     Eprim   Etot    LXe     LXeTime PMT     PMTTime Seed1           Seed2           First   Flags" 
0331            << G4endl;
0332       (*hitsfile) <<"#       MeV     MeV     hits    ns      hits    ns                                      hit"
0333            << G4endl
0334            << G4endl;
0335     }
0336 
0337   if(S_hits) {
0338 
0339     if(hitsfile->is_open()) {
0340 
0341 
0342       (*hitsfile) << std::setiosflags(std::ios::fixed)
0343           << std::setprecision(4)
0344           << std::setiosflags(std::ios::left)
0345           << std::setw(6)
0346           << event_id << "\t"
0347           << energy_pri/MeV << "\t" 
0348           << totEnergy/MeV << "\t"
0349           << S_hits  << "\t"
0350           << std::setiosflags(std::ios::scientific) 
0351           << std::setprecision(2)
0352           << firstLXeHitTime/nanosecond << "\t"
0353           << P_hits << "\t"
0354           << std::setiosflags(std::ios::fixed) 
0355           << std::setprecision(4)
0356           << aveTimePmtHits/nanosecond << "\t"
0357           << *seeds     << "\t"
0358           << *(seeds+1) << "\t"
0359           << firstParticleName << "\t"
0360           << (gamma_ev    ? "gamma " : "") 
0361           << (neutron_ev  ? "neutron " : "") 
0362           << (positron_ev ? "positron " : "") 
0363           << (electron_ev ? "electron " : "") 
0364           << (other_ev    ? "other " : "") 
0365           << G4endl;
0366 
0367       if (event_id%printModulo == 0)
0368     G4cout << "     Event summary in file " << filename << G4endl;  
0369     }
0370 
0371     G4AnalysisManager* man = G4AnalysisManager::Instance();
0372     G4int firstparticleIndex = 0;
0373     if(firstParticleName == "gamma") firstparticleIndex = 1;
0374     else if (firstParticleName == "neutron") firstparticleIndex = 2;
0375     else if(firstParticleName == "electron") firstparticleIndex = 3;
0376     else if(firstParticleName == "positron") firstparticleIndex = 4;
0377     else{
0378       firstparticleIndex = 5;
0379       man->FillH1(3,totEnergy);
0380     }
0381 
0382     man->FillH1(4,P_hits,10); //weight
0383     man->FillH1(5,P_hits);
0384 
0385     man->FillH1(1,energy_pri/keV);
0386     man->FillH1(2,totEnergy/keV);
0387     man->FillH1(6,aveTimePmtHits/ns);
0388 
0389     long seed1 = *seeds;
0390     long seed2 = *(seeds+1);    
0391     
0392     //Fill ntuple #2
0393     man->FillNtupleDColumn(2,0,event_id);
0394     man->FillNtupleDColumn(2,1,energy_pri/keV);
0395     man->FillNtupleDColumn(2,2,totEnergy);
0396     man->FillNtupleDColumn(2,3,S_hits);
0397     man->FillNtupleDColumn(2,4,firstLXeHitTime);
0398     man->FillNtupleDColumn(2,5,P_hits);
0399     man->FillNtupleDColumn(2,6,aveTimePmtHits);
0400     man->FillNtupleDColumn(2,7,firstparticleIndex);
0401     man->FillNtupleDColumn(2,8,firstParticleE);
0402     man->FillNtupleDColumn(2,9,gamma_ev);
0403     man->FillNtupleDColumn(2,10,neutron_ev);
0404     man->FillNtupleDColumn(2,11,positron_ev);
0405     man->FillNtupleDColumn(2,12,electron_ev);
0406     man->FillNtupleDColumn(2,13,other_ev);
0407     man->FillNtupleDColumn(2,14,seed1);
0408     man->FillNtupleDColumn(2,15,seed2);
0409     man->AddNtupleRow(2);
0410     
0411   }
0412 
0413 }
0414 
0415 
0416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0417 void DMXEventAction::writePmtHitsToFile(const DMXPmtHitsCollection* hits) 
0418 { 
0419   G4String filename="pmt.out";
0420   if (runAct)
0421     filename=runAct->GetsavepmtFile();
0422   
0423 
0424   //first time it is invoked: create it
0425   if (!pmtfile)
0426     {
0427       //check that we are in a worker: returns -1 in a master and -2 in sequential
0428       //one file per thread is produced ending with ".N", with N= thread number
0429       if (G4Threading::G4GetThreadId() >= 0)
0430     {
0431       std::stringstream sss;
0432       sss << filename.c_str() << "." << G4Threading::G4GetThreadId();    
0433       filename = sss.str();
0434       //G4cout << "Filename is: " << filename << G4endl;
0435     }
0436       pmtfile = new std::ofstream;
0437       pmtfile->open(filename);
0438     }
0439 
0440 
0441   G4double x; G4double y; G4double z;
0442   G4AnalysisManager* man = G4AnalysisManager::Instance();
0443 
0444   if(pmtfile->is_open()) {
0445     (*pmtfile) << "Hit#    X, mm   Y, mm   Z, mm" << G4endl;       
0446     (*pmtfile) << std::setiosflags(std::ios::fixed)
0447            << std::setprecision(3)
0448            << std::setiosflags(std::ios::left)
0449            << std::setw(6);
0450     for (G4int i=0; i<P_hits; i++)
0451       {
0452     x = ((*hits)[i]->GetPos()).x()/mm;
0453     y = ((*hits)[i]->GetPos()).y()/mm;
0454     z = ((*hits)[i]->GetPos()).z()/mm;
0455     (*pmtfile) << i << "\t"
0456            << x << "\t" 
0457            << y << "\t"
0458            << z << G4endl;
0459     
0460     man->FillH2(1,x/mm, y/mm);  // fill(x,y)
0461     if (event_id == 0 ) {
0462       man->FillH2(2,x,y);
0463     }
0464 
0465     //Fill ntuple #3
0466     man->FillNtupleDColumn(3,0,event_id);
0467     man->FillNtupleDColumn(3,1,(G4double) i);
0468     man->FillNtupleDColumn(3,2,x);
0469     man->FillNtupleDColumn(3,3,y);
0470     man->FillNtupleDColumn(3,4,z);
0471     man->AddNtupleRow(3);
0472 
0473       }
0474     if (event_id%printModulo == 0 && P_hits > 0) 
0475       G4cout << "     " << P_hits << " PMT hits in " << filename << G4endl;  
0476   }
0477   
0478 }
0479 
0480 
0481 
0482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0483 void DMXEventAction::drawTracks(const G4Event* evt) {
0484 
0485   if(G4VVisManager::GetConcreteInstance()) {
0486     G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers");    
0487     G4TrajectoryContainer* trajContainer = evt->GetTrajectoryContainer();
0488     G4int n_trajectories = 0;
0489 
0490     if(trajContainer) n_trajectories = trajContainer->entries();
0491     for (G4int i=0; i<n_trajectories; i++) {
0492       G4Trajectory* trj = (G4Trajectory*)(*trajContainer)[i];
0493       if (drawTrksFlag == "all") 
0494     trj->DrawTrajectory();
0495       else if ((drawTrksFlag == "charged") && (trj->GetCharge() != 0.))
0496     trj->DrawTrajectory();
0497       else if ((drawTrksFlag == "noscint") 
0498            && (trj->GetParticleName() != "opticalphoton"))
0499     trj->DrawTrajectory();
0500     }
0501     
0502     // G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");    
0503   } 
0504 
0505 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!