extracting c_beta information from obs mode

[Letizia Stanghellini]

Hi Christophe,

I am using PyNeb with the uncertainty/Montecarlo mode and the observations commands to calculate abundances. It works fine, the program is appended. Would it be possible to extract the cbeta (extinction constant at Hbeta) values in obs. mode? I could not find a hint in the manual, maybe it is obvious. 

Thank you!
Cheers,

Letizia

all_atoms = None  # Initialize variable to store unique atoms
H1 = pn.RecAtom('H', 1)
H1.name
scaler = StandardScaler()
default_T_O3 = np.nan
default_N_O2 = np.nan
obs = pn.Observation(obsFile= f'obs_{galaxy_id}_wholegalaxy.txt', fileFormat='lines_in_rows_err_cols')
obs.extinction.law = 'CCM89'
#obs.extinction.law = 'Cal00'
obs.def_EBV(label1="H1r_4341A", label2="H1r_4861A", r_theo=0.4721)
obs.addMonteCarloObs(N=100)
obs.correctData(normWave=4861.)
O3_4959_line = obs.getLine(label='O3_4959A')
cO3_4959_flux = O3_4959_line.corrIntens    
O3_5007_line = obs.getLine(label='O3_5007A')
cO3_5007_flux = O3_5007_line.corrIntens      
O3_4363_line = obs.getLine(label='O3_4363A')
cO3_4363_flux = O3_4363_line.corrIntens
T_ratio = cO3_4363_flux / cO3_5007_flux
O2_3726_line = obs.getLine(label='O2_3726A')
cO2_3726_flux = O2_3726_line.corrIntens      
O2_3729_line = obs.getLine(label='O2_3729A')
cO2_3729_flux = O2_3729_line.corrIntens
N_ratio = cO2_3726_flux / cO2_3729_flux
H1r_4341A_line = obs.getLine(label='H1r_4341A')
cH1r_4341A_flux=H1r_4341A_line.corrIntens
H1r_4341A_flux=H1r_4341A_line.obsIntens

T_ratio_min = 0
T_ratio_max = 100
N_ratio_min = 0.1
N_ratio_max = 3
diag = pn.Diagnostics()
diag.addDiagsFromObs(obs)
T_O3, N_O2 = diag.getCrossTemDen('[OIII] 4363/5007+', '[OII] 3726/3729', obs=obs, use_ANN=True)
if all_atoms is None:
    all_atoms = obs.getUniqueAtoms()
atomDict = pn.getAtomDict(atom_list=all_atoms)
ion_ab_dic = {}
for line in obs.getSortedLines():
    atom = atomDict.get(line.atom)
    if atom is not None:
        ion_ab_dic[line.label] = atom.getIonAbundance(line.corrIntens, T_O3, N_O2, to_eval=line.to_eval)
Alin_O3 = ion_ab_dic.get('O3_5007A', np.nan)
Alin_O2 = 0.5 * (ion_ab_dic.get('O2_3729A', np.nan)) + 0.5 * (ion_ab_dic.get('O2_3726A', np.nan))
Alin_OH = Alin_O2 + Alin_O3
Tratio = np.nanmedian(T_ratio) if not np.isnan(T_ratio).all() else np.nan
Nratio = np.nanmedian(N_ratio) if not np.isnan(N_ratio).all() else np.nan          
Te = np.nanmedian(T_O3) if not np.isnan(T_O3).all() else np.nan
eTe = np.nanstd(T_O3) if not np.isnan(T_O3).all() else np.nan
Ne = np.nanmedian(N_O2) if not np.isnan(N_O2).all() else np.nan
eNe = np.nanstd(N_O2) if not np.isnan(N_O2).all() else np.nan        
A_O2 = 12 + np.nanmedian(np.log10(Alin_O2)) if not np.isnan(Alin_O2).all() else np.nan
eA_O2 = np.nanstd(np.log10(Alin_O2)) if not np.isnan(Alin_O2).all() else np.nan            
A_O3 = 12 + np.nanmedian(np.log10(Alin_O3)) if not np.isnan(Alin_O3).all() else np.nan
eA_O3 = np.nanstd(np.log10(Alin_O3)) if not np.isnan(Alin_O3).all() else np.nan
A_OH = 12 + np.nanmedian(np.log10(Alin_OH)) if not np.isnan(Alin_OH).all() else np.nan
eA_OH = np.nanstd(np.log10(Alin_OH)) if not np.isnan(Alin_OH).all() else np.nan
eA_OH_plus = (12 + np.log10(np.nanmedian(Alin_OH) + np.nanstd(Alin_OH))) - A_OH
eA_OH_minus = (12 + np.log10(np.nanmedian(Alin_OH) - np.nanstd(Alin_OH))) - A_OH

print(Te,eTe,Ne,eNe,A_OH,eA_OH_plus,eA_OH_minus)

[Christophe Morisset]

Hi Letizia,

Sure, the obs object contains a RedCorr instantiation called extinction, which contains the reddening parameters. So you should be able to obtain cHbeta from:

obs.extinction.cHbeta

Hope it helps,

best saludos,

Ch.

[Letizia Stanghellini]

thank you!  I am still unsure. As you see in my program, I find the Te, Ne, and abundances for 100 fake observations around the flux uncertainties and then I get the median, but I can not do that with ext=obs.extinction.cHbeta, as this gives me just one value (even if I set this after the obs.addMonteCarloObs(N=100) command). Something I am missing. 

[Christophe Morisset]

Oh, yes, you are defining the reddening correction before adding the MC fake observations. So you have only one cHbeta. 

You could just move the obs.def_EBV(label1="H1r_4341A", label2="H1r_4861A", r_theo=0.4721) line after the obs.addMonteCarloObs(N=100)one and it should work.

Ch.

[Letizia Stanghellini]

thanks! I just could not see it without your help. Have  a great day