#!/tools32/bin/perl
#!/usr/local/bin/perl

BEGIN {
push @INC, qw( /net/db/src/SSRCalc/ssrcalc .  /tools32/lib/perl5/5.8.0/i386-linux-thread-multi /tools32/lib/perl5/5.8.0 /tools32/lib/perl5/site_perl/5.8.0/i386-linux-thread-multi /tools32/lib/perl5/site_perl/5.8.0 /tools32/lib/perl5/site_perl );
}

###############################################################################
# Program     : ShowOneSpectrum.cgi
#
# Description : This CGI program recieves consensus_library_spectrum_id and plots it.
#               Has feedback loop for changing plot.
#
# Based upon the ShowSpectrum.cgi in the Proteomics module by 
# Kerry & Eric Deutsch
#
###############################################################################


###############################################################################
# Basic SBEAMS setup
###############################################################################
use strict;
use FindBin;

use lib "$FindBin::Bin/../../lib/perl";
use lib '/tools32/lib/perl5/5.8.0/i386-linux-thread-multi';
use lib '/tools32/lib/perl5/5.8.0';
use lib '/tools32/lib/perl5/site_perl/5.8.0/i386-linux-thread-multi';
use lib '/tools32/lib/perl5/site_perl/5.8.0';
use lib '/tools32/lib/perl5/site_perl';

use vars qw ($q $sbeams $sbeamsMOD $PROG_NAME
             $current_username $modification_helper );

use SBEAMS::Connection qw($q $log);
use SBEAMS::Connection::Settings;
use SBEAMS::Connection::Tables;

use SBEAMS::PeptideAtlas;
use SBEAMS::PeptideAtlas::Settings;
use SBEAMS::PeptideAtlas::Tables;
use SBEAMS::PeptideAtlas::ModificationHelper;
use SBEAMS::PeptideAtlas::ConsensusSpectrum;

use PGPLOT;
use PDL;
use PDL::Graphics::PGPLOT;

use File::Basename;

my $denormalize = 1;
my $curveType = 'function';

my $error_state = 0;
$sbeams = new SBEAMS::Connection;
$sbeamsMOD = new SBEAMS::PeptideAtlas;
$sbeamsMOD->setSBEAMS($sbeams);

my $consensus = new SBEAMS::PeptideAtlas::ConsensusSpectrum;
$consensus->setSBEAMS($sbeams);

$PROG_NAME="ShowOneSpectrum";

###############################################################################
# Define global variables if any and execute main()
###############################################################################

my $ymax = 10;

main();


###############################################################################
# Main Program:
#
# If $sbeams->Authenticate() succeeds, print header, process the CGI request,
# print the footer, and end.
###############################################################################
sub main {


  #### Do the SBEAMS authentication and exit if a username is not returned
  exit unless ($current_username = $sbeams->Authenticate(
    permitted_work_groups_ref=>['PeptideAtlas_user','PeptideAtlas_admin',
      'PeptideAtlas_readonly', 'PeptideAtlas_exec'],
    #connect_read_only=>1,
    allow_anonymous_access=>1,
  ));

  #### Print the header, figure and do what the user want, and print footer
  my %parameters;
  $sbeams->parse_input_parameters(q=>$q,parameters_ref=>\%parameters);
  $sbeams->processStandardParameters(parameters_ref=>\%parameters);

  $denormalize = 0 if $parameters{normalize};
  $curveType = $parameters{curveType} if $parameters{curveType};
 
  my %results;
  my $ce_vals;
  my $has_median = 0;
  my $max_intensity = 0;
  # First, grab da peaks
  use Data::Dumper;
  for my $ce ( qw( low mlow medium mhigh high avg  ) ) {
    my $ce_results = mk_plot( %parameters, ce => $parameters{$ce} );
    next unless $ce_results->{max_intensity};

    # Going to base peaks to watch on this guy...
    $has_median++ if $ce eq 'medium';

    $ce_vals ||= $sbeamsMOD->get_Agilent_ce( %{$ce_results} );

    $ce_results->{ce} = $ce_vals->{$ce};
    $ce_results->{level} = ucfirst($ce);

    $parameters{modified_sequence} = $ce_results->{modified_sequence};
    $parameters{charge} = $ce_results->{charge};
    $max_intensity = $ce_results->{max_intensity} if $ce_results->{max_intensity} > $max_intensity; 

    $results{spectra} ||= {};
    $results{spectra}->{$ce} = $ce_results;
  }

  # This has to go into the header
  my $peak_plot = '';
  if ( $has_median ) {
    my $num_ions = get_num_ions( $parameters{modified_sequence} );
    my $chart_data = get_top_peaks( \%results, $num_ions );
    
    # Trying to get CE and spectra plots on same scale... 
    $parameters{ymax} = $ymax * 1.075;

    $peak_plot = get_chart( data => $chart_data,seq => $parameters{modified_sequence}, chg => $parameters{charge} );
  }

  # The plot them, now that we know the scale 
  for my $ce ( qw( low mlow medium mhigh high ) ) {
    next unless $results{spectra}->{$ce};
    my $plot_results = mk_plot( %parameters, ce => $parameters{$ce}, spectrum => $results{spectra}->{$ce} );
    $results{$ce}->{image} = $plot_results->{image};
    $results{$ce}->{lib_desc} = $plot_results->{lib_desc};
    $results{$ce}->{modified_sequence} = $plot_results->{modified_sequence};
    $results{$ce}->{charge} = $plot_results->{charge};
    $results{$ce}->{mz} = $plot_results->{precursor_mass};
  }


  my $content = "<TABLE>\n";
  $content .= "<TR><TD ALIGN=CENTER> <div id='chart_div'></div> </TD></TR>\n";
  $content .= "<TR><TD ALIGN=CENTER><H2>Spectra for $parameters{modified_sequence} <SUP>+$parameters{charge}</SUP></H2></TD></TR>\n";
  for my $ce ( qw( low mlow medium mhigh high avg  ) ) {
    next unless $results{$ce};
    $content .= "<TR><TD><IMG SRC='$results{$ce}->{image}' BORDER=0></TD></TR>\n";
    $content .= "<TR><TD ALIGN=CENTER><H4>CE: $results{spectra}->{$ce}->{ce}, Level: $results{spectra}->{$ce}->{level}, Library: $results{$ce}->{lib_desc}</TD></TR>\n";
  }
  $content .= "</TABLE>\n";


  $sbeamsMOD->display_page_header( header_info => $peak_plot );
  print $content;
  $sbeamsMOD->display_page_footer();

  $sbeams->display_page_footer(close_tables=>'YES',
    separator_bar=>'YES',display_footer=>'NO');

} # end main


sub mk_plot {

    my %parameters = @_;
	  $parameters{consensus_library_spectrum_id} = $parameters{ce};

    my %results = ( image => '', lib_desc => '' );
    return \%results unless $parameters{ce};

    #### Define some general variables
    my ($i,$element,$key,$value,$sql);


    my $apply_action  = $q->param('apply_action');

    my @charge;
    push (@charge, 1);
    push (@charge, 2);
    push (@charge, 3);

    $parameters{'ionlab'} = "Horizontal" unless $parameters{'ionlab'};
    my ($labangle,$fjust);
    if ($parameters{'ionlab'} eq "Vertical") {
      $labangle = 90;
      $fjust = 0;
    } else {
      $labangle = 0;
      $fjust = 0.5;
    }

    my $sql = qq~
       SELECT sequence,charge,mz_exact,protein_name,
       protein_name_alt,modified_sequence, consensus_library_name, file_path, entry_idx
       FROM $TBAT_CONSENSUS_LIBRARY_SPECTRUM CLS
       JOIN $TBAT_CONSENSUS_LIBRARY CL ON CL.consensus_library_id = CLS.consensus_library_id
       WHERE consensus_library_spectrum_id = '$parameters{consensus_library_spectrum_id}'
    ~;

    my @rows = $sbeams->selectSeveralColumns($sql);

    foreach my $row (@rows) {
      my ($seq, $chg, $mz, $pn, $pn_alt, $mod_seq, $lib_name, $file_path, $entry_idx ) = @{$row};
      $parameters{'peptide'} = $seq;
      $parameters{'charge'} = $chg;
      $parameters{'precursor_mass'} = $mz;
      $parameters{'protein_name'} = $pn;
      $parameters{'protein_name_alt'} = $pn_alt;
      $parameters{'modified_sequence'} = $mod_seq;
      $parameters{'library_name'} = $lib_name;
      $parameters{'file_path'} = $file_path;
      $parameters{'entry_idx'} = $entry_idx;
    }

    $results{modified_sequence} = $parameters{modified_sequence};
    $results{charge} = $parameters{charge};
    $results{mz} = $parameters{'precursor_mass'};

    $results{lib_desc} = "$parameters{library_name}";
    $results{lib_desc} =~ s/\.sptxt//;

    my $peptide = $parameters{peptide};
    $peptide =~ s/^.\.//;
    $peptide =~ s/\..$//;
    my $precursor_mass = $parameters{'precursor_mass'};
    my $charge = $parameters{'charge'};
    my $protein_name = $parameters{'protein_name'};
    my $protein_name_alt = $parameters{'protein_name_alt'};
    my $modified_sequence = $parameters{'modified_sequence'};

    $parameters{assumed_charge} = $charge;

    ## if no modified_sequence, use unmodified sequence
    unless ($modified_sequence) {
        $modified_sequence = $parameters{peptide};
    }

    #### Calculate peptide mass information
    my $masstype = $parameters{masstype} || 0;

    $modification_helper = new SBEAMS::PeptideAtlas::ModificationHelper();

    if ( !$parameters{spectrum} ) {
      my $peaks = $consensus->get_spectrum_peaks( %parameters, denormalize => 1 );
      $peaks->{mz} = $results{mz};
      $peaks->{modified_sequence} = $results{modified_sequence};
      $peaks->{charge} = $results{charge};
      return $peaks;
    }

    my %spectrum = %{$parameters{spectrum}};

    unless (%spectrum) {
			$error_state++;
      return \%results;
    }

    my ($i,$mass,$intensity,$massmin,$xticks);
    my ($massmax,$intenmax)=(0,0);
    my @spectrum_array;

    for ($i=0; $i<$spectrum{n_peaks}; $i++) {
      $mass = $spectrum{masses}->[$i];
      $intensity = $spectrum{intensities}->[$i];
      push(@spectrum_array,[$mass,$intensity]);
      $massmin = $mass if ($i == 0);
      $massmax = $mass if ($mass > $massmax);
      $intenmax = $intensity if ($intensity > $intenmax);
    }

    #### Compute data and plot bounds
    $parameters{xmin} = int($massmin/100)*100 unless $parameters{xmin};
    $parameters{xmax} = int($massmax/100)*100+100 unless $parameters{xmax};

    $parameters{ymax} ||= $intenmax + 600;
    $ymax = $parameters{ymax};

    if ($q->param("reset") eq "ZOOM OUT") {
      $parameters{xmin} = int($massmin/100)*100;
      $parameters{xmax} = int($massmax/100)*100+100;
    }

    my $maxval = $intenmax;
    my $interval = $parameters{ymax} / 20;
    my $interval_power = int( log($interval) / log(10) );
    my $roundval = 10**$interval_power;
    $parameters{ymax} = int( $parameters{ymax} /$roundval)*$roundval;
    my $ydiv = $parameters{ymax} / 5;

    #### Calculate fragment ions for the given peptide
    my @residues = Fragment($peptide);
    my $length = $#residues + 1;

    #### Initialize the plot environment
    my($progname)= basename $0;
    #my($tmpfile) = "$progname.$$.@{[time]}.gif";
    #### Reduce length because of PGPLOT 80 char limit??
    my($tmpfile) = "Spec.$$.@{[time]}.$parameters{ce}.gif";

    $parameters{gifwidth} = 800 unless $parameters{gifwidth};
    $parameters{gifheight} = 400 unless $parameters{gifheight};

    if ($apply_action eq "PRINTABLE FORMAT") {
      $parameters{gifwidth} = 480;
      $parameters{gifheight} = 384;
    }

    #print "Writing GIF to: $PHYSICAL_BASE_DIR/tmp/images/$tmpfile\n";
    my $win = pg_setup(Device=>"$PHYSICAL_BASE_DIR/tmp/images/$tmpfile/gif",
        title=>"$modified_sequence",
        xmin=>$parameters{xmin}, xmax=>$parameters{xmax},
        ymax=>$parameters{ymax}, ydiv=>$ydiv, nyticks=>5,
        gifwidth=>$parameters{gifwidth},
        gifheight=>$parameters{gifheight}
    );

    pgsch 0.9;
    pgmtext 'T',0.7,.01,0,"Peak value = $maxval";
#   pgmtext 'RV',-2,1.02,1,"Assumed charge = $parameters{assumed_charge}";
    pgmtext 'RV',-2,1.02,1,"Charge = $parameters{charge}";

    my @peakcolors;
    my @theoretical_spectrum;

    #### Loop over each desired charge
    my $charge;
    foreach $charge (@charge) {
      my ($masslist_ref) = CalcIons(Residues=>\@residues, Charge=>$charge,
                                    modified_sequence=>$modified_sequence);
      my ($B_ref,$Y_ref);

      #### Make the plot
      ($win,$B_ref,$Y_ref) = PlotPeaks(SpecData=>\%spectrum,
                                       Masslist=>$masslist_ref, Charge=>$charge,
                                       Win=>$win, Length=>$length,
                                       PeakColors=>\@peakcolors);

      #### Print out the peptide ion information
      printIons(masslist_ref=>$masslist_ref,color=>1,html=>1,
        charge=>$charge,length=>$length,
        theoretical_spectrum_ref=>\@theoretical_spectrum);

      #### Label the peaks on the plot
      LabelResidues(Ionmasses=>$masslist_ref, Binten=>$B_ref, Yinten=>$Y_ref,
                    Charge=>$charge, Win=>$win, Length=>$length,
                    Xmin=>$parameters{xmin}, Xmax=>$parameters{xmax},
                    Ymax=>$parameters{ymax}, Angle=>$labangle, Fjust=>$fjust);

    } # end foreach


    #### Draw Precursor mass symbol
    if ($parameters{precursor_mass}) {
      pgsci 2;    #### Default color
      pgsch 2;    #### Character size
      pgslw 5;    #### Line thinkness
      pgsclp(0);   #### Disable clipping
      pgpt(1,$parameters{precursor_mass},-0.34 * $interval,7);
      pgsch 1;    #### Character size
      pgpt(1,$parameters{precursor_mass},-0.34 * $interval,7);
    }


    #### Finish and close the plot
    $win->close();

	  $results{image} = "$HTML_BASE_DIR/tmp/images/$tmpfile";
    return \%results;

} # end printEntryForm

sub get_top_peaks {
  my $results = shift;
  my $num_peaks = shift || 7;
  return '' unless $results->{spectra}->{medium};
  my %top_peaks;
  $denormalize = 0;
  my @medium = $consensus->get_top_n_peaks( n_peaks => $num_peaks,
                                        denormalize => $denormalize,
                                            spectra => $results->{spectra}->{medium} );
                                                     
  $top_peaks{medium} = $medium[0];
  my $peak_list = $medium[1];
  $ymax = $top_peaks{medium}->{$peak_list->[0]};

  for my $ce ( qw( low mlow mhigh high ) ) {
    $top_peaks{$ce} = $consensus->get_top_n_peaks( n_peaks => $num_peaks,
                                               denormalize => $denormalize,
                                                   spectra => $results->{spectra}->{$ce},
                                                 peak_list => $peak_list );
    
    next unless $top_peaks{$ce};
    for my $ion ( keys( %{$top_peaks{$ce}} ) ) {
      next unless $ion;
      $ymax = $top_peaks{$ce}->{$ion} if $top_peaks{$ce}->{$ion} > $ymax;
    }
  }


  my $columns .= " data.addColumn( 'string', 'CE' )\n";
  for my $ion ( @{$peak_list} ) {
    $columns .= " data.addColumn( 'number', '$ion' )\n";
  }
#        data.addRows(4);
#        data.setValue(0, 0, '2004');
#        data.setValue(0, 1, 1000);
#        data.setValue(0, 2, 400);

  my $plot_data = '';
  my %plot_data;
  my $num_rows = 0;
  my $row = 0;
  for my $ce ( qw( low mlow medium mhigh high ) ) {
    next unless $results->{spectra}->{$ce};
    my $col = 0;
    $num_rows++;
    my $ce_label = ( $ce eq 'medium' ) ? "$results->{spectra}->{$ce}->{ce}*" : $results->{spectra}->{$ce}->{ce};
    $plot_data .= "data.setValue( $row, $col, '$ce_label' );\n";
    $col++;
    for my $ion ( @{$peak_list} ) {
      $top_peaks{$ce}->{$ion} = sprintf( "%i", $top_peaks{$ce}->{$ion} );

      $plot_data .= "data.setValue( $row, $col, $top_peaks{$ce}->{$ion} );\n";
      $col++;
    }
    $row++;
  }
  my $chart_data = qq~
    $columns
    data.addRows($num_rows);
    $plot_data
  ~;
  return $chart_data;
}

sub get_chart {

  my %args = @_;

  my $data = $args{data} || return '';
  my $ion_info = "$args{seq}+$args{chg}" || '';
  
  my $chart = qq~
    <script type="text/javascript" src="https://www.google.com/jsapi"></script>
    <script type="text/javascript">
      google.load("visualization", "1", {packages:["corechart"]});
      google.setOnLoadCallback(drawChart);
      function drawChart() {
        var data = new google.visualization.DataTable();
        $data
        var chart = new google.visualization.LineChart(document.getElementById('chart_div'));

//         chart.draw(data, {width: 800, height: 600, title: 'Intensity vs CE'});

        chart.draw(data, {width: 800, height: 600, title: 'Intensity vs CE, $ion_info',
                          hAxis: {title: 'CE value (Volts)             * is vendor optimum', showTextEvery: 1,
                                     textStyle: {color: 'black', fontSize: 16, fontWeight: 'bold' },
												          },
                          vAxis: {title: 'Intensity', minValue: 0, maxValue: $ymax },
                          pointSize: 7, curveType: '$curveType'
                         }
            );
      }
    </script>
  ~;
  return $chart;

#  <body>
#    <div id="chart_div"></div>
#  </body>


}

###############################################################################
# get_library_spectrum
###############################################################################
sub get_library_spectrum 
{
  my %args = @_;

  my $inputfile = $args{'inputfile'} || "";
  my $verbose = $args{'verbose'} || "";
  my $consensus_library_spectrum_id = $args{'consensus_library_spectrum_id'} || "";
	
  # For backwards compatibility with NIST tables (TIQAM). 
  $consensus_library_spectrum_id ||= $args{NIST_library_spectrum_id};


  #### Define some general variables
  my ($i,$element,$key,$value,$sql);
  my (@rows,$nrows);


  #### Define the data hash
  my %spectrum;
  my @mass_intensities;


  #### If we have a msms_spectrum_id, get the data from the database
  if ($consensus_library_spectrum_id) {

    #### Extract the actual mass,intensity pairs from database
    $sql = "SELECT mz, relative_intensity ".
           "  FROM $TBAT_CONSENSUS_LIBRARY_SPECTRUM_PEAK ".
           " WHERE consensus_library_spectrum_id = '$consensus_library_spectrum_id'";

    # For backwards compatibility with NIST tables (TIQAM). 
		$sql = $sbeams->evalSQL( $sql );			 
		if ( $args{NIST_library_spectrum_id} ) {
			$sql =~ s/consensus/NIST/gm;
 		  $log->info( "Translated to NIST tables: $sql" );
		}

    my @mass_intensities = $sbeams->selectSeveralColumns($sql);


    #### If we still have no spectrum data, then bail out
    unless (@mass_intensities) {
      print "\nERROR: Unable to get m/z,intensity pairs for ".
	"consensus_library_spectrum_id '$consensus_library_spectrum_id'.\n\n";
      $error_state++;
      return;
    }

    #### Extract rows into two arrays of masses and intensities
    my (@masses,@intensities);
    for ($i=0; $i<=$#mass_intensities; $i++) {
      push(@masses,$mass_intensities[$i]->[0]);
      push(@intensities,$mass_intensities[$i]->[1]);
    }

    $spectrum{n_peaks} = $#mass_intensities + 1;


    #### Put data into hash and return
    $spectrum{masses} = \@masses;
    $spectrum{intensities} = \@intensities;

    return %spectrum;

  #### Otherwise complain and return
  } else {
    print "\nERROR: Unable to determine which consensus_library_spectrum_id to load.\n\n";
    $error_state++;
    return;
  }

}

sub get_num_ions {
  my $sequence = shift;
  $log->info( "initial seq is $sequence" );
  $sequence =~ s/\[\d+\]//g;
  my $length = length $sequence;
  my $n_ions = int( $length/2 + 3.1 ) ;
  $n_ions = 12 if $n_ions > 12;
  $log->info( "mod seq is $sequence, length is $length, n_ions is $n_ions" );
  return $n_ions;
}


###############################################################################
# pg_setup
###############################################################################
sub pg_setup {
    my %args = @_;

    #### Default device is to screen (xserver)
    my $device = $args{'Device'} || "\/xs";
    #$device = "/xs";

    #### Plot title
    my $title = $args{'title'} || "";

    #### Default x limits are (0,2000)
    my $xmin = $args{'xmin'} || 0;
    my $xmax = $args{'xmax'} || 2000;

    #### Default y limits are (0,500000)
    my $ymin = $args{'ymin'} || 0;
    my $ymax = $args{'ymax'} || 500000;

    #### Default separation between ticks is 100000
    my $ytickdiv = $args{'ydiv'} || 100000;

    #### Default number of y ticks
    my $nyticks = ($args{'nyticks'}+1) || 4;

    #### Default image size is 640x480
    my $gifwidth = $args{'gifwidth'} || 640;
    my $gifheight = $args{'gifheight'} || 480;

    #### Set needed PGPLOT environment variables
    $ENV{"PGPLOT_GIF_WIDTH"} = $gifwidth;
    $ENV{"PGPLOT_GIF_HEIGHT"} = $gifheight;

    #### Try to set the proper location of rgb.txt
    if ($CONFIG_SETTING{PGPLOT_RGBTXT}) {
      $ENV{"PGPLOT_RGB"} = $CONFIG_SETTING{PGPLOT_RGBTXT};
    } elsif ( -e "/usr/share/X11/rgb.txt" ) {
      $ENV{"PGPLOT_RGB"} = "/usr/share/X11/rgb.txt";
    } else {
      $ENV{"PGPLOT_RGB"} = "/usr/X11R6/lib/X11/rgb.txt";
    }

    $ENV{"PGPLOT_FONT"} = $CONFIG_SETTING{PGPLOT_FONT} ||
      "/usr/local/lib/grfont.dat";
    $ENV{"PGPLOT_BACKGROUND"} = "white";

    #### Create a new graphics device
    #### WARNING: If the device is a filename, it apparently gets
    #### truncated at 80 characters!!!
    my $win = PDL::Graphics::PGPLOT::Window -> new({Device => "$device"});

    #### Set window limits
    pgswin $xmin, $xmax, 0, $ymax;

    #### Set viewport limits
    pgsvp .095,.9775,.065,.95;

    #### Set axis color to black (stealing lt. gray color)
    pgscr 15, 0,0,0;

    #### Set color index
    pgsci 15;

    #### Set character height
    pgsch 0.9;

    #### Set line width
    pgslw 1;

    #### Set character font (Normal)
    pgscf 1;

    #### Draw labeled frame around viewport: full frame (BC), labels on
    #### bottom and left of frame (N), major tick marks (T), y labels
    #### normal to y-axis (V), decimal labels instead of scientific
    #### notation (1), automatic x major ticks, $ytickdiv between y ticks,
    #### with $nyticks major divisions.
    pgbox 'BCNT',0,0,'BCNTV1',$ytickdiv,$nyticks;

    #### Reset character height (make labels larger)
    pgsch 1;

    #### Y label on left, centered vertically along axis
    pgmtxt 'L',3.8,.5,.5,'Normalized Intensity';

    #### X label on bottom, centered vertically along axis
    pgmtxt 'B',2.25,.5,.5,'m/z';

    #### Main title above, centered vertically along top
    pgmtxt 'T',1,.5,.5,"$title";

    #### Reset character height (want in-plot labels small)
    pgsch .8;

    #### Allow overplotting of this frame
    $win -> hold;

    return $win;
}

###############################################################################
# PlotPeaks
###############################################################################
sub PlotPeaks {
    my %args = @_;

    #### Spectrum data to be plotted
    my $specdata = $args{'SpecData'};

    #### Ions to be plotted
    my $masslist_ref = $args{'Masslist'};

    #### Charge
    my $charge = $args{'Charge'};

    #### Plot frame
    my $win = $args{'Win'};

    #### Peak Colors
    my $peakcolors_ref = $args{'PeakColors'};

    #### Peak finding window
    my $window = $args{'Window'} || 2;

    my $length = $args{'Length'};
    my @Bmaxinten = (0) x $length;
    my @Ymaxinten = (0) x $length;
    my @BYmaxinten = (0) x $length;

    my @Bmass;
    my @Ymass;
    my @BYmass;
    my @Rmass;

    my @Binten;
    my @Yinten;
    my @BYinten;
    my @Rinten;
    #my @Rinten = (0) x $specdata->{n_peaks};


    my ($redcol,$bluecol,$grcol);

    #### Define pink color to be lightcoral
    pgscr 6,0.94,0.5,0.5;

    #### Define lt. blue color to be navy
    pgscr 11,0,0,.5;

    #### Define lt. green color to be DarkSeaGreen
    pgscr 10,0.56,0.74,0.56;

    #### Define red color to be red
    pgscr 2,1,0,0;

    #### Define blue color to be blue
    pgscr 4,0,0,1;

    #### Define green color to be ForestGreen
    pgscr 3,0.13,0.55,0.13;

    if ($charge == 1) {
      $redcol = 2;
      $bluecol = 4;
      $grcol = 3;
    } else {
      $redcol = 6;
      $bluecol = 11;
      $grcol = 10;
    }


    #### Convert to piddle for easy sub-selecting
    my $bdata = pdl $masslist_ref->{Bions};
    my $ydata = pdl $masslist_ref->{Yions};

    #### Draw peaks
    my $i;
    my $lineclr;
    my ($mass, $intensity);

    for ($i=0; $i<$specdata->{n_peaks}; $i++) {
      $mass = $specdata->{masses}->[$i];
      $intensity = $specdata->{intensities}->[$i];

      #### Set default line color to Black
      $lineclr = $peakcolors_ref->[$i] || 14;

      my $mainx = pdl [$mass, $mass];
      my $mainy = pdl [0, $intensity];


      #### This kludge lets me not colorize the last B and/or
      #### first Y peaks found
## problem i'ere
      set $bdata, ($length-1),-9999;
      set $ydata, 0, -9999;

      my $Bind = which($bdata >= ($mass-$window) & $bdata <= ($mass+$window));
      my $Yind = which($ydata >= ($mass-$window) & $ydata <= ($mass+$window));

      #### If there are just B fragments with mass near enough this peak
      if (($Bind !~ 'Empty') && ($Yind =~ 'Empty')) {
        $lineclr = $redcol;
        push(@Binten,$intensity);
        push(@Bmass,$mass);
        if ($Bmaxinten[$Bind->at(0)] < $intensity) {
          $Bmaxinten[$Bind->at(0)] = $intensity;
          $masslist_ref->{Bcolor}->[$Bind->at(0)] = $lineclr
            unless ($masslist_ref->{Bcolor}->[$Bind->at(0)]);
        }

      #### Else if there are just Y fragments with mass near enough this peak
      } elsif (($Yind !~ 'Empty') && ($Bind =~ 'Empty')) {
        $lineclr = $bluecol;
        push(@Yinten,$intensity);
        push(@Ymass,$mass);
        if ($Ymaxinten[$Yind->at(0)] < $intensity) {
          $Ymaxinten[$Yind->at(0)] = $intensity;
          $masslist_ref->{Ycolor}->[$Yind->at(0)] = $lineclr;
        }

      #### Else if both B and Y fragments with mass near enough this peak
      } elsif (($Bind !~ 'Empty') && ($Yind !~ 'Empty')) {
        $lineclr = $grcol;
        push(@BYinten,$intensity);
        push(@BYmass,$mass);
        if ($Ymaxinten[$Yind->at(0)] < $intensity) {
          $Ymaxinten[$Yind->at(0)] = $intensity;
          $masslist_ref->{Ycolor}->[$Yind->at(0)] = $lineclr;
        }
        if ($Bmaxinten[$Bind->at(0)] < $intensity) {
          $Bmaxinten[$Bind->at(0)] = $intensity;
          $masslist_ref->{Bcolor}->[$Bind->at(0)] = $lineclr;
        }

      #### else if there are no fragments with mass near enough this peak
      } else {
        if (($peakcolors_ref->[$i] != 2) & ($peakcolors_ref->[$i] != 3) &
            ($peakcolors_ref->[$i] != 4) & ($peakcolors_ref->[$i] != 6) &
            ($peakcolors_ref->[$i] != 10) & ($peakcolors_ref->[$i] != 11)) {
          #$Rinten[$i] = $intensity
          push(@Rinten,$intensity);
          push(@Rmass,$mass);
        }

      }

      $peakcolors_ref->[$i] = $lineclr;
    }

    my ($mass2, $intensity2);
    $mass2 = $specdata->{masses};
    $intensity2 = $specdata->{intensities};

    #### Now we resort to plotting all peaks by "never lifting the pen"
    #### and drawing it all in a continuous line with line() because this
    #### is much faster
    my $rx = pdl (\@Rmass,\@Rmass,\@Rmass)->xchg(0,1)->clump(2);
    #my $rx = pdl ($mass2,$mass2,$mass2)->xchg(0,1)->clump(2);
    my $ra = [(0) x scalar(@Rinten)];
    my $ry = pdl ($ra,\@Rinten,$ra)->xchg(0,1)->clump(2);
    my $rh = {Color => 14};
    $win -> line ($rx,$ry,$rh);

    my $bx = pdl (\@Bmass,\@Bmass,\@Bmass)->xchg(0,1)->clump(2);
    my $ba = [(0) x scalar(@Binten)];
    my $by = pdl ($ba,\@Binten,$ba)->xchg(0,1)->clump(2);
    my $bh = {Color => $redcol};
    $win -> line ($bx,$by,$bh);

    my $yx = pdl (\@Ymass,\@Ymass,\@Ymass)->xchg(0,1)->clump(2);
    my $ya = [(0) x scalar(@Yinten)];
    my $yy = pdl ($ya,\@Yinten,$ya)->xchg(0,1)->clump(2);
    my $yh = {Color => $bluecol};
    $win -> line ($yx,$yy,$yh);

    my $byx = pdl (\@BYmass,\@BYmass,\@BYmass)->xchg(0,1)->clump(2);
    my $bya = [(0) x scalar(@BYinten)];
    my $byy = pdl ($bya,\@BYinten,$bya)->xchg(0,1)->clump(2);
    my $byh = {Color => $grcol};
    $win -> line ($byx,$byy,$byh);

    return ($win,\@Bmaxinten,\@Ymaxinten);
}


###############################################################################
# LabelResidues
###############################################################################
sub LabelResidues {
    my %args = @_;

    my $Ionmasses_ref = $args{'Ionmasses'};
    my $Bdata = pdl $Ionmasses_ref->{Bions};
    my $Ydata = pdl $Ionmasses_ref->{Yions};
    my $charge = $args{'Charge'};
    my $win = $args{'Win'};
    my $Binten_ref = $args{'Binten'};
    my @Binten = @$Binten_ref;
    my $Yinten_ref = $args{'Yinten'};
    my @Yinten = @$Yinten_ref;
    my $length = $args{'Length'};
    my $labht;
    my $angle = $args{'Angle'} || 0;
    my $fjust = $args{'Fjust'};
    my ($Bcol,$Ycol,$bothcol);
    my $i;
    my ($lineclr,$redcol,$bluecol,$grcol);
    my $Ymax = $args{'Ymax'};
    my $Xmin = $args{'Xmin'};
    my $Xmax = $args{'Xmax'};
    my $interval = $Ymax / 50.0;
    my $xshift = ($Xmax - $Xmin) / 200.0;

    #### Define pink color to be lightcoral
    pgscr 6,0.94,0.5,0.5;

    #### Define lt. blue color to be navy
    pgscr 11,0,0,.5;

    #### Define lt. green color to be DarkSeaGreen
    pgscr 10,0.56,0.74,0.56;

    #### Define red color to be lightcoral
    pgscr 2,1,0,0;

    #### Define blue color to be navy
    pgscr 4,0,0,1;

    #### Define green color to be DarkSeaGreen
    pgscr 3,0.13,0.55,0.13;

    if ($charge == 1) {
      $redcol = 2;
      $bluecol = 4;
      $grcol = 3;
    } else {
      $redcol = 6;
      $bluecol = 11;
      $grcol = 10;
    }

    for ($i=0; $i < $length; $i++) {
      if (($Binten[$i] != 0) && ($i != ($length-1))) {
        my $val = $Ionmasses_ref->{indices}->[$i];
        ++$val;
        my $index = "B$charge\-$val";
        my $mass = $Bdata->at($i);
        my $matchx = pdl [$mass, $mass];
        my $y = $Binten[$i];
        my $matchy = pdl [$y+($interval/3.), $y+4*($interval/3.)];
        my $Yind = which($Ydata >= ($mass-2) & $Ydata <= ($mass+2));
        if ($Yind !~ 'Empty') {
          #### Green text and line for both ion match
          pgsci $grcol;
          $lineclr = $grcol;

          #### Location of label above tick mark (moved up)
          $labht = $y+8.5*($interval/3.);
        } else {
          #### Red text and line for B ion match
          pgsci $redcol;
          $lineclr = $redcol;

          #### Location of label above tick mark
          $labht = $y+5*($interval/3.);
        }
        #### Plot ion marker line
        $win -> line($matchx, $matchy, {Color=>$lineclr});
        $win -> hold;

        #### Add ion label
        pgptext $mass+$xshift,$labht,$angle,$fjust,"$index"
          if (($labht < $Ymax) && ($mass > $Xmin) && ($mass < $Xmax));
      }
    }

    for ($i = $length; $i > 0; $i--) {
      if (($Yinten[$i] != 0) && ($i != 0)) {
        my $index = "Y$charge\-$i";
        my $mass = $Ydata->at($i);
        my $matchx = pdl [$mass, $mass];
        my $y = $Yinten[$i];
        my $matchy = pdl [$y+($interval/3.), $y+4*($interval/3.)];
        my $Bind = which($Bdata >= ($mass-2) & $Bdata <= ($mass+2));
        if ($Bind !~ 'Empty') {
          #### Green text and line for both ion match
          pgsci $grcol;
          $lineclr = $grcol;

          #### Location of label above tick mark (moved up)
          $labht = $y+8.5*($interval/3.);
        } else {
          #### Red text and line for B ion match
          pgsci $bluecol;
          $lineclr = $bluecol;

          #### Location of label above tick mark
          $labht = $y+5*($interval/3.);
        }
        #### Plot ion marker line
        $win -> line($matchx, $matchy, {Color=>$lineclr});
        $win -> hold;

        #### Add ion label
        pgptext $mass+$xshift,$labht,$angle,$fjust,"$index"
          if (($labht < $Ymax) && ($mass > $Xmin) && ($mass < $Xmax));
      }
    }

    return $win;
}

###############################################################################
# Fragment
###############################################################################
sub Fragment {
    my $peptide = shift;
    my $length = length($peptide);
    my @residues = ();
    my $i;

    for ($i=0; $i<$length; $i++) {
      if (substr($peptide,$i+1,1) eq '[') {
        push (@residues, substr($peptide,$i,6));
        $i = $i + 5;
      } elsif (substr($peptide,$i+1,1) =~ /\W/) {
        push (@residues, substr($peptide,$i,2));
        $i = $i + 1;
      } else {
        push (@residues, substr($peptide,$i,1));
      }
    }

    #### Return residue array
    return @residues;
}

###############################################################################
# CalcIons -- calculate theoretical ions (including modified masses)
###############################################################################
sub CalcIons {
    my %args = @_;
    my $i;

    my $residues_ref = $args{'Residues'};
    my @residues = @$residues_ref;
    my $charge = $args{'Charge'} || 1;
    my $length = scalar(@residues);
   
    my $modified_sequence = $args{'modified_sequence'};

    my @masses = $modification_helper->getMasses($modified_sequence);

    my $Nterm = 1.0078;
    my $Bion = 0.;
    my $Yion  = 19.0184;  ## H_2 + O

    my @Bcolor = (14) x $length;
    my @Ycolor = (14) x $length;

    my %masslist;
    my (@aminoacids, @indices, @rev_indices, @Bions, @Yions);


    #### Compute the ion masses
    for ($i = 0; $i<$length; $i++) {
      $Bion += $masses[$i];

      #### B index & Y index
      $indices[$i] = $i;
      $rev_indices[$i] = $length-$i;
      $Yion += $masses[ $rev_indices[$i] ]  if ($i > 0);

      #### B ion mass & Y ion mass
      $Bions[$i] = ($Bion + $charge*$Nterm)/$charge;
      $Yions[$i] = ($Yion + $charge*$Nterm)/$charge;
    }

    $masslist{residues} = \@residues;
    $masslist{indices} = \@indices;
    $masslist{Bions} = \@Bions;
    $masslist{Yions} = \@Yions;
    $masslist{rev_indices} = \@rev_indices;

    #### Return reference to a hash of array references
    return (\%masslist);
}

###############################################################################
# printIons
###############################################################################
sub printIons {
  my %args = @_;

  my $masslist_ref = $args{'masslist_ref'};
  my $color = $args{'color'} || 0;
  my $html = $args{'html'} || 0;
  my $charge = $args{'charge'};
  my $length = $args{'length'};
  my $theoretical_spectrum_ref = $args{'theoretical_spectrum_ref'};

# print "SEQ  #     B       Y    +$charge\n";
# print "--- --  ------  ------  --\n";

# my ($bcolbegin, $bcolend, $ycolbegin, $ycolend);

# my (%colors);
# $colors{2} = "#FF0000";
# $colors{4} = "#0000FF";
# $colors{3} = "#218D21";
# $colors{6} = "#F18080";
# $colors{11} = "#00088";
# $colors{10} = "#8FBE8F";


# #### Printing stuff
  for (my $i=0; $i < $length; $i++) {

#   #### If the output is in HTML, define the colorizing tags
#   if ($html) {

#     #### If a color for this B ion mass, set color tags
#     if ($masslist_ref->{Bcolor}->[$i] >= 2) {
#       $bcolbegin = "<FONT COLOR = $colors{$masslist_ref->{Bcolor}->[$i]}>";
#       $bcolend = "</FONT>";
#     #### else no color (default black)
#     } else {
#       $bcolbegin = "";
#       $bcolend = "";
#     }

#     #### If a color for this Y ion mass, set color tags
#     if ($masslist_ref->{Ycolor}->[$i] >= 2) {
#       $ycolbegin = "<FONT COLOR = $colors{$masslist_ref->{Ycolor}->[$i]}>";
#       $ycolend = "</FONT>";
#     #### else no color (default black)
#     } else {
#       $ycolbegin = "";
#       $ycolend = "";
#     }

#   }


    #### Define the m/z columns formats and values
    my $B_format = '%7.1f';
    my $Y_format = '%7.1f';
    my $B_value = $masslist_ref->{Bions}->[$i];
    my $Y_value = $masslist_ref->{Yions}->[$i];


    #### Special case --'s for first row
    if ($i == 0) {
      $Y_format = '%7s';
      $Y_value = '--  ';

    #### Special case --'s for last row
    } elsif ($i == ($length-1)) {
      $B_format = '%7s';
      $B_value = '--  ';
    }

#   #### Print out the data
#   printf "%3s %2d $bcolbegin$B_format$bcolend ".
#          "$ycolbegin$Y_format$ycolend %3d\n",
#     $masslist_ref->{residues}->[$i], $i+1,
#     $B_value, $Y_value, $length-$i;

    #### Fill the theoretical spectrum data in a different format
    #### (Residue,Index,Ion,Charge,m/z)
    $theoretical_spectrum_ref->[2*$length*($charge-1)+$i] =
      [$masslist_ref->{residues}->[$i],$i+1,'B',$charge,$B_value];
    $theoretical_spectrum_ref->[2*$length*($charge-1)+2*$length-1-$i] =
      [$masslist_ref->{residues}->[$i],$length-$i,'Y',
      $charge,$Y_value];


  } # end for

# print "\n";

} # end printIons