// Author: Andrew M. Rogers and Mark Wallace, NSCL 07/01/2006
//* Copyright (C) 2006-2008 Andrew M. Rogers and Mark Wallace

#include "TMcpDet.h"
#include "TMcp.h"

#include <fstream>
#include <math.h>
#include <stdio.h>

#include <TChain.h>
#include <TCutG.h>
#include <TF2.h>
#include <TFile.h>
#include <TObject.h>
#include <TPad.h>
#include <TMath.h>
#include <TRotation.h>
#include <TTree.h>


ClassImp(TMcpDet);

using namespace TMath;

////////////////////////////////////////////////////////////////////////////////
/* BEGIN_HTML
<p> The TMcp class inherits from TObject.  This Class provides the structure 
for MCP (MicroChannel Plate) detectors.  
</p>

<h3><a name="unpack">I. Unpacking</a></h3>

<h3><a name="vars">II. Important Data Members</a></h3>


END_HTML */
////////////////////////////////////////////////////////////////////////////////

const Double_t MASKDIFF = 5.000;    // mm

//extern Bool_t GetCentroid(TH2* histo,TCutG *cut,Double_t *centroid,Bool_t kPrint=kTRUE);


//______________________________________________________________________________
TMcpDet::~TMcpDet()
{
  // -- Destructor

  //  printf("%x\n",fPolyX0);
  //  for(Int_t i=0; i<=fMapOrder; i++) if(fPolyX[i]) delete fPolyX[i];
  //  for(Int_t i=0; i<=fMapOrder; i++) if(fPolyY[i]) delete fPolyY[i];

}

//______________________________________________________________________________
TMcpDet::TMcpDet(const TMcpDet &det) : TObject(det)
{
  // -- Copy constructor.

  ((TMcpDet&)det).Copy(*this);
}


//______________________________________________________________________________
void TMcpDet::Copy(TObject &det)const
{
  // -- Copy this method.
  //

  ((TMcpDet&)det).fRandom.SetSeed(0);
  ((TMcpDet&)det).fName   = fName;
  ((TMcpDet&)det).fTitle  = fTitle;
  ((TMcpDet&)det).fId     = fId;
  ((TMcpDet&)det).fNGainStages = fNGainStages;
  ((TMcpDet&)det).fMapOrder    = fMapOrder;
  ((TMcpDet&)det).fPosSumCorr  = fPosSumCorr;
  for(Int_t i=0; i<9; i++) ((TMcpDet&)det).fChargeMap[i] = fChargeMap[i];
  for(Int_t i=0; i<9; i++) ((TMcpDet&)det).fSigMap[i]    = fSigMap[i];
  ((TMcpDet&)det).fMethod = fMethod;
  for(Int_t i=0; i<4; i++){
    ((TMcpDet&)det).fHLSlope[i]  = fHLSlope[i];
    ((TMcpDet&)det).fHLOffset[i] = fHLOffset[i];
    for(Int_t j=0; j<2; j++) ((TMcpDet&)det).fHLLimits[i][j] = fHLLimits[i][j];
  }
  ((TMcpDet&)det).xScale      = xScale;
  ((TMcpDet&)det).yScale      = yScale;
  ((TMcpDet&)det).rotAlign    = rotAlign;
  ((TMcpDet&)det).xOffset     = xOffset;
  ((TMcpDet&)det).yOffset     = yOffset;
  ((TMcpDet&)det).zOffset     = zOffset;
  ((TMcpDet&)det).xShift      = xShift;
  ((TMcpDet&)det).yShift      = yShift;
  ((TMcpDet&)det).zShift      = zShift;
  ((TMcpDet&)det).foilAngle   = foilAngle;
  ((TMcpDet&)det).rotation    = rotation;
  ((TMcpDet&)det).fEulerPhi   = fEulerPhi;
  ((TMcpDet&)det).fEulerTheta = fEulerTheta;
  ((TMcpDet&)det).fEulerPsi   = fEulerPsi;
  ((TMcpDet&)det).sumThresh   = sumThresh;
  for(Int_t i=0; i<8; i++){
      ((TMcpDet&)det).cornerGain[i] = cornerGain[i];
      ((TMcpDet&)det).thresh[i]     = thresh[i];
      ((TMcpDet&)det).ped[i]        = ped[i];
  }

  for(Int_t i=0; i<=fMapOrder; i++){
    if(fPolyX[i]) ((TMcpDet&)det).fPolyX[i] = new TF2(*fPolyX[i]);
    else          ((TMcpDet&)det).fPolyX[i] = 0;
  }
  for(Int_t i=0; i<=fMapOrder; i++){
    if(fPolyY[i]) ((TMcpDet&)det).fPolyY[i] = new TF2(*fPolyY[i]);
    else          ((TMcpDet&)det).fPolyY[i] = 0;
  }
  
  for(Int_t i=0; i<2; i++){
    if(fPosSumX[i]) ((TMcpDet&)det).fPosSumX[i] = new TF2(*fPosSumX[i]);
    else            ((TMcpDet&)det).fPosSumX[i] = 0;
    if(fPosSumY[i]) ((TMcpDet&)det).fPosSumY[i] = new TF2(*fPosSumY[i]);
    else            ((TMcpDet&)det).fPosSumY[i] = 0;
  }

  //   if(fPolyX0) fPolyX0->Copy(*((TMcpDet&)det).fPolyX0);
  //   if(fPolyY0) fPolyY0->Copy(*((TMcpDet&)det).fPolyY0);
}


//______________________________________________________________________________
void TMcpDet::Clear(Option_t *) 
{
  // -- Clear MCP detector values.
  //

  posV.SetXYZ(Sqrt(-1.0),Sqrt(-1.0),Sqrt(-1.0));

  x      = sqrt(-1.0);
  y      = sqrt(-1.0);
  xHG    = sqrt(-1.0);
  yHG    = sqrt(-1.0);
  r      = sqrt(-1.0);
  z      = sqrt(-1.0);
  xRaw   = sqrt(-1.0);
  yRaw   = sqrt(-1.0);
  xRawHG = sqrt(-1.0);
  yRawHG = sqrt(-1.0);
  xRawM  = sqrt(-1.0);
  yRawM  = sqrt(-1.0);
  sum    = 0;
  sumHG  = 0;
  cMult  = 0;
  tSig   = sqrt(-1.0);

  for (Int_t z=0; z<12; z++) corner[z] = sqrt(-1.0);
}


//______________________________________________________________________________
Int_t TMcpDet::Calibrate(Long64_t entry)
{
  // -- Calibrate an entry.
  //
  // 

  // Check that a TTree has been initialized.
  if(!fChain){printf("TTree or TChain has not been initialized!!!\n"); return kFALSE;}

  Clear();
  fCon->b_charge->GetEntry(entry);
  fCon->b_time->GetEntry(entry);
  
  // Set the LOW GAIN corner values.
  if(fCon->charge[fChargeMap[0]][fSigMap[0]] < 3841 &&
     fCon->charge[fChargeMap[1]][fSigMap[1]] < 3841 &&
     fCon->charge[fChargeMap[2]][fSigMap[2]] < 3841 &&
     fCon->charge[fChargeMap[3]][fSigMap[3]] < 3841){
    for(Int_t i=0; i<4; i++){
      // If randomizing around bin center . . 
      if(parent->fAddRndm) corner[i] = (Double_t)((fCon->charge[fChargeMap[i]][fSigMap[i]]+fRandom.Rndm())-ped[i])*cornerGain[i];
      else                 corner[i] = (Double_t)((fCon->charge[fChargeMap[i]][fSigMap[i]])-ped[i])*cornerGain[i];
      if(corner[i]<0)      corner[i] = sqrt(-1.0);
    }
  }
    
  // Set the HIGH GAIN corner values.
  if(fNGainStages>1){
  if(fCon->charge[fChargeMap[4]][fSigMap[4]] < 3841 &&
     fCon->charge[fChargeMap[5]][fSigMap[5]] < 3841 &&
     fCon->charge[fChargeMap[6]][fSigMap[6]] < 3841 &&
     fCon->charge[fChargeMap[7]][fSigMap[7]] < 3841){
    for(Int_t i=4; i<8; i++){
      // If randomizing around bin center. . .
      if(parent->fAddRndm) corner[i] = (Double_t)((fCon->charge[fChargeMap[i]][fSigMap[i]]+fRandom.Rndm())-ped[i])*cornerGain[i];
      else                 corner[i] = (Double_t)((fCon->charge[fChargeMap[i]][fSigMap[i]])-ped[i])*cornerGain[i];
      if(corner[i]<0)      corner[i] = sqrt(-1.0);
    }
  }
  }

  // Set the time signal.
  Int_t cm = fChargeMap[8];
  Int_t st = fSigMap[8];
  if(fCon->charge[cm][st]<=4096 && parent->fAddRndm) tSig = (fCon->charge[cm][st]+fRandom.Rndm());
  else if(fCon->charge[cm][st]<=4096)                tSig = (fCon->charge[cm][st]);
  else                                               tSig = sqrt(-1.0);

  return 1;
}


//______________________________________________________________________________
Int_t TMcpDet::Calculate(Long64_t entry)
{
  // -- Calculate quantities related to the MCP detector for an entry.

  // Check that a TTree has been initialized.
  if(!fChain){ printf("TTree or TChain has not been initialized!!!\n"); return -1;}

  // Compute the corner multiplicy and the sum.
  for(Int_t z=0; z<4; z++){
    if(corner[z] > thresh[z])     cMult++;
    sum   += corner[z];
    if(corner[z+4] > thresh[z+4]) cMult++;
    sumHG += corner[z+4];
  }
  // If there is a signal from all 4 corners compute the RAW postion.
  if(cMult == 4){
    #ifdef EXPERIMENT_02023
    xRaw = (corner[1]+corner[3]-corner[0]-corner[2])/(corner[0]+corner[1]+corner[2]+corner[3]);
    yRaw = (corner[0]+corner[1]-corner[2]-corner[3])/(corner[0]+corner[1]+corner[2]+corner[3]);
    #else
    xRaw   = (corner[0]+corner[3]-corner[1]-corner[2])/(corner[0]+corner[1]+corner[2]+corner[3]);
    yRaw   = (corner[0]+corner[1]-corner[3]-corner[2])/(corner[0]+corner[1]+corner[2]+corner[3]);
    xRawHG = (corner[4]+corner[7]-corner[5]-corner[6])/(corner[4]+corner[5]+corner[6]+corner[7]);
    yRawHG = (corner[4]+corner[5]-corner[7]-corner[6])/(corner[4]+corner[5]+corner[6]+corner[7]);
    xRawM  = (corner[8]+corner[11]-corner[9]-corner[10])/(corner[8]+corner[9]+corner[10]+corner[11]);
    yRawM  = (corner[8]+corner[9]-corner[11]-corner[10])/(corner[8]+corner[9]+corner[10]+corner[11]);
    #endif

  }

  if(fNGainStages==2){
  // Match HIGH/LOW gain corners.
  Double_t corr;
  for(Int_t i=0; i<4; i++){
    corr = fHLSlope[i]*corner[i]+fHLOffset[i];
    if     (corner[i]>0 && corner[i]<fHLLimits[i][0]) corner[8+i] = corner[4+i];
    else if(corner[i]>0 && corner[i]>fHLLimits[i][1]) corner[8+i] = corr;
    else{
      corner[8+i] = (corner[4+i]*(fHLLimits[i][1]-corner[i])+corr*(corner[i]-fHLLimits[i][0]))/
	            (fHLLimits[i][1]-fHLLimits[i][0]);
      //corner[8+i] = fHLSlope[i]*corner[i]+fHLOffset[i];
    }
  }
  }

  posV.SetXYZ(xRaw+xOffset,yRaw+yOffset,0);
  posV.RotateZ(rotAlign*DegToRad());
  xRaw = posV.X();
  yRaw = posV.Y();
  posV.Clear();

  return 1;
}

// //______________________________________________________________________________
// Bool_t TMcpDet::Calculate()
// {
//   // -- Calculate quantities related to the MCP detector.
//   // Before the detector calculations can be done the fitting functions must be 
//   // constructed from a mask run.
//   //
//   // The following calculations are performed,
//   //  1. The x and y positions are evaluated using the xRaw and yRaw data along
//   //     with the rotation angle of the target.
//   //  2. A radial distance is also computed.

//   // Check that a TTree has been initialized.
//   if(!fChain){
//     printf("TTree or TChain has not been initialized!!!\n");
//     return kFALSE;
//   }

//   b_x->Reset();
//   b_y->Reset();
//   b_z->Reset();
//   b_r->Reset();

//   cout << fChain->GetEntries() << endl;
//   //  for(Long64_t i=0; i<fChain->GetEntries(); i++){

//   Int_t num = fChain->GetEntries();
//   TH2D *h_xRawTmp = new TH2D("h_xRawTmp","h_xRawTmp",600,-30,30,600,-30,30);
//   TH2D *h_yRawTmp = new TH2D("h_yRawTmp","h_yRawTmp",600,-30,30,600,-30,30);
//   Double_t step = 0.1;

//   // THIS METHOD IS CURRENTLY NOT OPERATIONAL.
//   if(fMethod==1){
//   cin >> num;
//   // Fill an array to search.
//   Double_t xx=-30,yy=30;
//   cout << "Total gird points: " << 600*600 << endl;
//   Int_t nGPoint =0;
//   while(xx<=30 && yy>=-30){
//     cout << nGPoint << "\r";
//     //    cout << "xx: " << xx << " yy: " << yy << " \r";
//     cout.flush();
//     nGPoint++;
//     Int_t    nRowSplines=0;
//     Double_t xTmp[13];
//     Double_t yTmp[13];
//     memset(xTmp,'0',sizeof(xTmp));
//     memset(yTmp,'0',sizeof(yTmp));
//     for(Int_t j=0; j<13; j++){
//       if(fRowX[j]){
// 	xTmp[nRowSplines] = fRowX[j]->Eval(yy);
// 	nRowSplines++;
//       }
//     }
//     fBicubicX = new TSpline3("BicubicX",GridX,xTmp,nRowSplines,"0",0,0);
//     h_xRawTmp->SetBinContent(h_xRawTmp->FindBin(xx,yy),fBicubicX->Eval(xx));
//     nRowSplines = 0;
//     for(Int_t j=0; j<13; j++){
//       if(fRowY[j]){
// 	yTmp[nRowSplines] = fRowY[j]->Eval(xx);
// 	nRowSplines++;
//       }
//     }
//     fBicubicY = new TSpline3("BicubicY",GridY,yTmp,nRowSplines,"0",0,0);
//     h_yRawTmp->SetBinContent(h_yRawTmp->FindBin(xx,yy),fBicubicY->Eval(yy));
    
//     delete fBicubicX;
//     delete fBicubicY;
//     xx += step;
//     if(xx>29.5) {
//       yy = yy-step; 
//       xx = -30;
//     }
//   }
//   printf("Filled Grid\n");
//   }


//   // Calculated 3D interaction position.
//   for(Long64_t i=0; i<num; i++){
//     Clear();
//     b_xRaw->GetEntry(i);
//     b_yRaw->GetEntry(i);

//     // Use a simple mixed polynomial fit method.
//     if(fMethod==0){
//       posV->SetXYZ(xRaw+xOffset,yRaw+yOffset,0);
//       // Rotation used to align beam mask holes to alpha mask holes.
//       // SHOULD WE DO THIS???  SHOULD IT BE DONE AFTER EVAL??? I think so now.
//       //      posV->RotateZ(rotation*DegToRad());
//       x = fPolyX0->Eval(posV->X(),posV->Y());
//       y = fPolyY0->Eval(posV->X(),posV->Y());
      

//       // Rotate the vector around the z-axis, this will change the frame of reference from
//       // the mask system to a system consitent with both MCP's.
//       // Rotate the vector around the vertical (y-axis) to acount for the foil angle.
//       posV->SetXYZ(x,y,0);
//       posV->RotateZ(rotAlign*DegToRad());
//       posV->RotateZ(rotation*DegToRad());
//       posV->RotateY(foilAngle*DegToRad());
//       // Apply shifts for target postion.
//       posV->SetX(posV->X()+xShift);
//       posV->SetY(posV->Y()+yShift);
//       posV->SetZ(posV->Z()+zShift);
//       // Apply Euler angle rotations to switch into standard coordinate system.
//       TRotation rot0,rot1,rot2,trans;
//       rot0.SetToIdentity();
//       rot1.SetToIdentity();
//       rot2.SetToIdentity();
//       trans.SetToIdentity();
//       rot0.SetXEulerAngles(fEulerPhi*DegToRad(),0,0);
//       trans = rot0.Inverse();
//       posV->Transform(trans);
//       trans.SetToIdentity();
//       rot1.SetXEulerAngles(0,fEulerTheta*DegToRad(),0);
//       trans = rot1.Inverse();
//       posV->Transform(trans);
//       trans.SetToIdentity();
//       rot2.SetXEulerAngles(0,0,fEulerPsi*DegToRad());
//       trans = rot2.Inverse();
//       posV->Transform(trans);

//       x = posV->X();
//       y = posV->Y();
//       z = posV->Z();
//     }
    
//     // Use the bicubic spline method. NOT IMPLEMENTED!!!
//     if(fMethod==1){
//       cout << i << "\r";
//       cout.flush();

//       x = -29;
//       y = -29;

//       // Search the histo.
//       Double_t xb,yb;
//       Double_t lowestX=0.05, lowestY=0.05;
//       if(!isnan(xRaw) && !isnan(yRaw)){
// 	for(Int_t jj=0; jj<600; jj++){
// 	  for(Int_t kk=0; kk<600; kk++){
// 	    xb = h_xRawTmp->GetBinContent(jj,kk);
// 	    yb = h_yRawTmp->GetBinContent(jj,kk);
// 	    if(Abs(xb-xRaw)<lowestX) {x = -30+jj*step; lowestX = Abs(xb-xRaw);}
// 	    if(Abs(xb-xRaw)<lowestY) {y =  30-kk*step; lowestY = Abs(yb-yRaw);}
// 	  }
// 	}
//       }else{
// 	x = sqrt(-1.0);
// 	y = sqrt(-1.0);
//       }
//     } 

//     r = Sqrt(x*x+y*y);

//     b_x->Fill();
//     b_y->Fill();
//     b_z->Fill();
//     b_r->Fill();
//     b_posV->Fill();
//   }

//   delete h_xRawTmp;
//   delete h_yRawTmp;

//   printf("%s is now calculated.\n",GetName());

//   return kTRUE;
// }


//______________________________________________________________________________
Bool_t TMcpDet::CreateMaskFile(TH2 *histo)
{
  // -- Generate the mask calibration file from TCutG objects.
  // This method assumes that the TCutG objects are named,
  // hole#_# where #_# are the row and column number from the 
  // matrix below,
  //            0 1 2 3 4 5 6 7 8 9 . . 12
  //            
  //      0     o o o o o o o o o o o o o 
  //      1     o o o o o o o o o o o o o
  //      2     o o o o o o o o o o o o o
  //      3     o o o o o o o o o o o o o
  //      4     o o o L L L L o o o o o o
  //      5     o o o o o o L o o o o o o
  //      6     o o o o o o L o o o o o o
  //      7     o o o o o o o o o o o o o
  //      8     o o o o o o o o o o o o o
  //      9     o o o o o o o o o o o o o
  //     10     o o o o o o o o o o o o o
  //     11     o o o o o o o o o o o o o
  //     12     o o o o o o o o o o o o o

  const Double_t diff = 5.000;    // mm

  Char_t fileName[100];
  sprintf(fileName,"%sMask.dat",GetName());
  ofstream file(fileName,ios::out);
  cout << "Creating File:  " << fileName << endl;
  if(!file.is_open()){cout << "File Creation FAILED!!!" << endl;  return kFALSE;}

  file << "Row \t" << "Col \t" 
       << "X \t" << "Y \t" 
       << "XRaw \t\t" << "YRaw \t\t" 
       << "Counts \t" << " sigmaX \t\t" << "sigmaY" << endl;

  TCutG  *cut;
  Char_t name[30],oBuf[2000];
  Double_t cen[5];

  for(Int_t i=0; i<13; i++){
    for(Int_t j=0; j<13; j++){
      sprintf(name,"hole%u_%u",i,j);
      if(gPad->FindObject(name)){
	cout << "Found TCutG Object: " << name << endl;
	cut = (TCutG*)gPad->FindObject(name);
	GetCentroid(histo,cut,cen,kTRUE);
	sprintf(oBuf,"%3i \t%3i \t% 6.2lf \t% 6.2lf \t% 12.8lf \t% 12.8lf \t%6i \t% 12.8lf \t% 12.8lf\n",
		i,j,(Double_t)(-30.0+j*diff),(Double_t)(30.0-i*diff),cen[0],cen[1],(Int_t)cen[2],cen[3],cen[4]);
	file << oBuf;
// 	file << i << " \t" << j << " \t"
// 	     << (-30.000 + j*diff) << " \t" << (30.000 - i*diff) << " \t"
// 	     << cen[0] << " \t" << cen[1] << " \t" << cen[2] << " \t"
// 	     << cen[3] << " \t" << cen[4] << endl;
      }
    }
  }
  file.close();
  return kTRUE;
}


//______________________________________________________________________________
Int_t TMcpDet::GenPrime(Long64_t entry)
{
  // -- Generate primary data for an entry.
  //

  // Check that a TTree has been initialized.
  if(!fChain){ printf("TTree or TChain has not been initialized!!!\n"); return -1;}

//   TH2D *h_xRawTmp = new TH2D("h_xRawTmp","h_xRawTmp",600,-30,30,600,-30,30);
//   TH2D *h_yRawTmp = new TH2D("h_yRawTmp","h_yRawTmp",600,-30,30,600,-30,30);
//   Double_t step = 0.1;
//
//   // THIS METHOD IS CURRENTLY NOT OPERATIONAL.
//   if(fMethod==1){
//     cin >> num;
//     // Fill an array to search.
//     Double_t xx=-30,yy=30;
//     cout << "Total gird points: " << 600*600 << endl;
//     Int_t nGPoint =0;
//     while(xx<=30 && yy>=-30){
//       cout << nGPoint << "\r";
//       //    cout << "xx: " << xx << " yy: " << yy << " \r";
//       cout.flush();
//       nGPoint++;
//       Int_t    nRowSplines=0;
//       Double_t xTmp[13];
//       Double_t yTmp[13];
//       memset(xTmp,'0',sizeof(xTmp));
//       memset(yTmp,'0',sizeof(yTmp));
//       for(Int_t j=0; j<13; j++){
// 	if(fRowX[j]){
// 	  xTmp[nRowSplines] = fRowX[j]->Eval(yy);
// 	  nRowSplines++;
// 	}
//       }
//       fBicubicX = new TSpline3("BicubicX",GridX,xTmp,nRowSplines,"0",0,0);
//       h_xRawTmp->SetBinContent(h_xRawTmp->FindBin(xx,yy),fBicubicX->Eval(xx));
//       nRowSplines = 0;
//       for(Int_t j=0; j<13; j++){
// 	if(fRowY[j]){
// 	  yTmp[nRowSplines] = fRowY[j]->Eval(xx);
// 	  nRowSplines++;
//       }
//       }
//       fBicubicY = new TSpline3("BicubicY",GridY,yTmp,nRowSplines,"0",0,0);
//       h_yRawTmp->SetBinContent(h_yRawTmp->FindBin(xx,yy),fBicubicY->Eval(yy));
      
//       delete fBicubicX;
//       delete fBicubicY;
//       xx += step;
//       if(xx>29.5) {
// 	yy = yy-step; 
// 	xx = -30;
//       }
//     }
//     printf("Filled Grid\n");
//   }


  // Calculated 3D interaction position.
  Clear();
  fConD->b_xRaw->GetEntry(entry);
  fConD->b_yRaw->GetEntry(entry);
  fConD->b_sum->GetEntry(entry);

  // Use a simple mixed polynomial fit method.
  if(fMethod==0){
    posV.SetXYZ(fConD->xRaw,fConD->yRaw,0);
    //    posV.RotateZ(rotAlign*DegToRad());
    // Rotation used to align beam mask holes to alpha mask holes.
    // SHOULD WE DO THIS???  SHOULD IT BE DONE AFTER EVAL??? I think so now.
    //      posV.RotateZ(rotation*DegToRad());
    x = fPolyX[0]->Eval(posV.X(),posV.Y());
    y = fPolyY[0]->Eval(posV.X(),posV.Y());
    posV.SetXYZ(x,y,0);
    
    // Apply shifts for target postion.
    posV.SetX(posV.X()+xShift);
    posV.SetY(posV.Y()+yShift);
    posV.SetZ(posV.Z()+zShift);

    // Rotate the vector around the z-axis, this will change the frame of reference from
    // the mask system to a system consitent with both MCP's.
    posV.RotateZ(rotation*DegToRad());

    // Apply second round of mask hole mapping.
    if(fMapOrder>0){
      posV.SetX(fPolyX[1]->Eval(posV.X(),posV.Y()));
      posV.SetY(fPolyY[1]->Eval(posV.X(),posV.Y()));
    }

    // Apply sum dependance correction.
    if(fPosSumCorr){
      Double_t xOff    = fPosSumX[0]->Eval(posV.Y(),posV.X());
      Double_t yOff    = fPosSumY[0]->Eval(posV.Y(),posV.X());
      Double_t xSlope  = fPosSumX[1]->Eval(posV.X(),posV.Y());
      Double_t ySlope  = fPosSumY[1]->Eval(posV.X(),posV.Y());

      if(fConD->sum!=0){
	posV.SetX(posV.X()-fConD->sum*xSlope);
	posV.SetY(posV.Y()-fConD->sum*ySlope);
      }
    }

    // Apply third round of mas hole mapping.
    if(fMapOrder>1){
      posV.SetX(fPolyX[2]->Eval(posV.X(),posV.Y()));
      posV.SetY(fPolyY[2]->Eval(posV.X(),posV.Y()));
    }

    // Rotate the vector around the vertical (y-axis) to acount for the foil angle.
    //    posV.RotateZ(rotAlign*DegToRad());
    posV.RotateY(foilAngle*DegToRad());

    // Apply Euler angle rotations to switch into standard coordinate system.
    TRotation rot0,rot1,rot2,trans;
    rot0.SetToIdentity();
    rot1.SetToIdentity();
    rot2.SetToIdentity();
    trans.SetToIdentity();
    rot0.SetXEulerAngles(fEulerPhi*DegToRad(),0,0);
    trans = rot0.Inverse();
    posV.Transform(trans);
    trans.SetToIdentity();
    rot1.SetXEulerAngles(0,fEulerTheta*DegToRad(),0);
    trans = rot1.Inverse();
    posV.Transform(trans);
    trans.SetToIdentity();
    rot2.SetXEulerAngles(0,0,fEulerPsi*DegToRad());
    trans = rot2.Inverse();
    posV.Transform(trans);
    
    x = posV.X();
    y = posV.Y();
    z = posV.Z();

    
    // TEMPORARY RECALIBRATION !!!!!!!!!!!!!
    //    if(strcmp(fName,"mcp0")==0) y = y+2.25;
//     Double_t a0[2] = {-0.9725, 1.063};
//     Double_t a1[2] = {0.07311, 1.186};

//     if(fName=="mcp0") x = a0[0]+a0[1]*x;
//     if(fName=="mcp1") x = a1[0]+a1[1]*x;
    // ---------------------------
  }

  r = Sqrt(x*x+y*y);
  return 1;
  
//   // Use the bicubic spline method. NOT IMPLEMENTED!!!
//   if(fMethod==1){
//     cout << i << "\r";
//     cout.flush();
    
//     x = -29;
//     y = -29;
    
//     // Search the histo.
//     Double_t xb,yb;
//     Double_t lowestX=0.05, lowestY=0.05;
//     if(!isnan(xRaw) && !isnan(yRaw)){
//       for(Int_t jj=0; jj<600; jj++){
// 	for(Int_t kk=0; kk<600; kk++){
// 	  xb = h_xRawTmp->GetBinContent(jj,kk);
// 	  yb = h_yRawTmp->GetBinContent(jj,kk);
// 	  if(Abs(xb-xRaw)<lowestX) {x = -30+jj*step; lowestX = Abs(xb-xRaw);}
// 	  if(Abs(xb-xRaw)<lowestY) {y =  30-kk*step; lowestY = Abs(yb-yRaw);}
// 	}
//       }
//     }else{
//       x = sqrt(-1.0);
//       y = sqrt(-1.0);
//     }
//   } 
//   delete h_xRawTmp;
//   delete h_yRawTmp;
}



//______________________________________________________________________________
void TMcpDet::InitClass(TString name, Int_t ID,Int_t nStages)
{
  // -- Initialize an MCP detector.
  //
  
  Char_t title[500];
  SetId(ID);
  sprintf(title,"MCP%i Tracking Detector",fId);
  SetNameTitle(name,title);

  fNGainStages = nStages;
  fMapOrder    = 0;
  fPosSumCorr  = kFALSE;

  fMethod     = 0;

  xScale      = 1;
  yScale      = 1;
  rotAlign    = 0;
  xOffset     = 0;
  yOffset     = 0;
  zOffset     = 0;
  xShift      = 0;
  yShift      = 0;
  zShift      = 0;
  foilAngle   = 0.00;
  fEulerPhi   = 0.00;
  fEulerTheta = 0.00;
  fEulerPsi   = 0.00;
  rotation    = 0.00;
  sumThresh   = 0;

  memset(fHLSlope,'1',sizeof(fHLSlope));
  memset(fHLOffset,'\0',sizeof(fHLOffset));
  for(Int_t i=0; i<4; i++){
    for(Int_t j=0; j<2; j++) fHLLimits[i][j] = 0;
  }

  for (Int_t i=0; i<8; i++){
    thresh[i]     = 0;
    cornerGain[i] = 1;
    ped[i]        = 0;
  }

  posV.SetXYZ(Sqrt(-1.0),Sqrt(-1.0),Sqrt(-1.0));

  Clear();
}


//______________________________________________________________________________
void TMcpDet::InitTree(TTree *tree)
{
  // -- Initialize the tree.
  //

  fChain   = tree;
  fCurrent = -1;

  b_xRaw   = fChain->GetBranch(fName+".xRaw");
  b_yRaw   = fChain->GetBranch(fName+".yRaw");
  b_xRawHG = fChain->GetBranch(fName+".xRawHG");
  b_yRawHG = fChain->GetBranch(fName+".yRawHG");
  b_xRawM  = fChain->GetBranch(fName+".xRawM");
  b_yRawM  = fChain->GetBranch(fName+".yRawM");
  b_x      = fChain->GetBranch(fName+".x");
  b_y      = fChain->GetBranch(fName+".y");
  b_xHG    = fChain->GetBranch(fName+".xHG");
  b_yHG    = fChain->GetBranch(fName+".yHG");
  b_z      = fChain->GetBranch(fName+".z");
  b_r      = fChain->GetBranch(fName+".r");
  b_sum    = fChain->GetBranch(fName+".sum");
  b_sumHG  = fChain->GetBranch(fName+".sumHG");
  b_cMult  = fChain->GetBranch(fName+".cMult");
  b_tSig   = fChain->GetBranch(fName+".tSig");
  b_corner = fChain->GetBranch(fName+".corner[12]");
  b_posVX  = fChain->GetBranch(fName+".posV.fX");
  b_posVY  = fChain->GetBranch(fName+".posV.fY");
  b_posVZ  = fChain->GetBranch(fName+".posV.fZ");
}


//______________________________________________________________________________
Bool_t TMcpDet::FitMask(Int_t mapOrder)
{
  // -- Generates the fitting functions from the mask holes.
  // 

  Int_t m = mapOrder;
  Char_t name[100];

  // Use a simple mixed polynomial fit.
  if(fMethod==0){
    Double_t rangelowX = g_maskHolesX[m].GetXmin();
    Double_t rangeupX  = g_maskHolesX[m].GetXmax();
    Double_t rangelowY = g_maskHolesY[m].GetYmin();
    Double_t rangeupY  = g_maskHolesY[m].GetYmax();

    if(fPolyX[m]) fPolyX[m]->Delete();
    if(fPolyY[m]) fPolyY[m]->Delete();

    sprintf(name,"f_%sMaskX%iFnc",GetName(),m);
    fPolyX[m] = new TF2(name, function, rangelowX-1, rangeupX+1, rangelowX-1, rangeupX+1, 10);
    sprintf(name,"f_%sMaskY%iFnc",GetName(),m);
    fPolyY[m] = new TF2(name, function, rangelowY-1, rangeupY+1, rangelowY-1, rangeupY+1, 10);

    fPolyX[m]->SetParNames("a0","a1","a2","a3","a4","a5","a6","a7","a8","a9");
    fPolyY[m]->SetParNames("b0","b1","b2","b3","b4","b5","b6","b7","b8","b9"); 

    fPolyX[m]->SetParameters(0,1,0,0,0,0,0,0,0,0);
    fPolyY[m]->SetParameters(0,0,1,0,0,0,0,0,0,0);
 
    //    if(mapOrder==0){
      g_maskHolesX[m].Fit(fPolyX[m],"QEN"); 
      g_maskHolesY[m].Fit(fPolyY[m],"QEN"); 
//     }else{
//       g_maskHolesX[m].Fit(fPolyX[m],"QENW"); 
//       g_maskHolesY[m].Fit(fPolyY[m],"QENW"); 
//     }

    printf("\t%s Coefficients \n",GetName());
    printf("X Polynomial: \t");
    for(Int_t i=0; i<10; i++) printf(" % 12.8lf",fPolyX[m]->GetParameter(i));
    printf("\nY Polynomial: \t");
    for(Int_t i=0; i<10; i++) printf(" % 12.8lf",fPolyY[m]->GetParameter(i));
    printf("\n");
    

//     Char_t *overW;
//     cout << "Overwrite fitting function in file? (y/n) ";
//     cin  >> overW;
//     TFile *curDir = (TFile*)gDirectory->GetFile();
//     TFile fFnc("mcpMaskFitFncs.root","UPDATE","MCP Mask Mapping Functions",0);
//     fFnc.cd();
//     if(strcmp(overW,"y")==0 || strcmp(overW,"Y")==0){
//       fPolyX[m]->Write("",TObject::kOverwrite);
//       fPolyY[m]->Write("",TObject::kOverwrite);  
//     }else{
//       fPolyX[m]->Write();
//       fPolyY[m]->Write();  
//     }

//    curDir->cd();
  }

  // Use the bicubic spline method.
  // The row splines which are used were created in the LoadMaskFile() method.
  if(fMethod==1){
    cout << "Functions are generated during Calculate Method." << endl;
  }
  
  return kTRUE;
}


//______________________________________________________________________________
Bool_t TMcpDet::LoadMaskFile(Char_t *fileName,Int_t mapOrder,Int_t fUseErrors)
{
  // -- Load the mask file for calibration.
  //

  Int_t  m = mapOrder;
  Char_t name[1000];

  // Clear the graphs.
  g_maskHolesX[m].Clear();
  g_maskHolesY[m].Clear();

  sprintf(name,"g_maskHolesX%i",m);
  g_maskHolesX[m].SetName(name);
  sprintf(name,"g_maskHolesY%i",m);
  g_maskHolesY[m].SetName(name);

  // Try to open the file.
  ifstream file(fileName,ios::in);
  if(!file.is_open()){ cout << "Could not open file." << endl; return kFALSE;}
  Char_t line[400];
  file.getline(line,400);

  // Make the grid for the bicubic splines.
  for(Int_t i=0; i<13; i++){ 
    GridX[i] = -30.00 + i*MASKDIFF;
    GridY[i] = -30.00 + i*MASKDIFF;
    for(Int_t j=0; j<13; j++){
      GridFX[i][j] = -1;   
      GridFY[i][j] = -1;
    }
  }

  Int_t row,            // Row of mask holes.
    col,                // Column of mask holes.
    counts,             // Counts in centroid of mask hole.
    nPoints=0;          // Number of holes used.
  Double_t x1Raw,       // X raw position of centroid.
    x2Raw,              // Y raw position of centroid.
    x1Mask,             // X position of mask hole.
    x2Mask,             // Y position of mask hole.
    x1Error,            // X position sigma error.
    x2Error;            // Y position sigma error.

  Double_t xA[1000],yA[1000],xrawA[1000],yrawA[1000];
  Double_t xF[1000],yF[1000],xrawF[1000],yrawF[1000];
  // Only Scale the RAW Calibration Function.
  Double_t XSCALE = 1;
  Double_t YSCALE = 1;
  if(mapOrder==0) {XSCALE=xScale; YSCALE=yScale;}

  while(!file.eof()){
    file >> row >> col >>  x1Mask >> x2Mask >> x1Raw >> x2Raw >> counts >> x1Error >> x2Error;
    //    printf("%i %i %lf %lf %lf %lf %i %lf %lf\n",row, col,  x1Mask, x2Mask, x1Raw, x2Raw, counts, x1Error, x2Error);
    
    if(fUseErrors==3){
      
    }else if(fUseErrors==2){
      g_maskHolesX[m].SetPoint(nPoints, x1Raw*XSCALE, x2Raw*YSCALE, x1Mask);
      g_maskHolesX[m].SetPointError(nPoints, x1Error/2, x2Error/2,0.0000001);
      g_maskHolesY[m].SetPoint(nPoints, x1Raw*XSCALE, x2Raw*YSCALE, x2Mask);
      g_maskHolesY[m].SetPointError(nPoints, x1Error/2, x2Error/2,0.0000001);
      // Create mask centroid graph.
      g_maskCentroids[m].SetPoint(nPoints,x1Raw*XSCALE,x2Raw*YSCALE);
      g_maskCentroids[m].SetPointError(nPoints,x1Error/2,x2Error/2);
    }else if(fUseErrors==1){
      g_maskHolesX[m].SetPoint(nPoints, x1Raw, x2Raw, x1Mask);
      g_maskHolesX[m].SetPointError(nPoints, counts, counts,0.0000001);
      g_maskHolesY[m].SetPoint(nPoints, x1Raw, x2Raw, x2Mask);
      g_maskHolesY[m].SetPointError(nPoints, counts, counts,0.0000001);
      // Create mask centroid graph.
      g_maskCentroids[m].SetPoint(nPoints,x1Raw,x2Raw);
      g_maskCentroids[m].SetPointError(nPoints,counts,counts);
    }else{
      g_maskHolesX[m].SetPoint(nPoints, x1Raw, x2Raw, x1Mask);
      g_maskHolesY[m].SetPoint(nPoints, x1Raw, x2Raw, x2Mask);
      // Create mask centroid graph.
      g_maskCentroids[m].SetPoint(nPoints,x1Raw,x2Raw);
      g_maskCentroids[m].SetPointError(nPoints,x1Error/2,x2Error/2);
    }
    
    GridFX[row][col] = x1Raw;
    GridFY[row][col] = x2Raw;

    nPoints++;
  }
  // Sort the arrays for input into TSpline3 constructor. (They must be acending)
//   Double_t x2RawSort[1000],xSorty[1000],x1RawSort[1000],ySortx[1000];
//   Double_t FX[1000][1000],FY[1000][1000];
//   Int_t positionY,equalY;
//   Int_t positionX,equalX;
//   for(Int_t i=0; i<nPoints; i++){
//     positionX = 0;
//     equalX    = -1;
//     positionY = 0;
//     equalY    = -1;
//     for(Int_t j=0; j<nPoints; j++){
//       if(xrawA[j]<xrawA[i]) positionX++;
//       if(xrawA[j]==xrawA[i]) equalX++;
//       if(yrawA[j]<yrawA[i]) positionY++;
//       if(yrawA[j]==yrawA[i]) equalY++;
//     }
//     positionX = positionX + equalX;
//     positionY = positionY + equalY;
//     x1RawSort[positionX] = xrawA[i];
//     ySortx[positionX]   = yA[i];
//     x2RawSort[positionY] = yrawA[i];
//     xSorty[positionY]   = xA[i];
//   }
//   // Fill the X and Y tables.
//   for(Int_t i=0; i<nPoints; i++){
//     for(Int_t j=0; j<nPoints; j++){
//     }
//   }
//   for(Int_t i=0; i<nPoints; i++){
//     cout << yrawA[i] << " \t" << x2RawSort[i] << " \t" << xA[i] << " \t" << xSorty[i] << endl;
//   }



    
  // Create the row 1D splines.
//   Char_t title[500];
//   TSpline3 *sp[13];
//   Int_t nRowPoints=0;
//   Double_t xRow[13],
//     yRow[13];
//   for(Int_t i=0; i<13; i++){
//     sprintf(title,"MCP X Row %i Cubic Spline",i);
//     memset(xRow,'0',sizeof(xRow));
//     memset(yRow,'0',sizeof(yRow));
//     nRowPoints = 0;
//     for(Int_t j=0; j<13; j++){
//       if(GridFX[i][j] != -1){
// 	xRow[nRowPoints] = GridY[j];
// 	yRow[nRowPoints] = GridFX[i][j];
// 	//	printf("%f \t%f \t%i\n",xRow[nRowPoints],yRow[nRowPoints],nRowPoints);
// 	nRowPoints++;

// 	//       	sp[i] = new TSpline3(title,xRow,yRow,nRowPoints,"0",0,0);
//        	fRowX[i] = new TSpline3(title,xRow,yRow,nRowPoints,"0",0,0);
// 	//	fRowX[i] = sp[i];
//       }
//     }
//     sprintf(title,"MCP Y Row %i Cubic Spline",i);
//     memset(xRow,'0',sizeof(xRow));
//     memset(yRow,'0',sizeof(yRow));
//     nRowPoints = 0;
//     for(Int_t j=0; j<13; j++){
//       if(GridFY[i][j] != -1){
// 	xRow[nRowPoints] = GridX[j];
// 	yRow[nRowPoints] = GridFY[i][j];
// 	//	printf("%f \t%f \t%i\n",xRow[nRowPoints],yRow[nRowPoints],nRowPoints);
// 	nRowPoints++;

// 	//       	sp[i] = new TSpline3(title,xRow,yRow,nRowPoints,"0",0,0);
//        	fRowY[i] = new TSpline3(title,xRow,yRow,nRowPoints,"0",0,0);
// 	//	fRowY[i] = sp[i];
//       }
//     }
//   }

  file.close();

  FitMask(m);
  
  return kTRUE;
}


//______________________________________________________________________________
Bool_t TMcpDet::LoadCorrections(Char_t *fileName)
{
  // --
  //

  TFile cFile(fileName,"READ");
  Char_t name[500];
  sprintf(name,"mcp%iPosSumXOffset",GetId());
  fPosSumX[0] = new TF2(*(TF2*)cFile.FindObjectAny(name));
  sprintf(name,"mcp%iPosSumYOffset",GetId());
  fPosSumY[0] = new TF2(*(TF2*)cFile.FindObjectAny(name));
  sprintf(name,"mcp%iPosSumXSlopes",GetId());
  fPosSumX[1] = new TF2(*(TF2*)cFile.FindObjectAny(name));
  sprintf(name,"mcp%iPosSumYSlopes",GetId());
  fPosSumY[1] = new TF2(*(TF2*)cFile.FindObjectAny(name));
  fPosSumCorr = kTRUE;

  printf("Loaded MCP Sum Correction Functions.\n");
  return kTRUE;
}



//______________________________________________________________________________
void TMcpDet::SetId(const Int_t ID)
{
  // -- Set the MCP ID number. In general 0 = Upstream, 1 = Target.
  //
  
  fId = ID;
}


//______________________________________________________________________________
void TMcpDet::SetName(const Char_t *name)
{
  // -- Change (i.e. set) the name of the TMcpDet.
  //

  fName = name;
}


//______________________________________________________________________________
void TMcpDet::SetNameTitle(const Char_t *name, const Char_t *title)
{
  // -- Change (i.e. set) all the TMcpDet parameters (name and title).
  //

  fName  = name;
  fTitle = title;
}


//______________________________________________________________________________
void TMcpDet::SetTitle(const Char_t *title)
{
  // -- Change (i.e. set) the title of the TMcpDet.
  //

  fTitle = title;
}


/******************* The fitting function *********************************/

Double_t function(Double_t *x, Double_t *a)
{
    Double_t val;

    val = a[0] 
      + a[1]*x[0] 
      + a[2]*x[1] 
      + a[3]*x[0]*x[0] 
      + a[4]*x[1]*x[1] 
      + a[5]*x[0]*x[0]*x[0] 
      + a[6]*x[1]*x[1]*x[1] 
      + a[7]*x[0]*x[1] 
      + a[8]*x[0]*x[0]*x[1] 
      + a[9]*x[0]*x[1]*x[1];

    return val;
}


Double_t fitCorners(Double_t *x, Double_t *par)
{
  Double_t val;

  val = par[3]*TMath::LogNormal(x[0],par[0],par[1],par[2]);

  return val;
}


Bool_t GetCentroid(TH2* histo,TCutG *cut,Double_t *centroid,Bool_t kPrint)
{
  // -- Calculate the x-y centroid of a TCutG object.
  // 

  Bool_t debugGC = kFALSE;

  Double_t val[2]   = {0.000,0.000};
  Double_t error[2] = {0.000,0.000};
  Double_t sumX = 0,
    sumY = 0;

  if (!histo) {printf("Cannot find histo!!!\n"); return kFALSE;}

  Int_t n = cut->GetN();
  Double_t *X = cut->GetX();
  Double_t *Y = cut->GetY();
  Double_t xmin = 1e200;
  Double_t xmax = -xmin;
  Double_t ymin = xmin;
  Double_t ymax = xmax;
  for (Int_t i=0;i<n;i++) {
    if (X[i] < xmin) xmin = X[i];
    if (X[i] > xmax) xmax = X[i];
    if (Y[i] < ymin) ymin = Y[i];
    if (Y[i] > ymax) ymax = Y[i];
  }
  TAxis *xaxis = histo->GetXaxis();
  TAxis *yaxis = histo->GetYaxis();
  Int_t binx1  = xaxis->FindBin(xmin);
  Int_t binx2  = xaxis->FindBin(xmax);
  Int_t biny1  = yaxis->FindBin(ymin);
  Int_t biny2  = yaxis->FindBin(ymax);
  Int_t nbinsx = histo->GetNbinsX();

  // Loop on bins for which the bin center is in the cut 
  Int_t bin, binx, biny, count=0;

  for (biny=biny1;biny<=biny2;biny++) {
    Double_t y = yaxis->GetBinCenter(biny);
    for (binx=binx1;binx<=binx2;binx++) {
      Double_t x = xaxis->GetBinCenter(binx);
      if (!cut->IsInside(x,y)) continue;
      //      bin = binx +(nbinsx+2)*biny;
      sumX  += histo->GetBinContent(binx,biny)*x;
      sumY  += histo->GetBinContent(binx,biny)*y;
      count += (Int_t)histo->GetBinContent(binx,biny);
      if(debugGC) printf("% 12.8lf \t % 12.8lf \t %i \n",sumX,sumY,count);
    }
  }

  if(count>0){
    val[0] = sumX/count;
    val[1] = sumY/count;
  }else{
    val[0] = sqrt(-1.0);
    val[1] = sqrt(-1.0);
  }

  //  cout << "HERE" << endl;
  // Loop again to get biased weighted SDOM error.
  Double_t sigmaX=0,sigmaY=0;
  for (biny=biny1;biny<=biny2;biny++) {
    Double_t y = yaxis->GetBinCenter(biny);
    for (binx=binx1;binx<=binx2;binx++) {
      Double_t x = xaxis->GetBinCenter(binx);
      // If we are in the cut and the bin count>0 add to sigmas.
      if (!cut->IsInside(x,y)) continue;
      bin = binx +(nbinsx+2)*biny;
      Int_t sCount =0;
      sCount = (Int_t)histo->GetBinContent(bin);
      if(sCount>0){
	sigmaX += sCount*(x-val[0])*(x-val[0]);
	sigmaY += sCount*(y-val[1])*(y-val[1]);
      }
    }
  }

  if(count>1){
    sigmaX = TMath::Sqrt(sigmaX/(count-1));
    sigmaY = TMath::Sqrt(sigmaY/(count-1));
  }else{
    sigmaX = TMath::Sqrt(-1.0);
    sigmaY = TMath::Sqrt(-1.0);
  }

  if(kPrint){
    printf("% 12.8lf \t% 12.8lf \t%i \t%12.8lf \t% 12.8f\n",val[0],val[1],count,sigmaX,sigmaY);
  }
  
  centroid[0] = val[0];
  centroid[1] = val[1];
  centroid[2] = count;
  centroid[3] = sigmaX;
  centroid[4] = sigmaY;
  
  return kTRUE;
}

//      if(!isnan(xRaw) && !isnan(yRaw)){
// 	Double_t step = 0.1;
// 	Double_t toler = 0.005;
// 	Double_t xx=-30,yy=30;
// 	Double_t xRawTmp, yRawTmp;
// 	while(xx<=30 && yy>=-30){
	  
// 	Int_t    nRowSplines=0;
// 	Double_t xTmp[13];
// 	Double_t yTmp[13];
// 	memset(xTmp,'0',sizeof(xTmp));
// 	memset(yTmp,'0',sizeof(yTmp));
// 	for(Int_t j=0; j<13; j++){
// 	  if(fRowX[j]){
// 	    xTmp[nRowSplines] = fRowX[j]->Eval(yy);
// 	    nRowSplines++;
// 	  }
// 	}
// 	fBicubicX = new TSpline3("BicubicX",GridX,xTmp,nRowSplines,"0",0,0);
// 	xRawTmp = fBicubicX->Eval(xx);
// 	nRowSplines = 0;
// 	for(Int_t j=0; j<13; j++){
// 	  if(fRowY[j]){
// 	    yTmp[nRowSplines] = fRowY[j]->Eval(xx);
// 	    nRowSplines++;
// 	  }
// 	}
// 	fBicubicY = new TSpline3("BicubicY",GridY,yTmp,nRowSplines,"0",0,0);
// 	yRawTmp = fBicubicY->Eval(yy);

// 	cout << xx << " \t" << yy << " \t" << xRawTmp << " \t" << yRawTmp << " \t" << xRaw << " \t" << yRaw 
// 	     << " \t" << Abs(xRawTmp-xRaw) << endl;

// 	if(Abs(xRawTmp-xRaw)<toler && Abs(yRawTmp-yRaw)<toler) {
// 	  x = xx;
// 	  y = yy;
// 	  // cout << x << " " << y << endl;
// 	  break;
// 	}

	
	
// 	delete fBicubicX;
// 	delete fBicubicY;
// 	xx += step;
// 	if(xx>29.95) {
// 	  yy = yy-step; 
// 	  xx = -30;
// 	  //	  cout << " HERE" << xx << " " << yy << " " << endl;
// 	}
// 	}
// 	//	cout <<endl;
// 	Char_t innn;
// 	//	cin >> innn;