// Progran program_2_outer_vertex.cpp - transformation of coordinates
// Created by Richard Hollis
// rholli3@uic.edu
// Summer 2000

////////////////
// This is the seventh transformation program.  Instructions;
//    1).  Find your new data file that contains all the points.
//    2).  Find the Template_2_outer_vertex.out file - this is the template.
//    3).  Create the new file newsurveybov_DATE.out - where DATE is today's date - add this
//         to line 13 of the program.
//    4).  Copy the new survey data into this file following the SAME REGIME as the template.
//    5).  Change the name of the output file to bottomouter1_DATE.out - to distinguish it -
//         line 12.
//    6).  Save, compile and run the program_2_outer_vertex.cpp file in ROOT.
//         ie root[0] .L program_2_outer_vertex.cpp
//            root[1] first()
//    7).  Find the Template_3_outer_vertex.out file - this is the template.
//    8).  Create the new file newsurveybov2_DATE.out - where DATE is today's date - add this
//         to line 13 of the program.
//    9).  Copy the new survey data into this file following the SAME REGIME as the template.
//    10). Change the name of the output file to bottomouter2_DATE.out - to distinguish it -
//         line 12.
//    11). Save, compile and run the .cpp file in ROOT.
//         ie root[0] .L program_2_outer_vertex.cpp
//            root[1] first()
//    12). This is the end of the ROOT programs - return to the InfoFile.txt information
////////////////



#include<iostream>
#include<fstream>
#include<iomanip>

using namespace std;

double modulus(double a[], double b[]);

int first()
{
  gROOT->Reset();
  ofstream outFile2("bottomouter1.out");
  ifstream inFileStart("newsurveybov2_1304.out");
  int number = 1;
  

  ha = new TH1F("ha","This is the distribution for bottom left (x)",100,-600,600);
  hb = new TH1F("hb","This is the distribution for bottom left (y)",100,-600,600);
  hc = new TH1F("hc","This is the distribution for bottom left (z)",100,-600,600);
  hd = new TH1F("hd","This is the distribution for bottom right (x)",100,-600,600);
  he = new TH1F("he","This is the distribution for bottom right (y)",100,-600,600);
  hf = new TH1F("hf","This is the distribution for bottom right (z)",100,-600,600);
  hg = new TH1F("hg","This is the distribution for top left (x)",100,-600,600);
  hh = new TH1F("hh","This is the distribution for top left (y)",100,-600,600);
  hk = new TH1F("hk","This is the distribution for top left (z)",100,-600,600);
  hl = new TH1F("hl","This is the distribution for top right (x)",100,-600,600);
  hm = new TH1F("hm","This is the distribution for top right (y)",100,-600,600);
  hn = new TH1F("hn","This is the distribution for top right (z)",100,-600,600);
      
  for(Int_t loop = 1 ; loop <= number ; loop++)
    {
      gROOT->Reset();
      gRandom->SetSeed();
      Float_t curve[12][3000] = {0};
      double x[4] = {0}, y[4] = {0}, z[4] = {0};
      double dx[4] = {0}, dy[4] = {0}, dz[4] = {0};
      char code[6] = {0}, waste[6] = {0};

      inFileStart >> waste
		  >> x[0] >> y[0] >> z[0]
		  >> dx[0] >> dy[0] >> dz[0];
      inFileStart >> code
		  >> x[1] >> y[1] >> z[1]
		  >> dx[1] >> dy[1] >> dz[1];
      for(Int_t i = 2 ; i <= 2 ; i++)
	{
	  inFileStart >> waste
		      >> x[i] >> y[i] >> z[i]
		      >> dx[i] >> dy[i] >> dz[i];
	};
      
      for ( Int_t i=0; i < 3000; i++)
	{
	  for(Int_t j = 0 ; j <= 2 ; j++)
	    {
	      curve[0+3*j][i]=gRandom->Gaus(x[j],dx[j]);
	      curve[1+3*j][i]=gRandom->Gaus(y[j],dy[j]);
	      curve[2+3*j][i]=gRandom->Gaus(z[j],dz[j]);
	    };
	};
      
      ofstream outFile("junky.out");
      for(Int_t i = 0 ; i < 3000 ; i++)
	{
	  for(Int_t j = 0 ; j <= 2 ; j++)
	    {
	      outFile << curve[0+3*j][i] << " " 
		      << curve[1+3*j][i] << " "
		      << curve[2+3*j][i] << " ";
	      outFile << endl;
	    };
	};
      outFile.close();
      
      cout << "Stage one complete - 3000 points formed"
	   << "\n\tabout a gaussian distribution for "
	   << code << "."
	   << endl;
      
      
      
      double zero[3] = {0}, one[3] = {0}, two[3] = {0};
      double three[3] = {0}, four[3] = {0};
      double v[3] = {0}, v1[3] = {0}, p[3] = {0};
      double modv[4] = {0}, xp[3] = {0}, yp[3] = {0}, zp[3] = {0};
      int padnum = 0, i2 = 0, inum2 = 0, iblcorn = 0, itot = 0;
      ifstream inFile("junky.out");
      ofstream outFile("total3.out");
	
      for(int inum = 1 ; inum <= 3000 ; inum++)
	{
	  inFile >> three[0] >> three[2] >> three[1];
	  inFile >> four[0] >> four[2] >> four[1];
	  inFile >> two[0] >> two[2] >> two[1];
	  
	  modv[1] = modulus(two, three);
	  for(int i = 0 ; i <= 2 ; i++)
	    {
	      zero[i] = 0;
	      xp[i] = (three[i] - two[i]) / modv[1];
	      v[i] = four[i] - two[i];
	    };
	  
	  modv[1] = modulus(v, zero);
	  for(int i = 0 ; i <= 2 ; i++)
	    {
	      v[i] /= modv[1];
	      v1[i]=0;
	    };
	    
	  zp[0] = xp[1]*v[2] - xp[2]*v[1];
	  zp[1] = xp[2]*v[0] - xp[0]*v[2];
	  zp[2] = xp[0]*v[1] - xp[1]*v[0];
	  yp[0] = zp[1]*xp[2] - zp[2]*xp[1];
	  yp[1] = zp[2]*xp[0] - zp[0]*xp[2];
	  yp[2] = zp[0]*xp[1] - zp[1]*xp[0];
	  
	  modv[1] = modulus(xp, zero);
	  modv[2] = modulus(yp, zero);
	  modv[3] = modulus(zp, zero);
	  
	  for(int i = 0 ; i <= 2 ; i++)
	    {
	      xp[i] /= modv[1];
	      yp[i] /= modv[2];
	      zp[i] /= modv[3];
	    };
	  double shift = 0.5;
	  double a[4] = { 1.017 - shift, 
			  12.815 - shift, 
			  36.885 - shift, 
			  61.087 - shift };
	  for(int i = 0 ; i <= 2 ; i++)
	    {
	      v1[i] = a[0]*xp[i] + a[1]*yp[i];
	      p[i]  = two[i] + v1[i];
	    };
	  outFile << setw(10) << p[0] 
		  << setw(10) << p[1] 
		  << setw(10) << p[2]
		  << endl;
	  for(int i = 0 ; i <= 2 ; i++)
	    {
	      v1[i] = a[3]*xp[i] + a[1]*yp[i];
	      p[i]  = two[i] + v1[i];
	    };
	  outFile << setw(10) << p[0] 
		  << setw(10) << p[1] 
		  << setw(10) << p[2]
		  << endl;
	  for(int i = 0 ; i <= 2 ; i++)
	    {
	      v1[i] = a[0]*xp[i] + a[2]*yp[i];
	      p[i]  = two[i] + v1[i];
	    };
	  outFile << setw(10) << p[0] 
		  << setw(10) << p[1] 
		  << setw(10) << p[2]
		  << endl;
	  for(int i = 0 ; i <= 2 ; i++)
	    {
	      v1[i] = a[3]*xp[i] + a[2]*yp[i];
	      p[i]  = two[i] + v1[i];
	    }
	  outFile << setw(10) << p[0] 
		  << setw(10) << p[1] 
		  << setw(10) << p[2]
		  << endl;
	};
      inFile.close();
      outFile.close();
	
	cout << "Stage two complete - points cycled around program."
	     << endl;
	
      ifstream inFile2("total3.out");
      
      Float_t ain, bin, cin, din, ein, fin, gin, hin, kin, lin, min, nin;
      ha->Reset();
      hb->Reset();
      hc->Reset();
      hd->Reset();
      he->Reset();
      hf->Reset();
      hg->Reset();
      hh->Reset();
      hk->Reset();
      hl->Reset();
      hm->Reset();
      hn->Reset();

      for(Int_t i = 1 ; i <= 3000 ; i++)
	{
	  inFile2 >> ain >> bin >> cin;
	  inFile2 >> din >> ein >> fin;
	  inFile2 >> gin >> hin >> kin;
	  inFile2 >> lin >> min >> nin;

	  ha->Fill(ain);
	  hb->Fill(bin);
	  hc->Fill(cin);
	  hd->Fill(din);
	  he->Fill(ein);
	  hf->Fill(fin);
	  hg->Fill(gin);
	  hh->Fill(hin);
	  hk->Fill(kin);
	  hl->Fill(lin);
	  hm->Fill(min);
	  hn->Fill(nin);
	};
      //      gBenchmark->Show("hsimple");
      outFile2 << code << " " << setw(10)
	       << ha->GetMean() << " " << setw(10)
	       << hb->GetMean() << " " << setw(10)
	       << hc->GetMean() << " " << setw(10)
	       << ha->GetRMS() << " " << setw(10)
	       << hb->GetRMS() << " " << setw(10)
	       << hc->GetRMS() << endl;
      outFile2 << code << " " << setw(10)
	       << hd->GetMean() << " " << setw(10)
	       << he->GetMean() << " " << setw(10)
	       << hf->GetMean() << " " << setw(10)
	       << hd->GetRMS() << " " << setw(10)
	       << he->GetRMS() << " " << setw(10)
	       << hf->GetRMS() << endl;
      outFile2 << code << " " << setw(10)
	       << hg->GetMean() << " " << setw(10)
	       << hh->GetMean() << " " << setw(10)
	       << hk->GetMean() << " " << setw(10)
	       << hg->GetRMS() << " " << setw(10)
	       << hh->GetRMS() << " " << setw(10)
	       << hk->GetRMS() << endl;
      outFile2 << code << " " << setw(10)
	       << hl->GetMean() << " " << setw(10)
	       << hm->GetMean() << " " << setw(10)
	       << hn->GetMean() << " " << setw(10)
	       << hl->GetRMS() << " " << setw(10)
	       << hm->GetRMS() << " " << setw(10)
	       << hn->GetRMS() << endl;
      cout << "Stage three complete - errors calculated from distributions."
	   << endl;
    };
  inFile.close();
  inFileStart.close();
}

double modulus(double a[],double b[])
{
  double c[3] = {0};
  for(int i = 0 ; i <= 2 ; i++)
    {
      c[i] = b[i] - a[i];
    };
  double answer = 0;
  answer = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
  
  return answer;
}