Tags:
create new tag
, view all tags

Script to run LHCb application (DaVinci, but can be used for others)

from ROOT import TFile,TH1F,TH2F,TCanvas,TBrowser,gStyle,TText,TMath,TF1,gROOT,MakeNullPointer,gSystem

# get the basic configuration from here                                                                                                     
#from LHCbConfig import *                                                                                                                   
from GaudiConf.Configuration import *
from Configurables import DstConf
from GaudiPython.Bindings import gbl, AppMgr, Helper
from AnalysisPython import Dir, Functors

import GaudiPython

gbl = GaudiPython.gbl
#ParticleID     = gbl.LHCb.ParticleID                                                                                                       
#MCTrackInfo    = gbl.MCTrackInfo                                                                                                           
#XYZVector      = gbl.ROOT.Math.XYZVector                                                                                                   
Particle       = gbl.LHCb.Particle
#LorentzVector  = gbl.Gaudi.LorentzVector                                                                                                   
#MCParticle     = gbl.LHCb.MCParticle                                                                                                       
Track          = gbl.LHCb.Track
RecVertex      = gbl.LHCb.RecVertex

appConf = ApplicationMgr( OutputLevel = INFO, AppName = 'myApp' )


lhcbApp = LHCbApp(
#                  DDDBtag = 'default',                                                                                                     
#                  CondDBtag = 'default',                                                                                                   
                  DataType = '2010',
                  Simulation = False)

dstConf = DstConf(EnableUnpack=True)

myPVLocation = "/Event/Rec/Vertex/myOffline"
myPV3DLocation = "/Event/Rec/Vertex/myP3D"
PVLocation =  "/Event/Rec/Vertex/Offline"
PVDummyLocation =  "/Event/Rec/Vertex/DummyPV"
trBestLocation = "/Event/Rec/Track/Best"
#                                                                                                                                           
# PV reconstruction                                                                                                                         
from Configurables import PatPVOffline, PVOfflineTool, LSAdaptPVFitter, PatPV3D, LSAdaptPV3DFitter, PVSeedTool

# dummy PV run to read tracks                                                                                                               
pvdummy = PatPVOffline("pvdummy")
pvdummy.OutputVertices = PVDummyLocation

# my PV offline                                                                                                                             
pvoffmy = PatPVOffline("pvoffmy")
pvoffmy.OutputVertices = myPVLocation
pvoffmy.addTool(PVOfflineTool,'PVOfflineTool')
pvoffmy.PVOfflineTool.addTool(LSAdaptPVFitter,'LSAdaptPVFitter')
pvoffmy.PVOfflineTool.PVSeedingName = "PVSeed3DTool"
pvoffmy.PVOfflineTool.LSAdaptPVFitter.OutputLevel = 2
pvoffmy.PVOfflineTool.LSAdaptPVFitter.trackMaxChi2 = 9.
pvoffmy.PVOfflineTool.LSAdaptPVFitter.trackMaxChi2Remove = 16.
pvoffmy.OutputLevel = 2
                                                                                  
# PV 3D                                                                                                                                     
pv3d = PatPV3D("pv3d")
pv3d.OutputVerticesName = myPV3DLocation
#pv3d.OutputLevel = 2                                                                                                                       
pv3d.addTool(PVOfflineTool,'PVOfflineTool')
pv3d.PVOfflineTool.addTool(LSAdaptPV3DFitter,'LSAdaptPV3DFitter')
pv3d.PVOfflineTool.PVSeedingName = "PVSeed3DTool"
pv3d.PVOfflineTool.LSAdaptPV3DFitter.trackMaxChi2 = 9.
pv3d.PVOfflineTool.LSAdaptPV3DFitter.trackMaxChi2Remove = 25.

appConf.TopAlg += [ pvdummy,pvoffmy, pv3d]                                                                                                 

appMgr = GaudiPython.AppMgr()
appMgr.config( files = ['$GAUDIPOOLDBROOT/options/GaudiPoolDbRoot.opts'])
#appMgr.initialize()                                                                                                                        

appMgr.ExtSvc += ['DataOnDemandSvc']

sel  = appMgr.evtsel()
sel.PrintFreq = 100
nextEvent = Functors.NextEvent(appMgr)
evt  = appMgr.evtsvc()
his  = appMgr.histsvc()
det  = appMgr.detsvc()
toolSvc = appMgr.toolsvc()
#partsvc = appMgr.ppSvc()                                                                                                                   
msg  = appMgr.service('MessageSvc', 'IMessageSvc')
dps  = appMgr.service('EventDataSvc', 'IDataProviderSvc')

#from Configurables import MyUtilsTool                                                                                                      
#myutils = toolSvc.create('MyUtilsTool', interface = "IGenericTool")                                                                        

MaxEvents = 2000
files = []
#file1  = 'PFN:/castor/cern.ch/grid/lhcb/data/2010/MINIBIAS.DST/00007022/0000/00007022_00000806_1.minibias.dst'                             
#file1  = 'PFN:/castor/cern.ch/grid/lhcb/data/2010/MINIBIAS.DST/00007022/0000/00007022_00000001_1.minibias.dst'                             
# probelmatic event wrt PV reco                                                                                                             
file1  = 'PFN:/afs/cern.ch/user/t/tblake/public/Stripping/000000.Bhadron.dst'
for i in range(10):
  varname = 'file'+str(i)
  if varname in vars():
    files.append(vars()[varname])

sel.open(files)

# histogram booking                                                                                                                         
h_1    = TH1F('h_1','Best size',100,0.0,500.)

h_nPV  = TH1F('h_nPV','n PVs', 11, -0.5, 10.5)
h_nPVmy  = TH1F('h_nPVmy','n PVs ,my', 11, -0.5, 10.5)

h_multiPV = TH1F('h_multiPV','pv multiplicity',51,-0.5,50.5)
h_multiPVmy = TH1F('h_multiPVmy','pv multiplicity my',51,-0.5,50.5)

h_multiFalsePV = TH1F('h_multiFalsePV','false pv multiplicity',21,-0.5,20.5)
h_multiFalsePVmy = TH1F('h_multiFalsePVmy','false pv multiplicity my',21,-0.5,20.5)

h_xPV    = TH1F('h_xPV',  'x PVs', 100,  0.0, 1.0)
h_yPV    = TH1F('h_yPV',  'y PVs', 100, -0.5, 0.5)
h_zPV    = TH1F('h_zPV',  'z PVs', 100, -200.0, 200.0)

h_xPVmy  = TH1F('h_xPVmy','x PVs', 100,  0.0, 1.0)
h_yPVmy  = TH1F('h_yPVmy','y PVs', 100, -0.5, 0.5)
h_zPVmy  = TH1F('h_zPVmy','z PVs', 100, -200.0, 200.0)

h_dzPV  = TH1F('h_dzPV','dz PVs', 100, -100., 100.)
h_dzPVmy  = TH1F('h_dzPVmy','dz PVs my', 100, -100., 100.)
#h_dzPVmy  = TH1F('h_dzPVmy','dz PVs my', 50, -9.5, 10.5)                                                                                   

####################################                                                                                                        
# loop over pvs and fill dz histo                                                                                                           
####################################                                                                                                        
def dz_pvs(PVs, hdzPV, hFalseMult):
  for i,ipv in enumerate(PVs):
    for k,kpv in enumerate(PVs):
      if k != i:
        dz = kpv.position().z()-ipv.position().z()
        hdzPV.Fill(dz)
        if abs(dz)<1.0:
           print 'i,k,dz ',i,k,dz,ipv.tracks().size(),kpv.tracks().size(),'\n'
           hFalseMult.Fill(kpv.tracks().size())


# excution                                                                                                                                  

pvLocation = "/Event/Rec/Vertex/Primary"
bestTrackLocation = "/Event/Rec/Track/Best"

nEvents=0
nRecEvents = 0
nPrimaryVertices = 0
while ( nextEvent() ) :
   nEvents+=1
#   print 'nev',nEvents                                                                                                                     
   if nEvents>MaxEvents: break
   # get the primary vertices                                                                                                               
#   evt.dump()                                                                                                                              
   best_trs = evt[bestTrackLocation]
   myPVs = evt[myPVLocation]
   if myPVs :
     dz_pvs(myPVs, h_dzPVmy, h_multiFalsePVmy)
     h_nPVmy.Fill(myPVs.size() )
     for PV in myPVs :
         if PV :
            h_xPVmy.Fill(PV.position().x() )
            h_yPVmy.Fill(PV.position().y() )
            h_zPVmy.Fill(PV.position().z() )
            pvTracks = PV.tracks()
            if pvTracks != None :
               h_multiPVmy.Fill(pvTracks.size() )
#            sinfo = "chi2"                                                                                                                 
#            for best_tr in best_trs:                                                                                                       
#              print myutils.trackToVertexInfo(best_tr,PV,sinfo)                                                                            

   PVs = evt[pvLocation]
   if PVs :
# dz of PVs histo                                                                                                                           
#     print 'len pvs ',len(PVs)                                                                                                             
     dz_pvs(PVs, h_dzPV, h_multiFalsePV)
##  z,mult PV                                                                                                                               
     nRecEvents += 1
     nPrimaryVertices += PVs.size()
     h_nPV.Fill(PVs.size() )
     for PV in PVs :
         if PV :
            h_xPV.Fill(PV.position().x() )
            h_yPV.Fill(PV.position().y() )
            h_zPV.Fill(PV.position().z() )
            pvTracks = PV.tracks()
            if pvTracks != None :
               h_multiPV.Fill(pvTracks.size() )


#  success = h_1.Fill(n)                                                                                                                    
#  evt.dump()                                                                                                                               
print "Analysed ", nEvents, " events"

f=TFile('myHistos.root','recreate')
for h in gROOT.GetList() :
  h.Write()
f.Close()
#gStyle.SetOptStat(1100011)                                                                                                                 
gStyle.SetOptFit(111)
tks = TCanvas('tks','my plots',1200,600)
tks.Divide(4,2)
tks.cd(1)
h_nPV.Draw()
tks.cd(2)
h_dzPV.Draw()
tks.cd(3)
h_multiFalsePV.Draw()
tks.cd(4)
h_multiPV.Draw()

tks.cd(5)
h_nPVmy.Draw()
tks.cd(6)
h_dzPVmy.Draw()
tks.cd(7)
h_multiFalsePVmy.Draw()
tks.cd(8)
h_multiPVmy.Draw()

tkres = TCanvas('tkres','res plots',900,600)
tkres.Divide(3,2)
tkres.cd(1)
h_xPV.Draw()
tkres.cd(2)
h_yPV.Draw()
tkres.cd(3)
h_zPV.Draw()
tkres.cd(4)
h_xPVmy.Draw()
tkres.cd(5)
h_yPVmy.Draw()
tkres.cd(6)
h_zPVmy.Draw()

tkpre1 = TCanvas('tkpre1','PV1 plots',300,600)
tkpre1.Divide(1,2)
tkpre1.cd(1)
h_nPV.Draw()
tkpre1.cd(2)
h_dzPV.Draw()

tkpre2 = TCanvas('tkpre2','PV2 plots',300,600)
tkpre2.Divide(1,2)
tkpre2.cd(1)
h_nPVmy.Draw()
tkpre2.cd(2)
h_dzPVmy.Draw()

a = raw_input('continue? : ')



-- MariuszWitek - 28 Apr 2011

Topic revision: r1 - 2011-04-28 - MariuszWitek
 
This site is powered by the TWiki collaboration platform Powered by Perl This site is powered by the TWiki collaboration platformCopyright © 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback