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

In This Package:

TestCoinc.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 #
00003 # Simple toy MC for testing IBD coincidence selection
00004 #
00005 # This script generates fake reconstructed data appropriate for input
00006 # to the muon tagging and coincidence algorithms.  After coincidence
00007 # identification and selection, the output can be cross-checked
00008 # against truth information.
00009 #
00010 #  Usage:
00011 #   nuwa.py -n 10000 -m"TestCoinc"
00012 #                    -m"MuonTagging.MuonTagLoose" \
00013 #                    -m"MuonTagging.MuonData" \
00014 #                    -m"ADCoincTagging.CoincidenceTight" \
00015 #                    -m"SimpleFilter.Keep /Event/Data/Physics/CoincidenceTight"\
00016 #                    -o coincOutput.root
00017 #
00018 #
00019 
00020 # Load DybPython
00021 from DybPython.DybPythonAlg import DybPythonAlg
00022 from GaudiPython import SUCCESS, FAILURE
00023 from GaudiPython import gbl
00024 from DybPython.Util import irange
00025 import GaudiKernel.SystemOfUnits as units
00026 
00027 # Loosen python garbage collection
00028 gbl.SetMemoryPolicy( gbl.kMemoryHeuristics )
00029 
00030 # Make shortcuts to any ROOT classes you want to use
00031 TimeStamp = gbl.TimeStamp
00032 Detector = gbl.DayaBay.Detector
00033 Site = gbl.Site
00034 DetectorId = gbl.DetectorId
00035 RecHeader = gbl.DayaBay.RecHeader
00036 rhe = RecHeader()
00037 TMath_Gaus = gbl.TMath.Gaus
00038 
00039 # Decorate TimeStamp class with proper comparison
00040 def compareTimeStamp(tsA, tsB):
00041     "Comparison function for timestamps"
00042     if tsA.GetSec() < tsB.GetSec():
00043         return -1
00044     elif tsA.GetSec() > tsB.GetSec():
00045         return 1
00046     else:
00047         if tsA.GetNanoSec() < tsB.GetNanoSec():
00048             return -1
00049         elif tsA.GetNanoSec() > tsB.GetNanoSec():
00050             return 1
00051     return 0
00052 
00053 TimeStamp.__cmp__ = compareTimeStamp
00054 
00055 # Detector Names
00056 detectorNames = {DetectorId.kAD1:'AD1',
00057                  DetectorId.kAD2:'AD2',
00058                  DetectorId.kAD3:'AD3',
00059                  DetectorId.kAD4:'AD4',
00060                  DetectorId.kIWS:'IWS',
00061                  DetectorId.kOWS:'OWS',
00062                  DetectorId.kRPC:'RPC'
00063                  }
00064 
00065 # Define source IDs
00066 sourceId = {"Unknown":0,
00067             "IBDPrompt":1,
00068             "IBDDelayed":2,
00069             "ADSingles":3,
00070             "IWS":4,
00071             "OWS":5,
00072             "OWS_IWS":6,
00073             "AD1_OWS_IWS":7,
00074             "AD2_OWS_IWS":8,
00075             "AD3_OWS_IWS":9,
00076             "AD4_OWS_IWS":10
00077             }
00078 
00079 class EventSource:
00080     "Source of Fake Events"
00081     def __init__(self,name,rate):
00082         self.name = name
00083         self.rate_Hz = rate
00084         self.site = gbl.Site.kDayaBay
00085         self.type = "Unknown"
00086         self.sourceTime = TimeStamp()
00087         self.triggers = []
00088         self.count = 0
00089         self.detectorLatency = {DetectorId.kAD1:0,
00090                                 DetectorId.kAD2:0,
00091                                 DetectorId.kAD3:0,
00092                                 DetectorId.kAD4:0,
00093                                 DetectorId.kIWS:100.0e-9,
00094                                 DetectorId.kOWS:150.0e-9,
00095                                 DetectorId.kRPC:50.0e-9
00096                                 }
00097         return
00098 
00099     def getNextTrigger(self):
00100         "Return the next trigger from this source"
00101         if len(self.triggers)<1:
00102             print "Error: No triggers from source",self.name
00103             return None
00104         nextTrigger = self.triggers[0]
00105         self.triggers = self.triggers[1:]
00106         return nextTrigger
00107 
00108     def currentTime(self):
00109         "Return the current source internal clock time"
00110         return self.sourceTime
00111 
00112     def queuedTriggers(self):
00113         "Return the number of triggers in the queue for this source"
00114         return len(self.triggers)
00115 
00116     def nextTriggerTime(self):
00117         "Return the time of the next trigger from this source"
00118         if len(self.triggers)<1:
00119             print "Warning: Requesting triggers from source, but none exist."
00120             return TimeStamp.GetEOT()
00121         return self.triggers[0]['triggerTime']
00122 
00123     def makeNextEvent(self):
00124         "Generate the next event from this source"
00125         nextEventTime = self.makeNextEventTime()
00126         status = self.makeOneEventAtTime(nextEventTime)
00127         if status.isFailure(): return status
00128         self.sortTriggers()
00129         self.sourceTime = nextEventTime
00130         self.count += 1
00131         return SUCCESS
00132 
00133     def sortTriggers(self):
00134         "Sort the current trigger by trigger time"
00135         if len(self.triggers)<1: return
00136         self.triggers.sort(key=lambda trigger: trigger['triggerTime'])
00137         return
00138 
00139     def makeNextEventTime(self):
00140         "Generate the time for the next event from this source"
00141         import random
00142         deltaT_s = random.expovariate(self.rate_Hz)
00143         nextTime = TimeStamp(self.sourceTime)
00144         nextTime.Add(deltaT_s)
00145         return nextTime
00146 
00147     def makeOneEventAtTime(self,time):
00148         "Generate one event at the specified time"
00149         print "Error: This is an abstract function!  Why are you calling this?"
00150         return FAILURE
00151 
00152     def makeCalibHeader(self,inputData):
00153         "Make a Calibrated data header given input data"
00154         calibHeader = gbl.DayaBay.CalibReadoutHeader()
00155         gbl.SetOwnership(calibHeader,False)
00156         calibHeader.setContext(inputData['context'])
00157         calibHeader.setExecNumber(inputData['execNumber'])
00158         calibHeader.setEarliest(inputData['earliest'])
00159         calibHeader.setLatest(inputData['latest'])
00160         detectorId = inputData['context'].GetDetId()
00161         calibReadout = None
00162         if detectorId == DetectorId.kRPC:
00163             calibReadout = gbl.DayaBay.CalibReadoutRpcCrate()
00164             gbl.SetOwnership(calibReadout,False)
00165         elif detectorId>=DetectorId.kAD1 and detectorId<=DetectorId.kOWS:
00166             calibReadout = gbl.DayaBay.CalibReadoutPmtCrate()
00167             gbl.SetOwnership(calibReadout,False)
00168         else:
00169             print "Error: Making calib header for an invalid detector:", detectorId
00170             return None
00171         calibReadout.setTriggerType(inputData['triggerType'])
00172         calibHeader.setCalibReadout(calibReadout)
00173         return calibHeader
00174 
00175     def makeCalibStats(self,inputData):
00176         "Make a CalibStats data header given input data"
00177         calibStats = gbl.DayaBay.UserDataHeader()
00178         gbl.SetOwnership(calibStats,False)
00179         calibStats.setContext(inputData['context'])
00180         calibStats.setExecNumber(inputData['execNumber'])
00181         calibStats.setEarliest(inputData['earliest'])
00182         calibStats.setLatest(inputData['latest'])
00183         for inputHdr in inputData['inputHeaders']:
00184             calibStats.addInputHeader(inputHdr)
00185         detectorId = inputData['context'].GetDetId()
00186         defaultIntVars = {'nHit':-1}
00187         defaultFloatVars = {'nPESum':-1,
00188                             'nPERMS':-1,
00189                             'nPulseSum':-1,
00190                             'nPulseRMS':-1,
00191                             'tEarliest':-1,
00192                             'tLatest':-1,
00193                             'tMean':-1,
00194                             'EarlyCharge':-1,
00195                             'NominalCharge':1,
00196                             'LateCharge':-1,
00197                             'MaxQ':0,
00198                             'Quadrant':0,
00199                             'flasher_flag':0,
00200                             'time_PSD':0,
00201                             'time_PSD1':0,
00202                             'MaxQ_2inchPMT':0,
00203                             'dtLastAD1_ms':-1,
00204                             'dtLastAD2_ms':-1,
00205                             'dtLastAD3_ms':-1,
00206                             'dtLastAD4_ms':-1,
00207                             'dtLastIWS_ms':-1,
00208                             'dtLastOWS_ms':-1,
00209                             'dtLastRPC_ms':-1,
00210                             'dtNextAD1_ms':-1,
00211                             'dtNextAD2_ms':-1,
00212                             'dtNextAD3_ms':-1,
00213                             'dtNextAD4_ms':-1,
00214                             'dtNextIWS_ms':-1,
00215                             'dtNextOWS_ms':-1,
00216                             'dtNextRPC_ms':-1
00217                             }
00218         for varName in defaultIntVars:
00219             if varName in inputData:
00220                 calibStats.setInt(varName,inputData[varName])
00221             else:
00222                 calibStats.setInt(varName,defaultIntVars[varName])
00223         for varName in defaultFloatVars:
00224             if varName in inputData:
00225                 calibStats.setFloat(varName,inputData[varName])
00226             else:
00227                 calibStats.setFloat(varName,defaultFloatVars[varName])
00228         return calibStats
00229 
00230     def makeRecHeader(self,inputData):
00231         "Make a Reconstructed data header given input data"
00232         recHeader = gbl.DayaBay.RecHeader()
00233         gbl.SetOwnership(recHeader,False)
00234         recHeader.setContext(inputData['context'])
00235         recHeader.setExecNumber(inputData['execNumber'])
00236         recHeader.setEarliest(inputData['earliest'])
00237         recHeader.setLatest(inputData['latest'])
00238         for inputHdr in inputData['inputHeaders']:
00239             recHeader.addInputHeader(inputHdr)
00240         defaultVars = {'triggerType':0,
00241                        'energy':0,
00242                        'rawEvis':0,
00243                        'x':0,
00244                        'y':0,
00245                        'z':0,
00246                        'energyStatus':1,
00247                        'positionStatus':1
00248                        }
00249         for varName in defaultVars:
00250             if varName in inputData:
00251                 defaultVars[varName] = inputData[varName]
00252         detectorId = inputData['context'].GetDetId()
00253         if detectorId == DetectorId.kRPC:
00254             pass
00255         elif detectorId>=DetectorId.kAD1 and detectorId<=DetectorId.kAD4:
00256             recTrigger = gbl.DayaBay.RecTrigger()
00257             recTrigger.setTriggerType(defaultVars['triggerType'])
00258             recTrigger.setEnergy(defaultVars['energy'])
00259             recTrigger.setEnergyStatus(defaultVars['energyStatus'])
00260             recTrigger.setRawEvis(defaultVars['rawEvis'])
00261             position = gbl.CLHEP.HepLorentzVector(defaultVars['x'],
00262                                                   defaultVars['y'],
00263                                                   defaultVars['z'])
00264             recTrigger.setPosition(position)
00265             recTrigger.setPositionStatus(defaultVars['positionStatus'])
00266             recHeader.setRecTrigger(recTrigger)
00267         elif detectorId>=DetectorId.kIWS and detectorId<=DetectorId.kOWS:
00268             pass
00269         else:
00270             print "Error: Making RecHeader for an invalid detector:", detectorId
00271             return None
00272         return recHeader
00273 
00274     def constructTrigger(self,inputData):
00275         "Construct the data objects needed for one trigger"
00276         triggerData = {}
00277         triggerData['triggerTime'] = inputData['time']
00278         triggerData['detectorId'] = inputData['detectorId']
00279         triggerData['site'] = inputData['site']
00280         headerData = {}
00281         # Calib Readout
00282         calibReadout = self.makeCalibHeader(inputData['calibHeader'])
00283         headerData['/Event/CalibReadout/CalibReadoutHeader'] = calibReadout
00284         # Calib Stats
00285         inputData['calibStats']['inputHeaders'] = [calibReadout]
00286         calibStats = self.makeCalibStats(inputData['calibStats'])
00287         headerData['/Event/Data/CalibStats'] = calibStats
00288         # Finish
00289         triggerData['headers'] = headerData
00290         return triggerData
00291 
00292     def constructAdTrigger(self,inputData):
00293         "Construct the data objects needed for one AD trigger"
00294         triggerData = self.constructTrigger(inputData)
00295         headerData = triggerData['headers']
00296         calibReadout = headerData['/Event/CalibReadout/CalibReadoutHeader']
00297         calibStats = headerData['/Event/Data/CalibStats']
00298         # AD Simple
00299         inputData['adSimple']['inputHeaders'] = [calibReadout,calibStats]
00300         adSimple = self.makeRecHeader(inputData['adSimple'])
00301         headerData['/Event/Rec/AdSimple'] = adSimple
00302         # AD Scaled
00303         inputData['adScaled']['inputHeaders'] = [calibReadout,calibStats]
00304         adScaled = self.makeRecHeader(inputData['adScaled'])
00305         headerData['/Event/Rec/AdScaled'] = adScaled
00306         # Finish
00307         return triggerData
00308 
00309     def constructPoolTrigger(self,inputData):
00310         "Construct the data objects needed for one AD trigger"
00311         triggerData = self.constructTrigger(inputData)
00312         headerData = triggerData['headers']
00313         calibReadout = headerData['/Event/CalibReadout/CalibReadoutHeader']
00314         calibStats = headerData['/Event/Data/CalibStats']
00315         # AD Simple
00316         inputData['adSimple']['inputHeaders'] = [calibReadout,calibStats]
00317         adSimple = self.makeRecHeader(inputData['adSimple'])
00318         headerData['/Event/Rec/AdSimple'] = adSimple
00319         # AD Scaled
00320         inputData['adScaled']['inputHeaders'] = [calibReadout,calibStats]
00321         adScaled = self.makeRecHeader(inputData['adScaled'])
00322         headerData['/Event/Rec/AdScaled'] = adScaled
00323         return triggerData
00324 
00325     def makeEmptyTrigger(self, time, site, detectorId):
00326         "Make an empty skeleton for AD trigger data"
00327         inputData = {}
00328         triggerTime = TimeStamp(time)
00329         triggerTime.Add(self.detectorLatency[detectorId])
00330         inputData['time']=triggerTime
00331         inputData['detectorId']=detectorId
00332         inputData['site']=site
00333         context = gbl.Context(site, 0, triggerTime, detectorId)
00334         earliest = TimeStamp(time)
00335         earliest.Add(-1.0e-6) # subtract 1us
00336         latest = TimeStamp(time)
00337         latest.Add(1.0e-6) # add 1us
00338 
00339         inputCalibHdr = {"context":context,
00340                          "earliest":earliest,
00341                          "latest":latest,
00342                          "execNumber":0,
00343                          "triggerType":0
00344                          }
00345 
00346         inputCalibStats = {"context":context,
00347                            "earliest":earliest,
00348                            "latest":latest,
00349                            "execNumber":0
00350                            }
00351         inputData['calibHeader'] = inputCalibHdr
00352         inputData['calibStats'] = inputCalibStats
00353         return inputData
00354 
00355     def makeEmptyAdTrigger(self, time, site, detectorId):
00356         "Make an empty skeleton for AD trigger data"
00357         inputData = self.makeEmptyTrigger(time, site, detectorId)
00358         inputAdSimple = {"context":inputData['calibHeader']['context'],
00359                          "earliest":inputData['calibHeader']['earliest'],
00360                          "latest":inputData['calibHeader']['latest'],
00361                          "execNumber":0,
00362                          "triggerType":0
00363                          }
00364 
00365         inputAdScaled = {"context":inputData['calibHeader']['context'],
00366                          "earliest":inputData['calibHeader']['earliest'],
00367                          "latest":inputData['calibHeader']['latest'],
00368                          "execNumber":0,
00369                          "triggerType":0
00370                          }
00371 
00372         inputData['adSimple'] = inputAdSimple
00373         inputData['adScaled'] = inputAdScaled
00374         return inputData
00375 
00376     def makeEmptyPoolTrigger(self, time, site, detectorId):
00377         "Make an empty skeleton for Water Pool trigger data"
00378         inputData = self.makeEmptyTrigger(time, site, detectorId)
00379         # Still need to make empty reconstructed data headers for Pool events
00380         inputAdSimple = {"context":inputData['calibHeader']['context'],
00381                          "earliest":inputData['calibHeader']['earliest'],
00382                          "latest":inputData['calibHeader']['latest'],
00383                          "execNumber":0,
00384                          "triggerType":0
00385                          }
00386 
00387         inputAdScaled = {"context":inputData['calibHeader']['context'],
00388                          "earliest":inputData['calibHeader']['earliest'],
00389                          "latest":inputData['calibHeader']['latest'],
00390                          "execNumber":0,
00391                          "triggerType":0
00392                          }
00393 
00394         inputData['adSimple'] = inputAdSimple
00395         inputData['adScaled'] = inputAdScaled
00396         return inputData
00397     
00398 class IBDSource(EventSource):
00399     "Generator for IBD triggers"""
00400     def __init__(self,name,rate):
00401         EventSource.__init__(self,name,rate)
00402         self.type = "IBD"
00403         self.detectorId = gbl.DetectorId.kAD1
00404         self.spillInFrac = 0.05
00405         self.captureTime = 30.0e-6
00406         self.captureTimeSpillIn = 200.0e-6
00407         self.thermalizationTime = 3.0e-6
00408         return
00409 
00410     def makeOneEventAtTime(self,time):
00411         "Generate one event at the specified time"
00412         dtIBD_s = self.makePromptDelayedDt()
00413         promptTrigger = self.makePromptTrigger(time)
00414         if not promptTrigger:
00415             print "Error: Failed to generate prompt IBD trigger"
00416             return FAILURE
00417         self.triggers.append(promptTrigger)
00418         delayedTime = TimeStamp(time)
00419         delayedTime.Add(dtIBD_s)
00420         delayedTrigger = self.makeDelayedTrigger(delayedTime)
00421         self.triggers.append(delayedTrigger)
00422         return SUCCESS
00423 
00424     def makePromptDelayedDt(self):
00425         "Generate the time between prompt/delayed triggers"
00426         # FIXME: This is very approximate; please improve based on data
00427         import random
00428         captureTime = self.captureTime
00429         if random.random() < self.spillInFrac:
00430             captureTime = self.captureTimeSpillIn
00431         deltaT_s = 0
00432         while True:
00433             deltaT_s = random.expovariate(1.0/captureTime)
00434             # Suppress short capture times to account for thermalization time
00435             if deltaT_s > (4*self.thermalizationTime): break
00436             suppressionFactor = 1-TMath_Gaus(deltaT_s,0,
00437                                              self.thermalizationTime)
00438             if random.random() < suppressionFactor:
00439                 break
00440         return deltaT_s
00441 
00442     def makePromptTrigger(self, time):
00443         "Make the data for one prompt IBD trigger"
00444         inputData = self.makeEmptyAdTrigger(time, self.site, self.detectorId)
00445         # Make event-specific modifications here
00446         inputData['adSimple']['energy'] = 3.0
00447         inputData['adSimple']['rawEvis'] = 3.0
00448         inputData['adScaled']['energy'] = 3.0
00449         inputData['adScaled']['rawEvis'] = 3.0
00450         # Hack: Put truth info into one of the unused variables... 
00451         inputData['calibStats']['nPulseSum'] = sourceId["IBDPrompt"]
00452         inputData['calibStats']['nPulseRMS'] = self.count 
00453         triggerData = self.constructAdTrigger(inputData)
00454         return triggerData
00455 
00456     def makeDelayedTrigger(self, time):
00457         "Make the data for one delayed IBD trigger"
00458         inputData = self.makeEmptyAdTrigger(time, self.site, self.detectorId)
00459         # Make special modifications here
00460         inputData['adSimple']['energy'] = 8.0
00461         inputData['adSimple']['rawEvis'] = 8.0
00462         inputData['adScaled']['energy'] = 8.0
00463         inputData['adScaled']['rawEvis'] = 8.0
00464         # Hack: Put truth info into one of the unused variables... 
00465         inputData['calibStats']['nPulseSum'] = sourceId["IBDDelayed"]
00466         inputData['calibStats']['nPulseRMS'] = self.count 
00467         triggerData = self.constructAdTrigger(inputData)
00468         return triggerData
00469 
00470 
00471 class ADSinglesSource(EventSource):
00472     "Generator for AD Singles triggers"""
00473     def __init__(self,name,rate):
00474         EventSource.__init__(self,name,rate)
00475         self.type = "ADSingles"
00476         self.detectorId = gbl.DetectorId.kAD1
00477         self.probabilityHighEnergy = 0.0001
00478         return
00479 
00480     def makeOneEventAtTime(self,time):
00481         "Generate one event at the specified time"
00482         inputData = self.makeEmptyAdTrigger(time, self.site, self.detectorId)
00483         # Make event-specific modifications here
00484         inputData['adSimple']['energy'] = 1.0
00485         inputData['adSimple']['rawEvis'] = 1.0
00486         inputData['adScaled']['energy'] = 1.0
00487         inputData['adScaled']['rawEvis'] = 1.0
00488         # Randomly make a high-energy singles
00489         import random
00490         if random.random() < self.probabilityHighEnergy:
00491             # Make a singles that satisfies delayed energy cut
00492             inputData['adSimple']['energy'] = 7.0
00493             inputData['adSimple']['rawEvis'] = 7.0
00494             inputData['adScaled']['energy'] = 7.0
00495             inputData['adScaled']['rawEvis'] = 7.0
00496         # Hack: Put truth info into one of the unused variables... 
00497         inputData['calibStats']['nPulseSum'] = sourceId["ADSingles"]
00498         inputData['calibStats']['nPulseRMS'] = self.count
00499         triggerData = self.constructAdTrigger(inputData)
00500         self.triggers.append(triggerData)
00501         return SUCCESS
00502 
00503 # FIXME: Muon source only generates prompt muon triggers
00504 #        (does not make retriggers, mu-decay, mu-capture, spallation products) 
00505 class MuonSource(EventSource):
00506     "Generator for Muon triggers"""
00507     def __init__(self,name,rate):
00508         EventSource.__init__(self,name,rate)
00509         self.type = "Muon"
00510         # FIXME: improve muon combinations and relative fractions
00511         self.combinations = {}
00512         self.combinations['OWS'] = [DetectorId.kOWS]
00513         self.combinations['IWS'] = [DetectorId.kIWS]
00514         self.combinations['OWS_IWS'] = [DetectorId.kOWS,
00515                                         DetectorId.kIWS]
00516         self.combinations['AD1_OWS_IWS'] = [DetectorId.kAD1,
00517                                             DetectorId.kIWS,
00518                                             DetectorId.kOWS]
00519         self.combinations['AD2_OWS_IWS'] = [DetectorId.kAD2,
00520                                             DetectorId.kIWS,
00521                                             DetectorId.kOWS]
00522         self.combinations['AD3_OWS_IWS'] = [DetectorId.kAD3,
00523                                             DetectorId.kIWS,
00524                                             DetectorId.kOWS]
00525         self.combinations['AD4_OWS_IWS'] = [DetectorId.kAD4,
00526                                             DetectorId.kIWS,
00527                                             DetectorId.kOWS]
00528         self.nearSiteFractions = {'OWS_IWS':0.66,
00529                                   'AD1_OWS_IWS':0.12,
00530                                   'AD2_OWS_IWS':0.12,
00531                                   'OWS':0.05,
00532                                   'IWS':0.05
00533                                   }
00534         # FIXME: Set proportions
00535         self.farSiteFractions = {'OWS_IWS':1.0,
00536                                  'AD1_OWS_IWS':0,
00537                                  'AD2_OWS_IWS':0,
00538                                  'AD3_OWS_IWS':0,
00539                                  'AD4_OWS_IWS':0,
00540                                  'OWS':0,
00541                                  'IWS':0
00542                                  }
00543         return
00544 
00545     def makeOneEventAtTime(self,time):
00546         "Generate one event at the specified time"
00547         muonType = self.pickMuonType()
00548         triggeredDetectors = self.combinations[muonType]
00549         for detectorId in triggeredDetectors:
00550             if DetectorId.kAD1 <= detectorId <= DetectorId.kAD4:
00551                 inputData = self.makeEmptyAdTrigger(time, self.site,
00552                                                     detectorId)
00553                 # Make event-specific modifications here
00554                 inputData['calibStats']['nPESum'] = 1.0e3 * 160
00555                 inputData['calibStats']['NominalCharge'] = 1.0e3 * 160
00556                 inputData['adSimple']['energy'] = 1.0e3
00557                 inputData['adSimple']['rawEvis'] = 1.0e3
00558                 inputData['adScaled']['energy'] = 1.0e3
00559                 inputData['adScaled']['rawEvis'] = 1.0e3
00560                 # Hack: Put truth info into one of the unused variables... 
00561                 inputData['calibStats']['nPulseSum'] = sourceId[muonType]
00562                 inputData['calibStats']['nPulseRMS'] = self.count
00563                 triggerData = self.constructAdTrigger(inputData)
00564                 self.triggers.append(triggerData)
00565             elif DetectorId.kIWS <= detectorId <= DetectorId.kOWS:
00566                 inputData = self.makeEmptyPoolTrigger(time, self.site,
00567                                                       detectorId)
00568                 # Make event-specific modifications here
00569                 inputData['calibStats']['nHit'] = 50
00570                 # Hack: Put truth info into one of the unused variables... 
00571                 inputData['calibStats']['nPulseSum'] = sourceId[muonType]
00572                 inputData['calibStats']['nPulseRMS'] = self.count
00573                 triggerData = self.constructPoolTrigger(inputData)
00574                 self.triggers.append(triggerData)
00575             else:
00576                 print "Error: Unhandled detector type:",detectorId
00577                 return FAILURE
00578         return SUCCESS
00579 
00580     def pickMuonType(self):
00581         "Pick a random muon type"
00582         import random
00583         muonFractions = self.nearSiteFractions
00584         if self.site == Site.kFar:
00585             muonFractions = self.farSiteFractions
00586         curVal = 0
00587         selectedType = None
00588         randVal = random.random()
00589         for muonType in muonFractions:
00590             selectedType = muonType
00591             curVal += muonFractions[muonType]
00592             if randVal < curVal: break
00593         return selectedType
00594 
00595 
00596 # Make your algorithm
00597 class TestCoincAlg(DybPythonAlg):
00598     "Coincidence Testing Algorithm"
00599     def __init__(self,name):
00600         DybPythonAlg.__init__(self,name)
00601         self.sources = []
00602         self.execNumber = 0
00603         self.minRetriggerTime = 750.0e-9 # minimum retrigger time in nanoseconds
00604         self.startTime = TimeStamp() # Set default start time to 'now'
00605         self.lastTriggerTime = {DetectorId.kAD1:TimeStamp.GetBOT(),
00606                                 DetectorId.kAD2:TimeStamp.GetBOT(),
00607                                 DetectorId.kAD3:TimeStamp.GetBOT(),
00608                                 DetectorId.kAD4:TimeStamp.GetBOT(),
00609                                 DetectorId.kIWS:TimeStamp.GetBOT(),
00610                                 DetectorId.kOWS:TimeStamp.GetBOT(),
00611                                 DetectorId.kRPC:TimeStamp.GetBOT()
00612                                 }
00613         self.triggerCount = {DetectorId.kAD1:0,
00614                              DetectorId.kAD2:0,
00615                              DetectorId.kAD3:0,
00616                              DetectorId.kAD4:0,
00617                              DetectorId.kIWS:0,
00618                              DetectorId.kOWS:0,
00619                              DetectorId.kRPC:0
00620                              }
00621         return
00622 
00623     def initialize(self):
00624         status = DybPythonAlg.initialize(self)
00625         if status.isFailure(): return status
00626         self.info("initializing")
00627         # Synchronize the start time for sources
00628         for source in self.sources:
00629             source.sourceTime = self.startTime
00630         return SUCCESS
00631 
00632     def execute(self):
00633         status = self.fillNextTrigger()
00634         if status.isFailure(): return status
00635         return SUCCESS
00636         
00637     def finalize(self):
00638         self.info("finalizing")
00639         status = DybPythonAlg.finalize(self)
00640         return status
00641 
00642     ######################################################
00643     # Helper functions
00644 
00645     def fillNextTrigger(self):
00646         """ Place the next fake trigger into the event store """
00647         nextTrigger = self.getNextTrigger()
00648         if not nextTrigger: return FAILURE
00649         evt = self.evtSvc()
00650         headerPaths = ["/Event/CalibReadout/CalibReadoutHeader",
00651                        "/Event/Data/CalibStats",
00652                        "/Event/Rec/AdSimple",
00653                        "/Event/Rec/AdScaled"
00654                        ]
00655         headers = nextTrigger['headers']
00656         for headerPath in headerPaths:
00657             # Add item to event store
00658             if headerPath in headers:
00659                 evt[headerPath] = headers[headerPath]
00660         return SUCCESS
00661 
00662     def getNextTrigger(self):
00663         """ Get the data headers for the next trigger """
00664         status = self.updateSources()
00665         if status.isFailure(): return None
00666         source = self.findNextSource()
00667         if not source: return None
00668         nextTrigger = source.getNextTrigger()
00669         if not self.ableToTrigger(nextTrigger):
00670             # If trigger system unable to issue trigger, try again
00671             nextTrigger = self.getNextTrigger()
00672         # Set exec number
00673         for header in nextTrigger['headers'].values():
00674             header.setExecNumber(self.execNumber)
00675         # Set dt last trigger
00676         status = self.setDtLastTrigger(nextTrigger)
00677         if status.isFailure(): return None
00678         # Set trigger number for this detector
00679         status = self.setTriggerNumber(nextTrigger)
00680         if status.isFailure(): return None
00681         # Update time since this detector last triggered
00682         triggerTime = TimeStamp(nextTrigger['triggerTime'])
00683         detectorId = nextTrigger['detectorId']
00684         self.lastTriggerTime[detectorId] = triggerTime
00685         self.execNumber += 1
00686         self.triggerCount[detectorId] += 1
00687         return nextTrigger
00688 
00689     def ableToTrigger(self, nextTrigger):
00690         "Check if trigger system is able to issue trigger"
00691         detectorId = nextTrigger['detectorId']
00692         lastTriggerTime = self.lastTriggerTime[detectorId]
00693         nextTriggerTime = nextTrigger['triggerTime']
00694         if lastTriggerTime >= nextTriggerTime:
00695             # Catch case where next trigger time == last trigger time
00696             return False
00697         dtTrigger = TimeStamp(nextTriggerTime)
00698         dtTrigger.Subtract(lastTriggerTime)
00699         if dtTrigger.GetSeconds() < self.minRetriggerTime:
00700             return False
00701         return True
00702 
00703     def setDtLastTrigger(self, nextTrigger):
00704         "Set the CalibStats time since last trigger"
00705         detectorId = nextTrigger['detectorId']
00706         site = nextTrigger['site']
00707         triggerTime = nextTrigger['triggerTime']
00708         calibStats = nextTrigger['headers']['/Event/Data/CalibStats']
00709         for detectorId in self.lastTriggerTime:
00710             if (site != Site.kFar
00711                 and (DetectorId.kAD3 <= detectorId <= DetectorId.kAD4)):
00712                 continue
00713             if self.lastTriggerTime[detectorId] > TimeStamp.GetBOT():
00714                 dtLast = TimeStamp(triggerTime)
00715                 dtLast.Subtract(self.lastTriggerTime[detectorId])
00716                 calibStats.setFloat('dtLast%s_ms'% detectorNames[detectorId],
00717                                     dtLast.GetSeconds()*1e3)
00718         return SUCCESS
00719 
00720     def setTriggerNumber(self, nextTrigger):
00721         "Set the trigger number"
00722         detectorId = nextTrigger['detectorId']
00723         headers = nextTrigger['headers']
00724         triggerNumber = self.triggerCount[detectorId]
00725         calibHeader = headers['/Event/CalibReadout/CalibReadoutHeader']
00726         calibHeader.calibReadout().setTriggerNumber(triggerNumber)
00727         adSimple = headers['/Event/Rec/AdSimple']
00728         adSimple.recTrigger().setTriggerNumber(triggerNumber)
00729         adScaled = headers['/Event/Rec/AdScaled']
00730         adScaled.recTrigger().setTriggerNumber(triggerNumber)
00731         return SUCCESS
00732 
00733     def updateSources(self):
00734         """ Ask each source to update to next trigger time """
00735         nextSourceTime = self.nextSourceTime()
00736         upToDate = False
00737         while not upToDate:
00738             upToDate = True
00739             for source in self.sources:
00740                 if source.currentTime() <= nextSourceTime:
00741                     status = source.makeNextEvent()
00742                     if status.isFailure(): return status
00743                     nextSourceTime = self.nextSourceTime()
00744                     upToDate = False
00745                     break
00746         return SUCCESS
00747 
00748     def findNextSource(self):
00749         """ Determine which source will provide the next trigger """
00750         if len(self.sources)<1:
00751             self.error("No Sources have been defined")
00752             return None
00753         curSource = self.sources[0]
00754         for source in self.sources[1:]:
00755             if source.nextTriggerTime() < curSource.nextTriggerTime():
00756                 curSource = source
00757         return curSource
00758 
00759     def nextSourceTime(self):
00760         """ Find the time of the next source """
00761         nextSource = self.findNextSource()
00762         if not nextSource:
00763             return TimeStamp.GetEOT()
00764         return nextSource.nextTriggerTime()
00765         
00766 #####  Job Configuration for nuwa.py ########################################
00767 
00768 def configure( argv=[] ):
00769     """ Configure coincidence testing algorithm """
00770     return
00771 
00772 def run(app):
00773     '''
00774     Configure and add coincidence testing algorithm to job
00775     '''
00776     myAlg = TestCoincAlg("MyTestCoincAlg")
00777     # IBD in DayaBay AD1
00778     ibdSrc_DayaBayAD1 = IBDSource("ibdSrc_DayaBayAD1",700./(24*60*60.))
00779     ibdSrc_DayaBayAD1.site = Site.kDayaBay
00780     ibdSrc_DayaBayAD1.detectorId = DetectorId.kAD1
00781     myAlg.sources += [ibdSrc_DayaBayAD1]
00782     # IBD in DayaBay AD2
00783     ibdSrc_DayaBayAD2 = IBDSource("ibdSrc_DayaBayAD2",700./(24*60*60.))
00784     ibdSrc_DayaBayAD2.site = Site.kDayaBay
00785     ibdSrc_DayaBayAD2.detectorId = DetectorId.kAD2
00786     myAlg.sources += [ibdSrc_DayaBayAD2]
00787     # Singles in DayaBay AD1
00788     singlesSrc_DayaBayAD1 = ADSinglesSource("singlesSrc_DayaBayAD1",200.0) #Hz
00789     singlesSrc_DayaBayAD1.site = Site.kDayaBay
00790     singlesSrc_DayaBayAD1.detectorId = DetectorId.kAD1
00791     myAlg.sources += [singlesSrc_DayaBayAD1]
00792     # Singles in DayaBay AD2
00793     singlesSrc_DayaBayAD2 = ADSinglesSource("singlesSrc_DayaBayAD2",200.0) #Hz
00794     singlesSrc_DayaBayAD2.site = Site.kDayaBay
00795     singlesSrc_DayaBayAD2.detectorId = DetectorId.kAD2
00796     myAlg.sources += [singlesSrc_DayaBayAD2]
00797     # Muons in DayaBay site
00798     muonSrc_DayaBay = MuonSource("muonSrc_DayaBay",200.0) #Hz
00799     muonSrc_DayaBay.site = Site.kDayaBay
00800     myAlg.sources += [muonSrc_DayaBay]
00801     # Add algorithm
00802     app.addAlgorithm(myAlg)
00803     pass
00804 
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:10:46 for ADCoincTagging by doxygen 1.7.4