In This Package:

DsFastMuonStackAction.cc

Go to the documentation of this file.

/*
A tool to weight photons based on the detector element.

Usage:
    from DetSim.DetSimConf import DsFastMuonStackAction
    saction = DsFastMuonStackAction("GiGa.DsFastMuonStackAction")
    saction.Detectors = [ '/dd/Structure/DayaBay/db-rock/db-ows', 
                          '/dd/Structure/DayaBay/db-rock/db-ows/db-curtain/db-iws',
                          '/dd/Structure/DayaBay/db-rock/db-ows/db-curtain/db-iws/db-ade1/db-sst1/db-oil1',
                          '/dd/Structure/DayaBay/db-rock/db-ows/db-curtain/db-iws/db-ade2/db-sst2/db-oil2' ]
    saction.Limits = [ 10000, 10000, 10000, 10000 ]
    saction.10.s = [ 10., 10., 10., 10. ]
*/

#include "DsFastMuonStackAction.h"

#include "DetDesc/IPVolume.h"
#include "DetHelpers/ICoordSysSvc.h"
#include "DetDesc/DetectorElement.h"
#include "DetDesc/IGeometryInfo.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/DeclareFactoryEntries.h"
#include "GaudiKernel/IRndmGenSvc.h"

#include <G4SDManager.hh>
#include <G4RunManager.hh>
#include <G4TrackStatus.hh>
#include <G4ThreeVector.hh>
#include <G4ParticleDefinition.hh>
#include <G4Event.hh>
#include <G4String.hh>
#include <G4ParticleTypes.hh>
#include <G4VTouchable.hh>
#include <G4Track.hh>
#include <G4ClassificationOfNewTrack.hh>

#ifdef DSFMSA_DEBUGPOSITION
#include <TFile.h>
#include <TTree.h>
#endif

DECLARE_TOOL_FACTORY(DsFastMuonStackAction);

DsFastMuonStackAction::DsFastMuonStackAction ( const std::string& type   , 
                                   const std::string& name   , 
                                   const IInterface*  parent ) 
    : GiGaStackActionBase( type, name, parent ) 
    , m_csvc(0)
    , fTrackRemainingPhotons(false)
{ 
    declareProperty("Detectors",  m_detectors, "A list with detector elements you are interested in.");
    declareProperty("Limits",     m_limits,    "A list with maximal number of photons to start weighting "
                                               "for each detector element respectively.");
    declareProperty("Weights",    m_weights,   "A list with photon weights for each detector respectively.");
}

StatusCode DsFastMuonStackAction::initialize() 
{
    if ( m_detectors.size()!=m_limits.size() || m_limits.size()!=m_weights.size() ){
        error() << "Specified parameter lists are inconsistent. Check sizes." << endreq;
        return StatusCode::FAILURE;
    }
    for (unsigned int i(0); i<m_detectors.size(); i++){
        IDetectorElement* de=getDet<IDetectorElement>(m_detectors[i]);;
        if ( !de ) {
            error() << "Detector '"<<m_detectors[i] <<"' not found." << endreq;
            return StatusCode::FAILURE;
        }
#ifdef DSFMSA_DEBUGPOSITION
        tmpMap[de]=int(i)+1;
#endif
        fCounter[de]=DECounter(m_limits[i], m_weights[i]);
        int size=m_limits[i]*sizeof(G4Track)/1024.;
        info()<<"Entry "<<i<<": photons are weighted with w="<<m_weights[i]
                           <<", after N="<<m_limits[i]
                           <<" (up to "<<size<<" KB) in detector element:"<<endreq;
        info() << m_detectors[i] << endreq;
    }

    if ( fCounter.empty() ) {
        error() << "Got no detector elements to look after" << endreq;
        return StatusCode::FAILURE;
    }

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

#ifdef DSFMSA_DEBUGPOSITION
    G4cerr<<"create file"<<G4endl;
    event=0;
    file = new TFile("testCoord.root", "RECREATE");
    tree = new TTree("coordtree", "coordinates tree");
    tree->Branch("x", &xx, "x/D");
    tree->Branch("y", &yy, "y/D");
    tree->Branch("z", &zz, "z/D");
    tree->Branch("i", &ii, "i/I");
    tree->Branch("e", &event, "e/I");
#endif

    return StatusCode::SUCCESS; 
}

StatusCode DsFastMuonStackAction::finalize() 
{
#ifdef DSFMSA_DEBUGPOSITION
  file->Write();
  file->Close();
#endif
  return  GiGaStackActionBase::finalize();
}

G4ClassificationOfNewTrack DsFastMuonStackAction::ClassifyNewTrack (const G4Track* aTrack) {
    //
    // Classify new track, and if it is a photon determine whether it should be weighted or not
    //
    if (aTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return fUrgent;
    //if ( fTrackRemainingPhotons )  info()<<"track remaining"<<endreq; //debug

    // 
    // Determine the whether the track is situated in interesting for us DE or not
    //
    IDetectorElement* de;
    const G4VTouchable* touchable=aTrack->GetTouchable();
    std::map<const G4VTouchable*, IDetectorElement*>::iterator cacheIt=fTouchablesCache.find(touchable);
    std::map<IDetectorElement*, DECounter>::iterator it, end=fCounter.end();
    if ( cacheIt==fTouchablesCache.end() ) {
        // 
        // cached DE is not found
        // determine it and save to cache
        //
        // info()<<"  cached DE is not found"<<endreq;//debug
        Gaudi::XYZPoint position(aTrack->GetPosition().x(),aTrack->GetPosition().y(),aTrack->GetPosition().z());
        IDetectorElement* de=m_csvc->belongsToDE(position);
        it=fCounter.find(de); 
        while ( it==end && de ){ 
            // iterate up through the hierarchy of detector elements, 
            // until we find the DE we are interested in,
            // or at least get to the top of the hierarchy
            de=de->parentIDetectorElement();    
            it=fCounter.find(de); 
        }
        //info()<<"  found "<<de<<" points to"<<(int)de<<" end:"<<(it==end)<<endreq;//debug
        if ( touchable ) fTouchablesCache[touchable]=de;
        else { 
            // G4VPhysicalVolume* v=aTrack->GetVolume();
            // G4String name=v?v->GetName():"no volume";
            // warning()<<"FastMuonStackAction got NULL touchable to '"<<name<<"'. Ignoring..."<<endreq;
            return fUrgent;
        }
        if ( !de ) {// point is not situated in any of DE we are interested in
            return fUrgent;
        }
    }
    else {
        // cached DE found
        // info()<<"  cache found"<<endreq;//debug
        de=cacheIt->second;
        if ( !de ) {
            // point is not situated in any of DE we are interested in
            return fUrgent;
        }
        it=fCounter.find(de);
    }

#ifdef DSFMSA_DEBUGPOSITION
    ii=tmpMap[it->first];
    if (ii==0){
        info()<<"found again strange de"<<endreq;
    }
    xx=aTrack->GetPosition().x();
    yy=aTrack->GetPosition().y();
    zz=aTrack->GetPosition().z();
    tree->Fill();
    if (ii==0) {
        G4VPhysicalVolume* v=aTrack->GetVolume();
        G4String name=v?v->GetName():"no volume";
        info()<<"strange volume"<<it->first->name()<<" "<<name.data()<<endreq;
    }
#endif

    DECounter& counter=it->second;
    counter++;
    G4double tweight=aTrack->GetWeight();
    if ( tweight>=counter.weight ) {
        // ignore the tracks with high weight
        // even if it's logically correct
        // large weights are bad
        return fUrgent;
    }

    if ( !counter.full()  ) {
        //number of photons for this DE is small, continue collecting
        return fTrackRemainingPhotons?fUrgent:fWaiting;
    }
      
    if( m_uni()>=1./counter.weight ) { 
        //info()<<"killed"<<endreq;//debug
        return fKill; 
    }
    
    // info()<<"push weighted track"<<endreq;//debug
    G4Track* track=const_cast<G4Track*>(aTrack); //yes, it's ugly, but should be safe
    track->SetWeight(tweight*counter.weight);
    counter.nTracks++;
    return fUrgent;
}

void DsFastMuonStackAction::NewStage()
{
    // 
    // The Urgent stack is empty. It means that all possible particles are simulated
    // except the photons, delayed fo several DEs.
    // info()<<"New stage, reclassifying"<<endreq;
    std::map<IDetectorElement*, DECounter>::iterator it, end=fCounter.end();
    for ( it=fCounter.begin(); it!=end; it++ ) it->second.reset();
    fTrackRemainingPhotons=true;
    stackManager->ReClassify();
    fTrackRemainingPhotons=false;
    //info()<<"Classification is finished"<<endreq;
}


void DsFastMuonStackAction::PrepareNewEvent()
{
#ifdef DSFMSA_DEBUGPOSITION
    event++; 
#endif
    fTrackRemainingPhotons=false;
    std::map<IDetectorElement*, DECounter>::iterator it, end=fCounter.end();
    for ( it=fCounter.begin(); it!=end; it++ ){
        DECounter& c=it->second;
        if ( c.full() ) 
            info()<<c.nTracks<<" tracks for "<<        
                    c.count-c.maxCount<<" photons are simulated in "<<
                    it->first->name()<<endreq;

        c.reset(true);
    }
}


//-----------------------END----------------------------------

---