/*

Standalone program that converts the mechanical database files
into a partial alignment

The location of the directories which contain the file 
are hardwired below.

Uses hardware layer numbering (0-4 for SVX-II)
The inputs are in various coordinate systems and mm 
but are converted on input to cm and the following coordinate system.
X points radially outward, y ponts along CDF phi, and z along CDF z.
The wafers are numbered 0-3 increasing along CDF z.

In the last step, where the numbers are inserted into an AlHolder
before that are written out, the odd number wedges have to be 
rotated 180deg around the radial vector.  This also effect the wafer 
numbering for those half-ladders.  This step bring it into the CDF
offline coordinate system.

Ray Culbertson 4/2001


7/16/2002 - added barrel ordering, Ritchie 1 2 3 -> offline position 1 2 0
and changed flipped barrel from offline #1 to #2
7/19/2002 - switch to Andy's electronic fingerprint method of 
finding where ladders are on the detector (ladderLocations.txt)
 */


#include <math.h>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <string>
#include <cstring>
#include <vector>
#include "TROOT.h"
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TMinuit.h"
#include "TH1.h"
#include "TString.h"
#include "TGraph.h"

#include "Alignment/AlHolder.hh"
#include "Alignment/AlDigiCode.hh"

extern void fcn(int& npar, double *gin, double& f, double* x, int iflag);

//this global vars contain the data for the Minuit fits
double xWarpWaf[200];
double yWarpWaf[200];
double zWarpWaf[200];
int useWarp[200];
int fitWafer;

using namespace std;

class Ladder {
public:
  int barrel;
  int wedge;
  int layer;
  bool flip;
  char name[5];
  int number;
  double width;  // ladder width
  // the transverse measurements, two positions for each wafer
  double yy[8];   // transverse measurements of the wafers
  double zz[8];
  double yygap[8];   // z measurements of the wafers on each
  double zzgap[8];   //    side of gap - not sure what these are
  double gap[3];   // nom gap between wafers - from CMM
  double sd[3];   // delta z of the wafers
  double yoff[4]; // output alignment of the wafers
  double xrot[4];
  double zoff[4];
  // from the warp measurements, the radial offset
  double xoff[4];
  double yrot[4];
  double zrot[4];
  // the wafer positions from the jig, 
  // before correcting to the barrel measurement
  double xoffJig[4];
  double yrotJig[4];

  // the wafer warp measurements
  int nWarp;
  double xWarp[200];
  double yWarp[200];
  double zWarp[200];
  double fitWarp[4][6]; //4=nwaf, 6 = "px","pay","paz","pyy","pyz","pzz"
  double chiWarp[4];
  double hiWarp[4];
  double betayy[4];
  double betayz[4];
  double betazz[4];

  // ladder bow, from fit wafer placement
  double bow;
  // from three radial measurements, center and ends
  double rawBow;
  double xofflad; // these are not used, they fall out of fit to radial
  double yrotlad; // measurements

  // measurements done while placing the ladder on the barrel
  //input numbers
  double barAx;  // rotation around radial vector
  double barDz;  // 
  double barDy;  // along +phi, opposite form file
  double barBow;
  double barRwaf;
  double barAy;  // unfilled
  double barAz;  // unfilled

  TGraph* gRaw;
  TGraph* gFit;

  Ladder(){};
  Ladder(const char* namea, int b, int w, int l, bool fl);
  void findXY();
  void findQuadratic();
  void adjustBow();
  void print(std::ostream& os);
  //  void store(AlHolder& al);
};

class FileReader {
public:
  FileReader(){};
  void setDir(const char* dira);
  void readBarrels(Ladder** det);
  char dir[200];
  void readFileXY(Ladder** det);         // transverse jig measurements
  void readFileQuadratic(Ladder** det);  // radial jig measurements
  void readFileBarrel(Ladder** det);     // measuremnts on the barrel
};

class Plotter {
public:
  Plotter();
  void makeTuple(Ladder& lad);
  void writeTuple();
  void fillAlignment(Ladder** det);

  Ladder plotLadder;
  TFile* hfile;
  TTree* trTree;

  TH1D *hdy,*hdyl[5];
  TH1D *hax,*haxl[5];
  TH1D *hdz,*hdzl[5];
  TH1D *hbo,*hbol[5];
  TH1D *hbochi,*hbochil[5];
};

       ////////////////////////////////////////// Main
int main(int argv, char* argc[]){

  TROOT simple("simple","Translate mechanical alignment");

  Ladder* det[5*12*3];
  for(int i=0; i<5*12*3; i++) det[i]=NULL;

  FileReader fr;

  // information of which ladders are on which barrels
  fr.setDir("mech/asbuilt");
  fr.readBarrels(det);

  for (int i=0; i<3*12*5; i++) {
    det[i]->print(std::cout);
  }

  // measurements of the transverse position of wafer on ladder
  fr.setDir("mech/xy");
  fr.readFileXY(det);
  // measurements of the radial position of wafer on ladder
  fr.setDir("mech/bow");
  fr.readFileQuadratic(det);
  // measurements of the position of the ladder on the barrel
  fr.setDir("mech/barrel");
  fr.readFileBarrel(det);

  Plotter* pPlo = new Plotter();

  for (int i=0; i<3*12*5; i++) {
  //for (int i=0; i<20; i++) {
    det[i]->findXY();
    det[i]->findQuadratic();
    pPlo->makeTuple(*det[i]);
  }

  pPlo->writeTuple();

  // now scale the jig bow shape to the barrel bow measurement
  for (int i=0; i<3*12*5; i++) {
    det[i]->adjustBow();
  }

  pPlo->fillAlignment(det);

  return 0;
}

                 ///////////// Ladder ///////////////////

Ladder::Ladder(const char* namea, int b, int w, int l, bool fl):
   barrel(b),wedge(w),layer(l),flip(fl)
{
  strncpy(name,namea,4);
  number = atoi(name+1);
  nWarp = 0;
  for(int i=0; i<8; i++) {
    yy[i]=-1001.0;
    zz[i]=-1001.0;
    yygap[i]=-1001.0;
    zzgap[i]=-1001.0;
  }
  for(int i=0; i<3; i++) {
    sd[i]=-1001.0;
    gap[i]=-1001.0;
  }
  if(l==0) width = 1.7140;
  if(l==1) width = 2.5594;
  if(l==2) width = 4.0300;
  if(l==3) width = 4.7860;
  if(l==4) width = 6.0170;
  barRwaf = -1001.0;
  bow=0;
  rawBow=0;
  gRaw=NULL;
  gFit=NULL;
}

void Ladder::findXY(){

  bool bad = false;
  bool prt = false;
  for(int i=0; i<8; i+=2) {
    if(yy[i]<-1000.0 || yy[i]<-1000.0) bad = true;
    yoff[i/2] = (yy[i] + yy[i+1])/2.0;
    // x rotation sense is the opposite of the slope in the y-z plane
    xrot[i/2] = -( (yy[i+1] - yy[i])/(zz[i+1]-zz[i]) );
  }

  //bad means there weren't a complete set of measurements
  if(bad) {
    for(int i=0; i<4; i++) {
      yoff[i] = 0;
      xrot[i] = 0;
    }
  }

  if(sd[0]<-1000.0 || sd[1]<-1000.0 || sd[2]<-1000.0) bad = true;
  zoff[0] = -sd[0] - sd[1]/2.0;
  zoff[1] = -sd[1]/2.0;
  zoff[2] =  sd[1]/2.0;
  zoff[3] =  sd[2] + sd[1]/2.0;
  if(bad) {
    for(int i=0; i<4; i++) {
      zoff[i] = 0;
    }
  }
  if(fabs(sd[0])>0.005 || fabs(sd[1])>0.005 || fabs(sd[2])>0.005  ) prt = true;
  if(prt) {
    std::cout << "found large delta Z " << std::endl;
    print(std::cout);
    for(int i=0;i<3;i++){
      std::cout << i << " " << sd[i] << std::endl;
    }
  }

}

void Ladder::findQuadratic(){

  float xg[1000];
  float yg[1000];
  char gname[100];

  for(int i=0; i<200; i++) { useWarp[i]=-1; }

  // assign points to each wafer
  for(int waf=0; waf<4; waf++) {
    double zmed = (-1.5 + waf)*7.43;
    double zlow = zmed - 7.43/2.0;
    double zhgh = zmed + 7.43/2.0;
    for(int i=0; i<nWarp; i++) {
      if(zWarp[i]>zlow && zWarp[i]<zhgh) {
	useWarp[i]=waf;
	xWarpWaf[i] = xWarp[i];
	yWarpWaf[i] = yWarp[i];
	zWarpWaf[i] = zWarp[i] - zmed;
      }
    }
  }  // end loop over wafers

  // find the simple bow from three extreme points
  float dy,dz,d;
  float dlowmn =99.0,dhghmn=99.0,dmedmn=99.0;
  int ilow=-1,imed=-1,ihgh=-1;

  for(int i=0; i<nWarp; i++) {
    dy = 0-yWarp[i];    dz = -14.0-zWarp[i];
    d = sqrt(dy*dy+dz*dz);
    if(d<dlowmn) { dlowmn=d; ilow=i; }
    dy = 0-yWarp[i];    dz = 0-zWarp[i];
    d = sqrt(dy*dy+dz*dz);
    if(d<dmedmn) { dmedmn=d; imed=i; }
    dy = 0-yWarp[i];    dz = 14.0-zWarp[i];
    d = sqrt(dy*dy+dz*dz);
    if(d<dhghmn) { dhghmn=d; ihgh=i; }
  }  // end loop over wafers
  if(ilow>-1 && imed>-1 && ihgh>-1) {
    //rawBow = xmed - (xlow+xhgh)/2.0;
    rawBow = xWarp[imed] - (xWarp[ilow]+xWarp[ihgh])/2.0;
  }

  // fill the TGraph for the raw measurements
  int ng=0;
  for(int i=0; i<nWarp; i++) {
    if(    ( layer<4  && fabs(yWarp[i])<0.4)
	|| ( layer==4 && fabs(yWarp[i]+0.8)<0.4)  ) {
      xg[ng]=zWarp[i];
      yg[ng]=xWarp[i];
      ng++;
    }
  }  // end loop over wafers
  if(ilow>-1 && ihgh>-1) {
    xg[ng] = zWarp[ilow];
    yg[ng] = xWarp[ilow];
    ng++;
    xg[ng] = zWarp[ihgh];
    yg[ng] = xWarp[ihgh];
    ng++;
  }
  if(ng>0) {
    gRaw = new TGraph(ng,xg,yg);
    sprintf(gname,"gRaw_%d",number);
    gRaw->SetName(gname);
    gRaw->SetMarkerStyle(8); gRaw->SetMarkerSize(0.8);
  }


  // start fit
  
  const char* pname[6] = {"px","pay","paz","pyy","pyz","pzz"};
  double strt[6] = {0,0,0,0,0,0};
  double step[6] = {0.0001,0.0001,0.0001,0.0001,0.0001,0.001};
  double lowb[6] = {0,0,0,0,0,0};
  double hghb[6] = {0,0,0,0,0,0};
  int ierflg;

  TMinuit *min = new TMinuit(6);
  min->mninit(5,6,7);
  min->SetFCN(fcn);

  for(int waf=0; waf<4; waf++) {

    int ndof=0;
    for(int i=0; i<nWarp; i++) { if(useWarp[i]==waf) ndof++; };
    // check that ere are enough points at the end of the wafer
    // so that warp can be reliably fit

    float zlocmin =  99.0;
    float zlocmax = -99.0;
    for(int i=0; i<nWarp; i++) { 
      if(useWarp[i]==waf) {
	if(zWarpWaf[i]<zlocmin) zlocmin = zWarpWaf[i];
	if(zWarpWaf[i]>zlocmax) zlocmax = zWarpWaf[i];
      }
    };

    if(ndof>6) {

      std::cout <<"Starting fit " << waf << " " ;
      print(cout); std::cout << std::endl;
      // pass to global variable for fcn
      fitWafer = waf;

      // external par number, name, start, step size, low bound, high bound
      for (int i=0; i<6; i++) {
	min->mnparm(i,pname[i],strt[i],step[i],lowb[i],hghb[i],ierflg);
      }
      // if there aren't points with large lever arm don't fit quadratic
      if(zlocmin>-2.0 || zlocmax<2.0) {
	std::cout << " turning off quadratic terms " << std::endl;
	double p4=4;  double p5=5;  double p6=6;
	min->mnexcm("FIX", &p4 ,1,ierflg); 
	min->mnexcm("FIX", &p5 ,1,ierflg);
	min->mnexcm("FIX", &p6 ,1,ierflg);
      }

      double p0=0;
      min->mnexcm("MIGRAD", &p0 ,0, ierflg);

      if(zlocmin>-2.0 || zlocmax<2.0) {

	std::cout << " turning on quadratic terms  " << std::endl;
	double p4=4;  double p5=5;  double p6=6;
	min->mnexcm("RELEASE", &p4 ,1,ierflg); 
	min->mnexcm("RELEASE", &p5 ,1,ierflg);
	min->mnexcm("RELEASE", &p6 ,1,ierflg);
      }
  
      double val[6], err, xlolim, xuplim;
      int iuint;

      for(int i=0; i<6; i++){
	TString ts(pname[i]);
	min->mnpout(i, ts, val[i], err, xlolim, xuplim, iuint); 
	//   std::cout<< val[i] << std::endl;
	fitWarp[waf][i] = val[i];
      }
      double f;
      double* gin =0;
      int npar =6;
      fcn(npar, gin, f, val, ierflg);
      // convert to chi**2/dof
      chiWarp[waf] = f/(ndof-6);

      xoff[waf] = val[0];
      yrot[waf] = val[1];
      zrot[waf] = val[2];
      betayy[waf]=val[3];
      betayz[waf]=val[4];
      betazz[waf]=val[5];

      hiWarp[waf]=0.0;
      for(double y0=-width/2.0; y0<width/2.0+0.1; y0+=width/2.0) {
	for(double z0=-7.43/2.0; z0<7.43/2.0+0.1; z0+=7.43/2.0) {
	  // hitemp sign is changed so the result is the height of 
	  // the center of the wafer compared to the edges
	  double hitemp = -( y0*y0*val[3] + y0*z0*val[4] + z0*z0*val[5]);
	  //if(fabs(hitemp) > fabs(hiWarp[waf])) hiWarp[waf] = hitemp;
	  hiWarp[waf] += hitemp;
	}
      }
      hiWarp[waf] = hiWarp[waf]/9.0;

      if(fabs(hiWarp[waf])>0.015 || chiWarp[waf]>10.0) {
	std::cout << "found large warp or chi2: " << hiWarp[waf] << " " 
	     << chiWarp[waf] << std::endl;
	for(int i=0; i<nWarp; i++) {
	  if(useWarp[i]==waf) {
	    std::cout << i << " " << xWarpWaf[i] << " " << 
	      yWarpWaf[i] << " " << zWarpWaf[i] << " " 
		 << zWarp[i] -(49.7 - 2*74.3)/10.0<< " " 
		 << xWarpWaf[i]/0.0005 << std::endl;
	  }
	}
	print(std::cout);
      }
    } else {  // end if has some data available
      
      std::cout << "setting a warp to zero " << waf << " " << nWarp<< std::endl;
      print(std::cout);

	for(int i=0; i<nWarp; i++) {
	  //if(useWarp[i]==waf) {
	    std::cout << i << " " << xWarpWaf[i] << " " << 
	      yWarpWaf[i] << " " << zWarpWaf[i] << " " 
		 << zWarp[i] -(49.7 - 2*74.3)/10.0<< " " 
		 << xWarpWaf[i]/0.0005 << " " <<useWarp[i] << std::endl;
	    //}
	}


      for(int i=0; i<6; i++){ fitWarp[waf][i] = 0.0; }
      xoff[waf] = 0;
      yrot[waf] = 0;
      zrot[waf] = 0;
      hiWarp[waf]=0.0;
      chiWarp[waf]=-1.0;
      bow=0.0;
    }  // end if section when no data avail

  }  // end loop over wafers

  // now correct for the overall ladder position in the jig
  double zpin = 15.16;  // pin positions
  double xll = xoff[0] + (1.5*7.43-zpin)*yrot[0];
  double xlh = xoff[3] + (zpin-1.5*7.43)*yrot[3];
  xofflad = (xlh + xll)/2.0;
  yrotlad = (xlh-xll)/(2*zpin);

  //correct the individual wafers
  for(int waf=0; waf<4; waf++) {
    xoff[waf] = xoff[waf] - ((-1.5+waf)*7.43*yrotlad + xofflad);
    yrot[waf] = yrot[waf] - yrotlad;
    // save these, they will be overwritten later
    xoffJig[waf] = xoff[waf];
    yrotJig[waf] = yrot[waf];
  }

  // now fill in bow, center height - pin
  bow = (xoff[1] + (7.43/2.0)*yrot[1] + (7.43/2.0)*(7.43/2.0)*betazz[1]);

  // fill the TGraph for the fit measurements
  ng=0;
  for(int i=0; i<201; i++) {
    xg[ng]=((i-100.0)/100.0)*2*7.43;
    int iwaf = i/50; if(iwaf>3) iwaf=3;
    float xloc = xg[ng] - (iwaf-1.5)*7.43;
    yg[ng]=(xoff[iwaf]+yrot[iwaf]*xloc+betazz[iwaf]*xloc*xloc
		   + xg[ng]*yrotlad + xofflad );
    ng++;
  }  // end loop over wafers
  if(ng>0) {
    gFit = new TGraph(ng,xg,yg);
    sprintf(gname,"gFit_%d",number);
    gFit->SetName(gname);
  }

  /*
  // this finds the largest deviation - can't be compared to barrel
  //measurements which select the center height
  double bowtemp;
  bow=0.0;
  for(int waf=0; waf<4; waf++) {
    bowtemp = (xoff[waf] + (7.43/2.0)*yrot[waf]);
    if(fabs(bowtemp)>fabs(bow)) bow = bowtemp;
    bowtemp = (xoff[waf] - (7.43/2.0)*yrot[waf]);
    if(fabs(bowtemp)>fabs(bow)) bow = bowtemp;
  }
  */

}

void Ladder::print(std::ostream& os){
  os << name << " l,w,b,flip: " << layer << setw(3) << wedge 
     << " " << barrel << "   " << flip 
     << "    warps: " << nWarp 
     << std::endl;
}

//void Ladder::store(AlHolder& al){}


        ////////////////////  fcn   ///////////////////
void fcn(int& npar, double *gin, double& f, double* xa, int iflag)
{
  f=0;
  int n=0;
  for(int i=0; i<200; i++){
    if(useWarp[i] == fitWafer) {
      n++;
      double x = xWarpWaf[i];
      double y = yWarpWaf[i];
      double z = zWarpWaf[i];
      //4/19/01 fixed bug by y->-z, z->y
      // a positve rotation of wafer about z axis causes x to decrease with y
      double xpred = xa[0] + xa[1]*z + xa[2]*(-y) 
	          + xa[3]*y*y + xa[4]*y*z +xa[5]*z*z;
      double xpara = x;
      f += (xpred-xpara)*(xpred-xpara);
    }
  }
  // approximate uncertainty in radial wafer measurements is 5um
  f=f/(0.0003*0.0003);

}

void Ladder::adjustBow() {
  // correct the jig bow measurement to the barrel overall scale
  // this will be the final position
  float d = barBow - bow;  //This is the difference in the two bows
  float x1, x2;
  for(int w=0; w<4; w++) {
    x1 = xoff[w] - yrot[w]*(7.43/2.0);
    x2 = xoff[w] + yrot[w]*(7.43/2.0);
    if(w<=1) {
      x1 = x1 + 0.5*d*w;
      x2 = x2 + 0.5*d*(w+1);
    } else {
      x1 = x1 + 0.5*d*(4-w);
      x2 = x2 + 0.5*d*(3-w);
    }
    xoff[w] = (x1+x2)/2.0;
    yrot[w] = (x2-x1)/7.43;
  }
  //float nw = xoff[1]+yrot[1]*(7.43/2.0) + (7.43/2.0)*(7.43/2.0)*betazz[1];
  //std::cout << "newandoldbow " << barBow << " " << nw << std::endl;
}



               ///////////// FileReader  /////////////////

void FileReader::setDir(const char* dira) 
{
  strncpy(dir,dira,200);
}

void FileReader::readBarrels(Ladder** det) 
{
  char file[100];
  char line[1000];

  if(false) {  // use Ritchie files (some known errors)

    for(int b=0; b<3; b++) {
      
      // added this on 7/16/02, reordering of barrels
      int filenum;
      if(b==0) filenum = 3;
      if(b==1) filenum = 1;
      if(b==2) filenum = 2;
      
      sprintf(file,"%s/ladinfBR%d.txt",dir,filenum);
      std::cout << "reading: " << file << std::endl;
      std::ifstream ifile(file);
      //skip header lines
      for(int i=0; i<4; i++) {
	ifile.getline(line,1000);
	//      std::cout << "|skip|" << line << std::endl;
      }
      
      char* cc, *name;
      while(!ifile.eof()) {
	ifile.getline(line,1000);
	if(strlen(line)>0) {
	  name = strtok(line," \t");
	  name[4]='\0';
	  int l = atoi(name+1)/100;
	  cc = strtok(NULL," \t");  // skip barrel number
	  cc = strtok(NULL," \t");  // xend
	  bool flip = false;
	  // if word is not east, then the ladder flipped
	  // end for end
	  if(strncmp(cc+4,"East",4)) flip = true;
	  cc = strtok(NULL," \t");  // zend
	  cc = strtok(NULL," \t");  // wedge
	  int w = atoi(cc);
	  int ind = l+w*5+b*12*5;
	  
	  // the entire barrel 2 (counting 1,2,3) was assembled backward
	  // so flip all those ladders.  Variable b counts 0,1,2
	  if(b==2) { 
	    if(flip) flip = false;
	    else flip = true;
	  }
	  
	  det[ind] = new Ladder(name,b,w,l,flip);
	  
	}
      }
    }
    
  } else { // use electornic fingerprint file

    sprintf(file,"%s/ladderLocations.txt",dir);
    std::cout << "reading: " << file << std::endl;
    std::ifstream ifile(file);
    //skip header lines
    for(int i=0; i<2; i++) {
      ifile.getline(line,1000);
    }
      
    char* cc, name[5];
    char temp[4];
    int ipt;
    while(!ifile.eof()) {
      ifile.getline(line,1000);
      //cout << "line" << line << endl;
      if(strlen(line)>0) {
	cc = strtok(line," \t");
	strncpy(temp,cc+2,1);   temp[1]='\0';
	int bh = atoi(temp);
	//cout << " found bh " << temp << " " << bh << endl;

	ipt = 2;
	if(!strncmp(cc+5,"L",1)) {
	  strncpy(temp,cc+4,1);   temp[1]='\0';
	  ipt=6;
	} else {
	  strncpy(temp,cc+4,2);   temp[2]='\0';
	  ipt=7;
	}
	int w = atoi(temp);
	strncpy(temp,cc+ipt,1);   temp[1]='\0';
	int l = atoi(temp);

	cc = strtok(NULL," \t");
	strncpy(name,cc,4);       name[4]='\0';
	strncpy(temp,cc+1,3);     temp[3]='\0';
	int id = atoi(temp);

	bool zend=false;
	if(!strncmp(cc+4,"z",1)) zend = true;
	bool flip = false;
	if(bh%2==0) {
	  if(!zend) flip = true;
	} else {
	  if(zend) flip = true;
	}

	int b = bh/2;
	int ind = l+w*5+b*12*5;

	if(det[ind]) { // check that it is consistent
	  if(flip!=det[ind]->flip) cout << "Found flip inconsistency  " 
					<< b << " " << w << " " << l << endl;
	  if(id!=det[ind]->number) cout << "Found number inconsistency  " 
					<< b << " " << w << " " << l << endl;
	} else {
	  det[ind] = new Ladder(name,b,w,l,flip);
	}
	
	
      } // if strlen
    } // while !eof

  } // use fingerprint file
  
}

void FileReader::readFileQuadratic(Ladder** det) 
{
  char file[100];
  char line[1000];
  //string str;

  for(int b=0; b<3 ; b++) {
  for(int w=0; w<12; w++) {
  for(int l=0; l<5 ; l++) {
    int ind = l+w*5+b*12*5;
    Ladder& lad = *det[ind];

    //if(!strncmp(lad.name,"F131",4)) {
    //  sprintf(file,"%s/Layer%d/f%3sflt1_v2.txt",dir,l,lad.name+1);
    //} else {
      sprintf(file,"%s/Layer%d/f%3sflt1.txt",dir,l,lad.name+1);
      //}

    std::cout << "reading Warp file: " << file << std::endl;

    std::ifstream ifile(file);
    //if(ifile.open(file)!=0) {
    if(!ifile.good()) {
      std::cout << "Couldn't open file!!!!!!!!!  " << file << std::endl;
    } else {
    //skip header lines
      for(int i=0; i<3; i++) {
	ifile.getline(line,1000);
      }
      
      char* cc;

      while(!ifile.eof()) {
	ifile.getline(line,1000);
	if(strlen(line)>0) {
	  cc = strtok(line," \t");  // "X"
	  if(strncmp(cc,"X",1)) {   // if this isn't "X"
	    // it could be a new series of measurements so restart counter
	    // and keep skipping while header lines are found
	    lad.nWarp = 0;
	    while(!ifile.eof() && strncmp(cc,"X",1)) {
	      ifile.getline(line,1000);
	      cc = strtok(line," \t");  // "X"
	    }
	  }
	  cc = strtok(NULL," \t");  // x
	  double x = atof(cc);
	  cc = strtok(NULL," \t");  // "Y"
	  cc = strtok(NULL," \t");  // y
	  double y = atof(cc);
	  cc = strtok(NULL," \t");  // "Z"
	  cc = strtok(NULL," \t");  // y
	  double z = atof(cc);
	  Ladder& lad = *det[ind];
	  
	  // x needs to be shifting by 50mm
	  // then put it in ladder coords
	  // this is an average offset but erors are second order
	  // in alignment
	  x += 49.7 - 2*74.3;
	  // layer 4 is taken with a different coordinate system
	  if(l==4) x-=48.45;

	  if(fabs(x)> 2*74.3 && fabs(x)< 2*74.3+5 && fabs(z)> 0.15 ) {
	    std::cout << "Removing Ledge points" << std::endl;
	  } else {
	    bool ok = true;
	    if(fabs(x)> 2*74.3 ) {
	      cout << "bad x in readQuad: " << x << std::endl; lad.print(std::cout);
	      ok = false;
	    }
	    if(fabs(y)> 32 ) {
	      cout << "bad y in readQuad: " << y << std::endl; lad.print(std::cout);
	      ok = false;
	    }
	    if(fabs(z)> 1.0 ) {
	      cout << "bad z in readQuad: " << z << std::endl; lad.print(std::cout);
	      ok = false;
	    }


	    //***********************
	    // remove the measuremnt outside the hybrid
	    // this seems to always be 90um off the others...
	    // - OK - an understood, real effect
	    //************************
	    //
	    //if(fabs(x)> 100.3 ) {
	    //  std::cout << "removing meas outside hybrid in readQuad: " 
	    //   << x << std::endl; lad.print(std::cout);
	    //  ok = false;
	    //}

	    
	    if(ok) {
	      lad.xWarp[lad.nWarp] =  z/10.0;  // convert to cm
	      lad.yWarp[lad.nWarp] = -y/10.0;  // change to standard coords
	      lad.zWarp[lad.nWarp] =  x/10.0;
	      
	      if(lad.flip) {
		lad.yWarp[lad.nWarp] = -lad.yWarp[lad.nWarp];
		lad.zWarp[lad.nWarp] = -lad.zWarp[lad.nWarp];
	      }
	      lad.nWarp++;
	    } // end OK

	  } // end ledge check
	}  // end if string has length
      } // end while

    } // end of else if on file found

  }  // b,w,l loops
  }
  }

}

void FileReader::readFileXY(Ladder** det) 
{
  char file[100];
  char line[1000];

  for(int b=0; b<3 ; b++) {
  for(int w=0; w<12; w++) {
  for(int l=0; l<5 ; l++) {
    int ind = l+w*5+b*12*5;
    Ladder& lad = *det[ind];

    sprintf(file,"%s/Layer%d/f%3s.txt",dir,l,lad.name+1);
    std::cout << "reading XY file: " << file << std::endl;

    std::ifstream ifile(file);

    char* cc;
    char temp[100];
    while(!ifile.eof()) {
      ifile.getline(line,1000);
      if(strlen(line)>0) {
	cc = strtok(line," \t");  // cc="NS1F1" or similar label

	// non-layer 4 files have the form NS[wafer]F[position 1 or 2]
	// some non-layer 4 files have NSnML and NSnMR, which are ignored
	// layer 4 files have NS[wafer]M[LorR] 
	// but need to reject NS[wafer]T[LorR]
	bool ok = false;
	int waf,pos;

	if(!strncmp(cc,"NS",2)) {
	  if(l==4) {
	    if(!strncmp(cc+3,"M",1)) {
	      ok = true;
	      strncpy(temp,cc+2,1);
	      waf = atoi(temp)-1;
	      pos=0;
	      if(!strncmp(cc+4,"R",1)) pos=1;
	    }
	  } else {
	    if(!strncmp(cc+3,"F",1)) {
	      ok = true;
	      strncpy(temp,cc+2,1);
	      waf = atoi(temp)-1;
	      strncpy(temp,cc+4,1);
	      pos = atoi(temp)-1;
	    }
	  }
	} 

	// if this is a transverse measurement we want
	if(ok) {

	  ifile.getline(line,1000);  // read X
	  cc = strtok(line," \t");  // cc="X"
	  cc = strtok(NULL," \t");  // x
	  double x = atof(cc);
	  ifile.getline(line,1000);  // read Y
	  cc = strtok(line," \t");  // cc="Y"
	  cc = strtok(NULL," \t");  // y
	  double y = atof(cc);

	  int ipt = 2*waf+pos;
	  if(lad.yy[ipt]>-1000.0) { 
	    std::cout << "found duplicate measurement in read XY " << ipt << std::endl;
	    lad.print(std::cout);
	  }
	  if(fabs(y)>32) { 
	    std::cout << "found bad y in read XY " << y << std::endl;
	    lad.print(std::cout);
	  }
	  if(fabs(x)>2*74.3) {
	    std::cout << "found bad x in read XY " << x << std::endl;
	    lad.print(std::cout);
	  }

	  // then put it in ladder coords
	  // when the CMM x axis is along the CDF z axis (flip=false)
	  // then the CMM y measurement is in the opposite of the 
	  // CDF phi direction and the CDF local y direction

	  lad.yy[ipt] = -y/10.0;  // convert to cm
	  lad.zz[ipt] =  x/10.0;  // convert to cm

	} // end of if for NS label


	// gaps, raw measurements 
	ok=false;
	if(!strncmp(cc,"S",1)) {
	  if(l==4) {
	    if(!strncmp(cc+2,"M",1)) {
	      ok = true;
	      strncpy(temp,cc+1,1);
	      waf = atoi(temp)-1;
	      pos=0;
	      if(!strncmp(cc+3,"R",1)) pos=1;
	    }
	  } else {
	    if(!strncmp(cc+3,"T",1)) {
	      ok = true;
	      strncpy(temp,cc+1,1);
	      waf = atoi(temp)-1;
	      strncpy(temp,cc+4,1);
	      if(!strncmp(cc+2,"R",1)) pos=1;
	    }
	  }
	} 

	// if this is a gap measurement we want
	if(ok) {

	  ifile.getline(line,1000);  // read X
	  cc = strtok(line," \t");  // cc="X"
	  cc = strtok(NULL," \t");  // x
	  double x = atof(cc);
	  ifile.getline(line,1000);  // read Y
	  cc = strtok(line," \t");  // cc="Y"
	  cc = strtok(NULL," \t");  // y
	  double y = atof(cc);

	  int ipt = 2*waf+pos;
	  if(lad.yygap[ipt]>-1000.0) { 
	    std::cout << "found duplicate measurement in read XY " << ipt << std::endl;
	    lad.print(std::cout);
	  }
	  if(fabs(y)>32) { 
	    std::cout << "found bad y in read XY " << y << std::endl;
	    lad.print(std::cout);
	  }
	  if(fabs(x)>2*74.3) {
	    std::cout << "found bad x in read XY " << x << std::endl;
	    lad.print(std::cout);
	  }

	  // then put it in ladder coords
	  // when the CMM x axis is along the CDF z axis (flip=false)
	  // then the CMM y measurement is in the opposite of the 
	  // CDF phi direction and the CDF local y direction

	  lad.yygap[ipt] = -y/10.0;  // convert to cm
	  lad.zzgap[ipt] =  x/10.0;  // convert to cm

	} // end of if for SxM label for z position of gaps



	// gaps,  derived measurements
	bool s12 = !strncmp(cc,"S1-S2D",6);
	bool s23 = !strncmp(cc,"S2-S3D",6);
	bool s34 = !strncmp(cc,"S3-S4D",6);
	if(s12 || s23 || s34) {
	  ifile.getline(line,1000);  // read X
	  cc = strtok(line," \t");  // cc="XD"
	  cc = strtok(NULL," \t");  // xmeas
	  cc = strtok(NULL," \t");  // xexp
	  double xexp = atof(cc);
	  cc = strtok(NULL," \t");  // xdiff
	  double xdiff = atof(cc);

	  int ipt;
	  if(s12) ipt=0;
	  if(s23) ipt=1;
	  if(s34) ipt=2;
	  if(lad.sd[ipt]>-1000.0) {
	    std::cout << "found duplicate delta z in read XY " << ipt << std::endl;
	    lad.print(std::cout);
	  }
	  lad.sd[ipt] = xdiff/10.0;  // convert to cm
	  lad.gap[ipt] = xexp/10.0;  // convert to cm

	} // end of if for S label


      }  // end if string has length
    } // end while reading file

    // now flip the ladder coordinates around if necessary
    // this will result in all cached measurements being 
    // increasing along CDF z

    if(lad.flip) {
      double tmp;
      int j;
      for(int i=0;i<4;i++) {
	j = 7-i;
	tmp = lad.yy[i];
	lad.yy[i] = -lad.yy[j];
	lad.yy[j] = -tmp;
	tmp = lad.zz[i];
	lad.zz[i] = -lad.zz[j];
	lad.zz[j] = -tmp;
	tmp = lad.yygap[i];
	lad.yygap[i] = -lad.yygap[j];
	lad.yygap[j] = -tmp;
	tmp = lad.zzgap[i];
	lad.zzgap[i] = -lad.zzgap[j];
	lad.zzgap[j] = -tmp;
      }
      tmp = lad.sd[0];
      lad.sd[0] = lad.sd[2];
      lad.sd[2] = tmp;
      tmp = lad.gap[0];
      lad.gap[0] = lad.gap[2];
      lad.gap[2] = tmp;

    }
    
    // check that all measurements were found
    bool problem = false;
    for(int i=0; i<8; i++) {
      if(lad.zz[i]<-1000.0) problem = true;
      if(lad.yy[i]<-1000.0) problem = true;
    }
    if(problem) {
      std::cout << "found a missing XY measurement ";
      lad.print(std::cout);
    }
    problem = false;
    for(int i=0; i<3; i++) {
      if(lad.sd[i]<-1000.0) problem = true;
    }
    if(problem) {
      std::cout << "found a missing delta z measurement ";
      lad.print(std::cout);
    }


  }  // b,w,l loops
  }
  }

}

void FileReader::readFileBarrel(Ladder** det) 
{
  char file[100];
  char line[1000];

  //  for(int b=0; b<3 ; b++) {
  //for(int b=1; b<2 ; b++) {

  //   sprintf(file,"%s/barrel%1d.txt",dir,b+1);
    sprintf(file,"%s/dball.txt",dir);
    std::cout << "reading barrel file: " << file << std::endl;

    std::ifstream ifile(file);

    char* cc;

    ifile.getline(line,1000); // read header line
    int ipt;

    for(int iline=0; iline<180; iline++) {
      ifile.getline(line,1000);
      //std::cout << line << std::endl;
      cc = strtok(line," \t");  // wafer number
      int num = atoi(cc);
      ipt = -1;
      for(int b=0; b<3; b++) {
      for(int w=0; w<12; w++) {
	for(int l=0; l<5 ; l++) {
	  int ind = l+w*5+b*12*5;
	  if((*det[ind]).number == num) ipt = ind;
	}
      }
      }

      if(ipt==-1) {
	std::cout << "Couldn't match barrel meas to existing ladder !!!!!" << std::endl;
	std::cout << " num " << num << std::endl;
	return;
      }

      //Ladder	Yaw	          Dz	 Dy	Bow	Rledge	Pitch	Roll
      Ladder& lad = *det[ipt];
      cc = strtok(NULL," \t");  // yaw
      lad.barAx  = atof(cc);
      cc = strtok(NULL," \t");  // Dz
      lad.barDz  = atof(cc);
      cc = strtok(NULL," \t");  // Dy
      lad.barDy  =  - atof(cc); // flip direction from -phi to phi
      cc = strtok(NULL," \t");  // Bow
      lad.barBow = atof(cc);
      cc = strtok(NULL," \t");  // Rledge
      lad.barRwaf = atof(cc);
      cc = strtok(NULL," \t");  // Pitch
      lad.barAy = atof(cc);
      cc = strtok(NULL," \t");  // Roll
      lad.barAz = atof(cc);
    } // end loop over 60 lines
    //  } // end loop over 3 barrels


  // check that all measurements were found
  for(int i=0; i<3*12*5; i++) {
    if(det[i]->barRwaf<-1000.0) {
      std::cout << "Found missing ladder in barrel files !!!!!!" << std::endl;
    }
  }

}



      //////////////////////  Plotter   //////////////////
Plotter::Plotter() {

 hfile = new 
   TFile("Translate.root","RECREATE","Translate Ntuple");

 trTree = new TTree("TranslateTree","Translate Tree");
 trTree->Branch("barrel",&(plotLadder.barrel),"barrel/I");
 trTree->Branch("layer", &(plotLadder.layer),"layer/I");
 trTree->Branch("wedge", &(plotLadder.wedge),"wedge/I");
 trTree->Branch("number", &(plotLadder.number),"number/I");
 trTree->Branch("width", &(plotLadder.width),"width/D");

 trTree->Branch("yy",&(plotLadder.yy),"yy[8]/D");
 trTree->Branch("zz",&(plotLadder.zz),"zz[8]/D");
 trTree->Branch("yygap",&(plotLadder.yygap),"yygap[8]/D");
 trTree->Branch("zzgap",&(plotLadder.zzgap),"zzgap[8]/D");
 trTree->Branch("sd",&(plotLadder.sd),"sd[3]/D");
 trTree->Branch("gap",&(plotLadder.gap),"gap[3]/D");
 trTree->Branch("yoff",&(plotLadder.yoff),"yoff[4]/D");
 trTree->Branch("xoff",&(plotLadder.xoff),"xoff[4]/D");
 trTree->Branch("zoff",&(plotLadder.zoff),"zoff[4]/D");
 trTree->Branch("xrot",&(plotLadder.xrot),"xrot[4]/D");
 trTree->Branch("yrot",&(plotLadder.yrot),"yrot[4]/D");
 trTree->Branch("zrot",&(plotLadder.zrot),"zrot[4]/D");

 trTree->Branch("betayy",&(plotLadder.betayy),"betayy[4]/D");
 trTree->Branch("betayz",&(plotLadder.betayz),"betayz[4]/D");
 trTree->Branch("betazz",&(plotLadder.betazz),"betazz[4]/D");

 trTree->Branch("bow",&(plotLadder.bow),"bow/D");
 trTree->Branch("rawBow",&(plotLadder.rawBow),"rawBow/D");
 trTree->Branch("xofflad",&(plotLadder.xofflad),"xofflad/D");
 trTree->Branch("yrotlad",&(plotLadder.yrotlad),"yrotlad/D");

 trTree->Branch("nWarp", &(plotLadder.nWarp),"nWarp/I");

 trTree->Branch("chiWarp",&(plotLadder.chiWarp),"chiWarp[4]/D");
 trTree->Branch("hiWarp",&(plotLadder.hiWarp),"hiWarp[4]/D");

 trTree->Branch("barAx",  &(plotLadder.barAx),  "barAx/D");
 trTree->Branch("barDz",  &(plotLadder.barDz),  "barDz/D");
 trTree->Branch("barDy",  &(plotLadder.barDy),  "barDy/D");
 trTree->Branch("barBow", &(plotLadder.barBow), "barBow/D");
 trTree->Branch("barRwaf",&(plotLadder.barRwaf),"barRwaf/D");
 trTree->Branch("barAy",  &(plotLadder.barAy),  "barAy/D");
 trTree->Branch("barAz",  &(plotLadder.barAz),  "barAz/D");

}

void Plotter::makeTuple(Ladder& lad) 
{
  plotLadder = lad;
  std::cout << " filling " << plotLadder.layer <<"  " 
       << plotLadder.hiWarp[0] << " "
       << plotLadder.hiWarp[1] << " "
       << plotLadder.hiWarp[2] << std::endl;
  trTree->Fill();
  hfile->cd();
  if(lad.gRaw!=NULL) lad.gRaw->Write();
  if(lad.gFit!=NULL) lad.gFit->Write();
}

void Plotter::writeTuple() 
{
  hfile->Write();
  hfile->Close(); 
}

void Plotter::fillAlignment(Ladder** det) 
{
  

AlHolder wafOnly;
AlHolder ladOnly;
AlHolder wafAndLad;
AlHolder jigWafer;

AlDigiCode al;
  for(int b=0; b<3 ; b++) {
  for(int w=0; w<12; w++) {
  for(int l=0; l<5 ; l++) {


    int ind = l+w*5+b*12*5;
    Ladder& lad = *det[ind];

    al.setPar(b,l+1,w);
    ladOnly.setLadder(al,0.0,lad.barDy,lad.barDz,lad.barAx,
                      lad.barAy,lad.barAz);

    for(int a=0; a<4; a++) {
      int hl = a/2;
      int aa = a%2;

      // if odd half ladder, then readout is in -phi direction
      if(hl==0) {
	al.setPar(b,l+1,w,hl,aa);
        wafOnly.setWafer(al,lad.xoff[a],lad.yoff[a],lad.zoff[a],lad.xrot[a],
                         lad.yrot[a],lad.zrot[a],0.0,0.0,0.0,0.0);
        wafOnly.setWarp(al,lad.fitWarp[a][3], lad.fitWarp[a][4],
                           lad.fitWarp[a][5]);
        jigWafer.setWafer(al,lad.xoffJig[a],lad.yoff[a],lad.zoff[a],
			  lad.xrot[a],lad.yrotJig[a],lad.zrot[a],
			  0.0,0.0,0.0,0.0);
        jigWafer.setWarp(al,lad.fitWarp[a][3], lad.fitWarp[a][4],
                           lad.fitWarp[a][5]);
      } else {
	al.setPar(b,l+1,w,hl,1-aa);
        wafOnly.setWafer(al,lad.xoff[a],-lad.yoff[a],-lad.zoff[a],lad.xrot[a],
                         -lad.yrot[a],-lad.zrot[a],0.0,0.0,0.0,0.0);
        wafOnly.setWarp(al,lad.fitWarp[a][3], lad.fitWarp[a][4],
                           lad.fitWarp[a][5]);
        jigWafer.setWafer(al,lad.xoffJig[a],-lad.yoff[a],-lad.zoff[a],
			  lad.xrot[a],-lad.yrotJig[a],-lad.zrot[a],
			  0.0,0.0,0.0,0.0);
        jigWafer.setWarp(al,lad.fitWarp[a][3], lad.fitWarp[a][4],
                           lad.fitWarp[a][5]);
      }
    }

  }
  }
  }
  wafAndLad = wafOnly + ladOnly;
  wafOnly.writeToFile("wafOnly.txt");
  ladOnly.writeToFile("ladOnly.txt");
  wafAndLad.writeToFile("wafAndLad.txt");
  jigWafer.writeToFile("jigWafer.txt");

  /*
Finally the above described parameters are insertied in the 
alignment database with the following final modification.  
Odd-numbered half-ladders (starting counting with 0) are read out in the 
negative phi direction, therefore the coordinate system for those wafers only,
are rotated again by 180 degrees around the wafer's local x axis.
This is done by negating the y and z offsets, and the y and z rotations
for these wafers.  Also the wafer number within a half-ladder needs to 
be reversed so for these half-ladders the wafer number counts up
while the wafer goes in the negative CDF z direction.

*/

}