/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Classes | Public Member Functions | Private Attributes
DsFastMuonStackAction Class Reference

#include <DsFastMuonStackAction.h>

Collaboration diagram for DsFastMuonStackAction:
Collaboration graph
[legend]

List of all members.

Classes

class  DECounter

Public Member Functions

 DsFastMuonStackAction (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~DsFastMuonStackAction ()
virtual StatusCode initialize ()
virtual StatusCode finalize ()
virtual G4ClassificationOfNewTrack ClassifyNewTrack (const G4Track *aTrack)
virtual void NewStage ()
virtual void PrepareNewEvent ()

Private Attributes

ICoordSysSvcm_csvc
bool m_useDynamicWeighting
std::map< IDetectorElement
*, DECounter
fCounter
std::map< const std::string,
IDetectorElement * > 
fTouchablesCache
std::map< IDetectorElement
*, IDetectorElement * > 
fInsideCache
G4bool fTrackRemainingPhotons
std::vector< std::string > m_detectors
std::vector< G4int > m_limits
std::vector< G4double > m_weights
bool m_applyTimeCut
double m_timeCut
double m_maxTimeForWeighting
Rndm::Numbers m_uni
bool m_debugme
std::map< const std::string,
IDetectorElement * >
::size_type 
m_TCsize

Detailed Description

Definition at line 28 of file DsFastMuonStackAction.h.


Constructor & Destructor Documentation

DsFastMuonStackAction::DsFastMuonStackAction ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 47 of file DsFastMuonStackAction.cc.

    : GiGaStackActionBase( type, name, parent ) 
    , m_csvc(0)
    , fTrackRemainingPhotons(false)
    , m_debugme(false) // for debug
    , m_TCsize(0) // for debug
{ 
  info() << "DsFastMuonStackingAction(" << type << "/" << name << ") created" << endreq;
    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.");
    declareProperty("ApplyTimeCut", m_applyTimeCut=false, "Apply time cut");
    declareProperty("TimeCut",    m_timeCut=300.*ns, "Set time cut value (default: 300 ns)");

    declareProperty("UseDynamicWeighting", m_useDynamicWeighting=false, "Use dynamic weighting of photons");
    declareProperty("MaxTimeForWeighting", m_maxTimeForWeighting=1.e16*ns,"Set maximum time cut for weighting (default 1year). Photons with time > cut will NOT be weighted.");
}
virtual DsFastMuonStackAction::~DsFastMuonStackAction ( ) [inline, virtual]

Definition at line 33 of file DsFastMuonStackAction.h.

{};

Member Function Documentation

StatusCode DsFastMuonStackAction::initialize ( ) [virtual]

Definition at line 68 of file DsFastMuonStackAction.cc.

{
    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], m_useDynamicWeighting);
        int size=m_limits[i]*sizeof(G4Track)/1024.;
        if ( m_useDynamicWeighting ) {
          info()<<"Dynamic weighting is applied for entry "<<i<<" after threshold N="<<m_limits[i];
        }
        else {
          info()<<"Entry "<<i<<": photons are weighted with w="<<m_weights[i] <<" after treshold N="<<m_limits[i] ;
        }
        info() <<" (up to "<<size<<" KB) in detector element:"<<endreq;
        info() << m_detectors[i] << endreq;

    }
    if (m_applyTimeCut) info()<<"Optical photons in detector elements of interest with time > "<<m_timeCut/ns
                              <<" ns will be killed rather than weighted, if the number of photons in detector elements exceeds the threshold."<<endreq;

    info() <<"Optical photons with time > " << m_maxTimeForWeighting/ns << " ns will never be weighted or counted towards weighting." << 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 ( ) [virtual]

Definition at line 140 of file DsFastMuonStackAction.cc.

{
#ifdef DSFMSA_DEBUGPOSITION
  file->Write();
  file->Close();
#endif
  return  GiGaStackActionBase::finalize();
}
G4ClassificationOfNewTrack DsFastMuonStackAction::ClassifyNewTrack ( const G4Track *  aTrack) [virtual]

Classify new track, and if it is an optical photon (OP) determine whether it should be weighted or not. Don't consider weighting or counting OP that are past the max time cut. Get and cache the touchable and detector element for each OP. If OP is in DE of interest, then apply acceptance/rejection and weighting as specified by user.

debug: print contents of touchables cache whenever it changes

end code for debugging

debug: compare DE from cache with DE from position. Note mismatches when comparing DE from position with DE from cache, ignore cases where DE from cache is NULL

potential mismatch, but must first check DE hierarchy

end of added debugging

Definition at line 166 of file DsFastMuonStackAction.cc.

                                                                                         {
  //
  if (aTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return fUrgent;

  if ( aTrack->GetGlobalTime() > m_maxTimeForWeighting ) {
    debug() << "track " << aTrack->GetTrackID() << " has time " << aTrack->GetGlobalTime()/ns << " (ns) > maxTimeForWeighting " 
            << m_maxTimeForWeighting/ns << " (ns) so do not weight it. Just put it on urgent stack" <<endreq;
    return fUrgent;
  }

  if ( fTrackRemainingPhotons )  debug()<<"track remaining optical photons"<<endreq; 

  if ( m_debugme ) {
    if ( m_TCsize != fTouchablesCache.size() && !fTouchablesCache.empty() )    {
      debug() << " fTouchablesCache map(touchable,DE) has changed. old " << m_TCsize << " new " << fTouchablesCache.size() << " sizes. key,values follow " << endreq;
      m_TCsize = fTouchablesCache.size();
      for ( std::map<const std::string, IDetectorElement*>::iterator cacheIt = fTouchablesCache.begin(); cacheIt!=fTouchablesCache.end(); cacheIt++) {
        debug() << "touchableName " << cacheIt->first ;
        if (cacheIt->second) debug() << " DE name " << cacheIt->second->name() << endreq;
        else                 debug() << " DE name " << cacheIt->second << endreq;
      }
    }
  }
  

  
  // 
  // Determine the whether the track is situated in interesting DE or not
  //   If there is no touchable, then DE is not interesting.
  //
  const G4TouchableHandle touchable=aTrack->GetTouchableHandle();
  std::string touchableName = GetTouchableName( touchable );
  if ( touchableName == "NULL" ) {
    debug() << "track " << aTrack->GetTrackID() << " TouchableName is " << touchableName << " so put track on urgent stack " << endreq;
    return fUrgent;
  }
  else {
    debug() << "track " << aTrack->GetTrackID() << " TouchableName is " << touchableName << endreq;
  }
  std::map<const std::string, IDetectorElement*>::iterator cacheIt=fTouchablesCache.find(touchableName);
  std::map<IDetectorElement*, DECounter>::iterator it, end=fCounter.end();
  if ( cacheIt==fTouchablesCache.end() ) {
    // 
    // cached DE is not found
    // determine it and save to cache
    //
    
    Gaudi::XYZPoint position(aTrack->GetPosition().x(),aTrack->GetPosition().y(),aTrack->GetPosition().z());
    IDetectorElement* de=m_csvc->belongsToDE(position);
    debug()<<"track " << aTrack->GetTrackID() << " in volume " << aTrack->GetVolume()->GetName() 
           << ". Cached DE is not found. DE of current position is "<< de->name() <<endreq; 
    
    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); 
      if ( m_debugme ) {
        if (de) debug() << "track " << aTrack->GetTrackID() << " iterate up to parent DE " << de->name() ;
        else    debug() << "track " << aTrack->GetTrackID() << " iterate up to parent DE " << de ;
        if ( it==end ) debug() << " This DE NOT found in fCounter " << endreq;
        else           debug() << " This DE found in fCounter " << endreq;
      }
    }


    fTouchablesCache[touchableName]=de;
    if ( m_debugme ) {
      if (de) debug() << "track " << aTrack->GetTrackID() << " add DE " << de->name() << " to cache of touchables " ;
      else    debug() << "track " << aTrack->GetTrackID() << " add DE " << de << " to cache of touchables " ;
      debug() << ". The touchable physical volume name " << touchableName << endreq;
    }
    if ( !de ) {// point is not situated in any of DE we are interested in
      debug() << "track " << aTrack->GetTrackID() << " is not in DE of interest. put it on urgent stack" << endreq;
      return fUrgent;
    }
  }
  else {
    // cached DE found

    debug() << "track " << aTrack->GetTrackID() << " current touchable physical volumne name " << touchableName <<" touchable for cached DE is "<< cacheIt->first << endreq;

    IDetectorElement* de = cacheIt->second;

    if ( m_debugme ) {
      Gaudi::XYZPoint position(aTrack->GetPosition().x(), aTrack->GetPosition().y(), aTrack->GetPosition().z() );
      IDetectorElement* deFromPos = m_csvc->belongsToDE( position );
      
      if ( deFromPos ) debug() << "track " << aTrack->GetTrackID() <<" as check, compare DE from the track position " << deFromPos->name() ;
      else             debug() << "track " << aTrack->GetTrackID() <<" as check, compare DE from the track position " << deFromPos         ;
      if ( de ) debug() << " with DE from cache " << de->name() << endreq;
      else      debug() << " with DE from cache " << de << endreq;
      while ( de && de != deFromPos && deFromPos ) {
        debug() << " iterating up DE hierarchy. current DE from track position " << deFromPos->name() << endreq;
        deFromPos = deFromPos->parentIDetectorElement();
      }
      if ( de && de != deFromPos ) {
        debug() <<  "track " << aTrack->GetTrackID() << " MISMATCH in DE from cache and (iterated) DE from position "  ;
        if ( deFromPos ) debug() << deFromPos->name() << endreq;
        else             debug() << deFromPos      << endreq;
      }
    }

    if ( !de ) {
      // point is not situated in any of DE we are interested in
      debug() << "track " << aTrack->GetTrackID() << ". Found cached DE, but it is not in DE of interest. put it on urgent stack" << endreq;
      return fUrgent;
    }
    it=fCounter.find(de);
  }

  // if DE is not a DE of interest, then the iterator in fCounter will be off the end
  if ( it==end ) {
    debug() <<"track " << aTrack->GetTrackID() << ". Not a DE of interest. fCounter iterator is off the end. put track on urgent stack" << endreq;
    return fUrgent;
  }
  
#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();
  debug()<<"track " << aTrack->GetTrackID() << " weight(before) "<<tweight<< " in " << it->first->name()<< endreq;
  if ( tweight>=counter.weight() ) {
    // ignore the tracks with high weight
    // even if it's logically correct
    // large weights are bad
    debug()<<"track " << aTrack->GetTrackID() << " weight= "<<tweight<<" same or larger than 1/prescale "<<counter.weight()<<" so weight remains unchanged "<<endreq;
    return fUrgent;
  }
  
  if ( !counter.full()  ) {
    //number of photons for this DE is small, continue collecting
    debug() << "track " << aTrack->GetTrackID() << " weight "<<tweight << " has small # of OP for this DE. Continue collecting" <<endreq;
    return fTrackRemainingPhotons?fUrgent:fWaiting;
  }
  
  // kill (weighted) photon if it is too late
  if( m_applyTimeCut && aTrack->GetGlobalTime()<m_timeCut ){
    debug() << "track " << aTrack->GetTrackID() << " weight "<<tweight<< " killed by timecut tTrack "<< aTrack->GetGlobalTime()/ns
            <<" > tCut "<<m_timeCut/ns<<endreq;
    return fKill;
  }

  // select photons for application of weighting
  if( m_uni()>=1./counter.weight() ) { 
    debug()<<"random>="<<1./counter.weight()<<" so track " << aTrack->GetTrackID() << " is killed"<<endreq;
    return fKill; 
  }
  
  debug() << "push weighted track " << aTrack->GetTrackID() <<endreq;
  G4Track* track=const_cast<G4Track*>(aTrack); //yes, it's ugly, but should be safe
  track->SetWeight(tweight*counter.weight());
  debug() << " track  " << aTrack->GetTrackID() << " weight set to " << tweight << " X " << counter.weight() << " should be equal to aTrack->GetWeight() " 
          << aTrack->GetWeight() << " Add to urgent stack" <<endreq;
  counter.nTracks++;
  return fUrgent;
}
void DsFastMuonStackAction::NewStage ( ) [virtual]

Definition at line 352 of file DsFastMuonStackAction.cc.

{
    // 
    // The Urgent stack is empty. It means that all possible particles are simulated
    // except the photons, delayed fo several DEs.
  debug()<<"New stage, reclassifying"<<endreq;
    std::map<IDetectorElement*, DECounter>::iterator it, end=fCounter.end();
    for ( it=fCounter.begin(); it!=end; it++ ) { 
      debug() << " volume "<< it->first->name() 
              << " filled " << it->second.filled 
              << " count " << it->second.count 
              << " nTracks " << it->second.nTracks << ". Will now be reset "<< endreq;
      it->second.reset(); 
    }
    fTrackRemainingPhotons=true;
    stackManager->ReClassify();
    fTrackRemainingPhotons=false;
    debug()<<"Classification is finished"<<endreq;
}
void DsFastMuonStackAction::PrepareNewEvent ( ) [virtual]

Definition at line 373 of file DsFastMuonStackAction.cc.

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

        c.reset(true);
    }
}

Member Data Documentation

Definition at line 44 of file DsFastMuonStackAction.h.

Definition at line 45 of file DsFastMuonStackAction.h.

std::map<IDetectorElement*, DECounter> DsFastMuonStackAction::fCounter [private]

Definition at line 62 of file DsFastMuonStackAction.h.

std::map<const std::string, IDetectorElement*> DsFastMuonStackAction::fTouchablesCache [private]

Definition at line 63 of file DsFastMuonStackAction.h.

std::map<IDetectorElement*, IDetectorElement*> DsFastMuonStackAction::fInsideCache [private]

Definition at line 64 of file DsFastMuonStackAction.h.

Definition at line 65 of file DsFastMuonStackAction.h.

std::vector<std::string> DsFastMuonStackAction::m_detectors [private]

Definition at line 67 of file DsFastMuonStackAction.h.

std::vector<G4int> DsFastMuonStackAction::m_limits [private]

Definition at line 68 of file DsFastMuonStackAction.h.

std::vector<G4double> DsFastMuonStackAction::m_weights [private]

Definition at line 69 of file DsFastMuonStackAction.h.

Definition at line 70 of file DsFastMuonStackAction.h.

Definition at line 71 of file DsFastMuonStackAction.h.

Definition at line 72 of file DsFastMuonStackAction.h.

Rndm::Numbers DsFastMuonStackAction::m_uni [private]

Definition at line 73 of file DsFastMuonStackAction.h.

Definition at line 75 of file DsFastMuonStackAction.h.

std::map<const std::string, IDetectorElement*>::size_type DsFastMuonStackAction::m_TCsize [private]

Definition at line 76 of file DsFastMuonStackAction.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:17:58 for DetSim by doxygen 1.7.4