AuxDetSD.cc

Go to the documentation of this file.

//
//               __        __ __  __  __
//   ____ ______/ /_____ _/ // / / /_/ /__
//  / __ `/ ___/ __/ __ `/ // /_/ __/ //_/
// / /_/ / /  / /_/ /_/ /__  __/ /_/ ,<
// \__,_/_/   \__/\__, /  /_/  \__/_/|_|
//               /____/
//
// artg4tk: art based Geant 4 Toolkit
//
//=============================================================================
// AuxDetSD.cc: Class representing a sensitive aux detector
// Author: Hans Wenzel (Fermilab)
//=============================================================================
#include "larg4/Services/AuxDetSD.h"
#include "Geant4/G4Event.hh"
#include "Geant4/G4EventManager.hh"
#include "Geant4/G4HCofThisEvent.hh"
#include "Geant4/G4SDManager.hh"
#include "Geant4/G4Step.hh"
#include "Geant4/G4SystemOfUnits.hh"
#include "Geant4/G4ThreeVector.hh"
#include "Geant4/G4UnitsTable.hh"
#include "Geant4/G4VSolid.hh"
#include "Geant4/G4VVisManager.hh"
#include "Geant4/G4ios.hh"
#include <algorithm>
#include <unordered_set>
//#define _verbose_ 1
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace larg4 {

  AuxDetSD::AuxDetSD(G4String name) : G4VSensitiveDetector(name)
  {
    hitCollection.clear();
  }

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

  AuxDetSD::~AuxDetSD() {}
  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  void AuxDetSD::Initialize(G4HCofThisEvent*)
  {
    hitCollection.clear();
    temphitCollection.clear();
  }
  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  G4bool AuxDetSD::ProcessHits(G4Step* step, G4TouchableHistory*)
  {
    G4double edep = step->GetTotalEnergyDeposit() / CLHEP::MeV;
    if (edep == 0.) return false;
    G4Track* track = step->GetTrack();
    const unsigned int trackID = track->GetTrackID();
    unsigned int ID = step->GetPreStepPoint()->GetPhysicalVolume()->GetCopyNo();
    TempHit tmpHit = TempHit(ID,
                             trackID,
                             track->GetParentID(),
                             step->IsFirstStepInVolume(),
                             step->IsLastStepInVolume(),
                             edep,
                             step->GetPreStepPoint()->GetPosition().getX() / CLHEP::cm,
                             step->GetPreStepPoint()->GetPosition().getY() / CLHEP::cm,
                             step->GetPreStepPoint()->GetPosition().getZ() / CLHEP::cm,
                             step->GetPreStepPoint()->GetGlobalTime() / CLHEP::ns,
                             step->GetPostStepPoint()->GetPosition().getX() / CLHEP::cm,
                             step->GetPostStepPoint()->GetPosition().getY() / CLHEP::cm,
                             step->GetPostStepPoint()->GetPosition().getZ() / CLHEP::cm,
                             step->GetPostStepPoint()->GetGlobalTime() / CLHEP::ns,
                             step->GetPostStepPoint()->GetMomentum().getX() / CLHEP::GeV,
                             step->GetPostStepPoint()->GetMomentum().getY() / CLHEP::GeV,
                             step->GetPostStepPoint()->GetMomentum().getZ() / CLHEP::GeV);
    temphitCollection.push_back(tmpHit);

    /*
  sim::AuxDetHit newHit = sim::AuxDetHit(ID,
                                    trackID,
                                    edep,
                                    step->GetPreStepPoint()->GetPosition().getX()/CLHEP::cm,
                                    step->GetPreStepPoint()->GetPosition().getY()/CLHEP::cm,
                                    step->GetPreStepPoint()->GetPosition().getZ()/CLHEP::cm,
                                    step->GetPreStepPoint()->GetGlobalTime()/CLHEP::ns,
                                    step->GetPostStepPoint()->GetPosition().getX()/CLHEP::cm,
                                    step->GetPostStepPoint()->GetPosition().getY()/CLHEP::cm,
                                    step->GetPostStepPoint()->GetPosition().getZ()/CLHEP::cm,
                                    step->GetPostStepPoint()->GetGlobalTime()/CLHEP::ns,
                                    step->GetPostStepPoint()->GetMomentum().getX()/CLHEP::GeV,
                                    step->GetPostStepPoint()->GetMomentum().getY()/CLHEP::GeV,
                                    step->GetPostStepPoint()->GetMomentum().getZ()/CLHEP::GeV
                                    );
  hitCollection.push_back(newHit);
  */
    return true;
  }

  void AuxDetSD::EndOfEvent(G4HCofThisEvent*)
  {
    if (temphitCollection.size() == 0) return; // No hits so nothing to do
#if defined _verbose_
    std::cout << " EndOfEvent number of temp hits: " << temphitCollection.size() << std::endl;
    std::cout << " EndOfEvent number of aux hits: " << hitCollection.size() << std::endl;
#endif
    std::sort(temphitCollection.begin(), temphitCollection.end());
    int geoId = -1;
    int trackId = -1;
    unsigned int counter = 0;
    std::unordered_set<unsigned int> setofIDs;

    for (auto it = temphitCollection.begin(); it != temphitCollection.end(); it++) {
#if defined _verbose_
      std::cout << "geoID: " << it->GetID() << " track ID: " << it->GetTrackID()
                << " Edep: " << it->GetEnergyDeposited() << " Parent Id: " << it->GetParentID()
                << " exit Time:  " << it->GetExitT() << " is first: " << it->IsIsfirstinVolume()
                << " is last:  " << it->IsIslastinVolume();
#endif
      if (it->GetID() == geoId && trackId == it->GetTrackID()) // trackid and detector didn't change
      {
#if defined _verbose_
        std::cout << " A" << std::endl;
#endif
        if (it->GetExitT()) // change exit vector and add to total charge
        {
          hitCollection[counter - 1].SetEnergyDeposited(
            hitCollection[counter - 1].GetEnergyDeposited() + it->GetEnergyDeposited());
          hitCollection[counter - 1].SetExitX(it->GetExitX());
          hitCollection[counter - 1].SetExitY(it->GetExitY());
          hitCollection[counter - 1].SetExitZ(it->GetExitZ());
          hitCollection[counter - 1].SetExitT(it->GetExitT());
        }
        else // just add to total charge
        {
          hitCollection[counter - 1].SetEnergyDeposited(
            hitCollection[counter - 1].GetEnergyDeposited() + it->GetEnergyDeposited());
        }
      }
      else if (setofIDs.find(it->GetParentID()) != setofIDs.end()) {
        setofIDs.insert(it->GetTrackID());
#if defined _verbose_
        std::cout << " A" << std::endl;
#endif
        hitCollection[counter - 1].SetEnergyDeposited(
          hitCollection[counter - 1].GetEnergyDeposited() + it->GetEnergyDeposited());
      }
      else if (it->GetID() != geoId) // new Detector
      {
        geoId = it->GetID();
        trackId = it->GetTrackID();
        setofIDs.clear();
        setofIDs.insert(it->GetTrackID());
#if defined _verbose_
        std::cout << " N" << std::endl;
#endif
        counter++;
        hitCollection.push_back(sim::AuxDetHit(it->GetID(),
                                               it->GetTrackID(),
                                               it->GetEnergyDeposited(),
                                               it->GetEntryX(),
                                               it->GetEntryY(),
                                               it->GetEntryZ(),
                                               it->GetEntryT(),
                                               it->GetExitX(),
                                               it->GetExitY(),
                                               it->GetExitZ(),
                                               it->GetExitT(),
                                               it->GetExitMomentumX(),
                                               it->GetExitMomentumY(),
                                               it->GetExitMomentumZ()));
      }
      else {
        trackId = it->GetTrackID();
#if defined _verbose_
        std::cout << " N" << std::endl;
#endif
        counter++;
        setofIDs.clear();
        setofIDs.insert(it->GetTrackID());
        hitCollection.push_back(sim::AuxDetHit(it->GetID(),
                                               it->GetTrackID(),
                                               it->GetEnergyDeposited(),
                                               it->GetEntryX(),
                                               it->GetEntryY(),
                                               it->GetEntryZ(),
                                               it->GetEntryT(),
                                               it->GetExitX(),
                                               it->GetExitY(),
                                               it->GetExitZ(),
                                               it->GetExitT(),
                                               it->GetExitMomentumX(),
                                               it->GetExitMomentumY(),
                                               it->GetExitMomentumZ()));
      }
    }
#if defined _verbose_
    std::cout << "Number of AuxDetHits: " << counter << std::endl;
#endif
  } // EndOfEvent
} // namespace sim