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

In This Package:

FilterAlgs.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 #
00003 # Data Filtering algorithms
00004 # 
00005 # Useful algorithms for filtering data written to output file
00006 #
00007 
00008 __all__ = [ "FilterAlg", "CoincidenceFilterAlg" ]
00009 
00010 # Load GaudiPython
00011 from DybPython.DybPythonAlg import DybPythonAlg
00012 from GaudiPython import SUCCESS, FAILURE
00013 from GaudiPython import gbl
00014 import GaudiKernel.SystemOfUnits as units
00015 
00016 # Shortcuts to DayaBay classes
00017 TimeStamp = gbl.TimeStamp
00018 Detector = gbl.DayaBay.Detector
00019 ReadoutHeader = gbl.DayaBay.ReadoutHeader
00020 RecHeader = gbl.DayaBay.RecHeader
00021 RegistrationSequence = gbl.DayaBay.RegistrationSequence
00022 Rndm = gbl.Rndm 
00023 
00024 # Basic Filtering algorithm
00025 class FilterAlg(DybPythonAlg):
00026     "Filter the reconstructed data from a general data file"
00027     def __init__(self,name):
00028         DybPythonAlg.__init__(self,name)
00029         # Properties
00030         self.RegSeqLocation = RegistrationSequence.defaultLocation()
00031         self.ClearStore = True
00032         self.StorePath = ""
00033         self.Prescale = None
00034         self.DetectorNames = []
00035         # Internal variables
00036         self.detectors = []
00037         self.uniform_random = Rndm.Numbers()
00038         return
00039 
00040     def initialize(self):
00041         status = DybPythonAlg.initialize(self)
00042         self.info("initializing")
00043         for detName in self.DetectorNames:
00044             self.detectors.append(Detector(detName))
00045         # Initialize random number service if applying prescale 
00046         if self.Prescale != None:
00047             self.randSvc = self.svc('IRndmGenSvc','RndmGenSvc')
00048             if self.randSvc == None:
00049                 self.error("Failed to initialize random number service.")
00050                 return FAILURE
00051             status= self.uniform_random.initialize(self.randSvc, Rndm.Flat(0,1))
00052             if not status.isSuccess():
00053                 self.error("Failed to initialize uniform random numbers.")
00054                 return status
00055         return status
00056 
00057     def execute(self):
00058         self.info("executing")
00059         evt = self.evtSvc()
00060         # Get Registration Sequence
00061         regSequence = evt[self.RegSeqLocation]
00062         if regSequence == None:
00063             self.error("Failed to retrieve Registration Sequence from"
00064                        + self.regseqLocation)
00065             return FAILURE
00066         # Get List of TES objects registered this execution cycle
00067         registrations = regSequence.registrations()
00068         for registration in registrations:
00069             hdrObj = registration.object()
00070             if hdrObj == None:
00071                 self.error("Failed to retrieve object at "+registration.path())
00072                 return FAILURE
00073             if self.ClearStore:
00074                 # Start by turning off storage of all items
00075                 regSequence.registration(hdrObj).setStore(False)
00076             if len(self.detectors)>0:
00077                 # Filter by detector if a detector list is defined
00078                 currentDetector = Detector(hdrObj.context().GetSite(),
00079                                            hdrObj.context().GetDetId())
00080                 if currentDetector not in self.detectors:
00081                     continue
00082             if registration.path() == self.StorePath:
00083                 # Store if it is at the selected TES path, check prescale
00084                 if self.Prescale == None:
00085                     regSequence.registration(hdrObj).setStore(True)
00086                     self.debug("Storing "+registration.path()+" in output file")
00087                 elif self.Prescale >= self.uniform_random():
00088                     regSequence.registration(hdrObj).setStore(True)
00089                     self.debug("Storing "+registration.path()+" in output file")
00090         return SUCCESS
00091         
00092     def finalize(self):
00093         self.info("finalizing")
00094         status = DybPythonAlg.finalize(self)
00095         return status
00096 
00097 
00098 # Time-correlation Filtering algorithm
00099 class CoincidenceFilterAlg(DybPythonAlg):
00100     "Filter stored data based on time between readouts"
00101     def __init__(self,name):
00102         DybPythonAlg.__init__(self,name)
00103         # Properties
00104         self.ReadoutLocation = ReadoutHeader.defaultLocation()
00105         self.RegSeqLocation = RegistrationSequence.defaultLocation()
00106         self.StorePaths = [ ReadoutHeader.defaultLocation() ]
00107         self.CoincidenceWindow = 1.0 * units.millisecond
00108         self.DetectorNames = []
00109         self.AdOnly = False
00110         self.SameDetectorOnly = False
00111         self.SaveIntermediate = False
00112         self.SaveNeighbors = False
00113         self.NeighborWindow = 1000.0 * units.nanosecond
00114         # Internal variables
00115         self.detectors = []
00116         return
00117 
00118     def initialize(self):
00119         status = DybPythonAlg.initialize(self)
00120         self.info("initializing")
00121         for detName in self.DetectorNames:
00122             self.detectors.append(Detector(detName))
00123         return status
00124 
00125     def execute(self):
00126         self.info("executing")
00127         # Get the readout headers in the archive
00128         readoutArchive = self.getAES(self.ReadoutLocation)
00129         if readoutArchive == None: return FAILURE
00130         if len(readoutArchive)<1: return SUCCESS # No ReadoutHeaders in archive
00131         # Get current readout (first in archive list)
00132         currentReadout = readoutArchive[0].readout()
00133         if currentReadout == None: return SUCCESS # No readout this cycle
00134         currentDet = currentReadout.detector()
00135         if self.AdOnly and not currentDet.isAD(): # Select AD coincidence only
00136             return SUCCESS
00137         if len(self.detectors)>0: # Use list of selected detectors
00138             if currentDet not in self.detectors:
00139                 return SUCCESS  # Detector not in selected list 
00140         # Get previous readouts from archive
00141         previousReadouts = []
00142         for readoutHdr in readoutArchive[1:]:
00143             readout = readoutHdr.readout()
00144             if readout != None: previousReadouts.append(readout)
00145         if len(previousReadouts)<1: return SUCCESS # No preceeding readouts
00146         # Look for coincidences with the current readout
00147         promptReadouts = self.findPromptReadouts(currentReadout,
00148                                                  previousReadouts)
00149         if len(promptReadouts)==0: return SUCCESS # No prompts
00150         # Make the list of readouts we want to store in the output file
00151         storageReadoutHeaders = self.makeStorageList(currentReadout,
00152                                                      promptReadouts,
00153                                                      previousReadouts)
00154         # Mark the readouts (and associated data) for storage
00155         status = self.markForStorage(storageReadoutHeaders)
00156         return status
00157         
00158     def finalize(self):
00159         self.info("finalizing")
00160         status = DybPythonAlg.finalize(self)
00161         return status
00162 
00163     def findPromptReadouts(self, currentReadout, previousReadouts):
00164         # Use the coincidence time to find the prompt readouts
00165         promptReadouts = []
00166         for previousReadout in previousReadouts:
00167             if ( (self.SameDetectorOnly or self.AdOnly) 
00168                 and currentReadout.detector() != previousReadout.detector()):
00169                 continue  # Limited to the same detector?            
00170             if len(self.detectors)>0: # Use list of selected detectors
00171                 if currentReadout.detector() not in self.detectors: continue
00172             # Check time between readouts
00173             deltaTime = TimeStamp(currentReadout.triggerTime().GetSec()
00174                                   - previousReadout.triggerTime().GetSec(),
00175                                   currentReadout.triggerTime().GetNanoSec()
00176                                   - previousReadout.triggerTime().GetNanoSec())
00177             self.verbose("dT readout: "+str(deltaTime.GetSeconds()))
00178             if deltaTime.GetSeconds()*units.second <= self.CoincidenceWindow:
00179                 # Found coincidence
00180                 self.debug("Found coincidence dT [s]: "
00181                            +str(deltaTime.GetSeconds()) )
00182                 promptReadouts.append(previousReadout)
00183         return promptReadouts
00184 
00185     def markForStorage(self, storageHeaders):
00186         # Mark the headers (and associated data) for storage
00187         # First, recursively add input headers to storage list
00188         for header in storageHeaders:
00189             self.addInputHeaders( header, storageHeaders )
00190         # Second, loop over RegistrationSequences, and find header Registrations
00191         registrations = []
00192         status = self.findRegistrations( storageHeaders, registrations )
00193         if not status.isSuccess(): return status
00194         for header, registration in zip( storageHeaders, registrations ):
00195             if registration == None:
00196                 self.error("Failed to find registration for header "
00197                            +header.name())
00198                 return FAILURE
00199             if registration.path() in self.StorePaths:
00200                 self.debug("Storing header at "+registration.path())
00201                 registration.setStore(True) # Mark header for storage
00202         return SUCCESS
00203 
00204     def findRegistrations(self, headers, registrations):
00205         # Find the ObjectReg registrations for each header in list
00206         # Initialize registration list with 'None'
00207         for header in headers:
00208             registrations.append(None)
00209         # Get Registration Archive
00210         regSeqArchive = self.getAES(self.RegSeqLocation)
00211         if regSeqArchive == None: return FAILURE
00212         for regSequence in regSeqArchive:
00213             for registration in regSequence.registrations():
00214                 hdrObj = registration.object()
00215                 if hdrObj == None:
00216                     self.error("Failed to retrieve object at "
00217                                +registration.path())
00218                     return FAILURE
00219                 if hdrObj in headers:
00220                     # Found the registration for this header
00221                     hdrIndex = headers.index(hdrObj)
00222                     registrations[hdrIndex] = regSequence.registration(hdrObj)
00223         return SUCCESS
00224     
00225     def addInputHeaders(self, header, headerList):
00226         # Recurse up the header.inputHeaders and add to headerList
00227         for inputHeader in header.inputHeaders():
00228             if inputHeader not in headerList:
00229                 headerList.append(inputHeader)
00230                 self.addInputHeaders(inputHeader, headerList)
00231         return
00232 
00233     def makeStorageList(self, delayedReadout, promptReadouts, previousReadouts):
00234         # Generate the list of readouts for storage
00235         storageHeaders = [delayedReadout.header()]
00236         storageHeaders += [readout.header() for readout in promptReadouts]
00237         if not self.SaveIntermediate and not self.SaveNeighbors:
00238             return storageHeaders
00239         # Save all readouts between or near prompt/delayed pair
00240         delayedTime = delayedReadout.triggerTime()
00241         earliestTime = TimeStamp(promptReadouts[0].triggerTime())
00242         # First, find earliest prompt readout
00243         for promptReadout in promptReadouts:
00244             if promptReadout.triggerTime() < earliestTime:
00245                 earliestTime = promptReadout.triggerTime()
00246         # Add time window to catch readouts just before prompt?
00247         if self.SaveNeighbors:
00248             earliestTime.Add(TimeStamp(0,
00249                                        int(-self.NeighborWindow
00250                                            /units.nanosecond)))
00251         # Catch intermediate and neigboring readouts
00252         for previousReadout in previousReadouts:
00253             prevTime = previousReadout.triggerTime() 
00254             deltaT_earliest = TimeStamp( prevTime.GetSec()
00255                                          - earliestTime.GetSec(),
00256                                          prevTime.GetNanoSec()
00257                                          - earliestTime.GetNanoSec())
00258             deltaT_delayed = TimeStamp( prevTime.GetSec()
00259                                         - delayedTime.GetSec(),
00260                                         prevTime.GetNanoSec()
00261                                         - delayedTime.GetNanoSec())
00262             if (deltaT_earliest.GetSeconds() >= 0
00263                 and deltaT_delayed.GetSeconds() <= 0):
00264                 if previousReadout.header() not in storageHeaders:
00265                     storageHeaders.append(previousReadout.header())
00266         return storageHeaders
00267 
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:09:43 for DybAlg by doxygen 1.7.4