ComboCode.cc.statistics.UnresoStats

1 # -*- coding: utf-8 -*- 2 3 """ 4 Performing simple peak-to-peak statistics. 5 6 Author: R. Lombaert 7 8 """ 9 10 import os 11 import scipy 12 from scipy import argmin,array,sqrt,log10 13 import numpy as np 14 import operator 15 16 import cc.path 17 from cc.tools.io import DataIO 18 from cc.modeling.objects import Transition 19 from cc.statistics.Statistics import Statistics 20 from cc.statistics import BasicStats as bs 21 from cc.plotting import Plotting2 22 23 24

25 -class UnresoStats(Statistics):

26 27 """ 28 Environment with several tools to perform statistics on peak-to-peak ratios 29 for unresolved lines. 30 31 """ 32

33 - def __init__(self,star_name,code='GASTRoNOoM',path_code='codeSep2010'):

34 35 """ 36 Initializing an instance of UnresoStats. 37 38 @param star_name: Star name from Star.dat 39 @type star_name: string 40 41 @keyword code: the code used for producing your output 42 43 (default: 'GASTRoNOoM') 44 @type code: string 45 @keyword path_code: Output folder in the code's home folder 46 47 (default: 'codeSep2010') 48 @type path_code: string 49 50 """ 51 52 super(UnresoStats,self).__init__(star_name=star_name,\ 53 code=code,path_code=path_code) 54 55 #-- Remember sample transitions per band and their central wavelengths. 56 # key: filename, value: list[Transition()] 57 self.sample_trans = dict() 58 self.central_mwav = dict() 59 60 #-- Remember the peak and integrated flux ratios, and the data errors on 61 # the ratios 62 # key: filename 63 # value: dict([instrument based id, list[float] or float]) 64 self.peak_ratios = dict() 65 self.int_ratios = dict() 66 self.int_ratios_err = dict() 67 68 #-- Remember the line strengths from observations and models per band 69 self.dint_bands = dict() 70 self.derr_bands = dict() 71 self.blends_bands = dict() 72 self.mint_bands = dict() 73 74 #-- Remember chi2 per model for integrated line strengths and straight 75 # up comparison between convolved model and observed spectrum. 76 # key: instrument based id 77 # value: float 78 self.chi2_inttot = dict() 79 self.chi2_con = dict()

80 81 82

83 - def setInstrument(self,instrument_name,*args,**kwargs):

84 85 ''' 86 Set an instrument, see Statistics.py. 87 88 @param instrument_name: The instrument (such as 'PACS', 'SPIRE') 89 @type instrument_name: string 90 91 ''' 92 93 if instrument_name == 'FREQ_RESO': 94 raise IOError('Unresolved data stats cannot be calculated for FREQ_RESO.') 95 super(UnresoStats,self).setInstrument(instrument_name=instrument_name,\ 96 *args,**kwargs)

97 98 99

100 - def setLineStrengths(self):

101 102 ''' 103 Find the peak to peak ratios of data versus model. 104 105 The result are saved in the self.peak_ratios dictionary, see 106 __init__.__doc__() 107 108 ''' 109 110 inst = self.instrument 111 print '***********************************' 112 print '** Gathering model and observed line strengths for %s.'\ 113 %inst.instrument 114 115 #-- Make a sample selection of Transition()s. 116 sample_trans = Transition.extractTransFromStars(self.star_grid,\ 117 dtype=inst.instrument) 118 119 #-- Set the tolerance: a factor that is multiplied with half the 120 # wavelength resolution of data when looking for the peak value 121 # around a central wavelength in the data. As a default, this is 122 # equal to the PACS_OVERSAMPLING. No point in changing this. 123 self.tolerance = inst.oversampling 124 125 for ifn,(fn,dwav) in enumerate(zip(inst.data_filenames,\ 126 inst.data_wave_list)): 127 #-- Create a list of sample transitions 128 self.sample_trans[fn] = [trans 129 for trans in sample_trans 130 if trans.wavelength*10**4 >= dwav[0]\ 131 and trans.wavelength*10**4 <= dwav[-2]] 132 #-- Get the central wavelength of the lines, corrected for 133 # Doppler shift due to vlsr of the central source. In micron. 134 self.central_mwav[fn] = [t.wavelength*10**4*1./(1-inst.vlsr/t.c) 135 for t in self.sample_trans[fn]] 136 137 self.__setPeakRatios(ifn,fn) 138 if not inst.linefit is None: 139 self.__setIntRatios(ifn,fn) 140 141 print '***********************************'

142 143 144

145 - def calcChiSquared(self,chi2_method='diff',ndf=0,filename=None):

146 147 ''' 148 Calculate the chi_squared value for given models and data. 149 150 If integrated fluxes are available, they are used without the line 151 blends. LINES FLAGGED AS A BLEND ARE EXCLUDED FROM THE CHI2 CALCULATION. 152 153 If not, the modeled, convolved spectrum is compared directly with the 154 data, where the model is not just 0 flux. 155 156 @keyword chi2_method: The type of chi-squared calculated for integrated 157 fluxes. 'diff' for standard chi^2 kind, 'log' for 158 redistributing the data/model ratios on an 159 absolute logarithmic scale before calculating the 160 chi2 161 162 (default: 'diff') 163 @type chi2_method: string 164 @keyword ndf: Number of degrees of freedom. Default in case of calculating 165 for one single model. Typically the number of variable grid 166 parameters in a grid calculation. 167 168 (default: 0) 169 @type ndf: int 170 @keyword filename: the filename for which you want to return the list, 171 if None, all filenames are used and the lists are 172 merged into one 173 174 (default: None) 175 @type filename: string 176 ''' 177 178 if chi2_method not in ['diff','log']: 179 chi2_method = 'diff' 180 print "WARNING: STAT_CHI2 not recognized. Using default 'diff'." 181 182 #-- For now excluding blends 183 inst = self.instrument 184 all_dints = self.getRatios(sel_type='dint_bands',data_type='dint_bands',\ 185 filename=filename) 186 all_derrs = self.getRatios(sel_type='dint_bands',data_type='derr_bands',\ 187 filename=filename) 188 189 for istar,star in enumerate(self.star_grid): 190 this_id = star['LAST_%s_MODEL'%inst.instrument.upper()] 191 192 #-- Get all integrated fluxes and errors for this particular model, 193 # no matter the filename of the data (ie no matter the band). 194 # If no integrated flux available, set chi2_inttot[this_id] to 0 195 all_mints = self.getRatios(this_id=this_id,sel_type='dint_bands',\ 196 data_type='mint_bands',filename=filename) 197 sel = np.isfinite(all_mints)*np.isfinite(all_dints) 198 199 if not all_mints[sel].size: 200 self.chi2_inttot[this_id] = 0 201 else: 202 chi2 = bs.calcChiSquared(data=all_dints,model=all_mints,\ 203 noise=all_derrs*all_dints,\ 204 mode=chi2_method,ndf=ndf) 205 self.chi2_inttot[this_id] = chi2 206 207 #-- Calculate chi2 based on the convolved model with respect to the 208 # continuum-subtracted data. 209 all_dflux = [] 210 all_mflux = [] 211 all_dstd = [] 212 if not filename: 213 fns = inst.data_filenames 214 fluxes = inst.data_flux_list 215 else: 216 fns = filename is None and inst.data_filenames or [filename] 217 fluxes = [inst.data_flux_list[inst.data_filenames.index(fn)] 218 for fn in fns] 219 220 for fn,dflux in zip(fns,fluxes): 221 dstd = self.data_stats[fn]['std'] 222 #-- Cannot return empty list as the selection of existing 223 # convolutions is done in Statistics.setModels() 224 mflux = inst.getSphinxConvolution(star,fn)[1] 225 #-- Make sure all points are positive. Some data points may be 226 # <0 due to continuum subtraction. 227 mflux = mflux[dflux>0] 228 dflux = dflux[dflux>0] 229 all_dflux.extend(dflux[mflux>0]) 230 all_mflux.extend(mflux[mflux>0]) 231 all_dstd.extend([dstd]*len(mflux[mflux>0])) 232 self.chi2_con[this_id] = bs.calcChiSquared(all_dflux,\ 233 all_mflux,all_dstd,\ 234 mode=chi2_method,ndf=ndf)

235 236

237 - def __setIntRatios(self,ifn,fn):

238 239 ''' 240 Calculate ratios of integrated intensities, if requested. 241 242 Assumes there is a linefit filename available in the Instrument() 243 object. Only done for those line present in this file, based on 244 Doppler shifted wavelength. 245 246 @param ifn: Index of the data band in self.instrument lists. 247 @type ifn: int 248 @param fn: The filename of the data set. Needed for book keeping. 249 @type fn: string 250 251 ''' 252 253 #-- Get some data properties, and extract data wavelength and flux 254 inst = self.instrument 255 self.dint_bands[fn] = [] 256 self.derr_bands[fn] = [] 257 self.blends_bands[fn] = [] 258 self.mint_bands[fn] = dict() 259 self.int_ratios[fn] = dict() 260 self.int_ratios_err[fn] = dict() 261 262 #-- Comparing integrated intensities between PACS and models. 263 # Comparisons only made per filename! 264 inst.intIntMatch(trans_list=self.sample_trans[fn],ifn=ifn) 265 266 for st in self.sample_trans[fn]: 267 dintint, dintinterr, blends = st.getIntIntUnresolved(fn) 268 if dintint is None or dintint == 'inblend': 269 self.dint_bands[fn].append(np.nan) 270 self.derr_bands[fn].append(np.nan) 271 self.blends_bands[fn].append(None) 272 else: 273 if not blends is None: 274 #-- Take abs in case dintint is already blend due to FWHM 275 dintint = -1*abs(dintint) 276 self.dint_bands[fn].append(dintint) 277 self.derr_bands[fn].append(dintinterr) 278 self.blends_bands[fn].append(blends) 279 280 self.dint_bands[fn] = array(self.dint_bands[fn]) 281 self.derr_bands[fn] = array(self.derr_bands[fn]) 282 283 for star in self.star_grid: 284 #-- From here on, we start extracting the model specific int ints. 285 this_id = star['LAST_%s_MODEL'%inst.instrument.upper()] 286 self.mint_bands[fn][this_id] = [] 287 mtrans = array([star.getTransition(t) 288 for t in self.sample_trans[fn]]) 289 290 for mt,dint,blend in zip(mtrans,self.dint_bands[fn],\ 291 self.blends_bands[fn]): 292 # 4) No trans == sample_trans found for this model, 293 # Note that model value is added even if PACS int int is 294 # not available 295 if mt is None: 296 self.mint_bands[fn][this_id].append(np.nan) 297 298 # 5) Match found with a wave_fit value. Get int ratio 299 # m/d. If dint is negative, it is a blend due to large 300 # FWHM! If blend is not None, multiple sample trans have 301 # been found in the wavelength resolution bin of the 302 # fitted line and also indicates a blend. 303 else: 304 if blend is None: 305 mintint = mt.getIntIntIntSphinx() 306 self.mint_bands[fn][this_id].append(mintint) 307 else: 308 #-- blend is a list of sample transitions that refers 309 # to the transitions involved in the blend, so get 310 # these from the model grid, and add them up. dintint 311 # was made negative to indicate blend presence. 312 blendlines = [star.getTransition(t) 313 for t in blend 314 if not star.getTransition(t) is None] 315 total_mint = sum([t.getIntIntIntSphinx() 316 for t in blendlines]) 317 self.mint_bands[fn][this_id].append(total_mint) 318 319 #-- Note that the error is relative, and the chi^2 wants 320 # an absolute value. 321 self.mint_bands[fn][this_id] = array(self.mint_bands[fn][this_id]) 322 self.int_ratios[fn][this_id] = self.mint_bands[fn][this_id]\ 323 /self.dint_bands[fn] 324 self.int_ratios_err[fn][this_id] = self.derr_bands[fn]*\ 325 abs(self.int_ratios[fn][this_id])

326 327 328

329 - def __setPeakRatios(self,ifn,fn):

330 331 ''' 332 Calculate Peak ratios for all models included in the star grid. 333 334 Done per filename. 335 336 @param ifn: Index of the data band in self.instrument lists. 337 @type ifn: int 338 @param fn: The filename of the data set. Needed for book keeping. 339 @type fn: string 340 341 ''' 342 343 #-- Get some data properties, and extract data wavelength and flux 344 inst = self.instrument 345 dwav = inst.data_wave_list[ifn] 346 dflux = inst.data_flux_list[ifn] 347 348 #-- Get some data statistics 349 d_mean = self.data_stats[fn]['mean'] 350 d_std = self.data_stats[fn]['std'] 351 #-- this sigma is used for determining whether a peak flux is 352 # significant compared to the std value of the band. 353 d_sigma = self.data_stats[fn]['sigma'] 354 355 self.peak_ratios[fn] = dict() 356 357 for star in self.star_grid: 358 #-- Read the convolved sphinx model 359 mwav, mflux = inst.getSphinxConvolution(star,fn) 360 if list(mflux[mflux < 0]) != []: 361 print 'There are negative sphinx flux values! They will '+\ 362 'not be taken into account.' 363 364 #-- Calculate the peak-to-peak ratios. 365 # 1) Central wavelengths of mtrans are set in previous method 366 # 2) Get the central model flux, which should coincide exactly 367 # with the Doppler shifted rest wavelength of the line. If the 368 # model flux is negative, the value is not used. 369 central_mflux = [mflux[argmin(abs(mwav-wav))] > 0 \ 370 and mflux[argmin(abs(mwav-wav))] \ 371 or None 372 for wav in self.central_mwav[fn]] 373 # 3) Get the central data flux, at the Doppler shifted central 374 # wavelength expected from the model. The maximum flux is 375 # taken in the wavelength bin tolerance*wav_resolution/2. 376 # allowing for small wave shifts due to instrumental effects. 377 central_dflux = [max(dflux[abs(dwav-wav)<= \ 378 self.tolerance/2.*( dwav[argmin(abs(dwav-wav))+1]\ 379 -dwav[argmin(abs(dwav-wav))])]) 380 for wav in self.central_mwav[fn]] 381 # 4) Check if the data flux point is actually significant 382 # compared to the noise in the spectrum. Compare with dstd, 383 # given d_sigma from path_combocode/usr/Data.dat . 384 # Insignificant values are multiplied by -1, to indicate they 385 # are upper limits in the data at that wavelength. 386 central_dflux = [d >= d_mean+(d_std*d_sigma) \ 387 and d \ 388 or -1*abs(d_mean+(d_std*d_sigma)) 389 for d in central_dflux] 390 # 5) Calculate the ratios, only if the model flux is not None 391 # (was a negative model flux value: We don't want that) 392 # Negative ratios are possible, in case of ratio lower limits 393 this_id = star['LAST_%s_MODEL'%inst.instrument.upper()] 394 self.peak_ratios[fn][this_id] = [not m is None and m/d or None 395 for m,d in zip(central_mflux,\ 396 central_dflux)]

397 398

399 - def getRatios(self,this_id=None,sel_type='peak_ratios',\ 400 data_type='peak_ratios',return_negative=0,filename=None):

401 402 ''' 403 Return all ratios for all filenames, including only the values that are 404 not None. 405 406 Selection type helps to filter out blended lines (negative values for 407 the ratios), or unavailable ratios (None). Can only be int_ratios or 408 peak_ratios when this_id is not None. 409 410 The data type that is returned can be any of the dictionaries tracking 411 stuff in this object. 412 413 @keyword this_id: The requested instrument id. Can be None in case both 414 sel_type and data_type point to dictionaries that are 415 not sorted per model. 416 417 (default: None) 418 @type this_id: string 419 @keyword sel_type: The type of data used to select values, one of 420 'peak_ratios' or 'int_ratios' 421 422 (default: 'peak_ratios') 423 @type sel_type: string 424 @keyword data_type: the type of data returned, one of 'peak_ratios', 425 'int_ratios', 'central_mwav', '*_bands' 426 427 (default: 'peak_ratios') 428 @type data_type: string 429 @keyword return_negative: only return the negative peak ratios or 430 equivalent (lower limits). 431 432 (default: 0) 433 @type return_negative: bool 434 @keyword filename: the filename for which you want to return the list, 435 if None, all filenames are used and the lists are 436 merged into one 437 438 (default: None) 439 @type filename: string 440 441 @return: The ratios requested 442 @rtype: array 443 444 ''' 445 446 inst = self.instrument 447 filenames = filename is None and inst.data_filenames or [filename] 448 449 bandsel = ['dint_bands','derr_bands','blends_bands','central_mwav',\ 450 'sample_trans'] 451 modelsel = ['int_ratios','peak_ratios','mint_bands','int_ratios_err'] 452 453 if sel_type in bandsel and data_type in bandsel: 454 this_id = None 455 456 #-- Cannot select values for model specific dicts when id not given 457 if this_id is None and (sel_type in modelsel or data_type in modelsel): 458 return None 459 460 elif this_id is None: 461 return array([v2 462 for fn in filenames 463 for v,v2 in zip(getattr(self,sel_type)[fn],\ 464 getattr(self,data_type)[fn]) 465 if not v is None \ 466 and ((return_negative and v < 0) \ 467 or (not return_negative and v > 0))]) 468 469 elif data_type in bandsel: 470 return array([v2 471 for fn in filenames 472 for v,v2 in zip(getattr(self,sel_type)[fn][this_id],\ 473 getattr(self,data_type)[fn]) 474 if not v is None \ 475 and ((return_negative and v < 0) \ 476 or (not return_negative and v > 0))]) 477 478 elif sel_type in bandsel: 479 return array([v2 480 for fn in filenames 481 for v,v2 in zip(getattr(self,sel_type)[fn],\ 482 getattr(self,data_type)[fn][this_id]) 483 if not v is None \ 484 and ((return_negative and v < 0) \ 485 or (not return_negative and v > 0))]) 486 487 else: 488 return array([v2 489 for fn in filenames 490 for v,v2 in zip(getattr(self,sel_type)[fn][this_id],\ 491 getattr(self,data_type)[fn][this_id]) 492 if not v is None \ 493 and ((return_negative and v < 0) \ 494 or (not return_negative and v > 0))])

495 496 497

498 - def plotRatioWav(self,inputfilename,no_peak=False):

499 500 ''' 501 Plot ratios as a function of their central wavelength. 502 503 In blue and green, the peak-to-peak ratios are plotted. Blue for normal 504 peak-to-peak ratios, green for ratios that involve a data point that is 505 in the noise. Plotting peak-to-peak ratios can be turned off with the 506 no_peak keyword. 507 508 In red and magenta, the integrated line strengths are plotted. Red for 509 isolated lines that are not tagged as a blend, magenta for lines that 510 are tagged as a blend due to 1) too wide lines with respect to the PACS 511 resolution, or 2) multiple sample transitions in the width of an 512 observed line. 513 514 @param inputfilename: the input filename for the grid, which will be 515 attached to the final plot filename 516 @type inputfilename: string 517 518 @keyword no_peak: Hide the peak ratios. 519 520 (default: False) 521 @type no_peak: bool 522 523 ''' 524 525 if not self.chi2_inttot: 526 print "No Sphinx models calculated. Aborting statistics plot. " 527 return 528 this_grid = self.sortStarGrid() 529 plot_filenames = [] 530 inst = self.instrument 531 for star in this_grid: 532 this_id = star['LAST_%s_MODEL'%inst.instrument.upper()] 533 lp = [] 534 waves = [] 535 ratios = [] 536 ratios_err = [] 537 538 #-- the ratios included in the statistics 539 if not no_peak: 540 this_wav_inc = self.getRatios(data_type='central_mwav',\ 541 this_id=this_id) 542 this_ratio_inc = self.getRatios(this_id=this_id) 543 if list(this_wav_inc): 544 waves.append(this_wav_inc) 545 ratios.append(this_ratio_inc) 546 ratios_err.append(None) 547 lp.append('ob') 548 549 #-- ratios replaced by lower limits if data point is in the noise 550 # ie data point can only be smaller than or equal to used value 551 this_wav_lower = self.getRatios(data_type='central_mwav',\ 552 this_id=this_id,\ 553 return_negative=1) 554 this_ratio_lower = self.getRatios(return_negative=1,\ 555 this_id=this_id) 556 this_ratio_lower = [abs(r) for r in this_ratio_lower] 557 if list(this_wav_lower): 558 waves.append(this_wav_lower) 559 ratios.append(this_ratio_lower) 560 ratios_err.append(None) 561 lp.append('dg') 562 563 if not inst.linefit is None: 564 #-- If integrated intensities are available for the instrument, get 565 # the integrated intensity ratios 566 this_wav_int = self.getRatios(sel_type='int_ratios',\ 567 data_type='central_mwav',\ 568 this_id=this_id) 569 this_ratio_int = self.getRatios(sel_type='int_ratios',\ 570 data_type='int_ratios',\ 571 this_id=this_id) 572 this_ratio_int_err = self.getRatios(sel_type='int_ratios',\ 573 data_type='int_ratios_err',\ 574 this_id=this_id) 575 if list(this_wav_int): 576 waves.append(this_wav_int) 577 ratios.append(this_ratio_int) 578 ratios_err.append(this_ratio_int_err) 579 lp.append('or') 580 581 #-- Get the ratios that are lower limits due to line blends. 582 # Line blends detected due to fitted FHWM/PACS FHWM > 120% 583 # ie model int can only be larger than or equal to used value 584 this_wav_lowerint = self.getRatios(sel_type='int_ratios',\ 585 data_type='central_mwav',\ 586 this_id=this_id,\ 587 return_negative=1) 588 this_ratio_lowerint = self.getRatios(sel_type='int_ratios',\ 589 data_type='int_ratios',\ 590 this_id=this_id,\ 591 return_negative=1) 592 this_ratio_lowerint_err = self.getRatios(sel_type='int_ratios',\ 593 data_type='int_ratios_err',\ 594 this_id=this_id,\ 595 return_negative=1) 596 this_ratio_lowerint = [abs(r) for r in this_ratio_lowerint] 597 if list(this_wav_lowerint): 598 waves.append(this_wav_lowerint) 599 ratios.append(this_ratio_lowerint) 600 ratios_err.append(this_ratio_lowerint_err) 601 lp.append('dm') 602 603 #- prepping input for the plot command 604 xmin = min([min(x) for x in waves]) 605 xmax = max([max(x) for x in waves]) 606 waves.extend([[0.5*xmin,1.5*xmax]]*3) 607 ratios.extend([[1,1],\ 608 [1-inst.absflux_err,\ 609 1-inst.absflux_err],\ 610 [1+inst.absflux_err,\ 611 1+inst.absflux_err]]) 612 ratios_err.extend([None,None,None]) 613 lp.extend(['-k','--k','--k']) 614 plot_filename = os.path.join(getattr(cc.path,self.code.lower()),\ 615 self.path_code,'stars',\ 616 self.star_name,\ 617 '%s_results_'%inst.instrument+\ 618 'ratio_wav_%s'%str(this_id)) 619 labels = [('Mdot = %.2e Msolar/yr'%star['MDOT_GAS'],0.05,0.05),\ 620 ('Teff = %.1f K'%star['T_STAR'],0.05,0.1),\ 621 ('$\psi$ = %0.2e'%star['DUST_TO_GAS_CHANGE_ML_SP'],0.05,\ 622 0.15),\ 623 ('A$_{H_2O}$/A$_{H_2}$ = %0.2e'%star['F_H2O'],0.05,0.2),\ 624 ('R$_(o,H_2O)$ = %i'%int(star['R_OUTER_H2O']),0.05,0.25)] 625 if star.getMolecule('1H1H16O') \ 626 and star.getMolecule('1H1H16O').set_keyword_change_abundance: 627 labels.append(('$H_2O$ profile = %s'\ 628 %os.path.split(star.getMolecule('1H1H16O')\ 629 .change_fraction_filename\ 630 .replace('_','\_'))[1],\ 631 0.05,0.30)) 632 plot_title = '%s: $\chi^2_\mathrm{con}$ %.4f'\ 633 %(str(this_id).replace('_','\_'),\ 634 self.chi2_con[this_id]) 635 if self.chi2_inttot[this_id]: 636 plot_title += ', $\chi^2_\mathrm{int}$ %.4f'\ 637 %(self.chi2_inttot[this_id]) 638 plot_filenames.append(Plotting2.plotCols(\ 639 filename=plot_filename,x=waves,y=ratios,yerr=ratios_err,\ 640 yaxis=r'$F_{\nu,p,m}/F_{\nu,p,d}$',\ 641 plot_title=plot_title,labels=labels,extension='pdf',\ 642 xlogscale=0,ylogscale=1,line_types=lp,xmin=xmin*0.9,\ 643 xmax=xmax*1.03,figsize=(10.*scipy.sqrt(2.), 10.),\ 644 linewidth=2,fontsize_title=20,fontsize_label=16)) 645 inputf_short = os.path.splitext(os.path.split(inputfilename)[1])[0] 646 new_filename = os.path.join(getattr(cc.path,self.code.lower()),\ 647 self.path_code,'stars',self.star_name,\ 648 '%s_results_'%inst.instrument+\ 649 'ratio_wav_%s.pdf'%inputf_short) 650 DataIO.joinPdf(old=plot_filenames,new=new_filename) 651 print '** Stat plots can be found at:' 652 print new_filename 653 print '***********************************'

654 655 656

657 - def sortStarGrid(self):

658 659 ''' 660 Return a sorted list of the star grid in this instance according to 661 the chi^2 value. If a chi^2 based on integrated line fluxes is 662 available, it is used. Otherwise the chi^2 determined from convolved 663 model vs data is taken. 664 665 @return: A sorted list of models according to chi^2 666 @rtype: list[Star] 667 668 ''' 669 670 if self.chi2_inttot.values()[0]: 671 styp = 'chi2_inttot' 672 else: 673 styp = 'chi2_con' 674 ikey = 'LAST_%s_MODEL'%self.instrument.instrument.upper() 675 return sorted(self.star_grid,key=lambda x: getattr(self,styp)[x[ikey]])

676

Source Code for Module ComboCode.cc.statistics.UnresoStats