| 1 |
import tarfile, re, math, ROOT, sys
|
| 2 |
sys.path.append("/afs/naf.desy.de/user/n/nowaf/UserCode/nowaf/PythonScripts/")
|
| 3 |
import Scan_cff as Scan
|
| 4 |
|
| 5 |
class PointLimit:
|
| 6 |
def __init__( self ):
|
| 7 |
self.m0 = 0
|
| 8 |
self.m12 = 0
|
| 9 |
self.tanb = 0
|
| 10 |
self.CLsObsAs = 0
|
| 11 |
self.CLsExpAs = 0
|
| 12 |
self.CLsExpAs_m2sigma = 0
|
| 13 |
self.CLsExpAs_m1sigma = 0
|
| 14 |
self.CLsExpAs_p2sigma = 0
|
| 15 |
self.CLsExpAs_p1sigma = 0
|
| 16 |
self.CLsObs = 0
|
| 17 |
self.CLsExp = 0
|
| 18 |
self.CLsExp_m2sigma = 0
|
| 19 |
self.CLsExp_m1sigma = 0
|
| 20 |
self.CLsExp_p2sigma = 0
|
| 21 |
self.CLsExp_p1sigma = 0
|
| 22 |
self.CLsObs_err = 0
|
| 23 |
self.CLsExp_err = 0
|
| 24 |
self.CLsExp_m2sigma_err = 0
|
| 25 |
self.CLsExp_m1sigma_err = 0
|
| 26 |
self.CLsExp_p2sigma_err = 0
|
| 27 |
self.CLsExp_p1sigma_err = 0
|
| 28 |
self.CLsExpPl = 0
|
| 29 |
self.CLsExpPl_m2sigma = 0
|
| 30 |
self.CLsExpPl_m1sigma = 0
|
| 31 |
self.CLsExpPl_p2sigma = 0
|
| 32 |
self.CLsExpPl_p1sigma = 0
|
| 33 |
self.acceptable = True
|
| 34 |
pass
|
| 35 |
pass
|
| 36 |
|
| 37 |
def readInLimits( tarDir, tarName ):
|
| 38 |
|
| 39 |
tar = tarfile.open( tarDir + tarName )
|
| 40 |
#print tar.getnames()
|
| 41 |
points = []
|
| 42 |
nInacceptablePoints = 0
|
| 43 |
nPoints = 0
|
| 44 |
for filename in tar.getnames():
|
| 45 |
if not re.search( "result.txt", filename ):
|
| 46 |
continue
|
| 47 |
#print filename
|
| 48 |
f = tar.extractfile( filename )
|
| 49 |
p = PointLimit()
|
| 50 |
filenameSnippets = filename.split( "_" )
|
| 51 |
#print filenameSnippets
|
| 52 |
p.m0 = int( filenameSnippets[ 1 ] )
|
| 53 |
p.m12 = int( filenameSnippets[ 2 ] )
|
| 54 |
filenameSnippetSnippets = filenameSnippets[ 3 ].split( "." )
|
| 55 |
#print filenameSnippetSnippets
|
| 56 |
p.tanb = int( filenameSnippetSnippets[ 0 ] )
|
| 57 |
pointAccepted = True
|
| 58 |
nPoints += 1
|
| 59 |
for line in f.readlines():
|
| 60 |
#print line
|
| 61 |
if not pointAccepted:
|
| 62 |
continue
|
| 63 |
if not re.search( "CLs", line ):
|
| 64 |
continue
|
| 65 |
if re.search( "scanned", line ):
|
| 66 |
continue
|
| 67 |
#print line
|
| 68 |
lineSnippets = line.split( "=" )
|
| 69 |
#print lineSnippets
|
| 70 |
value = lineSnippets[ 1 ].rstrip( " " )
|
| 71 |
try:
|
| 72 |
value = float( value.rstrip( "\n" ) )
|
| 73 |
except:
|
| 74 |
pointAccepted = False
|
| 75 |
value = -1
|
| 76 |
pass
|
| 77 |
if not pointAccepted:
|
| 78 |
continue
|
| 79 |
#print lineSnippets[ 0 ], "=", value
|
| 80 |
if re.search( "asymptotic", lineSnippets[ 0 ] ):
|
| 81 |
if re.search( "observed", lineSnippets[ 0 ] ):
|
| 82 |
p.CLsObsAs = value
|
| 83 |
elif re.search( "expected", lineSnippets[ 0 ] ):
|
| 84 |
if re.search( "m2sigma", lineSnippets[ 0 ] ) :
|
| 85 |
p.CLsExpAs_m2sigma = value
|
| 86 |
elif re.search( "m1sigma", lineSnippets[ 0 ] ) :
|
| 87 |
p.CLsExpAs_m1sigma = value
|
| 88 |
elif re.search( "p2sigma", lineSnippets[ 0 ] ) :
|
| 89 |
p.CLsExpAs_p2sigma = value
|
| 90 |
elif re.search( "p1sigma", lineSnippets[ 0 ] ) :
|
| 91 |
p.CLsExpAs_p1sigma = value
|
| 92 |
else:
|
| 93 |
p.CLsExpAs = value
|
| 94 |
pass
|
| 95 |
pass
|
| 96 |
pass
|
| 97 |
elif re.search( "likelihood", lineSnippets[ 0 ] ):
|
| 98 |
#if re.search( "observed", lineSnippets[ 0 ] ):
|
| 99 |
# p.CLsObsAs = value
|
| 100 |
if re.search( "expected", lineSnippets[ 0 ] ):
|
| 101 |
if re.search( "m2sigma", lineSnippets[ 0 ] ) :
|
| 102 |
p.CLsExpPl_m2sigma = value
|
| 103 |
elif re.search( "m1sigma", lineSnippets[ 0 ] ) :
|
| 104 |
p.CLsExpPl_m1sigma = value
|
| 105 |
elif re.search( "p2sigma", lineSnippets[ 0 ] ) :
|
| 106 |
p.CLsExpPl_p2sigma = value
|
| 107 |
elif re.search( "p1sigma", lineSnippets[ 0 ] ) :
|
| 108 |
p.CLsExpPl_p1sigma = value
|
| 109 |
else:
|
| 110 |
p.CLsExpPl = value
|
| 111 |
pass
|
| 112 |
pass
|
| 113 |
pass
|
| 114 |
else:
|
| 115 |
if re.search( "observed", lineSnippets[ 0 ] ):
|
| 116 |
if re.search( "error", lineSnippets[ 0 ] ):
|
| 117 |
p.CLsObs_err = value
|
| 118 |
else:
|
| 119 |
p.CLsObs = value
|
| 120 |
pass
|
| 121 |
pass
|
| 122 |
elif re.search( "expected", lineSnippets[ 0 ] ):
|
| 123 |
if not re.search( "error", lineSnippets[ 0 ] ):
|
| 124 |
if re.search( "m2sigma", lineSnippets[ 0 ] ) :
|
| 125 |
p.CLsExp_m2sigma = value
|
| 126 |
elif re.search( "m1sigma", lineSnippets[ 0 ] ) :
|
| 127 |
p.CLsExp_m1sigma = value
|
| 128 |
elif re.search( "p2sigma", lineSnippets[ 0 ] ) :
|
| 129 |
p.CLsExp_p2sigma = value
|
| 130 |
elif re.search( "p1sigma", lineSnippets[ 0 ] ) :
|
| 131 |
p.CLsExp_p1sigma = value
|
| 132 |
else:
|
| 133 |
p.CLsExp = value
|
| 134 |
pass
|
| 135 |
pass
|
| 136 |
else:
|
| 137 |
if re.search( "m2sigma", lineSnippets[ 0 ] ) :
|
| 138 |
p.CLsExp_m2sigma_err = value
|
| 139 |
elif re.search( "m1sigma", lineSnippets[ 0 ] ) :
|
| 140 |
p.CLsExp_m1sigma_err = value
|
| 141 |
elif re.search( "p2sigma", lineSnippets[ 0 ] ) :
|
| 142 |
p.CLsExp_p2sigma_err = value
|
| 143 |
elif re.search( "p1sigma", lineSnippets[ 0 ] ) :
|
| 144 |
p.CLsExp_p1sigma_err = value
|
| 145 |
else:
|
| 146 |
p.CLsExp_err = value
|
| 147 |
pass
|
| 148 |
pass
|
| 149 |
pass
|
| 150 |
pass
|
| 151 |
pass
|
| 152 |
#print p.m0
|
| 153 |
#print p.CLsObs
|
| 154 |
if pointAccepted:
|
| 155 |
points.append( p )
|
| 156 |
pass
|
| 157 |
else:
|
| 158 |
nInacceptablePoints += 1
|
| 159 |
print "inacceptable: m0=", p.m0, " m12=", p.m12
|
| 160 |
p.acceptable = False
|
| 161 |
points.append( p )
|
| 162 |
pass
|
| 163 |
print "Of ", nPoints, ", ", nInacceptablePoints, " were not accepted"
|
| 164 |
return points
|
| 165 |
|
| 166 |
def fillMissingPoints( points ):
|
| 167 |
lists = getAxisLists( points )
|
| 168 |
pointsPresent = getPointsPresent( points )
|
| 169 |
#print lists[ 0 ]
|
| 170 |
#print lists[ 1 ]
|
| 171 |
for x in lists[ 0 ]:
|
| 172 |
for y in lists[ 1 ]:
|
| 173 |
found = False
|
| 174 |
#for p in points:
|
| 175 |
# if p.m0 == x and p.m12 == y:
|
| 176 |
# found = True
|
| 177 |
# pass
|
| 178 |
# pass
|
| 179 |
if ( x,y ) in pointsPresent:
|
| 180 |
found = True
|
| 181 |
if not found:
|
| 182 |
#print "Filling up point with m0=", x, " m12=",y
|
| 183 |
#points.append( PointLimit() )
|
| 184 |
newPoint = PointLimit()
|
| 185 |
newPoint.m0 = x
|
| 186 |
newPoint.m12 = y
|
| 187 |
points.append( newPoint )
|
| 188 |
pass
|
| 189 |
pass
|
| 190 |
pass
|
| 191 |
return points
|
| 192 |
|
| 193 |
def getPointsPresent( points ):
|
| 194 |
pointsPresent = []
|
| 195 |
for p in points:
|
| 196 |
pointsPresent.append( ( p.m0, p.m12 ) )
|
| 197 |
pass
|
| 198 |
return pointsPresent
|
| 199 |
|
| 200 |
def getAxisLists( points ):
|
| 201 |
xwidth = getBinWidth( points, "X" )
|
| 202 |
ywidth = getBinWidth( points, "Y" )
|
| 203 |
xextrema = getExtrema( points, "X" )
|
| 204 |
yextrema = getExtrema( points, "Y" )
|
| 205 |
xlist = []
|
| 206 |
ylist = []
|
| 207 |
print "width x=", xwidth, " y=", ywidth
|
| 208 |
for x in range( xextrema[ 0 ], xextrema[ 1 ] + 1, int( xwidth ) ):
|
| 209 |
xlist.append( x )
|
| 210 |
pass
|
| 211 |
for y in range( yextrema[ 0 ], yextrema[ 1 ] + 1, int( ywidth ) ):
|
| 212 |
ylist.append( y )
|
| 213 |
pass
|
| 214 |
return ( xlist, ylist )
|
| 215 |
|
| 216 |
def getExtrema( points, axis ):
|
| 217 |
lowest = 1000000
|
| 218 |
highest = -10000
|
| 219 |
for p in points:
|
| 220 |
value = p.m0
|
| 221 |
if axis == "Y":
|
| 222 |
value = p.m12
|
| 223 |
pass
|
| 224 |
if value < lowest:
|
| 225 |
lowest = value
|
| 226 |
pass
|
| 227 |
if value > highest:
|
| 228 |
highest = value
|
| 229 |
pass
|
| 230 |
pass
|
| 231 |
return ( lowest, highest )
|
| 232 |
|
| 233 |
def getBinWidth( points, axis ):
|
| 234 |
delta = 100000.
|
| 235 |
for p1 in points:
|
| 236 |
for p2 in points:
|
| 237 |
if p1.m0 == p2.m0:
|
| 238 |
continue
|
| 239 |
if p1.m12 == p2.m12:
|
| 240 |
continue
|
| 241 |
if axis == "X":
|
| 242 |
if math.fabs( p1.m0 - p2.m0 ) < delta:
|
| 243 |
delta = math.fabs( p1.m0 - p2.m0 )
|
| 244 |
pass
|
| 245 |
pass
|
| 246 |
elif axis == "Y":
|
| 247 |
if math.fabs( p1.m12 - p2.m12 ) < delta:
|
| 248 |
delta = math.fabs( p1.m12 - p2.m12 )
|
| 249 |
pass
|
| 250 |
pass
|
| 251 |
pass
|
| 252 |
pass
|
| 253 |
return delta
|
| 254 |
|
| 255 |
def findNeighbours( p, points, axis ):
|
| 256 |
lower = None
|
| 257 |
higher = None
|
| 258 |
binWidth = getBinWidth( points, "X" )
|
| 259 |
if axis == "Y":
|
| 260 |
binWidth = getBinWidth( points, "Y" )
|
| 261 |
pass
|
| 262 |
for n in points:
|
| 263 |
if axis == "X":
|
| 264 |
if not n.m12 == p.m12:
|
| 265 |
continue
|
| 266 |
if n.m0 == p.m0 - binWidth:
|
| 267 |
lower = n
|
| 268 |
elif n.m0 == p.m0 + binWidth:
|
| 269 |
higher = n
|
| 270 |
pass
|
| 271 |
pass
|
| 272 |
if axis == "Y":
|
| 273 |
if not n.m0 == p.m0:
|
| 274 |
continue
|
| 275 |
if n.m12 == p.m12 - binWidth:
|
| 276 |
lower = n
|
| 277 |
elif n.m12 == p.m12 + binWidth:
|
| 278 |
higher = n
|
| 279 |
pass
|
| 280 |
pass
|
| 281 |
pass
|
| 282 |
### neighbour stays None if point is at the edge
|
| 283 |
return ( lower, higher )
|
| 284 |
|
| 285 |
def interpolatePoints( points ):
|
| 286 |
points = fillMissingPoints( points )
|
| 287 |
#interpolatePointsInY( points )
|
| 288 |
|
| 289 |
for p in points:
|
| 290 |
if not p.CLsObs == 0:
|
| 291 |
continue
|
| 292 |
noNeighboursY = False
|
| 293 |
noNeighboursX = False
|
| 294 |
neighboursY = findNeighbours( p, points, "Y" )
|
| 295 |
neighboursX = findNeighbours( p, points, "X" )
|
| 296 |
if neighboursY[ 0 ] == None or \
|
| 297 |
neighboursY[ 1 ] == None:
|
| 298 |
noNeighboursY = True
|
| 299 |
pass
|
| 300 |
if neighboursX[ 0 ] == None or \
|
| 301 |
neighboursX[ 1 ] == None:
|
| 302 |
noNeighboursX = True
|
| 303 |
pass
|
| 304 |
if not noNeighboursY and \
|
| 305 |
( neighboursY[ 0 ].CLsObs == 0 or \
|
| 306 |
neighboursY[ 1 ].CLsObs == 0 ):
|
| 307 |
noNeighboursY = True
|
| 308 |
pass
|
| 309 |
if not noNeighboursX and \
|
| 310 |
( neighboursX[ 0 ].CLsObs == 0 or \
|
| 311 |
neighboursX[ 1 ].CLsObs == 0 ):
|
| 312 |
noNeighboursX = True
|
| 313 |
pass
|
| 314 |
if noNeighboursX and noNeighboursY:
|
| 315 |
continue
|
| 316 |
print "will interpolate point m0=", p.m0, " m12=", p.m12
|
| 317 |
if noNeighboursX:
|
| 318 |
print "\tIn Y"
|
| 319 |
interpolatePoints1D( p, neighboursY )
|
| 320 |
pass
|
| 321 |
elif noNeighboursY:
|
| 322 |
print "\tIn X"
|
| 323 |
interpolatePoints1D( p, neighboursX )
|
| 324 |
else:
|
| 325 |
print "\tIn X and Y"
|
| 326 |
### interpolate in both
|
| 327 |
interpolatePoints2D( p, neighboursX, neighboursY )
|
| 328 |
pass
|
| 329 |
|
| 330 |
#print "\tResults:"
|
| 331 |
#print "\t\tCLsObs=", p.CLsObs
|
| 332 |
#print "\t\tCLsExp=", p.CLsExp
|
| 333 |
|
| 334 |
#ind = points.index( p )
|
| 335 |
#points[ ind ] = pnew
|
| 336 |
|
| 337 |
pass
|
| 338 |
return points
|
| 339 |
|
| 340 |
def interpolatePoints1D( p, neighbour ):
|
| 341 |
### only one point at the moment
|
| 342 |
p.CLsObs = neighbour[ 0 ].CLsObs + 1/2. * ( neighbour[ 1 ].CLsObs - neighbour[ 0 ].CLsObs )
|
| 343 |
p.CLsObsAs = neighbour[ 0 ].CLsObsAs + 1/2. * ( neighbour[ 1 ].CLsObsAs - neighbour[ 0 ].CLsObsAs )
|
| 344 |
p.CLsExpAs = neighbour[ 0 ].CLsExpAs + 1/2. * ( neighbour[ 1 ].CLsExpAs - neighbour[ 0 ].CLsExpAs )
|
| 345 |
p.CLsExpAs_m2sigma = neighbour[ 0 ].CLsExpAs_m2sigma + 1/2. * ( neighbour[ 1 ].CLsExpAs_m2sigma - neighbour[ 0 ].CLsExpAs_m2sigma )
|
| 346 |
p.CLsExpAs_m1sigma = neighbour[ 0 ].CLsExpAs_m1sigma + 1/2. * ( neighbour[ 1 ].CLsExpAs_m1sigma - neighbour[ 0 ].CLsExpAs_m1sigma )
|
| 347 |
p.CLsExpAs_p2sigma = neighbour[ 0 ].CLsExpAs_p2sigma + 1/2. * ( neighbour[ 1 ].CLsExpAs_p2sigma - neighbour[ 0 ].CLsExpAs_p2sigma )
|
| 348 |
p.CLsExpAs_p1sigma = neighbour[ 0 ].CLsExpAs_p1sigma + 1/2. * ( neighbour[ 1 ].CLsExpAs_p1sigma - neighbour[ 0 ].CLsExpAs_p1sigma )
|
| 349 |
p.CLsExp = neighbour[ 0 ].CLsExp + 1/2. * ( neighbour[ 1 ].CLsExp - neighbour[ 0 ].CLsExp )
|
| 350 |
p.CLsExp_m2sigma = neighbour[ 0 ].CLsExp_m2sigma + 1/2. * ( neighbour[ 1 ].CLsExp_m2sigma - neighbour[ 0 ].CLsExp_m2sigma )
|
| 351 |
p.CLsExp_m1sigma = neighbour[ 0 ].CLsExp_m1sigma + 1/2. * ( neighbour[ 1 ].CLsExp_m1sigma - neighbour[ 0 ].CLsExp_m1sigma )
|
| 352 |
p.CLsExp_p2sigma = neighbour[ 0 ].CLsExp_p2sigma + 1/2. * ( neighbour[ 1 ].CLsExp_p2sigma - neighbour[ 0 ].CLsExp_p2sigma )
|
| 353 |
p.CLsExp_p1sigma = neighbour[ 0 ].CLsExp_p1sigma + 1/2. * ( neighbour[ 1 ].CLsExp_p1sigma - neighbour[ 0 ].CLsExp_p1sigma )
|
| 354 |
p.CLsObs_err = neighbour[ 0 ].CLsObs_err + 1/2. * ( neighbour[ 1 ].CLsObs_err - neighbour[ 0 ].CLsObs_err )
|
| 355 |
p.CLsExp_err = neighbour[ 0 ].CLsExp_err + 1/2. * ( neighbour[ 1 ].CLsExp_err - neighbour[ 0 ].CLsExp_err )
|
| 356 |
p.CLsExp_m2sigma_err = neighbour[ 0 ].CLsExp_m2sigma_err + 1/2. * ( neighbour[ 1 ].CLsExp_m2sigma_err - neighbour[ 0 ].CLsExp_m2sigma_err )
|
| 357 |
p.CLsExp_m1sigma_err = neighbour[ 0 ].CLsExp_m1sigma_err + 1/2. * ( neighbour[ 1 ].CLsExp_m1sigma_err - neighbour[ 0 ].CLsExp_m1sigma_err )
|
| 358 |
p.CLsExp_p2sigma_err = neighbour[ 0 ].CLsExp_p2sigma_err + 1/2. * ( neighbour[ 1 ].CLsExp_p2sigma_err - neighbour[ 0 ].CLsExp_p2sigma_err )
|
| 359 |
p.CLsExp_p1sigma_err = neighbour[ 0 ].CLsExp_p1sigma_err + 1/2. * ( neighbour[ 1 ].CLsExp_p1sigma_err - neighbour[ 0 ].CLsExp_p1sigma_err )
|
| 360 |
pass
|
| 361 |
|
| 362 |
def interpolatePoints2D( p, neighbours1, neighbours2 ):
|
| 363 |
p.CLsObs = ( neighbours1[ 0 ].CLsObs + 1/2. * ( neighbours1[ 1 ].CLsObs - neighbours1[ 0 ].CLsObs ) + \
|
| 364 |
neighbours2[ 0 ].CLsObs + 1/2. * ( neighbours2[ 1 ].CLsObs - neighbours2[ 0 ].CLsObs ) ) / 2.
|
| 365 |
p.CLsObsAs = ( neighbours1[ 0 ].CLsObsAs + 1/2. * ( neighbours1[ 1 ].CLsObsAs - neighbours1[ 0 ].CLsObsAs ) + \
|
| 366 |
neighbours2[ 0 ].CLsObsAs + 1/2. * ( neighbours1[ 1 ].CLsObsAs - neighbours2[ 0 ].CLsObsAs ) ) / 2.
|
| 367 |
p.CLsExpAs = ( neighbours1[ 0 ].CLsExpAs + 1/2. * ( neighbours1[ 1 ].CLsExpAs - neighbours1[ 0 ].CLsExpAs ) + \
|
| 368 |
neighbours2[ 0 ].CLsExpAs + 1/2. * ( neighbours2[ 1 ].CLsExpAs - neighbours2[ 0 ].CLsExpAs ) ) / 2.
|
| 369 |
p.CLsExpAs_m2sigma = ( neighbours1[ 0 ].CLsExpAs_m2sigma + 1/2. * ( neighbours1[ 1 ].CLsExpAs_m2sigma - neighbours1[ 0 ].CLsExpAs_m2sigma ) + \
|
| 370 |
neighbours2[ 0 ].CLsExpAs_m2sigma + 1/2. * ( neighbours2[ 1 ].CLsExpAs_m2sigma - neighbours2[ 0 ].CLsExpAs_m2sigma ) ) / 2.
|
| 371 |
p.CLsExpAs_m1sigma = ( neighbours1[ 0 ].CLsExpAs_m1sigma + 1/2. * ( neighbours1[ 1 ].CLsExpAs_m1sigma - neighbours1[ 0 ].CLsExpAs_m1sigma ) + \
|
| 372 |
neighbours2[ 0 ].CLsExpAs_m1sigma + 1/2. * ( neighbours2[ 1 ].CLsExpAs_m1sigma - neighbours2[ 0 ].CLsExpAs_m1sigma ) ) / 2.
|
| 373 |
p.CLsExpAs_p2sigma = ( neighbours1[ 0 ].CLsExpAs_p2sigma + 1/2. * ( neighbours1[ 1 ].CLsExpAs_p2sigma - neighbours1[ 0 ].CLsExpAs_p2sigma ) + \
|
| 374 |
neighbours2[ 0 ].CLsExpAs_p2sigma + 1/2. * ( neighbours2[ 1 ].CLsExpAs_p2sigma - neighbours2[ 0 ].CLsExpAs_p2sigma ) ) / 2.
|
| 375 |
p.CLsExpAs_p1sigma = ( neighbours1[ 0 ].CLsExpAs_p1sigma + 1/2. * ( neighbours1[ 1 ].CLsExpAs_p1sigma - neighbours1[ 0 ].CLsExpAs_p1sigma ) + \
|
| 376 |
neighbours2[ 0 ].CLsExpAs_p1sigma + 1/2. * ( neighbours2[ 1 ].CLsExpAs_p1sigma - neighbours2[ 0 ].CLsExpAs_p1sigma ) ) / 2.
|
| 377 |
p.CLsExp = ( neighbours1[ 0 ].CLsExp + 1/2. * ( neighbours1[ 1 ].CLsExp - neighbours1[ 0 ].CLsExp ) + \
|
| 378 |
neighbours2[ 0 ].CLsExp + 1/2. * ( neighbours2[ 1 ].CLsExp - neighbours2[ 0 ].CLsExp ) ) / 2.
|
| 379 |
p.CLsExp_m2sigma = ( neighbours1[ 0 ].CLsExp_m2sigma + 1/2. * ( neighbours1[ 1 ].CLsExp_m2sigma - neighbours1[ 0 ].CLsExp_m2sigma ) + \
|
| 380 |
neighbours2[ 0 ].CLsExp_m2sigma + 1/2. * ( neighbours2[ 1 ].CLsExp_m2sigma - neighbours2[ 0 ].CLsExp_m2sigma ) ) / 2.
|
| 381 |
p.CLsExp_m1sigma = ( neighbours1[ 0 ].CLsExp_m1sigma + 1/2. * ( neighbours1[ 1 ].CLsExp_m1sigma - neighbours1[ 0 ].CLsExp_m1sigma ) + \
|
| 382 |
neighbours2[ 0 ].CLsExp_m1sigma + 1/2. * ( neighbours2[ 1 ].CLsExp_m1sigma - neighbours2[ 0 ].CLsExp_m1sigma ) ) / 2.
|
| 383 |
p.CLsExp_p2sigma = ( neighbours1[ 0 ].CLsExp_p2sigma + 1/2. * ( neighbours1[ 1 ].CLsExp_p2sigma - neighbours1[ 0 ].CLsExp_p2sigma ) + \
|
| 384 |
neighbours2[ 0 ].CLsExp_p2sigma + 1/2. * ( neighbours2[ 1 ].CLsExp_p2sigma - neighbours2[ 0 ].CLsExp_p2sigma ) ) /2.
|
| 385 |
p.CLsExp_p1sigma = ( neighbours1[ 0 ].CLsExp_p1sigma + 1/2. * ( neighbours1[ 1 ].CLsExp_p1sigma - neighbours1[ 0 ].CLsExp_p1sigma ) + \
|
| 386 |
neighbours2[ 0 ].CLsExp_p1sigma + 1/2. * ( neighbours2[ 1 ].CLsExp_p1sigma - neighbours2[ 0 ].CLsExp_p1sigma ) ) / 2.
|
| 387 |
p.CLsObs_err = ( neighbours1[ 0 ].CLsObs_err + 1/2. * ( neighbours1[ 1 ].CLsObs_err - neighbours1[ 0 ].CLsObs_err ) + \
|
| 388 |
neighbours2[ 0 ].CLsObs_err + 1/2. * ( neighbours2[ 1 ].CLsObs_err - neighbours2[ 0 ].CLsObs_err ) ) / 2.
|
| 389 |
p.CLsExp_err = ( neighbours1[ 0 ].CLsExp_err + 1/2. * ( neighbours1[ 1 ].CLsExp_err - neighbours1[ 0 ].CLsExp_err ) + \
|
| 390 |
neighbours2[ 0 ].CLsExp_err + 1/2. * ( neighbours2[ 1 ].CLsExp_err - neighbours2[ 0 ].CLsExp_err ) ) / 2.
|
| 391 |
p.CLsExp_m2sigma_err = ( neighbours1[ 0 ].CLsExp_m2sigma_err + 1/2. * ( neighbours1[ 1 ].CLsExp_m2sigma_err - neighbours1[ 0 ].CLsExp_m2sigma_err ) + \
|
| 392 |
neighbours2[ 0 ].CLsExp_m2sigma_err + 1/2. * ( neighbours2[ 1 ].CLsExp_m2sigma_err - neighbours2[ 0 ].CLsExp_m2sigma_err ) ) / 2.
|
| 393 |
p.CLsExp_m1sigma_err = ( neighbours1[ 0 ].CLsExp_m1sigma_err + 1/2. * ( neighbours1[ 1 ].CLsExp_m1sigma_err - neighbours1[ 0 ].CLsExp_m1sigma_err ) + \
|
| 394 |
neighbours2[ 0 ].CLsExp_m1sigma_err + 1/2. * ( neighbours2[ 1 ].CLsExp_m1sigma_err - neighbours2[ 0 ].CLsExp_m1sigma_err ) ) / 2.
|
| 395 |
p.CLsExp_p2sigma_err = ( neighbours1[ 0 ].CLsExp_p2sigma_err + 1/2. * ( neighbours1[ 1 ].CLsExp_p2sigma_err - neighbours1[ 0 ].CLsExp_p2sigma_err ) + \
|
| 396 |
neighbours2[ 0 ].CLsExp_p2sigma_err + 1/2. * ( neighbours2[ 1 ].CLsExp_p2sigma_err - neighbours2[ 0 ].CLsExp_p2sigma_err ) ) / 2.
|
| 397 |
p.CLsExp_p1sigma_err = ( neighbours1[ 0 ].CLsExp_p1sigma_err + 1/2. * ( neighbours1[ 1 ].CLsExp_p1sigma_err - neighbours1[ 0 ].CLsExp_p1sigma_err ) + \
|
| 398 |
neighbours2[ 0 ].CLsExp_p1sigma_err + 1/2. * ( neighbours2[ 1 ].CLsExp_p1sigma_err - neighbours2[ 0 ].CLsExp_p1sigma_err ) ) / 2.
|
| 399 |
pass
|
| 400 |
|
| 401 |
def readInSanjaysUncerts( histDict, them0, them12, filename ):
|
| 402 |
|
| 403 |
#xbin = histDict[ "Process" ].GetXaxis().FindBin( them0 )
|
| 404 |
#ybin = histDict[ "Process" ].GetYaxis().FindBin( them12 )
|
| 405 |
#zSlice = histDict[ "Process" ].ProjectionZ( "_pz", xbin, xbin, ybin, ybin )
|
| 406 |
|
| 407 |
|
| 408 |
### 0 1 2 3 4 5 6 7 8 9
|
| 409 |
### ng ns nn ll sb ss tb bb gg sg
|
| 410 |
#pDict = {}
|
| 411 |
#sumbin = histDict[ "ScanAnalyzerAllMC" ].FindBin( them0, them12 )
|
| 412 |
#sum = histDict[ "ScanAnalyzerAllMC" ].GetBinContent( sumbin )
|
| 413 |
#for zbin in range( 0, zSlice.GetNbinsX() ):
|
| 414 |
# pDict[ zbin ] = zSlice.GetBinContent( zbin )
|
| 415 |
|
| 416 |
# #### percentage of decays
|
| 417 |
# pDict[ zbin ] /= sum
|
| 418 |
# pass
|
| 419 |
|
| 420 |
#### read in uncertainties
|
| 421 |
uDict = {}
|
| 422 |
f = open( filename, "r" )
|
| 423 |
sum = 0
|
| 424 |
uncert = 0
|
| 425 |
notfound = True
|
| 426 |
for line in f.readlines():
|
| 427 |
### ignore first line
|
| 428 |
#if re.search( "Sub-processes", line ):
|
| 429 |
if re.search( "Interactions", line ):
|
| 430 |
continue
|
| 431 |
lineSnippets = line.split( "|" )
|
| 432 |
pointDefs = lineSnippets[ 1 ].split( " " )
|
| 433 |
m0 = int( pointDefs[ 6 ] )
|
| 434 |
if not m0 == them0:
|
| 435 |
continue
|
| 436 |
m12 = int( pointDefs[ 9 ] )
|
| 437 |
if not m12 == them12:
|
| 438 |
continue
|
| 439 |
notfound = False
|
| 440 |
### 2 3 4 5 6 7 8 9 10 11
|
| 441 |
### ng ns nn ll sb ss tb bb gg sg
|
| 442 |
for u in range( 2, 12 ):
|
| 443 |
vals = lineSnippets[ u ].split( " pm " )
|
| 444 |
#uDict[ u ] = ( float( vals[ 1 ] )/ float( vals [ 0 ] * pDict[ u - 2 ] * float( vals [ 0 ] )**2
|
| 445 |
#uncert += ( float( vals[ 1 ] )/ float( vals [ 0 ] ) * pDict[ u - 2 ] * float( vals [ 0 ] ) )**2
|
| 446 |
#sum += pDict[ u - 2 ] * float( vals [ 0 ] )
|
| 447 |
uDict[ u ] = float( vals[ 1 ] )
|
| 448 |
sum += float( vals [ 0 ] )
|
| 449 |
pass
|
| 450 |
pass
|
| 451 |
f.close()
|
| 452 |
uncert = 0
|
| 453 |
### add up lineraly because of correlation
|
| 454 |
if not notfound:
|
| 455 |
for i in range( 2, 12 ):
|
| 456 |
uncert += uDict[ i ]
|
| 457 |
pass
|
| 458 |
|
| 459 |
return uncert / sum
|
| 460 |
return None
|
| 461 |
|
| 462 |
def getRelTheoUncert( m0, m12, histDict ):
|
| 463 |
#### do something here
|
| 464 |
|
| 465 |
|
| 466 |
### toy uncerts
|
| 467 |
#nloUp = 0.2
|
| 468 |
#nloDn = 0.3
|
| 469 |
#pdfUp = 0.12
|
| 470 |
#pdfDn = 0.12
|
| 471 |
|
| 472 |
#### nlo
|
| 473 |
#bin = histDict[ "RateRel0NLO2" ].FindBin( m0, m12 )
|
| 474 |
#nloUp = histDict[ "RateRel0NLO2" ].GetBinContent( bin )
|
| 475 |
#nloDn = histDict[ "RateRel0NLO05" ].GetBinContent( bin )
|
| 476 |
#pdfUp = histDict[ "RateRel0PDFxsecUncert" ].GetBinContent( bin )
|
| 477 |
#pdfDn = pdfUp
|
| 478 |
|
| 479 |
|
| 480 |
#print "m0=", m0, "m12=", m12, " nloUp=", nloUp, " pdfUp=", pdfUp
|
| 481 |
|
| 482 |
#up = math.sqrt( nloUp**2 + pdfUp**2 )
|
| 483 |
#dn = math.sqrt( nloDn**2 + pdfDn**2 )
|
| 484 |
|
| 485 |
filename = "./theoUncerts/combined_cross_section_Errmsugra_m0_m12_10_0_1.txt"
|
| 486 |
up = readInSanjaysUncerts( histDict, m0, m12, filename )
|
| 487 |
dn = up
|
| 488 |
|
| 489 |
#print "m0=", m0, "m12=", m12, " Up=", up
|
| 490 |
|
| 491 |
return up, dn
|
| 492 |
|
| 493 |
def fillTheoErrPlot( hist ):
|
| 494 |
filename = "./theoUncerts/combined_cross_section_Errmsugra_m0_m12_10_0_1.txt"
|
| 495 |
for m0 in range( 220, 3020, 20 ):
|
| 496 |
for m12 in range( 100, 1020, 20 ):
|
| 497 |
val = readInSanjaysUncerts( {}, m0, m12, filename )
|
| 498 |
bin = hist.FindBin( m0, m12 )
|
| 499 |
if not val == None:
|
| 500 |
hist.SetBinContent( bin, val )
|
| 501 |
pass
|
| 502 |
pass
|
| 503 |
pass
|
| 504 |
pass
|
| 505 |
|
| 506 |
def makeStandardPlots( points, file, lumi ):
|
| 507 |
|
| 508 |
histDict = {}
|
| 509 |
|
| 510 |
#### nlo uncert
|
| 511 |
histDict.update( Scan.readInScanHists( file, list=[ "ScanAnalyzerNLO2",
|
| 512 |
"ScanAnalyzerPreNLO2",
|
| 513 |
"ScanAnalyzerNLO05",
|
| 514 |
"ScanAnalyzerPreNLO05",
|
| 515 |
"ScanAnalyzerNLO",
|
| 516 |
"ScanAnalyzerPreNLO",
|
| 517 |
"ScanAnalyzerAllNLO",
|
| 518 |
"ScanAnalyzerAllMC",
|
| 519 |
"ScanAnalyzerPDFxsec",
|
| 520 |
"ScanAnalyzerPDFxsecUncert",
|
| 521 |
] ) )
|
| 522 |
Scan.getRateHist( histDict, "NLO", lumi=lumi, Type="" )
|
| 523 |
Scan.getRateHist( histDict, "NLO05", lumi=lumi, Type="" )
|
| 524 |
Scan.getRateHist( histDict, "NLO2", lumi=lumi, Type="" )
|
| 525 |
Scan.getRelRate( histDict, "NLO05", "NLO", Type="" )
|
| 526 |
Scan.getRelRate( histDict, "NLO2", "NLO", Type="" )
|
| 527 |
Scan.getRateHist( histDict, "PDFxsec", lumi=lumi, Type="" )
|
| 528 |
Scan.getRateHist( histDict, "PDFxsecUncert", lumi=lumi, Type="" )
|
| 529 |
Scan.getRelRate( histDict, "PDFxsecUncert", "PDFxsec", Type="" )
|
| 530 |
|
| 531 |
histDict[ "Process" ] = file.Get( "ScanAnalyzerAllNLO/process" )
|
| 532 |
histDict[ "Process" ].SetDirectory( 0 )
|
| 533 |
|
| 534 |
### now create first plots ( cMSSM )
|
| 535 |
### evaluate later
|
| 536 |
nBinsX = 140
|
| 537 |
nBinsY = 46
|
| 538 |
xmin = 210.
|
| 539 |
xmax = 3010.
|
| 540 |
ymin = 90.
|
| 541 |
ymax = 1010.
|
| 542 |
h_obs = ROOT.TH2F( "obs", "obs", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 543 |
h_obs_eval = ROOT.TH2F( "obsEv", "obsEv", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 544 |
h_exp = ROOT.TH2F( "exp", "exp", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 545 |
h_exp_eval = ROOT.TH2F( "expEv", "expEv", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 546 |
h_exp_p1 = ROOT.TH2F( "expP1", "expP1", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 547 |
h_exp_m1 = ROOT.TH2F( "expM1", "expM1", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 548 |
h_exp_p1_eval = ROOT.TH2F( "expP1Ev", "expP1Ev", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 549 |
h_exp_m1_eval = ROOT.TH2F( "expM1Ev", "expM1Ev", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 550 |
|
| 551 |
h_inacc = ROOT.TH2F( "inacc", "inacc", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 552 |
h_inacc.SetFillColor( ROOT.kRed )
|
| 553 |
|
| 554 |
#### for the theretical uncertainties
|
| 555 |
h_obs_theoUp_eval = ROOT.TH2F( "obsTheoUpEv", "obsTheoUpEv", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 556 |
h_obs_theoDn_eval = ROOT.TH2F( "obsTheoDnEv", "obsTheoDnEv", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 557 |
h_exp_theoUp_eval = ROOT.TH2F( "expTheoUpEv", "expTheoUpEv", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 558 |
h_exp_theoDn_eval = ROOT.TH2F( "expTheoDnEv", "expTheoDnEv", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 559 |
h_obs_theoUp = ROOT.TH2F( "obsTheoUp", "obsTheoUp", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 560 |
h_obs_theoDn = ROOT.TH2F( "obsTheoDn", "obsTheoDn", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 561 |
h_exp_theoUp = ROOT.TH2F( "expTheoUp", "expTheoUp", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 562 |
h_exp_theoDn = ROOT.TH2F( "expTheoDn", "expTheoDn", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 563 |
|
| 564 |
h_theorel = ROOT.TH2F( "reltheo", "reltheo", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 565 |
|
| 566 |
#for xbin in range( 1, h_obs.GetNbinsX() + 1 ):
|
| 567 |
# for ybin in range( 1, h_obs.GetNbinsY() + 1 ):
|
| 568 |
for p in points:
|
| 569 |
|
| 570 |
if p.CLsObs == 0:
|
| 571 |
continue
|
| 572 |
|
| 573 |
bin = h_obs.FindBin( p.m0, p.m12 )
|
| 574 |
|
| 575 |
#if p.acceptable:
|
| 576 |
|
| 577 |
#### theoretical uncerts
|
| 578 |
relUp, relDn = getRelTheoUncert( p.m0, p.m12, histDict )
|
| 579 |
|
| 580 |
h_obs.SetBinContent( bin, p.CLsObs )
|
| 581 |
h_exp.SetBinContent( bin, p.CLsExp )
|
| 582 |
h_exp_m1.SetBinContent( bin, p.CLsExp_m1sigma )
|
| 583 |
h_exp_p1.SetBinContent( bin, p.CLsExp_p1sigma )
|
| 584 |
h_obs_theoUp.SetBinContent( bin, p.CLsObs * ( 1 + relUp ) )
|
| 585 |
h_obs_theoDn.SetBinContent( bin, p.CLsObs * ( 1 - relUp ) )
|
| 586 |
h_exp_theoUp.SetBinContent( bin, p.CLsExp * ( 1 + relUp ) )
|
| 587 |
h_exp_theoDn.SetBinContent( bin, p.CLsExp * ( 1 - relUp ) )
|
| 588 |
#h_theorel.SetBinContent( bin, relUp )
|
| 589 |
|
| 590 |
####### temporary Asymptotoc!
|
| 591 |
#print "######################### WARNING! Using Asymptotic !! ############################"
|
| 592 |
|
| 593 |
val = 0.01
|
| 594 |
if p.CLsObs > 1:
|
| 595 |
val = 1.
|
| 596 |
pass
|
| 597 |
if not p.CLsObs == 0:
|
| 598 |
h_obs_eval.SetBinContent( bin, val )
|
| 599 |
pass
|
| 600 |
val = 0.01
|
| 601 |
if p.CLsExp > 1:
|
| 602 |
val = 1.
|
| 603 |
pass
|
| 604 |
if not p.CLsExp == 0:
|
| 605 |
h_exp_eval.SetBinContent( bin, val )
|
| 606 |
pass
|
| 607 |
val = 0.01
|
| 608 |
if p.CLsExp_m1sigma > 1:
|
| 609 |
val = 1
|
| 610 |
pass
|
| 611 |
if not p.CLsExp_m1sigma == 0:
|
| 612 |
h_exp_m1_eval.SetBinContent( bin, val )
|
| 613 |
pass
|
| 614 |
val = 0.01
|
| 615 |
if p.CLsExp_p1sigma > 1:
|
| 616 |
val = 1
|
| 617 |
pass
|
| 618 |
if not p.CLsExp_p1sigma == 0:
|
| 619 |
h_exp_p1_eval.SetBinContent( bin, val )
|
| 620 |
pass
|
| 621 |
|
| 622 |
if not p.acceptable:
|
| 623 |
h_inacc.SetBinContent( bin, 40 )
|
| 624 |
pass
|
| 625 |
|
| 626 |
###### now, the theory errs
|
| 627 |
val = 0.01
|
| 628 |
if p.CLsObs * ( 1 + relUp ) > 1:
|
| 629 |
val = 1.
|
| 630 |
pass
|
| 631 |
if not p.CLsObs * ( 1 + relUp ) == 0:
|
| 632 |
h_obs_theoUp_eval.SetBinContent( bin, val )
|
| 633 |
pass
|
| 634 |
val = 0.01
|
| 635 |
if p.CLsObs * ( 1 - relDn ) > 1:
|
| 636 |
val = 1.
|
| 637 |
pass
|
| 638 |
if not p.CLsObs * ( 1 - relDn ) == 0:
|
| 639 |
h_obs_theoDn_eval.SetBinContent( bin, val )
|
| 640 |
pass
|
| 641 |
val = 0.01
|
| 642 |
if p.CLsExp * ( 1 + relUp ) > 1:
|
| 643 |
val = 1.
|
| 644 |
pass
|
| 645 |
if not p.CLsExp * ( 1 + relUp ) == 0:
|
| 646 |
h_exp_theoUp_eval.SetBinContent( bin, val )
|
| 647 |
pass
|
| 648 |
val = 0.01
|
| 649 |
if p.CLsExp * ( 1 - relDn ) > 1:
|
| 650 |
val = 1.
|
| 651 |
pass
|
| 652 |
if not p.CLsExp * ( 1 - relDn ) == 0:
|
| 653 |
h_exp_theoDn_eval.SetBinContent( bin, val )
|
| 654 |
pass
|
| 655 |
|
| 656 |
pass
|
| 657 |
|
| 658 |
histDict[ "Obs" ] = h_obs.Clone()
|
| 659 |
histDict[ "ObsEval" ] = h_obs_eval.Clone()
|
| 660 |
histDict[ "Exp" ] = h_exp.Clone()
|
| 661 |
histDict[ "ExpEval" ] = h_exp_eval.Clone()
|
| 662 |
histDict[ "Inacc" ] = h_inacc.Clone()
|
| 663 |
histDict[ "ExpM1" ] = h_exp_m1.Clone()
|
| 664 |
histDict[ "ExpP1" ] = h_exp_p1.Clone()
|
| 665 |
histDict[ "ExpM1Eval" ] = h_exp_m1_eval.Clone()
|
| 666 |
histDict[ "ExpP1Eval" ] = h_exp_p1_eval.Clone()
|
| 667 |
histDict[ "ObsTheoUpEval" ] = h_obs_theoUp_eval.Clone()
|
| 668 |
histDict[ "ObsTheoDnEval" ] = h_obs_theoDn_eval.Clone()
|
| 669 |
histDict[ "ExpTheoUpEval" ] = h_exp_theoUp_eval.Clone()
|
| 670 |
histDict[ "ExpTheoDnEval" ] = h_exp_theoDn_eval.Clone()
|
| 671 |
histDict[ "ObsTheoUp" ] = h_obs_theoUp.Clone()
|
| 672 |
histDict[ "ObsTheoDn" ] = h_obs_theoDn.Clone()
|
| 673 |
histDict[ "ExpTheoUp" ] = h_exp_theoUp.Clone()
|
| 674 |
histDict[ "ExpTheoDn" ] = h_exp_theoDn.Clone()
|
| 675 |
histDict[ "TheoRel" ] = h_theorel.Clone()
|
| 676 |
#fillTheoErrPlot( histDict[ "TheoRel" ] )
|
| 677 |
|
| 678 |
#### now smooth the 2D non evaluated plots
|
| 679 |
smoo = 3
|
| 680 |
histDict[ "ObsFullSmooth" ] = h_obs.Clone()
|
| 681 |
histDict[ "ObsFullSmooth" ].Smooth( smoo )
|
| 682 |
histDict[ "ObsFullSmooth" ].SetName( "ObsFullSmooth" )
|
| 683 |
histDict[ "ExpFullSmooth" ] = h_exp.Clone()
|
| 684 |
histDict[ "ExpFullSmooth" ].Smooth( smoo )
|
| 685 |
histDict[ "ExpFullSmooth" ].SetName( "ExpFullSmooth" )
|
| 686 |
histDict[ "ExpM1FullSmooth" ] = h_exp_m1.Clone()
|
| 687 |
histDict[ "ExpM1FullSmooth" ].Smooth( smoo )
|
| 688 |
histDict[ "ExpM1FullSmooth" ].SetName( "ExpM1FullSmooth" )
|
| 689 |
histDict[ "ExpP1FullSmooth" ] = h_exp_p1.Clone()
|
| 690 |
histDict[ "ExpP1FullSmooth" ].Smooth( smoo )
|
| 691 |
histDict[ "ExpP1FullSmooth" ].SetName( "ExpP1FullSmooth" )
|
| 692 |
histDict[ "ObsTheoUpFullSmooth" ] = h_obs_theoUp.Clone()
|
| 693 |
histDict[ "ObsTheoUpFullSmooth" ].Smooth( smoo )
|
| 694 |
histDict[ "ObsTheoUpFullSmooth" ].SetName( "ObsTheoUpFullSmooth" )
|
| 695 |
histDict[ "ObsTheoDnFullSmooth" ] = h_obs_theoDn.Clone()
|
| 696 |
histDict[ "ObsTheoDnFullSmooth" ].Smooth( smoo )
|
| 697 |
histDict[ "ObsTheoDnFullSmooth" ].SetName( "ObsTheoDnFullSmooth" )
|
| 698 |
histDict[ "ExpTheoUpFullSmooth" ] = h_exp_theoUp.Clone()
|
| 699 |
histDict[ "ExpTheoUpFullSmooth" ].Smooth( smoo )
|
| 700 |
histDict[ "ExpTheoUpFullSmooth" ].SetName( "ExpTheoUpFullSmooth" )
|
| 701 |
histDict[ "ExpTheoDnFullSmooth" ] = h_exp_theoDn.Clone()
|
| 702 |
histDict[ "ExpTheoDnFullSmooth" ].Smooth( smoo )
|
| 703 |
histDict[ "ExpTheoDnFullSmooth" ].SetName( "ExpTheoDnFullSmooth" )
|
| 704 |
|
| 705 |
|
| 706 |
h_obs_smoo = ROOT.TH2F( "obsSmooth", "obsSmooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 707 |
h_obs_smoo_eval = ROOT.TH2F( "obsSmoothEval", "obsSmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 708 |
h_exp_smoo = ROOT.TH2F( "expSmooth", "expSmooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 709 |
h_exp_smoo_eval = ROOT.TH2F( "expSmoothEval", "expSmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 710 |
h_exp_m1_smoo = ROOT.TH2F( "expM1Smooth", "expM1Smooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 711 |
h_exp_m1_smoo_eval = ROOT.TH2F( "expM1SmoothEval", "expM1SmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 712 |
h_exp_p1_smoo = ROOT.TH2F( "expP1Smooth", "expP1Smooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 713 |
h_exp_p1_smoo_eval = ROOT.TH2F( "expP1SmoothEval", "expP1SmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 714 |
h_obs_theoUp_smoo = ROOT.TH2F( "obsTheoUpSmooth", "obsTheoUpSmooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 715 |
h_obs_theoUp_smoo_eval = ROOT.TH2F( "obsTheoUpSmoothEval", "obsTheoUpSmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 716 |
h_obs_theoDn_smoo = ROOT.TH2F( "obsTheoDnSmooth", "obsTheoDnSmooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 717 |
h_obs_theoDn_smoo_eval = ROOT.TH2F( "obsTheoDnSmoothEval", "obsTheoDnSmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 718 |
h_exp_theoUp_smoo = ROOT.TH2F( "expTheoUpSmooth", "expTheoUpSmooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 719 |
h_exp_theoUp_smoo_eval = ROOT.TH2F( "expTheoUpSmoothEval", "expTheoUpSmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 720 |
h_exp_theoDn_smoo = ROOT.TH2F( "expTheoDnSmooth", "expTheoDnSmooth", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 721 |
h_exp_theoDn_smoo_eval = ROOT.TH2F( "expTheoDnSmoothEval", "expTheoDnSmoothEval", nBinsX, xmin, xmax, nBinsY, ymin, ymax )
|
| 722 |
#### refill smoothing plots
|
| 723 |
#nBinsX = 140
|
| 724 |
#nBinsY = 46
|
| 725 |
#xmin = 210.
|
| 726 |
#xmax = 3010.
|
| 727 |
#ymin = 90.
|
| 728 |
#ymax = 1010.
|
| 729 |
for xval in range( xmin, xmax + ( xmax - xmin )/ ( nBinsX * 2 ), ( xmax - xmin )/ ( nBinsX ) ):
|
| 730 |
for yval in range( ymin, ymax + ( ymax - ymin )/ ( nBinsY * 2 ), ( ymax - ymin )/ ( nBinsY ) ):
|
| 731 |
bin = histDict[ "ObsFullSmooth" ].FindBin( xval, yval )
|
| 732 |
if yval <= 340 or xval > 460:
|
| 733 |
h_obs_smoo.SetBinContent( bin, histDict[ "ObsFullSmooth" ].GetBinContent( bin ) )
|
| 734 |
h_exp_smoo.SetBinContent( bin, histDict[ "ExpFullSmooth" ].GetBinContent( bin ) )
|
| 735 |
h_exp_m1_smoo.SetBinContent( bin, histDict[ "ExpM1FullSmooth" ].GetBinContent( bin ) )
|
| 736 |
h_exp_p1_smoo.SetBinContent( bin, histDict[ "ExpP1FullSmooth" ].GetBinContent( bin ) )
|
| 737 |
h_obs_theoUp_smoo.SetBinContent( bin, histDict[ "ObsTheoUpFullSmooth" ].GetBinContent( bin ) )
|
| 738 |
h_obs_theoDn_smoo.SetBinContent( bin, histDict[ "ObsTheoDnFullSmooth" ].GetBinContent( bin ) )
|
| 739 |
h_exp_theoUp_smoo.SetBinContent( bin, histDict[ "ExpTheoUpFullSmooth" ].GetBinContent( bin ) )
|
| 740 |
h_exp_theoDn_smoo.SetBinContent( bin, histDict[ "ExpTheoDnFullSmooth" ].GetBinContent( bin ) )
|
| 741 |
|
| 742 |
#####
|
| 743 |
val = 0.01
|
| 744 |
if histDict[ "ObsFullSmooth" ].GetBinContent( bin ) > 1:
|
| 745 |
val = 1
|
| 746 |
pass
|
| 747 |
h_obs_smoo_eval.SetBinContent( bin, val )
|
| 748 |
val = 0.01
|
| 749 |
if histDict[ "ExpFullSmooth" ].GetBinContent( bin ) > 1:
|
| 750 |
val = 1
|
| 751 |
pass
|
| 752 |
h_exp_smoo_eval.SetBinContent( bin, val )
|
| 753 |
val = 0.01
|
| 754 |
if histDict[ "ExpM1FullSmooth" ].GetBinContent( bin ) > 1:
|
| 755 |
val = 1
|
| 756 |
pass
|
| 757 |
h_exp_m1_smoo_eval.SetBinContent( bin, val )
|
| 758 |
val = 0.01
|
| 759 |
if histDict[ "ExpP1FullSmooth" ].GetBinContent( bin ) > 1:
|
| 760 |
val = 1
|
| 761 |
pass
|
| 762 |
h_exp_p1_smoo_eval.SetBinContent( bin, val )
|
| 763 |
val = 0.01
|
| 764 |
if histDict[ "ObsTheoUpFullSmooth" ].GetBinContent( bin ) > 1:
|
| 765 |
val = 1
|
| 766 |
pass
|
| 767 |
h_obs_theoUp_smoo_eval.SetBinContent( bin, val )
|
| 768 |
val = 0.01
|
| 769 |
if histDict[ "ObsTheoDnFullSmooth" ].GetBinContent( bin ) > 1:
|
| 770 |
val = 1
|
| 771 |
pass
|
| 772 |
h_obs_theoDn_smoo_eval.SetBinContent( bin, val )
|
| 773 |
val = 0.01
|
| 774 |
if histDict[ "ExpTheoUpFullSmooth" ].GetBinContent( bin ) > 1:
|
| 775 |
val = 1
|
| 776 |
pass
|
| 777 |
h_exp_theoUp_smoo_eval.SetBinContent( bin, val )
|
| 778 |
val = 0.01
|
| 779 |
if histDict[ "ExpTheoDnFullSmooth" ].GetBinContent( bin ) > 1:
|
| 780 |
val = 1
|
| 781 |
pass
|
| 782 |
h_exp_theoDn_smoo_eval.SetBinContent( bin, val )
|
| 783 |
|
| 784 |
else:
|
| 785 |
h_obs_smoo.SetBinContent( bin, histDict[ "Obs" ].GetBinContent( bin ) )
|
| 786 |
h_exp_smoo.SetBinContent( bin, histDict[ "Exp" ].GetBinContent( bin ) )
|
| 787 |
h_exp_m1_smoo.SetBinContent( bin, histDict[ "ExpM1" ].GetBinContent( bin ) )
|
| 788 |
h_exp_p1_smoo.SetBinContent( bin, histDict[ "ExpP1" ].GetBinContent( bin ) )
|
| 789 |
h_obs_theoUp_smoo.SetBinContent( bin, histDict[ "ObsTheoUp" ].GetBinContent( bin ) )
|
| 790 |
h_obs_theoDn_smoo.SetBinContent( bin, histDict[ "ObsTheoDn" ].GetBinContent( bin ) )
|
| 791 |
h_exp_theoUp_smoo.SetBinContent( bin, histDict[ "ExpTheoUp" ].GetBinContent( bin ) )
|
| 792 |
h_exp_theoDn_smoo.SetBinContent( bin, histDict[ "ExpTheoDn" ].GetBinContent( bin ) )
|
| 793 |
|
| 794 |
####
|
| 795 |
val = 0.01
|
| 796 |
if histDict[ "Obs" ].GetBinContent( bin ) > 1:
|
| 797 |
val = 1
|
| 798 |
pass
|
| 799 |
h_obs_smoo_eval.SetBinContent( bin, val )
|
| 800 |
val = 0.01
|
| 801 |
if histDict[ "Exp" ].GetBinContent( bin ) > 1:
|
| 802 |
val = 1
|
| 803 |
pass
|
| 804 |
h_exp_smoo_eval.SetBinContent( bin, val )
|
| 805 |
val = 0.01
|
| 806 |
if histDict[ "ExpM1" ].GetBinContent( bin ) > 1:
|
| 807 |
val = 1
|
| 808 |
pass
|
| 809 |
h_exp_m1_smoo_eval.SetBinContent( bin, val )
|
| 810 |
val = 0.01
|
| 811 |
if histDict[ "ExpP1" ].GetBinContent( bin ) > 1:
|
| 812 |
val = 1
|
| 813 |
pass
|
| 814 |
h_exp_p1_smoo_eval.SetBinContent( bin, val )
|
| 815 |
val = 0.01
|
| 816 |
if histDict[ "ObsTheoUp" ].GetBinContent( bin ) > 1:
|
| 817 |
val = 1
|
| 818 |
pass
|
| 819 |
h_obs_theoUp_smoo_eval.SetBinContent( bin, val )
|
| 820 |
val = 0.01
|
| 821 |
if histDict[ "ObsTheoDn" ].GetBinContent( bin ) > 1:
|
| 822 |
val = 1
|
| 823 |
pass
|
| 824 |
h_obs_theoDn_smoo_eval.SetBinContent( bin, val )
|
| 825 |
val = 0.01
|
| 826 |
if histDict[ "ExpTheoUp" ].GetBinContent( bin ) > 1:
|
| 827 |
val = 1
|
| 828 |
pass
|
| 829 |
h_exp_theoUp_smoo_eval.SetBinContent( bin, val )
|
| 830 |
val = 0.01
|
| 831 |
if histDict[ "ExpTheoDn" ].GetBinContent( bin ) > 1:
|
| 832 |
val = 1
|
| 833 |
pass
|
| 834 |
h_exp_theoDn_smoo_eval.SetBinContent( bin, val )
|
| 835 |
pass
|
| 836 |
pass
|
| 837 |
pass
|
| 838 |
|
| 839 |
histDict[ "ObsSmooth" ] = h_obs_smoo.Clone()
|
| 840 |
histDict[ "ObsSmoothEval" ] = h_obs_smoo_eval.Clone()
|
| 841 |
histDict[ "ExpSmooth" ] = h_exp_smoo.Clone()
|
| 842 |
histDict[ "ExpSmoothEval" ] = h_exp_smoo_eval.Clone()
|
| 843 |
histDict[ "ExpM1Smooth" ] = h_exp_m1_smoo.Clone()
|
| 844 |
histDict[ "ExpM1SmoothEval" ] = h_exp_m1_smoo_eval.Clone()
|
| 845 |
histDict[ "ExpP1Smooth" ] = h_exp_p1_smoo.Clone()
|
| 846 |
histDict[ "ExpP1SmoothEval" ] = h_exp_p1_smoo_eval.Clone()
|
| 847 |
histDict[ "ObsTheoUpSmooth" ] = h_obs_theoUp_smoo.Clone()
|
| 848 |
histDict[ "ObsTheoUpSmoothEval" ] = h_obs_theoUp_smoo_eval.Clone()
|
| 849 |
histDict[ "ObsTheoDnSmooth" ] = h_obs_theoDn_smoo.Clone()
|
| 850 |
histDict[ "ObsTheoDnSmoothEval" ] = h_obs_theoDn_smoo_eval.Clone()
|
| 851 |
histDict[ "ExpTheoUpSmooth" ] = h_exp_theoUp_smoo.Clone()
|
| 852 |
histDict[ "ExpTheoUpSmoothEval" ] = h_exp_theoUp_smoo_eval.Clone()
|
| 853 |
histDict[ "ExpTheoDnSmooth" ] = h_exp_theoDn_smoo.Clone()
|
| 854 |
histDict[ "ExpTheoDnSmoothEval" ] = h_exp_theoDn_smoo_eval.Clone()
|
| 855 |
|
| 856 |
return histDict
|
| 857 |
|
| 858 |
def getContours( histDict, hist="ObsEval" ):
|
| 859 |
obsClone = histDict[ hist ].Clone()
|
| 860 |
obsClone.SetContour( 3 )
|
| 861 |
### needed
|
| 862 |
obsClone.Draw( "CONT Z List" )
|
| 863 |
### needed, else the contous object is not found
|
| 864 |
ROOT.gPad.Update()
|
| 865 |
|
| 866 |
contours = ROOT.gROOT.GetListOfSpecials().FindObject( "contours" )
|
| 867 |
#print "contours size: ", contours.GetSize()
|
| 868 |
|
| 869 |
#print "---------------", hist, "------------------"
|
| 870 |
resultingGraphs = []
|
| 871 |
for i in range( 0, contours.GetSize() ):
|
| 872 |
cList = contours.At( i )
|
| 873 |
#print "i=", i, "cList size=", cList.GetSize()
|
| 874 |
curve = cList.First()
|
| 875 |
if curve:
|
| 876 |
#print type( curve )
|
| 877 |
resultingGraphs.append( curve.Clone() )
|
| 878 |
pass
|
| 879 |
for j in range( 0, cList.GetSize() ):
|
| 880 |
curve = cList.After( curve )
|
| 881 |
if curve:
|
| 882 |
#print type( curve )
|
| 883 |
resultingGraphs.append( curve.Clone() )
|
| 884 |
pass
|
| 885 |
pass
|
| 886 |
pass
|
| 887 |
#curv.Draw( "AP" )
|
| 888 |
#histDict[ "test" ] = obsClone
|
| 889 |
### sorting ...
|
| 890 |
## not done yet
|
| 891 |
#print resultingGraphs
|
| 892 |
histDict[ "contourList" + hist ] = resultingGraphs
|
| 893 |
#print "===== hist=", hist, " ", type( histDict[ "contourList" + hist ][ 0 ] )
|
| 894 |
pass
|
| 895 |
|
| 896 |
def drawContourList( histDict, cDict, hist="ObsEval" ):
|
| 897 |
for i in range( 0, len( histDict[ "contourList" + hist ] ) ):
|
| 898 |
cDict[ "cList" + str( i ) ] = ROOT.TCanvas( "cList" + str( i ), "cList" + str( i ) )
|
| 899 |
cDict[ "cList" + str( i ) ].cd()
|
| 900 |
histDict[ "contourList" + hist ][ i ].Draw( "Al" )
|
| 901 |
pass
|
| 902 |
pass
|
| 903 |
|
| 904 |
def drawAllContours( histDict, cDict, hist="ObsEval" ):
|
| 905 |
cDict[ "cList" ] = ROOT.TCanvas( "cList" , "cList" )
|
| 906 |
cDict[ "cList" ].cd()
|
| 907 |
### find graph with most points:
|
| 908 |
indexPoints = 0
|
| 909 |
nPoints = 0
|
| 910 |
for j in range( 0, len( histDict[ "contourList" + hist ] ) ):
|
| 911 |
if histDict[ "contourList" + hist ][ j ].GetN() > nPoints:
|
| 912 |
nPoints = histDict[ "contourList" + hist ][ j ].GetN()
|
| 913 |
indexPoints = j
|
| 914 |
pass
|
| 915 |
pass
|
| 916 |
histDict[ "contourList" + hist ][ indexPoints ].Draw( "Al" )
|
| 917 |
histDict[ "contourList" + hist ][ indexPoints ].SetLineColor( ROOT.kAzure - 6 )
|
| 918 |
histDict[ "contourList" + hist ][ indexPoints ].SetLineWidth( 3 )
|
| 919 |
for i in range( 0, len( histDict[ "contourList" + hist ] ) ):
|
| 920 |
histDict[ "contourList" + hist ][ i ].SetLineColor( ROOT.kAzure - 6 )
|
| 921 |
histDict[ "contourList" + hist ][ i ].SetLineWidth( 3 )
|
| 922 |
if i == indexPoints:
|
| 923 |
continue
|
| 924 |
else:
|
| 925 |
histDict[ "contourList" + hist ][ i ].Draw( "samel" )
|
| 926 |
pass
|
| 927 |
pass
|
| 928 |
pass
|
| 929 |
|
| 930 |
## def findPointIndex( graph, tuple ):
|
| 931 |
## for i in range( 0, graph.GetN() ):
|
| 932 |
## x = ROOT.Double( 0 )
|
| 933 |
## y = ROOT.Double( 0 )
|
| 934 |
## graph.GetPoint( i, x, y )
|
| 935 |
## if ( x, y ) == tuple:
|
| 936 |
## return i
|
| 937 |
## pass
|
| 938 |
## return -1
|
| 939 |
|
| 940 |
def smooth( graph, N, ymin = 280 ):
|
| 941 |
old = graph.Clone()
|
| 942 |
if N > 2 * graph.GetN():
|
| 943 |
N = 2 * graph.GetN() - 1
|
| 944 |
pass
|
| 945 |
#print "N=", N
|
| 946 |
|
| 947 |
gauss = [ 0 ] *N
|
| 948 |
sigma = N / 4.
|
| 949 |
|
| 950 |
#print "sigma=", sigma
|
| 951 |
|
| 952 |
sum = 0
|
| 953 |
lim = N / 2.
|
| 954 |
|
| 955 |
#print "lim=", lim
|
| 956 |
|
| 957 |
fb = ROOT.TF1( "fb", "gaus(0)", -lim, lim )
|
| 958 |
fb.SetParameter( 0, 1. / ( math.sqrt( 2 * 3.1415 ) * sigma ) )
|
| 959 |
fb.SetParameter( 1, 0 )
|
| 960 |
fb.SetParameter( 2, sigma )
|
| 961 |
for i in range( 0, N ):
|
| 962 |
gauss[i] = fb.Integral( -lim + i, -lim + i + 1 )
|
| 963 |
sum += gauss[ i ]
|
| 964 |
pass
|
| 965 |
#print "sum=", sum
|
| 966 |
|
| 967 |
for i in range( 0, N ):
|
| 968 |
gauss[ i ] /= sum
|
| 969 |
#print "gauss[ i ]=", gauss[ i ]
|
| 970 |
pass
|
| 971 |
|
| 972 |
#print "i-loop"
|
| 973 |
for i in range( 0, graph.GetN() ):
|
| 974 |
#print "\ti=", i
|
| 975 |
avy = 0.
|
| 976 |
avx = 0.
|
| 977 |
|
| 978 |
### use only 3 points to each side
|
| 979 |
### if m0 < 450
|
| 980 |
x1 = ROOT.Double( 0 )
|
| 981 |
y1 = ROOT.Double( 0 )
|
| 982 |
old.GetPoint( i, x1, y1 )
|
| 983 |
##if x1 == 0:
|
| 984 |
## continue
|
| 985 |
#if y1 > 460:
|
| 986 |
# r = 6
|
| 987 |
#elif x1 < 450:
|
| 988 |
# r = 6
|
| 989 |
# pass
|
| 990 |
#else :
|
| 991 |
r = N
|
| 992 |
|
| 993 |
x = ROOT.Double( 0 )
|
| 994 |
x0 = 0
|
| 995 |
y = ROOT.Double( 0 )
|
| 996 |
points = 0
|
| 997 |
#print "\tj-loop"
|
| 998 |
for j in range( i - r / 2,i + r / 2 ):
|
| 999 |
#print "\t\tj=", j
|
| 1000 |
|
| 1001 |
if x1 == 0:
|
| 1002 |
continue
|
| 1003 |
|
| 1004 |
if j < 0:
|
| 1005 |
old.GetPoint( 0, x, y )
|
| 1006 |
#old.GetPoint( 1, x, y )
|
| 1007 |
elif j >= graph.GetN():
|
| 1008 |
old.GetPoint( old.GetN() - 1, x, y )
|
| 1009 |
else:
|
| 1010 |
old.GetPoint( j, x, y )
|
| 1011 |
pass
|
| 1012 |
#print "\t\tx=", x, " y=", y
|
| 1013 |
|
| 1014 |
if i == j:
|
| 1015 |
x0 = x1
|
| 1016 |
#y0 = y ####?
|
| 1017 |
pass
|
| 1018 |
avy += y * gauss[ points ]
|
| 1019 |
avx += x * gauss[ points ]
|
| 1020 |
points += 1
|
| 1021 |
|
| 1022 |
#if x == 0:
|
| 1023 |
# break
|
| 1024 |
|
| 1025 |
pass
|
| 1026 |
|
| 1027 |
#print "\tx=", x1, " y=", y1
|
| 1028 |
#print "\tavx=", avx, " avy=", avy
|
| 1029 |
if x1 == 0:
|
| 1030 |
graph.SetPoint( i, x1, y1 )
|
| 1031 |
#print "SetPoint(", x1, y1, ")"
|
| 1032 |
|
| 1033 |
elif y1 == ymin:
|
| 1034 |
graph.SetPoint( i, avx, y1 )
|
| 1035 |
#elif avy < ymin:
|
| 1036 |
# graph.SetPoint( i, avx, y1 )
|
| 1037 |
|
| 1038 |
#elif x > 426. and x < 427. and y > 399. and y < 460.: ####temporary?
|
| 1039 |
# graph.SetPoint( i, x, y )
|
| 1040 |
elif i - r / 2 < 0 or i + r / 2 >= graph.GetN():
|
| 1041 |
#graph.SetPoint( i, x0, avy )
|
| 1042 |
graph.SetPoint( i, avx, avy )
|
| 1043 |
#print "SetPoint(", x0, avy, ")"
|
| 1044 |
#graph.SetPoint( i, x0, y0 ) ###?
|
| 1045 |
#elif avx < 250:
|
| 1046 |
# graph.SetPoint( i, x0, avy )
|
| 1047 |
else:
|
| 1048 |
graph.SetPoint( i, avx, avy )
|
| 1049 |
#print "SetPoint(", avx, avy, ")"
|
| 1050 |
pass
|
| 1051 |
pass
|
| 1052 |
pass
|
| 1053 |
|
| 1054 |
|
| 1055 |
def getFullContour( histDict, hist="ObsEval", color=ROOT.kAzure - 6, xmax=790, ymax=510, ymin = 280 ):
|
| 1056 |
#cDict[ "cList" ] = ROOT.TCanvas( "cList" , "cList" )
|
| 1057 |
#cDict[ "cList" ].cd()
|
| 1058 |
### find graph with most points:
|
| 1059 |
print "-------------", hist, "--------------------"
|
| 1060 |
indexPoints = 0
|
| 1061 |
nPoints = 0
|
| 1062 |
for j in range( 0, len( histDict[ "contourList" + hist ] ) ):
|
| 1063 |
#print type( histDict[ "contourList" + hist ][ j ] )
|
| 1064 |
if histDict[ "contourList" + hist ][ j ].GetN() > nPoints:
|
| 1065 |
nPoints = histDict[ "contourList" + hist ][ j ].GetN()
|
| 1066 |
indexPoints = j
|
| 1067 |
pass
|
| 1068 |
pass
|
| 1069 |
### make new graph
|
| 1070 |
#get the x,y tuple first
|
| 1071 |
pointTuples = []
|
| 1072 |
for i in range( 0, histDict[ "contourList" + hist ][ indexPoints ].GetN() ):
|
| 1073 |
x = ROOT.Double( 0 )
|
| 1074 |
y = ROOT.Double( 0 )
|
| 1075 |
histDict[ "contourList" + hist ][ indexPoints ].GetPoint( i, x, y )
|
| 1076 |
pointTuples.append( ( x,y ) )
|
| 1077 |
pass
|
| 1078 |
## now find the index of the list, where y is greatest and x is smallest
|
| 1079 |
# now find the index of the list, where x is greatest and y is smallest
|
| 1080 |
#y = -1
|
| 1081 |
#ind = 0
|
| 1082 |
#x = 100000
|
| 1083 |
#for pt in pointTuples:
|
| 1084 |
# if pt[ 1 ] > y:
|
| 1085 |
# #x = pt[ 0 ]
|
| 1086 |
# y = pt[ 1 ]
|
| 1087 |
# pass
|
| 1088 |
# pass
|
| 1089 |
#for pt in pointTuples:
|
| 1090 |
# if pt[ 1 ] < y:
|
| 1091 |
# continue
|
| 1092 |
# if pt[ 0 ] <= x:
|
| 1093 |
# x = pt[ 0 ]
|
| 1094 |
# ind = pointTuples.index( pt )
|
| 1095 |
# pass
|
| 1096 |
# pass
|
| 1097 |
y = 10000
|
| 1098 |
ind = 0
|
| 1099 |
x = -1
|
| 1100 |
for pt in pointTuples:
|
| 1101 |
if pt[ 0 ] > x:
|
| 1102 |
#x = pt[ 0 ]
|
| 1103 |
x = pt[ 0 ]
|
| 1104 |
pass
|
| 1105 |
pass
|
| 1106 |
for pt in pointTuples:
|
| 1107 |
if pt[ 0 ] < x:
|
| 1108 |
continue
|
| 1109 |
if pt[ 1 ] <= y:
|
| 1110 |
y = pt[ 1 ]
|
| 1111 |
ind = pointTuples.index( pt )
|
| 1112 |
pass
|
| 1113 |
pass
|
| 1114 |
print hist
|
| 1115 |
#print "greatest x=", x
|
| 1116 |
#print "smallest y=", y
|
| 1117 |
|
| 1118 |
### split tuple List in list before ind and after ind:
|
| 1119 |
list1 = pointTuples[ : ind ]
|
| 1120 |
list2 = pointTuples[ ind : ]
|
| 1121 |
|
| 1122 |
#print "------------", hist, "------------------"
|
| 1123 |
|
| 1124 |
#print "list1: ", list1
|
| 1125 |
#print "list2: ", list2
|
| 1126 |
|
| 1127 |
### now fill new graph
|
| 1128 |
cont = ROOT.TGraph()
|
| 1129 |
p = 0
|
| 1130 |
x = 100000
|
| 1131 |
lastTuple = ( None, None )
|
| 1132 |
start = False
|
| 1133 |
#print "List 2"
|
| 1134 |
for pt2 in list2:
|
| 1135 |
if pt2[ 0 ] >= xmax:
|
| 1136 |
continue
|
| 1137 |
elif pt2[ 1 ] >= ymax:
|
| 1138 |
continue
|
| 1139 |
elif pt2[ 1 ] <= ymin:
|
| 1140 |
continue
|
| 1141 |
if pt2[ 0 ] == 0 and pt2[ 1 ] == 0:
|
| 1142 |
continue
|
| 1143 |
if not start:
|
| 1144 |
cont.SetPoint( p, pt2[ 0 ], ymin )
|
| 1145 |
#print "\tAdding Point (", pt2[ 0 ], ", ", ymin, ")"
|
| 1146 |
p += 1
|
| 1147 |
start = True
|
| 1148 |
pass
|
| 1149 |
|
| 1150 |
|
| 1151 |
|
| 1152 |
cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] )
|
| 1153 |
#print "\tAdding Point (", pt2[ 0 ], ", ", pt2[ 1 ], ")"
|
| 1154 |
|
| 1155 |
|
| 1156 |
if pt2[ 0 ] < 480 and pt2[ 0 ] > 420:
|
| 1157 |
#print "setting additional points"
|
| 1158 |
p += 1
|
| 1159 |
cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] + 0.0001 )
|
| 1160 |
p +=1
|
| 1161 |
cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] + 0.0002 )
|
| 1162 |
p +=1
|
| 1163 |
cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] + 0.0003 )
|
| 1164 |
#p +=1
|
| 1165 |
#cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] + 0.0004 )
|
| 1166 |
#p +=1
|
| 1167 |
#cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] + 0.0005 )
|
| 1168 |
#p +=1
|
| 1169 |
#cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] + 0.0006 )
|
| 1170 |
pass
|
| 1171 |
|
| 1172 |
p +=1
|
| 1173 |
lastTuple = ( pt2[ 0 ], pt2[ 1 ] )
|
| 1174 |
pass
|
| 1175 |
#print "List1"
|
| 1176 |
start = False
|
| 1177 |
for pt1 in list1:
|
| 1178 |
if pt1[ 0 ] >= xmax:
|
| 1179 |
continue
|
| 1180 |
elif pt1[ 1 ] >= ymax:
|
| 1181 |
continue
|
| 1182 |
elif pt1[ 1 ] <= ymin:
|
| 1183 |
continue
|
| 1184 |
if pt1[ 0 ] == 0 and pt1[ 1 ] == 0:
|
| 1185 |
continue
|
| 1186 |
|
| 1187 |
if not start:
|
| 1188 |
cont.SetPoint( p, pt1[ 0 ], ymin )
|
| 1189 |
print "\tAdding Point (", pt1[ 0 ], ", ", ymin, ")"
|
| 1190 |
|
| 1191 |
p += 1
|
| 1192 |
start = True
|
| 1193 |
pass
|
| 1194 |
|
| 1195 |
cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] )
|
| 1196 |
#print "\tAdding Point (", pt1[ 0 ], ", ", pt1[ 1 ], ")"
|
| 1197 |
|
| 1198 |
if pt1[ 0 ] < 480 and pt1[ 0 ] > 420:
|
| 1199 |
#print "setting additional points"
|
| 1200 |
p += 1
|
| 1201 |
cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] + 0.0001 )
|
| 1202 |
p +=1
|
| 1203 |
cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] + 0.0002 )
|
| 1204 |
p +=1
|
| 1205 |
cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] + 0.0003 )
|
| 1206 |
#p +=1
|
| 1207 |
#cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] + 0.0004 )
|
| 1208 |
#p +=1
|
| 1209 |
#cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] + 0.0005 )
|
| 1210 |
#p +=1
|
| 1211 |
#cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] + 0.0006 )
|
| 1212 |
pass
|
| 1213 |
|
| 1214 |
|
| 1215 |
|
| 1216 |
p +=1
|
| 1217 |
lastTuple = ( pt1[ 0 ], pt1[ 1 ] )
|
| 1218 |
pass
|
| 1219 |
#### extrapolate last tuple to m0 = 0
|
| 1220 |
#print lastTuple
|
| 1221 |
#if re.search( "1", hist ):
|
| 1222 |
cont.SetPoint( p, 0 , lastTuple[ 1 ] )
|
| 1223 |
# pass
|
| 1224 |
|
| 1225 |
cont.SetLineColor( color )
|
| 1226 |
cont.SetLineWidth( 2 )
|
| 1227 |
#smooth( cont, 15 )
|
| 1228 |
smooth( cont, 15, ymin )
|
| 1229 |
histDict[ "cont" + hist ] = cont
|
| 1230 |
#cont.Draw( "AL" )
|
| 1231 |
pass
|
| 1232 |
|
| 1233 |
def getContourSnippet( histDict, hist="ObsEval", color=ROOT.kAzure -6, xmax=640, ymax=240 ):
|
| 1234 |
cont = ROOT.TGraph()
|
| 1235 |
startPointFound = False
|
| 1236 |
p = 0
|
| 1237 |
#print "---------------contSnippet----------------"
|
| 1238 |
for i in range( 0, histDict[ "cont" + hist ].GetN() ):
|
| 1239 |
x = ROOT.Double( 0 )
|
| 1240 |
y = ROOT.Double( 0 )
|
| 1241 |
histDict[ "cont" + hist ].GetPoint( i, x, y )
|
| 1242 |
if x <= xmax and y >= ymax and not startPointFound:
|
| 1243 |
#if y>ymax:
|
| 1244 |
startPointFound = True
|
| 1245 |
pass
|
| 1246 |
if startPointFound:
|
| 1247 |
#print i, " x=", x, " y=", y
|
| 1248 |
cont.SetPoint( p, x, y )
|
| 1249 |
p +=1
|
| 1250 |
pass
|
| 1251 |
pass
|
| 1252 |
cont.SetLineColor( color )
|
| 1253 |
cont.SetLineWidth( 2 )
|
| 1254 |
histDict[ "contSnippet" + hist ] = cont
|
| 1255 |
pass
|
| 1256 |
|
| 1257 |
def getFullContoursWithBand( histDict ):
|
| 1258 |
Exp = histDict[ "contExpSmoothEval" ].Clone()
|
| 1259 |
Exp.SetLineWidth( 2 )
|
| 1260 |
Exp.SetLineColor( ROOT.kAzure - 6 )
|
| 1261 |
#Exp.SetLineStyle( 7 )
|
| 1262 |
histDict[ "FullBandExp" ] = Exp
|
| 1263 |
Obs = histDict[ "contObsSmoothEval" ].Clone()
|
| 1264 |
Obs.SetLineWidth( 2 )
|
| 1265 |
Obs.SetLineColor( 1 )
|
| 1266 |
histDict[ "FullBandObs" ] = Obs
|
| 1267 |
|
| 1268 |
### now make the band
|
| 1269 |
band = ROOT.TGraph( histDict[ "contExpM1SmoothEval" ].GetN() + histDict[ "contExpP1SmoothEval" ].GetN() + 1 )
|
| 1270 |
|
| 1271 |
print "number of expected points=", histDict[ "contExpM1SmoothEval" ].GetN() + histDict[ "contExpP1SmoothEval" ].GetN() + 1
|
| 1272 |
p = 0;
|
| 1273 |
firstx = ROOT.Double( 0 )
|
| 1274 |
firsty = ROOT.Double( 0 )
|
| 1275 |
for i in range( 0, histDict[ "contExpM1SmoothEval" ].GetN() + 1 ):
|
| 1276 |
x = ROOT.Double( 0 )
|
| 1277 |
y = ROOT.Double( 0 )
|
| 1278 |
histDict[ "contExpM1SmoothEval" ].GetPoint( i, x, y )
|
| 1279 |
if x == 0 and y == 0:
|
| 1280 |
continue
|
| 1281 |
band.SetPoint( p, x, y )
|
| 1282 |
print p, " x=",x, " y=",y
|
| 1283 |
p += 1
|
| 1284 |
|
| 1285 |
if i == 0:
|
| 1286 |
firstx = x
|
| 1287 |
firsty = y
|
| 1288 |
pass
|
| 1289 |
pass
|
| 1290 |
for i in range( histDict[ "contExpP1SmoothEval" ].GetN() , 0 - 1, -1 ):
|
| 1291 |
x = ROOT.Double( 0 )
|
| 1292 |
y = ROOT.Double( 0 )
|
| 1293 |
histDict[ "contExpP1SmoothEval" ].GetPoint( i, x, y )
|
| 1294 |
if x == 0 and y == 0:
|
| 1295 |
continue
|
| 1296 |
band.SetPoint( p, x, y )
|
| 1297 |
print p, " x=",x, " y=",y
|
| 1298 |
p += 1
|
| 1299 |
pass
|
| 1300 |
band.SetPoint( p, firstx, firsty )
|
| 1301 |
print "used points=", p
|
| 1302 |
print "firstx=", firstx, " firsty=", firsty
|
| 1303 |
p += 1
|
| 1304 |
band.SetLineColor( ROOT.kAzure - 9 )
|
| 1305 |
band.SetFillColor( ROOT.kAzure - 9 )
|
| 1306 |
band.SetFillStyle( 1001 )
|
| 1307 |
histDict[ "FullBand" ] = band
|
| 1308 |
pass
|
| 1309 |
|
| 1310 |
def getSnippetContoursWithBand( histDict ):
|
| 1311 |
Exp = histDict[ "contSnippetExpEval" ].Clone()
|
| 1312 |
Exp.SetLineWidth( 2 )
|
| 1313 |
Exp.SetLineColor( 1 )
|
| 1314 |
Exp.SetLineStyle( 7 )
|
| 1315 |
histDict[ "SnippetBandExp" ] = Exp
|
| 1316 |
Obs = histDict[ "contSnippetObsEval" ].Clone()
|
| 1317 |
Obs.SetLineWidth( 2 )
|
| 1318 |
Obs.SetLineColor( 1 )
|
| 1319 |
histDict[ "SnippetBandObs" ] = Obs
|
| 1320 |
|
| 1321 |
### now make the band
|
| 1322 |
band = ROOT.TGraph( histDict[ "contSnippetExpM1Eval" ].GetN() + histDict[ "contSnippetExpP1Eval" ].GetN() + 1 )
|
| 1323 |
|
| 1324 |
print "number of expected points=", histDict[ "contSnippetExpM1Eval" ].GetN() + histDict[ "contSnippetExpP1Eval" ].GetN() + 1
|
| 1325 |
p = 0;
|
| 1326 |
firstx = ROOT.Double( 0 )
|
| 1327 |
firsty = ROOT.Double( 0 )
|
| 1328 |
for i in range( 0, histDict[ "contSnippetExpM1Eval" ].GetN() ):
|
| 1329 |
x = ROOT.Double( 0 )
|
| 1330 |
y = ROOT.Double( 0 )
|
| 1331 |
histDict[ "contSnippetExpM1Eval" ].GetPoint( i, x, y )
|
| 1332 |
if x == 0 and y == 0:
|
| 1333 |
continue
|
| 1334 |
band.SetPoint( p, x, y )
|
| 1335 |
print p, " x=",x, " y=",y
|
| 1336 |
p += 1
|
| 1337 |
|
| 1338 |
if i == 0:
|
| 1339 |
firstx = x
|
| 1340 |
firsty = y
|
| 1341 |
pass
|
| 1342 |
pass
|
| 1343 |
for i in range( histDict[ "contSnippetExpP1Eval" ].GetN() + 1 , 0, -1 ):
|
| 1344 |
x = ROOT.Double( 0 )
|
| 1345 |
y = ROOT.Double( 0 )
|
| 1346 |
histDict[ "contSnippetExpP1Eval" ].GetPoint( i, x, y )
|
| 1347 |
if x == 0 and y == 0:
|
| 1348 |
continue
|
| 1349 |
band.SetPoint( p, x, y )
|
| 1350 |
print p, " x=",x, " y=",y
|
| 1351 |
p += 1
|
| 1352 |
pass
|
| 1353 |
band.SetPoint( p, firstx, firsty )
|
| 1354 |
print "used points=", p
|
| 1355 |
print "firstx=", firstx, " firsty=", firsty
|
| 1356 |
p += 1
|
| 1357 |
band.SetLineColor( 1 )
|
| 1358 |
band.SetFillColor( ROOT.kYellow )
|
| 1359 |
band.SetFillStyle( 1001 )
|
| 1360 |
histDict[ "SnippetBand" ] = band
|
| 1361 |
pass
|
| 1362 |
|
| 1363 |
## def getContourSnippet( histDict, hist="ObsEval", color=ROOT.kAzure + 7 ):
|
| 1364 |
## #cDict[ "cList" ] = ROOT.TCanvas( "cList" , "cList" )
|
| 1365 |
## #cDict[ "cList" ].cd()
|
| 1366 |
## ### find graph with most points:
|
| 1367 |
## #print "-------------", hist, "--------------------"
|
| 1368 |
## indexPoints = 0
|
| 1369 |
## nPoints = 0
|
| 1370 |
## for j in range( 0, len( histDict[ "contourList" + hist ] ) ):
|
| 1371 |
## #print type( histDict[ "contourList" + hist ][ j ] )
|
| 1372 |
## if histDict[ "contourList" + hist ][ j ].GetN() > nPoints:
|
| 1373 |
## nPoints = histDict[ "contourList" + hist ][ j ].GetN()
|
| 1374 |
## indexPoints = j
|
| 1375 |
## pass
|
| 1376 |
## pass
|
| 1377 |
## ### make new graph
|
| 1378 |
## #get the x,y tuple first
|
| 1379 |
## pointTuples = []
|
| 1380 |
## for i in range( 0, histDict[ "contourList" + hist ][ indexPoints ].GetN() ):
|
| 1381 |
## x = ROOT.Double( 0 )
|
| 1382 |
## y = ROOT.Double( 0 )
|
| 1383 |
## histDict[ "contourList" + hist ][ indexPoints ].GetPoint( i, x, y )
|
| 1384 |
## pointTuples.append( ( x,y ) )
|
| 1385 |
## pass
|
| 1386 |
## # now find the index of the list, where y is greatest and x is smallest
|
| 1387 |
## y = -1
|
| 1388 |
## ind = 0
|
| 1389 |
## x = 100000
|
| 1390 |
## for pt in pointTuples:
|
| 1391 |
## if pt[ 1 ] > y:
|
| 1392 |
## #x = pt[ 0 ]
|
| 1393 |
## y = pt[ 1 ]
|
| 1394 |
## pass
|
| 1395 |
## pass
|
| 1396 |
## for pt in pointTuples:
|
| 1397 |
## if pt[ 1 ] < y:
|
| 1398 |
## continue
|
| 1399 |
## if pt[ 0 ] <= x:
|
| 1400 |
## x = pt[ 0 ]
|
| 1401 |
## ind = pointTuples.index( pt )
|
| 1402 |
## pass
|
| 1403 |
## pass
|
| 1404 |
## ### split tuple List in list before ind and after ind:
|
| 1405 |
## list1 = pointTuples[ : ind ]
|
| 1406 |
## list2 = pointTuples[ ind : ]
|
| 1407 |
|
| 1408 |
## #print list2
|
| 1409 |
## #print list1
|
| 1410 |
|
| 1411 |
## ### now fill new graph
|
| 1412 |
## cont = ROOT.TGraph()
|
| 1413 |
## p = 0
|
| 1414 |
## startPointFound = False
|
| 1415 |
## #print "List 2"
|
| 1416 |
|
| 1417 |
## print "=========== Snippet:"
|
| 1418 |
## for pt2 in list2:
|
| 1419 |
## if pt2[ 0 ] >= 700:
|
| 1420 |
## continue
|
| 1421 |
## elif pt2[ 1 ] >= 520:
|
| 1422 |
## continue
|
| 1423 |
## if pt2[ 0 ] <= 410 and pt2[ 1 ] >= 370:
|
| 1424 |
## startPointFound = True
|
| 1425 |
## pass
|
| 1426 |
## if startPointFound:
|
| 1427 |
## print "point added: x=", pt2[ 0 ], " y=", pt2[ 1 ]
|
| 1428 |
## cont.SetPoint( p, pt2[ 0 ], pt2[ 1 ] )
|
| 1429 |
## #print p, ": x=", pt2[ 0 ], "y=", pt2[ 1 ]
|
| 1430 |
## p +=1
|
| 1431 |
## pass
|
| 1432 |
## pass
|
| 1433 |
## #print "List1"
|
| 1434 |
## for pt1 in list1:
|
| 1435 |
## if pt1[ 0 ] >= 700:
|
| 1436 |
## continue
|
| 1437 |
## elif pt1[ 1 ] >= 520:
|
| 1438 |
## continue
|
| 1439 |
## if pt1[ 0 ] <= 410 and pt1[ 1 ] >= 370:
|
| 1440 |
## startPointFound = True
|
| 1441 |
## pass
|
| 1442 |
## if startPointFound:
|
| 1443 |
## print "point added: x=", pt1[ 0 ], " y=", pt1[ 1 ]
|
| 1444 |
## cont.SetPoint( p, pt1[ 0 ], pt1[ 1 ] )
|
| 1445 |
## #print p, ": x=", pt1[ 0 ], "y=", pt1[ 1 ]
|
| 1446 |
## p +=1
|
| 1447 |
## pass
|
| 1448 |
## pass
|
| 1449 |
## cont.SetLineColor( color )
|
| 1450 |
## cont.SetLineWidth( 3 )
|
| 1451 |
## smooth( cont, 15 )
|
| 1452 |
## histDict[ "contSnippet" + hist ] = cont
|
| 1453 |
## #cont.Draw( "AL" )
|
| 1454 |
## pass
|
| 1455 |
|
| 1456 |
## def drawContour( histDict, cDict, hist="ObsEval" ):
|
| 1457 |
## cDict[ "cList" ] = ROOT.TCanvas( "cList" , "cList" )
|
| 1458 |
## cDict[ "cList" ].cd()
|
| 1459 |
## ### find graph with most points:
|
| 1460 |
## indexPoints = 0
|
| 1461 |
## nPoints = 0
|
| 1462 |
## for j in range( 0, len( histDict[ "contourList" + hist ] ) ):
|
| 1463 |
## if histDict[ "contourList" + hist ][ j ].GetN() > nPoints:
|
| 1464 |
## nPoints = histDict[ "contourList" + hist ][ j ].GetN()
|
| 1465 |
## indexPoints = j
|
| 1466 |
## pass
|
| 1467 |
## pass
|
| 1468 |
## ### make new graph
|
| 1469 |
## cont = ROOT.TGraph()
|
| 1470 |
## removePointTuples = []
|
| 1471 |
## cont = histDict[ "contourList" + hist ][ indexPoints ].Clone()
|
| 1472 |
## for p in range( 0, histDict[ "contourList" + hist ][ indexPoints ].GetN() ):
|
| 1473 |
## x = ROOT.Double( 0 )
|
| 1474 |
## y = ROOT.Double( 0 )
|
| 1475 |
## histDict[ "contourList" + hist ][ indexPoints ].GetPoint( p, x, y )
|
| 1476 |
## #print p, ": x=", x, "y=", y
|
| 1477 |
## if x >= 800:
|
| 1478 |
## #counter += 1
|
| 1479 |
## removePointTuples.append( ( x,y ) )
|
| 1480 |
## continue
|
| 1481 |
## elif y >= 600:
|
| 1482 |
## removePointTuples.append( ( x,y ) )
|
| 1483 |
## #counter += 1
|
| 1484 |
## continue
|
| 1485 |
## ##cont.SetPoint( p - counter, x, y )
|
| 1486 |
## #pointTuples.append( ( x, y ) )
|
| 1487 |
## #pass
|
| 1488 |
## #sorted( pointTuples, key=lambda point: point[ 0 ])
|
| 1489 |
## for p in range( 0, histDict[ "contourList" + hist ][ indexPoints ].GetN() ):
|
| 1490 |
## x = ROOT.Double( 0 )
|
| 1491 |
## y = ROOT.Double( 0 )
|
| 1492 |
## histDict[ "contourList" + hist ][ indexPoints ].GetPoint( p, x, y )
|
| 1493 |
## if ( x, y ) in removePointTuples:
|
| 1494 |
## #print p, ": x=", pointTuples[ p ][ 0 ], "y=", pointTuples[ p ][ 1 ]
|
| 1495 |
## ind = findPointIndex( cont, ( x,y ) )
|
| 1496 |
## #print ind
|
| 1497 |
## if not ind == -1:
|
| 1498 |
## cont.RemovePoint( ind )
|
| 1499 |
## pass
|
| 1500 |
## pass
|
| 1501 |
## pass
|
| 1502 |
## cont.SetLineColor( ROOT.kAzure + 7 )
|
| 1503 |
## cont.SetLineWidth( 3 )
|
| 1504 |
## histDict[ "cont" + hist ] = cont
|
| 1505 |
## cont.Draw( "AL" )
|
| 1506 |
## pass
|
| 1507 |
|
