#include <sstream>
#include <vector>
#include "eventDisplay.h"
using namespace std;

#define PI 3.1415962
#define MAXHITSFDC 2000
#define MAXHITSCDC 2000
#define MAXPOINTS 20000

/////////////////////////
// default constructor //
/////////////////////////
eventDisplay::eventDisplay() : fm("fitter.hddm") {

  // create a canvas
  universe = new TCanvas( "universe", "track", 100, 100, 1000, 800);
  trackTypeColors[0] = kGray;
  trackTypeColors[1] = kRed;
  trackTypeColors[2] = kGreen;
  trackTypeColors[3] = kBlue;
}

#define NPARAMS 5

void eventDisplay::readAndDraw(void){
  int run, event, status;
  bool search = true;
  while (search) {
  status = fm.getNextEvent(run, event);
  if (status) {
    cout << "no more events" << endl;
    return;
  }
  cout << "got next event, run=" << run << " event=" << event << endl;

  eventInfo_t info;
  fm.readInfo(&info);
  cout << "and the number is " << info.event << endl;
  for (int ip = 0; ip < NPARAMS; ip++) {
    cout << "param " << ip << " " << info.paramsTrue[ip] << " " << info.paramsFit[ip] << endl;
  }
  double thetaDeg = info.paramsTrue[2]*180/3.1415926;
  double ptot = 1.0/info.paramsTrue[4]/sin(info.paramsTrue[2]);
  cout << "theta = " << thetaDeg << endl;
  cout << "chisq = " << info.chisq << endl;
  cout << "ptot = " << ptot << endl;
  // if (info.chisq > 100 && thetaDeg > 20. && ptot > 2.0) search = false;
  if (true) search = false;
  }

  // define the geometry

  universe->cd();
  new TGeoManager("ggm", "the global geometry manager"); // define global
						     // geometry
						     // manager
  mat = new TGeoMaterial("Vacuum", 0, 0,0); // define
							  // vacuum
							  // material
  med = new TGeoMedium("Vacuum", 1, mat); // define medium
						      // made of
						      // vacuum
  // define a box, destined to be the top volume
  top = gGeoManager->MakeBox("Top", med, 100.0, 100.0, 500.0);

  // define sphere shape
  sphere = new TGeoSphere(0.0, 0.1, 0.0, 180.0, 0.0, 360.0);
  // define the volume bubble as a sphere made of vacuum
  sphereVacRaw = new TGeoVolume("bubbleRaw", sphere, med);
  sphereVacRaw->SetLineColor(4);
  sphereVacCorr = new TGeoVolume("bubbleCorr", sphere, med);
  sphereVacCorr->SetLineColor(2);
  // define bead
  bead = new TGeoSphere(0.0, .02, 0.0, 180, 0.0, 360);
  beadVacFDC = new TGeoVolume("beadFDC", bead, med);
  beadVacFDC->SetLineColor(1);
  beadVacCDC = new TGeoVolume("beadCDC", bead, med);
  beadVacCDC->SetLineColor(1);
  TGeoTranslation *trans[100]; // define an array of translations for hits
  Double_t x, y, z, theta, phi, dist; // to receive data from file
  Double_t thetadeg, phideg; // to hold angles in degrees
  Int_t index; // ditto
  Double_t xmin = 1.e10; // initialize the bounding parameters
  Double_t xmax = -1.e10;
  Double_t ymin = 1.e10;
  Double_t ymax = -1.e10;
  Double_t zmin = 1.e10;
  Double_t zmax = -1.e10;
  Int_t nhitsFDC = 0; // zero hit counter
  hitInfoFDC_t hitsFDC[MAXHITSFDC];
  // readFDCData(hitsFDC, nhitsFDC);
  fm.readFDCData(hitsFDC, nhitsFDC);
  hitInfoFDC_t* thisHitFDCPtr;
  for (int i = 0; i < nhitsFDC; i++) { // loop through hits
    thisHitFDCPtr = &(hitsFDC[i]);
    index = thisHitFDCPtr->index;
    x = thisHitFDCPtr->x;
    y = thisHitFDCPtr->y;
    z = thisHitFDCPtr->z;
    if (x < xmin) xmin = x; // test for new bounding parameters
    if (x > xmax) xmax = x;
    if (y < ymin) ymin = y;
    if (y > ymax) ymax = y;
    if (z < zmin) zmin = z;
    if (z > zmax) zmax = z;
    // store the translation for this hit
    trans[i] = new TGeoTranslation(x, y, z);
    top->AddNode(sphereVacRaw, i, trans[i]); // add the ball to
						  // the world volume
  }

  createCDCHits();
  createFDCHits();

  Int_t nresidsCDC = 0; // zero hit counter
  residInfoCDC_t residsCDC[MAXHITSCDC];
  Int_t nresidsFDC = 0; // zero hit counter
  residInfoFDC_t residsFDC[MAXHITSCDC];
  fm.readResiduals(residsCDC, nresidsCDC, residsFDC, nresidsFDC);
  createResidualsCDC(residsCDC, nresidsCDC, xmin, xmax, ymin, ymax, zmin, zmax);
  createResidualsFDC(residsFDC, nresidsFDC, xmin, xmax, ymin, ymax, zmin, zmax);

  Double_t rMinDef[3], rMaxDef[3]; // store bounds into an array
  Double_t xbar, dx, ybar, dy, zbar, dz;
  dx = 0.5*(xmax - xmin);
  xbar = 0.5*(xmax + xmin);
  dy = 0.5*(ymax - ymin);
  ybar = 0.5*(ymax + ymin);
  dz = 0.5*(zmax - zmin);
  zbar = 0.5*(zmax + zmin);
  Double_t delta = dx;
  if (dy > delta) delta = dy;
  if (dz > delta) delta = dz;
  rMinDef[0] = xbar - delta - 1.0;
  rMaxDef[0] = xbar + delta + 1.0;
  rMinDef[1] = ybar - delta - 1.0;
  rMaxDef[1] = ybar + delta + 1.0;
  rMinDef[2] = zbar - delta - 1.0;
  rMaxDef[2] = zbar + delta + 1.0;

  for (int i = 0; i < 3; i++) {
    cout << i << " rMinDef=" << rMinDef[i] << " rMaxDef=" << rMaxDef[i] << endl;
  }

  TVirtualGeoTrack* track[3];
  Int_t track_index, track_id;
  for (track_id = 0; track_id < 4; track_id++) {
    track_index = gGeoManager->AddTrack(track_id,0,0); // create a track
    track[track_id] = gGeoManager->GetTrack(track_index); // get its pointer
  }

  //read in the track
  Int_t npoints = 0; // zero point counter
  Int_t id; // track id (starting or final)
  pointInfo_t points[MAXPOINTS];
  //readTrajectoryData(points, npoints);
  fm.readTrajectoryData(points, npoints);
  for (int i = 0; i < npoints; i++) {
    track_id = points[i].trajId;
    index = points[i].index;
    x = points[i].x;
    y = points[i].y;
    z = points[i].z;
    // store the point on the track
    track[track_id]->AddPoint(x, y, z, (Double_t)(index)); // use index as a
						 // surrogate for time
  }

  // define virtual track bounding box for visual effect only
  TGeoBBox *trackBox = new TGeoBBox("trackBox", 0.5*(xmax - xmin), 0.5*(ymax - ymin), 0.5*(zmax - zmin));
  TGeoVolume *trackVol = new TGeoVolume("track_volume", trackBox, med);
  TGeoTranslation *trackTrans = new TGeoTranslation((xmax + xmin)*0.5, (ymax + ymin)*0.5, (zmax + zmin)*0.5);
  //top->AddNode(trackVol, 0, trackTrans);

  gGeoManager->SetTopVolume(top); // make top the top volume
  gGeoManager->CloseGeometry(); // we are finished with creating the geometry
  top->SetLineColor(kMagenta); // set color of the universe
  //gGeoManager->SetTopVisible(); // allow top volume to be seen

  top->Draw(); // draw it
  for (int it = 0; it < 4; it++) {
    track[it]->SetLineColor(trackTypeColors[it]);
    track[it]->Draw();
  }

  // draw some axes
  TAxis3D *axes = new TAxis3D();
  axes->Draw();

  // zoom in on the track
  viewPtr = TView::CreateView(1, rMinDef, rMaxDef);
  //top->Raytrace();

  universe->Update();

}

void eventDisplay::createResidualsCDC(residInfoCDC_t* residsCDC, int nresidsCDC, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, double &zmax) {
  TGeoTranslation transCDC[1024]; // translations for hits
  TGeoRotation rotCDC[1024]; // rotation
  TGeoCombiTrans *combiCDC[1024]; // an array of combo tranlations and rotations
  TGeoTranslation *transCDCTraj[1024]; // trajectory point of closest approach
  residInfoCDC_t *thisResidCDCPtr;
  char strawName[128];
  Int_t index;
  Double_t x, y, z, theta, phi, dist, thetadeg, phideg;
  for (int i = 0; i < nresidsCDC; i++) {
    thisResidCDCPtr = &(residsCDC[i]);
    index = thisResidCDCPtr->index;
    x = thisResidCDCPtr->x;
    y = thisResidCDCPtr->y;
    z = thisResidCDCPtr->z;
    theta = thisResidCDCPtr->theta;
    phi = thisResidCDCPtr->phi;
    dist = thisResidCDCPtr->dist;
    if (x < xmin) xmin = x; // test for new bounding parameters
    if (x > xmax) xmax = x;
    if (y < ymin) ymin = y;
    if (y > ymax) ymax = y;
    if (z < zmin) zmin = z;
    if (z > zmax) zmax = z;
    // store the translation for this hit
    transCDC[i].SetTranslation(x, y, z);
    phideg = phi*180.0/PI + 90;
    thetadeg = theta*180.0/PI;
    rotCDC[i].SetAngles(phideg, thetadeg, 0.0);
    // define straw shape as a cylinder
    TGeoTube *straw = new TGeoTube(0.0, fabs(dist), 0.005);
    // define the straw as a cylinder made of vacuum
    sprintf(strawName, "straw%d",i);
    TGeoVolume *strawVac = new TGeoVolume("straw", straw, med);
    int idTraj = thisResidCDCPtr->idTraj;
    int color;
    if (fabs(thetadeg) < 0.1) { // axial
      color = kMagenta;
    } else { // stereo
      color = kCyan;
    }
    if (idTraj != 1) color -= 5;
    strawVac->SetLineColor(color);
    combiCDC[i] = new TGeoCombiTrans(transCDC[i], rotCDC[i]);
    combiCDC[i]->RegisterYourself();
    top->AddNode(strawVac, i, combiCDC[i]); // add the straw to
						  // the world volume
    transCDCTraj[i] = new TGeoTranslation;
    transCDCTraj[i]->SetTranslation(thisResidCDCPtr->xTraj, thisResidCDCPtr->yTraj, thisResidCDCPtr->zTraj);
    top->AddNode(beadVacCDC, i, transCDCTraj[i]);
  }
}

void eventDisplay::createResidualsFDC(residInfoFDC_t* residsFDC, int nresidsFDC, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, double &zmax) {
  TGeoTranslation *transFDC[1024]; // translations for hits
  TGeoTranslation *transFDCTraj[1024]; // trajectory point of closest approach
  residInfoFDC_t *thisResidFDCPtr;
  char strawName[128];
  Int_t index;
  Double_t x, y, z;
  for (int i = 0; i < nresidsFDC; i++) {
    thisResidFDCPtr = &(residsFDC[i]);
    index = thisResidFDCPtr->index;
    x = thisResidFDCPtr->x;
    y = thisResidFDCPtr->y;
    z = thisResidFDCPtr->z;
    if (x < xmin) xmin = x; // test for new bounding parameters
    if (x > xmax) xmax = x;
    if (y < ymin) ymin = y;
    if (y > ymax) ymax = y;
    if (z < zmin) zmin = z;
    if (z > zmax) zmax = z;
    // store the translation for this hit
    transFDC[i] = new TGeoTranslation;
    transFDC[i]->SetTranslation(x, y, z);
    top->AddNode(sphereVacCorr, i, transFDC[i]);
    transFDCTraj[i] = new TGeoTranslation;
    transFDCTraj[i]->SetTranslation(thisResidFDCPtr->xTraj,
				    thisResidFDCPtr->yTraj,
				    thisResidFDCPtr->zTraj);
    top->AddNode(beadVacFDC, i, transFDCTraj[i]);
  }
}

void eventDisplay::createCDCHits() {
  TGeoTranslation trans[1024]; // translations for hits
  TGeoRotation rot[1024]; // rotation
  TGeoCombiTrans *combi[1024]; // an array of combo tranlations and rotations
  Int_t nCDC = 0; // zero hit counter
  hitInfoCDC_t hits[MAXHITSCDC];
  fm.readCDCData(hits, nCDC);
  cout << "number of CDC hits=" << nCDC << endl;
  hitInfoCDC_t *thisHitPtr;
  char strawName[128];
  Int_t index;
  Double_t x, y, z, theta, phi, L, tdrift, thetadeg, phideg;
  for (int i = 0; i < nCDC; i++) {
    thisHitPtr = &(hits[i]);
    index = thisHitPtr->index;
    x = thisHitPtr->x;
    y = thisHitPtr->y;
    z = thisHitPtr->z;
    theta = thisHitPtr->theta;
    phi = thisHitPtr->phi;
    L = thisHitPtr->L;
    tdrift = thisHitPtr->tdrift;
    // store the translation for this hit
    trans[i].SetTranslation(x, y, z);
    phideg = phi*180.0/PI + 90;
    thetadeg = theta*180.0/PI;
    rot[i].SetAngles(phideg, thetadeg, 0.0);
    // define straw shape as a cylinder
    TGeoTube *straw_again = new TGeoTube(0.0, 0.8, 0.5*L);
    // define the straw as a cylinder made of vacuum
    sprintf(strawName, "straw%d",i);
    TGeoVolume *strawVac_again = new TGeoVolume("straw_again", straw_again,
						med);
    int color = kGreen - 9;
    strawVac_again->SetLineColor(color);
    combi[i] = new TGeoCombiTrans(trans[i], rot[i]);
    combi[i]->RegisterYourself();
    top->AddNode(strawVac_again, i, combi[i]); // add the straw to
						  // the world volume
  }
}

void eventDisplay::createFDCHits() {
  TGeoTranslation trans[1024]; // translations for hits
  TGeoRotation rot[1024]; // rotation
  TGeoCombiTrans *combi[1024]; // an array of combo tranlations and rotations
  Int_t nFDC = 0; // zero hit counter
  hitInfoFDC_t hits[MAXHITSFDC];
  fm.readFDCData(hits, nFDC);
  cout << "number of FDC hits=" << nFDC << endl;
  hitInfoFDC_t *thisHitPtr;
  char wireName[128];
  Int_t index;
  Double_t x, y, z, theta, phi, L, tdrift, thetadeg, phideg;
  for (int i = 0; i < nFDC; i++) {
    thisHitPtr = &(hits[i]);
    index = thisHitPtr->index;
    x = thisHitPtr->xWire;
    y = thisHitPtr->yWire;
    z = thisHitPtr->zWire;
    theta = thisHitPtr->theta;
    phi = thisHitPtr->phi;
    L = thisHitPtr->L;
    // store the translation for this hit
    trans[i].SetTranslation(x, y, z);
    phideg = phi*180.0/PI + 90;
    thetadeg = theta*180.0/PI;
    rot[i].SetAngles(phideg, thetadeg, 0.0);
    // define wire shape as a cylinder
    TGeoTube *wire = new TGeoTube(0.0, 0.01, 0.5*L);
    // define the wire as a cylinder made of vacuum
    sprintf(wireName, "wire%d",i);
    TGeoVolume *wireVac = new TGeoVolume("wire", wire, med);
    int color = kBlack;
    wireVac->SetLineColor(color);
    combi[i] = new TGeoCombiTrans(trans[i], rot[i]);
    combi[i]->RegisterYourself();
    top->AddNode(wireVac, i, combi[i]); // add the straw to
                                        // the world volume
  }
}

void eventDisplay::doZoomIn() {
  viewPtr->Zoom();
  universe->Modified();
  universe->Update();
}

void eventDisplay::doZoomOut() {
  viewPtr->UnZoom();
  universe->Modified();
  universe->Update();
}

void eventDisplay::doPanXP() {
  Double_t rmin[3], rmax[3], fraction = 0.25, distance;
  viewPtr->GetRange(rmin, rmax);
  distance = fraction*(rmax[0] - rmin[0]);
  rmin[0] += distance;
  rmax[0] += distance;
  viewPtr->SetRange(rmin, rmax);
  universe->Modified();
  universe->Update();
  return;
}

void eventDisplay::doPanXM() {
  Double_t rmin[3], rmax[3], fraction = -0.25, distance;
  viewPtr->GetRange(rmin, rmax);
  distance = fraction*(rmax[0] - rmin[0]);
  rmin[0] += distance;
  rmax[0] += distance;
  viewPtr->SetRange(rmin, rmax);
  universe->Modified();
  universe->Update();
  return;
}

void eventDisplay::doPanYP() {
  Double_t rmin[3], rmax[3], fraction = 0.25, distance;
  viewPtr->GetRange(rmin, rmax);
  distance = fraction*(rmax[1] - rmin[1]);
  rmin[1] += distance;
  rmax[1] += distance;
  viewPtr->SetRange(rmin, rmax);
  universe->Modified();
  universe->Update();
  return;
}

void eventDisplay::doPanYM() {
  Double_t rmin[3], rmax[3], fraction = -0.25, distance;
  viewPtr->GetRange(rmin, rmax);
  distance = fraction*(rmax[1] - rmin[1]);
  rmin[1] += distance;
  rmax[1] += distance;
  viewPtr->SetRange(rmin, rmax);
  universe->Modified();
  universe->Update();
  return;
}

void eventDisplay::doPanZP() {
  Double_t rmin[3], rmax[3], fraction = 0.25, distance;
  viewPtr->GetRange(rmin, rmax);
  distance = fraction*(rmax[2] - rmin[2]);
  rmin[2] += distance;
  rmax[2] += distance;
  viewPtr->SetRange(rmin, rmax);
  universe->Modified();
  universe->Update();
  return;
}

void eventDisplay::doPanZM() {
  Double_t rmin[3], rmax[3], fraction = -0.25, distance;
  viewPtr->GetRange(rmin, rmax);
  distance = fraction*(rmax[2] - rmin[2]);
  rmin[2] += distance;
  rmax[2] += distance;
  viewPtr->SetRange(rmin, rmax);
  universe->Modified();
  universe->Update();
  return;
}