AmCTag.py

Go to the documentation of this file.

#!/usr/bin/env python
'''
AmC Tagger 
Usage:
    nuwa.py --no-history -A "0.1*s" -n -1 -m"Quickstart.EH1Tables" -m"Quickstart.Calibrate" -m"Quickstart.CalculateCalibStats"  -m"Quickstart.Reconstruct" -m"MuonTagging.MuonTagLoose" -m"CalibrationTagging.AmCTag" -o output.root input.root
Help:
    nuwa.py -m"CalibrationTagging.AmCTag -h"
    
tsang@caltech.edu, 2011/08/22
'''

from UserTagging.UserTaggingAlg import UserTaggingAlg
from GaudiPython import SUCCESS, FAILURE
from GaudiPython import gbl
import GaudiKernel.SystemOfUnits as units
import math

# ===========================================
class AmCTag(UserTaggingAlg):
    'Coincidence Tagger for AD events'

    def __init__(self, name):
        UserTaggingAlg.__init__(self, name)
        
        self.dTCut = 200*units.microsecond
        self.vertexCut = 1*units.meter
        self.ratioQCut = 0.25
        self.location = "/Event/Rec/AdSimple"
        self.muonVetoTime = 1000*units.microsecond
        self.adMuonEnergyCut = 20*units.MeV
        self.promptLowCut   = 0.5*units.MeV
        self.promptHighCut  = 2.0*units.MeV
        self.delayedLowCut  = 6.5*units.MeV
        self.delayedHighCut = 9.5*units.MeV
        self.timeWindowLow  = -1640 #ns
        self.timeWindowHigh = -1450 #ns
        self.tRmsCut = 30 
            
    # -------------------------------------------
    def initTagList(self):
        self.detId = [
            gbl.DetectorId.kAD1,
            gbl.DetectorId.kAD2,
            gbl.DetectorId.kAD3,
            gbl.DetectorId.kAD4]
        self.inputHeaders = {}
        self.lastTriggerTime = {}  # dict
        self.lastPromptTime = {}  # dict
        self.lastDelayedTime = {} # dict
        self.lastMuonTime = {}
        self.lastWasPrompt = {}
        self.stage = {}  # 0: looking for prompt
                         # 1: prompt found, looking for delayed
                         # 2: delayed found, ensuring clearance
                         # 3: pair found!
        self.freshStart = {}
        for id in self.detId:
            self.inputHeaders[id] = []
            self.freshStart[id] = True
            self.lastWasPrompt[id] = False
            self.stage[id] = 0

        self.addTag('AmC' , '/Event/Tag/Calib/AmC')
        return
        
    # -------------------------------------------
    def check(self, evt):
        print '='*80
        recHeader = evt['/Event/Rec/AdSimple']
        if not recHeader:
            self.warning('Failed to get current recon header')
            return FAILURE

        # Access the Calib Readout Header.
        # This is a container for calibrated data
        calibHdr = evt["/Event/CalibReadout/CalibReadoutHeader"]
        if calibHdr == None:
            self.error("Failed to get current calib readout header")
            return FAILURE

        # Access the Readout.  This is the calibrated data from one trigger.
        calibReadout = calibHdr.calibReadout()
        if calibReadout == None:
            self.error("Failed to get calibrated readout from header")
            return FAILURE

        recResult = recHeader.recTrigger()
        triggerTime = recResult.triggerTime()
        detectorId = recHeader.context().GetDetId()

        # Fresh start: assume first event is muon
        if not detectorId in self.lastMuonTime.keys(): #self.detId and self.freshStart[detectorId]:
            self.lastMuonTime[detectorId] = triggerTime
            self.lastTriggerTime[detectorId] = triggerTime
            self.freshStart[detectorId] = False
            print 'Fresh Start!',detectorId
            return SUCCESS

        # Treat all non-AD triggers as muons
        if not detectorId in self.detId:
#            for id in self.detId:
#              self.lastMuonTime[id] = triggerTime
#            print 'Treat all non-AD triggers as muons.'
            print 'Ignore all non-AD triggers.'
            return SUCCESS

        print triggerTime.AsString(), detectorId
        print 'Last trigger time:',self.lastTriggerTime[detectorId].AsString()
        print 'Last muon time:',self.lastMuonTime[detectorId].AsString()
        timeSinceMu = gbl.TimeStamp( recResult.triggerTime() )
        timeSinceMu.Subtract(self.lastMuonTime[detectorId])
        timeSinceTrig = gbl.TimeStamp( recResult.triggerTime() )
        timeSinceTrig.Subtract(self.lastTriggerTime[detectorId])

        reject = False
        # Get Recon info
        recES = recResult.energyStatus() 
        recE = recResult.energy()
        recPS = recResult.positionStatus()
        recX = recResult.position().x()
        recY = recResult.position().y()
        recZ = recResult.position().z()

        # Fix for wrong energy scale
        recE *= 184./160.

        # Ignore triggers with failed recon
        if not reject and recPS != 1:
          print 'Position recon failed. Trigger ignored'
          reject = True
        if not reject and recES != 1:
          print 'Energy recon failed. Trigger ignored'
          reject = True

        # Calculate charge sum and other quantities
        maxQ = None
        sumQ = 0 
        nChan = 0 
        meanT = 0 
        meanT2 = 0 
        for channel in calibReadout.channelReadout():
          q = 0 
          for hitIdx in range( channel.size() ):
            hitTime = channel.time( hitIdx )
            hitCharge = channel.charge( hitIdx )
            if hitTime > self.timeWindowLow and hitTime < self.timeWindowHigh:
              nChan += 1
              q = hitCharge
              meanT += hitTime
              meanT2 += hitTime * hitTime
              break
          if not maxQ or maxQ < q: maxQ = q 
          sumQ += q
        if sumQ != 0 and nChan > 0:
          ratioQ = maxQ/sumQ
          meanT /= nChan
          meanT2 /= nChan
          Trms = math.sqrt(meanT2-meanT*meanT)
        else:
          ratioQ = meanT = meanT2 = Trms = -1

#        # Get ratio Max(Q)/Sum(Q)
#        #ratioMaxQ < 0.25
#        maxQ = None
#        sumQ = 0
#        for channel in calibReadout.channelReadout():
#          q = 0
#          for hitIdx in range( channel.size() ):
#            hitTime = channel.time( hitIdx )
#            hitCharge = channel.charge( hitIdx )
#            if hitTime > -1640 and hitTime < -1450:
#              q = hitCharge
#              break
#          if not maxQ or maxQ < q: maxQ = q
#          sumQ += q
#        if sumQ != 0:
#          ratioQ = maxQ/sumQ
#        else:
#          ratioQ = 0

#        # Max/Sum cut
#        if not reject and (ratioQ > self.ratioQCut or sumQ == 0):
#          print 'Max(Q)/Sum(Q)','=', ratioQ,'>', self.ratioQCut, '. Trigger ignored.'
#          reject = True

        # Muon trigger? If so, start veto
        isMuon = False
        #isMuon = evt['/Event/Tag/Physics/MuonLoose']
        if recE > self.adMuonEnergyCut:
          isMuon = True
          
        if isMuon:
          print 'Loose Muon in',detectorId, '. Start veto:', self.muonVetoTime
          reject = True

        # Muon veto
        if not reject:
          if timeSinceMu.GetSeconds() * units.second < self.muonVetoTime:
            print 'Muon veto on. Trigger ignored', detectorId
            reject = True

        # Discard non-physics events
        if not reject and nChan == 0:
          print 'No hits selected in hit time window. Ignored.'
          reject = True
        if not reject and ratioQ > self.ratioQCut:
          print 'Max(Q)/Sum(Q)','=', ratioQ,'>', self.ratioQCut, '. Ignored.'
          reject = True
        if not reject and Trms > self.tRmsCut:
          print 'tRms =', Trms,'>',self.tRmsCut, '. Ignored.'
          reject = True


        # Energy cut for identifying prompt or delayed signal
        isPrompt = False
        isDelayed = False
        if not reject:
          if recE > self.promptLowCut and recE < self.promptHighCut:
            isPrompt = True
            print 'Prompt candidate'
          if recE > self.delayedLowCut and recE < self.delayedHighCut: 
            if self.stage[detectorId] == 1:
              promptRecResult = self.inputHeaders[detectorId][0].recTrigger() 
              pRecX = promptRecResult.position().x()
              pRecY = promptRecResult.position().y()
              pRecZ = promptRecResult.position().z()
              deltaX = (pRecX - recX)
              deltaY = (pRecY - recY)
              deltaZ = (pRecZ - recZ)
              deltaVertex2 = deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ
              if deltaVertex2 < self.vertexCut*self.vertexCut:
                isDelayed = True
              else:
                print 'rejected by vertex cut'
            else:
              isDelayed = True
            if isDelayed: print 'Delayed candidate'

        # Timing coincidence 
        print 'Stage', self.stage
        if reject:
          if self.stage[detectorId] == 1:
            self.stage[detectorId] = 0
            self.inputHeaders[detectorId] = []
            print '1->0'
        else:
#          print 'Current trigger time:',triggerTime.AsString()
#          print 'Last trigger time:',self.lastTriggerTime[detectorId].AsString()
#          print 'Last muon time:',self.lastMuonTime[detectorId].AsString()
          print 'dT =', timeSinceTrig.GetSeconds(), 'second'
          # Stage 1->2/0
          if self.stage[detectorId] == 1:
            if isDelayed and timeSinceTrig.GetSeconds() * units.second < self.dTCut:
              self.inputHeaders[detectorId].append(recHeader)
              self.stage[detectorId] = 2
              print '1->2'
            else:
              self.stage[detectorId] = 0
              self.inputHeaders[detectorId] = []
              print '1->0'
          # Stage 2->3/0 
          elif self.stage[detectorId] == 2:
            if timeSinceTrig.GetSeconds() * units.second > self.dTCut:
              self.stage[detectorId] = 3
              print '2->3'
              # save tag
              myTag = self.getTag('AmC')
              print 'Length of self.inputheaders', len(self.inputHeaders[detectorId])
              myTag.setInputHeaders(self.inputHeaders[detectorId])
              myTag.tagIt()
              self.inputHeaders[detectorId] = []
              self.stage[detectorId] = 0
              print 'Tag saved!'
              print '3->0'
            else:
              self.stage[detectorId] = 0
              self.inputHeaders[detectorId] = []
              print '2->0'
          # Stage 0->1
          if isPrompt and self.stage[detectorId] == 0 and timeSinceTrig.GetSeconds() * units.second > self.dTCut:
            self.inputHeaders[detectorId].append(recHeader)
            self.stage[detectorId] = 1
            print '0->1'


        if isMuon:
            self.lastMuonTime[detectorId] = triggerTime
        self.lastTriggerTime[detectorId] = recResult.triggerTime()
        return SUCCESS
        
#####  Job Configuration for nuwa.py ########################################
options = None

def configure(argv=[]):
    """Configuration with command line arguments"""   
    global options 
    from optparse import OptionParser
    
    parser = OptionParser()
    parser.add_option("-c", "--coincidence-dt", type="string",
        default="300*microsecond",
        help="AD coincidence time window")
    parser.add_option("-l", "--location", type="string",
        default="/Event/Rec/AdSimple",
        help="Target for AD coincidence tagging")
    (options, args) = parser.parse_args(args=argv)

def run(app):
    from DybPython.Tools import unitify
    
    myAlg = AmCTag("AmCTag")
    myAlg.dTCut = unitify(options.coincidence_dt)
    myAlg.location = options.location
    app.addAlgorithm(myAlg)