#!/usr/bin/perl
# db_create: converts decay and gamma transitions data from ENSDF
# database to binary form of interlinked hashes.
# TODO: draw a database scheme.

use strict;
use Storable;
my $ensdf;
for my $v qw(/Volumes/Archiv/ensdf /home/hatta/work/ensdf ensdf) {
  if (-d $v) {
    $ensdf = $v;
    last;
  }
}

# convert energy error from ensdf format to a real number
sub decipher_de {
  my ($e, $de) = @_;
  my $e1 = $e;
  $de = 0.5 unless ($de =~ /\d/);
  return $de unless ($e =~ /\./);
  $e =~ s/(\d+)\.(\d+)/$2/;
  my $n_acc = length($e);
  my $r = $de * 10**(-$n_acc);
  return $r;
}

# filter out garbage from number values
sub numfilter {
  my $s = $_[0];
  if ($s =~ /([\d\.]+)(E[\+-]?\d+)?/i) {
#    print "notice:numfilter: '$s' => '$1$2'.\n";
    return "$1$2";
  } else {
    print "error:invalid number value $s.\n" if $s =~ /\S/;
  }
}

# convert ENSDF time field to seconds
sub atohl {
  my $s = shift;
  unless ($s =~ /^\s*([\d\.]+E?\+?-?\d?\d?)\s*(Y|D|H|M|S|MS|US|NS|PS|FS|AS|EV|KEV|MEV)\s*\??\s*$/i) {
    return;
  }
  my $num = $1;
  my $unit = $2;
  my %conversion_time = (Y => 31536000, D => 86400, H => 3600, M => 60,
			 S => 1, MS => 1e-3, US => 1e-6, NS => 1e-9,
			 PS => 1e-12, FS => 1e-15, AS => 1e-18);
  my %conversion_energy = (EV => 1, KEV => 1000, MEV => 1000000);
  if (exists $conversion_time{uc($unit)}) {
    return $conversion_time{uc($unit)} * $num;
  } elsif (exists $conversion_energy{uc($unit)}) {
    return 4.13566733e-15 / ($num * $conversion_energy{uc($unit)});
  } else {
    return;
  }
}

unless(defined($ensdf)) {
  print "error: ENSDF not found.\n";
  exit(1);
}
print "$0 started.\n";
print '$Id: db_create 57 2009-12-28 10:17:54Z hatta $', "\n";
print "reading ENSDF source from $ensdf.\n";
my %ztonam;
my %namtoz;
open(ELEMENTS, "< elements.txt");
while (<ELEMENTS>) {
  if (/(\d+)\s*(\w+)/) {
    $ztonam{$1} = lc($2);
    $namtoz{lc($2)} = $1;
  }
}
close(ELEMENTS);

opendir(ENSDFDIR, $ensdf) || die "can't open ensdf\n";
my @ensdf_files = readdir(ENSDFDIR);
closedir(ENSDFDIR);
my %decay_table;

# look for decay data files and build lookup table
for my $file (@ensdf_files) {
  open(FILE, "< $ensdf/$file");
  $_=<FILE>;			# read first line
  /(\S+)\s*DECAY/ || goto cleanup;	# should contain "DECAY" substring
  my $dtype1 = $1;
  my $dtype = "unknown";
  if (/(\S+)\s*DECAY/) {
    $dtype = $1;
  }
  if (/\sB[\+-]\s+DECAY/ || /\sEC\s+DECAY/) {
#    $dtype = "beta";
  } elsif (/\sA\s+DECAY/) {
#    $dtype = "alpha";
  } elsif (/\sIT\s+DECAY/) {
    $dtype = "isomeric";
  } elsif (/\sSF\s+DECAY/) {
#    print "fission decays are not processed\n";
    goto cleanup;
  } else {
    /(\S+)\s*DECAY/;
    $dtype = $1;
#    print "unknown decay mode $1!\n";
  }
  $_ =~ /(\(.*\))/;
  my $header_halflife = $1;
  (my $final_nucl = substr($_, 0, 5)) =~ s/\s//g;
  $final_nucl =~ /^\s*(\d+)(\w+)\s*$/;
  my $final_a = $1;
  my $final_z = $namtoz{lc($2)};
  /(\d+)(\w+)\s+\S+\s+DECAY/;
  my $parent_nucl = "$1$2";
  my $parent_a = $1;
  my $parent_z = $namtoz{lc($2)};
  my %p_lines;			# parent levels 
  my @l_lines;			# final levels
  my $current_l_line;
  my %n_lines;			# normalizations
  my %b_lines;			# beta emissions
  my %e_lines;			# e+ / EC
  my %g_lines;			# photon emissions
  my %a_lines;			# alpha emissions
  my %particle_lines;		# other immediate particle emissions
  while (<FILE>) {
    if (/^.....  P( |\d)/) {
      $p_lines{$1} = $_;
    } elsif (/^.....  N( |\d)/) {
      $n_lines{$1} = $_;
    } elsif (/^.....  L /) {
      $current_l_line = $_;
      push @l_lines, $_;
    } elsif (/^.....  B /) {
      my $ib = numfilter(substr($_, 21, 8));
      (my $un = substr($_, 77, 2)) =~ s/^\s*(\S+)\s*$/$1/;
      if (defined($current_l_line)) {
	if (exists $b_lines{$current_l_line}) {
	  print "warning:$file: skipping multiple betas for one final level (-).\n";
	  next;
	}
	$b_lines{$current_l_line} = {ib => $ib, un => $un};
      }
    } elsif (/^.....  E /) {
      my $ib = numfilter(substr($_, 21, 8));
      my $ie = numfilter(substr($_, 31, 8));
      my $ti = numfilter(substr($_, 64, 10));
      (my $un = substr($_, 77, 2)) =~ s/^\s*(\S+)\s*$/$1/;
      if (defined($current_l_line)) {
	if (exists $e_lines{$current_l_line}) {
	  print "warning:$file: skip multiple betas for one final level (+).\n";
	  next;
	}
	$e_lines{$current_l_line} = {ib => $ib, un => $un, ie => $ie, ti => $ti};
      }
    } elsif (/^.....  A /) {
      my $ia = numfilter(substr($_, 21, 8));
      if (defined($current_l_line)) {
	if (exists $a_lines{$current_l_line}) {
	  print "warning:$file: skip multiple alphas for one final level.\n";
	  next;
	}
	$a_lines{$current_l_line} = {ia => $ia};
      }
    } elsif (/^.....   (P|N|A)/) {
      my $particle = $1;
      my $ip = numfilter(substr($_, 21, 8));
      if (defined($current_l_line)) {
	if (exists $particle_lines{$current_l_line}) {
	  print "warning:$file: skip multiple particles ($particle) for one final level.\n";
	  next;
	}
	$particle_lines{$current_l_line} = {ip => $ip};
      }
    } elsif (/^.....  G /) {
      if (defined($current_l_line) && exists($g_lines{$current_l_line})) {
	push @{$g_lines{$current_l_line}}, $_;
      } else {
	$g_lines{$current_l_line} = [$_];
      }
    }
  }
  unless (%p_lines) {
    print "warning:$file: no parent record, add a one.\n";
    $p_lines{" "} = "       P 0.0                            $header_halflife                        ";
  }
  print "warning:$file: N problem\n"
    if (keys(%p_lines) > 1 && keys(%n_lines) != keys(%p_lines));
  for my $p_record (values %p_lines) {
    my $p_energy = numfilter(substr($p_record, 9, 10));
    my $q_value = numfilter(substr($p_record, 64, 10));
    (my $spin_parity = substr($p_record, 21, 18)) =~ s/^\s*(\S+)\s*$/$1/;
    (my $halflife = substr($p_record, 39, 10)) =~ s/^\s*(\S+)\s*$/$1/;
    $halflife = "*" unless ($halflife =~ /\d/);
    my $n_record;
    foreach (values %n_lines) {
      $n_record = $_ if(substr($_, 8, 1) eq substr($p_record, 8, 1));
    }
    unless ($n_record) {
      if (scalar(values(%n_lines)) == 1) {
	$n_record = (values(%n_lines))[0];
      }
    }
    my ($nr, $nt, $br, $nb);
    if ($n_record) {
      $nr = numfilter(substr($n_record, 9, 10));
      $nt = numfilter(substr($n_record, 21, 8));
      $br = numfilter(substr($n_record, 31, 8));
      $nb = numfilter(substr($n_record, 41, 8));
    }
    $nr = 1 unless $nr > 0;
    $nt = 1 unless $nt > 0;
    # $br = 1 unless $br > 0;
    # $nb = 0.01 unless ($br > 0 || $nb > 0);
    if ($dtype1 =~ /(B\+|B-|EC)/ && !($nb > 0) && $br > 0) {
      print "warning:$file: beta decay and no NB information, restoring it manually.\n";
      my $sum_ib = 0;
      if ($dtype1 =~ /(B\+)/) {
	foreach (keys %e_lines) {
	  my $iii = $e_lines{$_}->{ib} > 0 ? $e_lines{$_}->{ib} : $e_lines{$_}->{ti};
	  $sum_ib += $iii;
	}
      } elsif ($dtype1 =~ /EC/) {
	foreach (keys %e_lines) {
	  my $iii = $e_lines{$_}->{ie} > 0 ? $e_lines{$_}->{ie} : $e_lines{$_}->{ti};
	  $sum_ib += $iii;
	}
      } elsif ($dtype1 =~ /B-/) {
	$sum_ib += $b_lines{$_}->{ib} foreach keys %b_lines;
      }
      if ($sum_ib > 0) {
	$nb = 100 / $sum_ib;
      } else {
	$nb = 1;
      }
    }
    my @final_levels;
    my %final_level_intensities;
    foreach (@l_lines) {
      my $f_energy = numfilter(substr($_, 9, 10));
      if ($dtype1 =~ /B-/) {
	if (exists $b_lines{$_}) {
	  push @final_levels, $f_energy;
	  $final_level_intensities{$f_energy} = $b_lines{$_}->{ib} * $nb / 100;
	}
      } elsif ($dtype1 =~ /B\+/) {
	if (exists $e_lines{$_}) {
	  push @final_levels, $f_energy;
	  my $iii = $e_lines{$_}->{ib} > 0 ? $e_lines{$_}->{ib} : $e_lines{$_}->{ti};
	  $final_level_intensities{$f_energy} = $iii * $nb / 100;
	}
      } elsif ($dtype1 =~ /EC/) {
	if (exists $e_lines{$_}) {
	  push @final_levels, $f_energy;
	  my $iii = $e_lines{$_}->{ie} > 0 ? $e_lines{$_}->{ie} : $e_lines{$_}->{ti};
	  $final_level_intensities{$f_energy} = $iii * $nb / 100;
	}
      } elsif ($dtype1 =~ /A/) {
	if (exists $a_lines{$_}) {
	  push @final_levels, $f_energy;
	  $final_level_intensities{$f_energy} = $a_lines{$_}->{ia} / 100;
	}
      } else {
	push @final_levels, $f_energy;
	$final_level_intensities{$f_energy} = $particle_lines{$_}->{ip};
      }
    }
    # make sure that there is at least one final level (g.s.)
    if (scalar(@final_levels) == 0) {
      print "warning:$file: add fake 0.0 final level.\n";
      push @final_levels, "0.0";
      $final_level_intensities{"0.0"} = 1;
    }
    # XXX: merely a placeholder of real decay lines processing
    my @gammas;
    foreach my $g (values %g_lines) {
      foreach (@$g) {
	my $g_energy = numfilter(substr($_, 9, 10));
	my $ti = numfilter(substr($_, 64, 10));
	my $ri = numfilter(substr($_, 21, 8));
	push @gammas, {energy => $g_energy, intensity => $ri};
      }
    }
    my $decay =
      {file => $file,
       parent_a => $parent_a,
       parent_z => $parent_z,
       parent_name => $parent_nucl,
       final_a => $final_a,
       final_z => $final_z,
       final_name => $final_nucl,
       type => $dtype,
       P => $p_energy,
       halflife => $halflife,
       halflife_sec => ($halflife eq "*" ? "*" : atohl($halflife)),
       spin_parity => $spin_parity,
       G => \@gammas,
       L => \@final_levels,
       I => \%final_level_intensities,
       NR => $nr,
       NT => $nt,
       BR => $br,
       NB => $nb};
    print "$$decay{type} decay of $$decay{parent_name} to $$decay{final_name}, levels: ";
    print "$_ ", $$decay{I}->{$_}, "; " foreach keys %{$$decay{I}};
    print "BR = $$decay{BR}, NB = $$decay{NB}.\n";
    push @{$decay_table{"$parent_a-$parent_z"}}, $decay;
  }
  print "addding data from $file ", scalar(@{$decay_table{"$parent_a-$parent_z"}}), " decays for $parent_a",
    uc($ztonam{lc($parent_z)}), ".\n";
 cleanup:
  close(FILE);
}

# now build table of nuclear level information
my %gamma_transitions;
for my $file (@ensdf_files) {
  open(LFILE, "< $ensdf/$file");
  $_=<LFILE>;
  goto cleanup unless (/ADOPTED LEVELS/);
  my $nucl = substr($_, 0, 5);
  $nucl =~ /^\s*(\d+)(\w+)\s*$/;
  my $a = $1;
  my $z = $namtoz{lc($2)};
  my %levels;
  my %gammas;
  my %all_gammas;
  my %t12;
  my %t12_sec;
  my %metastable;
  my $current_level;
  my %spin_parity;
  my %normalization;
  while (my $line = <LFILE>) {
    if ($line =~ /^.....  L /) { # read level record
      # read and filter energy field
      my $e = numfilter(substr($line, 9, 10));
      # read and filter de field
      my $de = numfilter(substr($line, 19, 2));
      # read and filter metastable flag
      (my $m_flag = substr($line, 77, 2)) =~ s/^\s*(\S+)\s*$/$1/;
      # read and filter level half-life time
      (my $t12 = substr($line, 39, 10)) =~ s/^\s*(\S+)\s*$/$1/;
      unless ($t12 =~ /\S/) {
	$t12 = "*";
      }
      (my $spin_parity = substr($line, 21, 18)) =~ s/^\s*(\S+)\s*$/$1/;
      $current_level = $e;
      $levels{$current_level} = decipher_de($current_level, $de);
      $gammas{$current_level} = []; # an array of gammas
      $metastable{$current_level} = $m_flag if ($m_flag =~ /M/);
      $t12{$current_level} = $t12;
      $t12_sec{$current_level} = atohl($t12) unless ($t12 eq "*" || $t12 =~ /STABLE/i);
      $spin_parity{$current_level} = $spin_parity;
    } elsif ($line =~ /^..... PN/) {
      my $nrbr = numfilter(substr($line, 9, 10));
      $normalization{NRBR} = $nrbr if $nrbr > 0;
#      print $line;
    } elsif ($line =~ /^.....  N( |\d)/) {
#      print $line;
      my ($nr, $nt, $br, $nb);
      $nr = numfilter(substr($line, 9, 10));
      $nt = numfilter(substr($line, 21, 8));
      $br = numfilter(substr($line, 31, 8));
      $nb = numfilter(substr($line, 41, 8));
      $normalization{NR} = $nr;
      $normalization{BR} = $br;
      unless ($normalization{NRBR} > 0) {
	$normalization{NRBR} = $nr * $br;
      }
    } elsif ($line =~ /^.....  G / && defined($current_level)) { # gamma record
      # read and filter energy field
      my $e = numfilter(substr($line, 9, 10));
      # read and filter de field
      my $de = numfilter(substr($line, 19, 2));
      # read and filter intensity field
      my $ri = numfilter(substr($line, 21, 8));
      my $ti = numfilter(substr($line, 64, 10));
      if ($ri > 0 || $ti > 0) {
	push @{$gammas{$current_level}},
	  {energy => $e, ti => $ti, ri => $ri, de => decipher_de($e, $de)};
      }
      # XXX: also special case: T12 > 0 and I is 0 and # of transitions is 1
      push @{$all_gammas{$current_level}},
	{energy => $e, ti => $ti, ri => $ri, de => decipher_de($e, $de)};
    }
  }
  for my $l1 (keys %all_gammas) {
    my @g1 = @{$all_gammas{$l1}};
    if (scalar(@g1) == 1 && $t12_sec{$l1} > 0 && $t12_sec{$l1} < 0.1
	&& !($g1[0]->{ti} > 0 || $g1[0]->{ri} > 0)) {
      my %t1 = %{$g1[0]};
      print "warning:$file: only one gamma level from $l1, changing I to 1.\n";
      $t1{ti} = 1;
      $t1{ri} = 1;
      push @{$gammas{$l1}}, \%t1;
    }
  }
  $gamma_transitions{"$a-$z"} =
    {levels => \%levels,
     gammas => \%gammas,
     halflife => \%t12,
     halflife_sec => \%t12_sec,
     metastable => \%metastable,
     spin_parity => \%spin_parity,
     normalization => \%normalization,
     file => $file
    };
  print "adding data for ", uc("$a$ztonam{$z}"), ", ", scalar(keys %levels), " levels, ", 
    scalar keys %gammas, " transitions (file $file)\n";
 cleanup:
  close FILE;
}

print "writing gamma table.\n";
store \%decay_table, 'ensdf-decays.dat';
print "writing decays table.\n";
store \%gamma_transitions, 'ensdf-gamma-transitions.dat';