170 lines
5.0 KiB
C++
170 lines
5.0 KiB
C++
|
|
#include "CloverEventAction.hh"
|
|
#include "CloverCrystalSD.hh"
|
|
#include "CloverCrystalHit.hh"
|
|
// #include "CloverAnalysis.hh"
|
|
#include "Detector_Parameters.hh"
|
|
|
|
|
|
#include "G4RunManager.hh"
|
|
#include "G4Event.hh"
|
|
#include "G4SDManager.hh"
|
|
#include "G4HCofThisEvent.hh"
|
|
#include "G4UnitsTable.hh"
|
|
#include "G4SystemOfUnits.hh"
|
|
|
|
//#include <cmath.h>
|
|
#include <stdlib.h>
|
|
|
|
#include "Randomize.hh"
|
|
#include <iomanip>
|
|
#include <vector>
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
CloverEventAction::CloverEventAction()
|
|
: G4UserEventAction(),
|
|
fCrystalHCID(-1),
|
|
fNDet(Crystal_Num) // number of crystal
|
|
{
|
|
fdEList.clear();
|
|
fdLList.clear();
|
|
for(int i = 0 ; i < fNDet; i++){
|
|
fdEList.push_back(0.);
|
|
fdLList.push_back(0.);
|
|
}
|
|
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
CloverEventAction::~CloverEventAction()
|
|
{}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
CloverCrystalHitsCollection* CloverEventAction::GetHitsCollection(G4int hcID, const G4Event* event) const
|
|
{
|
|
auto hitsCollection = static_cast<CloverCrystalHitsCollection*>(event->GetHCofThisEvent()->GetHC(hcID));
|
|
|
|
if ( ! hitsCollection ) {
|
|
G4ExceptionDescription msg;
|
|
msg << "Cannot access hitsCollection ID " << hcID;
|
|
G4Exception("CloverEventAction::GetHitsCollection()",
|
|
"MyCode0003", FatalException, msg);
|
|
}
|
|
|
|
return hitsCollection;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void CloverEventAction::PrintEventStatistics( G4double absoEdep, G4double absoTrackLength) const
|
|
{
|
|
// print event statistics
|
|
G4cout
|
|
<< "CloverEventAction::PrintEventStatistics" << G4endl
|
|
<< " Absorber: total energy: "
|
|
<< std::setw(7) << G4BestUnit(absoEdep, "Energy")
|
|
<< " total track length: "
|
|
<< std::setw(7) << G4BestUnit(absoTrackLength, "Length")
|
|
<< G4endl;
|
|
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void CloverEventAction::BeginOfEventAction(const G4Event* /*event*/)
|
|
{}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void CloverEventAction::EndOfEventAction(const G4Event* event)
|
|
{
|
|
|
|
G4TrajectoryContainer* trajectoryContainer = event->GetTrajectoryContainer();
|
|
G4int n_trajectories = 0;
|
|
if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();
|
|
|
|
// periodic printing
|
|
//G4cout << "######################################### CloverEventAction::EndOfEventAction" << G4endl;
|
|
//G4int eventID = event->GetEventID();
|
|
//if ( eventID < 100 || eventID % 1000 == 0) {
|
|
// G4cout << ">>>>>>>> Event: " << eventID << G4endl;
|
|
// if ( trajectoryContainer ) {
|
|
// G4cout << " " << n_trajectories
|
|
// << " trajectories stored in this event." << G4endl;
|
|
// }
|
|
// G4VHitsCollection* hc = event->GetHCofThisEvent()->GetHC(0);
|
|
// G4cout << " "
|
|
// << hc->GetSize() << " hits stored in this event" << G4endl;
|
|
//}
|
|
|
|
G4PrimaryVertex * pv = event->GetPrimaryVertex();
|
|
G4PrimaryParticle * pp = pv->GetPrimary();
|
|
|
|
G4double beamEnergy = pp->GetKineticEnergy() * 1000; // to keV
|
|
|
|
G4ThreeVector pmomt = pp->GetMomentumDirection();
|
|
G4double beamTheta = pmomt.theta() / degree;
|
|
G4double beamPhi = pmomt.phi() / degree;
|
|
|
|
//G4cout << "Beam Energy : " << beamEnergy << G4endl;
|
|
// << "Beam theta : " << beamTheta << G4endl
|
|
// << "Beam phi : " << beamPhi << G4endl;
|
|
|
|
// Get hits collections IDs (only once)
|
|
if ( fCrystalHCID == -1 ) {
|
|
fCrystalHCID = G4SDManager::GetSDMpointer()->GetCollectionID("CrystalHitsCollection");
|
|
}
|
|
|
|
// Get hits collections
|
|
auto crystalHC = GetHitsCollection(fCrystalHCID, event); //this is G4VHitsCollection
|
|
|
|
for( int i = 0 ; i < fNDet ; i++){
|
|
CloverCrystalHit * crystalHit = (*crystalHC)[i]; //this is CloverCrystalHit :: G4VHit
|
|
|
|
//add detector resolution of 1%
|
|
|
|
G4double edep = crystalHit->GetEdep();
|
|
|
|
|
|
G4double resol = G4RandGauss::shoot(1, 0.001);
|
|
|
|
//Energy Resolution of the detectors is Energy dependent.
|
|
/*
|
|
double Fano = 0.112;
|
|
double electron_hole_work_function_GE = .003; //3 ev in HPGe.
|
|
//
|
|
|
|
double energy_width = 2*std::pow((2* 0.693 * Fano * edep * electron_hole_work_function_GE),0.5);
|
|
|
|
double electron_width = 0.001;
|
|
|
|
double e_res = std::pow( std::pow(energy_width,2) + std::pow(electron_width,2),0.5);
|
|
|
|
|
|
G4double resol = G4RandGauss::shoot(1,e_res);
|
|
*/
|
|
fdEList[i] = edep * resol * 1000; // to keV
|
|
fdLList[i] = crystalHit->GetStepLength();
|
|
|
|
}
|
|
|
|
|
|
// get analysis manager
|
|
auto analysisManager = G4AnalysisManager::Instance();
|
|
|
|
// fill ntuple
|
|
analysisManager->FillNtupleIColumn(0, n_trajectories);
|
|
|
|
analysisManager->FillNtupleDColumn(3, beamEnergy);
|
|
analysisManager->FillNtupleDColumn(4, beamTheta);
|
|
analysisManager->FillNtupleDColumn(5, beamPhi);
|
|
|
|
analysisManager->AddNtupleRow();
|
|
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|