# ============================================================================= __author__ = "Daan van Eijk" __date__ = "2010-07-29" __version__ = "CVS Tag \$Name: v14r3 \$, version \$Revision: 1.1 \$ " # ============================================================================= ## import everything from bender #from Bender.Awesome import * from Bender.Main import * from LoKiCore.math import * from ROOT import Double, GaudiAlg from PartProp.Nodes import * from LoKiPhys.decorators import * #from GaudiPython.GaudiAlgs import HistoAlgo import math # ============================================================================= class HitEfficiency( Algo ): def uniqueLayer( self,channelID ): return (channelID.station()-1) * 4 + channelID.layer() def uniqueQuarter( self,channelID ): return self.uniqueLayer(channelID) * 4 + channelID.quarter() def uniqueModule( self,channelID ): return self.uniqueQuarter(channelID)*9 + channelID.module() - 1 def uniqueOtis( self,channelID ): return self.uniqueModule(channelID)*4 + (channelID.straw()-1)/32 """ Algorithm to fill NTuple for hit efficiency """ ## initialize the algorithm def initialize ( self ) : """ Initialize the algorithm """ ## initialize the base sc = Algo.initialize ( self ) ## initialize the base if sc.isFailure() : return sc #Get the OT self._detOT = self.getDet('/dd/Structure/LHCb/AfterMagnetRegion/T/OT') #Create the linear extrapolator tool self._linearextrapolator = self.tool(cpp.ITrackExtrapolator,'TrackLinearExtrapolator') #Create the interpolator tool self._interpolator = self.tool(cpp.ITrackInterpolator,'TrackInterpolator') #Create the decoder tool self._decoder = self.tool(cpp.IOTRawBankDecoder,'OTRawBankDecoder') #Define thresholds, book histograms self._mydict = {} beginstep = 1 endstep = 120 nthresholdperlayer = 10 firstthreshold = 1 layer = -1 #Will increase every time the counter modulo number of steps is equal to 1 for i in range(beginstep,endstep+1): if i%nthresholdperlayer == firstthreshold: layer += 1 #print 'i = ', i #print 'layer = ', layer self._mydict[i]= layer string = 'threshold_'+str(i) self.bookProfile2D(GaudiAlg.ID(string),'Histograms',-3157,3157,74,-2800.,2800.,100) print 'This is the dictionary of the steps and layers: ', self._mydict self._minOTHitsPerTrack = 20 self._maxChi2PerDoF = 2 self._maxDistError = 0.2 self._numUniqueStation = 3 self._numUniqueLayer =12 self._numUniqueModule = 432 self._numUniqueOtis = 432*4 self._maxChan = 128 return SUCCESS ## Main method def analyse ( self ) : """ Main method """ self._tuple = self.nTuple( 'Efficiency') self._odin = self.get( 'DAQ/ODIN' ) self._tracks = self.get( '/Event/Rec/Track/Best' ) runNum = self._odin.runNumber() evtNum = self._odin.eventNumber() #clear the hitmap self._moduleHitMap = [] for i in range(0,self._numUniqueModule): listOfchannels = [] for j in range(0,self._maxChan): listOfchannels.append(0.0) self._moduleHitMap.append(listOfchannels) #decode all OT hits ottimes = LHCb.OTLiteTimeContainer() sc = self._decoder.decode( ottimes ) #fill all the hits in the hit map if( sc ): for ihit in ottimes: self._moduleHitMap[self.uniqueModule(ihit.channel())][ihit.channel().straw()-1] = ihit.channel() if self._tracks: for track in self._tracks: if (( track.fitStatus() == track.Fitted ) and ( track.nDoF()>1 ) and (track.chi2PerDoF() < self._maxChi2PerDoF) ): nodes = track.nodes() _nTrackOTHits = 0 for node in nodes: if ( ( node.type() == node.HitOnTrack ) and ( node.measurement().type() == node.measurement().OT )): _nTrackOTHits += 1 #if _nTrackOTHits > 24: #print 'Found ', _nTrackOTHits, ' HitOnTrack nodes which are OT measurements on a track!' if _nTrackOTHits >= self._minOTHitsPerTrack: self.fillEfficiency(track,_nTrackOTHits) return SUCCESS def fillEfficiency(self, track,ntrackothits): modulesOnTrack = {} nodes = track.nodes() for node in nodes: if ( ( node.type() == node.HitOnTrack ) and ( node.measurement().type() == node.measurement().OT )): meas = node.measurement() module = meas.module() mono = ( meas.channel().straw()-1 ) / ( module.nChannels()/2 ) #FIXME: Find the node with the smallest residual if you have 2 nodes in one monolayer (Wouter told me this is a small effect...) if module.uniqueModule() in modulesOnTrack: if mono == 0: modulesOnTrack[module.uniqueModule()][0] = module modulesOnTrack[module.uniqueModule()][1][0] = node elif mono == 1: modulesOnTrack[module.uniqueModule()][0] = module modulesOnTrack[module.uniqueModule()][1][1] = node else: modulesOnTrack[module.uniqueModule()] = [0,[0,0]] if mono == 0: modulesOnTrack[module.uniqueModule()][0] = module modulesOnTrack[module.uniqueModule()][1][0] = node elif mono == 1: modulesOnTrack[module.uniqueModule()][0] = module modulesOnTrack[module.uniqueModule()][1][1] = node #FOR THE FIXME: ## // only take the closest node in a mono-layer. we can anyway only unbias 1. ## int mono = (meas->channel().straw()-1) / (module->nChannels()/2) ; ## const LHCb::Node*& nodeInModule = mono == 0 ? ## modulesOnTrack[ module ].first : modulesOnTrack[ module ].second ; ## if( nodeInModule == 0 || ## std::abs(nodeInModule->residual()) > ## std::abs((*inode)->residual()) ) nodeInModule = *inode ; ## //modulesOnTrack[ module ].push_back( *inode ) ; ## } #Process the efficiency for these modules. Make sure to unbias the state when necessary. keep track of the layers that we have seen, such that we can skip them below... layersOnTrack = [] for mod in modulesOnTrack: layersOnTrack.append(self._detOT.findLayer(modulesOnTrack[mod][0].elementID())) nodeMonoA = modulesOnTrack[mod][1][0] nodeMonoB = modulesOnTrack[mod][1][1] self.fillEfficiency2(track, modulesOnTrack[mod][0], nodeMonoA, nodeMonoB,ntrackothits) #Store the z-position of the layers and find the missing planes state = LHCb.State() modulesNotOnTrack = [] for ilay in self._detOT.layers(): zpos = ilay.geometry().toGlobal(Gaudi.XYZPoint(0.0,0.0,0.0)).z() if ilay not in layersOnTrack: self._interpolator.interpolate(track,zpos,state).ignore() if (math.sqrt(state.covariance()(0,0))< self._maxDistError): tolerance =0.1 #to get a more precies answer, to propagate the hit to the xy plane (modules can be tilted out of this plane) self._linearextrapolator.propagate(state,ilay.plane(),tolerance).ignore() module = self._detOT.findModule(state.position()) if module: modulesNotOnTrack.append(module) self.fillEfficiency4(track, module, state,ntrackothits) def fillEfficiency2( self, track, module, nodeMonoA, nodeMonoB,ntrackothits): states = [0,0] strawpos = [Double(),Double()] yfrac = [Double(),Double()] if nodeMonoA: states[0] = nodeMonoA.unbiasedState() else: states[0] = nodeMonoB.state() if nodeMonoB: states[1] = nodeMonoB.unbiasedState() else: states[1] = nodeMonoA.state() for imono in range(0,2): module.monoLayerIntersection(imono, states[imono].position(), states[imono].tx(), states[imono].ty(), strawpos[imono], yfrac[imono] ) refstate = states[1] #WHY? Because it's only used to calculate the slopes in fillefficiency3, and this change for monolayer A and B very little! self.fillEfficiency3(track, module, strawpos, yfrac, refstate,ntrackothits) def fillEfficiency3(self, track, module, strawpos, yfrac, refstate,ntrackothits): ##make a temporary structure that tells which hits are in this module. modid = module.elementID() pitch = module.xPitch() uniquemodule = self.uniqueModule(modid) nstraws = module.nChannels() #compute the direction in the local frame. localslopes = module.geometry().toLocal(refstate.slopes()) localTx = localslopes.x()/localslopes.z() cosalpha = 1/math.sqrt( 1+localTx*localTx) #find the hits in a +- 3 straw window for imono in range(0,2): monooffset = imono * nstraws/2 -1 minstraw = max(int(strawpos[imono])-2,1) maxstraw = min(int(strawpos[imono]) +3, int(nstraws/2)) for istraw in range(minstraw,maxstraw+1): channel = self._moduleHitMap[uniquemodule][istraw+monooffset] foundhit = False if channel: foundhit = True #Efficiency unbiasing procedure if (foundhit): #print "foundhit was true, so in if-statement" #print "_nTrackOTHits = ", ntrackothits #print "_minOTHitsPerTrack", self._minOTHitsPerTrack if (ntrackothits == self._minOTHitsPerTrack): #print "_nTrackOTHits == _minOTHitsPerTrack, continue !" continue #print "foundhit is ", foundhit, ". Filling...." monocoord = strawpos[imono]-istraw #print "monocoord = ", monocoord dstraw = monocoord*pitch*cosalpha pvcontainer = self.get('Rec/Vertex/Primary') step = self._odin.calibrationStep() layer = module.elementID().uniqueLayer()-4 #print 'step =', step #print 'layer =', layer if abs(dstraw) < 1.25 and pvcontainer.size() == 1 and self._tracks.size() > 0 and self._tracks.size() < 1000 and track.type() == 3: if step in self._mydict: if layer == self._mydict[step]: self.profile2D(refstate.x(), refstate.y(),foundhit,GaudiAlg.ID('threshold_'+str(step)),'Histograms',-3157,3157,-2800,2800,74,100,1) self.setFilterPassed ( True ) def fillEfficiency4(self, track, module, state,ntrackothits): strawpos = [Double(),Double()] yfrac = [Double(),Double()] for imono in range(0,2): module.monoLayerIntersection(imono, state.position(), state.tx(), state.ty(), strawpos[imono], yfrac[imono] ) self.fillEfficiency3(track, module, strawpos, yfrac, state,ntrackothits) # ============================================================================= ## configure the job # @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl # @date 2010-05-31 def configure ( datafiles , catalogs = [] ) : """ Configure the job """ ## ## Static Configuration: ## ## DaVinci Configurable from Configurables import Brunel brunel = Brunel ( DataType = '2011', DDDBtag = 'head-20110914', CondDBtag = 'head-20110914', OutputType = 'None' VetoHltErrorEvents = False ) from Configurables import RecMoniConf, RecSysConf, GaudiSequencer #RecSysConf().RecoSequence = ["Decoding","VELO","TT","IT","OT","Tr","Vertex","CALO","MUON"] RecSysConf().RecoSequence = ["Decoding","VELO","TT","IT","OT","Tr","Vertex"] RecMoniConf().setProp("MoniSequence",["OT"]) GaudiSequencer('LumiSeq').Enable = False GaudiSequencer('CaloBanksHandler').Enable = False #importOptions('\$APPCONFIGOPTS/DisableLFC.py') #importOptions('\$APPCONFIGOPTS/UseOracle.py') from Configurables import ( CondDBAccessSvc, CondDB ) AlignmentCondition = CondDBAccessSvc("AlignmentCondition") AlignmentCondition.ConnectionString = "sqlite_file:2011-05.v2-scan.db/LHCBCOND" CondDB().addLayer(AlignmentCondition) #CondDB().IgnoreHeartBeat = True from Configurables import HistogramPersistencySvc HistogramPersistencySvc ( OutputFile = 'histos.root' ) from Configurables import NTupleSvc NTupleSvc().Output += [ "ACCEPT DATAFILE='tuples.root' TYPE='ROOT' OPT='NEW'"] NTupleSvc().OutputLevel = 1 ## define input data setData ( datafiles ) ## ## Dynamic Configuration: Jump into the wonderful world of GaudiPython ## ## get the actual application manager (create if needed) gaudi = appMgr() ## create local algorithm: hitefficiency = HitEfficiency ('HitEfficiency', ## N-tuple LUN NTupleLUN = "ACCEPT" ) gaudi.addAlgorithm( hitefficiency ) return SUCCESS # ============================================================================= ## job steering # @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl # @date 2010-05-31 if __name__ == '__main__' : ## make printout of the own documentations print '*'*120 print __doc__ print ' Author : %s ' % __author__ print ' Version : %s ' % __version__ print ' Date : %s ' % __date__ print ' dir(%s) : %s ' % ( __name__ , dir() ) print '*'*120 ## configure the job: files = [ # "DATAFILE='castor:/castor/cern.ch/grid/lhcb/data/2011/RAW/FULL/LHCb/TEST/103331/103331_0000000002.raw' SVC='LHCb::MDFSelector'" "DATAFILE='/project/bfys/dveijk/ThresholdScans/TestResubmit/103332_0000000100.raw' SVC='LHCb::MDFSelector'" ] configure( files ) ## run the job run( 100 ) # ============================================================================= # The END # =============================================================================