#! /usr/local/bin/perl

##############################################
# sub routine "input" : fetch FFF file name  #
#                                            #
# argument = FFF file name;                  #
# set FFF ($fff) and HK  ($hk) file names    #
##############################################
sub input {
    my $status = 0;

    $fffile     = `pquery2 hxdmkgainhist_gso input_name $invokestring`;
    chop($fffile);
    $hkfile     = `pquery2 hxdmkgainhist_gso hk_name $invokestring`;
    chop($hkfile);
    ${fitresults_gd} = `pquery2 hxdmkgainhist_gso gso_gd_fitlog_name $invokestring`;
    chop(${fitresults_gd});
    ${fitresults_511} = `pquery2 hxdmkgainhist_gso gso_511_fitlog_name $invokestring`;
    chop(${fitresults_511});
    ${fitresults_152gd} = `pquery2 hxdmkgainhist_gso gso_152gd_fitlog_name $invokestring`;
    chop(${fitresults_152gd});

    ${HV_W0_VOLT_MIN} = `pquery2 hxdmkgainhist_gso hv_w0_volt_min $invokestring`;
    chop(${HV_W0_VOLT_MIN});
    ${HV_W1_VOLT_MIN} = `pquery2 hxdmkgainhist_gso hv_w1_volt_min $invokestring`;
    chop(${HV_W1_VOLT_MIN});
    ${HV_W2_VOLT_MIN} = `pquery2 hxdmkgainhist_gso hv_w2_volt_min $invokestring`;
    chop(${HV_W2_VOLT_MIN});
    ${HV_W3_VOLT_MIN} = `pquery2 hxdmkgainhist_gso hv_w3_volt_min $invokestring`;
    chop(${HV_W3_VOLT_MIN});

    ${fit_start_gd_s} = `pquery2 hxdmkgainhist_gso fit_start_gd_s $invokestring`;
    chop(${fit_start_gd_s});
    ${fit_end_gd_s} = `pquery2 hxdmkgainhist_gso fit_end_gd_s $invokestring`;
    chop(${fit_end_gd_s});
    ${fit_start_gd_f} = `pquery2 hxdmkgainhist_gso fit_start_gd_f $invokestring`;
    chop(${fit_start_gd_f});
    ${fit_end_gd_f} = `pquery2 hxdmkgainhist_gso fit_end_gd_f $invokestring`;
    chop(${fit_end_gd_f});
    
    ${fit_start_511_s} = `pquery2 hxdmkgainhist_gso fit_start_511_s $invokestring`;
    chop(${fit_start_511_s});
    ${fit_end_511_s} = `pquery2 hxdmkgainhist_gso fit_end_511_s $invokestring`;
    chop(${fit_end_511_s});
    ${fit_start_511_f} = `pquery2 hxdmkgainhist_gso fit_start_511_f $invokestring`;
    chop(${fit_start_511_f});
    ${fit_end_511_f} = `pquery2 hxdmkgainhist_gso fit_end_511_f $invokestring`;
    chop(${fit_end_511_f});
    
    ${fit_start_152gd_s} = `pquery2 hxdmkgainhist_gso fit_start_152gd_s $invokestring`;
    chop(${fit_start_152gd_s});
    ${fit_end_152gd_s} = `pquery2 hxdmkgainhist_gso fit_end_152gd_s $invokestring`;
    chop(${fit_end_152gd_s});
    ${fit_start_152gd_f} = `pquery2 hxdmkgainhist_gso fit_start_152gd_f $invokestring`;
    chop(${fit_start_152gd_f});
    ${fit_end_152gd_f} = `pquery2 hxdmkgainhist_gso fit_end_152gd_f $invokestring`;
    chop(${fit_end_152gd_f});

    $EXPOSURE   = `pquery2 hxdmkgainhist_gso exposure $invokestring`;
    chop($EXPOSURE);

    $process_id        = `pquery2 hxdmkgainhist_gso process_id $invokestring`;
    chop(${process_id});
     
    $save_tmpfiles   = `pquery2 hxdmkgainhist_gso save_tmpfiles $invokestring`;
    chop($save_tmpfiles);

    return $status;
}

###############################################
# sub routine "maketime" :                    #
#       process "hk" file to  "maketime"      #
#                                             #
# input HK file ($hk)                         #
###############################################
sub maketime {

    local($hk_file)="${hkfile}";
  ## make HK dump ##
    local($total_gti)="${totalgti}";
    local($maketime) = "maketime infile=\"${hk_file}+1\" outfile=\"${total_gti}\" expr=\"(HXD_HV_W0_CAL >= ${HV_W0_VOLT_MIN}) && (HXD_HV_W1_CAL >= ${HV_W1_VOLT_MIN}) && (HXD_HV_W2_CAL >= ${HV_W2_VOLT_MIN}) && (HXD_HV_W3_CAL >= ${HV_W3_VOLT_MIN})\" name=anything value=anything time=TIME compact=no prefr=0.5 postfr=0.5\n";
    print $maketime;
    system($maketime);
}

###############################################
# sub routine "makefdump" :                   #
#       process "hk" file to  "makegti"       #
#                                             #
# input HK file ($hk)                         #
# make fdump log file name ($fdump_log)       #
###############################################
sub makefdump {
    local($gti)="${totalgti}";
  ## make HK dump ##

    local($fdump) = "fdump infile=\"${gti}\" outfile=\"${fdump_log}\" columns='START,STOP' rows=- prhead=no pagewidth=100 showrow=no\n";
    print $fdump;
    system($fdump);
}

########################
##  maketime routine  ##
########################
sub makegti{
#### read HK dump ####
local(${start_time})=0.0;         # initialize gti start time
local(${end_time})=0.0;           # initialize gti end time
local(${rest_time})=$EXPOSURE;    # initialize exposure time counter
local(${g_num})=0;                # initialize gti file ID
local(${g_flg})=0;
    
    open(FDUMP,"${fdump_log}")
	|| die "Cannot open fdump.log (${fdump_log}).: $!\n";
    while(<FDUMP>){     # read new line from fdump.log
      if(/\d+/){
	  ${g_flg} = 1;
	  chomp;
	  s/^\s+//g;

	  ($start,$stop)
	      = split(/\s+/,"$_");

	  ${start_time}=${start};

	while((${start_time}+${rest_time}) <= $stop){
	    ${end_time}=${start_time}+${rest_time};

	open(GTI,">> ${working_dir}/${process_id}_gti_${g_num}.gti")
	    || die "Cannot open gti_${g_num}(${process_id}_gti_${g_num}.gti).: $!\n";
	    print GTI "${start_time} ${end_time}\n";
	close GTI;

	open(GTILIST,">> ${gti_list}") 
	    || die "Cannot open gti_log(${gti_list}).: $!\n";
	print GTILIST "${working_dir}/${process_id}_gti_${g_num}.gti\n";
	close GTILIST;
	    ${g_num}++;
	    ${start_time}=${end_time};
	    ${rest_time} = ${EXPOSURE};
	}
	open(GTI,">> ${working_dir}/${process_id}_gti_${g_num}.gti")
	    || die "Cannot open gti_${g_num}(${process_id}_gti_${g_num}.gti).: $!\n";
	    print GTI "${start_time} ${stop}\n";
	    ${rest_time} -= $stop - ${start_time};
	}
      else{;}
  }
if(${g_flg} == 1){
    open(GTILIST,">> ${gti_list}") 
	|| die "Cannot open gti_log(${gti_list}).: $!\n";
    print GTILIST "${working_dir}/${process_id}_gti_${g_num}.gti\n";
    close GTILIST;
}
else{;}
}

######################
#   SORT OUT UNITs   #
######################
sub sort {
    local($wpu);  local($wid);  local($unit);
  for($unit=0;$unit<16;$unit++){
      $wpu = int($unit/4);
      $wid = $unit - 4*$wpu;
      $unit_s = $unit;
      $unit_f = $unit + 16;

      $phagsow[$unit_s] = "${process_id}_gso_w${wpu}${wid}_gti_${loop_gti}_slow.pha";
      $phagsow[$unit_f] = "${process_id}_gso_w${wpu}${wid}_gti_${loop_gti}_fast.pha";
      ${phagsow_hitpat8[$unit_s]} = "${process_id}_gso_w${wpu}${wid}_gti_${loop_gti}_hitpat8_slow.pha";
      ${phagsow_hitpat8[$unit_f]} = "${process_id}_gso_w${wpu}${wid}_gti_${loop_gti}_hitpat8_fast.pha";

      ${phagsow_rbn16[$unit_s]} = "${process_id}_gso_w${wpu}${wid}_reb16_gti_${loop_gti}_slow.pha";
      ${phagsow_rbn16[$unit_f]} = "${process_id}_gso_w${wpu}${wid}_reb16_gti_${loop_gti}_fast.pha";
      ${phagsow_hitpat8_rbn16[$unit_s]} = "${process_id}_gso_w${wpu}${wid}_reb16_gti_${loop_gti}_hitpat8_slow.pha";
      ${phagsow_hitpat8_rbn16[$unit_f]} = "${process_id}_gso_w${wpu}${wid}_reb16_gti_${loop_gti}_hitpat8_fast.pha";

  }
}

#####################################
#    sub rutine "makeselect"        #
#                                   #
#    select PMT events every unit   #
#    make "fits.pha"                #
#####################################
sub eventselect{

    local($unit); local($wpu); local($wid);

#### Non Selection ####

    print "\nMaking Non-selection fits file now.\n";

    for($unit=0;$unit<16;$unit++){
	$wpu = int($unit/4);
	$wid = $unit - 4*$wpu;	

        local($fselect[$unit]) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_${unit}.fff\" expr=\"(UNITID==\"${unit}\")&&(TRIG==b0100000)&&(QUALITY_FLAGS==bxxxxxx1)\"\n";

        print $fselect[$unit];
        system($fselect[$unit]);

    }
    
#### Hit Pattern 8 ####

    print "\nMaking Hit Pattern 8 fits file now.\n";
    
    local(${fselect_hitpat8[0]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_0_hitpat8.fff\" expr=\"(UNITID==0) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0111000000000000) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b11100000000000000011) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[1]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_1_hitpat8.fff\" expr=\"(UNITID==1) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b1011100100000000) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b01110000000000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[2]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_2_hitpat8.fff\" expr=\"(UNITID==2) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b1101100110001100) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[3]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_3_hitpat8.fff\" expr=\"(UNITID==3) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b1110000000001100) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000000000111) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[4]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_4_hitpat8.fff\" expr=\"(UNITID==4) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0110011100000000) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00111000000000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[5]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_5_hitpat8.fff\" expr=\"(UNITID==5) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0000101100000000) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00011111000000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[6]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_6_hitpat8.fff\" expr=\"(UNITID==6) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0000110111000000) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000011100000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[7]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_7_hitpat8.fff\" expr=\"(UNITID==7) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0110111011000100) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[8]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_8_hitpat8.fff\" expr=\"(UNITID==8) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0010001101110110) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[9]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_9_hitpat8.fff\" expr=\"(UNITID==9) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0000001110110000) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000001110000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[10]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_10_hitpat8.fff\" expr=\"(UNITID==10) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0000000011010000) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000111110000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[11]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_11_hitpat8.fff\" expr=\"(UNITID==11) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0000000011100110) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000111000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[12]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_12_hitpat8.fff\" expr=\"(UNITID==12) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0011000000000111) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000000001110) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[13]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_13_hitpat8.fff\" expr=\"(UNITID==13) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0011000110011011) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000000000000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[14]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_14_hitpat8.fff\" expr=\"(UNITID==14) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0000000010011101) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000011100000) == (b00000000000000000000)))\"\n";
    local(${fselect_hitpat8[15]}) = "fselect infile=\"${fffile}\" outfile=\"${working_dir}/tmp_${process_id}_15_hitpat8.fff\" expr=\"(UNITID==15) && (TRIG==b0100000) && (QUALITY_FLAGS==bxxxxxx1) && (((HIT_PATTERN_WELL & b0000000000001110) == (b0000000000000000)) && ((HIT_PATTERN_ANTI & b00000000000001111100) == (b00000000000000000000)))\"\n";

        print  ${fselect_hitpat8[0]};
        system(${fselect_hitpat8[0]});
        print  ${fselect_hitpat8[1]};
        system(${fselect_hitpat8[1]});
        print  ${fselect_hitpat8[2]};
        system(${fselect_hitpat8[2]});
        print  ${fselect_hitpat8[3]};
        system(${fselect_hitpat8[3]});
        print  ${fselect_hitpat8[4]};
        system(${fselect_hitpat8[4]});
        print  ${fselect_hitpat8[5]};
        system(${fselect_hitpat8[5]});
        print  ${fselect_hitpat8[6]};
        system(${fselect_hitpat8[6]});
        print  ${fselect_hitpat8[7]};
        system(${fselect_hitpat8[7]});
        print  ${fselect_hitpat8[8]};
        system(${fselect_hitpat8[8]});
        print  ${fselect_hitpat8[9]};
        system(${fselect_hitpat8[9]});
        print  ${fselect_hitpat8[10]};
        system(${fselect_hitpat8[10]});
        print  ${fselect_hitpat8[11]};
        system(${fselect_hitpat8[11]});
        print  ${fselect_hitpat8[12]};
        system(${fselect_hitpat8[12]});
        print  ${fselect_hitpat8[13]};
        system(${fselect_hitpat8[13]});
        print  ${fselect_hitpat8[14]};
        system(${fselect_hitpat8[14]});
        print  ${fselect_hitpat8[15]};
        system(${fselect_hitpat8[15]});
   
}

####################################
#   make slow and fast pha files   #
####################################  
sub extractor {
    local($unit);
    for($unit=0;$unit<16;$unit++){
        local($unit_s) = $unit;
        local($unit_f) = $unit + 16;
        local($extractor[$unit_s]) = "extractor filename=\"${working_dir}/tmp_${process_id}_${unit}.fff\" eventsout=\"NONE\" imgfile=\"NONE\" fitsbinlc=\"NONE\" regionfile=\"NONE\" phafile=\"${working_dir}/${phagsow[$unit_s]}\" timefile=\"${working_dir}/${process_id}_gti_${loop_gti}.gti\" gtinam=\"GTI\" tcol=\"TIME\" ecol=\"PHA_SLOW\" xcolh=\"UNITID\" ycolh=\"DET_TYPE\" xcolf=\"UNITID\" ycolf=\"UNITID\"\n";
        local($extractor[$unit_f]) = "extractor filename=\"${working_dir}/tmp_${process_id}_${unit}.fff\" eventsout=\"NONE\" imgfile=\"NONE\" fitsbinlc=\"NONE\" regionfile=\"NONE\" phafile=\"${working_dir}/${phagsow[$unit_f]}\" timefile=\"${working_dir}/${process_id}_gti_${loop_gti}.gti\" gtinam=\"GTI\" tcol=\"TIME\" ecol=\"PHA_FAST\" xcolh=\"UNITID\" ycolh=\"DET_TYPE\" xcolf=\"UNITID\" ycolf=\"UNITID\"\n";

	print $extractor[$unit_s];
        system($extractor[$unit_s]);
	print $extractor[$unit_f];
        system($extractor[$unit_f]);

        local(${extractor_hitpat8[$unit_s]}) = "extractor filename=\"${working_dir}/tmp_${process_id}_${unit}_hitpat8.fff\" eventsout=\"NONE\" imgfile=\"NONE\" fitsbinlc=\"NONE\" regionfile=\"NONE\" phafile=\"${working_dir}/${phagsow_hitpat8[$unit_s]}\" timefile=\"${working_dir}/${process_id}_gti_${loop_gti}.gti\" gtinam=\"GTI\" tcol=\"TIME\" ecol=\"PHA_SLOW\" xcolh=\"UNITID\" ycolh=\"DET_TYPE\" xcolf=\"UNITID\" ycolf=\"UNITID\"\n";
        local(${extractor_hitpat8[$unit_f]}) = "extractor filename=\"${working_dir}/tmp_${process_id}_${unit}_hitpat8.fff\" eventsout=\"NONE\" imgfile=\"NONE\" fitsbinlc=\"NONE\" regionfile=\"NONE\" phafile=\"${working_dir}/${phagsow_hitpat8[$unit_f]}\" timefile=\"${working_dir}/${process_id}_gti_${loop_gti}.gti\" gtinam=\"GTI\" tcol=\"TIME\" ecol=\"PHA_FAST\" xcolh=\"UNITID\" ycolh=\"DET_TYPE\" xcolf=\"UNITID\" ycolf=\"UNITID\"\n";

        print ${extractor_hitpat8[$unit_s]};
        system(${extractor_hitpat8[$unit_s]});
        print ${extractor_hitpat8[$unit_f]};
        system(${extractor_hitpat8[$unit_f]});

    } 
}

########################################
# rebin pha files                      #
########################################
sub rebin{
    local($unit);
    for($unit=0;$unit<32;$unit++){

        local($rbnpha[$unit]) = "rbnpha infile=\"${working_dir}/${phagsow[$unit]}\" outfile=\"${working_dir}/${phagsow_rbn16[$unit]}\" cmpmode=linear finchan=256\n";

        local(${rbnpha_hitpat8[$unit]}) = "rbnpha infile=\"${working_dir}/${phagsow_hitpat8[$unit]}\" outfile=\"${working_dir}/${phagsow_hitpat8_rbn16[$unit]}\" cmpmode=linear finchan=256\n";
        print $rbnpha[$unit];
        system($rbnpha[$unit]);
        print ${rbnpha_hitpat8[$unit]};
        system(${rbnpha_hitpat8[$unit]});

    }
}

###############################################
# sub routine " makedummyrsp" :               #
#       make dummy response with 256 ch       #
#   no input                                  #
###############################################
sub makedummyrsp {
##### HEADAS 6.0.x version
#    $makedummyrsp_exe = "genrsp inrfil=\"none\" rmffil=\"${rmf}\" res_reln=\"constant\" fwhm=0 disperse=no tlscpe=\"SUZAKU\" instrm=\"HXD\" low=-0.5 high=255.5 delta=1.0 break=-1.0 delta2=0.1 chan_low=-0.5 chan_high=255.5 chan_delta=1.0 chan_break=-1.0 chan_delta2=0.1 efffil=\"none\" detfil=\"none\" filfil=\"none\" max_elements=100000 rsp_min=1.e-6 clobber=yes mode=ql \n";
##### HEADAS new version
    $makedummyrsp_exe = "genrsp inrfil=\"none\" rmffil=\"${rmf}\" resol_reln=\"constant\" resol_file=\"none\" fwhm=0 disperse=no tlscpe=\"SUZAKU\" instrm=\"HXD\" resp_reln=\"linear\" resp_file=\"none\" resp_low=-0.5 resp_high=255.5 resp_number=256 resp_break=-1.0 resp_bnumber=0 chan_reln=\"linear\" chan_low=-0.5 chan_high=255.5 chan_number=256 chan_break=-1.0 chan_bnumber=0 efffil=\"none\" detfil=\"none\" filfil=\"none\" max_elements=100000 rsp_min=1.e-6 clobber=yes mode=ql \n";
    print $makedummyrsp_exe;
    system($makedummyrsp_exe);
}

########################################
# drive XSPEC version 11.xx            #
########################################
sub drivexspec{
    local(${unit_id});
    $fit_model_id = 1;                  # Fit mode ID : gauss + powerlaw
    for ($unit_id=0;$unit_id<16;$unit_id++){
	local($wpu) = int($unit_id/4);
	local($wid) = ($unit_id - 4*$wpu);
	local($num_s) = ${unit_id};
	local($num_f) = ${unit_id} + 16;
	local(${log_gd}) = "${working_dir}/${process_id}_xspec_w${wpu}${wid}_Gd_gti_${loop_gti}.log";
	local(${log_511}) = "${working_dir}/${process_id}_xspec_w${wpu}${wid}_511_gti_${loop_gti}.log";
	local(${log_152gd}) = "${working_dir}/${process_id}_xspec_w${wpu}${wid}_152gd_gti_${loop_gti}.log";
	local(${fit_gd_ps}) = "${working_dir}/${process_id}_hxd_w${wpu}${wid}_Gd_gti_${loop_gti}.ps";
	local(${fit_511_ps}) = "${working_dir}/${process_id}_hxd_w${wpu}${wid}_511_gti_${loop_gti}.ps";
	local(${fit_152gd_ps}) = "${working_dir}/${process_id}_hxd_w${wpu}${wid}_152gd_gti_${loop_gti}.ps";
      open(COM20, "> xspec_gd_${wpu}${wid}_${loop_gti}.xcm")
	  || die "Cannot open xspec_gd.xcm.: $!\n";
        print COM20 "query yes\n";
        print COM20 "data 1:1 ${working_dir}/${phagsow_hitpat8_rbn16[$num_s]} 2:2 ${working_dir}/${phagsow_hitpat8_rbn16[$num_f]}\n";
        print COM20 "response 1 ${rmf} 2 ${rmf}\n";
      	print COM20 "setplot energy\n";
	print COM20 "ignore bad\n";
        print COM20 "ignore 1:**-${fit_start_gd_s} 1:${fit_end_gd_s}-** 2:**-${fit_start_gd_f} 2:${fit_end_gd_f}-**\n";
	print COM20 "model gauss + powerlaw\n"; # $fit_model_id = 1
        print COM20 "90.0\n";      # Gd slow:initial K alpha line energy in ch
	print COM20 "14.0\n";      # Gd slow:initial sigma in ch
	print COM20 "\n";          # Gd slow:gauss normalization
	print COM20 "0.9\n";       # Gd slow:powerlaw PhoIndex
        print COM20 "\n";          # Gd slow:powerlaw normalization
        print COM20 "90.0\n";      # Gd fast:initial K alpha line energy in ch
	print COM20 "14.0\n";      # Gd fast:initial sigma in ch
	print COM20 "\n";          # Gd fast:gauss normalization
	print COM20 "0.9\n";       # Gd fast:powerlaw PhoIndex
        print COM20 "\n";          # Gd fast:powerlaw normalization
        print COM20 "log ${log_gd}\n";
	print COM20 "fit 200\n";
	print COM20 "chatter 0 0\nfreeze 2 3 4 5 6 7 9\nchatter 10 10\n";
        print COM20 "chatter 0 0\nfit\nchatter 10 10\n";
        print COM20 "query no\n";
	print COM20 "error 1\n";
        print COM20 "query yes\n";
	print COM20 "chatter 0 0\nfreeze 1\nchatter 10 10\n";
        print COM20 "chatter 0 0\nthaw 6\nchatter 10 10\n";
        print COM20 "chatter 0 0\nfit\nchatter 10 10\n";
        print COM20 "query no\n";
        print COM20 "error 6\n";
        print COM20 "query yes\n";
        print COM20 "chatter 0 0\nthaw 1 2 3 4 5 7 9\nchatter 10 10\n";
        print COM20 "chatter 0 0\nfit\nchatter 10 10\n";
	print COM20 "flux ${fit_start_gd_s} ${fit_end_gd_s}\n";
        print COM20 "puts \"exposure_time=[format \"\%E\" [tcloutr expos 1]]\"\n";
	print COM20 "show files\n";
	print COM20 "setplot add\n";
	print COM20 "plot ldata delchi\n";
        print COM20 "iplot\n";
	print COM20 "r y 0.001 0.025\n";
        print COM20 "hardcopy ${fit_gd_ps}/cps\n";
  	print COM20 "plot\n";
	print COM20 "quit\n";
	print COM20 "exit\n";
	print COM20 "\n";
      close COM20;

    ${xspec_gd} = "xspec < xspec_gd_${wpu}${wid}_${loop_gti}.xcm";
    print ${xspec_gd};
    system(${xspec_gd});

    #unlink COM20;

      open(COM21, "> xspec_511_${wpu}${wid}_${loop_gti}.xcm")
	  || die "Cannot open xspec_511.xcm.: $!\n";
        print COM21 "query yes\n";
        print COM21 "data 1:1 ${working_dir}/${phagsow_rbn16[$num_s]} 2:2  ${working_dir}/${phagsow_rbn16[$num_f]}\n";
        print COM21 "backgrnd 1:1 ${working_dir}/${phagsow_hitpat8_rbn16[$num_s]} 2:2 ${working_dir}/${phagsow_hitpat8_rbn16[$num_f]}\n";
        print COM21 "response 1 ${rmf} 2 ${rmf}\n";
      	print COM21 "setplot energy\n";
	print COM21 "ignore bad\n";
        print COM21 "ignore 1:**-${fit_start_511_s} 1:${fit_end_511_s}-** 2:**-${fit_start_511_f} 2:${fit_end_511_f}-**\n";
	print COM21 "model gauss + powerlaw\n"; # $fit_model_id = 1
        print COM21 "130.0\n";     # 511 slow:initial K alpha line energy in ch
	print COM21 "5.0\n";       # 511 slow:initial sigma in ch  
	print COM21 "\n";          # 511 slow:gauss normalization
	print COM21 "0.9\n";       # 511 slow:powerlaw PhoIndex
        print COM21 "\n";          # 511 slow:powerlaw normalization
        print COM21 "130.0\n";     # 511 fast:initial K alpha line energy in ch
	print COM21 "5.0\n";       # 511 fast:initial sigma in ch
	print COM21 "\n";          # 511 fast:gauss normalization
	print COM21 "0.9\n";       # 511 fast:powerlaw PhoIndex
        print COM21 "\n";          # 511 fast:powerlaw normalization	      
        print COM21 "log ${log_511}\n";
	print COM21 "fit 200\n";
	print COM21 "chatter 0 0\nfreeze 2 3 4 5 6 7 9\nchatter 10 10\n";
        print COM21 "chatter 0 0\nfit\nchatter 10 10\n";
        print COM21 "query no\n";
	print COM21 "error 1\n";
        print COM21 "query yes\n";
	print COM21 "chatter 0 0\nfreeze 1\nchatter 10 10\n";
        print COM21 "chatter 0 0\nthaw 6\nchatter 10 10\n";
        print COM21 "chatter 0 0\nfit\nchatter 10 10\n";
        print COM21 "query no\n";
        print COM21 "error 6\n";
        print COM21 "query yes\n";
        print COM21 "chatter 0 0\nthaw 1 2 3 4 5 7 9\nchatter 10 10\n";
        print COM21 "chatter 0 0\nfit\nchatter 10 10\n";
	print COM21 "flux ${fit_start_511_s} ${fit_end_511_s}\n";
        print COM21 "puts \"exposure_time=[format \"\%E\" [tcloutr expos 1]]\"\n";
	print COM21 "show files\n";
  	print COM21 "setplot add\n";
	print COM21 "plot ldata delchi\n";
        print COM21 "iplot\n";
	print COM21 "r y 0.005 0.03\n";
        print COM21 "hardcopy ${fit_511_ps}/cps\n";
  	print COM21 "plot\n";
	print COM21 "quit\n";
	print COM21 "exit\n";						       
	print COM21 "\n";
      close COM21;

    ${xspec_511} = "xspec < xspec_511_${wpu}${wid}_${loop_gti}.xcm";
    print  ${xspec_511};
    system(${xspec_511});

#	unlink COM21;

      open(COM22, "> xspec_152gd_${wpu}${wid}_${loop_gti}.xcm")
	  || die "Cannot open xspec_gd.xcm.: $!\n";
        print COM22 "query yes\n";
        print COM22 "data 1:1 ${working_dir}/${phagsow_hitpat8_rbn16[$num_s]} 2:2 ${working_dir}/${phagsow_hitpat8_rbn16[$num_f]}\n";
        print COM22 "response 1 ${rmf} 2 ${rmf}\n";
      	print COM22 "setplot energy\n";
	print COM22 "ignore bad\n";
        print COM22 "ignore 1:**-${fit_start_152gd_s} 1:${fit_end_152gd_s}-** 2:**-${fit_start_152gd_f} 2:${fit_end_152gd_f}-**\n";
	print COM22 "model gauss + powerlaw\n"; # $fit_model_id = 1
        print COM22 "90.0\n";    # 152Gd slow:initial K alpha line energy in ch
	print COM22 "10.0\n";    # 152Gd slow:initial sigma in ch
	print COM22 "\n";        # 152Gd slow:gauss normalization
	print COM22 "0.9\n";     # 152Gd slow:powerlaw PhoIndex
        print COM22 "\n";        # 152Gd slow:powerlaw normalization
        print COM22 "90.0\n";    # 152Gd fast:initial K alpha line energy in ch
	print COM22 "10.0\n";    # 152Gd fast:initial sigma in ch
	print COM22 "\n";        # 152Gd fast:gauss normalization
	print COM22 "0.9\n";     # 152Gd fast:powerlaw PhoIndex
        print COM22 "\n";        # 152Gd fast:powerlaw normalization
        print COM22 "log ${log_152gd}\n";
	print COM22 "fit 200\n";
	print COM22 "chatter 0 0\nfreeze 2 3 4 5 6 7 9\nchatter 10 10\n";
        print COM22 "chatter 0 0\nfit\nchatter 10 10\n";
        print COM22 "query no\n";
	print COM22 "error 1\n";
        print COM22 "query yes\n";
	print COM22 "chatter 0 0\nfreeze 1\nchatter 10 10\n";
        print COM22 "chatter 0 0\nthaw 6\nchatter 10 10\n";
        print COM22 "chatter 0 0\nfit\nchatter 10 10\n";
        print COM22 "query no\n";
        print COM22 "error 6\n";
        print COM22 "query yes\n";
        print COM22 "chatter 0 0\nthaw 1 2 3 4 5 7 9\nchatter 10 10\n";
        print COM22 "chatter 0 0\nfit\nchatter 10 10\n";
	print COM22 "flux ${fit_start_152gd_s} ${fit_end_152gd_s}\n";
        print COM22 "puts \"exposure_time=[format \"\%E\" [tcloutr expos 1]]\"\n";
	print COM22 "show files\n";
  	print COM22 "setplot add\n";
  	print COM22 "plot ldata delchi\n";
        print COM22 "iplot\n";
	print COM22 "r y 0.001 0.02\n";
        print COM22 "hardcopy ${fit_152gd_ps}/cps\n";
  	print COM22 "plot\n";
	print COM22 "quit\n";
	print COM22 "exit\n";
	print COM22 "\n";
      close COM22;

    ${xspec_152gd} = "xspec < xspec_152gd_${wpu}${wid}_${loop_gti}.xcm";
    print ${xspec_152gd};
    system(${xspec_152gd});

#	unlink COM22;

        #######################################
        # Parse the sdtout to get fit results #
        #######################################

#### Write gti ####
	local($start);
	local($end);
open(RGTI,"${working_dir}/${process_id}_gti_${loop_gti}.gti")
    || die "Cannot open ${process_id}_gti_${loop_gti}.gti: $!\n";
while(<RGTI>){   #read gti log file
    if(/\d+/){
        s/^\s+//g;
        chomp;
        ($start,$end) = split(" +",$_);
	push(@START, $start);
    }
}
	$start=@START[0];
	@RES_GD = ($start, $end, ${unit_id}, ${fit_model_id});
	@RES_511 = ($start, $end, ${unit_id}, ${fit_model_id});
	@RES_152GD = ($start, $end, ${unit_id}, ${fit_model_id});
#	print "=====gti1=====\n@RES \n @START\n";
	print "=====gti1=====\n", (join ' ',@RES) , "\n", (join ' ',@START), "\n";
	undef @START;
        close RGTI;
#	print "=====gti=====\n@RES \n @START\n";
	print "=====gti=====\n", (join ' ',@RES) , "\n", (join ' ',@START), "\n";
#### Write Gd fit result ####

        $dtset_2_flg = 0;

open(FLOG1,"${log_gd}")
    || die "Cannot open ${log_gd}: $!\n";
while(<FLOG1>) {   #read xspec fit log file
    open(RES1,">>${fitresults_gd}");
    if(/LineE/){
        s/^\s+//g;
        chomp;
        ($dum00,$dum01,$dum02,$dum03,$dum04,$dum05,${best_line},$dum06,${line_error}) = split(" +",$_);
       print "=====best line=====\n ${best_line} ${line_error}\n";
	${best_line_16} = ${best_line} * 16;
	${line_error_16} = ${line_error} * 16;
        push(@RES_GD, ${best_line_16}, ${line_error_16});
    }
    elsif(/Sigma/){
        s/^\s+//g;
        chomp;
        ($dum10,$dum11,$dum12,$dum13,$dum14,$dum15,${best_sigma},$dum16,${sigma_error}) = split(" +",$_);
        print "=====best sigma=====\n ${best_sigma} ${sigma_error}\n";
        ${best_sigma_16} = ${best_sigma} * 16;
	${sigma_error_16} = ${sigma_error} * 16;
	push(@RES_GD, ${best_sigma}, ${sigma_error});
    }
    elsif(/gaussian/ && /norm/){
        s/^\s+//g;
        chomp;
        ($dum20,$dum21,$dum22,$dum23,$dum24,${best_norm},$dum25,${norm_error}) = split(" +",$_);
       print "=====norm=====\n ${best_norm} ${norm_error}\n";
        push(@RES_GD, ${best_norm}, ${norm_error});
    }
    elsif(/PhoIndex/){
        s/^\s+//g;
        chomp;
        ($dum30,$dum31,$dum32,$dum33,$dum34,${phoindex},$dum35,${pho_error}) = split(" +",$_);
       print "=====phoindx=====\n${phoindex} ${pho_error}\n";
        push(@RES_GD, ${phoindex}, ${pho_error});
    }
    elsif(/powerlaw/ && /norm/){
        s/^\s+//g;
        chomp;
        ($dum40,$dum41,$dum42,$dum43,$dum44,${pow_norm},$dum45,${pow_norm_error}) = split(" +",$_);
       print "=====pow_norm=====\n ${pow_norm} ${pow_norm_error}\n";
        push(@RES_GD, ${pow_norm}, ${pow_norm_error});
    }
    elsif(/Reduced/ && /chi-squared/){
        s/^\s+//g;
        chomp;
        ($dum50,$dum51,$dum52,$dum52a,${redu_chi},$dum53,$dum54,$dum55,$dum56,$dum57) = split(" +",$_);
       print "=====redu_chi=====\n ${redu_chi}\n";
    }
    elsif(/Model/ && /Flux/ && $dtset_2_flg == 0){
        $dtset_2_flg = 1;
        s/^\s+//g;
        chomp;
        ($dum60,$dum61,$dum62,${area_flux_s},$dum63,$dum64,$dum65,$dum66,$dum67,$dum68,$dum69) = split(" +",$_);
       print "=====area_flux=====\n ${area_flux_s}\n ";
        splice(@RES_GD, 14, 0, ${area_flux_s});
    }
    elsif(/Model/ && /Flux/){
        s/^\s+//g;
        chomp;
       ($dum70,$dum71,$dum72,${area_flux_f},$dum73,$dum74,$dum75,$dum76,$dum77,$dum78,$dum79) = split(" +",$_);
       print "=====area_flux_f=====\n ${area_flux_f}\n ";
        push(@RES_GD, ${area_flux_f});
    }
    elsif(/XSPEC12> error 1/){
	$err_flg=1;
    }
    elsif(/XSPEC12> error 6/){
	$err_flg=6;
    }
    elsif(/XSPEC12>/){
	$err_flg=0;
    }
    elsif($err_flg==1 && /^#\s+1\s+\S+\s+\S+\s+\(\s*(\S+?)\s*,\s*(\S+?)\s*\)/ ){
        ${slow_err_low} = $1;
        ${slow_err_hi}  = $2;
	print "=====slow error=====\n${slow_err_low} ${slow_err_hi}\n";
     }
    elsif($err_flg==6 && /^#\s+6\s+\S+\s+\S+\s+\(\s*(\S+?)\s*,\s*(\S+?)\s*\)/ ){
        ${fast_err_low} = $1;
        ${fast_err_hi}  = $2;
	print "=====fast error=====\n ${fast_err_low} ${fast_err_hi}\n";
    }
    elsif(/^#exposure_time=/){
        s/^\s+//g;
        chomp;
        ($duma0,${exp}) = split("=",$_);
	print "=====exp=====\n ${exp}\n";
    }
}
        
        # these were the default bad errors in xspec11
        if ( !defined $slow_err_low && !defined $slow_err_hi ) {
            $slow_err_low = -90.0000;
            $slow_err_hi  = 999910.;
            print "=====slow error=====\n${slow_err_low} ${slow_err_hi}\n";
        }
        if ( !defined $fast_err_low && !defined $fast_err_hi ) {
            $fast_err_low = -90.0000;
            $fast_err_hi  = 999910.;
            print "=====fast error=====\n ${fast_err_low} ${fast_err_hi}\n";
        }

	${slow_err} = ((${slow_err_hi} - ${slow_err_low})/2) * 16;
	${fast_err} = ((${fast_err_hi} - ${fast_err_low})/2) * 16;
	print "=====slow_fast error=====\n ${slow_err} ${fast_err} \n";

        # undefine the slow/fast error
        $slow_err_low = $slow_err_hi = $fast_err_low = $fast_err_hi = undef;

	splice(@RES_GD, 5, 1, ${slow_err});
	splice(@RES_GD, 16, 1, ${fast_err});
	${fit_start_gd_s_16} = ${fit_start_gd_s} * 16;
	${fit_start_gd_f_16} = ${fit_start_gd_f} * 16;
	${fit_end_gd_s_16} = ${fit_end_gd_s} * 16;
	${fit_end_gd_f_16} = ${fit_end_gd_f} * 16;
	splice(@RES_GD,  4, 0, ${fit_start_gd_s_16}, ${fit_end_gd_s_16});
	splice(@RES_GD, 17, 0, ${fit_start_gd_f_16}, ${fit_end_gd_f_16});
	push(@RES_GD, ${redu_chi});
	splice(@RES_GD,2, 0, $exp);
	
#	print "=====pre RES_GDULT=====\n@RES_GD \n";
	@SLOW_SIGMA = splice(@RES_GD, 9, 2 );
	splice(@RES_GD, 11, 0, @SLOW_SIGMA);
	@FAST_SIGMA = splice(@RES_GD, 22, 2);
	splice(@RES_GD, 24, 0, @FAST_SIGMA);
	splice(@RES_GD, 23, 1, $RES_GD[10]);
	splice(@RES_GD, 29, 1, $RES_GD[16]);
	
	if($RES_GD[7]>$RES_GD[5] && $RES_GD[7]<$RES_GD[6] && $RES_GD[8]<($RES_GD[6]-$RES_GD[5]) && $RES_GD[11]<($RES_GD[7]*0.5/2.35) && $RES_GD[20]>$RES_GD[18] && $RES_GD[20]<$RES_GD[19] && $RES_GD[21]<($RES_GD[19]-$RES_GD[18]) && $RES_GD[24]<($RES_GD[20]*0.5/2.35) && (0 < $RES_GD[31] && $RES_GD[31] < 3.0)){
	    ${quality_flg} = 0;
	}
	else{
	    print "GDULT:Bad Quality\n";
	    ${quality_flg} = 1;
	}

	push(@RES_GD,${quality_flg});

#	print "=====RES_GDULT=====\n@RES_GD \n";
	print "=====RES_GDULT=====\n", (join ' ',@RES_GD), "\n";
        print RES1 join(" ", @RES_GD),"\n";
	close RES1;
	close FLOG1;


#### Write 511keV fit result ####
        $dtset_2_flg = 0;

open(FLOG2,"${log_511}")
    || die "Cannot open ${log_511}: $!\n";
while(<FLOG2>) {   #read xspec fit log file
    open(RES2,">>${fitresults_511}");
    if(/LineE/){
        s/^\s+//g;
        chomp;
        ($dum00,$dum01,$dum02,$dum03,$dum04,$dum05,${best_line},$dum06,${line_error}) = split(" +",$_);
       print "=====best line=====\n ${best_line} ${line_error}\n";
	${best_line_16} = ${best_line} * 16;
	${line_error_16} = ${line_error} * 16;
        push(@RES_511, ${best_line_16}, ${line_error_16});
    }
    elsif(/Sigma/){
        s/^\s+//g;
        chomp;
        ($dum10,$dum11,$dum12,$dum13,$dum14,$dum15,${best_sigma},$dum16,${sigma_error}) = split(" +",$_);
       print "=====best sigma=====\n ${best_sigma} ${sigma_error}\n";
	${best_sigma_16} = ${best_sigma} * 16;
	${sigma_error_16} = ${sigma_error}* 16;
        push(@RES_511, ${best_sigma_16}, ${sigma_error_16});
    }
    elsif(/gaussian/ && /norm/){
        s/^\s+//g;
        chomp;
        ($dum20,$dum21,$dum22,$dum23,$dum24,${best_norm},$dum25,${norm_error}) = split(" +",$_);
       print "=====norm=====\n ${best_norm} ${norm_error}\n";
        push(@RES_511, ${best_norm}, ${norm_error});
    }
    elsif(/PhoIndex/){
        s/^\s+//g;
        chomp;
        ($dum30,$dum31,$dum32,$dum33,$dum34,${phoindex},$dum35,${pho_error}) = split(" +",$_);
       print "=====phoindx=====\n${phoindex} ${pho_error}\n";
        push(@RES_511, ${phoindex}, ${pho_error});
    }
    elsif(/powerlaw/ && /norm/){
        s/^\s+//g;
        chomp;
        ($dum40,$dum41,$dum42,$dum43,$dum44,${pow_norm},$dum45,${pow_norm_error}) = split(" +",$_);
       print "=====pow_norm=====\n ${pow_norm} ${pow_norm_error}\n";
        push(@RES_511, ${pow_norm}, ${pow_norm_error});
    }
    elsif(/Reduced/ && /chi-squared/){
        s/^\s+//g;
        chomp;
        ($dum50,$dum51,$dum52,$dum52a,${redu_chi},$dum53,$dum54,$dum55,$dum56,$dum57) = split(" +",$_);
       print "=====redu_chi=====\n ${redu_chi}\n";
    }
    elsif(/Model/ && /Flux/ && $dtset_2_flg == 0){
        $dtset_2_flg = 1;
        s/^\s+//g;
        chomp;
        ($dum60,$dum61,$dum62,${area_flux_s},$dum63,$dum64,$dum65,$dum66,$dum67,$dum68,$dum69) = split(" +",$_);
       print "=====area_flux=====\n ${area_flux_s}\n ";
        splice(@RES_511, 14, 0, ${area_flux_s});
    }
    elsif(/Model/ && /Flux/){
        s/^\s+//g;
        chomp;
       ($dum70,$dum71,$dum72,${area_flux_f},$dum73,$dum74,$dum75,$dum76,$dum77,$dum78,$dum79) = split(" +",$_);
       print "=====area_flux_f=====\n ${area_flux_f}\n ";
        push(@RES_511, ${area_flux_f});
    }
    elsif(/XSPEC12> error 1/){
	$err_flg=1;
    }
    elsif(/XSPEC12> error 6/){
	$err_flg=6;
    }
    elsif(/XSPEC12>/){
	$err_flg=0;
    }
    elsif($err_flg==1 && /^#\s+1\s+\S+\s+\S+\s+\(\s*(\S+?)\s*,\s*(\S+?)\s*\)/ ){
        ${slow_err_low} = $1;
        ${slow_err_hi}  = $2;
	print "=====slow error=====\n${slow_err_low} ${slow_err_hi}\n";
     }
    elsif($err_flg==6 && /^#\s+6\s+\S+\s+\S+\s+\(\s*(\S+?)\s*,\s*(\S+?)\s*\)/ ){
        ${fast_err_low} = $1;
        ${fast_err_hi}  = $2;
	print "=====fast error=====\n ${fast_err_low} ${fast_err_hi}\n";
    }
    elsif(/^#exposure_time=/){
        s/^\s+//g;
        chomp;
        ($duma0,${exp}) = split("=",$_);
	print "=====exp=====\n ${exp}\n";
    }
}    

        # these were the default bad errors in xspec11
        if ( !defined $slow_err_low && !defined $slow_err_hi ) {
            $slow_err_low = -130.000;
            $slow_err_hi  = 999870.;
            print "=====slow error=====\n${slow_err_low} ${slow_err_hi}\n";
        }
        if ( !defined $fast_err_low && !defined $fast_err_hi ) {
            $fast_err_low = -130.000;
            $fast_err_hi  = 999870.;
            print "=====fast error=====\n ${fast_err_low} ${fast_err_hi}\n";
        }

	${slow_err} = ((${slow_err_hi} - ${slow_err_low})/2) * 16;
	${fast_err} = ((${fast_err_hi} - ${fast_err_low})/2) * 16;
	print "=====slow_fast error=====\n ${slow_err} ${fast_err} \n";

        # undefine the slow/fast error
        $slow_err_low = $slow_err_hi = $fast_err_low = $fast_err_hi = undef;

	splice(@RES_511, 5, 1, ${slow_err});
	splice(@RES_511, 16, 1, ${fast_err});
	${fit_start_511_s_16} = ${fit_start_511_s} * 16;
	${fit_start_511_f_16} = ${fit_start_511_f} * 16;
	${fit_end_511_s_16} = ${fit_end_511_s} * 16;
	${fit_end_511_f_16} = ${fit_end_511_f} * 16;
	splice(@RES_511,  4, 0, ${fit_start_511_s_16}, ${fit_end_511_s_16});
	splice(@RES_511, 17, 0, ${fit_start_511_f_16}, ${fit_end_511_f_16});
	push(@RES_511, ${redu_chi});
	splice(@RES_511,2, 0, $exp);

#	print "=====pre RES_511ULT=====\n@RES_511 \n";
	@SLOW_SIGMA = splice(@RES_511, 9, 2 );
	splice(@RES_511, 11, 0, @SLOW_SIGMA);
	@FAST_SIGMA = splice(@RES_511, 22, 2);
	splice(@RES_511, 24, 0, @FAST_SIGMA);
	splice(@RES_511, 23, 1, $RES_511[10]);
	splice(@RES_511, 29, 1, $RES_511[16]);
	
	if($RES_511[7]>$RES_511[5] && $RES_511[7]<$RES_511[6] && $RES_511[8]<($RES_511[6]-$RES_511[5]) && $RES_511[11]<($RES_511[7]*0.5/2.35) && $RES_511[20]>$RES_511[18] && $RES_511[20]<$RES_511[19] && $RES_511[21]<($RES_511[19]-$RES_511[18]) && $RES_511[24]<($RES_511[20]*0.5/2.35) && (0 < $RES_511[31] && $RES_511[31] < 3.0)){
	    ${quality_flg} = 0;
	}
	else{
	    print "511ULT:Bad Quality\n";
	    ${quality_flg} = 1;
	}

	push(@RES_511,${quality_flg});

#	print "=====RES_511ULT=====\n@RES_511 \n";
	print "=====RES_511ULT=====\n", (join ' ',@RES_511), "\n";
        print RES2 join(" ", @RES_511),"\n";
	close RES2;
	close FLOG2;


#### Write 152Gd fit result ####
        $dtset_2_flg = 0;

open(FLOG3,"${log_152gd}")
    || die "Cannot open ${log_152gd}: $!\n";
while(<FLOG3>) {   #read xspec fit log file
    open(RES3,">>${fitresults_152gd}");
    if(/LineE/){
        s/^\s+//g;
        chomp;
        ($dum00,$dum01,$dum02,$dum03,$dum04,$dum05,${best_line},$dum06,${line_error}) = split(" +",$_);
       print "=====best line=====\n ${best_line} ${line_error}\n";
	${best_line_16} = ${best_line} * 16;
	${line_error_16} = ${line_error} * 16;
        push(@RES_152GD, ${best_line_16}, ${line_error_16});
    }
    elsif(/Sigma/){
        s/^\s+//g;
        chomp;
        ($dum10,$dum11,$dum12,$dum13,$dum14,$dum15,${best_sigma},$dum16,${sigma_error}) = split(" +",$_);
       print "=====best sigma=====\n ${best_sigma} ${sigma_error}\n";
	${best_sigma_16} = ${best_sigma} * 16;
	${sigma_error_16} = ${sigma_error}* 16;
        push(@RES_152GD, ${best_sigma}, ${sigma_error});
    }
    elsif(/gaussian/ && /norm/){
        s/^\s+//g;
        chomp;
        ($dum20,$dum21,$dum22,$dum23,$dum24,${best_norm},$dum25,${norm_error}) = split(" +",$_);
       print "=====norm=====\n ${best_norm} ${norm_error}\n";
        push(@RES_152GD, ${best_norm}, ${norm_error});
    }
    elsif(/PhoIndex/){
        s/^\s+//g;
        chomp;
        ($dum30,$dum31,$dum32,$dum33,$dum34,${phoindex},$dum35,${pho_error}) = split(" +",$_);
       print "=====phoindx=====\n${phoindex} ${pho_error}\n";
        push(@RES_152GD, ${phoindex}, ${pho_error});
    }
    elsif(/powerlaw/ && /norm/){
        s/^\s+//g;
        chomp;
        ($dum40,$dum41,$dum42,$dum43,$dum44,${pow_norm},$dum45,${pow_norm_error}) = split(" +",$_);
       print "=====pow_norm=====\n ${pow_norm} ${pow_norm_error}\n";
        push(@RES_152GD, ${pow_norm}, ${pow_norm_error});
    }
    elsif(/Reduced/ && /chi-squared/){
        s/^\s+//g;
        chomp;
        ($dum50,$dum51,$dum52,$dum52a,${redu_chi},$dum53,$dum54,$dum55,$dum56,$dum57) = split(" +",$_);
       print "=====redu_chi=====\n ${redu_chi}\n";
    }
    elsif(/Model/ && /Flux/ && $dtset_2_flg == 0){
        $dtset_2_flg = 1;
        s/^\s+//g;
        chomp;
        ($dum60,$dum61,$dum62,${area_flux_s},$dum63,$dum64,$dum65,$dum66,$dum67,$dum68,$dum69) = split(" +",$_);
       print "=====area_flux=====\n ${area_flux_s}\n ";
        splice(@RES_152GD, 14, 0, ${area_flux_s});
    }
    elsif(/Model/ && /Flux/){
        s/^\s+//g;
        chomp;
       ($dum70,$dum71,$dum72,${area_flux_f},$dum73,$dum74,$dum75,$dum76,$dum77,$dum78,$dum79) = split(" +",$_);
       print "=====area_flux_f=====\n ${area_flux_f}\n ";
        push(@RES_152GD, ${area_flux_f});
    }
    elsif(/XSPEC12> error 1/){
	$err_flg=1;
    }
    elsif(/XSPEC12> error 6/){
	$err_flg=6;
    }
    elsif(/XSPEC12>/){
	$err_flg=0;
    }
    elsif($err_flg==1 && /^#\s+1\s+\S+\s+\S+\s+\(\s*(\S+?)\s*,\s*(\S+?)\s*\)/ ){
        ${slow_err_low} = $1;
        ${slow_err_hi}  = $2;
	print "=====slow error=====\n${slow_err_low} ${slow_err_hi}\n";
     }
    elsif($err_flg==6 && /^#\s+6\s+\S+\s+\S+\s+\(\s*(\S+?)\s*,\s*(\S+?)\s*\)/ ){
        ${fast_err_low} = $1;
        ${fast_err_hi}  = $2;
	print "=====fast error=====\n ${fast_err_low} ${fast_err_hi}\n";
    }
    elsif(/^#exposure_time=/){
        s/^\s+//g;
        chomp;
        ($duma0,${exp}) = split("=",$_);
	print "=====exp=====\n ${exp}\n";
    }
}

        # these were the default bad errors in xspec11
        if ( !defined $slow_err_low && !defined $slow_err_hi ) {
            $slow_err_low = -90.0000;
            $slow_err_hi  = 999910.;
            print "=====slow error=====\n${slow_err_low} ${slow_err_hi}\n";
        }
        if ( !defined $fast_err_low && !defined $fast_err_hi ) {
            $fast_err_low = -90.0000;
            $fast_err_hi  = 999910.;
            print "=====fast error=====\n ${fast_err_low} ${fast_err_hi}\n";
        }

	${slow_err} = ((${slow_err_hi} - ${slow_err_low})/2) * 16;
	${fast_err} = ((${fast_err_hi} - ${fast_err_low})/2) * 16;
	print "=====slow_fast error=====\n ${slow_err} ${fast_err} \n";

        # undefine the slow/fast error
        $slow_err_low = $slow_err_hi = $fast_err_low = $fast_err_hi = undef;

	splice(@RES_152GD, 5, 1, ${slow_err});
	splice(@RES_152GD, 16, 1, ${fast_err});
	${fit_start_152gd_s_16} = ${fit_start_152gd_s} * 16;
	${fit_start_152gd_f_16} = ${fit_start_152gd_f} * 16;
	${fit_end_152gd_s_16} = ${fit_end_152gd_s} * 16;
	${fit_end_152gd_f_16} = ${fit_end_152gd_f} * 16;
	splice(@RES_152GD,  4, 0, ${fit_start_152gd_s_16}, ${fit_end_152gd_s_16});
	splice(@RES_152GD, 17, 0, ${fit_start_152gd_f_16}, ${fit_end_152gd_f_16});
	push(@RES_152GD, ${redu_chi});
	splice(@RES_152GD,2, 0, $exp);

#	print "=====pre RES_152GDULT=====\n@RES_152GD \n";
	@SLOW_SIGMA = splice(@RES_152GD, 9, 2 );
	splice(@RES_152GD, 11, 0, @SLOW_SIGMA);
	@FAST_SIGMA = splice(@RES_152GD, 22, 2);
	splice(@RES_152GD, 24, 0, @FAST_SIGMA);
	splice(@RES_152GD, 23, 1, $RES_152GD[10]);
	splice(@RES_152GD, 29, 1, $RES_152GD[16]);
	
	if($RES_152GD[7]>$RES_152GD[5] && $RES_152GD[7]<$RES_152GD[6] && $RES_152GD[8]<($RES_152GD[6]-$RES_152GD[5]) && $RES_152GD[11]<($RES_152GD[7]*0.5/2.35) && $RES_152GD[20]>$RES_152GD[18] && $RES_152GD[20]<$RES_152GD[19] && $RES_152GD[21]<($RES_152GD[19]-$RES_152GD[18]) && $RES_152GD[24]<($RES_152GD[20]*0.5/2.35) && (0 < $RES_GD[31] && $RES_GD[31] < 3.0)){
	    ${quality_flg} = 0;
	}
	else{
	    print "152GdULT:Bad Quality\n";
	    ${quality_flg} = 1;
	}

	push(@RES_152GD,${quality_flg});

#	print "=====RES_152GDULT=====\n@RES_152GD \n";
	print "=====RES_152GDULT=====\n", (join ' ',@RES_152GD),"\n";
        print RES3 join(" ", @RES_152GD),"\n";
	close RES3;
    close FLOG3;

    }
}

##########
## MAIN ##
##########
if(defined(@ARGV)) {
    foreach $i (0 .. $#ARGV) {
        $ARGV[$i] = "\"".$ARGV[$i]."\"";
    }
    $invokestring = join(' ',@ARGV);
}

&input($ARGV[0]);  print "FFF = ${fffile}, HK = ${hkfile}\n";

## set file name ##
${working_dir}="${process_id}_hxdmkgainhist_tmp";
${totalgti}="${working_dir}/total.gti";       # original GTI,maketime out put
${fdump_log}="${working_dir}/fdump.log";
${gti_list}="${working_dir}/gtilist";         # List of GTI
${rmf}="${working_dir}/${process_id}dmy.rsp"; # response file

#### data dir ####
$mkworkdir = "rm -rf ${working_dir}; mkdir ${working_dir} \n";
print $mkworkdir;
system($mkworkdir);

&maketime;         print "original GTI = ${totalgti}\n";
&makefdump;        print "fdump log = ${fdump_log}\n";
&makegti;          print "GTI LIST = ${gti_list}\n";
                   print "Exposure = $EXPOSURE\n";
&eventselect;      print "EVENT SELECTION\n";
&makedummyrsp;

${loop_gti}=0;
open(LOOP,"${gti_list}")
    || warn "Cannot open gtilist (${gti_list}): $!\n";
while(<LOOP>) {
    &sort;
    &extractor;
    &rebin;
    &drivexspec;
    ${loop_gti}++;
}
close LOOP;

$rmxspeclog = "rm -rf ${process_id}_xspec*.log xspec*.xcm xautosav.xcm ${rmf} \n";
print $rmxspeclog;
system($rmxspeclog);

if($save_tmpfiles eq 'no') {
    $rmtmpfiles = "rm -rf ${working_dir}\n";
    print $rmtmpfiles;
    system($rmtmpfiles);
}