#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)
#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[2] ; TH1F * Mreco44;  TH1F * Mreco42;    
TH1F * UTAstyle;        TH1F * DXi      ; TH1F * xP[20] ;      
   
TH2F * DPvP[51][2] ;     TH2F * DPvx[51][2] ; TH2F * DPvy[51][2] ; 
TH1F * Magfail ;        
TH1F * yrap44 ;  TH1F * yrap42 ; 
TH1F * XMax44 ;  TH1F * XMax4with2 ; TH1F * XMax2with4 ;
TH1F * YXMax44 ;  TH1F * YXMax4with2 ; TH1F * YXMax2with4 ;
TH1F * DelM44 ; TH1F * DelM42 ; TH1F * DelMall ;
TH1F * Delyall ;

TRandom3 gau ;

// 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[35][2], Plane_nmag[51][2], Magnet_aperclass[35][2], 
       ncoll=0, Coll_nmag[51][2],  
       symmetrise=0, ivec=0, nev, coltinue, masscan, icali ;
double Magnet_x[35][2],Magnet_y[35][2],Magnet_z[35][2],Magnet_length[35][2], 
       Magnet_strength[35][2], grad, 
       Magnet_aperA[35][4][2],  
       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,
       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 
       sigdxdz = 0., sigx = 0.,
       apermb = 0. ,
       xcol1=999., xcol2=999., Coll_z[5], Coll_xap[5], xcoll0[5][2],
       meanm, ximean,ximean2,mmean,mmean2,ymean,ymean2, ymmean;
float  cut[20], ptval, dptval, mscale=1000.,  // for milliradians 
       aside[2]={-1.,1.},
       Higgsmassscale = 1. ,  p0_subgev=0. , 
       xinner0= 0.003, xinnerp[51] = {-1., -1., -1., -1., -1.},
       x0IP=0., y0IP=0., ax0IP=0., ay0IP=0., xmax[10][2], yxmax[10][2],
       pscan[3100][5][2], xscan[3100][5][2], axscan[3100][5][2],
       yscan[3100][5][2], ayscan[3100][5][2], pmom[2], p0[2],
       sumrvar42, sumrvar44, sumrvar22, sumvar42, sumvar44, sumvar22 ,
       massvec[20],a42vec[20],a44vec[20], a22vec[20], dm44vec[20],dm42vec[20], 
       dm22vec[20], dmAllvec[20],sigma_xi[20], sigma_m[20],sigma_y[20],
       sigma_M44[20],sigma_M42[20], sigma_mall[20],sigma_yall[20],
       corr_my[20], readplane[51]; 
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,
       useaper, drawhists, psmear=false, x0smear=false,  donexP=false,
       plotcoin = false,  initmess=true, initaccep=true, outputdat=false ;
Int_t  ncand, nmagnet=0, ntype;
int    isay=1, irtype, nplane=0, iplane, naccep,naccep0,naccep1, nmean, 
       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 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};
std::ofstream outputdatafile ;
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

  Preco[0]    = new TH1F("P-reco+","P-reco+",60,0.,120.);
  Preco[1]    = new TH1F("P-reco-","P-reco-",60,0.,120.);
  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 y, float dirx, float diry)
{
  double pnew,disa,pprev,c1calc,c2calc, precon, xval,aval;
  double asmear, xsmear;

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

  if(y+diry==-99.) cout<<endl; // dummy commmand to remove warnings

    //-- apply smearing if required ---
    asmear = 0. ; xsmear = 0.;
    if(sigdxdz>0) asmear = sigdxdz*gau.Gaus();  
    if(sigx>0)    xsmear = sigx*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;

    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);
      }

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

    if( !initialising && !calibrating) Preco[side]->Fill(pnew);
    return precon ;
}
//=======================================================================
void FPTrack::ReadData()
{
  // Read data cards from Data.txt file 

  std::string dataname ;
  std::ifstream mydata;
  cout<<"Open data file"<<endl; 
  mydata.open("FPData.txt");
  bool readon=true ;
  int item ;
  ireadp=0;

  IP = 0 ;
 
  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") xinnerp[1] = float(item)*0.000001;
     else if(dataname == "SIM2") xinnerp[2] = 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 == "CALI") icali =  item ;
     else if(dataname == "PSIG") psmear = (item>0) ;
     else if(dataname == "DXIP") x0smear = (item>0) ;

     else if(dataname == "ASIG") sigdxdz = item*0.000001 ;
     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 == "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") ;
  if(outputdat) cout<<" Opened fptraj "<<endl;

  if(irtype<3) masscan=0 ;
  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);

}
//=======================================================================
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[35][2],magcen[35][2],magstrength[35][2],
         L,K0L,K1L,K2L,K3L, HKICK,VKICK, 
         dx, dy, dist, magaper[35][4][2] ;
  int    magnet=0, maxmag=0, max1mag=0, minmag=0, magver, 
         side, magtype[35][2], magapertype[35][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

   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.tfs");
       if(ifile==0 && IP==5) MyMags.open("LHCB2IR5.tfs");
       if(ifile==1 && IP==1) MyMags.open("LHCB1IR1.tfs");
       if(ifile==1 && IP==5) MyMags.open("LHCB1IR5.tfs");
      }             // Also version 6.5 at present (Sep 06)


    cout<<"Open magnets file"<<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)
    {
      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;

      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(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. ; 
        }
      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 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 ;

     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
      
      MagnetPlace(mag,is,magtype[minmag+mag][i], dx,dy,
                aside[j]*(magcen[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()
{
  CollPlace(149., xcol1) ;
  CollPlace(184., xcol2) ;
}
//===========================================================================

void FPTrack::Run()
{
//---------- Steering routine for the analysis. ------------    
   float h1=0.,h2=0.,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 ;
                  }   

   //-- Calibrate the momentum measurement if asked for.
   if(irtype ==3 && docali) { cout<<"Calibrating"<<endl;
                   calibrating = true ; 
                   Tracks(2) ;
                   calibrating = false ;
                  }       

   //-- Do tracking for input data file. Possible mass scan.
   if(irtype!=3) masscan = 0 ;
   if(irtype==2) calibrating = docali;
   if(masscan==0) Tracks(irtype) ;
   else 
     {if(masscan==1){h1=0.4; h2=2.0; 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++ ;
         }
       cout<<"Double_t mass["<<ivec<<"]={";
       for(int i=0; i<ivec; i++) cout<<massvec[i]<<", "; 
          cout<<endl; 
       if(masscan!=2){
          cout<<"Double_t a44["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<a44vec[i]<<", "; 
	  cout<<endl;}
       cout<<"Double_t a42["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<a42vec[i]<<", "; 
          cout<<endl; 
       cout<<"Double_t a22["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<a22vec[i]<<", "; 
          cout<<endl; 
       cout<<"Double_t dm42["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<dm42vec[i]<<", "; 
          cout<<endl; 
       if(masscan!=2){
          cout<<"Double_t dm44["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<dm44vec[i]<<", ";
          cout<<endl;}
       cout<<"Double_t dm22["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<dm22vec[i]<<", ";           
          cout<<endl; 
       if(masscan!=2){
          cout<<"Double_t sigma_M44["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<sigma_M44[i]<<", "; 
          cout<<endl; } 
       cout<<"Double_t sigma_M42["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<sigma_M42[i]<<", "; 
          cout<<endl; 
       cout<<"Double_t sigma_mall["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<sigma_mall[i]<<", "; 
          cout<<endl; 
       cout<<"Double_t sigma_yall["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<sigma_yall[i]<<", "; 
          cout<<endl;       
       cout<<"Double_t corr_my["<<ivec<<"]={"; 
          for(int i=0; i<ivec; i++) cout<<corr_my[i]<<", "; 
          cout<<endl;       
       cout<<"Xp error widths "<<endl;
       for (int i=0; i<20; 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,np,nptx,npty,iptx,ipty, 
         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, 
         p,tandir,xc,yc,zc,dxyz[3],dir[3],
         dzc,xp,yp,zp, dx0ip, theta=0., phi=0.;
  float  ppx[2],ppy[2],ppz[2], xangle, yangle, s;
  float  dphi=0., delp=0., xdisp, ydisp, themr, thx,thy; 
  float  xposition, yposition, xwrtbeam=-999.; double zposition;
  bool   outside, coutside; 

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

  p0[0] = pbeam0 ; p0[1] = pbeam0 ;
   
// Define the events for iruncon<=2 case.
   if(iruncon==1) numevents=1;
   if(iruncon==2)                   // calibration and scans
     {
     dptval  = 1.5/30.;
     nptx     = 30;     //<-- should be no more than 30
     npty     = 30 ;
     npt=nptx;
     if(docali) {npty=1; nptx=2*nptx +1;}  
                            // Calibration needs p:ptx grid only
     nthphi  = nptx*npty;
     ndelp   = 50;    
     delp    = 20. ;
     if(masscan==1) delp = 20. ;   
     if(masscan==2) delp = 20.; 
     dphi    = twopi/nphi;
     numevents=nthphi*ndelp;

     if(numevents >3100) {cout<<"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. ;
       }}
     }   
   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"<<endl; 
     MyInput.open("p1p2data");
     if(MyInput.eof()) cout<<" Empty Input File!"<<endl ;
     numevents = 100000000 ;
     initmess = false ;
     }
   
   bool p0smear = !initialising && !calibrating && psmear ;
   bool ipsmear = !initialising && !calibrating && x0smear ;
//------- Process events according to chosen method:-
   for (nev=0; nev <numevents; nev++)
     {
     if(nev > 0) initialising = false ; // first event did this.
     if(isay>0 && iruncon==3) cout<<"Event "<<nev<<" "<<endl; 
     if(p0smear) {
                  p0[0] = pbeam0 + gau.Gaus()*sigbeam0 ;
                  p0[1] = pbeam0 + gau.Gaus()*sigbeam0 ;
                 } 
     dx0ip = 0.;
     if(ipsmear) dx0ip =  sigx0*gau.Gaus() ;  // common for both sides


     if(iruncon == 1) 
       {
       pmom[0] = p0[0] - p0_subgev;
       pmom[1] = p0[1] - p0_subgev ; 
       cout<<"Beam momentum "<<pmom[0]<<endl; 
       theta=0;
       phi=0;
       }
     else if(iruncon == 2 || iruncon == 4 || iruncon == 5)
       {
        if(iruncon ==2 )
        {nnn  =  nev/npty;
         np   =  nnn/nptx; 
         iptx = nnn - nptx*np;
         ipty =  nev - npty*nnn; 
         thx = (npt - iptx)*dptval; 
         thy = (float(npty)*0.5  - 0.5 - float(ipty))*dptval;
         phi = atan2(thy,thx);
         ptval = sqrt(thy*thy + thx*thx);
         ipt = iptx ;
        }
        else  
       {
        nnn  =  nev/nphi;
        np   =  nnn/npt; 
        ipt = nnn - npt*np;
        iphi =  nev - nphi*nnn; 
        phi     = dphi*iphi  ;
        ptval   = dptval*ipt ; 
        nnn     = nev/nthphi + ipt ;
        if(ipt == nphi/2) nnn=nnn+npt*ndelp ;
       }
       if(calibrating) 
            {pmom[0] = p0[0] - pgrid(np,delp) ;
	     pmom[1] = p0[1] - pgrid(np,delp) ; }
       else {pmom[0]= p0[0] -  (np+1)*delp;
	     pmom[1]= p0[1] -  (np+1)*delp;
            }
       }
       else if(iruncon == 3)
       { 
       if(MyInput.eof()) break ;
       MyInput >>ppx[0]>>ppy[0]>>ppz[0]>>ppx[1]>>ppy[1]>>ppz[1];
       if(isay>1) cout<<"Input "<<ppx[0]<<" "<<ppy[0]<<" "<<ppz[0]<<" "
                             <<ppx[1]<<" "<<ppy[1]<<" "<<ppz[1]<<endl; 

//--- Rescale the mass.  But we want to keep PT the same.
       for(int side=0; side<2; side++)
            {ppx[side]=ppx[side]*1.;
             ppy[side]=ppy[side]*1.;
             ppz[side]=p0[side]-(pbeam0-fabs(ppz[side]))*Higgsmassscale ;
	    } 
       }   
//----- Loop over the two sides
     for(int side = 0; side<2; side++)
       { 
	 Precoval[1][side] = pbeam0;   //adjust first argument if > 1 track.
       // set initial position, angle, momentum.  This x is positive inwards! 

       if(iruncon == 2 || iruncon == 4 || iruncon == 5)
	 theta      = ptval/pmom[side] ; 
       if(iruncon == 3 && !initialising)
         { 
         ptval      = sqrt( sq(ppx[side]) + sq(ppy[side])) ;
         pmom[side] = sqrt( sq(ptval)     + sq(ppz[side])) ; 
         phi        = atan2(ppy[side],ppx[side]);  
         theta      = ptval/pmom[side] ;  
      }   

       idetec[1][side]=0; idetec[2][side]=0;
       xp = -x0IP  + dx0ip;   // second term can allow for beam spot size.
       yp =  y0IP ;
       zp =  0. ;
       xmax[1][side] = 0. ;  xmax[2][side] = 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.
       dir[0]  = dir[0] - ax0IP ;
       dir[1]  = dir[1] + ay0IP ;

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

       //if(iruncon==3 && nev==26) abort(); //FOR NOW

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

       outside = false ;  coutside  = false;
       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) cout<<magnet<<" IP "
	     <<dxyz[0]<<" "<<dxyz[1]<<" "<<dir[0]<<" "<<dir[1]<<" p="<<p<<endl;

         if(isay>2)cout<<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] ;

	// 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) cout<<magnet <<" EN "<<dxyz[0]<<" "
                            <<dxyz[1]<<" "<<zc+dzc+dxyz[2]<<endl;
        if(isay > 1 && outside)
          {cout<<" - outside - "<<magnet<<endl;
	   if(isay > 2)
	    cout<<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.

          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) cout<<magnet<<" DE " ;
          if(isay>1 && i==1) cout<<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.
	  }
        }
	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)cout<<magnet<<" QE ";

	    }     // end of quadrupole treatment 
        if(isay>1)cout<<dxyz[0]<<" "<<dxyz[1]<<" "
		    <<0.5* Magnet_length[magnet][side]*aside[side]+zc
                    <<" "<<dir[0]<<" "<<dir[1]<<endl;
  
        dxyz[2]  = 0.5* Magnet_length[magnet][side]*aside[side] ;   

//  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)
          {cout<<" - outside - "<<magnet<<endl;
	   if(isay > 2)
	    cout<<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) cout<<magnet<<" "<<zp<<" "<<xp<<" "<<yp<<endl;

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

        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(initialising) 
             {
              xcoll0[ic][side]  = xposition;
	      if(isay>0)cout<<"INIC "<<ic<<" "<<side<<" "<<xposition<<endl;
	     } 
           else
	     {
	      xwrtbeam  = xposition - xcoll0[ic][side] ;
              coutside  = coutside || fabs(xwrtbeam) > Coll_xap[ic] ;
              if(coltinue==0) outside = coutside ;
	     }
	  } 
        if(outside) 
	  {if(isay>1)cout<<"- hit collimator - "<<magnet<<" "<<xwrtbeam<<endl;
	   break;}

// 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];
           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] ; 
	       if(isay>0) cout<<"INIP "<<ip<<" "<<side<<" "<<xwrtbeam<<" "
                                  <<dir[0]<<endl;
	       goto ENDPLANES ;
             } //------------------------end of initialising.

           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.) {xdisp  = -xdisp  ;
	                             xangle = -xangle ;}  // +ve into Si detector.
           if( (iruncon==4 && ip==1) || (iruncon==5 && ip==2))
	      {
	       if(xdisp > xinnerp[ip])   iaccep[np][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 = nnn ; 
               if(n>=3100) goto ENDPLANES; 
	      }
           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;
           if(xdisp > xinnerp[ip] ) // Si width.  PARTICLE ACCEPTED.
	     {     
             Precoval[1][side]=Pformula(IP,ip,side,xdisp,ydisp,xangle,yangle);
             idetec[ip][side]=1 ;
             DPvP[ip][side]->Fill(p0[side]-Precoval[1][side],p0[side]-pmom[side]);
             DPvx[ip][side]->Fill(p0[side]-Precoval[1][side],xdisp*s);
             DPvy[ip][side]->Fill(p0[side]-Precoval[1][side],ydisp*s);
             }
	   // 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>0 && (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,ptval,themr,phi,xdisp,ydisp,xposition,yposition,zposition,xangle,p,Precoval[1][side], np);

	    ENDPLANES:;
	   
	    //cout<<xplane0[1][side]<<" "<<xplane0[2][side]<<endl;
       }  // ---plane loop
       if(outside) break;
     }   // ---magnet loop
   }  // --- side loop

    if(isay>0 && iruncon==3) cout<<" 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)) Output(iruncon,n) ; 
 }  //event loop 
  nev--;  

  if(iruncon==2 && !initialising) Calibrate(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++)
     {cout<<"IP="<<IP<<" xinner(mm)="<<xinnerp[iruncon-3]
                  <<" Side="<<isi<<" Plane="<<iruncon-3<<endl;
         for(imom=0; imom<200; imom++)
          {for(ipt=0; ipt<10; ipt++)cout<<iaccep[imom][ipt][isi]<<" ";
	    cout<<endl;
    } } } 
  MyInput.close() ;   // Close the input file. 
}
//==============================================================
void FPTrack::Calibrate(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 ;   

   nptx = 2*npt + 1 ;
   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++) {
     
     //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 (j=npt; j <ndelp*nptx; j=j+nptx)  
        {
	float a, p, x;
        np = j/nptx ;
        pvals[np]  = p0[side] - pscan[j][ip][side];
        if(pscan[j][ip][side]==0.) pvals[np] = 0.; 
        if(pvals[np]==0.) continue ;  // wasn't filled    
 
        a = axscan[j][ip][side] ;
        p = pvals[np] ;
        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 <<" Bad x  determinant  "<<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; j <(np+1)*nptx; j++)
       {  
        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 + npt) cout<< "centre ";
        if(isay>0) cout<<"np/p/disp/disa "<<np<<" "<<pval<<" "<<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<<"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 <<" 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; j++)
        {
        if(pscan[j][ip][side]==0.) goto LOOP3 ;  // wasn't filled    
        np  = j/nptx ; 
        if(np==npt) continue ;
        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);
          }
        cout<<" 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 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;   
}
//===============================================================
double FPTrack::pgrid(int np, float delp)
{ double p = (np+1)*delp - 0.08*(np+1)*(np+1) ;
  return p;
  }
//===============================================================
double FPTrack::fitf(int ip, int k, int i, double pval) 
{
  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; 

  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(xscan[n][ip][side]>0.003) 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]); 
      if(ip<3 && xscan[n][ip][side] > xmax[ip][side]) 
	            {xmax[ip][side]  = xscan[n][ip][side] ;
	             yxmax[ip][side] = yscan[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(isay>0&&Mrecoval>0.)cout<<"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])>0)
    || idetec[1][0]*idetec[2][1]>0
    || idetec[1][1]*idetec[2][0]>0
    || idetec[2][0]*idetec[2][1]>0 
     )
     { for(int side=0; side<2 ; side++)
	{

      for(int ip=1; ip<=nplane; ip++)
	if(idetec[ip][side]>0) 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]>0 || idetec[1][1]*idetec[2][0]>0)
	  DelM42->Fill(dm);
	if(idetec[2][0]*idetec[2][1] >0 )  DelM44->Fill(dm);
        if( (idetec[1][0]+idetec[2][0])*(idetec[1][1]+idetec[2][1]) >0)
	  {DelMall->Fill(dm);
	    Delyall->Fill(dy);
	    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]);
		//int iplot = int(fabs(pbeam0-pmom[iside])/20. -1.);
		//int iplot = int(fabs(pbeam0-pmom[iside])/10. -2.);
		int iplot = int(fabs(pbeam0-pmom[iside])/5. -4.);
               if(iplot<0 || iplot>20) continue;
	       //if(iplot >=10) iplot=10+int((fabs(pbeam0-pmom[iside]) -200.)/100.) ;
	       if(iplot >=20) iplot=19;
               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(idetec[1][0] > 0) idetec[2][0] = 0;    
   if(idetec[1][1] > 0) idetec[2][1] = 0;
   if(idetec[2][0]>0) naccep0++ ;
   if(idetec[2][1]>0) naccep1++ ;
   if(idetec[2][0]*idetec[2][1]>0) 
       {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]>0) naccep00++ ;
   if(idetec[1][1]>0) naccep01++ ;
   if(idetec[1][0]*idetec[1][1]) 
       {naccep22++ ;
        sumrvar22 = sumrvar22 + 1./Del2M ;
        sumvar22  = sumvar22 +      Del2M ;
       }
   if(idetec[1][0]*idetec[2][1] + idetec[2][0]*idetec[1][1]>0) 
       {naccep2++ ; 
	 if(idetec[1][0]*idetec[2][1] > 0){
	   XMax2with4->Fill(s*xmax[1][0]);
	   XMax4with2->Fill(s*xmax[2][1]); }
	 if(idetec[1][1]*idetec[2][0] > 0){
           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() ;
  cout<<"Higgs Mass = "<<Higgsmassscale*120.<<" " ;
  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 length, double strength)  

   { Magnet_x[n][side] = x;
     Magnet_y[n][side] = y;    
     Magnet_z[n][side] = z;
     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[20]  ={"0","1","2","3","4","5","6","7","8","9",
                    "10","11","12","13","14","15","16","17","18","19"};
  std::string PM[2] ={"+","-"};
  std::string Xp    ="Xplane";
  std::string Yp    ="Yplane";
  std::string XYp   ="XYplanes"; 

  int nh=50, Nh=100, i, nn, nm;

  nplane++; int n = nplane;
  float  d1=0.00005*mscale, d2=0., pmax, pmin;
  if(n==2)  {d1=0.0005*mscale;}
  float s=lengthscale; // scale, to convert from metres used in calculation.

  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.) xinnerp[n] = xinner0;

  cout<<nplane<<" = New plane at z = "<<z<<" at x distance "<<
                     xinnerp[n]<<" 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 = 2.0; ymax = 1.;}
  else
    {pmax = 130. ; xmax = 2.5; ymax = 0.3;}

  UTAstyle     = new TH1F("UTAstyle","UTAstyle",nh,xmin*s,xmax*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(),
                                                nh,xmin,xmax);
  Yplane[n][i] = new TH1F((Yp+PM[i]+ia[n]).c_str(),(Yp+PM[i]+ia[n]).c_str(),
                                                nh,-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(),
                                                nh,xmin,xmax,Nh,d2,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(donexP) return;
  for(i=0; i<20; 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 ;
      cout<<ncoll<<"Collimator at z = "<<z<<" at x aperture +- "<<xap<<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("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]*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[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(4,5,0.005,0.01) ;
   for(k=1; k<=20; k++) {page->cd(k); xP[k-1]->Draw();}  
   page->Update() ;
   //ps.NewPage();
}
//****************************************************************************
void FPTrack::Drawhist2()
{
  float z;
  TLatex tex;             TLatex tex2;
  tex.SetTextAlign(11);   tex2.SetTextAlign(11);
  tex.SetTextSize(0.05);  tex.SetTextSize(0.03);
  
  if(!drawhists || calibrating) 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");

  TH1F * Hnew2 = new TH1F(*Xplane[1][0]);
  Hnew2->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew2->Draw("AXIS");  Hnew2->Draw("SAME") ;  
  page2->Print("IP1X10.eps");
  TH1F * Hnew2a = new TH1F(*Xplane[1][1]);
  Hnew2a->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew2a->Draw("AXIS");  Hnew2a->Draw("SAME") ;  
  page2->Print("IP1X11.eps");
  TH1F * Hnew2b= new TH1F(*Xplane[2][0]);
  Hnew2b->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew2b->Draw("AXIS");  Hnew2b->Draw("SAME") ;  
  page2->Print("IP1X20.eps");

  TH1F * Hnew2c = new TH1F(*Xplane[2][1]);
  Hnew2c->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->SetLogy(0) ;
  Hnew2c->GetXaxis()->SetNdivisions(1603,kTRUE) ; 
  Hnew2c->GetXaxis()->SetTitle("X from beam line (m)") ;
  float maxbin = 0;
  for(Int_t ibin=1; ibin<=16; ibin++){
    int nbi =  (int)Hnew2c->GetBinContent(ibin); //cout<<nbi<<endl;
    if (nbi > maxbin) maxbin = nbi ;
  }  
  page2->cd(1); Hnew2c->Draw("AXIS");  Hnew2c->Draw("SAME") ;  
  page2->Print("IP1X21.eps");
  page2->SetLogy(0) ;

  TH2F * Hnew4 = new TH2F(*XYplane[1][0]);
  Hnew4->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew4->Draw("AXIS");  Hnew4->Draw("SAME") ;  
  page2->Print("IP1XY10.eps");
  TH2F * Hnew4a = new TH2F(*XYplane[1][1]);
  Hnew4a->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew4a->Draw("AXIS");  Hnew4a->Draw("SAME") ;  
  page2->Print("IP1XY11.eps");
  TH2F * Hnew4b= new TH2F(*XYplane[2][0]);
  Hnew4b->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew4b->Draw("AXIS");  Hnew4b->Draw("SAME") ;  
  page2->Print("IP1XY20.eps");
  TH2F * Hnew4c = new TH2F(*XYplane[2][1]);
  Hnew4c->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew4c->Draw("AXIS");  Hnew4c->Draw("SAME") ;  
  page2->Print("IP1XY21.eps");

  TH2F * Hnew3 = new TH2F(*XdX[2][0]);
  Hnew3->SetMarkerSize(3.0);
  page2->Divide(1,1,0.01,0.02); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew3->Draw("AXIS");  Hnew3->Draw("SAME") ;  
  z=xinnerp[2]*lengthscale;
  Redline(Hnew3,z,3 );
  page2->Print("IP1XdX0.eps");
  TH2F * Hnew3a = new TH2F(*XdX[2][1]);
  Hnew3a->SetMarkerSize(3.0);
  page2->Divide(1,1,0.01,0.02); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew3a->Draw("AXIS");  Hnew3a->Draw("SAME") ;  
  Redline(Hnew3a,z,3 );
  page2->Print("IP1XdX1.eps");

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

  TH1F * Hnew6 = new TH1F(*XMax44);
  Hnew2->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew6->Draw("AXIS");  Hnew6->Draw("SAME") ;  
  page2->Print("Xmaxplot44.eps");
  TH1F * Hnew6a = new TH1F(*XMax4with2);
  Hnew2->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew6a->Draw("AXIS");  Hnew6a->Draw("SAME") ;  
  page2->Print("Xmaxplot4with2.eps");
  TH1F * Hnew6b = new TH1F(*XMax2with4);
  Hnew2->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew6b->Draw("AXIS");  Hnew6b->Draw("SAME") ;  
  page2->Print("Xmaxplot2with4.eps");
  TH1F * Hnew7 = new TH1F(*YXMax44);
  Hnew7->SetMarkerSize(3.0);
  page2->Divide(1,1,0.005,0.015); page2->Clear();page2->UseCurrentStyle();
  page2->cd(1); Hnew7->Draw("AXIS");  Hnew7->Draw("SAME") ;  
  page2->Print("YXmaxplot44.eps");

  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->UseCurrentStyle();
  page2->Divide(2,2,0.005,0.015);
  //page2->cd(1)->SetLogy(1) ;
  page2->SetLogy(1); 
  Hnew2b->GetYaxis()->SetTitle("Events "); 
  Hnew2b->GetXaxis()->SetTitle("Horiz. distance from beam line (cm)");
  *Hnew2b = *Hnew2b + *Hnew2c ; 
  TH1F * Hnew6u = new TH1F(*UTAstyle);
  
  Hnew2b->Draw("AXIS");  Hnew2b->Draw("SAME") ;
  tex.DrawLatex(2. ,1550.," (c) ");
  //page2->SetLogy(0) ;
  Hnew6u->GetYaxis()->SetTitle("Events "); 
  Hnew6u->GetXaxis()->SetTitle("Horiz. distance from beam line (cm)"); 

  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->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();
}