readCo60data.py

Go to the documentation of this file.

#!/usr/bin/env python
#
# Usage:
#  nuwa.py -n -1 --dbconf=tmp_offline_db:offline_db -m "Quickstart.Calibrate" -m "Quickstart.CalculateCalibStats" -m "getCo60ES -d 2 -r 39329" <rawdatafile>
#
# jpochoa, June 2013

# load stuff
import os,math
from DybPython.DybPythonAlg import DybPythonAlg
from GaudiPython import SUCCESS, FAILURE
from GaudiPython import gbl

from ROOT import TF1
from ROOT import TObject

Detector = gbl.DayaBay.Detector
ServiceMode = gbl.ServiceMode
DetectorId = gbl.DetectorId
FeeChannelId = gbl.DayaBay.FeeChannelId
AdPmtSensor = gbl.DayaBay.AdPmtSensor
Context = gbl.Context 
TH1F = gbl.TH1F
TParameter = gbl.TParameter
TH2F = gbl.TH2F
Trigger = gbl.DayaBay.Trigger

# Algorithm
class getCo60ESAlg(DybPythonAlg):
    "Template Python Algorithm"
    def __init__(self,name):
        DybPythonAlg.__init__(self,name)
        self.det = None
        self.runno = None
        return

    def initialize(self):
        status = DybPythonAlg.initialize(self)
        if status.isFailure(): return status
        self.info("initializing")
        
        # global variables
        self.Nevts=0;
        self.prevTimeSec=0;
        self.prevTimeNanoSec=0;
        self.firstTimeSec=0;
        self.firstTimeNanoSec=0;
        self.lastTimeSec=0;
        self.lastTimeNanoSec=0;

        #  Statistics Service
        self.statsSvc = self.svc('IStatisticsSvc','StatisticsSvc')
        if self.statsSvc == None:
            self.error("Failed to initialize statistics service.")
            return FAILURE
        
        # energy histogram
        self.stats["/file0/Co60energy/h_e"] = TH1F("h_e","",100,0,740);
        self.stats["/file0/Co60energy/h_e"].GetXaxis().SetTitle("Q_{corr}^{nominal} (PEs)");

        self.stats["/file0/Co60energy/h_n"] = TH1F("h_n","",193,-0.5,192.5);
        self.stats["/file0/Co60energy/h_n"].GetXaxis().SetTitle("N_{hit}");
        
        return SUCCESS

    def execute(self):
        #self.info("executing")
        evt = self.evtSvc()

        # 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

        #  Remove events that are too close to other events in the same hall
        triggerTime = calibReadout.triggerTime()
        triggerTimeSec=triggerTime.GetSec()
        triggerTimeNanoSec=triggerTime.GetNanoSec()

        deltat=(triggerTimeSec-self.prevTimeSec)+1e-9*(triggerTimeNanoSec-self.prevTimeNanoSec)
        
        self.prevTimeSec=triggerTimeSec
        self.prevTimeNanoSec=triggerTimeNanoSec

        if deltat<200e-6:
            return SUCCESS;

        # Avoid irrelevant triggers (from other detectors and/or non-Physics)
        triggerType = calibReadout.triggerType();
        if triggerType & Trigger.kCalib != Trigger.kNone:
            return SUCCESS;
        if triggerType & Trigger.kRandom != Trigger.kNone:
            return SUCCESS;
        detId = Detector(calibHdr.context().GetSite(),calibHdr.context().GetDetId());
        if detId.isAD()!=1:
            return SUCCESS;
        if calibHdr.context().GetDetId()!=self.det:
            return SUCCESS;

        # Increment counter and record times
        self.Nevts+=1
        self.lastTimeSec=triggerTimeSec
        self.lastTimeNanoSec=triggerTimeNanoSec;
        if self.Nevts==1:
            self.firstTimeSec=triggerTimeSec;
            self.firstTimeNanoSec=triggerTimeNanoSec;
                    
        # Access the Calibrated Statistics Data Header.
        calibStats = evt["/Event/Data/CalibStats"]
        if calibStats == None:
            self.debug("No calibrated statistics!")
            return FAILURE
        
        # Access variables from CalibStats
        Qnominal = calibStats.get('NominalCharge').value()                
        nPMT = calibStats.get('nActivePMTs').value()
        nHit = calibStats.get('nHit').value()
        
        # fill histogram
        corrE=Qnominal*192/nPMT
        corrN=nHit*192/nPMT
        self.stats["/file0/Co60energy/h_e"].Fill(corrE);
        self.stats["/file0/Co60energy/h_n"].Fill(corrN);
        
        return SUCCESS

    def finalize(self):
        self.info("finalizing")

        # save first and last times (for normalization)
        print "trying to save estas ondas"
        self.stats["/file0/Co60energy/firstTimeSec"] = TParameter(int)("firstTimeSec",self.firstTimeSec)
        self.stats["/file0/Co60energy/firstTimeNanoSec"] = TParameter(int)("firstTimeNanoSec",self.firstTimeNanoSec)
        self.stats["/file0/Co60energy/lastTimeSec"] = TParameter(int)("lastTimeSec",self.lastTimeSec)
        self.stats["/file0/Co60energy/lastTimeNanoSec"] = TParameter(int)("lastTimeNanoSec",self.lastTimeNanoSec)
        
        # finalize--
        status = DybPythonAlg.finalize(self)
        return status


#####  Job Configuration for nuwa.py ########################################

det = None
runno = None

def configure( argv=[] ):
    """ Template job module """
    # Process module command line arguments here
    from optparse import OptionParser
    parser = OptionParser() 

    # Example module option
    parser.add_option("-d","--det",
                      default=1,
                      type="int",
                      help="DetectorId of AD with Co60 [default: %default]")

    parser.add_option("-r","--runno",
                      default=0,
                      type="int",
                      help="Run number [default: %default]")
    
    (options,args) = parser.parse_args(args=argv)

    # Must save option so it can be use in 'run()' function
    global det,runno
    det = options.det
    runno = options.runno
    print " Set DetectorId to ", det
    print " Set run to ", runno
    
    return

def run(app):
    '''
    Configure and add the algorithm to job
    '''
    # Add Python algorithm to job here
    myAlg = getCo60ESAlg("MyCo60ESAlg")
    # Retrieve module option, and set in algorithm
    global det
    global runno
    myAlg.det = det
    myAlg.runno = runno
    app.ExtSvc += ["StatisticsSvc","DybPmtCalibSvc"]
    # Add algorithm to job
    app.addAlgorithm(myAlg)
    pass