TModel *refa;
TModel *sig1a;
TModel *sig2a;
TModel *sig1ph;
TModel *sig2ph;
TCanvas *can1;
TCanvas *can2;

void compare(char* solfile) {
   TBitmap *data = new TBitmap(50,50);
   ifstream in;
   in.open("data.ser");
   in >> data;
   in.close();
   can1 = new TCanvas("can1","graph window",200,10,700,500);
   data->Draw();
   in.open(solfile);
   int ndim;
   double bestCost, nextCost;
   in >> ndim >> bestCost >> nextCost;
   refa = new TModel;
   sig1a = new TModel;
   sig1ph = new TModel;
   sig2ph = new TModel;
   in >> refa >> sig1a >> sig1ph >> sig2ph;
   in.close();
   sig2a = new TModel(*sig1a);
   *sig2a *= 0.83;
   TBitmap *model = new TBitmap(*data);
   model->Precompute(5);
   model->Interfere(*refa,*sig1a,*sig2a,*sig1ph,*sig2ph);
   std::cout << "Solution claims ChiSquare of " << bestCost
             << ", I find " << data->ChiSquared(*model) << std::endl;
   can2 = new TCanvas("can2","graph window",200,10,700,500);
   model->Draw();
}

void compare2(char* sol1, char* sol2, int sel) {
   TBitmap *model1 = new TBitmap(50,50);
   TBitmap *model2 = new TBitmap(50,50);
   refa = new TModel;
   sig1a = new TModel;
   sig2a = new TModel;
   sig1ph = new TModel;
   sig2ph = new TModel;
   ifstream in;
   in.open(sol1);
   int ndim;
   double bestCost, nextCost;
   in >> ndim >> bestCost >> nextCost;
   in >> refa >> sig1a >> sig1ph >> sig2ph;
   in.close();
   *sig2a = *sig1a;
   *sig2a *= 0.83;
   can1 = new TCanvas("can1","graph window",200,10,700,500);
   model1->Precompute(5);
   switch (sel) {
     case 0:
       refa->Print();
       model1->Compute(*refa);
       break;
     case 1:
       sig1a->Print();
       sig2a->Print();
       model1->Compute(*sig1a);
       break;
     case 2:
       sig1ph->Print();
       model1->Compute(*sig1ph);
       break;
     case 3:
       sig2ph->Print();
       model1->Compute(*sig2ph);
       break;
     default:
       model1->Interfere(*refa,*sig1a,*sig2a,*sig1ph,*sig2ph);
   }
   model1->Draw();
   in.open(sol2);
   in >> ndim >> bestCost >> nextCost;
   in >> refa >> sig1a >> sig1ph >> sig2ph;
   in.close();
   *sig2a = *sig1a;
   *sig2a *= 0.83;
   can2 = new TCanvas("can2","graph window",200,10,700,500);
   model2->Precompute(5);
   switch (sel) {
     case 0:
       refa->Print();
       model2->Compute(*refa);
       break;
     case 1:
       sig1a->Print();
       sig2a->Print();
       model2->Compute(*sig1a);
       break;
     case 2:
       sig1ph->Print();
       model2->Compute(*sig1ph);
       break;
     case 3:
       sig2ph->Print();
       model2->Compute(*sig2ph);
       break;
     default:
       model2->Interfere(*refa,*sig1a,*sig2a,*sig1ph,*sig2ph);
   }
   model2->Draw();
}

void plot1(char* sol1, int sel) {
   refa = new TModel;
   sig1a = new TModel;
   sig2a = new TModel;
   sig1ph = new TModel;
   sig2ph = new TModel;
   ifstream in;
   in.open(sol1);
   int ndim;
   double bestCost, nextCost;
   in >> ndim >> bestCost >> nextCost;
   in >> refa >> sig1a >> sig1ph >> sig2ph;
   in.close();
   *sig2a = *sig1a;
   *sig2a *= 0.83;
   can1 = new TCanvas("can1","graph window",200,10,700,500);
   TF2 *fplot = new TF2("fplot",modf,-0.5,0.5,-0.5,0.5,1);
   switch (sel) {
     case 0:
       refa->Print();
       fplot->SetParameter(0,0.);
       break;
     case 1:
       sig1a->Print();
       fplot->SetParameter(0,1.);
       break;
     case 2:
       sig1ph->Print();
       fplot->SetParameter(0,2.);
       break;
     case 3:
       sig2ph->Print();
       fplot->SetParameter(0,3.);
       break;
     default:
       cout << "plot1 error: second argument must be in range 0-3" << endl;
   }
   fplot->Draw("surf4");
}

void plot2(char* sol1, int sel) {
   refa = new TModel;
   sig1a = new TModel;
   sig2a = new TModel;
   sig1ph = new TModel;
   sig2ph = new TModel;
   ifstream in;
   in.open(sol1);
   int ndim;
   double bestCost, nextCost;
   in >> ndim >> bestCost >> nextCost;
   in >> refa >> sig1a >> sig1ph >> sig2ph;
   in.close();
   *sig2a = *sig1a;
   *sig2a *= 0.83;
   TF2 *fplot = new TF2("fplot",modf,-0.5,0.5,-0.5,0.5,1);
   switch (sel) {
     case 0:
       refa->Print();
       fplot->SetParameter(0,0.);
       break;
     case 1:
       sig1a->Print();
       fplot->SetParameter(0,1.);
       break;
     case 2:
       sig1ph->Print();
       fplot->SetParameter(0,2.);
       break;
     case 3:
       sig2ph->Print();
       fplot->SetParameter(0,3.);
       break;
     default:
       cout << "plot1 error: second argument must be in range 0-3" << endl;
   }
   fplot->DrawCopy("surf4same");
}

Double_t modf(Double_t *var, Double_t *par)
{
   switch (par[0]) {
      case 0:
         return refa(var[0],var[1]);
      case 1:
         return sig1a(var[0],var[1]);
      case 2:
         return sig1ph(var[0],var[1])*0.650/(2*TMath::Pi());
      case 3:
         return sig2ph(var[0],var[1])*0.650/(2*TMath::Pi());
      default:
         return 0;
   }
}