/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

histograms.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 '''
00003 Make some generator level histograms.
00004 
00005 Eg
00006 
00007 nuwa.py --output-stats="{'file1':'stats.root'}" \
00008         [-G $XMLDETDESCROOT/DDDB/dayabay-radslabs.xml]
00009         -n100 -o blah.root \
00010         -m 'GenTools.Opts --helper=radslab' \
00011         -m GenTools.histograms
00012 
00013 
00014 <bv@bnl.gov> Sat Jul 16 02:04:47 2011
00015 '''
00016 
00017 import PyCintex
00018 Gaudi = PyCintex.makeNamespace('Gaudi')
00019 
00020 from DybPython.DybPythonAlg import DybPythonAlg
00021 from DybPython.Util import hists
00022 from GaudiPython import SUCCESS, FAILURE, loaddict
00023 
00024 from DybPython.Util import irange
00025 import GaudiKernel.SystemOfUnits as units
00026 import math
00027 
00028 from ROOT import TH2F
00029 
00030 class HistBase(DybPythonAlg):
00031     '''
00032     Basic histogram handling common to all algs
00033     '''
00034 
00035     class PerSite:
00036         def __init__(self,nick,site_de,ad_des,keeper):
00037             self.nick = nick
00038             self.site_de = site_de
00039             self.ad_de = ad_des
00040             self.keeper = keeper
00041             return
00042         pass
00043 
00044     def __init__(self,name,prefix):
00045         DybPythonAlg.__init__(self,name)
00046         self.prefix = prefix
00047         self.sites = []         # PerSite infor
00048         return
00049 
00050     def initialize(self):
00051         print 'Calling DybPythonAlg.initialize(self)'
00052         sc = DybPythonAlg.initialize(self)
00053         if sc.isFailure(): return sc
00054 
00055         print 'HistBase.initialize()'
00056 
00057         statSvc = self.svc('IStatisticsSvc','StatisticsSvc')
00058 
00059         for nads,nick in [(2,'db'),(2,'la'),(4,'far')]:
00060             keeper = hists.Keeper(statSvc,self.prefix+nick)
00061             site_de_name = '/dd/Structure/Sites/%s-rock' % nick
00062             site_de = self.getDet(site_de_name)
00063             ad_des = []
00064             for iad in range(nads):
00065                 iad += 1
00066                 ad_de_name = '/dd/Structure/AD/%s-oil%d' % (nick,iad)
00067                 ad_de = self.getDet(ad_de_name)
00068                 ad_des.append(ad_de)
00069                 continue                
00070             self.sites.append(HistBase.PerSite(nick,site_de,ad_des,keeper))
00071             continue
00072 
00073         return SUCCESS
00074 
00075     def site_ad_by_point(self, gaudi_point):
00076         '''
00077         Return tuple of (site,ad_de) that contains the given global gaudi point.
00078         If not in a site, both are None.  If not in an AD ad_de is None.
00079         '''
00080 
00081         xyz  = [gaudi_point.x(),gaudi_point.y(),gaudi_point.z()]
00082         xyz_m = map(lambda x: x/units.meter, xyz)
00083         for site in self.sites:
00084             site_gi = site.site_de.geometry()
00085             local_point = site_gi.toLocal(gaudi_point)
00086             if not site_gi.isInside(gaudi_point):
00087                 continue
00088             for ad_de in site.ad_de:
00089                 if ad_de.geometry().isInside(gaudi_point):
00090                     return (site,ad_de)
00091                 continue
00092             return (site,None)
00093             continue
00094         return (None,None)
00095 
00096 class HistVerts(HistBase):
00097     def __init__(self,name='GenVertex',prefix='/file1/gen'):
00098         HistBase.__init__(self,name,prefix)
00099         return
00100 
00101     def initialize(self):
00102         print 'Calling HistBase.initialize(self)'
00103         sc = HistBase.initialize(self)
00104         if sc.isFailure(): return sc
00105 
00106         binsPerMeter = 10
00107         for site in self.sites:
00108             k = site.keeper
00109             nick = site.nick
00110             
00111             #print 'Booking histograms for "%s"' % nick
00112             # Pool
00113             pool_height = 11.0
00114             pool_xsize = 17.0
00115             pool_ysize = 11.0
00116             nads = 2
00117             if nick == 'far':
00118                 pool_ysize = 17.0
00119                 nads = 4
00120             k.poolXY = TH2F('poolXY','%s Pool X-Y Vertices (m)'%nick.upper(),
00121                             pool_xsize*binsPerMeter,-pool_xsize/2,pool_xsize/2,
00122                             pool_ysize*binsPerMeter,-pool_ysize/2,pool_ysize/2)
00123 
00124             k.poolXZ = TH2F('poolXZ','%s Pool X-Z Vertices (m)'%nick.upper(),
00125                             pool_xsize*binsPerMeter,-pool_xsize/2,pool_xsize/2,
00126                             pool_height*binsPerMeter,0,-pool_height)
00127 
00128             k.poolYZ = TH2F('poolYZ','%s Pool Y-Z Vertices (m)'%nick.upper(),
00129                             pool_ysize*binsPerMeter,-pool_ysize/2,pool_ysize/2,
00130                             pool_height*binsPerMeter,0,-pool_height)
00131 
00132             # ADs
00133             ad_radius = 5.0
00134             ad_height = 5.0
00135             for iad in range(nads):
00136                 iad += 1
00137                 ad = 'AD%d' % iad
00138                 #print 'Booking AD #%d/%d' % (iad,nads)
00139 
00140                 nam = ad+'_XY'
00141                 nbinsx = 2*ad_radius*binsPerMeter
00142                 nbinsy = 2*ad_radius*binsPerMeter
00143                 #print 'TH2F %s with %d * %d = %d' % (nam,nbinsx,nbinsy,nbinsx*nbinsy)
00144                 h1 = TH2F(nam,'%s %s X-Y Vertices (m)'%(nick.upper(), ad.upper()),
00145                          nbinsx,-ad_radius,ad_radius,
00146                          nbinsy,-ad_radius,ad_radius)
00147                 k.__setattr__(nam,h1)
00148 
00149                 nam = ad+'_R2Z'
00150                 nbinsx = 2*ad_radius**2*binsPerMeter/10
00151                 nbinsy = ad_height*binsPerMeter
00152                 #print 'TH2F %s with %d * %d = %d' % (nam,nbinsx,nbinsy,nbinsx*nbinsy)
00153                 h2 = TH2F(nam,'%s %s R2Z Vertices (m2 v m)'%(nick.upper(), ad.upper()),
00154                          nbinsx,0,ad_radius**2,
00155                          nbinsy,-ad_height/2.0,ad_height/2.0)
00156                 k.__setattr__(nam,h2)
00157 
00158                 continue        # loop over ADs
00159             continue            # loop over site
00160         return SUCCESS          # initialize()
00161 
00162     def execute(self):
00163         evt = self.evtSvc()
00164         path = '/Event/Gen/GenHeader'
00165         gh = evt[path]
00166         event = gh.event()
00167 
00168         if event is None:
00169             self.error('Failed to get GenHeader from "%s"' % path)
00170             return FAILURE
00171 
00172 
00173         for vtx in irange(event.vertices_begin(),event.vertices_end()):
00174             
00175             pos = vtx.position()
00176             gpoint = Gaudi.XYZPoint(pos.x(),pos.y(),pos.z())
00177             site,ad = self.site_ad_by_point(gpoint)
00178             if not site:
00179                 self.warning('Found vertex not in site: (%s,%s,%s) (meter)' % \
00180                                  (pos.x()/units.meter,
00181                                   pos.y()/units.meter,
00182                                   pos.z()/units.meter))
00183                 continue
00184 
00185             self.fill_vertex(gpoint,site,ad)
00186             continue
00187         return SUCCESS
00188 
00189     def fill_vertex(self,gpoint,site,ad_de = None):
00190         '''
00191         Fill vertex histograms using given gpoint using the site and ad_de
00192         '''
00193         k = site.keeper
00194         site_point = site.site_de.geometry().toLocal(gpoint)
00195         sx_m = site_point.x()/units.meter
00196         sy_m = site_point.y()/units.meter
00197         sz_m = site_point.z()/units.meter
00198         k.poolXY.Fill(sx_m,sy_m)
00199         k.poolXZ.Fill(sx_m,sz_m)
00200         k.poolYZ.Fill(sy_m,sz_m)
00201 
00202         if not ad_de: return
00203 
00204         ad_point = ad_de.geometry().toLocal(gpoint)
00205 
00206         adx_m = ad_point.x()/units.meter
00207         ady_m = ad_point.y()/units.meter
00208         adz_m = ad_point.z()/units.meter
00209 
00210         iad = site.ad_de.index(ad_de) + 1
00211         k.__getattr__('AD%d_XY'%iad).Fill(adx_m,ady_m)
00212         k.__getattr__('AD%d_R2Z'%iad).Fill(adx_m**2+ady_m**2,adz_m)
00213         return
00214 
00215 def configure(argv = None):
00216     return
00217 
00218 def run(app):
00219     app.addAlgorithm(HistVerts())
00220     return
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:18:50 for GenTools by doxygen 1.7.4