#include "analysis.h"
#include "scripts.h"

// ******************************************************
//
//  This program makes C-shell script files which run
//
//       pot_two_exp_fit_run
//
//  on one of the SU(3) data sets for one channel--one level.
//  The location of the refined data files is hardwired into
//  the program.  Also, the file containing the fitting
//  ranges is also hardwired into the program, so these
//  files better exist.
//
//  Input:    (a) one choice of data set
//            (b) one choice for symmetry channel
//            (c) one choice for level in that channel
//
//  All appropriate quark-antiquark orientations in the
//  two_exp_fit_ranges.dat file will be processed.
//
//  Fitting range file specification:
//
//     For each data set, there should be ONE file
//     named "two_exp_fit_ranges.dat"
//
//     Each line in this file should look like
//      channel level     R     t_start  t_stop
//         Sg     1    (3,0,0)     1       11
//
//  If the line contains the symbol '*', that line is ignored.
//  If the line contains the symbol '!', reoptimiziation using 4/0 is performed.
//
//  The following assumptions concerning the quark-antiquark
//  orientations are made:
//              on-axis   (R,0,0)  ->  index R     orientation = 0
//      planar-diagonal   (R,R,0)  ->  index R     orientation = 1
//       cubic-diagonal   (R,R,R)  ->  index R     orientation = 2
//
// *****************************************************


int main(){

charString symchannel,fname,str,prompt,fname2,binfile,emfile,fitrangefile;
int i,j,Rindex,channel,level,tstart,tstop,
    run_choice,orientation=0,plural_fuzzings=0,eff_mass_flag=1;
int maxlevel=4;
int max_line_length=8;
int Nspace_dim=3;
int reoptimize_flag,omit_flag;
int reoptimize_tD=5;
int reoptimize_t0=0;
int block_scale=1;
int em_plot_flag=2;

int Nboot=1024;
read_int("  Nboot? ",Nboot);

run_choice=-1;
while ((run_choice<0)||(run_choice>2)){
   cout << "  Which run?"<<endl;
   cout << "   [0] beta = 3.0, xi=3"<<endl;
   cout << "   [1] beta = 2.5, xi=5"<<endl;
   cout << "   [2] beta = 2.2, xi=5"<<endl;
   read_int("  Run? ",run_choice,max_line_length);
   }

if (run_choice==0){
   binfile="/home1/REFINED_DATA/su3_lgt4/B300_AR3/B30_AR3_REFBINS";
   emfile="/home1/REFINED_DATA/su3_lgt4/B300_AR3/B30_AR3_effmass";
   fitrangefile="/home1/STATIC_POTENTIAL_RESULTS/SU3_4D/B30_AR3/two_exp_fit_ranges.dat";
   plural_fuzzings=0;
   eff_mass_flag=1;
   orientation=0;}
else if (run_choice==1){
   binfile="/home1/REFINED_DATA/su3_lgt4/B250_AR5/B25_AR5_REFBINS";
   emfile="/home1/REFINED_DATA/su3_lgt4/B250_AR5/B25_AR5_effmass";
   fitrangefile="/home1/STATIC_POTENTIAL_RESULTS/SU3_4D/B25_AR5/two_exp_fit_ranges.dat";
   plural_fuzzings=0;
   eff_mass_flag=1;
   orientation=0;}
else if (run_choice==2){
   binfile="/home1/REFINED_DATA/su3_lgt4/B220_AR5/B22_AR5_REFBINS";
   emfile="/home1/REFINED_DATA/su3_lgt4/B220_AR5/B22_AR5_effmass";
   fitrangefile="/home1/STATIC_POTENTIAL_RESULTS/SU3_4D/B22_AR5/two_exp_fit_ranges.dat";
   plural_fuzzings=0;
   eff_mass_flag=1;
   orientation=0;}

   
 if (orientation==0){  // on-axis
    Nirreps=10;
    labels.reset(Nirreps);
    labels(1) = "A1g";
    labels(2) = "A1u";
    labels(3) = "A2g";
    labels(4) = "A2u";
    labels(5) = "B1g";
    labels(6) = "B1u";
    labels(7) = "B2g";
    labels(8) = "B2u";
    labels(9) = "Eg";
    labels(10)= "Eu";}
 else if (orientation==2){  // cubic-diagonal
    Nirreps=6;
    labels.reset(Nirreps);
    labels(1) = "A1g";
    labels(2) = "A1u";
    labels(3) = "A2g";
    labels(4) = "A2u";
    labels(5) = "Eg";
    labels(6) = "Eu";}
 else if (orientation==1){  // planar-diagonal
    Nirreps=8;
    labels.reset(Nirreps);
    labels(1) = "A1g";
    labels(2) = "A1u";
    labels(3) = "A2g";
    labels(4) = "A2u";
    labels(5) = "B1g";
    labels(6) = "B1u";
    labels(7) = "B2g";
    labels(8) = "B2u";}

channel=0;
prompt="    Channel to analyze? [1-10] ";
while ((channel<1) || (channel>Nirreps)){
   cout << endl<<"                      CHANNELS:"<<endl<<endl;
   j=0; i=0;
   while (j<Nirreps){
      j++; i++;
      cout <<"     ["<<int_to_string(j,2)<<"]  "<<labels(j);
      if (!(i%4)){
         cout << endl; i=0;}
      }
   cout <<endl<<endl;
   channel=0;
   read_int(prompt,channel,max_line_length);
   if ((channel<1) || (channel>Nirreps)){
      cout << endl<<endl<<" Value 1-10 expected! "<<endl;
      read_string("  Hit <return> to continue .... ",str);
      channel=0;}
   }
level=0;
while ((level<1)||(level>maxlevel)){
   read_int("  Choose which operator? [1-"
      +int_to_string(maxlevel)+"]",level,max_line_length);
   if ((level<1)||(level>maxlevel))
      cout << "Bad value for level"<<endl;
   }

symchannel=labels(channel);
for (j=2;j<=level;j++) symchannel+="p";
if (orientation==1) symchannel=symchannel+"_PD";
else if (orientation==2) symchannel=symchannel+"_CD";

ifstream in(convert_to_charstar(fitrangefile));
if (!in){
   cerr << " could not open the fitting range file <" << fitrangefile <<"> "<<endl;
   exit(1);
   }


fname="fit_"+symchannel+".csh";
cout << " Creating file: "<<fname<<endl;
ofstream out(convert_to_charstar(fname));

out << "#!/bin/csh"<<endl<<endl;
out << "rm -f "<<symchannel<<"_R*.gif"<<endl;
out << "rm -f "<<symchannel<<"_R*.ps"<<endl;
out << "rm -f "<<symchannel<<"_R*.sm"<<endl;
out << "/home/colin/research/analysis/RUN/pot_two_exp_fit_run << END"<<endl;
out << binfile <<endl;
out << "fits_"+symchannel<<".log"<<endl;
if (eff_mass_flag){
   out << "y"<<endl;
   out << emfile<<endl;
   }
else
   out << "n"<<endl;
out << orientation<<endl;
out << channel <<endl;
out << level << endl;
if (plural_fuzzings){
   out << "3"<<endl;
   out << "n"<<endl;
   out << plural_fuzzings<<endl;
   }

while (!in.eof()){ 
   if (parse_next_line(in,channel,level,orientation,Nspace_dim,Rindex,
                       tstart,tstop,reoptimize_flag,omit_flag,block_scale)){
      add_node(Rindex);
      out << "2"<<endl;       // choose R value
      out << Rindex << endl;
      if (!reoptimize_flag){
         out << "4"<<endl;       // read the data
         out << "1"<<endl;       // rebin data
         out << endl;
         if (eff_mass_flag)
            out << em_plot_flag<<endl;       // eff_mass  flag
         out << "7"<<endl;       // do 2-exponential fit
         out << tstart<<endl;
         out << tstop<<endl;
         out << "1"<<endl;       // automatic plot limits
         out << Nboot<<endl;     // number of bootstraps
         if (omit_flag) out << "y" <<endl;
         else out << "n" << endl;
         out << "n"<<endl;       // tweak plot?
         out << "y"<<endl;       // save plot?
         out << "4"<<endl;       // output choice: .sm and .ps
         out << symchannel<<"_R"<<int_to_string(Rindex*block_scale)<<endl;
         out << "n"<<endl;       // re-tweak?
         out << endl;
         }
      else{
         out << "10"<<endl;       // read the data and reoptimize
         out << reoptimize_tD <<endl;
         out << reoptimize_t0 <<endl;
         out << "1"<<endl;       // rebin data
         out << endl;
         if (eff_mass_flag)
            out << em_plot_flag<<endl;       // eff_mass  flag
         out << "7"<<endl;       // do 2-exponential fit
         out << tstart<<endl;
         out << tstop<<endl;
         out << "1"<<endl;       // automatic plot limits
         out << Nboot<<endl;     // number of bootstraps
         if (omit_flag) out << "y" <<endl;
         else out << "n" << endl;
         out << "n"<<endl;       // tweak plot?
         out << "y"<<endl;       // save plot?
         out << "4"<<endl;       // output choice: .sm and .ps
         out << symchannel<<"_R"<<int_to_string(Rindex*block_scale)<<endl;
         out << "n"<<endl;       // re-tweak?
         out << endl;
         }
      }
   }

in.close();
in.clear();

out << "9"<<endl;
out << "END"<<endl;
fname2=symchannel+"_res.dat";

node *pnode=head;
while (pnode!=NULL){
   out << "/usr/local/bin/ps2gif "<<symchannel
       <<"_R"<<int_to_string(pnode->value*block_scale)<<endl;
   pnode=pnode->next;
   }
if (head!=NULL){
   out << "grep \"R=\" "<<symchannel<<"_R"
       << int_to_string(head->value*block_scale)<<".sm >! "<<fname2<<endl;
   out << "grep \"E=\" "<<symchannel<<"_R"
       << int_to_string(head->value*block_scale)<<".sm >> "<<fname2<<endl;
   out << "grep \"Q=\" "<<symchannel<<"_R"
       << int_to_string(head->value*block_scale)<<".sm >> "<<fname2<<endl;
   pnode=head->next;
   while (pnode!=NULL){
      out << "grep \"R=\" "<<symchannel<<"_R"
          << int_to_string(pnode->value*block_scale)<<".sm >> "<<fname2<<endl;   
      out << "grep \"E=\" "<<symchannel<<"_R"
          << int_to_string(pnode->value*block_scale)<<".sm >> "<<fname2<<endl;   
      out << "grep \"Q=\" "<<symchannel<<"_R"
          << int_to_string(pnode->value*block_scale)<<".sm >> "<<fname2<<endl;
      pnode=pnode->next;
      }
   out << "echo \" \" >> "<<fname2<<endl;
   out << "grep \"*** Final fit result\" fits_"
       << symchannel<<".log >> "<<fname2<<endl;
   }
out.close();
str="chmod +x "+fname;
system(convert_to_charstar(str)); // makes script file executable

delete_nodes();

return 0;
}