CLARION2_GEANT4/CloverEventAction.cc

#include "CloverEventAction.hh"
#include "CloverCrystalSD.hh"
#include "CloverCrystalHit.hh"
#include "CloverAnalysis.hh"
#include "Detector_Parameters.hh"

#include "BGOSD.hh"
#include "BGOHIT.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),
   fBGOHCID(-1),
   fNBGODet(Crystal_Num/4), // number of bgo
   fNDet(Crystal_Num) // number of crystal   
{
  fdEList.clear();
  fdLList.clear();
  fBEList.clear();
  fBLList.clear();
  for(int i = 0 ; i < fNDet; i++){
    fdEList.push_back(0.);
    fdLList.push_back(0.);
    fBEList.push_back(0.);
    fBLList.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;
}    

BGOHitsCollection* CloverEventAction::GetHitsBGOCollection(G4int hcID, const G4Event* event) const
{
  auto hitsCollection  = static_cast<BGOHitsCollection*>(event->GetHCofThisEvent()->GetHC(hcID));
  
  if ( ! hitsCollection ) {
    G4ExceptionDescription msg;
    msg << "Cannot access hitsCollection ID " << hcID; 
    G4Exception("CloverEventAction::GetHitsCollection() BGOHIT",
      "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");
  }

  if ( fBGOHCID == -1 ) {
    fBGOHCID = G4SDManager::GetSDMpointer()->GetCollectionID("BGOHitsCollection");
  }
  // 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();
  //  std::cout<< "crystalhere "<< i << ", " << fdEList[i] << std::endl;
  }

  auto BGOHC = GetHitsBGOCollection(fBGOHCID, event); //this is G4VHitsCollection

  for( int i = 0 ; i < fNBGODet ; i++){
    BGOHit * bgoHit = (*BGOHC)[i];  //this is CloverCrystalHit :: G4VHit

    //add detector resolution of 1%

    G4double edep = bgoHit->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);
  */  
    fBEList[i] = edep * resol * 1000;  // to keV
    fBLList[i] = bgoHit->GetStepLength();

   // std::cout<< "bgohere "<< i << ", " << fBEList[i] << std::endl;
  }
  
  // get analysis manager
  auto analysisManager = G4AnalysisManager::Instance();
 
  // fill ntuple
  analysisManager->FillNtupleIColumn(0, n_trajectories);

  analysisManager->FillNtupleDColumn(5, beamEnergy);
  analysisManager->FillNtupleDColumn(6, beamTheta);
  analysisManager->FillNtupleDColumn(7, beamPhi);
  
  analysisManager->AddNtupleRow();

}  

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......