In This Package:

DsPmtSensDet.cc

Go to the documentation of this file.

#include "DsPmtSensDet.h"
#include "DsPmtModel.h"

#include "Event/SimPmtHit.h"

#include "Conventions/Detectors.h"

#include "GiGaCnv/IGiGaGeomCnvSvc.h"
//#include "GiGaCnv/GiGaVolumeUtils.h"
#include "GaudiKernel/ToolHandle.h"

#include "DetDesc/DetectorElement.h"

#include "GaudiKernel/IRndmGenSvc.h"

#include "G4DataHelpers/ITouchableToDetectorElement.h"

#include "G4Step.hh"
#include "G4HCofThisEvent.hh"
#include "G4SDManager.hh"
#include "G4StepPoint.hh"
#include "G4Region.hh"
#include "G4LogicalVolume.hh"
#include "G4VFastSimulationModel.hh"
#include "G4ProductionCuts.hh"
#include "G4ThreeVector.hh"
#include "G4TouchableHistory.hh"

#include "G4OpticalPhoton.hh"

#include "CLHEP/Units/PhysicalConstants.h"

#include <string>
#include <vector>
#include <cstdio>
using namespace std;

DetectorId::DetectorId_t detector_ids[] = 
{  DetectorId::kAD1, DetectorId::kAD2, DetectorId::kAD3,
   DetectorId::kAD4, DetectorId::kIWS, DetectorId::kOWS,
   (DetectorId::DetectorId_t)-1 };

Site::Site_t site_ids[] = 
{ Site::kDayaBay, Site::kLingAo, Site::kFar, (Site::Site_t)-1 };


DsPmtSensDet::DsPmtSensDet(const std::string& type,
                           const std::string& name, 
                           const IInterface*  parent)
    : G4VSensitiveDetector(name)
    , GiGaSensDetBase(type,name,parent)
    , m_t2de(0)
{
    info() << "DsPmtSensDet (" << type << "/" << name << ") created" << endreq;

    declareProperty("CathodeLogicalVolume",
                    m_cathodeLogVols,
                    "Photo-Cathode logical volume to which this SD is attached.");

    declareProperty("TouchableToDetelem", m_t2deName = "TH2DE",
                    "The ITouchableToDetectorElement to use to resolve sensor.");

    declareProperty("SensorStructures",m_sensorStructures,
                    "TDS Paths in which to look for sensor detector elements"
                    " using this sensitive detector");

    declareProperty("PackedIdPropertyName",m_idParameter="PmtID",
                    "The name of the user property holding the PMT ID.");

    declareProperty("QEffParameterName",m_qeffParamName="EFFICIENCY",
                    "name of user parameter in the photo cathode volume that"
                    " holds the quantum efficiency tabproperty");

    declareProperty("QEScale",m_qeScale=1.0 / 0.9,
                    "Upward scaling of the quantum efficiency by inverse of mean PMT-to-PMT efficiency in electronics simulation.");
    
    declareProperty("PreQE", m_preQE=0.32,
                    "Miximal possible PMT QE value to be applied on the stage of "
                    "generating photons in scintillation/Cerenkov processes");

    declareProperty("ApplyPreQE",m_applyPreQE=true,
                    "Apply preQE on stage of generating photons of not.");

    m_cathodeLogVols.push_back("/dd/Geometry/PMT/lvPmtHemiCathode");
    m_cathodeLogVols.push_back("/dd/Geometry/PMT/lvHeadonPmtCathode");
}

DsPmtSensDet::~DsPmtSensDet()
{
}

StatusCode DsPmtSensDet::initialize()
{
    this->GiGaSensDetBase::initialize();
    info() << "DsPmtSensDet initialize" << endreq;

    m_t2de = tool<ITouchableToDetectorElement>(m_t2deName);

    // define collections.

    // first a catch-all
    collectionName.insert("unknown");
    for (int isite=0; site_ids[isite] >= 0; ++isite) {
        Site::Site_t site = site_ids[isite];

        for (int idet=0; detector_ids[idet] >= 0; ++idet) {
            DetectorId::DetectorId_t detid = detector_ids[idet];

            DayaBay::Detector det(site,detid);

            if (det.bogus()) continue;

            string name=det.detName();
            collectionName.insert(name.c_str());
            debug() << "Inserted collectionName " << collectionName.size()-1 
                    << " = \"" << name << "\"" 
                    << "(" << isite << "," << idet << "): "
                    << "(" << site << "," << detid << ")"
                    << endreq;
        }
    }

    // Set up a fast sim model for PMT (see ExN05).  Note, this is
    // kind of circular, but I find no other way to attach a fast sim
    // to a sens det.

    //const G4LogicalVolume* lv = GiGaVolumeUtils::findLVolume(m_cathodeLogVol);
    IGiGaGeomCnvSvc *geoSvc = 0;
    if (service("GiGaGeo",geoSvc,true).isFailure()) {
        fatal() << "Failed to get GiGa Geometry service" << endreq;
        return StatusCode::FAILURE;        
    }
    std::vector<std::string>::iterator volumeIter, 
      volumeEnd = m_cathodeLogVols.end();
    int PMTii = 1;
    char str[10] = "";

    for(volumeIter = m_cathodeLogVols.begin(); volumeIter != volumeEnd; 
        volumeIter++){
      std::string cathodeLogVol = *volumeIter;
      G4LogicalVolume* lv = geoSvc->g4LVolume( cathodeLogVol );
      if (!lv) {
        error() << "Failed to get " << cathodeLogVol << endreq;
        return StatusCode::FAILURE;
      }
      std::string PMT_region_name = "PMT_region";
      std::string PmtModel_name = "PmtModel";
      sprintf(str,"%d",PMTii);
      PMT_region_name += str;
      PmtModel_name += str;
      PMTii++;

      G4Region* pmt_reg = new G4Region(PMT_region_name.c_str());
      pmt_reg->AddRootLogicalVolume(const_cast<G4LogicalVolume*>(lv));
      
      vector<double> cuts; 
      cuts.push_back(1.0*mm);
      cuts.push_back(1.0*mm);
      cuts.push_back(1.0*mm);
      pmt_reg->SetProductionCuts(new G4ProductionCuts());
      pmt_reg->GetProductionCuts()->SetProductionCuts(cuts);
      
      DsPmtModel* pmtModel = new DsPmtModel(PmtModel_name.c_str(),pmt_reg);
      pmtModel = 0;
    }

    // set up random numbers
    IRndmGenSvc *rgs = 0;
    if (service("RndmGenSvc",rgs,true).isFailure()) {
        fatal() << "Failed to get random service" << endreq;
        return StatusCode::FAILURE;        
    }
    StatusCode sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
    if (sc.isFailure()) {
        fatal() << "Failed to initialize uniform random numbers" << endreq;
        return StatusCode::FAILURE;
    }

    info()<<"Photons prescaling is "<<( m_applyPreQE?"on":"off" )
          <<". Preliminary applied efficiency is "<<m_preQE<<endreq;
    info()<<"Make sure the same values are set for the DsPhysConsOptical"<<endreq;

    return StatusCode::SUCCESS;
}

StatusCode DsPmtSensDet::finalize()
{
    info() << "DsPmtSensDet finalize" << endreq;
    return this->GiGaSensDetBase::finalize();
}


void DsPmtSensDet::Initialize(G4HCofThisEvent* hce)
{
    m_hc.clear();

    G4DhHitCollection* hc = new G4DhHitCollection(SensitiveDetectorName,collectionName[0]);
    m_hc[0] = hc;
    int hcid = G4SDManager::GetSDMpointer()->GetCollectionID(hc);
    hce->AddHitsCollection(hcid,hc);

    for (int isite=0; site_ids[isite] >= 0; ++isite) {
        for (int idet=0; detector_ids[idet] >= 0; ++idet) {
            DayaBay::Detector det(site_ids[isite],detector_ids[idet]);

            if (det.bogus()) continue;

            string name=det.detName();
            G4DhHitCollection* hc = new G4DhHitCollection(SensitiveDetectorName,name.c_str());
            short int id = det.siteDetPackedData();
            m_hc[id] = hc;

            int hcid = G4SDManager::GetSDMpointer()->GetCollectionID(hc);
            hce->AddHitsCollection(hcid,hc);
            debug() << "Add hit collection with hcid=" << hcid << ", cached ID=" 
                    << (void*)id 
                    << " name= \"" << SensitiveDetectorName << "/" << name << "\"" 
                    << endreq;
        }
    }

    debug() << "DsPmtSensDet Initialize, made "
           << hce->GetNumberOfCollections() << " collections"
           << endreq;
    
}

// Return the SensitiveDetector ID for the given touchable history.
const DetectorElement* DsPmtSensDet::SensDetElem(const G4TouchableHistory& hist)
{
    const IDetectorElement* idetelem = 0;
    int steps=0;

    if (!hist.GetHistoryDepth()) {
        error() << "DsPmtSensDet::SensDetElem given empty touchable history" << endreq;
        return 0;
    }

    StatusCode sc = 
        m_t2de->GetBestDetectorElement(&hist,m_sensorStructures,idetelem,steps);
    if (sc.isFailure()) {      // verbose warning
        warning() << "Failed to find detector element in:\n";
        for (size_t ind=0; ind<m_sensorStructures.size(); ++ind) {
            warning() << "\t\t" << m_sensorStructures[ind] << "\n";
        }
        warning() << "\tfor touchable history:\n";
        for (int ind=0; ind < hist.GetHistoryDepth(); ++ind) {
            warning() << "\t (" << ind << ") " 
                      << hist.GetVolume(ind)->GetName() << "\n";
        }
        warning() << endreq;
        return 0;
    }
    
    return dynamic_cast<const DetectorElement*>(idetelem);
}

int  DsPmtSensDet::SensDetId(const DetectorElement& de)
{
    const DetectorElement* detelem = &de;

    while (detelem) {
        if (detelem->params()->exists(m_idParameter)) {
            break;
        }
        detelem = dynamic_cast<const DetectorElement*>(detelem->parentIDetectorElement());
    }
    if (!detelem) {
        warning() << "Could not get PMT detector element starting from " << de << endreq;
        return 0;
    }

    return detelem->params()->param<int>(m_idParameter);
}
double DsPmtSensDet::SensDetQE(G4LogicalVolume* logvol, double energy)
{
    G4Material* mat = logvol->GetMaterial();
    if (!mat) {
        warning () << "No material for " << logvol->GetName() << endreq;
        return -1;
    }
    

     G4MaterialPropertiesTable* mattab = mat->GetMaterialPropertiesTable();
    if (mattab) {
        G4MaterialPropertyVector* qevec = mattab->GetProperty(m_qeffParamName.c_str());
        if (qevec) {
        
          debug() << m_qeffParamName << ":("
                  << qevec->GetMinPhotonEnergy()/CLHEP::eV << " eV," 
                  << qevec->GetMinProperty() << ")-->("
                  << qevec->GetMaxPhotonEnergy()/CLHEP::eV << " eV," 
                  << qevec->GetMaxProperty() << ")"
                  << " particle energy is " << energy/CLHEP::eV
                  << endreq;

          return qevec->GetProperty(energy);

        }
    }
    else {
        debug () << "No material properties in " << logvol->GetName() << endreq;
    }

    int ndaught = logvol->GetNoDaughters();
    for (int ind=0; ind < ndaught; ++ind) {
        G4VPhysicalVolume* physvol = logvol->GetDaughter(ind);
        double qe = this->SensDetQE(physvol->GetLogicalVolume(),energy);
        if (qe < 0) return qe;
    }
    warning() << "All attempts failed to find " << m_qeffParamName 
              << " in " << logvol->GetName() << endreq;
    return -1;
}

bool DsPmtSensDet::ProcessHits(G4Step* step,
                               G4TouchableHistory* /*history*/)
{
    //if (!step) return false; just crash for now if not defined

    // Find out what detector we are in (ADx, IWS or OWS)
    G4StepPoint* preStepPoint = step->GetPreStepPoint();

    double energyDep = step->GetTotalEnergyDeposit();

    if (energyDep <= 0.0) {
        //debug() << "Hit energy too low: " << energyDep/CLHEP::eV << endreq;
        return false;
    }

    const G4TouchableHistory* hist = 
        dynamic_cast<const G4TouchableHistory*>(preStepPoint->GetTouchable());
    if (!hist or !hist->GetHistoryDepth()) {
        error() << "ProcessHits: step has no or empty touchable history" << endreq;
        return false;
    }

    const DetectorElement* de = this->SensDetElem(*hist);
    if (!de) return false;

    // wangzhe QE calculation starts here.
    int pmtid = this->SensDetId(*de);
    DayaBay::Detector detector(pmtid);
    G4Track* track = step->GetTrack();
    double weight = track->GetWeight();

    // Make sure the energy deposit is caused by an optical photon.
    // In the future we may decide that energy deposits in the photocathode by other particles
    // are capable of generating a PMT signal. In that case, we will need a model for the
    // effects. djaffe

    if ( track->GetDefinition() != G4OpticalPhoton::Definition() ) 
      {
        debug() << "Track definition " << track->GetDefinition()->GetParticleName() 
                << "This track is not an optical photon. abort further processing " << endreq;
        return false;
      }

    double qe = this->SensDetQE(preStepPoint->GetPhysicalVolume()->GetLogicalVolume(),energyDep);

        // Scale total efficiency upward to allow for PMT-to-PMT
        // efficiency variation applied in the electronics simulation.
        // The QEScale used here should be the inverse of the mean
        // PMT efficiency applied in the electronics simulation.
        qe *= m_qeScale;
      
        // For now, hard code "extra" decrease in efficiency for water shield PMTs to match G4dyb.
        if (detector.detectorId() == DetectorId::kIWS || detector.detectorId() == DetectorId::kOWS) {
          qe *= 0.6;
        }
    if (qe < 0.) return false;
    if ( m_applyPreQE ) qe/=m_preQE; 
    if (qe > 1.) {
            warning() << "Non physical QE > 1." << endreq;
        }

        // apply QE
        double looser = m_uni();
        // debug()<<"qe= "<<qe<<"   looser= "<<looser<<endreq;
        if (looser > qe) {
          //verbose() << "Failed quantum efficiency " << looser << " > " << qe << endreq;
          return true;
        }
    // wz QE calculation and sampling end here. PE is going to be created next.

    double wavelength = CLHEP::twopi*CLHEP::hbarc/energyDep;

#if 1
    if (!pmtid) {
        warning () << "Dumping TouchableHistory:\n";
        for (int ind=0; ind<hist->GetHistoryDepth(); ++ind) {
            warning () << "\t (" << ind << ") " << hist->GetVolume(ind)->GetName() 
                    << " #" << hist->GetReplicaNumber(ind)
                    << " physvol @ " << (void*)hist->GetVolume(ind)
                    << "\n";
        }
        warning() << endreq;
    }

    verbose() << "Dumping TouchableHistory:\n";
    for (int ind=0; ind<hist->GetHistoryDepth(); ++ind) {
        verbose() << "\t (" << ind << ") " << hist->GetVolume(ind)->GetName() 
                  << " #" << hist->GetReplicaNumber(ind)
                  << " physvol @ " << (void*)hist->GetVolume(ind)
                  << "\n";
    }
    verbose() << endreq;
#endif


    DayaBay::SimPmtHit* sphit = new DayaBay::SimPmtHit();

    // base hit

    // Time since event created
    sphit->setHitTime(preStepPoint->GetGlobalTime());

    //#include "G4NavigationHistory.hh"

    const G4AffineTransform& trans = hist->GetHistory()->GetTopTransform();
    const G4ThreeVector& global_pos = preStepPoint->GetPosition();
    G4ThreeVector pos = trans.TransformPoint(global_pos);
    sphit->setLocalPos(pos);
    sphit->setSensDetId(pmtid);
    
    // pmt hit
    // sphit->setDir(...);       // for now
    G4ThreeVector pol = trans.TransformAxis(track->GetPolarization());
    pol = pol.unit();
    G4ThreeVector dir = trans.TransformAxis(track->GetMomentum());
    dir = dir.unit();
    sphit->setPol(pol);
    sphit->setDir(dir);
    sphit->setWavelength(wavelength);
    sphit->setType(0);
    // G4cerr<<"PMT: set hit weight "<<weight<<G4endl; //gonchar
    sphit->setWeight(weight);

    // Dump info on the hit and touchable history
#if 0
    debug() << "hit: id=" << (void*)pmtid << ", "
            << wavelength/CLHEP::nm << " nm, "
            << "pol=[" << pol[0] << "," << pol[1] << "," << pol[2] << "] "
            << "pos=[" << pos[0] << "," << pos[1] << "," << pos[2] << "] "
            << endreq;

    for (int ind=0; ind<hist->GetHistoryDepth(); ++ind) {
        verbose () << "\t (" << ind << ") " << hist->GetVolume(ind)->GetName() 
                << " #" << hist->GetReplicaNumber(ind)
                << " physvol @ " << (void*)hist->GetVolume(ind)
                << "\n";
    }
    verbose() << endreq;
#endif


    int trackid = track->GetTrackID();
    this->StoreHit(sphit,trackid);

    return true;
}

void DsPmtSensDet::StoreHit(DayaBay::SimPmtHit* hit, int trackid)
{
    int did = hit->sensDetId();
    DayaBay::Detector det(did);
    short int sdid = det.siteDetPackedData();

    G4DhHitCollection* hc = m_hc[sdid];
    if (!hc) {
        warning() << "Got hit with no hit collection.  ID = " << (void*)did
                  << " which is detector: \"" << DayaBay::Detector(did).detName()
                  << "\". Storing to the " << collectionName[0] << " collection"
                  << endreq;
        sdid = 0;
        hc = m_hc[sdid];
    }

#if 1
    verbose() << "Storing hit PMT: " << (void*)did 
              << " from " << DayaBay::Detector(did).detName()
              << " in hc #"<<  sdid << " = "
              << hit->hitTime()/CLHEP::ns << "[ns] "
              << hit->localPos()/CLHEP::cm << "[cm] "
              << hit->wavelength()/CLHEP::nm << "[nm]"
              << endreq;
#endif

    hc->insert(new G4DhHit(hit,trackid));
}

void DsPmtSensDet::EndOfEvent(G4HCofThisEvent* hce)
{
    debug() << "Cache has " << m_hc.size() << " collections" << endreq;

    int ncols = hce->GetNumberOfCollections();
    debug() << "DsPmtSensDet EndOfEvent " << ncols << " collections.";
    //    for (int ind=0; ind<ncols; ++ind) {
    //  G4VHitsCollection* hc = hce->GetHC(ind);
    //  info () << ind << ": " 
    //          << hc->GetSDname() << "//" << hc->GetName() << " has " 
    //          << hc->GetSize() << " hits\n";
    //}
    //info() << endreq;


    int tothits = 0;
    for (int ind=0; ind<ncols; ++ind) {
      G4VHitsCollection* hc = hce->GetHC(ind);
      if ( hc->GetSize() > 0) 
        {
          if ( tothits == 0) debug() << endreq;
          debug() << ind << ": " 
                  << hc->GetSDname() << "//" << hc->GetName() << " has " 
                  << hc->GetSize() << " hits" << endreq;
        }
      tothits += hc->GetSize() ;
    }
    if ( tothits == 0 ) debug() << " No hits found in " << ncols << " collections."  << endreq;
}

---