DsRpcSensDet.cc

Go to the documentation of this file.

#include "DsRpcSensDet.h"
#include "DsRpcModel.h"

#include "Event/SimRpcHit.h"

#include "Conventions/Detectors.h"

#include "GiGaCnv/GiGaVolumeUtils.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>
using namespace std;


DetectorId::DetectorId_t detector_ids1[] = 
{  DetectorId::kRPC,
   (DetectorId::DetectorId_t)-1 };

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


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

    declareProperty("StripLogicalVolume",
                    m_stripLogVol="/dd/Geometry/RPC/lvRPCStrip",
                    "RPC read-out strip 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="RpcID",
                    "The name of the user property holding the RPC ID.");


}

DsRpcSensDet::~DsRpcSensDet()
{
}

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

    m_t2de = tool<ITouchableToDetectorElement>(m_t2deName);

    // define collections.

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

        for (int idet=0; detector_ids1[idet] >= 0; ++idet) {
            DetectorId::DetectorId_t detid = detector_ids1[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 RPC.  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_stripLogVol);
    if (!lv) {
        error() << "+++++++++++++++++++++++++++++++Failed to get " << m_stripLogVol << endreq;
        return StatusCode::FAILURE;
    }

    G4Region* rpc_reg = new G4Region("RPC_region");
    rpc_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);
    rpc_reg->SetProductionCuts(new G4ProductionCuts());
    rpc_reg->GetProductionCuts()->SetProductionCuts(cuts);

    DsRpcModel* rpcModel = new DsRpcModel("RpcModel",rpc_reg);
    rpcModel = 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;
    }

    return StatusCode::SUCCESS;
}

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


void DsRpcSensDet::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_ids1[isite] >= 0; ++isite) {
        for (int idet=0; detector_ids1[idet] >= 0; ++idet) {
            DayaBay::Detector det(site_ids1[isite],detector_ids1[idet]);

            if (det.bogus()) continue;

            string name=det.detName();
            debug()<<"det.detName(): "<<name<<"*********************"<<endreq;  // was info()
            G4DhHitCollection* hc = new G4DhHitCollection(SensitiveDetectorName,name.c_str());
            short int id = det.siteDetPackedData();
            debug()<<"det.siteDetPackedData(): "<<id<<"**************"<<endreq; // was info()
            m_hc[id] = hc;

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

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

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

    if (!hist.GetHistoryDepth()) {
        error() << "DsRpcSensDet::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  DsRpcSensDet::SensDetId(const DetectorElement& de)
{
    const DetectorElement* detelem = &de;
    //info()<<"+++++++++++++++++++DetectorElement+++++++++"<<detelem<<endreq; 
    while (detelem) {
       // info()<<"detelem = dynamic_cast<const DetectorElement*>(detelem->parentIDetectorElement())::::::::"<<detelem<<endreq;
        if (detelem->params()->exists(m_idParameter)) {
            break;
        }
        detelem = dynamic_cast<const DetectorElement*>(detelem->parentIDetectorElement());
    }
    if (!detelem) {
        warning() << "Could not get RPC detector element starting from " << de << endreq;
        return 0;
    }

    return detelem->params()->param<int>(m_idParameter);
}


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

    // Find out what detector we are in (RPC)
    G4StepPoint* preStepPoint = step->GetPreStepPoint();

    double energyDep = step->GetTotalEnergyDeposit();
 
    // Do not process non-positive energy deposits 
    if (energyDep <= 0.0) {
        debug()<< "Hit energy too low: " << energyDep/CLHEP::eV << endreq;
        return true;
    }




    const G4TouchableHistory* hist = 
        dynamic_cast<const G4TouchableHistory*>(preStepPoint->GetTouchable());
    if (!hist or !hist->GetHistoryDepth()) {
      error() << "ProcessHits: step has no or empty touchable history." 
              << " Particle name for track assoc. with this step is " << step->GetTrack()->GetDefinition()->GetParticleName() << endreq;
        return false;
    }

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

    int rpcid = this->SensDetId(*de);
    DayaBay::Detector detector(rpcid);
 
    G4Track* track = step->GetTrack();
    double weight = track->GetWeight();
    weight = 1.;


#if 1
    if (!rpcid) {
        warning () << "Dumping TouchableHistory: (Particle name for track assoc. with this step is " 
                   << step->GetTrack()->GetDefinition()->GetParticleName() <<" .)\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::SimRpcHit* srhit = new DayaBay::SimRpcHit();
    srhit->setEnergyDep(energyDep);

    // base hit

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

    //#include "G4NavigationHistory.hh"

    const G4AffineTransform& trans = hist->GetHistory()->GetTopTransform();
    const G4ThreeVector& global_pos = preStepPoint->GetPosition();
    G4ThreeVector pos = trans.TransformPoint(global_pos);
    srhit->setLocalPos(pos);
    srhit->setSensDetId(rpcid);
    
    // rpc hit
    // sprit->setDir(...);       // for now
    //G4ThreeVector pol = trans.TransformAxis(track->GetPolarization());
   // pol = pol.unit();
   // G4ThreeVector dir = trans.TransformAxis(track->GetMomentum());
    //dir = dir.unit();
    //srhit->setPol(pol);
    //srhit->setDir(dir);
    //srhit->setWavelength(wavelength);
    //srhit->setType(0);
    //srhit->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(srhit,trackid);

    return true;
}

void DsRpcSensDet::StoreHit(DayaBay::SimRpcHit* 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 RPC: " << (void*)did 
              << " from " << DayaBay::Detector(did).detName()
              << " in hc #"<<  sdid << " = "
              << hit->hitTime()/CLHEP::ns << "[ns] "
              << hit->localPos()/CLHEP::cm << "[cm] "
              <<"Storing to the "<<collectionName[0]<<"collection****************"
              << endreq;
#endif

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

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

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