testCosmicT0Finder.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
 
 
 
# ---------------------------------------------------------------------------------
# Test of performance and validation of TOPCosmicT0Finder with simulation
# Input file can be prepared with top/examples/simGCRData.py
# Output ntuple to cosmicT0FinderNtuple.root
# ---------------------------------------------------------------------------------
 
import basf2 as b2
from ROOT import Belle2
from ROOT import TH1F
from ROOT import gROOT, AddressOf
import math
import ROOT
 
timeShift = 100  # ns (can be arbitrary chosen)
 
 
class ShiftDigits(b2.Module):
    ''' shift digits in time by timeShift'''
 
    def event(self):
        ''' event processing '''
 
        digits = Belle2.PyStoreArray('TOPDigits')
        for digit in digits:
            digit.setTime(digit.getTime() + timeShift)
 
 
gROOT.ProcessLine('struct TreeStruct {\
   int slot;        /* slot number */ \
   int nfot;        /* number of photons in the slot */ \
   int pdg;         /* PDG code from barHit */ \
   float x;         /* local x from barHit */ \
   float y;         /* local y from barHit */ \
   float z;         /* local z from barHit */ \
   float t;         /* time from barHit */ \
   float p;         /* momentum from barHit */ \
   float theta;     /* momentum theta [deg] from barHit */ \
   float phi;       /* momentum local phi [deg] from barHit */ \
   float dx;        /* extHit - barHit difference in local x */ \
   float dy;        /* extHit - barHit difference in local y */ \
   float dz;        /* extHit - barHit difference in local z */ \
   float dt;        /* extHit - barHit difference in time */ \
   float dp;        /* extHit - barHit difference in momentum */ \
   float dtheta;    /* extHit - barHit difference in theta [deg] */ \
   float dphi;      /* extHit - barHit difference in local phi [deg] */ \
   float t0;        /* t0 determined by TOPCosmicT0Finder */ \
   float t0err;     /* error on t0 */ \
   float timeShift; /* time shift applied to TOPDigits */ \
   float dt0;       /* difference: dt + t0 - timeShift */ \
};')
 
from ROOT import TreeStruct  # noqa
 
 
class Ntuple(b2.Module):
    ''' makes ntuple for MC studies of TOPCosmicT0Finder performance '''
 
    def initialize(self):
        ''' initialize '''
 
        
        self.filefile = ROOT.TFile('cosmicT0FinderNtuple.root', 'recreate')
        
        self.treetree = ROOT.TTree('tree', '')
        
        self.datadata = TreeStruct()
 
        for key in TreeStruct.__dict__.keys():
            if '__' not in key:
                formstring = '/F'
                if isinstance(self.datadata.__getattribute__(key), int):
                    formstring = '/I'
                self.treetree.Branch(key, AddressOf(self.datadata, key), key + formstring)
 
        
        self.resoreso = TH1F('reso', 'T0 residuals', 200, -1.0, 1.0)
        self.resoreso.SetXTitle('#Delta t0 [ns]')
        
        self.pullpull = TH1F('pull', 'T0 pulls', 200, -10.0, 10.0)
 
    def event(self):
        ''' event processing: fill ntuple and histograms '''
 
        mcParticles = Belle2.PyStoreArray('MCParticles')
        timeZeros = Belle2.PyStoreArray('TOPTimeZeros')
        for timeZero in timeZeros:
            extHit = timeZero.getRelated('ExtHits')
            if not extHit:
                b2.B2WARNING('no related extHits')
                continue
            barHit = None
            for barhit in mcParticles[0].getRelationsWith('TOPBarHits'):
                if barhit.getModuleID() == extHit.getCopyID():
                    barHit = barhit
            if not barHit:
                b2.B2WARNING('no corresponding TOPBarHit')
                continue
 
            moduleID = timeZero.getModuleID()
            if extHit.getCopyID() != moduleID:
                b2.B2ERROR('moduleID in extHit differs: ' + str(extHit.getCopyID()) +
                           ' ' + str(moduleID))
            if barHit.getModuleID() != moduleID:
                b2.B2ERROR('moduleID in barHit differs: ' + str(barHit.getModuleID()) +
                           ' ' + str(moduleID))
 
            self.datadata.slot = moduleID
            self.datadata.nfot = timeZero.getNumPhotons()
            self.datadata.pdg = barHit.getPDG()
            pos = barHit.getLocalPosition()
            self.datadata.x = pos.x()
            self.datadata.y = pos.y()
            self.datadata.z = pos.z()
            self.datadata.t = barHit.getTime()
            phiBar = math.pi / 2 - math.pi / 8 * (moduleID - 0.5)
            mom = barHit.getMomentum()
            mom.RotateZ(phiBar)
            self.datadata.p = mom.Mag()
            self.datadata.theta = math.degrees(mom.Theta())
            self.datadata.phi = math.degrees(mom.Phi())
            dr = extHit.getPosition() - barHit.getPosition()
            dr.RotateZ(phiBar)
            self.datadata.dx = dr.x()
            self.datadata.dy = dr.y()
            self.datadata.dz = dr.z()
            self.datadata.dt = extHit.getTOF() - self.datadata.t
            momExt = extHit.getMomentum()
            momExt.RotateZ(phiBar)
            self.datadata.dp = momExt.Mag() - self.datadata.p
            self.datadata.dtheta = math.degrees(momExt.Theta()) - self.datadata.theta
            self.datadata.dphi = math.degrees(momExt.Phi()) - self.datadata.phi
            self.datadata.t0 = timeZero.getTime()
            self.datadata.t0err = timeZero.getError()
            self.datadata.timeShift = timeShift
            self.datadata.dt0 = self.datadata.dt + self.datadata.t0 - self.datadata.timeShift
            self.filefile.cd()
            self.treetree.Fill()
            self.resoreso.Fill(self.datadata.dt0)
            self.pullpull.Fill(self.datadata.dt0/self.datadata.t0err)
            if timeZero.getChi2().GetEntries() > 0:
                timeZero.getChi2().Write()
                timeZero.getPDF().Write()
                timeZero.getHits().Write()
 
    def terminate(self):
        ''' terminate: write and close root file '''
 
        self.filefile.cd()
        self.filefile.Write()
        self.filefile.Close()
 
 
# Suppress messages and warnings during processing:
b2.set_log_level(b2.LogLevel.WARNING)
 
# Define a global tag (note: the one given bellow will become out-dated!)
b2.use_central_database('data_reprocessing_proc8')
 
# Create path
main = b2.create_path()
 
# input
roinput = b2.register_module('RootInput')
main.add_module(roinput)
 
# Initialize TOP geometry parameters (creation of Geant geometry is not needed)
main.add_module('TOPGeometryParInitializer')
 
# shift digits in time
main.add_module(ShiftDigits())
 
# t0 finder
finder = b2.register_module('TOPCosmicT0Finder')
finder.param('useIncomingTrack', False)  # or True
# finder.param('saveHistograms', True)
main.add_module(finder)
 
# ntuple
main.add_module(Ntuple())
 
# Print progress
main.add_module('Progress')
 
# Process events
b2.process(main)
 
# Print statistics
print(b2.statistics)