#define FPTrack_cxx
#include "FPTrack.h"
#include "TCanvas.h"
#include "TLine.h"
#include "Getline.h"
#include "TRandom3.h"
#include "TF1.h"
#include <string>  
//------- Deal with the eccentricities of gcc 2.95.2  (From Ludovic)
//---------------NEW for pt ---------
#if defined(__GNUC__) && (__GNUC__ < 3)
#include <strstream>
#else
#include <sstream>
#endif
using std::cout;
using std::endl;
//------------------ look also in .h file

//--- Some histogram declarations
TCanvas * page2 = new TCanvas("page2"," TThadronic ",600,800);
TCanvas * page = new TCanvas("page"," TThadronic ",600,800);

TH1F * Xplane[51][2];    TH1F * ThisHist ; TH1F * Yplane[51][2] ; 
TH2F * XYplane[51][2];   TH2F * Xt[51][2] ; TH2F * XdX[51][2] ;
TH1F * Preco[3][2] ; TH1F * Mreco44;  TH1F * Mreco42;    
TH1F * UTAstyle;        TH1F * DXi      ; TH1F * xP[30] ;      
TH1F * ArcL[51][2] ;   
TH2F * DPvP[51][2] ;     TH2F * DPvx[51][2] ; TH2F * DPvy[51][2] ; 
TH1F * Magfail ;        
TH1F * yrap44 ;  TH1F * yrap42 ;  TH1F * yrap_all ; 
TH1F * XMax44 ;  TH1F * XMax4with2 ; TH1F * XMax2with4 ;
TH1F * YXMax44 ;  TH1F * YXMax4with2 ; TH1F * YXMax2with4 ;
TH1F * DelM44 ; TH1F * DelM42 ; TH1F * DelMall ;
TH1F * Delyall ; TH1F * DelAzi ; TH1F * DelAzi2 ; TH1F * DelPx[3][2];
TH1F * DelPy[3][2];TH1F * DelPt[3][2];
TH2F * y12[2] ; TH2F * x12[2]; TH2F * xy12[2];
TH1F * DelP[3][2];  TH1F * DXbyX[3][2]; 

TRandom3 gau ;
#define Nm 36
#define NH 50
#define NC 12
#define NY 40

// Magnet type:  0 = dipole, 1 = horiz focussing quad, 2 = vert focussing quad
// Magnet Strength is given by:
// B  = bend angle for dipole.  Obtained from beamline specs.  
//     Dipole assumed aligned with correct field for this.                  
// G  = quadrupole field gradient.  
// The quadrupole is assumed aligned in the relevant portion of beamline.

int    Magnet_type[Nm][2], Plane_nmag[51][2], Magnet_aperclass[Nm][2], 
       ncoll=0, Coll_nmag[NC][2], 
       symmetrise=0, ivec=0, nev, coltinue, masscan, icali ;
double Magnet_x[Nm][2],Magnet_y[Nm][2],Magnet_z[Nm][2],Magnet_length[Nm][2], 
       Magnet_strength[Nm][2], grad, Magnet_xb[Nm][2], arclen0[51][2], arclen, 
       Magnet_aperA[Nm][4][2],UTAmin,UTAdx,  
       Plane_z[51],Precoval[10][2], xplane0[51][2],yplane0[51][2], aplane0[51][2],
       ayplane0[51][2], Mrecoval , xmindis, xmaxdis, ymindis,ymaxdis,
       pi=3.141592654, twopi=6.283185308, delp=0.,
       Brho = 23349.499,   // this is 7000 GeV/0.3  
       pbeam0 = 7000.,  // Don't alter, Magnet formulae depend on this value
       sigbeam0 = 0.77 , // 1.1E-4 * 7000
       sigx0    = 0.0000118, // 16.8/sqrt(2) microns smearing at I P 
       siga0    = 0.21, //  0.21GeV from beam divergence in x and in y.
       sigdxdz = 0., sigx = 0., sigy=0.,
       apermb = 0. ,
       Si_edge = 0.,
       xcol1=999.,xcol2=999.,xcol3=999.,Coll_z[NC],Coll_xap[NC],xcoll0[NC][2],
       meanm, ximean,ximean2,mmean,mmean2,ymean,ymean2, ymmean, ptest;
float  cut[NH], ptval0[2], ptval[2], dptval, mscale=1000., // for milliradians 
       aside[2]={-1.,1.},
       Higgsmassscale = 1. ,  p0_subgev=0. , hima = 2.0,
  xinner0= 0.003, xinnerp[5][2], collhit,
       x0IP=0., y0IP=0., ax0IP=0., ay0IP=0., xmax[10][2], yxmax[10][2],
       pscan[50000][5][2], xscan[50000][5][2], axscan[50000][5][2],
       yscan[50000][5][2], ayscan[50000][5][2], pmom[2], p0[2],
       sumrvar42, sumrvar44, sumrvar22, sumvar42, sumvar44, sumvar22 ,
       massvec[NH],a42vec[NH],a44vec[NH], a22vec[NH], dm44vec[NH],dm42vec[NH], 
       dm22vec[NH], dmAllvec[NH],sigma_xi[NH], sigma_m[NH],sigma_y[NH],
       sigma_M44[NH],sigma_M42[NH], sigma_mall[NH],sigma_yall[NH],
       corr_my[NH], readplane[51], pxval,pyval,ptmean,azival[2],azival0[2]; 
double xdis = 0.097,  // half the beamline sepn. in the parallel region, COLD.
       centreIP = 0,
       beamangle =-0.0001425, // Angle of beam should be -142.5 murad
       lengthscale = 100.; // Convert metres into cm for plotting.
bool   initialising = false, docali=false, calibrating=false, 
       chrom=false,
       // -- ALTER HERE!
       useaper,drawhists, psmear=false,x0smear=false,a0smear=false, donexP=false,
       plotcoin = false,  initmess=true, initaccep=true, outputdat=false, keepreading ;
Int_t  ncand, nmagnet=0, ntype;
int    isay=1, irtype, nplane=0, iplane, naccep,naccep0,naccep1, nmean, UTAplane, 
       naccep00,naccep01,naccep2,naccep22,idetec[10][2],iaccep[200][10][2],
       ireadp;
int IP ; 

// These were calibrations up to April 2006
/*
  float cal_alpha[8]={-0.000636377,0.0168016,0.00082655, 0.016661,
                      -0.00621056, 0.01752, -0.00067558, 0.0174556};
  float cal_beta[8] ={ -0.0438394,  -0.00288172, -0.00350358,  -0.00544231, 
                       -0.0877092,  -0.0017645, 0.0153147, -0.00251031};
  float cal_a[8]={  60023.9, 3826.22, 57578.4, 4744.18, 
                    87021.3, 3851.32, 87068.5, 4724.73};
  float cal_b[8]={ -503759, 5239.91, -508822, 11199.6,  
                   -1.27205e+06, 4354.5, -1.3496e+06, 9029.42};
  float cal_c_0[8]={ -113378, 141421,   257963, -57956.5, 
                     25611, 62937.5, 541949, -163025};
  float cal_c_1[8] ={6605.83, -3.30494e+07, 6834.66, -1.85099e+07, 
           	   8823.45, -2.60729e+07, 8830.04, -1.0767e+07};
*/
//These calibrations apply using better plane placing. But less good
//IP1
/*
  float cal_alpha[8]={ -0.00466448, 0.0168316, 0.000826551, 0.0166905,
  		     -0.00621057, 0.01752, -0.000675582, 0.0174555};
  float cal_beta[8] ={-0.0712839, -0.00396075, -0.00350538, -0.00659849,
  		-0.0877053, -0.00176401, 0.0153153, -0.00251032};
  float cal_a[8]={ 57980.9, 3826.55, 57578.4, 4744.09,
  		 87021.4, 3851.33, 87068.7, 4724.73};
  float cal_b[8]={  -483181, 5213.83, -508822, 11353.5,
  		  -1.27207e+06, 4354.5, -1.34966e+06, 9029.01};
  float cal_c_0[8]={ -80202.9, 497087, 257966, -145843,
  		   25621.5, 62918.2, 541989, -163025};
  float cal_c_1[8] ={ 6708.58, -4.99871e+07, 6834.66, -1.36849e+07,
  		    8823.4, -2.60717e+07, 8829.81, -1.0767e+07};
*/

//Calibrations through summer 2006 amended.
/*
  float cal_alpha[8]={-0.00466448, 0.0168016,0.00082655, 0.016661,
                      -0.00621056, 0.01752, -0.00067558, 0.0174556};
  float cal_beta[8] ={ -0.0712839,  -0.00288172, -0.00350358,  -0.00544231, 
                       -0.0877092,  -0.0017645, 0.0153147, -0.00251031};
  float cal_gamma[8] ={0., 0., 0., 0., 0.,  0., 0., 0.};
  
  float cal_a[8]={  57980.9, 3826.22, 57578.4, 4744.18, 
                    87021.3, 3851.32, 87068.5, 4724.73};
  float cal_b[8]={ -483181, 5239.91, -508822, 11199.6,  
                   -1.27205e+06, 4354.5, -1.3496e+06, 9029.42};
  float cal_c[8] ={0., 0., 0., 0., 0.,  0., 0., 0.};
  float cal_c_0[8]={ -114740, 141421,   257963, -57956.5, 
                     25611, 62937.5, 541949, -163025};
  float cal_c_1[8] ={6916.8, -3.30494e+07, 6834.66, -1.85099e+07, 
           	   8823.45, -2.60729e+07, 8830.04, -1.0767e+07};
  float cal_c_2[8] ={0., 0., 0., 0., 0.,  0., 0., 0.};
  float cal_c2_1[8]={0., 0., 0., 0., 0.,  0., 0., 0.};
*/

//  Calibrations from final romp version.
 float cal_alpha[8]={ -0.00466905, 0.0168327, 0.000825495, 0.0166901,
                      -0.00621193, 0.0175206, -0.000676976, 0.0174554};
 float cal_beta[8] ={ -0.0692079, -0.00402188, -0.00300072, -0.00654559,
                      -0.0867181, -0.00179805, 0.0162908, -0.0024871};
 float cal_gamma[8] ={ -0.212972, 0, -0.0533164, 0,
                       -0.158438, 0, -0.15199, 0};
 float cal_a[8]={ 58050.8, 3825.3, 57653.7, 4744.7,
                  87151.1, 3850.38, 87210.2, 4725.18};
 float cal_b[8]={ -518976, 5287.21, -547203, 11270.8,
                  -1.37174e+06, 4410.24, -1.45922e+06, 8969.88};
 float cal_c[8]={ 3.9612e+06, 0, 4.21307e+06, 0,
                  1.64612e+07, 0, 1.81171e+07, 0};
 float cal_c_0[8] ={ -1.2417e+06, 138883, -324841, -42743.9,
                     -1.66331e+06, 121214, -317845, -62902}; 
 float cal_c_1[8] ={ 10054.3, -3.28456e+07, 10438.1, -1.88503e+07,
                    -1.66331e+06, 121214, -317845, -62902};
 float cal_c_2[8] ={ -33.6882, 0, -32.1523, 0,  -50.8771, 0, -50.7701, 0};
 float cal_c2_1[8]={ -205696, 0, -185161, 0,     -361955, 0, -346118, 0};


//  Calibrations using precise 6.5 optics file.
/*
float cal_alpha[8]={-0.00466916, 0.0168326, 0.000825148, 0.0166902,
                     -0.0062142, 0.0175223, -0.000674388, 0.0174571};
float cal_beta[8] ={ -0.0692103, -0.0040225, -0.00300311, -0.00654314,
                     -0.0866598, -0.00179245, 0.016336, -0.00247823};
float cal_gamma[8]={-0.212964, 0, -0.053308, 0,
                      -0.157977, 0, -0.152639, 0};
float cal_a[8]={ 58050.6, 3825.41, 57654.9, 4744.52,
                 87144.1, 3849.13, 87209.5, 4723};
float cal_b[8]={ -518970, 5287.84, -547229, 11268.7,
                 -1.37221e+06, 4402.74, -1.45997e+06, 8950.99};
float cal_c[8]={ 3.96104e+06, 0, 4.21363e+06, 0,
                  1.64732e+07, 0, 1.8137e+07, 0}; 
float cal_c_0[8]={ -1.24512e+06, 143292, -327934, -42214.7,
                   -1.66749e+06, 81229, -322183, -62580.5};
float cal_c_1[8]={ 10055.8, -3.31551e+07, 10438.7, -1.89274e+07,
                   14200.6, -2.58038e+07, 14851.6, -1.57679e+07};
float cal_c_2[8]={ -33.7, 0, -32.1575, 0,
                   -50.8734, 0, -50.7598, 0};
float cal_c2_1[8]={ -205847, 0, -185240, 0,
                   -361712, 0, -346231, 0};
*/

float calpx_alpha[8], calpx_beta[8], calpx_gamma[8],
      calpay_alpha[8], calpay_beta[8], calpay_gamma[8],
      calpx_cfit[8][3][3],  calpy_cyfit[8][3][3], calpay_cayfit[8][3][3],
      calpy_aa[8], calpy_bb[8], calpy_cc[8];

std::ofstream outputdatafile, debugfile, cfile ;
std::string thisIP;
bool message=true;
// ============================================================================
void FPTrack::SetTitles()
{
  int   nh    = 80;
  float ys    = 0.01*lengthscale; 
  float xs    = 0.025*lengthscale;
  // nh = 20; xs = 3.0 ; //  REMOVE

  Magfail     = new TH1F("outside+-","outside+-",nh,0.,80.);
  XMax44      = new TH1F("XMax44","XMax440+440",nh,0.,xs);
  XMax4with2  = new TH1F("XMax42","XMax440+220",nh,0.,xs);
  XMax2with4  = new TH1F("XMax24","XMax220+440",nh,0.,xs);
  YXMax44     = new TH1F("YXMax44","YXMax440+440",40,-ys,ys);
  YXMax4with2 = new TH1F("YXMax42","YXMax440+220",40,-ys,ys);
  YXMax2with4 = new TH1F("YXMax24","YXMax220+440",40,-ys,ys);
  yrap44      = new TH1F("yrap44","Rapidity 440+440",40,-2.,2.);
  yrap42      = new TH1F("yrap42","Rapidity 440+220",40,-2.,2.);

  // bool errorbars = false;
}
// *****************************************************************************
double FPTrack::sq(double x)
{ return x*x;}
float  FPTrack::sq(float x)
{ return x*x;}
// =============================================================================
float FPTrack::Pformula(int IP, int plane,int side, 
			float xdisp, float ydisp, float dirx, float diry)
{
  double pnew,disa,disx,pprev,c1calc,c2calc, precon, xval,aval, a,f0,f1,f2,
         disy,yval,asmear, xsmear, ysmear, aysmear, wy,way,py_y,py_ay,
         cycalc[3],caycalc[3];

  // -- Formula to reconstruct momentum from the silicon measurements
  // -- This is the "pomeron momentum", i.e. momentum lost by proton.


    // -- apply smearing if required ---
    asmear = 0. ; xsmear = 0.; ysmear = 0.;
    if(sigdxdz>0) asmear = sigdxdz*gau.Gaus();  
    if(sigx>0)    xsmear = sigx*gau.Gaus();  
    if(sigy>0)    ysmear = sigy*gau.Gaus();  
    xval = xdisp+ xsmear ;
    aval = dirx + asmear ;
    
    int index= IP-1 + plane-1 +2*side; 
    pnew  = cal_a[index]*xval + cal_b[index]*sq(xval) + cal_c[index]*xval*sq(xval) ; 
    disa  = cal_alpha[index]*xval  + cal_beta[index]*sq(xval)
           +cal_gamma[index]*xval*sq(xval);
    disa  = aval - disa;
    c1calc=0.; 
    pprev=pnew;
   
    //  Now apply the angular correction term.
    bool ang=true;  
    if(!ang) goto FINISH;    

    for(int iter=0; iter<4; iter++){
      c1calc =  cal_c_0[index]*fitf(plane,1,0,pnew) 
	      + cal_c_1[index]*fitf(plane,1,1,pnew)  
              + cal_c_2[index]*fitf(plane,1,2,pnew);
      c2calc =  cal_c2_1[index]*fitf(plane,2,1,pnew) ;
      pnew =  pprev - c1calc*disa -c2calc*sq(disa);
      }

// -- Now the px material.

     a = aval*1000. ;
     disx  = xval 
     - (calpx_alpha[index]*a + calpx_beta[index]*sq(a) +calpx_gamma[index]*sq(a)*a);
     f0   = fitpxf(plane,side,0,pnew) ; 
     f1   = fitpxf(plane,side,1,pnew) ; 
     f2   = fitpxf(plane,side,2,pnew) ; 
     c1calc =  calpx_cfit[index][1][0]*f0  + calpx_cfit[index][1][1]*f1 
             + calpx_cfit[index][1][2]*f2;
     c2calc = calpx_cfit[index][2][0]*f0  + calpx_cfit[index][2][1]*f1 
             + calpx_cfit[index][2][2]*f2;
     pxval =  c1calc*disx +c2calc*sq(disx);

// -- Now the py material.

     aysmear=ysmear*1.414/8.; //  8 metres between detectors
     yval = ydisp+ ysmear ;
     aval = diry + aysmear ;
     a = aval*1000. ;
     disa  = a - (calpay_alpha[index]*pnew + calpay_beta[index]*sq(pnew)
                                           +calpay_gamma[index]*sq(pnew)*pnew);
     disy  = yval - (calpy_aa[index]*pnew + calpy_bb[index]*sq(pnew) 
                                     +calpy_cc[index]*sq(pnew)*pnew);
     f0   = fitpyf(plane,side,0,pnew) ; 
     f1   = fitpyf(plane,side,1,pnew) ; 
     f2   = fitpyf(plane,side,2,pnew) ; 
     cycalc[1] = calpy_cyfit[index][1][0]*f0  + calpy_cyfit[index][1][1]*f1 
                                              + calpy_cyfit[index][1][2]*f2;
     caycalc[2] = calpy_cyfit[index][2][0]*f0 + calpy_cyfit[index][2][1]*f1 
                                              + calpy_cyfit[index][2][2]*f2;     
     caycalc[1] = calpay_cayfit[index][1][0]*f0  + calpay_cayfit[index][1][1]*f1 
                                                 + calpay_cayfit[index][1][2]*f2;
     caycalc[2] = calpay_cayfit[index][2][0]*f0  + calpay_cayfit[index][2][1]*f1 
                                                 + calpay_cayfit[index][2][2]*f2; 
     py_y =  (1./cycalc[1])*disy + (1./cycalc[2])*sq(disy);
     py_ay =  caycalc[1]*disa + caycalc[2]*sq(disa);
      wy = sq(cycalc[1]) ;
      way= 1000000./(32.*sq(caycalc[1])); //  (sqrt(2)/8 m)^2.   we worked here in mrad.
      pyval = (wy*py_y + way*py_ay)/(wy + way );


      // debugfile<<pnew<<" "<<side<<" "<<plane<<" "<<ptest<<" "<<pyval<<" "<<py_y<<" "<<py_ay<<endl;


 FINISH: // if(pnew<0.) pnew=0.;  // through upward fluct, this can occur.
         precon = pbeam0 - pnew ;

      return precon ;
}
// =======================================================================
void FPTrack::ReadData()
{
  //  Read data cards from Data.txt file 

  std::string dataname ;
  std::ifstream mydata;
  cout<<"Open data cards file"<<endl; 
  mydata.open("FPData.txt");
  if(chrom) cfile.open("chromdata");
  bool readon=true ;
  int item,i,j ;
  ireadp=0;

  IP = 0 ; 
  for(i=0;i<5;i++) {for(j=0;j<2;j++) xinnerp[i][j]=-1.;}
 
  while(readon)
   { 
     mydata>>dataname;
     if(dataname == "END") break ;
     mydata>>item;
     mydata.ignore(1000,'\n');
          if(dataname == "IP")   IP = item;
     else if(dataname == "ISAY") isay = item;
     else if(dataname == "RTYP") irtype = item;
     else if(dataname == "HIGG") Higgsmassscale = float(item)*0.01 ;
     else if(dataname == "PGEV") p0_subgev = float(item);
     else if(dataname == "SYMM") symmetrise = item ;
     else if(dataname == "PLNE") {ireadp++; readplane[ireadp] = 0.01*float(item);}
     else if(dataname == "SIMA") xinner0 = float(item)*0.000001;
     else if(dataname == "SIM1") {if(xinnerp[1][0] < 0.)    //  SM10,SM11 have precedence
                                 {xinnerp[1][0] = float(item)*0.000001;
                                  xinnerp[1][1] = float(item)*0.000001;}}
     else if(dataname == "SIM2") {xinnerp[2][0] = float(item)*0.000001;
                                  xinnerp[2][1] = float(item)*0.000001;}
     else if(dataname == "SM10")  xinnerp[1][0] = float(item)*0.000001;
     else if(dataname == "SM11")  xinnerp[1][1] = float(item)*0.000001;
     else if(dataname == "EDGE") Si_edge = float(item)*0.000001;
     else if(dataname == "APER") useaper = (item>0) ;
     else if(dataname == "MBAP") apermb = float(item)*0.001 ;
     else if(dataname == "HIST") drawhists = (item>0) ;
     else if(dataname == "ODAT") outputdat = (item>0) ;
     else if(dataname == "MSCA") masscan = item ;
     else if(dataname == "HIMA") hima = float(item*0.1) ;
     else if(dataname == "CALI") icali =  item ;
     else if(dataname == "PSIG") psmear = (item>0) ;
     else if(dataname == "DXIP") x0smear = (item>0) ;
     else if(dataname == "DAIP") a0smear = (item>0) ;
     else if(dataname == "YSIG") sigy = item*0.000001    ;
     else if(dataname == "ASIG") {sigdxdz = item*0.000001 ;
                    if(item>=100) sigdxdz*=0.001;}  //  *1000 for fractions of murad
     else if(dataname == "XSIG") sigx = item*0.000001    ;
     else if(dataname == "COL1") xcol1 = float(item)*0.0001;
     else if(dataname == "COL2") xcol2 = float(item)*0.0001;
     else if(dataname == "COL3") xcol3 = float(item)*0.0001;
     else if(dataname == "COLC") coltinue = item;
   }

  mydata.close() ;   

  cout<<"IP = "<<IP<<" Run type "<<irtype<<endl;
  cout<<" Use magnet apertures? "<<useaper<<endl;
  if(useaper && apermb>0) 
    cout<<" Set MB apertures to new value "<<apermb<<endl;  
  cout<<" Draw plots? "<<drawhists<<endl;  
  if(outputdat) outputdatafile.open("FPtraj.txt") ;
  debugfile.open("debugf.txt") ;
  if(outputdat) cout<<" Opened fptraj "<<endl;

  if(irtype<3) masscan=0 ;
  if(irtype==3) chrom=false;
  docali = icali>0 ;
  meanm = 120.*Higgsmassscale ;
  int nh = 100;  if(psmear || x0smear ) nh = 50 ;
  float ww = 15.;
  Mreco44     = new TH1F("M-reco44","M-reco44",nh,0.95*meanm,1.05*meanm);
  Mreco42     = new TH1F("M-reco42","M-reco42",nh,0.9*meanm,1.1*meanm);
  DelM44      = new TH1F("DelM44","DelM44",100,-ww,ww);
  DelM42      = new TH1F("DelM42","DelM42",100,-ww,ww);
  DelMall     = new TH1F("DelMall","DelMall",100,-ww,ww);
  Delyall     = new TH1F("Delyall","Delyall",100,-0.03,0.03);
  DelAzi      = new TH1F("DelAzi","DelAzi", 40,-2., 2.);
  DelAzi2     = new TH1F("DelAzi_cut","DelAzi_cut",40,-2., 2.);
}
// =======================================================================
void FPTrack::MagnetSet()
{
//   We count the magnets from 1  not zero !
//   Coordinate system: take z = forwards, x= to left , y = up.
//   For quadrupoles, input the quantity K1L.  type = 1,2 if K1L is +,-ve.
//   Input type = 3 here to let the program decide.
//   HKICK and VKICK are read and treated as bending magnets. 

  double endpos,maglength[Nm][2],magcen[Nm][2],magstrength[Nm][2],
    L,K0L,K1L,K2L,K3L, HKICK,VKICK, magxb[Nm][2],
         dx, dy, dist, magaper[Nm][4][2] ;
  int    magnet=0, maxmag=0, max1mag=0, minmag=0, magver, 
         side, magtype[Nm][2], magapertype[Nm][2] ;
  float  BETX,ALFX,MUX,DX,DPX,X,PX,BETY,ALFY,MUY,DY,DPY,Y,PY,A1,A2,A3,A4;

// -------- Arrange to read the magnet data from CERN files.
//          twiss_b1.txt and twiss_b2.txt for beam 1 and beam 2
//          beam 2 is on side=0 here and beam 1 is on side=1
//          Read through the files to pick up the requested magnets.

   std::string magname,magsort, aperture, basicmag ;
   std::ifstream MyMags ;

   magver = 1;   //  6.5 preliminary version
   magver = 2;   //  6.5 summer 06  version
   magver = 3;   //  6.503 summer 09  version

   magver=2;

   for(int ifile=0; ifile<2; ifile++)
    {

    if(magver==1)
      {
       if(IP == 1 && ifile==0)   centreIP = 26658.88320  ;
       if(IP == 1 && ifile==1)   centreIP = 0.;
       if(IP == 5 && ifile==0)   centreIP = 13329.593967 ;
       if(IP == 5 && ifile==1)   centreIP = 13329.289233 ;
       if(ifile==0) MyMags.open("twiss_65_b2_old.txt");
       if(ifile==1) MyMags.open("twiss_65_b1_old.txt");
      }
    if(magver==2)
      {
       if(IP == 1 && ifile==0)   centreIP = 13329.28923 ;
       if(IP == 1 && ifile==1)   centreIP = 13329.59397 ;
       if(IP == 5 && ifile==0)   centreIP = 13329.59397 ;
       if(IP == 5 && ifile==1)   centreIP = 13329.28923 ;
       if(ifile==0 && IP==1) MyMags.open("LHCB2IR1_6500");
       if(ifile==0 && IP==5) MyMags.open("LHCB2IR5_6500");
       if(ifile==1 && IP==1) MyMags.open("LHCB1IR1_6500");
       if(ifile==1 && IP==5) MyMags.open("LHCB1IR5_6500");
      }             //  Also version 6.500 at present (Sep 06)
    if(magver==3)
      {
       if(IP == 1 && ifile==0)   centreIP = 13329.28923 ;
       if(IP == 1 && ifile==1)   centreIP = 13329.59397 ;
       if(IP == 5 && ifile==0)   centreIP = 13329.59397 ;
       if(IP == 5 && ifile==1)   centreIP = 13329.28923 ;
       if(ifile==0 && IP==1) MyMags.open("v6503_IR1_twiss.thick.b2.tfs");
       if(ifile==0 && IP==5) {std::cout<<" NO OPTICS FILE"; abort();}
       if(ifile==1 && IP==1) MyMags.open("v6503_IR1_twiss.thick.b1.tfs");
       if(ifile==1 && IP==5) {std::cout<<" NO OPTICS FILE"; abort();}
      }             //  Version 6.503 at present (Sep 06)


    cout<<"Open magnets file for beam "<<ifile<<endl; 
// --- Skip introductory material, indicated by the 1st character of a record.
    while (MyMags.peek() =='@' || MyMags.peek() =='*' || MyMags.peek() =='$') 
       MyMags.ignore(1000,'\n');

    //  side 1 == direction to +ve s of IP  == beam 1.
    //  NB "dist" is -Delta(s) but magcen, i.e. z,  is then -dist.

    side = ifile ;
    maxmag = 35 ;
    if(side==0) magnet = maxmag;
    if(side==1) magnet = 1;
    bool readon=true;

    while(readon)
    {
      aperture = " ";
      if(magver==1)
      MyMags>>magname>>magsort>>endpos>>L>>K0L>>K1L>>K2L>>K3L>>HKICK>>VKICK
	    >>BETX>>ALFX>>MUX>>DX>>DPX>>X>>PX>>BETY>>ALFY>>MUY>>DY>>DPY>>Y>>PY
            >>aperture>>A1>>A2>>A3>>A4;
      if(magver==2)
	MyMags>>magname>>magsort>>basicmag>>endpos>>L>>HKICK>>VKICK>>K0L
              >>K1L>>K2L>>K3L>>X>>Y>>PX>>PY>>BETX>>BETY>>ALFX>>ALFY
              >>MUX>>MUY>>DX>>DY>>DPX>>DPY>>aperture>>A1>>A2>>A3>>A4;
      if(magver==3)
	MyMags>>magname>>magsort>>basicmag>>endpos>>L>>HKICK>>VKICK>>K0L
              >>K1L>>K2L>>K3L>>X>>Y>>PX>>PY>>BETX>>BETY>>ALFX>>ALFY
              >>MUX>>MUY>>DX>>DY>>DPX>>DPY>>A1>>A2>>A3>>A4;

      MyMags.ignore(1000,'\n') ;             //  skip rest of line if necessary
 
      if( (magname == "\"IP1\"" && IP==1)  ||      
          (magname == "\"IP5\"" && IP==5))       //  this is the IP
        {
	 if(magver==1)
	    {PY = PY * 142.5/142. ; PX = PX * 142.5/142. ;}
	                   // ** a fudge because the file is only to 3 sig fig. 
          y0IP = Y; ay0IP = PY;   
          x0IP = X; ax0IP = PX;
          if(IP==1) thisIP= "IP1"; else thisIP = "IP5";
       }
      if(magver==1)       //  Some more fudges to give more precision
	{
        if(fabs(K0L)==0.000188) K0L=K0L*188.175/188.0  ;   
        if(fabs(K0L)==0.001129) K0L=K0L*1129.049/1129. ;   
        if(fabs(K1L)==0.055577) K0L=K0L*55576.561/55577. ; 
        if(fabs(K1L)==0.047986) K0L=K0L*47986.041/47986. ;
        if(fabs(VKICK)==0.00020) VKICK=VKICK*20.171/20. ; 
        if(fabs(VKICK)==0.00019) VKICK=VKICK*19.17/19. ; 
        }

      if(magver==3)
        {if(ifile==0) K0L = - K0L ; } // fudge for 6.503 bends in beam 2 *!*!

      dist = centreIP - endpos + 0.5*L ;
      float adist=fabs(dist) ;

     if(adist > 437. && side==0)  continue ;     
     if(side==1 && dist > 0.)     continue ;
     if(side==0 && dist <= 0.)     break    ;
     if(adist > 437. &&  side==1)  break    ;  
     if(side==0) minmag = magnet;
 
     magtype[magnet][side] = -1 ;

     if     (K0L != 0.)  { magtype[magnet][side]     =  0; 
                           if(L>10.) magstrength[magnet][side] =  aside[side]*K0L; 
                           else      magstrength[magnet][side] =  aside[side]*K0L; }
     else if(K1L != 0.)  { magtype[magnet][side]     = 3; 
                           magstrength[magnet][side] =  K1L; }
     else if(HKICK != 0.){ magtype[magnet][side]     = 0; 
                           magstrength[magnet][side] =  -HKICK; }
     else if(VKICK != 0.){ magtype[magnet][side]     = 4; 
                           magstrength[magnet][side] =  VKICK;} 

     if(magtype[magnet][side] < 0) continue ;  //  --(not a magnet) 

     maglength[magnet][side] = L   ;
     magcen[magnet][side]    = -dist;

     cout<<magname<<" "<<endpos<<" "<<L<<" "<<K0L<<" "<<K1L<<" "
              <<magnet<<" --> "<<magcen[magnet][side]<<" "<<maglength[magnet][side]
              <<" "<<magstrength[magnet][side]<<endl; 

     //  CLassify apertures.
          if(aperture == "\"NONE\"") magapertype[magnet][side] = 0 ;
     else if(aperture == "\"CIRCLE\"" && A1 == 0.) magapertype[magnet][side] =0;
     else if(aperture == "\"CIRCLE\"" && A1 >  0.) magapertype[magnet][side] =1;
     else if(aperture == "\"RECTELLIPSE\"")        magapertype[magnet][side] =2;
     else if(magver==3 && A1==0) magapertype[magnet][side] =0;  
     else if(magver==3 && A1!=0) magapertype[magnet][side] =2;  

     else cout<<" Unknown magnet aperture type "<<aperture<<endl ;

     if(apermb > 0)      //  -------alter the MB apertures
       { int pos = magname.find("MB.");
         if(pos==1){ A3=apermb; A4=apermb; }
       }

     magaper[magnet][0][side] = A1 ;
     magaper[magnet][1][side] = A2 ;
     magaper[magnet][2][side] = A3 ;
     magaper[magnet][3][side] = A4 ;
     magxb[magnet][side]      = X  ;

     if(side==1)  magnet++;
     if(side==0)  magnet-- ; 
     if(magnet>maxmag || magnet<1 ) 
        {cout<<"INCREASE MAXMAG VALUE"<<endl; break;}
    }
    MyMags.close();
  }

  //  Load up the settings.  "magnet" contains final side=1 setting + 1

  for (int is=0; is<2; is++) 
   {
    int i =is;     if(symmetrise==1) i = 0  ; 
    int j =1 ;     if(symmetrise==1) j = is ; 
    if(i == 1) minmag = 0 ;
    if(i == 0) max1mag = maxmag - minmag ;     
    if(i == 1) max1mag = magnet - 1 ;     
    for (int mag=1; mag<=max1mag; mag++)
      {
      dx = 0.;
      dy = 0.;

      if(fabs(magcen[minmag+mag][i]) > 155. ) dx = -xdis;     //  note the sign convention
      
      // CR SHIFT!!!
      int mz = int(fabs(magcen[minmag+mag][i]));
      if(mz ==228 ) 
	{
           dx = dx+0.001;   
        }

      MagnetPlace(mag,is,magtype[minmag+mag][i], dx,dy,
		  aside[j]*(magcen[minmag+mag][i]), magxb[minmag+mag][i],
                         maglength[minmag+mag][i],
		  // -aside[i]*magstrength[minmag+mag][i]) ;
                     magstrength[minmag+mag][i]) ;
      Magnet_aperclass[mag][is] = magapertype[minmag+mag][i] ;
      for(int k=0; k<4; k++)  Magnet_aperA[mag][k][is] = magaper[minmag+mag][k][i] ;
      }
      cout<<endl;
   }
      nmagnet=35 ;
}
// =======================================================================
void FPTrack::PlaneSet()
{
  float siliconwidth = 0.025;
  // float xmar    = 0.001;
  float xclose  = -0.001 ;             // allow a margin the other side
  float xfar ;    

  nplane  = 0;
  ymindis = -0.2*siliconwidth;
  ymaxdis = +0.2*siliconwidth;
  if(!initialising) {xclose = xplane0[0][0]; 
                     if(xplane0[0][1] < xclose) xclose = xplane0[0][1] ;
		                            xclose = xclose - 0.001;}
 
  // This all needs tidying up if many planes are used.

  xfar= xclose + siliconwidth;
  if(xclose < xfar){ xmindis = xclose; xmaxdis = xfar;} 
  if(xclose > xfar){ xmindis = xfar  ; xmaxdis = xclose;} 
  if(ireadp==0) 
    {PlanePlace(220., xmindis, xmaxdis, ymindis, ymaxdis) ;
     PlanePlace(420., xmindis, xmaxdis, ymindis, ymaxdis) ;
    }
  else for(int i = 1; i<=ireadp; i++) 
      PlanePlace(readplane[i], xmindis, xmaxdis, ymindis, ymaxdis) ;
}
// ===========================================================================
void FPTrack::CollSet()
{
  // Place an effective collimator at each end of the nominal one. 
  float hlength = 0.25;   // half the length of 0.5m
  CollPlace(150.345 - hlength, xcol1) ;
  CollPlace(150.345 + hlength, xcol1) ;
  CollPlace(184.857 - hlength, xcol2) ;
  CollPlace(184.857 + hlength, xcol2) ;
  CollPlace(224.8   - hlength, xcol3) ;
  CollPlace(224.8   + hlength, xcol3) ;
}
// ===========================================================================

void FPTrack::Run()
{
// ---------- Steering routine for the analysis. ------------    
   float h1=1.,h2=1.,h3=1.;

   ReadData() ;
   SetTitles();
   MagnetSet();
   if(irtype<1) return;
   PlaneSet();
   CollSet() ;

   // -- Run a single track to calibrate beamline  
   if(irtype > 1) { cout<<"Initialising"<<endl;
                   initialising = true ;
                   Tracks(1) ;
                   // PlaneSet();  
                   initialising = false ;
                  }   

   if(isay>0) debugfile<<" Start calibration "<<std::endl;
   // -- Calibrate the momentum measurement if asked for.
   if(irtype ==3 && docali) { cout<<"Calibrating"<<endl;
                   calibrating = true ; 
                   Tracks(2) ;
                   calibrating = false ;
                  }       
   if(isay>0) debugfile<<" Start tracking "<<std::endl;
   // -- Do tracking for input data file. Possible mass scan.
   if(irtype!=3) masscan = 0 ;
   if(irtype==2) calibrating = docali;
   if(masscan==0 || masscan>=50) Tracks(irtype) ;
   else 
    { if(masscan==1){h1=0.4; h2=hima; h3=0.1;}
      if(masscan==2){h1=2.0; h2=6.0; h3=0.25;} 
      // if(masscan==2){h1=6.6; h2=16.6; h3=0.6;}
      // ----- Types of mass scan ---------------
      for(Higgsmassscale = h1; 
            Higgsmassscale<= h2; 
              Higgsmassscale = Higgsmassscale + h3)
         {initaccep = true;
          DelM42->Reset(); DelM44->Reset() ; 
          DelMall->Reset(); Delyall->Reset();
          Tracks(irtype) ; 
          ivec++ ;
         }
     for(int ip=1; ip<3; ip++)
     for(int side=0; side<2; side++)
     {
     DelP[ip][side]->Fit("gaus","Q");
     TF1 *fit42 = DelP[ip][side]->GetFunction("gaus");
     float par2 = fit42->GetParameter(2);
     cout<<"Plane "<<ip<<" side "<<side<<" Presn "<<par2<<endl ;
     }

       cout<<"Double_t mass["<<ivec<<"]={";
       for(int i=0; i<ivec; i++) {cout<<massvec[i];
	 if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}
       if(masscan!=2){
          cout<<"Double_t a44["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<a44vec[i];
	  if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}}
       cout<<"Double_t a42["<<ivec<<"]={"; 
         for(int i=0; i<ivec; i++) {cout<<a42vec[i];
	 if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}
       cout<<"Double_t a22["<<ivec<<"]={"; 
         for(int i=0; i<ivec; i++) {cout<<a22vec[i];
	 if(i<ivec-1) cout<<", "; else cout<<"};"<<endl; }
       cout<<"Double_t dm42["<<ivec<<"]={"; 
         for(int i=0; i<ivec; i++) {cout<<dm42vec[i];
	 if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}
       if(masscan!=2){
          cout<<"Double_t dm44["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<dm44vec[i];
          if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}}
       cout<<"Double_t dm22["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<dm22vec[i];
	  if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}
       if(masscan!=2){
          cout<<"Double_t sigma_M44["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<sigma_M44[i];
	  if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}}
       cout<<"Double_t sigma_M42["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<sigma_M42[i];
	  if(i<ivec-1) cout<<", "; else cout<<"};"<<endl; }
       cout<<"Double_t sigma_mall["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<sigma_mall[i];
          if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;}
       cout<<"Double_t sigma_yall["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<sigma_yall[i];
          if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;} 
       cout<<"Double_t corr_my["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) {cout<<corr_my[i];
          if(i<ivec-1) cout<<", "; else cout<<"};"<<endl;} 
       cout<<"Xp error widths "<<endl;
       for (int i=0; i<30; i++)
        {
        Double_t par2;
        xP[i]->Fit("gaus","Q");
        TF1 *fitx = xP[i]->GetFunction("gaus");
        par2 = fitx->GetParameter(2);
        cout<<" "<<par2;}
        cout<<endl;
     }
   cout<<" Smear beam momentum? "<<psmear
           <<" Smear x at interaction point? "<<x0smear
	   <<endl<<" sig-x = "<<sigx<<" m. "
	   <<" sig-angle = "<<sigdxdz<<" rad. "<<endl;
}
// ======================================================================
void FPTrack::Tracks(int iruncon) 
{
  int    numevents=0,n=0,nnn,i_p=0,nptx=0,npty=0,iptx=0,ipty=0, 
         ipchoice,sidechoice,
         iphi=0,ipt=0,imom=0,ndelp=0,nthphi=0,npt=0,nphi=0;
  double qk,qkl,qkc,xe=0,xae=0,ye=0,yae=0,bend,bend0,anglemax=0.,ptx,pty, 
         p=0.,tandir,xc,yc,zc,dxyz[3],dir[3], xa1,x1,ya1,
         dzc,xp,yp,zp, dx0ip,dy0ip,dpx0ip,dpy0ip, theta=0., phi=0.;
  double ppx[2],ppy[2],ppz[2], ppzin[2], xangle, yangle, s=0., zposition,
    ppx0[2],ppy0[2], ptset=0., pxlast, pylast;
  float  dphi=0., xdisp, ydisp, themr, xip,thip,xamur; 
  float  xposition, yposition, xwrtbeam=-999., ptscale, xshift=0.;
  bool   outside, coutside, chromprin, accepted; 

  double radius =0.2; //  this default beampipe radius should be properly set 

  p0[0] = pbeam0 ; p0[1] = pbeam0 ;  pxlast=-99999.;pylast=99999.;
   
//  Define the events for iruncon<=2 case.
   if(iruncon==1) numevents=1;
   if(iruncon==2)                   //  calibration and scans
     {
     nphi     = 1;
     npt      = 29;
     nptx     = npt;       //  should be no more than 31
     npty     = npt ;
     if(npt>31) std::cout<<" CALIB ERROR npt is "<<npt<<std::endl; 
     delp     = 20. ; 
     // delp    = 10. ; // ================ TEMP for plane 2
     if(masscan==1) delp = 20. ;   
     if(masscan==2) delp = 20.; 
     ndelp    = 50;       //  step in longitudinal momentum
     docali   = true;
     // anglemax = 0.0004;   //  This gives +-200 murad
     anglemax = 0.004;   //  This gives +-2000 murad  FOR CHR. ROYON.

     if(chrom) {ndelp=11; nptx=50; anglemax=0.0005; delp=10.;} // 10 or 100

     dptval  = pbeam0*anglemax/float(npt);  
     if(chrom) {npty=1; nptx=2*nptx +1;}  
                            //  Calibration needs p:ptx grid only
     nthphi  = nptx*npty;
     dphi    = twopi/nphi;
     numevents=nthphi*ndelp;

     if(numevents >50000) {cout<<numevents<<" PSCAN INDEXING TOO SMALL!"<<endl;
                          abort();}
     for(n = 0; n < numevents; n++){
     for(int ip = 1; ip <= nplane; ip++){ 
         pscan[n][ip][0]    = 0. ; pscan[n][ip][1]    = 0. ;
         axscan[n][ip][0]   = 0. ; axscan[n][ip][1]   = 0. ;
       }}
     }   

   if(iruncon==4 || iruncon==5)
     {            // -->acceptances for Marta
     dptval  = 0.2;
      npt     = 10;     
     nphi    = 100;
     nthphi  = nphi*npt;
     ndelp   = 200;    
     delp    = 0.001*pbeam0 ;
     if(iruncon==5) delp    = 0.0001*pbeam0 ;
     dphi    = twopi/nphi;
     numevents=nthphi*ndelp;
     if(iruncon==1) numevents=1;
     // -------- iaccept operates over Del(p) and Del(pt) bins for given plane
     //          (iruncon=4 -> 220m, iruncon=5 -> 420m) 
     for(n = 0; n < ndelp; n++){
     for(int ipt = 0; ipt < npt; ipt++){ 
         iaccep[n][ipt][0]    = 0 ;
         iaccep[n][ipt][1]    = 0 ;
       }}
     }   
// --------- Prepare to read input file for iruncon=3 case
   std::ifstream MyInput;
   if(iruncon == 3)
     { 
     if(initmess) cout<<"Open input data file -"; 
     if(masscan<5)  MyInput.open("p1p2data");
     if(masscan==11)  MyInput.open("PileUpData");
     if(masscan==12)  MyInput.open("DiffData1");   // Prague from Vojta
     if(masscan==41)  MyInput.open("PileUpData");
     if(masscan==54) MyInput.open("lpair_mumupt4.dat");
     if(masscan==56) MyInput.open("lpair_mumupt6.dat");
     if(MyInput.eof()) cout<<" Empty Input File!"<<endl ;
     else if(initmess) cout<<" done. "<<endl;
     numevents = 100000000 ;
     initmess = false ;
     }
   
   bool p0smear = !initialising && !calibrating && psmear && iruncon>2;
   bool ipsmear = !initialising && !calibrating && (x0smear || a0smear) 
                  && iruncon && 2;

// ------- EVENT LOOP --------------------------------
// ------- Process events according to chosen method:-
   for (nev=0; nev <numevents; nev++)
     {

     accepted = false; 
     if(nev==1000 && !calibrating) isay=0; //  stop writing
     if(nev > 0) initialising = false ; //  first event did initialising.
     if(isay>0 && iruncon==3) debugfile<<"Event "<<nev<<" "<<endl; 
     if(p0smear) {
                  p0[0] = pbeam0 + gau.Gaus()*sigbeam0 ;
                  p0[1] = pbeam0 + gau.Gaus()*sigbeam0 ;
                 } 
     dx0ip = 0.; dy0ip = 0.; dpx0ip=0.; dpy0ip=0.;
     if(ipsmear) 
       { if(x0smear) dx0ip  =  sigx0*gau.Gaus() ;  //  common for both sides
         if(a0smear) dpx0ip =  siga0*gau.Gaus() ;  //  common for both sides
         if(a0smear) dpy0ip =  siga0*gau.Gaus() ;  //  common for both sides
         if(x0smear) dy0ip  =  sigx0*gau.Gaus() ;  //  common for both sides
       }
     if(iruncon == 1) 
       {
       pmom[0] = p0[0] - p0_subgev;
       pmom[1] = p0[1] - p0_subgev ; 
       debugfile<<"Beam momentum "<<pmom[0]<<endl; 
       theta=0;
       phi=0;
       }
     else if(iruncon == 2 || iruncon == 4 || iruncon == 5)
       {
        if(iruncon ==2 )
        {nnn  =  nev/npty;
         i_p   =  nnn/nptx; 
         iptx = nnn - nptx*i_p;
         ipty =  nev - npty*nnn; 
         ptx = (float(nptx)*0.5  - 0.5 - float(iptx))*dptval;
         pty = (float(npty)*0.5  - 0.5 - float(ipty))*dptval;
         phi = atan2(pty,ptx);
         ptset = sqrt(pty*pty + ptx*ptx);
         ipt = iptx ;
	         if(chrom) i_p=i_p-1;
        }
        else  
       {
        nnn  =  nev/nphi;
        i_p   =  nnn/npt; 
        ipt = nnn - npt*i_p;
        iphi =  nev - nphi*nnn; 
        phi     = dphi*iphi  ;
        ptset   = dptval*ipt ; 
        nnn     = nev/nthphi + ipt ;
        if(ipt == nphi/2) nnn=nnn+npt*ndelp ;
       }
       if(calibrating) 
             {pmom[0] = p0[0] - pgrid(i_p,delp) ;
	      pmom[1] = p0[1] - pgrid(i_p,delp) ; }
        else {pmom[0]= p0[0] -  (i_p+1)*delp;
	      pmom[1]= p0[1] -  (i_p+1)*delp;
             }
       }
       else if(iruncon == 3)
       { 
       if(MyInput.eof()) break ;
       if(masscan<12)
           MyInput >>ppx0[0]>>ppy0[0]>>ppzin[0]>>ppx0[1]>>ppy0[1]>>ppzin[1];
       else if(masscan==12) //  Prague diffractive data
	 {keepreading=true;
          while(keepreading)   //  Many irrelevant lines.
	    {
            if(MyInput.eof()) break ;
            MyInput >>ppx0[0]>>ppy0[0]>>ppzin[0];
            MyInput.ignore(1000,'\n') ;             //  skip rest of line if necessary
            if(ppx0[0]==0. || ppy0[0]==0.) continue;
            if((ppx0[0]==pxlast) || (ppy0[0]==pylast)) continue;  // repeated lines
            ppzin[0] = fabs(ppzin[0]); 
            pxlast=ppx0[0]; pylast=ppy0[0];
            if(ppzin[0] < 6000000.) continue;   //  it lost too much energy
            ppzin[1] = -ppzin[0]; ppx0[1]= -ppx0[0]; ppy0[1]= -ppy0[0];
            keepreading = false;
            }
	    if(keepreading && MyInput.eof()) break ;  //  final event was unusable 
            for(int iii=0; iii<2; iii++) {ppx0[iii]=0.001*ppx0[iii];
	      ppy0[iii]=ppy0[iii]*0.001; ppzin[iii]=ppzin[iii]*0.001;}  //  MeV -> GeV
	 } 
	else if(masscan>=40) //  forget why it was done like this    
	 {for(int skipline=0;skipline<8;skipline++) MyInput.ignore(1000,'\n');
	   MyInput>>ppx0[0]>>ppy0[0]>>ppzin[0];
           for(int skipline=0;skipline<6;skipline++) MyInput.ignore(1000,'\n');
	   MyInput>>ppx0[1]>>ppy0[1]>>ppzin[1];
           for(int skipline=0;skipline<8;skipline++) MyInput.ignore(1000,'\n');
         } 
       //    LPAIR ascii file.
       if(isay>1) debugfile<<"Input "<<ppx0[0]<<" "<<ppy0[0]<<" "<<ppzin[0]
                          <<" "<<ppx0[1]<<" "<<ppy0[1]<<" "<<ppzin[1]<<endl; 
       }  //  iruncon==3 section   
// --- Rescale the mass.  But we want to keep PT the same.
       ptscale = 1.0;
       for(int side=0; side<2; side++)
            {ppx[side]=ppx0[side]*ptscale + dpx0ip;
             ppy[side]=ppy0[side]*ptscale + dpy0ip;
             ppz[side]=p0[side]-(pbeam0-fabs(ppzin[side]))*Higgsmassscale ;
	    } 
// ----- Loop over the two sides
     for(int side = 0; side<2; side++)
      { 
	 Precoval[1][side] = pbeam0;   // Default.  Adjust argument 1 if > 1 track.
       //  set initial position, angle, momentum.  This x is positive inwards! 

       if(iruncon == 2 || iruncon == 4 || iruncon == 5)
	 {theta      = ptset/pmom[side] ; 
	   if(chrom) {theta = anglemax - float(iptx)*anglemax/float(npt);
	     phi = 0.;}
         }
       if(iruncon == 3 && !initialising)
         { 
         ptval0[side]     = sqrt( sq(ppx0[side]) + sq(ppy0[side])) ;
         ptval[side]      = sqrt( sq(ppx[side]) + sq(ppy[side])) ;
         pmom[side] = sqrt( sq(ptval[side])     + sq(ppz[side])) ; 
         phi        = atan2(ppy[side],ppx[side]);  
         theta      = ptval[side]/pmom[side] ;  
         }   
       idetec[1][side]=0; idetec[2][side]=0;
       xp = -x0IP  + dx0ip;   //  second term can allow for beam spot size.
       yp =  y0IP  + dy0ip;
       zp =  0. ;

       if(iruncon==3 && isay>0) 
       debugfile<<" Input "<<pmom[side]<<" "<<iruncon<<" "<<initialising<<std::endl;

       xmax[1][side] = 0. ;  xmax[2][side] = 0. ;  arclen = 0.;
                          //  NB definitions. These are dx/ds etc not dx/dz.
       dir[0]  = sin(theta)*cos(phi);   
       dir[1]  = sin(theta)*sin(phi);          
       dir[2]  = aside[side]*cos(theta); 
       p       = pmom[side];

       //  Rotate by beam rotation angle.  A very small horiz. rotation (IP5)
       //                                  or vert. (IP1) Set from optics file.
       if(! (chrom && iruncon ==2 && IP==5))  //  IP1 use correct angle
       { dir[0]  = dir[0] - ax0IP ;
         dir[1]  = dir[1] + ay0IP ; }
       // if(chrom) {xp=0.; yp=0.;}for test comparisons.

       //  store the initial parameters
       xip = xp; thip = 1000000.*dir[0]; 

       //  dir[1] =  - ay0IP ;  //  to test out the reverse beam focussing.  


   // **** LOOP OVER THE BEAM LINE ITEMS.
   //          Magnet labels start at 1.

       outside = false ;  coutside  = false;
       collhit = 99999.;  
       for(int magnet = 1; magnet <= nmagnet; magnet++)
         { 
	 //  can by hand switch off magnet apertures:
	 // Magnet_aperclass[23][side] = 0 ;
	 // Magnet_aperclass[24][side] = 0 ;
	 // Magnet_aperclass[26][0] = 0 ;
	 // Magnet_aperclass[28][1] = 0 ;
	 // Magnet_aperclass[27][side] = 0 ;
	 // Magnet_aperclass[30][side] = 0 ;
	 // Magnet_aperclass[31][side] = 0 ;

	 // -- track to next magnet
         // -- track through next magnet if not the last one. 
	   if(fabs(zp) > 430.) break ;      //  Terminate when far enough.

	 zc = Magnet_z[magnet][side] ;
	 yc = Magnet_y[magnet][side] ;
	 xc = Magnet_x[magnet][side] ;

  //   Get positions and angles.  We have xp,yp,zp wrt Hall. 
  //   Take relative to centre of magnet aperture.
	 dxyz[0] = xp - xc ;
	 dxyz[1] = yp - yc ;
         dxyz[2] = zp - zc ; 
          if(isay>2 && magnet == 1) debugfile<<magnet<<" IP "
	     <<dxyz[0]<<" "<<dxyz[1]<<" "<<dir[0]<<" "<<dir[1]<<" p="<<p<<endl;

         if(isay>2)debugfile<<magnet<<" IN "<<dxyz[0]<<" "<<dxyz[1]<<" "<<zp<<endl;
  //  take trajectory  to the start of this magnet.
  
        dzc     = -0.5*Magnet_length[magnet][side]*aside[side] - dxyz[2] ;
        dxyz[0] = dxyz[0] + dzc*dir[0]/dir[2] ;
        dxyz[1] = dxyz[1] + dzc*dir[1]/dir[2] ;
        arclen  = arclen  + aside[side]*dzc/dir[2] ;   // added aside
        xa1 = dir[0]; ya1 = dir[1];

// std::cout<<" A " <<side<<" "<<zp<<" "<<arclen<<std::endl;

	// Are we still in the beampipe? 
        if(Magnet_aperclass[magnet][side] == 1) 
           { radius  = Magnet_aperA[magnet][0][side] ;
	     outside = outside || sq(dxyz[0])  +sq(dxyz[1]) > sq(radius);}
        if(Magnet_aperclass[magnet][side] == 2) 
           {
	    outside = outside ||
                       sq(dxyz[0]/Magnet_aperA[magnet][2][side])  
	            +  sq(dxyz[1]/Magnet_aperA[magnet][3][side])  > 1. ;  
            for(int k = 0; k<2; k++) if(Magnet_aperA[magnet][k][side] > 0.) 
              outside = outside || fabs(dxyz[k]) > Magnet_aperA[magnet][k][side] ; 
	   }

        if(isay>2) debugfile<<magnet <<" EN "<<dxyz[0]<<" "
                            <<dxyz[1]<<" "<<zc+dzc+dxyz[2]<<endl;
        if(isay > 1 && outside)
          {debugfile<<" - outside - "<<magnet<<endl;
	   if(isay > 2)
	    debugfile<<dxyz[0]<<" "<<Magnet_aperA[magnet][2][side]<<" "<<dxyz[1]<<" "
	    <<Magnet_aperA[magnet][3][side]<<" "<<Magnet_aperA[magnet][0][side]
	    <<" "<<Magnet_aperA[magnet][1][side]<<endl;
	  } 

        if(useaper && outside) {Magfail->Fill(float(magnet+40*side)); break;} 

//  -- Find what kind of element it is and take appropriate action. 

	if(Magnet_type[magnet][side] == 0 ||  Magnet_type[magnet][side] == 4)
	 { 
	  int i = 0;                                   //  -- Horiz BEND
          if( Magnet_type[magnet][side] == 4) i=1 ;    //  -- Vert BEND        

//  --                       DIPOLE MAGNET

	  //  tandir is in x,y,z bearing in mind the sign of z/
          //  But the magnet strength is stored as its effect on dx/ds or dy/ds.
          //  bend =  magnet strength, corrected for momentum, = bend for this proton. 
          //  dzc  =  length of magnet, signed
	  //  dxyz[i] =  New displ. = old + effect of gradient + effect of bend.
          //  dxyz[i-1] = old displacement + effect of gradient in non-bending plane.

          x1      =  dxyz[0];
          tandir  =  dir[i]/dir[2] ;
          bend    =  Magnet_strength[magnet][side]*pbeam0/p;
          dzc     =  Magnet_length[magnet][side]*aside[side] ;
          dxyz[i] =  dxyz[i] + dzc*tandir +(dzc/aside[side])*tan(0.5*bend) ;
          dxyz[1-i] =  dxyz[1-i] + dzc*dir[1-i]/dir[2];
          dir[i]  =  (tandir*aside[side] +tan(bend))/(1.-aside[side]*tandir*tan(bend));
          dir[i]  =  dir[i]/sqrt(1. + sq(dir[i]))  ;

          //  tan(new angle) = tan(sum of orig. angle + angle of bend)
          //  Then we turn it into the sine, e.g. dx/ds (perhaps unncessarily)

          if(isay>1 && i==0) debugfile<<magnet<<" DE " ;
          if(isay>1 && i==1) debugfile<<magnet<<" VE " ;

	  //   If main bending magnet, the local coord system must be redefined.
          //   Rotate wrt the local lab coords if necessary, and return to lab frame.
          //   The main bends are curved, so we must do this before evaluating the aperture.
         if(Magnet_length[magnet][side] > 10.) //  Identifies MB magnet.
          {
            bend0   =  Magnet_strength[magnet][side] ;
            dir[0]  = dir[0]  - bend0;         //  Rotates the coords.
	    dxyz[0] = dxyz[0] - 0.5*bend0*Magnet_length[magnet][side] ;
                                               //  Shifts coord system laterally.
            dxyz[2] = dxyz[2] - dxyz[0]*bend0*aside[side] ;
            arclen  = arclen + 0.5*(x1 + dxyz[0])*bend0*aside[side] ;  // added aside

	  }
        }
	else
//                                          QUADRUPOLE        
	{   
	     qk  = 0.2997925* Magnet_strength[magnet][side]/p ;
             //                                    ---- basic formula is m Te GeV 
	     qk = sqrt(qk);
	     qkl = qk * Magnet_length[magnet][side]*aside[side]  ;       
	     qkc = qk *dir[2] ;

	     if(Magnet_type[magnet][side] == 1) 
	     {
              //                                  ----- horizontally focussing 
              xe  =     cos(qkl)* dxyz[0] + sin(qkl)*dir[0]/qkc ; 
              xae = -qk*sin(qkl)* dxyz[0] + cos(qkl)*dir[0]/dir[2] ; 
              ye  =    cosh(qkl)* dxyz[1] + sinh(qkl)*dir[1]/qkc  ;
              yae = qk*sinh(qkl)* dxyz[1] + cosh(qkl)*dir[1]/dir[2] ; 
	     }
	     else if (Magnet_type[magnet][side] == 2)
	     {
	      //                                  ----- vertically focussing
	      ye  =     cos(qkl)* dxyz[1] + sin(qkl)*dir[1]/qkc ;
	      yae = -qk*sin(qkl)* dxyz[1] + cos(qkl)*dir[1]/dir[2]  ;
	      xe  =    cosh(qkl)* dxyz[0] + sinh(qkl)*dir[0]/qkc ;
	      xae = qk*sinh(qkl)* dxyz[0] + cosh(qkl)*dir[0]/dir[2] ;
             } 
             else{cout<<"Bad magnet type "<<Magnet_type[magnet][side]
			   <<endl;}
        dxyz[0] = xe;
        dir[0]  = xae;
        dxyz[1] = ye;
        dir[1]  = yae;
        dir[2]  = sqrt(1. - sq(dir[0]) - sq(dir[1]))*aside[side] ;
        dir[0]  = dir[0]*dir[2] ;
        dir[1]  = dir[1]*dir[2] ;

        if(isay>1)debugfile<<magnet<<" QE ";

	    }     //  end of quadrupole treatment 
        if(isay>1)debugfile<<dxyz[0]<<" "<<dxyz[1]<<" "
		    <<0.5* Magnet_length[magnet][side]*aside[side]+zc
		    <<" "<<dir[0]<<" "<<dir[1]<<"  DXB "<<
                    dxyz[0] +Magnet_xb[magnet][side]<< endl;
  
        dxyz[2] = 0.5* Magnet_length[magnet][side]*aside[side] ;   
        arclen  = arclen + 2.*dxyz[2]*(1. + 0.5*(xa1*dir[0] + ya1*dir[1]) 
			 + 0.1667*((xa1-dir[0])*(xa1-dir[0]) +(ya1-dir[1])*(ya1-dir[1])));      

// std::cout<<" B " <<side<<" "<<zp<<" "<<arclen<<std::endl;

//   Are we still in the beampipe? 
        if(Magnet_aperclass[magnet][side] == 1) 
          {radius = Magnet_aperA[magnet][0][side] ;
	   outside = outside || sq(dxyz[0])  + sq(dxyz[1])  > sq(radius);}
        if(Magnet_aperclass[magnet][side] == 2) 
           {
	    outside = outside ||  sq(dxyz[0]/Magnet_aperA[magnet][2][side])  
	                        + sq(dxyz[1]/Magnet_aperA[magnet][3][side])  > 1.;  
            for(int k = 0; k<2; k++)       if(Magnet_aperA[magnet][k][side] > 0.) 
	      outside = outside || fabs(dxyz[k]) > Magnet_aperA[magnet][k][side] ; 
           }
        if(isay > 1 && outside)
          {debugfile<<" - outside - "<<magnet<<endl;
	   if(isay > 2)
	    debugfile<<dxyz[0]<<" "<<Magnet_aperA[magnet][2][side]<<" "<<dxyz[1]<<" "
	    <<Magnet_aperA[magnet][3][side]<<" "<<Magnet_aperA[magnet][0][side]
	    <<" "<<Magnet_aperA[magnet][1][side]<<endl;
	  } 

        if(useaper && outside) {Magfail->Fill(float(magnet+40*side)); break;} 
	// --------  reject particles that are outside the apertures.

        xp = dxyz[0] + xc ;
        yp = dxyz[1] + yc ;
        zp = dxyz[2] + zc ;

        // if(isay>1)debugfile<<"HALL "<<magnet<<" "<<zp<<" "<<xp<<" "<<yp<<endl;

        float M=1000000.;
	if(irtype==1)debugfile<<magnet<<" "<<zp<<" "<<int(xp*M)<<" "<<int(dxyz[0]*M)<<endl;

//  See if we are entering a COLLIMATOR and if necessary remove the particle.


        coutside = false;
        collhit = 99999999.;  // z position of first collimator hit.
        for(int ic=1; ic<=ncoll; ic++)
	  {if(magnet != Coll_nmag[ic][side]) continue;  // ---- in right region?
	   zposition = aside[side]*Coll_z[ic];
           xposition = xp + (zposition - zp)*dir[0]/dir[2]; 
           if(iruncon==1) 
             {
              xcoll0[ic][side]  = xposition;
	      if(isay>0)debugfile<<"INIC "<<ic<<" "<<side<<" "<<zposition<<" "<<xposition<<endl;
	     } 
           else
	     {
	      xwrtbeam  = xposition - xcoll0[ic][side] ;
              if (fabs(xwrtbeam) > Coll_xap[ic])
                {coutside = true;
                if(fabs(zposition) < collhit) collhit = fabs(zposition);}
	     }
	  } 
          if(coltinue==0) outside = outside||coutside ;  //  Normal case, flag outside collimator.
	  if(coutside && isay>1)debugfile<<"- hit collimator - "<<magnet<<" "<<xwrtbeam<<endl;

//  See if we are entering a ROMAN POTS area and if so store info, do plots.
 
        for(int ip=1; ip<=nplane; ip++)
	  {if(magnet != Plane_nmag[ip][side]) continue;  // ---- in right region?
	   zposition = aside[side]*Plane_z[ip];  // z of plane
                            
           if(fabs(zposition) > collhit && iruncon==3 && collhit < 9999.)  
                     continue;   // It hit a collimator first.
           xposition = xp + (zposition - zp)*dir[0]/dir[2]; 
           yposition = yp + (zposition - zp)*dir[1]/dir[2]; 
           xwrtbeam  = xposition + xdis;
           xdisp     = 0. ;

           if(initialising) 
             { xplane0[ip][side]  = xwrtbeam;
               aplane0[ip][side]  = dir[0] ;
	       yplane0[ip][side]  = yposition;
               ayplane0[ip][side] = dir[1] ; 
	       arclen0[ip][side]  = arclen ;
               if(isay>0) debugfile<<"INIP "<<ip<<" "<<side<<" "<<xwrtbeam<<" "
                                  <<dir[0]<<endl;
	       goto ENDPLANES ;
             } // ------------------------end of initialising.

           // Evaluate the displacements wrt calibrated beamline.
           idetec[ip][side]=0;
           xdisp = xwrtbeam  - xplane0[ip][side] ; // correct for initial beam.
           ydisp = yposition - yplane0[ip][side] ;
           xangle = dir[0]   - aplane0[ip][side] ;
           yangle = dir[1]   - ayplane0[ip][side] ;
           if(fabs(zposition)<270.)               // my sign convention in region. 
                               {xdisp  = -xdisp  ;
	                        xangle = -xangle ;}  //  +ve into Si detector.
           if( (iruncon==4 && ip==1) || (iruncon==5 && ip==2))
	      {
	       if(xdisp > xinnerp[ip][side])   iaccep[i_p][ipt][side]++;
               goto ENDPLANES;
              }
     	   // --- Now fill up arrays for the calibration -----------. 
           n=1;                 //  fills a few numbers for plotting.
           s=lengthscale;

           if(iruncon==2)        //  calibration always uses this setting 
              {                  //  Fill for phi = 0. and pi.
	       n = nev; 
               if(n>=50000) goto ENDPLANES; 
	      }
           accepted = true;
           xscan[n][ip][side]   = xdisp ; 
           axscan[n][ip][side]  = xangle ;
           pscan[n][ip][side]   = p ;
           yscan[n][ip][side]   = ydisp ;
           ayscan[n][ip][side]  = yangle ;	
	   Precoval[1][side]    = pbeam0;
           if(!initialising && !calibrating && ip==ireadp && outputdat)    //  Final plane 
	     outputdatafile<<nev<<" "<<p<<" "<<zposition<<" "<<xwrtbeam
                           <<" "<<yposition<<" "<<xangle<<" "<<yangle<<endl;
           //  Calibration shift (silicon has moved)
	   if(!calibrating)
               {xshift = 0.;
		 xdisp = xdisp + xshift ;
                 // xangle = xangle + 2.0e-06;
               }
                //  "If it moves in compared to nominal position, we accept more.
                //   At the same time, the measured value increases compared to physical."  

            if(xdisp > xinnerp[ip][side] - Si_edge) idetec[ip][side]=1; 
                                                 //  PARTICLE HITS DETECTOR MATERIAL       
            if(xdisp > xinnerp[ip][side] ) //  Si width.  PARTICLE ACCEPTED.
	    if(!calibrating)
             {
             //  ---- CALL PFORMULA -----
             Precoval[1][side]=Pformula(IP,ip,side,xdisp,ydisp,xangle,yangle);
             idetec[ip][side]=2 ;
             DPvP[ip][side]->Fill(pbeam0-Precoval[1][side],p0[side]-pmom[side]);
             DPvx[ip][side]->Fill(pbeam0-Precoval[1][side],xdisp*s);
             DPvy[ip][side]->Fill(pbeam0-Precoval[1][side],ydisp*s);
             if(ip<3) Preco[ip][side]->Fill(pbeam0-Precoval[1][side]);
             if(ip<3 && xdisp>0.) DXbyX[ip][side]->Fill((xangle - 0.0167*xdisp)*1000000.);
             //float xxxx = ((xangle - 0.0167*xdisp)*1000000.);
	     DelPx[ip][side]->Fill(pxval-ppx0[side]);       
             DelPy[ip][side]->Fill(pyval-ppy0[side]);       
             DelP[ip][side]->Fill(Precoval[1][side]-pmom[side]);       
             ptval[side]   =  sqrt(sq(pxval)+sq(pyval));
             DelPt[ip][side]->Fill(ptval[side] -ptval0[side]);
             azival[side]  = atan2(pyval,pxval); 
             azival0[side] = atan2(ppy0[side],ppx0[side]); 
	     // arclen     = arclen + (zposition-zp)/dir[2] ; //ignore this bit.
             ArcL[ip][side]->Fill((arclen - arclen0[ip][side])*1000.*aside[side]) ;  

            //  The Precoval value already has effects of smearing included.
             }
	   //  Plane 1 momentum becomes replaced by Plane 2 value later if relevant.

           if(iruncon==3 and coltinue==ip) outside=coutside ;
	                 // -- make the collimator effect active from now on. 

           themr = theta*1000.;
	   if(isay>1 && (iruncon!=3 || nev<100)) printf(
	   "Plane%i pt=%5.3f themr=%6.3f phi=%4.2f x,ydisp=%5.3e %5.3e x,y,zposn=%5.3f %5.3e %5.1f ax=%7.2e p,precon=%7.2f %7.2f %i \n"
	   ,ip,ptval0[side],themr,phi,xdisp,ydisp,xposition,yposition,zposition,xangle,p,Precoval[1][side],i_p);

           chromprin=false;
	   if(chrom && chromprin)
            {ipchoice=2; sidechoice=1; xamur=xangle*1000000.;
              if(ipchoice==ip && sidechoice==side) 
                printf("side %i plane %i %7.1f %8.5f %7.1f %7.1f %7.1f \n" ,
                side,ip,pmom[side],xip,thip, 1000.*xdisp, xamur) ;
            }
	    ENDPLANES:;

       }  //  ---plane loop
	if(useaper && outside) break;  //  Terminate if outside an aperture or collimator.
     }   //  ---magnet loop


    if(accepted && iruncon==2) //  This is for Christophe R.
      {if(fabs(yscan[n][2][side])>xinnerp[2][side])  
             x12[side]->Fill(xscan[n][1][side]*s,xscan[n][2][side]*s);
       if(yscan[n][1][side]*yscan[n][2][side]!=0.)  
          y12[side]->Fill(yscan[n][1][side]*s,yscan[n][2][side]*s);
       xy12[side]->Fill(xscan[n][1][side]*s,yscan[n][2][side]*s);
       // if(yscan[n][1][side]*yscan[n][2][side]!=0.)  std::cout<<side<<" "<<xscan[n][2][side]<<" "<<yscan[n][2][side]<<" "<<zposition<<" "<<ptset<<std::endl;
       
      }
   }  //  --- side loop

    if(accepted && iruncon==3) debugfile<<"Acc "; 
//    if(isay>0 && iruncon==3) debugfile<<" Evt "<<nev<<" P1S0,P2S0,P1S1,P2S1 = "<<
    if(accepted && iruncon==3) debugfile<<" Evt "<<nev<<" P1S0,P2S0,P1S1,P2S1 = "<<
    idetec[1][0]<<idetec[2][0]<< idetec[1][1]<< idetec[2][1]
                <<" "<<Precoval[1][0]<<" "<<Precoval[1][1]<<endl;

    if(!initialising && (!calibrating || iruncon==2) && accepted) 
         Output(iruncon,n) ; 
 }  // event loop 
  nev--;  

  if(iruncon==2 && !initialising && !chrom)
               {
                Calibrate_p(npt, ndelp); 
                Calibrate_px(npt, ndelp); 
                Calibrate_py(npt, ndelp); 
               }
                //   always calibrate with this iruncon setting. 
                //   fits to evaluate the momentum formula parameters.
                //   first for each p value do an x and a theta_0 fit 

  if(iruncon==3) Statistics();
 
  if(iruncon==4 || iruncon==5) 
    {for(int isi=0; isi<2; isi++)
     {debugfile<<"IP="<<IP<<" xinner(m)="<<xinnerp[iruncon-3][isi]
                  <<" Side="<<isi<<" Plane="<<iruncon-3<<endl;
         for(imom=0; imom<200; imom++)
          {for(ipt=0; ipt<10; ipt++)debugfile<<iaccep[imom][ipt][isi]<<" ";
	    debugfile<<endl;
    } } } 
  MyInput.close() ;   //  Close the input file. 
  if(chrom) cfile.close();
}
// ==============================================================
void FPTrack::Calibrate_p(int npt, int ndelp)
{
 //  Do fits to establish calibration constants.
    double c1calc=0., c2calc=0., pprev, pnew=0., pval, pvals[100], 
    f0, f1, f2, alpha=0,beta=0,gamma=0, aa=0,bb=0,cc=0;
   float  c1_0 =0., c1_1=0., c1_2=0., c2_1=0.;
   int    np, j, nptx,i_p, n0,nnn ;   
   float a, p, x;

   nptx = npt   ;
   n0   = npt/2 ;
   if(isay>0) cout<<"side,plane alpha,beta, a,b, c_0,c_1 "<<endl; 

   for (int side=0; side<2; side++) {
   for (int ip=1; ip<=nplane ; ip++) {

   if(isay>0) debugfile<<" side, ipla "<<side<<" "<<ip<<endl;
     
     // select planes according to icali    
     if(ip==1 && icali==2) continue;
     if(ip==2 && icali==1) continue;

     float sx=0.,    sx2=0., sx3=0., sx4=0.,    
	   sdpx=0., sdpx2=0., disp, disa, D ;
     double sf01=0.,sf02=0.,sf12=0., sf00=0., sf11=0., sf22=0.,saf0=0.,
            saf1=0., saf2=0.,spf2=0.,spf1=0.,spf0=0.;

     //  The first fit looks at beamline particles at series of p values
     // 
     for(i_p=0; i_p<ndelp; i_p++) pvals[i_p]=0.;
     for (j=0; j <ndelp*nptx*nptx; j++)  
        {
        if(pscan[j][ip][side]==0.) continue; //  wasn't filled    
	 nnn  =  j/nptx;     //  = npty 
         i_p   =  nnn/nptx; 
         if(  j - nptx*nnn != n0) continue;    //  only use one pty scan
         if(nnn - i_p*nptx != n0) continue;    //  only use ptx=0
 
          pvals[i_p]  = p0[side] - pscan[j][ip][side];
          a = axscan[j][ip][side] ;
          p = pvals[i_p] ;
          x = xscan[j][ip][side] ;
          f0   = x ; 
          f1   = x*x ;
          f2   = x*x*x ;
          sf01  = sf01  + f0*f1 ;
          sf02  = sf02  + f0*f2; 
          sf12  = sf12  + f1*f2; 
          sf00  = sf00  + sq(f0);      
          sf11  = sf11  + sq(f1); 
          sf22  = sf22  + sq(f2); 
          saf0  = saf0   + a*f0 ; 
          saf1  = saf1   + a*f1;  
          saf2  = saf2   + a*f2;  
          spf0  = spf0   + p*f0; 
          spf1  = spf1   + p*f1;  
          spf2  = spf2   + p*f2;  
        }
        //  To deal with singular cases, reduce matrix:

       if(ip==2){sf22=0.;sf02=0.;sf12=0.;saf2=0.; spf2=0.;} //turn off cubic

        if(sf22==0.) sf22=1.; if(sf11==0.) sf11=1.; if(sf00==0.) sf00=1.;

	  D = sf00*sf11*sf22 -sf00*sf12*sf12 -sf22*sf01*sf01 -sf11*sf02*sf02 
               +2.*sf01*sf02*sf12;
           if(D !=0)
            {
	    alpha = saf0*(sf11*sf22-sf12*sf12) -saf1*(sf01*sf22-sf02*sf12) 
                   +saf2*(sf01*sf12 -sf02*sf11);
            alpha = alpha/D ; 
	    beta  = saf1*(sf22*sf00-sf02*sf02) -saf2*(sf12*sf00-sf01*sf02) 
                   +saf0*(sf12*sf02 -sf01*sf22);
            beta  = beta/D ;
 	    gamma = saf2*(sf00*sf11-sf01*sf01) -saf0*(sf02*sf11-sf12*sf01) 
                   +saf1*(sf02*sf01 -sf12*sf00);
            gamma = gamma/D ;
	    aa    = spf0*(sf11*sf22-sf12*sf12) -spf1*(sf01*sf22-sf02*sf12) 
                   +spf2*(sf01*sf12 -sf02*sf11);
            aa    = aa/D ;
	    bb    = spf1*(sf22*sf00-sf02*sf02) -spf2*(sf12*sf00-sf01*sf02) 
                   +spf0*(sf12*sf02 -sf01*sf22);
            bb    = bb/D ;
 	    cc    = spf2*(sf00*sf11-sf01*sf01) -spf0*(sf02*sf11-sf12*sf01) 
                   +spf1*(sf02*sf01 -sf12*sf00);
            cc    = cc/D ;
            }
	   else cout <<"Calib-p Bad x  determinant  "<<ndelp<<" "<<nptx<<" "<<endl;
 
    //  Now examine the coefficients with angle.  Do this for each momentum.
     double c[3][50];

     for (np = 0; np < ndelp; np++)
      { 
      sx=0.; sx2=0.; sx3=0.; sx4=0.; sdpx = 0.; sdpx2 = 0.;
      pval = pvals[np] ;
      if(pval!=0.) for (j=np*nptx*nptx + n0; j <(np+1)*nptx*nptx + n0; j=j+nptx)
       {  
        disp  =  aa*xscan[j][ip][side] + bb*sq(xscan[j][ip][side]) 
	       + cc*sq(xscan[j][ip][side])*xscan[j][ip][side] - pval; 
        disa  =  alpha*xscan[j][ip][side] + beta*sq(xscan[j][ip][side])
        	+gamma*sq(xscan[j][ip][side])*xscan[j][ip][side] ;
        disa  = axscan[j][ip][side] - disa;
        sdpx  = sdpx  + disp*disa ; 
        sdpx2 = sdpx2 + disp*sq(disa) ; 
        sx    = sx    + disa ;
        sx2   = sx2   + sq(disa) ;
        sx3   = sx3   + sq(disa)*disa ;
        sx4   = sx4   + sq(sq(disa)) ;

        if(isay>0 && j == np*nptx*nptx + npt*nptx +n0) debugfile<< "centre ";
        if(isay>0) debugfile<<"np/p/a/disp/disa "<<np<<" "<<pval<<" "<<axscan[j][ip][side]<<" "<<disp<<" "<<disa<<" "<<j<<endl;
      } 

      //  The linear coefficient is made to vary as a function of momentum
      //  with two or three  coefficients.  
      //  For ip = 1  f1 = pval dominates
      //  For ip = 2, it is sufficient to have f0 = 1. and f1 = 1./pval
      //  First obtain the coefficients at each momentum. 
      c[1][np]=0.;  
      c[2][np]=0. ;
      if(sx2==0.) continue;
      D = sx2*sx4 - sq(sx3) ;
      if(D != 0.) {c[1][np] =  (sdpx*sx4 - sdpx2*sx3)/D ;
	           c[2][np] = -(sdpx*sx3 - sdpx2*sx2)/D ;}
      else {cout<<"Calib-p Bad determinant for calculating c[]["<<np<<"] ip="<<ip<<endl;
            c[1][np] =  sdpx/sx2 ;}
      }  // --------------------------- end of np loop

     //  Now do the summations to fit the coefficients 
     c1_0=0.; c1_1=0.; c1_2 = 0.; c2_1=0.;
     for(int ic=1; ic<=2; ic++)
       {
        double sf01=0.,sf02=0.,sf12=0., sf00=0., sf11=0., sf22=0.,scf0=0.,
               scf1=0., scf2=0.,scp2=0.,sp=0.,sp2=0.,sp3=0.,sp4=0.,scp=0.;
        for (np = 0; np < ndelp; np++)
	  {
          pval = pvals[np]  ;
          if(pval == 0.) continue; 
          if(c[ic][np]==0.) continue;  // wasn't filled
          f0   = fitf(ip,ic,0,pval) ; 
          f1   = fitf(ip,ic,1,pval) ;
          f2   = fitf(ip,ic,2,pval) ;
          sf01  = sf01  + f0*f1 ;
          sf02  = sf02  + f0*f2; 
          sf12  = sf12  + f1*f2; 
          sf00  = sf00  + sq(f0);      
          sf11  = sf11  + sq(f1); 
          sf22  = sf22  + sq(f2); 
          scf0  = scf0   + c[ic][np]*f0 ; 
          scf1  = scf1   + c[ic][np]*f1;  
          scf2  = scf2   + c[ic][np]*f2;  
          // 
          if(ip!=5) continue;
          scp  = scp   + pval*c[ic][np] ; 
          scp2 = scp2  + c[ic][np]*sq(pval) ; 
          sp   = sp    + pval ;
          sp2  = sp2   + sq(pval) ;
          sp3  = sp3   + sq(pval)*pval ;
          sp4  = sp4   + sq(sq(pval)) ;
         }                             

        //  To deal with some singular cases, reduce matrix:
        if(sf22==0.) sf22=1.;   if(sf11==0.) sf11=1.;   if(sf00==0.) sf00=1.;
        if(ic==1 && (ip==2 || ip==1)) 
	  {D = sf00*sf11*sf22 -sf00*sf12*sf12 -sf22*sf01*sf01 -sf11*sf02*sf02 
               +2.*sf01*sf02*sf12;
           if(D !=0)
            {
	    c1_0 = scf0*(sf11*sf22-sf12*sf12) -scf1*(sf01*sf22-sf02*sf12) 
                  +scf2*(sf01*sf12 -sf02*sf11);
            c1_0 = c1_0/D ; 
	    c1_1 = scf1*(sf22*sf00-sf02*sf02) -scf2*(sf12*sf00-sf01*sf02) 
                  +scf0*(sf12*sf02 -sf01*sf22);
            c1_1 = c1_1/D ;
 	    c1_2 = scf2*(sf00*sf11-sf01*sf01) -scf0*(sf02*sf11-sf12*sf01) 
                  +scf1*(sf02*sf01 -sf12*sf00);
            c1_2 = c1_2/D ;
            }
	    else cout <<"Calib-p Bad final c1 determinant  "<<endl;
          }  // --- end of ic==1 ip==1||2 condition -------
 	  
          if(ic==1 && ip==5)
             {D = sp2*sp4 - sq(sp3);
              if(D==0.) {cout<<" Bad Det!! side/ip"<<side<<" "<<ip<<endl;}
              else {
                    c1_0 =   (scp*sp4 - scp2*sp3)/D ;
                    c1_1 =  -(scp*sp3 - scp2*sp2)/D ;
               }}
          if(ic==2 && (ip==2 || ip==1))
           { if(sf11 !=0.) c2_1 = scf1/sf11;	    
           }    
      }

     for (np = 0; np < ndelp; np++)                 
       {  
       pval=pvals[np] ; 
       if(pval==0) continue;
       c2calc = c2_1*fitf(ip,2,1,pval) ;
       c1calc = c1_0*fitf(ip,1,0,pval) + c1_1*fitf(ip,1,1,pval) +c1_2*fitf(ip,1,2,pval);
       if(c[1][np]!=0.  && isay>0)
       cout<<"ip np c1:meas/fit "<<ip<<" "<<np<<" "<<c[1][np]<<" "<<c1calc<<" "
                                                   <<c[2][np]<<" "<<c2calc<<endl;
       }
     cout<<"coeffs "<<side<<" "<<ip<<" "<<alpha<<" "<<beta<<" "<<      
       aa<<" "<<bb<<" "<<cc<<" "<<c1_0<<" "<<c1_1<<" "<<c1_2<<" "<<c2_1<<endl; 

      j = IP-1 + ip-1 + 2*side;
      cal_a[j]     = aa;     cal_b[j]    = bb;    cal_c[j]    = cc; 
      cal_alpha[j] = alpha;  cal_beta[j] = beta;  cal_gamma[j] = gamma;
      cal_c_0[j]   = c1_0;   cal_c_1[j]  = c1_1;  cal_c_2[j]  = c1_2; 
      cal_c2_1[j]  = c2_1;

      if(isay>0)      // -- work out momentum to print it out --
      for (int j=0; j <ndelp*nptx*nptx; j++)
        {
        if(pscan[j][ip][side]==0.) goto LOOP3 ;  //  wasn't filled    
        np  = j/(nptx*nptx) ; 

        pnew  = aa*xscan[j][ip][side] + bb*sq(xscan[j][ip][side]) 
                +cc*xscan[j][ip][side]*sq(xscan[j][ip][side]) ; 
        disa  = alpha*xscan[j][ip][side] + beta*sq(xscan[j][ip][side])
        	+gamma*sq(xscan[j][ip][side])*xscan[j][ip][side] ;
        disa  = axscan[j][ip][side] - disa;

        pprev=pnew ;
        for(int iter=0; iter< 4  ; iter++){
	  c1calc = c1_0*fitf(ip,1,0,pnew) +c1_1*fitf(ip,1,1,pnew) +c1_2*fitf(ip,1,2,pnew);
	  c2calc = c2_1*fitf(ip,2,1,pnew) ;
           pnew =  pprev - c1calc*disa -c2calc*sq(disa);
          }
        if(isay>1)
        debugfile<<" FitPl "<<ip<<" side"<<side<<" "<<pvals[np]<<" "<<j<<" "<<pnew<<
	  " "<<c1_0<<" "<<c1_1<<" "<<c1calc<<
          " x,a "<<xscan[j][ip][side]<<" "<<axscan[j][ip][side]<<endl;
        LOOP3:;
        }
     } } 
     cout<<"Calibrations-p set "<<endl ;
     cout<<"alpha ";  
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_alpha[j]<<", ";}
     cout<<endl<<"beta ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_beta[j]<<", ";}
     cout<<endl<<"gamma ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_gamma[j]<<", ";}
     cout<<endl<<"a ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_a[j]<<", ";}
     cout<<endl<<"b ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_b[j]<<", ";}
     cout<<endl<<"c ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_c[j]<<", ";}
     cout<<endl<<"c_0 ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_c_0[j]<<", ";}
     cout<<endl<<"c_1 ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_c_1[j]<<", ";}
     cout<<endl<<"c_2 ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_c_2[j]<<", ";}
     cout<<endl<<"c_21 ";
     for(int j = IP-1; j<IP+3; j++) {cout<<cal_c2_1[j]<<", ";}
     cout<<endl;   
}
// ==============================================================
void FPTrack::Calibrate_px(int npt, int ndelp)
{
 //  Do fits to establish calibration constants.
    double pxnew=0., pxvals[100], 
      f0, f1, f2, alpha=0,beta=0,gamma=0, x0vals[100], a0vals[100],
      a, p, x, ccalc[3];
   float  cfit[5][3][3];
   int    np, j, nptx,i_p, n0,nnn ;   

   nptx = npt   ;
   n0   = npt/2 ;
   if(isay>0) cout<<"side,plane alpha,beta, a,b, c_0,c_1 "<<endl; 
   for(i_p=0; i_p<nptx; i_p++) pxvals[i_p]= (float(nptx)*0.5  -0.5 -float(i_p))*dptval;

   for (int side=0; side<2; side++) {
   for (int ip=1; ip<=nplane ; ip++) {
   for(np=0; np<ndelp; np++) x0vals[np]= 0.;
   for(np=0; np<ndelp; np++) a0vals[np]= -999.;


   if(isay>0) debugfile<<" side, ipla "<<side<<" "<<ip<<endl;
     
     // select planes according to icali    
     if(ip==1 && icali==2) continue;
     if(ip==2 && icali==1) continue;

     float sx=0.,    sx2=0., sx3=0., sx4=0.,    
       sdpx=0., sdpx2=0., D, dispx, disx ;
     double sf01=0.,sf02=0.,sf12=0., sf00=0., sf11=0., sf22=0.,saf0=0.,
            saf1=0., saf2=0.,spf2=0.,spf1=0.,spf0=0.;

     //  The first fit looks at beamline particles at series of p values
     //  and asks what is the expected x value at a given angle.

     for (j=0; j <ndelp*nptx*nptx; j++)  
        {
        if(pscan[j][ip][side]==0.) continue; //  wasn't filled    
	 nnn  =  j/nptx;     //  = npty 
         np   =  nnn/nptx; 
         if(  j - nptx*nnn != n0) continue;    //  only use one pty scan
         if(nnn - np*nptx != n0) continue;    //  only use ptx=0

       	  a = axscan[j][ip][side]*1000. ;   //  work in mrad
          if(a==0.) continue;
          p = pxvals[np] ;
          x = xscan[j][ip][side] ;
          x0vals[np] = x;
          a0vals[np] = a;
          f0   = a ; 
          f1   = a*a;
          f2   = a*a*a ;
          sf01  = sf01  + f0*f1 ;
          sf02  = sf02  + f0*f2; 
          sf12  = sf12  + f1*f2; 
          sf00  = sf00  + sq(f0);      
          sf11  = sf11  + sq(f1); 
          sf22  = sf22  + sq(f2); 
          saf0  = saf0   + x*f0 ; 
          saf1  = saf1   + x*f1;  
          saf2  = saf2   + x*f2;  
          spf0  = spf0   + p*f0; 
          spf1  = spf1   + p*f1;  
          spf2  = spf2   + p*f2;  
        }
        //  To deal with singular cases, reduce matrix:

        if(ip==2){sf22=0.;sf02=0.;sf12=0.;saf2=0.; spf2=0.;} //turn off cubic

        if(sf22==0.) sf22=1.; if(sf11==0.) sf11=1.; if(sf00==0.) sf00=1.;

	D = sf00*sf11*sf22 -sf00*sf12*sf12 -sf22*sf01*sf01 -sf11*sf02*sf02 
               +2.*sf01*sf02*sf12;
           if(D !=0)
            {
	    alpha = saf0*(sf11*sf22-sf12*sf12) -saf1*(sf01*sf22-sf02*sf12) 
                   +saf2*(sf01*sf12 -sf02*sf11);
            alpha = alpha/D ; 
	    beta  = saf1*(sf22*sf00-sf02*sf02) -saf2*(sf12*sf00-sf01*sf02) 
                   +saf0*(sf12*sf02 -sf01*sf22);
            beta  = beta/D ;
 	    gamma = saf2*(sf00*sf11-sf01*sf01) -saf0*(sf02*sf11-sf12*sf01) 
                   +saf1*(sf02*sf01 -sf12*sf00);
            gamma = gamma/D ;


            }
	    else cout <<"Calib-px Bad x  determinant  "<<endl;
 
    //  Now examine the coefficients with position x. 
    //  Do this for each momentum value

     double c[3][50];

     for (np=0 ; np< ndelp; np++)
       {c[1][np]=0.;  
        c[2][np]=0.;
        int j0 = np*nptx*nptx + n0*nptx + n0 ;
	if(a0vals[np] < -998.) continue; // wasn't filled
        if(pscan[j0][ip][side]==0.) continue; //  wasn't filled           
        sx=0.; sx2=0.; sx3=0.; sx4=0.; sdpx = 0.; sdpx2 = 0.;
        int npoints = 0;
        for (i_p = 0; i_p < nptx; i_p++)
        { 
        j =  np*nptx*nptx + n0 + i_p*nptx;  
        if(pscan[j][ip][side]==0.) continue; //  wasn't filled    
        a = axscan[j][ip][side]*1000. ;   //  work in mrad
	dispx =  pxvals[i_p];
        disx  =  alpha*a + beta*sq(a) + gamma*sq(a)*a ;
        disx  = xscan[j][ip][side] - disx;
        sdpx  = sdpx  + dispx*disx ; 
        sdpx2 = sdpx2 + dispx*sq(disx) ; 
        sx    = sx    + disx ;
        sx2   = sx2   + sq(disx) ;
        sx3   = sx3   + sq(disx)*disx ;
        sx4   = sx4   + sq(sq(disx)) ;
        npoints++;
        if(isay>0) debugfile<<"i_p/dispx/x/a/disx "<<i_p<<" "<<dispx<<" "
		   <<xscan[j][ip][side]<<" "<<a<<" "<<disx<<" "<<j<<endl;
      } // --------------------------- end of i_p loop 
      //  The coefficients are made to vary as a function of p 
      //  with two or three  coefficients.  
      //  First obtain the coefficients at each p. 

      if(npoints<5) continue;
      D = sx2*sx4 - sq(sx3) ;

      // D  = 0.; // ========= TEST
      if(D != 0.) {c[1][np] =  (sdpx*sx4 - sdpx2*sx3)/D ;
	           c[2][np] = -(sdpx*sx3 - sdpx2*sx2)/D ;}
      else {cout<<"Calib-px Bad determinant for c[]["<<np<<"] ip="<<ip<<endl;
            c[1][np] =  sdpx/sx2 ;} 


      }  // --------------------------- end of np loop

     //  Now do the summations to fit the coefficients as a function of p 
     //  The planes must be treated differently.

   for(int ix=0; ix<2; ix++)  for(int ia=0; ia<2; ia++) cfit[ip][ix][ia]=0.;  
   for(int ix=1; ix<=2; ix++)
       {
        double sf01=0.,sf02=0.,sf12=0.,sf00=0.,sf11=0.,sf22=0.,scf0=0.,scf1=0.,scf2=0.;
       for (np = 0; np < ndelp; np++)
	  {
	  if(a0vals[np] < -998.) continue; // wasn't filled
          a = pgrid(np,delp);
          f0   = fitpxf(ip,side,0,a) ; 
          f1   = fitpxf(ip,side,1,a) ; 
          f2   = fitpxf(ip,side,2,a) ; 
          sf01  = sf01  + f0*f1 ;
          sf02  = sf02  + f0*f2; 
          sf12  = sf12  + f1*f2; 
          sf00  = sf00  + sq(f0);      
          sf11  = sf11  + sq(f1); 
          sf22  = sf22  + sq(f2); 
          scf0  = scf0   + c[ix][np]*f0 ; 
          scf1  = scf1   + c[ix][np]*f1;  
          scf2  = scf2   + c[ix][np]*f2;  
          }
        //  To deal with some singular cases, reduce matrix:
        if(sf22==0.) sf22=1.;   if(sf11==0.) sf11=1.;   if(sf00==0.) sf00=1.;
        if(ip==2 || ip==1) 
	  {D = sf00*sf11*sf22 -sf00*sf12*sf12 -sf22*sf01*sf01 -sf11*sf02*sf02 
               +2.*sf01*sf02*sf12;
           if(D !=0)
            {
	    cfit[ip][ix][0] = ( scf0*(sf11*sf22-sf12*sf12) -scf1*(sf01*sf22-sf02*sf12) 
                  +scf2*(sf01*sf12 -sf02*sf11) )/D;
	    cfit[ip][ix][1] = ( scf1*(sf22*sf00-sf02*sf02) -scf2*(sf12*sf00-sf01*sf02) 
                  +scf0*(sf12*sf02 -sf01*sf22) )/D;
 	    cfit[ip][ix][2] = ( scf2*(sf00*sf11-sf01*sf01) -scf0*(sf02*sf11-sf12*sf01) 
                  +scf1*(sf02*sf01 -sf12*sf00) )/D;
            }
	    else cout <<"Calib-px Bad final c1 determinant  "<<endl;
          }  // --- ip==1||2 condition -------

       }  // -- end of ix loop

     // for(i_p=0; i_p<nptx; i_p++) debugfile<<pxvals[i_p]<<" "<<endl;

   if(isay>0) for (np = 0; np < ndelp; np++)       // -- Check the fit.             
       {  
       if(a0vals[np] < -998.) continue; // wasn't filled
           a = pgrid(np,delp);
	   f0   = fitpxf(ip,side,0,a) ; 
           f1   = fitpxf(ip,side,1,a) ; 
           f2   = fitpxf(ip,side,2,a) ; 
       ccalc[1] = cfit[ip][1][0]*f0  + cfit[ip][1][1]*f1 + cfit[ip][1][2]*f2;
       ccalc[2] = cfit[ip][2][0]*f0  + cfit[ip][2][1]*f1 + cfit[ip][2][2]*f2;	
       if(c[1][np]!=0.)
       cout<<"ip np a/c1,2:meas/fit "<<ip<<" "<<np<<" "<<a<<" "<<c[1][np]
        <<"/"<<ccalc[1]<<" "<<c[2][np]<<"/"<<ccalc[2]<<endl;
       }

      j = IP-1 + ip-1 + 2*side;
      calpx_alpha[j] = alpha;  calpx_beta[j] = beta;  calpx_gamma[j] = gamma;
      calpx_cfit[j][1][0] = cfit[ip][1][0]; calpx_cfit[j][2][0] = cfit[ip][2][0]; 
      calpx_cfit[j][1][1] = cfit[ip][1][1]; calpx_cfit[j][2][1] = cfit[ip][2][1]; 
      calpx_cfit[j][1][2] = cfit[ip][1][2]; calpx_cfit[j][2][2] = cfit[ip][2][2]; 

      if(isay>0)      // -- work out momentum to print it out --
      for (np=0; np<ndelp; np++)
      for(i_p=0; i_p<nptx; i_p++)
        {
        j =  np*nptx*nptx + n0 + i_p*nptx;  
        if(pscan[j][ip][side]==0.) goto LOOP3 ;  //  wasn't filled    
        a = axscan[j][ip][side]*1000. ;
        if(a<-998.) goto LOOP3;
        disx  = xscan[j][ip][side] - (alpha*a + beta*sq(a) +gamma*sq(a)*a);
           a = pgrid(np,delp);
	   f0   = fitpxf(ip,side,0,a) ; 
           f1   = fitpxf(ip,side,1,a) ; 
           f2   = fitpxf(ip,side,2,a) ; 
       ccalc[1] = cfit[ip][1][0]*f0  + cfit[ip][1][1]*f1 + cfit[ip][1][2]*f2;
       ccalc[2] = cfit[ip][2][0]*f0  + cfit[ip][2][1]*f1 + cfit[ip][2][2]*f2;
        pxnew =  ccalc[1]*disx +ccalc[2]*sq(disx);

        if(isay>0)
        debugfile<<" FitPlax"<<ip<<" side"<<side<<" "<<pxvals[i_p]<<" "<<np<<" "<<i_p
        <<" "<<pxnew<< " "<<cfit[ip][1][0]<<" "<<cfit[ip][1][1]<<" "<<ccalc[1]<<
          " x,a "<<xscan[j][ip][side]<<" "<<axscan[j][ip][side]<<endl;

       LOOP3:;
        }
     } } 
     
}
// ==============================================================
void FPTrack::Calibrate_py(int npt, int ndelp)
{
 //  Do fits to establish calibration constants.
    double pynew=0., pyvals[100], 
      f0, f1, f2, alpha=0,beta=0,gamma=0, y0vals[100], a0vals[100], aa=0.,bb=0.,cc=0.,
      a, p=0., y, cycalc[3], caycalc[3], wy,way,py_y,py_ay, pvals[100],
      cyfit[5][3][3],cayfit[5][3][3],cy[3][50],cay[3][50];
   int    np, j, npty,i_p, n0,nnn ;   

   npty = npt   ;
   n0   = npt/2 ;
   if(isay>0) cout<<"y side,plane alpha,beta, a,b, c_0,c_1 "<<endl; 
   for(i_p=0; i_p<npty; i_p++) pyvals[i_p]= (float(npty)*0.5  -0.5 -float(i_p))*dptval;

   for (int side=0; side<2; side++) {
   for (int ip=1; ip<=nplane ; ip++) {
   for(np=0; np<ndelp; np++) y0vals[np]= 0.;
   for(np=0; np<ndelp; np++) a0vals[np]= -999.;


   if(isay>0) debugfile<<"y side, ipla "<<side<<" "<<ip<<endl;
     
     // select planes according to icali    
     if(ip==1 && icali==2) continue;
     if(ip==2 && icali==1) continue;

     double sy=0.,    sy2=0., sy3=0., sy4=0.,sdpy=0., sdpy2=0., sdpa2=0.,
            disa, D, dispy, disy,sa=0.,sa2=0., sa3=0., sa4=0.,sdpa=0.,
            sf01=0.,sf02=0.,sf12=0., sf00=0., sf11=0., sf22=0.,saf0=0.,
            saf1=0., saf2=0.,syf2=0.,syf1=0.,syf0=0.;

     //  The first fit looks at beamline particles at series of p values
     //  and asks what is the expected y value and angle value at given p

     for (j=0; j <ndelp*npty*npty; j++)  
        {
        if(pscan[j][ip][side]==0.) continue; //  wasn't filled    
	 nnn  =  j/npty;     //  = nptx
         np   =  nnn/npty; 
         if( j - npty*nnn != n0) continue;    //  only use pty=0
         if(nnn - np*npty != n0) continue;    //  only use ptx=0
       	  a = ayscan[j][ip][side]*1000. ;   //  work in mrad
          p = pgrid(np,delp);
          y = yscan[j][ip][side] ;
          if(pscan[j][ip][side]==0.) continue; //  wasn't filled.
          y0vals[np] = y;
          a0vals[np] = a;
          pvals[np]  = p;
          f0   = p ; 
          f1   = p*p;
          f2   = p*p*p ;
          sf01  = sf01  + f0*f1 ;
          sf02  = sf02  + f0*f2; 
          sf12  = sf12  + f1*f2; 
          sf00  = sf00  + sq(f0);      
          sf11  = sf11  + sq(f1); 
          sf22  = sf22  + sq(f2); 
          saf0  = saf0   + a*f0 ; 
          saf1  = saf1   + a*f1;  
          saf2  = saf2   + a*f2;  
          syf0  = syf0   + y*f0; 
          syf1  = syf1   + y*f1;  
          syf2  = syf2   + y*f2;  
        }
        //  To deal with singular cases, reduce matrix:

        if(ip==2){sf22=0.;sf02=0.;sf12=0.;saf2=0.; syf2=0.;} //turn off cubic

        if(sf22==0.) sf22=1.; if(sf11==0.) sf11=1.; if(sf00==0.) sf00=1.;

	D = sf00*sf11*sf22 -sf00*sf12*sf12 -sf22*sf01*sf01 -sf11*sf02*sf02 
               +2.*sf01*sf02*sf12;
           if(D !=0)
            {
	    alpha = saf0*(sf11*sf22-sf12*sf12) -saf1*(sf01*sf22-sf02*sf12) 
                   +saf2*(sf01*sf12 -sf02*sf11);
            alpha = alpha/D ; 
	    beta  = saf1*(sf22*sf00-sf02*sf02) -saf2*(sf12*sf00-sf01*sf02) 
                   +saf0*(sf12*sf02 -sf01*sf22);
            beta  = beta/D ;
 	    gamma = saf2*(sf00*sf11-sf01*sf01) -saf0*(sf02*sf11-sf12*sf01) 
                   +saf1*(sf02*sf01 -sf12*sf00);
            gamma = gamma/D ;

	    aa    = syf0*(sf11*sf22-sf12*sf12) -syf1*(sf01*sf22-sf02*sf12) 
                   +syf2*(sf01*sf12 -sf02*sf11);
            aa    = aa/D ;
	    bb    = syf1*(sf22*sf00-sf02*sf02) -syf2*(sf12*sf00-sf01*sf02) 
                   +syf0*(sf12*sf02 -sf01*sf22);
            bb    = bb/D ;
 	    cc    = syf2*(sf00*sf11-sf01*sf01) -syf0*(sf02*sf11-sf12*sf01) 
                   +syf1*(sf02*sf01 -sf12*sf00);
            cc    = cc/D ;
            }
	    else cout <<"Calib-py Bad y  determinant  "<<endl;

     //  Now examine the coefficients with position y and angle ay
    //  Do this for each momentum value

     double c[3][50];

     for (np=0 ; np< ndelp; np++)
       {cy[1][np]=0.;  cy[2][np]=0.; cay[1][np]=0.;  cay[2][np]=0.;
        int j0 = np*npty*npty + n0*npty + n0 ;
        if(pscan[j0][ip][side]==0.) continue; //  wasn't filled           
        sy=0.; sy2=0.; sy3=0.; sy4=0.; sdpy = 0.; sdpy2 = 0.;
        sa=0.; sa2=0.; sa3=0.; sa4=0.; sdpa = 0.; sdpa2 = 0.;
        int npoints = 0;
        p = pgrid(np,delp);
        for (i_p = 0; i_p < npty; i_p++)
        { 
        j =  np*npty*npty + n0*npty + i_p;  
        if(pscan[j][ip][side]==0.) continue; //  wasn't filled    
        a = ayscan[j][ip][side]*1000. ;   //  work in mrad
	dispy =  pyvals[i_p];
        disa  =  a -(alpha*p + beta*sq(p) + gamma*sq(p)*p) ;
        disy  =   yscan[j][ip][side] -(aa*p + bb*sq(p) + cc*sq(p)*p) ;
        sdpa  = sdpa  + dispy*disa ; 
        sdpa2 = sdpa2 + dispy*sq(disa) ; 
        sdpy  = sdpy  + dispy*disy ; 
        sdpy2 = sdpy2 + dispy*sq(disy) ; 
        sa    = sa    + disa ;
        sa2   = sa2   + sq(disa) ;
        sa3   = sa3   + sq(disa)*disa ;
        sa4   = sa4   + sq(sq(disa)) ;     
        sy    = sy    + disy ;
        sy2   = sy2   + sq(disy) ;
        sy3   = sy3   + sq(disy)*disy ;
        sy4   = sy4   + sq(sq(disy)) ;
        npoints++;
        if(isay>0) debugfile<<"i_p/dispy/y/a/disy/disa "<<i_p<<" "<<dispy<<" "
		   <<yscan[j][ip][side]<<" "<<a<<" "<<disy<<" "<<disa<<" "<<j<<endl;
      } // --------------------------- end of i_p loop 
      //  The linear coefficient is made to vary as a function of disy 
      //  with two or three  coefficients.  
      //  First obtain the coefficients at each p. 
      if(npoints<5) continue; 
      D = sa2*sa4 - sq(sa3) ;
      if(D != 0.) {cay[1][np] =  (sdpa*sa4 - sdpa2*sa3)/D ;
	           cay[2][np] = -(sdpa*sa3 - sdpa2*sa2)/D ;}
      else {cout<<"Bad determinant for ca[]["<<np<<"] ip="<<ip<<endl;
            cay[1][np] =  sdpa/sa2 ;} 

      D = sy2*sy4 - sq(sy3) ;
      if(D != 0.) {cy[1][np] =  (sdpy*sy4 - sdpy2*sy3)/D ;
	           cy[2][np] = -(sdpy*sy3 - sdpy2*sy2)/D ;}
      else {cout<<"Calib-py Bad determinant for cay[]["<<np<<"] ip="<<ip<<endl;
            cy[1][np] =  sdpy/sy2 ;} 

      }  // --------------------------- end of np loop


/*
     for (np=0 ; np< ndelp; np++)
       {
        int j0 = np*npty*npty + n0*npty + n0 ;
        if(pscan[j0][ip][side]==0.) continue; //  wasn't filled           
                p = pgrid(np,delp);
        for (i_p = 0; i_p < npty; i_p++)
        { 
        j =  np*npty*npty + n0*npty + i_p;  
        if(pscan[j][ip][side]==0.) continue; //  wasn't filled    
        a = ayscan[j][ip][side]*1000. ;   //  work in mrad
	dispy =  pyvals[i_p];
        disa  =  a -(alpha*p + beta*sq(p) + gamma*sq(p)*p) ;
        disy  =   yscan[j][ip][side] -(aa*p + bb*sq(p) + cc*sq(p)*p) ;
        float pay = cay[1][np]*disa + cay[2][np]*disa*disa ;
        float py  = cy[1][np]*disy + cy[2][np]*disy*disy ;
 
       if(isay>0) debugfile<<"np/i_p/dispy/disy/disa/cy[1]/cay[1]/py/pya "
       <<np<<" "<<i_p<<" "<<dispy<<" "<<disy<<" "<<disa<<" "<<cy[1][np]<<" "<<cay[1][np]
	<<" "<<py<<" "<<pay<<" "<<j<<endl;
      } // --------------------------- end of i_p loop 
      }  // --------------------------- end of np loop
*/


     //  Now do the summations to fit the coefficients as a function of momentum
     //  The planes must be treated differently.
     //  There is a singularity in the y coefficient, so use its reciprocal.

   for(int ic=0; ic<2; ic++)  for(int ia=0; ia<2; ia++) 
         {cyfit[ip][ic][ia]=0.;cayfit[ip][ic][ia]=0.;}   
   for(int ic=1; ic<=2; ic++)
       {
        double sf01=0.,sf02=0.,sf12=0.,sf00=0.,sf11=0.,sf22=0.,
               scayf1=0.,scayf2=0.,scayf0=0.,scyf0=0.,scyf1=0.,scyf2=0.;
       for (np = 0; np < ndelp; np++)
	  {
	  if(a0vals[np] < -998.) continue; // wasn't filled
          p = pvals[np];
          cy[ic][np] = 1./ cy[ic][np];
          f0   = fitpyf(ip,side,0,p) ; 
          f1   = fitpyf(ip,side,1,p) ; 
          f2   = fitpyf(ip,side,2,p) ; 
          sf01  = sf01  + f0*f1 ;
          sf02  = sf02  + f0*f2; 
          sf12  = sf12  + f1*f2; 
          sf00  = sf00  + sq(f0);      
          sf11  = sf11  + sq(f1); 
          sf22  = sf22  + sq(f2); 
          scyf0  = scyf0   + cy[ic][np]*f0 ; 
          scyf1  = scyf1   + cy[ic][np]*f1;  
          scyf2  = scyf2   + cy[ic][np]*f2;  
          scayf0  = scayf0   + cay[ic][np]*f0 ; 
          scayf1  = scayf1   + cay[ic][np]*f1;  
          scayf2  = scayf2   + cay[ic][np]*f2;  
          }
        //  To deal with some singular cases, reduce matrix:
        if(sf22==0.) sf22=1.;   if(sf11==0.) sf11=1.;   if(sf00==0.) sf00=1.;
        if(ip==2 || ip==1) 
	  {D = sf00*sf11*sf22 -sf00*sf12*sf12 -sf22*sf01*sf01 -sf11*sf02*sf02 
               +2.*sf01*sf02*sf12;
           if(D !=0)
            {
	    cyfit[ip][ic][0] = ( scyf0*(sf11*sf22-sf12*sf12) -scyf1*(sf01*sf22-sf02*sf12) 
                  +scyf2*(sf01*sf12 -sf02*sf11) )/D;
	    cyfit[ip][ic][1] = ( scyf1*(sf22*sf00-sf02*sf02) -scyf2*(sf12*sf00-sf01*sf02) 
                  +scyf0*(sf12*sf02 -sf01*sf22) )/D;
 	    cyfit[ip][ic][2] = ( scyf2*(sf00*sf11-sf01*sf01) -scyf0*(sf02*sf11-sf12*sf01) 
                  +scyf1*(sf02*sf01 -sf12*sf00) )/D;
	    cayfit[ip][ic][0] = ( scayf0*(sf11*sf22-sf12*sf12) -scayf1*(sf01*sf22-sf02*sf12) 
                  +scayf2*(sf01*sf12 -sf02*sf11) )/D;
	    cayfit[ip][ic][1] = ( scayf1*(sf22*sf00-sf02*sf02) -scayf2*(sf12*sf00-sf01*sf02) 
                  +scayf0*(sf12*sf02 -sf01*sf22) )/D;
 	    cayfit[ip][ic][2] = ( scayf2*(sf00*sf11-sf01*sf01) -scayf0*(sf02*sf11-sf12*sf01) 
                  +scayf1*(sf02*sf01 -sf12*sf00) )/D;
            }
	    else cout <<"Calib-py Bad final c1 determinant  "<<endl;
          }  // --- ip==1||2 condition -------

       }  // -- end of ic loop

     // for(i_p=0; i_p<npty; i_p++) debugfile<<pyvals[i_p]<<" "<<endl;

   if(isay>0) for (np = 0; np < ndelp; np++)       // -- Check the fit.             
       {  
       if(a0vals[np] < -998.) continue; // wasn't filled
           a =   pgrid(np,delp);
	   f0   = fitpyf(ip,side,0,a) ; 
           f1   = fitpyf(ip,side,1,a) ; 
           f2   = fitpyf(ip,side,2,a) ; 
       cycalc[1] = cyfit[ip][1][0]*f0  + cyfit[ip][1][1]*f1 + cyfit[ip][1][2]*f2;
       cycalc[2] = cyfit[ip][2][0]*f0  + cyfit[ip][2][1]*f1 + cyfit[ip][2][2]*f2;	
       caycalc[1] = cayfit[ip][1][0]*f0  + cayfit[ip][1][1]*f1 + cayfit[ip][1][2]*f2;
       caycalc[2] = cayfit[ip][2][0]*f0  + cayfit[ip][2][1]*f1 + cayfit[ip][2][2]*f2;	
       if(c[1][np]!=0.)
       cout<<"y ip np a/c1,2:meas/fit "<<ip<<" "<<np<<" "<<a<<" "
        <<" "<<cy[1][np]<<"/"<<cycalc[1]<<" "<<cy[2][np]<<"/"<<cycalc[2]
        <<" "<<cay[1][np]<<"/"<<caycalc[1]<<" "<<cay[2][np]<<"/"<<caycalc[2]<<endl;
       }

      j = IP-1 + ip-1 + 2*side;
      calpay_alpha[j] = alpha;  calpay_beta[j] = beta;  calpay_gamma[j] = gamma;
      calpy_aa[j] = aa;  calpy_bb[j] = bb;  calpy_cc[j] = cc;
      calpy_cyfit[j][1][0] = cyfit[ip][1][0]; calpy_cyfit[j][2][0] = cyfit[ip][2][0]; 
      calpy_cyfit[j][1][1] = cyfit[ip][1][1]; calpy_cyfit[j][2][1] = cyfit[ip][2][1]; 
      calpy_cyfit[j][1][2] = cyfit[ip][1][2]; calpy_cyfit[j][2][2] = cyfit[ip][2][2]; 
      calpay_cayfit[j][1][0] = cayfit[ip][1][0]; calpay_cayfit[j][2][0] = cayfit[ip][2][0]; 
      calpay_cayfit[j][1][1] = cayfit[ip][1][1]; calpay_cayfit[j][2][1] = cayfit[ip][2][1]; 
      calpay_cayfit[j][1][2] = cayfit[ip][1][2]; calpay_cayfit[j][2][2] = cayfit[ip][2][2]; 

      if(isay>0)      // -- work out momentum to print it out --
      for (np=0; np<ndelp; np++)
        p = pvals[np]; if(p==0.) continue; //  wasn't filled
      for(i_p=0; i_p<npty; i_p++)
        {
        j =  np*npty*npty + npty*n0 + i_p;  
        if(pscan[j][ip][side]==0.) goto LOOP3 ;  //  wasn't filled    
        a = ayscan[j][ip][side]*1000. ;
        y = yscan[j][ip][side];
        if(a<-998.) goto LOOP3;
        disa  = a - (alpha*p + beta*sq(p) +gamma*sq(p)*p);
        disy  = yscan[j][ip][side] -  (aa*p + bb*sq(p) + cc*sq(p)*p);

	   f0   = fitpyf(ip,side,0,p) ; 
           f1   = fitpyf(ip,side,1,p) ; 
           f2   = fitpyf(ip,side,2,p) ; 
       cycalc[1] = cyfit[ip][1][0]*f0  + cyfit[ip][1][1]*f1 + cyfit[ip][1][2]*f2;
       cycalc[2] = cyfit[ip][2][0]*f0  + cyfit[ip][2][1]*f1 + cyfit[ip][2][2]*f2;	
       caycalc[1] = cayfit[ip][1][0]*f0  + cayfit[ip][1][1]*f1 + cayfit[ip][1][2]*f2;
       caycalc[2] = cayfit[ip][2][0]*f0  + cayfit[ip][2][1]*f1 + cayfit[ip][2][2]*f2;	
        py_y =  (1./cycalc[1])*disy + (1./cycalc[2])*sq(disy);
        py_ay =  caycalc[1]*disa + caycalc[2]*sq(disa);
        wy = sq(cycalc[1]) ;
        way= 1000000./(32.*sq(caycalc[1])); //  (sqrt(2)/8 m)^2.   working here in mrad.
        pynew = (wy*py_y + way*py_ay)/(wy + way );

        if(isay>0)
        debugfile<<" FitPlay"<<ip<<" side"<<side<<" "<<pyvals[i_p]<<" "<<np<<" "<<i_p
		 <<" "<<py_y<<" "<<py_ay<<" "<<pynew<< " "<<cyfit[ip][1][0]<<" "<<cyfit[ip][1][1]<<" "<<cycalc[1]<<" "<<caycalc[1]<<
          " y,a "<<yscan[j][ip][side]<<" "<<ayscan[j][ip][side]<<endl;

       LOOP3:;
        }
     } } 

     
}
// ===============================================================
double FPTrack::pgrid(int np, float delp)
{ double p = (np+1)*delp - 0.08*(np+1)*(np+1) ;
  return p;
  }
// ===============================================================
double FPTrack::fitpxf(int ip, int side, int i, double a) 
{
  double f;
  f = 0.;
    if(ip==1)
    {
      if(side==0)
      { 
      if(i==0) f=1. ;
      if(i==1) f=a;
      if(i==2) f=a*a;      
      }
      if(side==1)
      { 
      if(i==0) f=1. ;
      if(i==1) f=log(a);
      if(i==2) f=0.;      
      } 
    }
    if(ip==2)
      {
      if(side==0){    
                  if(i==0) f=1.;   
                  if(i==1 && a!=0.) f=1./(a) ;
                  if(i==2 && a!=0.) f=1./(a*a*sqrt(a)) ;
                 }
      if(side==1){    
                  if(i==0) f=1.;   
                  if(i==1 && a>0.) f=log(a) ;
                  if(i==2 && a!=0.) f= 1/a ;
                 }
      }
  return f;  
}
// ===============================================================
double FPTrack::fitpyf(int ip, int side, int i, double a) 
{
  double f;
  f = 0.;
    if(ip==1)
    { 
      if(i==0) f=1. ;
      if(i==1) f=a;
      if(i==2) f=a*a;      
    };    
    if(ip==2)
      {
      if(side==0){    
                  if(i==0) f=1.;   
                  if(i==1 && a!=0.) f=a ;
                  if(i==2 && a!=0.) f=a*a ;
                 }
      if(side==1){    
                  if(i==0) f=1.;   
                  if(i==1 && a>0.) f=a ;
                  if(i==2 && a!=0.) f= 1/(a*sqrt(a)) ;
                 }
      }

  return f;  
}
//  =============================================================
double FPTrack::fitf(int ip, int k, int i, double pval) 
{
  //  The function used in fitting the calibrations.

  double f;
  f = 0.;
  if(k==1)   //  refers to the order of the angle term.
   {
    if(ip==1)
    { 
      if(i==0) f=1./sqrt(sqrt(pval))  ;
      if(i==1) f=pval;
      if(i==2) f=pval*pval/sqrt(sqrt(pval));      
    };    
    if(ip==2)
     {    
     if(i==0) f=1.;   
     if(i==1) f=1./pval ;
     }
   }
  if(k==2 && ip==1 && i==1 ) f=pval*sqrt(pval)  ;
  return f;  
}
//  =============================================================
void FPTrack::Output(int iruncon, int n)
{ 
  float s=lengthscale, yrap, Mtrue, Del2M, fdxi, dm, dy, yrec, xi,delxi; 

  if(initaccep)
   {
     //   cout<<" initialise "<<mmean2<<endl;
    naccep  =0; naccep0=0;  naccep1=0;  naccep00=0; 
    naccep01=0; naccep2=0; naccep22=0; 
    sumvar22=0.;sumvar42=0.;sumvar44=0;sumrvar22=0.;sumrvar42=0.;sumrvar44=0; 
    nmean=0; ximean=0.; ximean2=0.; mmean=0.; mmean2=0.; ymean=0. ;ymean2=0.;
    ymmean=0.;
    initaccep = false;
   }
  if(iruncon<0||iruncon>5) cout<<"Bad run constant ="<<iruncon<<
                           " What was intended?"<<endl;

  for(int side=0; side<2; side++)
    {
    for(int ip=1; ip<=nplane; ip++)
    { if(pscan[n][ip][side] == 0.) continue;  //  didn't reach this plane.
      if(ip<3 && xscan[n][ip][side] > xmax[ip][side]) 
	            {xmax[ip][side]  = xscan[n][ip][side] ;
	             yxmax[ip][side] = yscan[n][ip][side] ;}
      if(iruncon==3 && (initialising || calibrating)) continue;
      if(ip==UTAplane) 
         if(xscan[n][ip][side]>0.001) UTAstyle->Fill(xscan[n][ip][side]*s) ;
      Xplane[ip][side]->Fill(xscan[n][ip][side]*s) ; 
      // Yplane[ip][side]->Fill(yscan[n][ip][side]*s) ;
      XYplane[ip][side]->Fill(xscan[n][ip][side]*s,yscan[n][ip][side]*s) ;
      XdX[ip][side]->Fill(xscan[n][ip][side]*s,mscale*axscan[n][ip][side]); 
    

     }}
  bool badmom=false;
   if(pmom[0]>=pbeam0 || pmom[1]>=pbeam0 
      || Precoval[1][0]>=pbeam0 ||Precoval[1][1]>=pbeam0)
     {yrap = 999., yrec = 999.; Mrecoval=0.; Mtrue=0.; Del2M=0.;badmom=true;}
   else{
   yrap     = 0.5*log((pbeam0-pmom[0])/(pbeam0-pmom[1])) ;
   yrec     = 0.5*log((pbeam0-Precoval[1][0])/(pbeam0-Precoval[1][1])) ;
   Mrecoval = sqrt(4.*(pbeam0-Precoval[1][0])*(pbeam0-Precoval[1][1])) ;
   Mtrue    = sqrt(4.*(p0[0]-pmom[0])*(p0[1]-pmom[1]));
   Del2M    = 2.*sq(sigbeam0)*(1. + 2.*sq((pmom[0]-pmom[1])/Mtrue)) ;
   // -- this formula treats the two beam uncertainties independently. 
  if(iruncon==3 && isay>0&&Mrecoval>0.) debugfile<<"Mass reconstructed "<<Mrecoval<<endl;
   
  // cout<<pmom[0]<<" "<<Precoval[1][0]<<" "<<pmom[1]<<" "<<Precoval[1][1]<<endl;

  //  Evaluate the rms differences fitted wrt input. 
  //  iplane = 1 for 215 m and 2 for 420 m
  if(
      (masscan >0 && (idetec[1][0]+idetec[2][0])*(idetec[1][1]+idetec[2][1])>3)
    || idetec[1][0]*idetec[2][1]>3
    || idetec[1][1]*idetec[2][0]>3
    || idetec[2][0]*idetec[2][1]>3 
      )                                  // i.e. we have two silicon planes hit
     { for(int side=0; side<2 ; side++)
	{
        for(int ip=1; ip<=nplane; ip++)
	 if(idetec[ip][side]>1) Yplane[ip][side]->Fill(yscan[n][ip][side]*s) ;

        //fdxi  = (pbeam0 - p0[side])/(p0[side]-pmom[side]) ;
        fdxi    = (Precoval[1][side]-pmom[side])/(p0[side]-pmom[side]) ;
	ximean  = ximean  + fdxi;	ximean2 = ximean2 + fdxi*fdxi;
        }
        nmean++ ;
        dm      = (Mrecoval - Mtrue) ;
        dy      = yrec - yrap ;
        mmean   = mmean + dm;           mmean2  = mmean2 + dm*dm;
        ymean   = ymean + dy;           ymean2  = ymean2 + dy*dy;
        ymmean  = ymmean + dm*dy;

	if(isay>1) cout<<" nmean "<<nmean<<" "<<dm<<endl;

	if(isay>1) cout<<idetec[1][0]<<idetec[2][0]<<idetec[1][1]<<idetec[2][1]<<endl;

        if( idetec[1][0]*idetec[2][1]>3 || idetec[1][1]*idetec[2][0]>3)
	  DelM42->Fill(dm);
	if(idetec[2][0]*idetec[2][1] >3 )  DelM44->Fill(dm);
        if( (idetec[1][0]+idetec[2][0])*(idetec[1][1]+idetec[2][1]) >3)
	  {DelMall->Fill(dm);
	    Delyall->Fill(dy);
            float daz = azival[0]-azival[1]-azival0[0]+azival0[1];
            if(daz>pi) daz=daz-2.*pi;  if(daz>pi) daz=daz-2.*pi;  
            if(daz<-pi) daz=daz+2.*pi;  if(daz<-pi) daz=daz+2.*pi;
            DelAzi->Fill(daz);
            if(ptval[0]>0.300 && ptval[1]>0.300)  DelAzi2->Fill(daz);
                                         //  Cut (GeV) for better azi meas.
	    for(int iside=0; iside<2; iside++) 
              {
	      if(isay>1) cout<<" fillM44 "<<pbeam0<<" "<<pmom[iside]<<endl;

               fdxi=(Precoval[1][iside]-pmom[iside])/(pbeam0 -pmom[iside]);
               xi = fabs(pbeam0 - pmom[iside])/pbeam0;
               delxi = 5./pbeam0;
               int iplot = int(xi/delxi -4.);
                    { delxi=0.01; 
                      if(xi<0.02) continue;
		      iplot = int((xi-0.02)/delxi);
		    }
               if(iplot<0) continue;
	       if(iplot >=30) iplot=29;
               xP[iplot]->Fill(fdxi);
	      } }
	//if(fabs(yrap)>1.2  && fabs(yrap)<1.5)     Delyall->Fill(dy);}
       if(nev<0)
	cout<<idetec[1][0]<<idetec[2][0]<<idetec[1][1]<<idetec[2][1]<<
	  " "<<fdxi<<" "<<nmean<<" "<<Mrecoval<<" "<<Mtrue<<" "<<dy<<" "<<mmean<<endl;
     }
   } // -- badmom else loop.

   //  If we hit one plane, the next one is now cancelled for the acceptances.
   //  This could be altered.      
   if(idetec[1][0] > 0) idetec[2][0] = 0;    
   if(idetec[1][1] > 0) idetec[2][1] = 0;
   if(idetec[2][0]>1) naccep0++ ;
   if(idetec[2][1]>1) naccep1++ ;
   if(idetec[2][0]*idetec[2][1]>3) 
       {naccep++ ; 
        if(xmax[2][0]>0.0001)XMax44->Fill(s*xmax[2][0]);
	if(xmax[2][1]>0.0001)XMax44->Fill(s*xmax[2][1]);
        if(xmax[2][0]>0.0001)YXMax44->Fill(s*yxmax[2][0]);
	if(xmax[2][1]>0.0001)YXMax44->Fill(s*yxmax[2][1]);
        yrap44->Fill(yrap);
        Mreco44->Fill(Mrecoval) ;
        sumrvar44 = sumrvar44 + 1./Del2M ;
        sumvar44  = sumvar44  +   Del2M ;
       }
   if(idetec[1][0]>1) naccep00++ ;
   if(idetec[1][1]>1) naccep01++ ;
   if(idetec[1][0]*idetec[1][1]>3) 
       {naccep22++ ;
        sumrvar22 = sumrvar22 + 1./Del2M ;
        sumvar22  = sumvar22 +      Del2M ;
       }

   if(idetec[1][0]*idetec[2][1] + idetec[2][0]*idetec[1][1]>3) 
     // if(idetec[1][0]*idetec[2][1]>3)  //  just one side combination
       {naccep2++ ; 
	 if(idetec[1][0]*idetec[2][1] > 3){
	   XMax2with4->Fill(s*xmax[1][0]);
	   XMax4with2->Fill(s*xmax[2][1]); }
	 if(idetec[1][1]*idetec[2][0] > 3){
           XMax2with4->Fill(s*xmax[1][1]);
	   XMax4with2->Fill(s*xmax[2][0]); }
          Mreco42->Fill(Mrecoval) ; 
          yrap42->Fill(yrap);
          sumrvar42 = sumrvar42 + 1./Del2M ;
          sumvar42  = sumvar42  + Del2M ;
        }
}
// ===================================================================
void FPTrack::Statistics()  
{ 
  if(outputdat) outputdatafile.close() ;
  debugfile.close() ;
  if(masscan <20) cout<<"Higgs Mass = "<<Higgsmassscale*120.<<" " ;
  if(masscan >=50) cout<<"LPAIR Mode = "<<masscan<<" " ;
  cout<<"No. events = "<<nev
           <<" Accepted at 420 = "<<naccep0<<","<<naccep1
           <<"-->coincidences = "<<naccep<<endl
           <<" Accepted at 220 = "<<naccep00<<","<<naccep01
           <<"-->coincidences ="<<naccep2<<endl;
  massvec[ivec] = Higgsmassscale*120.;
  a44vec[ivec] = float(naccep)/float(nev) ;
  a42vec[ivec] = float(naccep2)/float(nev) ;
  a22vec[ivec] = float(naccep22)/float(nev) ;
  // There are two possibilities: a weighted or unweighted mass width. 
  bool weighted = false;

  if(weighted)  
    {if(sumrvar44==0) sumrvar44=0.000001;
     if(sumrvar42==0) sumrvar42=0.000001;
     if(sumrvar22==0) sumrvar22=0.000001;
     dm44vec[ivec]   = sqrt(float(naccep)/sumrvar44) ;
     dm42vec[ivec]   = sqrt(float(naccep2)/sumrvar42) ;
     dm22vec[ivec]   = sqrt(float(naccep22)/sumrvar22) ;
     dmAllvec[ivec]  = sqrt(float(naccep+naccep2+naccep22)
			   /(sumrvar44+sumrvar42+sumrvar22)) ;
    }
  if(!weighted)  
    {
     dm44vec[ivec]   = sqrt(sumvar44/(float(naccep)+0.0001)) ;
     dm42vec[ivec]   = sqrt(sumvar42/(float(naccep2)+0.0001)) ;
     dm22vec[ivec]   = sqrt(sumvar22/(float(naccep22)+0.0001)) ;
     dmAllvec[ivec]  = sqrt((sumvar44+sumvar42+sumvar22)
			   /(float(naccep+naccep2+naccep22)+0.0001)) ;    
    }
  double sm,sy;
  corr_my[ivec] = 0.;  sigma_M44[ivec]=0.; sigma_M42[ivec]=0.; 
  sigma_mall[ivec]=0.; sigma_yall[ivec]=0.;

  float fnmean=nmean + 0.00001;
  sigma_xi[ivec]= sqrt((0.5*ximean2/fnmean) -sq(0.5*ximean/fnmean)) ;  
  sm  = sqrt( (mmean2/fnmean) -sq(mmean/fnmean)) ;  
  sy = sqrt( (ymean2/fnmean) -sq(ymean/fnmean)) ;
  if(sm*sy != 0.)
  corr_my[ivec]  = (ymmean/fnmean - mmean*ymean/sq(fnmean))/(sm*sy) ;

  Double_t par2;
  if(a44vec[ivec] > 0.) {
  DelM44->Fit("gaus","Q");
  TF1 *fit44 = DelM44->GetFunction("gaus");
  par2 = fit44->GetParameter(2);
  sigma_M44[ivec] = par2 ;}
 
  if(a42vec[ivec] > 0.) {
    DelM42->Fit("gaus","Q");
    TF1 *fit42 = DelM42->GetFunction("gaus");
    par2 = fit42->GetParameter(2);
    sigma_M42[ivec] = par2 ;}

  if(a44vec[ivec] +a42vec[ivec] + a22vec[ivec]> 0.) {
    DelMall->Fit("gaus","Q");
    TF1 *fitall = DelMall->GetFunction("gaus");
    par2 = fitall->GetParameter(2);
    sigma_mall[ivec] = par2 ;
   // 
   Delyall->Fit("gaus","Q");
   TF1 *fityall = Delyall->GetFunction("gaus");
   par2 = fityall->GetParameter(2);
   sigma_yall[ivec] = par2 ; 
  }

  // DelM42->Reset(); DelM44->Reset() ; DelMall->Reset(); Delyall->Reset(); 
}
//  =============================================================
void FPTrack::MagnetPlace(int n, int side, int type, 
			  double x, double y, double z, double xb, 
                         double length, double strength)  

   { Magnet_x[n][side] = x;
     Magnet_y[n][side] = y;    
     Magnet_z[n][side] = z;
     Magnet_xb[n][side]= xb;
     Magnet_type[n][side] = type;
     Magnet_length[n][side] = length;

     //  calculate qpole gradient from formula g = K1L Brho /L
     //                                            K1L from CERN database
     //                                            Brho = 23349. = 3.345*7000
     //                                            L = mag length

     ntype = type;
     if(type ==1) grad =  Brho*strength/length;
     if(type ==2) grad = -Brho*strength/length;
     if(type ==3 && strength >= 0.) {ntype = 1; grad =  Brho*strength/length;}
     if(type ==3 && strength <  0.) {ntype = 2; grad = -Brho*strength/length;}
 
     if(type>0) cout<<" "<<n<<" Grad "<<grad;

     Magnet_type[n][side] = ntype;
     Magnet_strength[n][side] = grad;  
     if(type==0 || type==4) Magnet_strength[n][side] = strength;  
 

     //  strength = bend angle for dipoles, -> gradient for quadrupoles
     //  type = 0 for dipole,
     //         1 for hor focussing qpoles,
     //         2 for vert focussing qpoles.
   }
//  =============================================================
void FPTrack::PlanePlace( double z, float xmin, float xmax, float ymin, float ymax)  
{ 
  std::string ia[30]  ={"0","1","2","3","4","5","6","7","8","9",
			"10","11","12","13","14","15","16","17","18","19",
                        "20","21","22","23","24","25","26","27","28","29"};
  std::string PM[2] ={"+","-"};
  std::string Xp    ="Xplane";
  std::string Yp    ="Yplane";
  std::string XYp   ="XYplanes"; 
  float xmax2;

  int nh=50, Nh=100, i, nn, nm;
  
  ymin=0.;
  nplane++; int n = nplane;
  float  d1=0.0005*mscale, d0=0., pmax, pmin;
  float s=lengthscale; //  scale, to convert from metres used in calculation.

  int nnh=nh; int nNh=Nh; if(chrom){nnh=5*nh; nNh=5*Nh;}
  Plane_z[n]= z;
  for (i=0; i<2;  i++)
    {
    nn = 0;
    for(nm=0; nm<nmagnet; nm++) {if(fabs(Magnet_z[nm][0]) > z) break; nn=nm;}
      Plane_nmag[n][i] = nn ;
    }

  if(xinnerp[n][0]<0.) xinnerp[n][0] = xinner0;
  if(xinnerp[n][1]<0.) xinnerp[n][1] = xinner0;

  cout<<nplane<<" = New plane at z = "<<z<<" at x distances "<<
                     xinnerp[n][0]<<","<<xinnerp[n][1]<<" m, after magnet "<<
                     Plane_nmag[n][0]<<", "<<Plane_nmag[n][1]<<endl; 

  if(initialising || calibrating) return ;
  if(z < 0.) cout<<" ERROR! z HAS BEEN SET NEGATIVE!"<<endl;
  pmin=0.; xmin = 0.;
  if(fabs(z) < 350.)
    {pmax = 1000. ; xmax = 4.0; ymax = 1.; d1 = 0.0003*mscale; d0=-0.0001*mscale;}
  else
    {pmax = 130. ; xmax = 2.5; ymax = 0.4;}
  

  //  Choose UTA plane
  UTAplane = 1;
  if(n==UTAplane){
  UTAmin=xinnerp[UTAplane][0]; UTAdx=0.003;
  if(UTAplane==1) UTAmin=0.003; 
  UTAstyle     = new TH1F("UTAstyle","UTAstyle",8,UTAmin*s,(UTAmin+8*UTAdx)*s);
  DXi          = new TH1F("DXi/Xi","DXi/Xi",50,0.75,1.25);
  }
  for(i=0; i<2; i++)
  {
  Xplane[n][i] = new TH1F((Xp+PM[i]+ia[n]).c_str(),(Xp+PM[i]+ia[n]).c_str(),
                                                25,xmin,xmax);
  Yplane[n][i] = new TH1F((Yp+PM[i]+ia[n]).c_str(),(Yp+PM[i]+ia[n]).c_str(),
                                                NY,-ymax,ymax);
  XYplane[n][i]= new TH2F((XYp+PM[i]+ia[n]).c_str(),(XYp+PM[i]+ia[n]).c_str(), 
                                                nh,xmin,xmax,nh,-ymax,ymax);
  Xt[n][i]     = new TH2F(("Xt"+PM[i]+ia[n]).c_str(),("Xt"+PM[i]+ia[n]).c_str(),
                                                nh,xmin,xmax,nh,0.,3.);
  XdX[n][i]   = new TH2F(("XdX"+PM[i]+ia[n]).c_str(),("XdX"+PM[i]+ia[n]).c_str(),
                                                nnh,xmin,xmax,nNh,d0,d1);
  DPvx[n][i] = new TH2F(("DPvX"+PM[i]+ia[n]).c_str(),("DPvX"+PM[i]+ia[n]).c_str(),
                                                Nh,pmin,pmax,nh,xmin,xmax);
  DPvy[n][i] = new TH2F(("DPvY"+PM[i]+ia[n]).c_str(),("DPvY"+PM[i]+ia[n]).c_str(),
			                        Nh,pmin,pmax, nh,-ymax,ymax);
  DPvP[n][i] = new TH2F(("DPvP"+PM[i]+ia[n]).c_str(),("DPvP"+PM[i]+ia[n]).c_str(),
                                                Nh,pmin,pmax,nh,pmin,pmax);
  if(n==2) Preco[n][i] = new TH1F(("Preco"+PM[i]+ia[n]).c_str(),
				 ("Preco"+PM[i]+ia[n]).c_str(),24,0.,120.);
  if(n==1) Preco[n][i] = new TH1F(("Preco"+PM[i]+ia[n]).c_str(),
				 ("Preco"+PM[i]+ia[n]).c_str(),80,0.,400.);
  DelP[n][i]  = new TH1F(("DelP"+PM[i]+ia[n]).c_str(),("DelP"+PM[i]+ia[n]).c_str(),
			                        40,-10.,10.) ;
  DelPx[n][i] = new TH1F(("DelPx"+PM[i]+ia[n]).c_str(),("DelPx"+PM[i]+ia[n]).c_str(),
                                                40,-0.4,0.4);
  DelPy[n][i] = new TH1F(("DelPy"+PM[i]+ia[n]).c_str(),("DelPy"+PM[i]+ia[n]).c_str(),
                                                40,-0.4,0.4);
  DelPt[n][i] = new TH1F(("DelPt"+PM[i]+ia[n]).c_str(),("DelPt"+PM[i]+ia[n]).c_str(),
                                                40,-0.4,0.4);
  ArcL[n][i]  = new TH1F(("ArcLength"+PM[i]+ia[n]).c_str(),
                         ("ArcLength"+PM[i]+ia[n]).c_str(),40,-2.,2.);
  DXbyX[n][i]  = new TH1F(("Angle-Mean"+PM[i]+ia[n]).c_str(),("Angle-Mean"+PM[i]+ia[n]).c_str(),
                          40,-20.,20.);

  }

  if(n==2) for(i=0; i<2; i++)
    {
    xmax2=0.4; if(z>390.) xmax2=2.5;
    x12[i] = new TH2F(("x12"+PM[i]).c_str(),("x12"+PM[i]).c_str(),
		    100,-0.25,0.25,100,-xmax2,xmax2);
    y12[i] = new TH2F(("y12"+PM[i]).c_str(),("y12"+PM[i]).c_str(),
		   100,-1.5*ymax,1.5*ymax,100,-1.5*ymax,1.5*ymax);
    xy12[i] = new TH2F(("xy12"+PM[i]).c_str(),("xy12"+PM[i]).c_str(),
		    100,-0.2,0.2,100,-1.5*ymax,1.5*ymax);

    }

  if(donexP) return;
  for(i=0; i<30; i++) 
  {xP[i]=new TH1F(("X-p"+ia[i]).c_str(),("X-p"+ia[i]).c_str(),100,-0.10,0.10);}
  donexP=true; 
}
//===========================================================================
void FPTrack::CollPlace( double z, double xap)  
{ 
 int nn,nm,i;
 ncoll++; 
 Coll_z[ncoll]   = z;
 Coll_xap[ncoll] = xap;
 for (i=0; i<2;  i++)
  {
  nn = 0;
  for(nm=0; nm<nmagnet; nm++) {if(fabs(Magnet_z[nm][0]) > z) break; nn=nm;}
    Coll_nmag[ncoll][i] = nn ;
    if(i==0)
    cout<<ncoll<<" Collimator at z = "<<z<<", x aperture +- "<<xap<<" m"<<endl; 
  }
}
//===========================================================================
void FPTrack::Sethist()
{
 page->SetFillColor(0);
 TStyle * mystyle = new TStyle("mystyle","mystyle");
 mystyle->SetCanvasBorderMode(0);  
 mystyle->SetCanvasBorderSize(0); 
 mystyle->SetPadBorderMode(0); 
 mystyle->SetStatBorderSize(1); 
 mystyle->SetPadColor(0); 
 // mystyle->SetPadLeftMargin(0.05);
 // mystyle->SetPadRightMargin(0.05);
 mystyle->SetCanvasColor(0); 
 mystyle->SetStatColor(0); 
 mystyle->SetOptStat(110010);
 mystyle->SetOptDate(11); 
 mystyle->SetDateX(0.0); 
 mystyle->SetDateY(0.0); 
 mystyle->GetAttDate()->SetTextSize(0.02); 
 mystyle->SetFrameBorderSize(1);  // -- seems to refer to the histograms
 mystyle->SetFrameFillColor(0);   // -- coloured background in histograms!
 mystyle->SetLabelSize(0.07, "X");
 mystyle->SetTitleColor(3); 
 mystyle->SetTitleBorderSize(0);
 mystyle->SetTitleFontSize(0.1);
 mystyle->SetTitleOffset(0.5);
 mystyle->SetTitleSize(1.5);
 // mystyle->SetTitleFillColor(5);   // --- requires a very recent ROOT version 
 gROOT->SetStyle("mystyle"); 

page->Divide(4,4,0.005,0.015);

page->SetTitle("LHC tracking");
 page->UseCurrentStyle();
 if(isay>0) cout<<"Set the histograms"<<endl; 
}
//======================================================================
void FPTrack::Sethist2()
{
 TStyle * mystyle2 = new TStyle("mystyle2","mystyle2");
 mystyle2->SetCanvasBorderMode(0);  
 mystyle2->SetCanvasBorderSize(0); 
 mystyle2->SetPadBorderMode(0); 
 mystyle2->SetStatBorderSize(0);   //  was 1 
 mystyle2->SetPadColor(0); 
 // mystyle2->SetPadLeftMargin(0.05);
 // mystyle2->SetPadRightMargin(0.05);
 mystyle2->SetCanvasColor(0); 
 mystyle2->SetStatColor(10); 
 // mystyle2->SetOptStat(110010);
 mystyle2->SetOptStat(0);
 mystyle2->SetOptDate(0);  
 mystyle2->SetDateX(0.0); 
 mystyle2->SetDateY(0.0); 
 mystyle2->GetAttDate()->SetTextSize(0.02); 
 mystyle2->SetFrameBorderSize(1);  // -- seems to refer to the histograms
 mystyle2->SetFrameFillColor(0);   // -- coloured background in histograms!
 mystyle2->SetLabelSize(0.05, "X");
 mystyle2->SetTitleColor(3); 
 mystyle2->SetTitleBorderSize(0);
 mystyle2->SetTitleOffset(0.75);
 mystyle2->SetTitleSize(1.5);
 // mystyle2->SetTitleStyle(3016);
 // mystyle2->SetTitleFillColor(2);   // --- requires a very recent ROOT version 
 gROOT->SetStyle("mystyle2");
 //  gROOT->ForceStyle();
page2->SetFillColor(0);
page2->UseCurrentStyle();
page2->Divide(1,1,0.005,0.005);
//page2->SetTitle(" D analysis (5.3.3)");

}
//========================================================================
void FPTrack::Drawhist()
{

  float z;
  TLatex tex;
  tex.SetTextAlign(22);
  tex.SetTextSize(0.04);

  if(!drawhists) return ;
  char * mytitle = Getline(" Input a histogram title! ");
  tex.DrawLatex(0.5,1.02,mytitle);
 
  TPostScript * ps = new TPostScript("FPTracks.ps",111);
  ps->NewPage();
  int k = 0;

  cout<<" draw hists"<<endl;
  page->Clear() ;
  page->Divide(4,4,0.005,0.015);
  for (int side=0;side<2; side++)
    {
    for(int i=1; i<=nplane; i++)
    {    
    z=xinnerp[i][side]*lengthscale;
    k++; page->cd(k) ; Xplane[i][side]->Draw(); Redline(Xplane[i][side],z,1);
    k++; page->cd(k) ; Yplane[i][side]->Draw(); 
    k++; page->cd(k) ; XYplane[i][side]->Draw();Redline(XYplane[i][side],z,1);
    k++; page->cd(k) ; XdX[i][side]->Draw(); Redline(XdX[i][side],z,1 );
    // k++; page->cd(k) ; DPvP[i][side]->Draw();

    k++; page->cd(k) ; DPvx[i][side]->Draw();
    k++; page->cd(k) ; DPvy[i][side]->Draw(); 
    // k++; page->cd(k) ; Preco[i][side]->Draw();
    }
    k++; page->cd(k) ; DelM42->Draw();
    k++; page->cd(k) ; DelMall->Draw();

   if(side==0) {k++; page->cd(k); Magfail->Draw();}  // just draw once
   if(side==0) {k++; page->cd(k); Mreco44->Draw();}
   if(side==1) {k++; page->cd(k); Mreco42->Draw();}
   page->Update();
   // if(side==1) continue ;
   ps->NewPage(); k=0;
   page->UseCurrentStyle();
   }
    page->Clear(); 
   page->Divide(6,5,0.005,0.01) ;
   for(k=1; k<=30; k++) {page->cd(k); xP[k-1]->Draw();}  
   page->Update() ;
   // ps.NewPage();
   ps->Close();
}
//****************************************************************************
void FPTrack::Drawhist2()
{
  int kk;
  // float z;
  TLatex tex;             TLatex tex2;
  tex.SetTextAlign(11);   tex2.SetTextAlign(11);
  tex.SetTextSize(0.05);  tex.SetTextSize(0.03);
  
  if(!drawhists || (calibrating && !chrom)) return ;
  gROOT->SetStyle("mystyle2");
  gROOT->ForceStyle()  ;
  TH2F * Hnew1 = new TH2F(*XYplane[1][0]);
  Hnew1->SetMarkerSize(3.0);
  page2->Divide(2,1,0.005,0.015); page2->Clear();
  page2->cd(1); Hnew1->Draw("AXIS");  Hnew1->Draw("SAME") ;
  // page2->Print("XYPlane1.eps");  
  TH2F * Hnew1a = new TH2F(*XYplane[2][1]);
  page2->cd(1); Hnew1a->Draw("AXIS");  Hnew1a->Draw("SAME") ;
   page2->Print("XYPlane2.eps");

   float lpairxsec=1.;
   if(masscan==56) lpairxsec = 4.12;    //  mumu 6 GeV pt
   if(masscan==54) lpairxsec = 10.73;   //   mumu 4 GeV pt 


  for(int ip=1; ip<=2; ip++) for(int side=0; side<2; side++)
    {TH1F * Hnew2x = new TH1F(*Xplane[ip][side]);
     TH1F * Hnew2y = new TH1F(*Yplane[ip][side]);
      if(masscan>=50) *Hnew2x = *Hnew2x * (lpairxsec*1000./float(nev));
      if(masscan>=50) *Hnew2y = *Hnew2y * (lpairxsec*1000./float(nev));
     page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
     page2->SetLogy(1) ;
     Hnew2x->GetXaxis()->SetNdivisions(1603,kTRUE) ; 
     Hnew2x->GetXaxis()->SetTitle("X from beam line (cm)") ;
     float maxbin = 0;
     for(Int_t ibin=1; ibin<=16; ibin++){
       int nbi =  (int)Hnew2x->GetBinContent(ibin); // cout<<nbi<<endl;
       if (nbi > maxbin) maxbin = nbi ;
        }  
     page2->cd(1); Hnew2x->Draw("AXIS");  Hnew2x->Draw("SAME") ;
     float z=xinnerp[ip][side]*lengthscale;
     Redline(Hnew2x, z, 2)   ;
     if(ip==1&&side==0) page2->Print((thisIP+"X10.eps").c_str());
     if(ip==1&&side==1) page2->Print((thisIP+"X11.eps").c_str());
     if(ip==2&&side==0) page2->Print((thisIP+"X20.eps").c_str());
     if(ip==2&&side==1) page2->Print((thisIP+"X21.eps").c_str());
     page2->SetLogy(0) ;

     page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
     Hnew2y->GetXaxis()->SetNdivisions(808,kTRUE) ; 
     Hnew2y->GetXaxis()->SetTitle("Y from beam line (cm)") ;
           maxbin = 0;
     for(Int_t ibin=1; ibin<=16; ibin++){
       int nbi =  (int)Hnew2y->GetBinContent(ibin); // cout<<nbi<<endl;
       if (nbi > maxbin) maxbin = nbi ;
        }  
     page2->cd(1); Hnew2y->Draw("AXIS");  Hnew2y->Draw("SAME") ;
     if(ip==1&&side==0) page2->Print((thisIP+"Y10.eps").c_str());
     if(ip==1&&side==1) page2->Print((thisIP+"Y11.eps").c_str());
     if(ip==2&&side==0) page2->Print((thisIP+"Y20.eps").c_str());
     if(ip==2&&side==1) page2->Print((thisIP+"Y21.eps").c_str());
     page2->SetLogy(0) ;

     std::cout<<" Y plane "<<ip<<" side "<<side<<std::endl;
     for(int ibin=1; ibin<=NY; ibin++) std::cout<<Hnew2y->GetBinContent(ibin)<<" ";
     std::cout<<endl; 


    }


/*
  for(int ip=1; ip<=2; ip++) for(int side=0; side<2; side++)
    {
     TH2F * Hnew4 = new TH2F(*XYplane[ip][side]);
     Hnew4->SetMarkerSize(3.0);
     page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
     page2->cd(1); Hnew4->Draw("BOX");  
     Hnew4->GetXaxis()->SetNdivisions(1604,kTRUE) ; 
     Hnew4->GetXaxis()->SetTitle("X from beam line (cm)") ;
     Hnew4->GetYaxis()->SetNdivisions(1604,kTRUE) ; 
     Hnew4->GetYaxis()->SetTitle("Y from beam line (cm)") ;
     float z=xinnerp[ip][side]*lengthscale;
     Redline(Hnew4, z, 2)   ;
     if(ip==1&&side==0) page2->Print((thisIP+"XY10.eps").c_str());
     if(ip==1&&side==1) page2->Print((thisIP+"XY11.eps").c_str());
     if(ip==2&&side==0) page2->Print((thisIP+"XY20.eps").c_str());
     if(ip==2&&side==1) page2->Print((thisIP+"XY21.eps").c_str());
    }
*/
    
    kk=0; 
    page2->Clear(); 
    page2->UseCurrentStyle();
    page2->SetGrid();
    page2->Divide(2,2,0.01,0.02);
    for(int ip=1; ip<=2; ip++) for(int side=0; side<2; side++)
    {
     TH2F * Hnew3 = new TH2F(*XdX[ip][side]);
     Hnew3->SetMarkerSize(3.0);
     kk++;
     page2->cd(kk); 
     Hnew3->GetXaxis()->SetNdivisions(1603,kTRUE) ; 
     Hnew3->GetXaxis()->SetTitle("X from beam line (cm)") ;
     Hnew3->GetYaxis()->SetNdivisions(1605,kTRUE) ; 
     Hnew3->GetYaxis()->SetTitle("dx/dz (mrad)") ;
     if(chrom) Hnew3->Draw() ; else Hnew3->Draw("BOX");
     float z=xinnerp[ip][side]*lengthscale;
     Redline(Hnew3,z,3 );
    } 
    page2->Update();
    page2->Print("XdX.eps");

/*
   for(int ip=1; ip<=2; ip++) for(int side=0; side<2; side++)
    {
     TH1F * Hnew6 = new TH1F(*Preco[ip][side]);
     if(masscan>=50) *Hnew6 = *Hnew6 * (lpairxsec*1000./float(nev));
     Hnew6->SetMarkerSize(3.0);
     page2->Divide(1,1,0.01,0.02); page2->Clear();page2->UseCurrentStyle();
     page2->cd(1); Hnew6->Draw("AXIS"); Hnew6->Draw("SAME") ;  
     Hnew6->GetXaxis()->SetNdivisions(1603,kTRUE) ; 
     Hnew6->GetXaxis()->SetTitle("Reconstructed Momentum (GeV)") ;
     if(ip==1 && side==0) page2->Print("Prec10.eps");
     if(ip==2 && side==0) page2->Print("Prec20.eps");
     if(ip==1 && side==1) page2->Print("Prec11.eps");
     if(ip==2 && side==1) page2->Print("Prec21.eps");
    } 

  // TH1F * Hnew5 = new TH1F(*Mreco);

  TH1F * Hnew9a  = new TH1F(*Mreco44);
  TH1F * Hnew9b  = new TH1F(*Mreco42);
  Hnew9a->SetTitle("Reconstructed mass, 420+420 ") ;  
  Hnew9b->SetTitle("Reconstructed mass, 420+220 ") ;  
  page2->Clear(); page2->Divide(2,2,0.005,0.015);page2->UseCurrentStyle();
  Hnew9a->SetLabelSize(0.06,"X");
  Hnew9b->SetLabelSize(0.06,"X");
  page2->cd(1); Hnew9a->Draw("AXIS");  Hnew9a->Draw("SAME") ;
  // double xx = page2->cd(1)->GetUxmin() , yy = gPad->GetUymax();
  tex2.DrawLatex(meanm,1.06*Hnew9a->GetMaximum()," 420+420 ");
  page2->cd(2); 
  Hnew9b->Draw("AXIS");  Hnew9b->Draw("SAME") ;
  tex2.DrawLatex(meanm,1.06*Hnew9b->GetMaximum()," 420+220 ");
  page2->Print("MassRecon.eps");
  if(masscan==0){cout<<" Fit mass plot "<<endl;
    Hnew9a->Fit("gaus"); Hnew9b->Fit("gaus");
  }

  page2->Clear(); page2->Divide(2,2,0.005,0.015);page2->UseCurrentStyle();  
  kk=0;
  for(int is=0; is<2; is++) 
  for(int ip=1; ip<=2; ip++)
  {
  TH1F * HnewA  = new TH1F(* DelPx[ip][is]);
  kk++;
  page2->cd(kk);
  HnewA->Draw();
  HnewA->GetXaxis()->SetTitle("Del(pX) (GeV)") ;
  }
  page2->Print("DelpX.eps");
  page2->Clear(); page2->Divide(2,2,0.005,0.015);page2->UseCurrentStyle();  
  kk=0;
  for(int is=0; is<2; is++) 
  for(int ip=1; ip<=2; ip++)
  {
  TH1F * HnewA  = new TH1F(* DelPy[ip][is]);
  kk++;  page2->cd(kk);
  HnewA->Draw();
  HnewA->GetXaxis()->SetTitle("Del(py) (GeV)") ;
  }
  page2->Print("DelpY.eps")
*/
  page2->Clear(); page2->Divide(2,2,0.005,0.015);page2->UseCurrentStyle();  
  kk=0;
  for(int is=0; is<2; is++) 
  for(int ip=1; ip<=2; ip++)
  {
  TH1F * HnewA  = new TH1F(* Preco[ip][is]);
  kk++;  page2->cd(kk);
  HnewA->Draw();
  HnewA->GetXaxis()->SetTitle("Dp = p0.xi (GeV)") ;
  }
  page2->Print("P0xi.eps");

  page2->Clear(); page2->Divide(1,2,0.005,0.015);page2->UseCurrentStyle();  
  kk=0;
  for(int is=0; is<2; is++) 
  {
  int ip=2;
  TH1F * HnewA  = new TH1F(* DXbyX[ip][is]);
  kk++;  page2->cd(kk);
  HnewA->Draw();
  HnewA->GetXaxis()->SetTitle("Dx (murad) - mean ") ;
  }
  page2->Print("DX-mean.eps");

/*
  page2->Clear(); page2->Divide(2,2,0.005,0.015);page2->UseCurrentStyle();  
  kk=0;
  for(int is=0; is<2; is++) 
  for(int ip=1; ip<=2; ip++)
  {
  TH1F * HnewA  = new TH1F(* DelPt[ip][is]);
  kk++;
  page2->cd(kk);
  HnewA->Draw();
  HnewA->GetXaxis()->SetTitle("Del(pt) (GeV)") ;
  }
  page2->Print("DelpT.eps");
*/

  page2->Clear(); page2->Divide(1,1,0.005,0.015);page2->UseCurrentStyle();  
  TH1F * HnewY  = new TH1F(*yrap42);
  page2->cd(1);
  HnewY->Draw();
  HnewY->GetXaxis()->SetTitle("Higgs rapidity") ;
  page2->Print("HiggsY.eps");



  page2->Clear(); page2->Divide(2,2,0.005,0.015);page2->UseCurrentStyle();  
  kk=0;
  for(int is=0; is<2; is++) 
  for(int ip=1; ip<=2; ip++)
  {
  TH1F * HnewA  = new TH1F(* ArcL[ip][is]);
  kk++;
  page2->cd(kk);
  HnewA->Draw();
  HnewA->GetXaxis()->SetTitle("Arc Variation (mm)") ;
  }
  page2->Print("ArcVariation.eps");

  page2->Clear(); page2->Divide(1,1,0.005,0.015);page2->UseCurrentStyle();  
  TH1F * HnewA  = new TH1F(* DelAzi);
  page2->cd(1);
  HnewA->SetMinimum(0.);
  HnewA->Draw();
  HnewA->GetXaxis()->SetTitle("Delta(Azimuth) ") ;
  page2->Print("DelAzi.eps");
  page2->Clear(); page2->Divide(1,1,0.005,0.015);page2->UseCurrentStyle();  
  TH1F * HnewB  = new TH1F(* DelAzi2);
  page2->cd(1);
  HnewB->SetMinimum(0.);
  HnewB->Draw();
  HnewB->GetXaxis()->SetTitle("Delta(Azimuth) ") ;
  page2->Print("DelAzi2.eps");
  page2->Clear(); page2->Divide(1,1,0.005,0.015);page2->UseCurrentStyle();  
 TH1F * HnewC  = new TH1F(* UTAstyle);
  page2->cd(1);
  HnewC->SetMinimum(0.);
  HnewC->Draw();
  HnewC->GetXaxis()->SetTitle("cm ") ;
  page2->Print("UTAstyle.eps");

  /*
  Hnew9b->SetLabelSize(0.06,"X");
  %page2->cd(1); Hnew9a->Draw("AXIS");  Hnew9a->Draw("SAME") ;
  // double xx = page2->cd(1)->GetUxmin() , yy = gPad->GetUymax();
  tex2.DrawLatex(meanm,1.06*Hnew9a->GetMaximum()," 420+420 ");
  page2->cd(2); 
  Hnew9b->Draw("AXIS");  Hnew9b->Draw("SAME") ;
  tex2.DrawLatex(meanm,1.06*Hnew9b->GetMaximum()," 420+220 ");
  page2->Print("MassRecon2.eps");
 
  page2->SetLogy(0) ;
  DelM42->GetXaxis()->SetTitle("(GeV)") ;
  page2->cd(2); DelM42->Draw();  
  tex.DrawLatex(0. ,10.,"420+220 ");
  page2->Update() ;
  page2->Print("DelM42.eps");
  */

  page2->Clear();page2->UseCurrentStyle();
  page2->Divide(1,2,0.005,0.025); 
  for(kk=1; kk<=2; kk++)
    {
  page2->cd(kk);

  x12[kk-1]->GetXaxis()->SetTitle("First plane     cm") ;
  x12[kk-1]->GetYaxis()->SetTitle("Second plane    cm") ;
  x12[kk-1]->SetLabelSize(0.04,"X");
  x12[kk-1]->SetMarkerStyle(20);
  x12[kk-1]->SetMarkerSize(0.3);
  x12[kk-1]->Draw();
    }
  page2->Update();
  page2->Print("X2_vs_X1.eps");
  page2->Clear();
  page2->Divide(1,2,0.005,0.025); 
  for(kk=1; kk<=2; kk++)
    {
  page2->cd(kk);
  y12[kk-1]->SetLabelSize(0.04,"X");
  y12[kk-1]->GetXaxis()->SetTitle("First plane  y  cm") ;
  y12[kk-1]->GetYaxis()->SetTitle("Second plane y  cm") ;
 
  y12[kk-1]->SetMarkerStyle(20);
  y12[kk-1]->SetMarkerSize(0.3);
  y12[kk-1]->Draw();
    }
  page2->Update();
  page2->Print("Y2_vs_Y1.eps");
  page2->Clear();
  page2->Divide(1,2,0.005,0.025); 
  for(kk=1; kk<=2; kk++)
    {
  page2->cd(kk);
  xy12[kk-1]->SetLabelSize(0.04,"X");
  xy12[kk-1]->GetXaxis()->SetTitle("First plane  x  cm") ;
  xy12[kk-1]->GetYaxis()->SetTitle("Second plane y  cm") ;
  xy12[kk-1]->SetMarkerStyle(20) ;
  xy12[kk-1]->SetMarkerSize(0.3);
  xy12[kk-1]->Draw("scat");
    }
  page2->Update();
  page2->Print("Y2_vs_X1.eps");

  page2->Close();
}
//==========================================================
void FPTrack::Redline(TH1F & hist, float & cut, int kcol)
{
  float ymax = 1.05*hist.GetMaximum();
  float ymin = 1.05*hist.GetMinimum();
   TLine* line = new TLine(cut, ymin, cut, ymax ) ;
   line->SetLineColor(kcol);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
}
//==========================================================
void FPTrack::Redline(TH1F * hist, float & cut, int kcol)
{
  float ymax = 1.0*hist->GetMaximum();
  float ymin = 1.0*hist->GetMinimum();
   TLine* line = new TLine(cut, ymin, cut, ymax ) ;
   line->SetLineColor(kcol);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
}
//==========================================================
void FPTrack::Redline(TH2F & hist, float & cut, int kcol)
{
  float ymax = 1.0*hist.GetMaximum();
  float ymin = 1.0*hist.GetMinimum();
   TLine* line = new TLine(cut,ymin, cut, ymax ) ;
   line->SetLineColor(kcol);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
}
//==========================================================
void FPTrack::Redline(TH2F * hist, float & cut, int kcol)
{
  float ymax = 1.0*hist->GetMaximum();
  float ymin = 1.0*hist->GetMinimum();
   TLine* line = new TLine(cut, ymin, cut, ymax ) ;
   line->SetLineColor(kcol);
   line->SetLineWidth(2);
   line->SetLineStyle(2);
   line->Draw();
}