/// Sujeewa Kumaratunga's code for FGD-TPC alignment
/// Comparison between two different types of runs, ex: B=0 and B!=0
/// To be used with the output of TPCFGD_Align.cxx 
/// Place outputfiles from one type (ex. B!=0) in something similar to RunListB.log
/// Place outputfiles from the other type (ex. B=0) in something similar to RunListNoB.log
/// several *dataalignconst.eps output files; you can change the syntax of this last part of the string by changing oFile

#include <iostream>
#include <iomanip>
#include <vector>
#include <string>
#include <cmath>
#include <vector>
#include <stdlib.h>
#include <fstream>
#include <sstream>
#include <stdio.h>

#include "TGraph.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPaveText.h"
#include "TLegend.h"
#include "TH2.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TH1.h"
#include "TH3.h"
#include "TMath.h"
#include "TPostScript.h"
#include <TTree.h>
#include "TAxis.h"
#include "TFile.h"
#include "TChain.h"
#include "TROOT.h"
#include "TEllipse.h"
#include "TLine.h"
#include "TChain.h"
#include "TTree.h"
#include "TF1.h"
#include "TMultiGraph.h"
#include "TGraphErrors.h"
#include <TSpline.h>
#include "TMVA/VariablePCATransform.h" //for tmva
#include "TMVA/TSpline1.h"
#include "TGraphAsymmErrors.h"
#include "TLegend.h"
#include "TSystem.h"
#include "TROOT.h"

#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TF1.h>
#include <TProfile.h>
#include <TTree.h>
#include <unistd.h>
#include "t2kalignvar.C"
#include "t2kalignvart.C"

using namespace std;

void tpcfgdalign();


int FGD(int z){

  if( z > 100. && z < 300. ) return -1;  // begin fgd 1 

  if( z > 300.  && z < 500 ) return 1;  // end fgd 1

  if( z > 1400. && z < 1600 ) return -2; // begin fgd 2 

  if( z > 1600. && z < 1900 ) return 2; // end fgd 2 

}

int TPC(int z){

  if( z < 100. ) return 1; 

  if( z < 1400. ) return 2; 

  return 3; 

}

TH1F *colorhists(TH1F *hist, int i){
  hist->SetLineColor(i+1);
  hist->SetMarkerColor(i+1);
  hist->SetMarkerStyle(i+19);

  return hist;
}

TH1F *fit1d(TH1F *hist, int i, char func[256], double lolim, double hilim)
{
  hist->SetLineColor(i+1);
  hist->SetMarkerColor(i+1);
  hist->SetMarkerStyle(i+19);

  char fitname[256];
  strcpy (fitname,hist->GetName());
  strcat (fitname,"Fit");



  TF1 *fitfunc = new TF1 (fitname,func, lolim, hilim);
  fitfunc->SetLineColor(i+1);
  hist->Fit(fitname,"0EMR");
  if(fitfunc->GetNDF()>0) hist->GetFunction(fitname)->ResetBit(TF1::kNotDraw); // make "fitname" visible
  
  return hist;

}

TH1F  *fitslices(TH2F *hist, int i){

  hist->FitSlicesY(); 

  char namemean[256];
  char namerms[256];

  sprintf(namemean,"%s_1",hist->GetName()); 
  sprintf(namerms,"%s_2",hist->GetName()); 

  
  TH1F *hfit =  (TH1F*)gDirectory->Get(namemean);
  TH1F *hfitrms =  (TH1F*)gDirectory->Get(namerms);

  hfit->SetLineColor(i+1);
  hfit->SetMarkerColor(i+1);
  (hfit->GetYaxis())->SetRangeUser(-5.,5.); 


  for(Int_t bi=1; bi <= hfit->GetNbinsX(); bi++)
    {
      //to set the rms of the gaussian fit as the error
      //hfit->SetBinError(bi,hfitrms->GetBinContent(bi));
      //to set sqrt(n) as the percentage error
      Double_t binxcont = 0;
      for(Int_t by=0; by<=hist->GetNbinsY(); by++)
   {
     binxcont += hist->GetBinContent(bi,by);
   }
      hfit->SetBinError(bi,((hfitrms->GetBinContent(bi))/sqrt(binxcont)));
      cout<<"see here "<<hist->GetName()<<"\t"<<binxcont<<"\t"<<hfitrms->GetBinContent(bi)<<endl;
    }

  return hfit;

}

TH1F  *fitslicesAng(TH2F *hist, int i){

  hist->FitSlicesY(); 

  char namemean[256];
  char namerms[256];

  sprintf(namemean,"%s_1",hist->GetName()); 
  sprintf(namerms,"%s_2",hist->GetName()); 

  
  TH1F *hfit =  (TH1F*)gDirectory->Get(namemean);
  TH1F *hfitrms =  (TH1F*)gDirectory->Get(namerms);

  hfit->SetLineColor(i+1);
  hfit->SetMarkerColor(i+1);
  (hfit->GetYaxis())->SetRangeUser(-1.,1.); 


  for(Int_t bi=1; bi <= hfit->GetNbinsX(); bi++)
    {
      //to set the rms of the gaussian fit as the error
      //hfit->SetBinError(bi,hfitrms->GetBinContent(bi));
      //to set sqrt(n) as the percentage error
      Double_t binxcont = 0;
      for(Int_t by=0; by<=hist->GetNbinsY(); by++)
   {
     binxcont += hist->GetBinContent(bi,by);
   }
      hfit->SetBinError(bi,((hfitrms->GetBinContent(bi))/sqrt(binxcont)));
      cout<<"see here "<<hist->GetName()<<"\t"<<binxcont<<"\t"<<hfitrms->GetBinContent(bi)<<endl;
    }

  return hfit;

}

TH1F  *fitslices2(TH2F *hist, int i){

  hist->FitSlicesY(); 

  char name[256];

  sprintf(name,"%s_1",hist->GetName()); 

  
  TH1F *hfit =  (TH1F*)gDirectory->Get(name);
  hfit->SetLineColor(i+1);
  (hfit->GetYaxis())->SetRangeUser(-1000.,1000.); 

  return hfit;

}

bool FileExists(const char *name){
  struct stat stFileInfo;
  int intStat;

  // Attempt to get the file attributes
  intStat = stat(name,&stFileInfo); 
  if(intStat == 0) 
    return true;
  else{
    return false;
  }
}


TLegend *legshape1(int nhist, TH1F **hh, std::string hisname[]){
 
  TLegend *legorig = new TLegend(0.5,0.8,0.8,1.);
  legorig->SetBorderSize(0);
  legorig->SetFillColor(0);
  legorig->SetTextSize(0.035);

  legorig->AddEntry(hh[0]->GetName());
  for(int hi=0; hi<nhist; hi++)
    {
      
      legorig->AddEntry(hh[hi],
      Form((hisname[hi]).c_str()),
      "l");
    }

  return legorig;
}


TLegend *legshape2(int nhist, TH1F **hh, std::string hisname[]){
 
  TLegend *legorig = new TLegend(0.1,0.8,0.9,0.9);
  legorig->SetBorderSize(0);
  legorig->SetFillColor(0);
  legorig->SetTextSize(0.035);

  legorig->AddEntry(hh[0]->GetName());
  for(int hi=0; hi<nhist; hi++)
    {
      legorig->AddEntry(hh[hi],
      Form((hisname[hi]+" Entries=%-6.0f  #mu=%-6.2f #pm %-3.2f  #sigma=%4.2f #pm %-3.2f").c_str(),
           hh[hi]->GetEntries(),
           hh[hi]->GetMean(),hh[hi]->GetMeanError(),
           hh[hi]->GetRMS(),hh[hi]->GetRMSError()),
      "l");
    }

  return legorig;
}





TLegend *legshape3(int nhist, TH1F **hh, std::string hisname[]){
 
  TLegend *legorig = new TLegend(0.1,0.8,0.994,0.994);
  legorig->SetBorderSize(0);
  legorig->SetFillColor(0);
  //legorig->SetTextSize(0.035);
  legorig->SetTextSize(0.0485523);
  legorig->SetTextFont(42);

  //legorig->AddEntry(hh[0]->GetName());
  for(int hi=0; hi<nhist; hi++)
    {
      char fitname[256];
      strcpy (fitname,hh[hi]->GetName());
      strcat (fitname,"Fit");
      //cout<<"and the fitnames from hscale "<<fitname<<endl;

      legorig->AddEntry(hh[hi],
      Form((hisname[hi]+" incpt=%-6.2f #pm %-3.2f slope=%-6.4f #pm %-3.4f #frac{#chi^{2}}{ndf}=#frac{%-3.2f}{%-3.2i}  ").c_str(),
           hh[hi]->GetFunction(fitname)->GetParameter(0),
           hh[hi]->GetFunction(fitname)->GetParError(0),
            hh[hi]->GetFunction(fitname)->GetParameter(1),
            hh[hi]->GetFunction(fitname)->GetParError(1),
           hh[hi]->GetFunction(fitname)->GetChisquare(),
           hh[hi]->GetFunction(fitname)->GetNDF()),
      "l");

    }

  return legorig;
}

TH1F *hscale(int nhist, TH1F **hh)
{
  Double_t maxhist = 0.;
  for(int hi=0; hi<nhist; hi++)
    {
      if(maxhist < hh[hi]->GetMaximum()) maxhist=hh[hi]->GetMaximum();
    }

  
  for(int hi=0; hi<nhist; hi++)
    {
      if(hh[hi]->GetMaximum()>0.)
  {
    //hh[hi]->Scale(maxhist/hh[hi]->GetMaximum());
    char fitname[256];
    strcpy (fitname,hh[hi]->GetName());
    strcat (fitname,"Fit");
  }
      hh[hi]->GetYaxis()->SetRangeUser(1.,maxhist*1.15);
    } 
  return *hh;  
}

void tpcfgdalign(){

  gROOT->ProcessLine(".L t2kalignvar.C++"); 
  gROOT->ProcessLine(".L t2kalignvart.C++"); 
    
  //gStyle->SetOptFit(11111111); 
  gStyle->SetOptFit(11111111); 
  //gStyle->SetOptStat(11111111); 


  TString oFile("dataalignconst");
  

#define ANG 1.0
#define BIN 10

  // Define histograms

  std::string tpcname[3];
  for (Int_t tii= 0; tii<3; tii++)
    {
      std::ostringstream eee;
      eee << tii;
      tpcname[tii]=eee.str();
    }
  std::string fgdname[2];
  for (Int_t fii= 0; fii<2; fii++)
    {
      std::ostringstream eee;
      eee << fii;
      fgdname[fii]=eee.str();
    }
  std::string detname[6];
  //   for (Int_t dii= 0; dii<5; dii++)
  //     {
  //       std::ostringstream eee;
  //       eee << dii;
  //       detname[dii]=eee.str();
  //     } 
  detname[0]="TPC1";
  detname[1]="FGD1";
  detname[2]="TPC2";
  detname[3]="FGD2";
  detname[4]="TPC3";
  detname[5]="DSECAL";

  TFile * out;  // canvas output file
  out = new TFile("tpcfgdalignhisto.root","recreate");

  vector <TH1F* > redchitpctrk; 
  vector <TH1F* > redchitpctrk2; 

  for (Int_t tii=0; tii<3; tii++)
    {
      redchitpctrk.push_back(new TH1F (("redchitrk"+tpcname[tii]).c_str(), ("reduced chi of tracks in TPC"+tpcname[tii]).c_str(), 100,0.,10.));
      redchitpctrk2.push_back(new TH1F (("redchitrk2"+tpcname[tii]).c_str(), ("reduced chi of tracks in TPC"+tpcname[tii]).c_str(), 100,0.,10.));

      //temp vectors for 2d arrays
      //vector <TH1F* > nytemp;


    }

  //loop over 6 detector parts, creating 5 pairs

  //for first file
  vector<TH1F* > ny;
  vector<TH1F* > nx;
  vector<TH1F* > dy;
  vector<TH1F* > dx;
  vector<TH2F* > dyvsy;
  vector<TH2F* > dyvsyM;
  vector<TH2F* > dyvsyP;
  vector<TH2F* > dyvsx;
  vector<TH2F* > constdyvsy;
  vector<TH2F* > constdyvsx;
  vector<TH2F* > mupulldyvsy;
  vector<TH2F* > mupulldyvsyty;
  vector<TH2F* > mupulldyvsx;
  vector<TH2F* > dxvsy;
  vector<TH2F* > dxvsx;
  vector<TH2F* > constdxvsy;
  vector<TH2F* > constdxvsx;
  vector<TH2F* > mupulldxvsy;
  vector<TH2F* > mupulldxvsx;
  vector<TH2F* > mupulldxvsxtx;
  vector<TH2F* > tyvsy;
  vector<TH2F* > tyvsyA;
  vector<TH2F* > tyvsyPrim;
  vector<TH2F* > tyvsx;
  vector<TH2F* > txvsy;
  vector<TH2F* > txvsx;
  vector<TH2F* > dyvsty;
  vector<TH2F* > dyvstx;
  vector<TH2F* > dxvstx;
  vector<TH2F* > dxvsty;
  vector<TH1F* > fgdy;
  vector<TH1F* > fgdx;
  vector<TH1F* > tyTotDet;
  //vector<TH2F* > tysy;

  for (Int_t dii=0; dii<5; dii++)
    {
      ny.push_back(new TH1F (("ny"+detname[dii]+detname[dii+1]).c_str(), ("ny ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 11,0.,10.));
      nx.push_back(new TH1F (("nx"+detname[dii]+detname[dii+1]).c_str(), ("nx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 11,0.,10.));
      dy.push_back(new TH1F (("dy"+detname[dii]+detname[dii+1]).c_str(), ("dy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 100,-50.,50.));
      dx.push_back(new TH1F (("dx"+detname[dii]+detname[dii+1]).c_str(), ("dx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 100,-50.,50.));
  
      dyvsy.push_back(new TH2F (("dyvsy"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dyvsyM.push_back(new TH2F (("dyvsyM"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y for -X("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dyvsyP.push_back(new TH2F (("dyvsyP"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y for +X("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dyvsx.push_back(new TH2F (("dyvsx"+detname[dii]+detname[dii+1]).c_str(), ("dy vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdyvsy.push_back(new TH2F (("constdyvsy"+detname[dii]+detname[dii+1]).c_str(), ("dyvsy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdyvsx.push_back(new TH2F (("constdyvsx"+detname[dii]+detname[dii+1]).c_str(), ("dy vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldyvsy.push_back(new TH2F (("mupulldyvsy"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldyvsyty.push_back(new TH2F (("mupulldyvsyty"+detname[dii]+detname[dii+1]).c_str(), ("dy vs yty ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-700.,700.,20,-20.,20.));
      mupulldyvsx.push_back(new TH2F (("mupulldyvsx"+detname[dii]+detname[dii+1]).c_str(), ("dy vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      
      dxvsy.push_back(new TH2F (("dxvsy"+detname[dii]+detname[dii+1]).c_str(), ("dx vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dxvsx.push_back(new TH2F (("dxvsx"+detname[dii]+detname[dii+1]).c_str(), ("dx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdxvsy.push_back(new TH2F (("constdxvsy"+detname[dii]+detname[dii+1]).c_str(), ("dxvsy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdxvsx.push_back(new TH2F (("constdxvsx"+detname[dii]+detname[dii+1]).c_str(), ("dx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldxvsy.push_back(new TH2F (("mupulldxvsy"+detname[dii]+detname[dii+1]).c_str(), ("dx vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldxvsx.push_back(new TH2F (("mupulldxvsx"+detname[dii]+detname[dii+1]).c_str(), ("dx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldxvsxtx.push_back(new TH2F (("mupulldxvsxtx"+detname[dii]+detname[dii+1]).c_str(), ("dx vs xtx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));

      tyvsy.push_back(new TH2F (("tyvsy"+detname[dii]+detname[dii+1]).c_str(), ("ty vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 35,-1200.,1200.,20,-4.,4.));
      tyvsyA.push_back(new TH2F (("tyvsyA"+detname[dii]+detname[dii+1]).c_str(), ("All ty vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 35,-1200.,1200.,20,-1.,1.));
      tyvsyPrim.push_back(new TH2F (("tyvsyPrim"+detname[dii]+detname[dii+1]).c_str(), ("Primary: ty vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 35,-1200.,1200.,20,-4.,4.));
      tyvsx.push_back(new TH2F (("tyvsx"+detname[dii]+detname[dii+1]).c_str(), ("ty vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-1.,1.));
      txvsy.push_back(new TH2F (("txvsy"+detname[dii]+detname[dii+1]).c_str(), ("tx vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-1.,1.));
      txvsx.push_back(new TH2F (("txvsx"+detname[dii]+detname[dii+1]).c_str(), ("tx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-1.,1.));

      dyvsty.push_back(new TH2F (("dyvsty"+detname[dii]+detname[dii+1]).c_str(), ("dy vs ty ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));
      dxvstx.push_back(new TH2F (("dxvstx"+detname[dii]+detname[dii+1]).c_str(), ("dx vs tx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));
      dxvsty.push_back(new TH2F (("dxvsty"+detname[dii]+detname[dii+1]).c_str(), ("dx vs ty ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));
      dyvstx.push_back(new TH2F (("dyvstx"+detname[dii]+detname[dii+1]).c_str(), ("dy vs tx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));

      fgdy.push_back(new TH1F (("fgdy"+detname[dii]+detname[dii+1]).c_str(), ("fgdy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.));
      fgdx.push_back(new TH1F (("fgdx"+detname[dii]+detname[dii+1]).c_str(), ("fgdx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.));

      tyTotDet.push_back(new TH1F (("tyTotDet"+detname[dii]+detname[dii+1]).c_str(), ("total ty ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 25,-1.,1.));

      //tysy.push_back(new TH2F (("tyvsSy"+detname[dii]+detname[dii+1]).c_str(), ("ty vs Sy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-1200.,1200.));

    }

  TH1F *hmom = new TH1F("hmom","hmom",1000, -1000., 1000.);
  TH1F *hmom2 = new TH1F("hmom2","hmom2",1000, -1000., 1000.);


  TH1F *nyT = new TH1F("nyT","ny",21, 0, 20);
  
  TH1F *nxT = new TH1F("nxT","nx",21, 0, 20);
  
  TH1F *dyT = new TH1F("dyT","dy",100.,-50.,50.);
  TH1F *SyTf1 = new TH1F("SyTf1","Sy FGD1",200,-961.,961.);
  TH1F *SxTf1 = new TH1F("SxTf1","Sx FGD1",200,-961.,961.);
  TH1F *SyTf2 = new TH1F("SyTf2","Sy FGD2",200,-961.,961.);
  TH1F *SxTf2 = new TH1F("SxTf2","Sx FGD2",200,-961.,961.);
  TH1F *SyzT = new TH1F("SyzT","Syz",3600.,0.,3600.);
  TH1F *SxzT = new TH1F("SxzT","Sxz",3600.,0.,3600.);
  TH2F *newstddy = new TH2F("newstddy","dy std vs mine B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest1 = new TH2F("newstddytest1","dy : std vs mine B>0 dir-y>0, rho>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest2 = new TH2F("newstddytest2","dy : std vs mine B>0 dir-y>0, rho<=0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest3 = new TH2F("newstddytest3","dy : std vs mine B>0 dir-y<0, rho>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest4 = new TH2F("newstddytest4","dy : std vs mine B>0 dir-y<0, rho<=0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest5 = new TH2F("newstddytest5","dy : std vs mine B>0 ^z>0, rho>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest6 = new TH2F("newstddytest6","dy : std vs mine B>0 ^z>0, rho<=0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest7 = new TH2F("newstddytest7","dy : std vs mine B>0 ^z<0, rho>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest8 = new TH2F("newstddytest8","dy : std vs mine B>0 ^z<0, rho<=0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest9 = new TH2F("newstddytest9","dy : std vs mine B>0 yPositiveRoot",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddytest10 = new TH2F("newstddytest10","dy : std vs mine B>0 yNegativeRoot",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddy2 = new TH2F("newstddy2","dy : std vs mine B=0",1800,-900.,900.,1800,-900.,900.); 
  TH1F *deltady = new TH1F("deltady","extrapolated dy mine-std B>0",180,-50.,90.); 
  TH1F *deltadx = new TH1F("deltadx","extrapolated dx mine-std B>0",180,-50.,90.); 
  TH1F *deltady2 = new TH1F("deltady2","extrapolated dy mine-std B=0",180,-50.,50.); 
  TH1F *deltadx2 = new TH1F("deltadx2","extrapolated dx mine-std B=0",180,-50.,50.); 
  TH2F *deltadyy = new TH2F("deltadyy","extrapolated dy mine-std VS std y B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *deltadyx = new TH2F("deltadyx","extrapolated dy mine-std VS std x B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *deltadxy = new TH2F("deltadxy","extrapolated dx mine-std VS std y B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *deltadxx = new TH2F("deltadxx","extrapolated dx mine-std VS std x B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *deltadyy2 = new TH2F("deltadyy2","extrapolated dy mine-std VS std y B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *deltadyx2 = new TH2F("deltadyx2","extrapolated dy mine-std VS std x B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *deltadxy2 = new TH2F("deltadxy2","extrapolated dx mine-std VS std y B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *deltadxx2 = new TH2F("deltadxx2","extrapolated dx mine-std VS std x B>0",1800,-900.,900.,1800,-900.,900.); 
  newstddy->GetXaxis()->SetTitle("dy: my calculation (mm)");
  newstddy->GetYaxis()->SetTitle("dy: std calculation (mm)");
  newstddy2->GetXaxis()->SetTitle("dy: my calculation (mm)");
  newstddy2->GetYaxis()->SetTitle("dy: std calculation (mm)");
  TH2F *newstddx = new TH2F("newstddx","dx std vs mine B>0",1800,-900.,900.,1800,-900.,900.); 
  TH2F *newstddx2 = new TH2F("newstddx2","dx std vs mine B=0",1800,-900.,900.,1800,-900.,900.); 
  newstddx->GetXaxis()->SetTitle("dx: my calculation (mm)");
  newstddx->GetYaxis()->SetTitle("dx: std calculation (mm)");
  newstddx2->GetXaxis()->SetTitle("dx: my calculation (mm)");
  newstddx2->GetYaxis()->SetTitle("dx: std calculation (mm)");

  TH2F *hndofty = new TH2F("hndofty","y-angle VS ndof in track ",100,-2.,2.,100,0.,100.); 

  TH1F *dxT = new TH1F("dxT","dx",100.,-50.,50.);
  TH1F *dxN = new TH1F("dxN","dx (x<0)",100.,-50.,50.);
  TH1F *dxP = new TH1F("dxP","dx (x>0)",100.,-50.,50.);

  TH1F *dx3eN = new TH1F("dx3eN","dx (TPC3-DSECAL) (x<0)",100.,-100.,100.);
  TH1F *dx3eP = new TH1F("dx3eP","dx (TPC3-DSECAL) (x>0)",100.,-100.,100.);
  
  TH2F *dxvsxT = new TH2F("dxvsxT","dx vs x",10,-900.,900.,20,-20.,20.); 
  TH2F *dxvsyT = new TH2F("dxvsyT","dx vs y",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsxT = new TH2F("dyvsxT","dy vs x",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsyT = new TH2F("dyvsyT","dy vs y",10,-900.,900.,20,-20.,20.); 
  
  TH1F *tyT = new TH1F("tyT","ty",BIN,-ANG,ANG);
  TH1F *txT = new TH1F("txT","tx",BIN,-ANG,ANG);
  TH1F *txN = new TH1F("txN","tx (x<0)",100.,-ANG,ANG);
  TH1F *txP = new TH1F("txP","tx (x>0)",100.,-ANG,ANG);
  TH1F *tyeT = new TH1F("tyeT","ty (DSECAL)",100.,-ANG,ANG);
  TH1F *txeT = new TH1F("txeT","tx (DSECAL)",100.,-ANG,ANG);
  TH1F *txeN = new TH1F("txeN","tx (DSECAL) (x<0)",100.,-ANG,ANG);
  TH1F *txeP = new TH1F("txeP","tx (DSECAL) (x>0)",100.,-ANG,ANG);
  
  TH2F *txvsxT = new TH2F("txvsxT","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *txvsyT = new TH2F("txvsyT","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *tyvsxT = new TH2F("tyvsxT","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *tyvsyT = new TH2F("tyvsyT","  ",10,-900.,900.,20,-20.,20.); 
 
  TH2F *dyvstxT = new TH2F("dyvstxT","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dyvstyT = new TH2F("dyvstyT","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstxT = new TH2F("dxvstxT","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstyT = new TH2F("dxvstyT","  ",20,-ANG,ANG,20,-20.,20.);

  TH2F *dyvstx3e = new TH2F("dyvstx3e","  ",BIN,-ANG,ANG,20,-20.,20.);
  
  TH2F *dxvsty3e = new TH2F("dxvsty3e","  ",BIN,-ANG,ANG,20,-20.,20.);

  TH1F *hqx1 = new TH1F("hqx1","hqx1",50,0.,100.); 
  TH1F *hqx2 = new TH1F("hqx2","hqx2",50,0.,100.); 
  TH1F *hqy1 = new TH1F("hqy1","hqy1",50,0.,100.); 
  TH1F *hqy2 = new TH1F("hqy2","hqy2",50,0.,100.); 


  //for second file

  vector<TH1F* > ny2;
  vector<TH1F* > nx2;
  vector<TH1F* > dy2;
  vector<TH1F* > dx2;
  vector<TH2F* > dyvsy2;
  vector<TH2F* > dyvsyM2;
  vector<TH2F* > dyvsyP2;
  vector<TH2F* > dyvsx2;
  vector<TH2F* > constdyvsy2;
  vector<TH2F* > constdyvsx2;
  vector<TH2F* > mupulldyvsy2;
  vector<TH2F* > mupulldyvsyty2;
  vector<TH2F* > mupulldyvsx2;
  vector<TH2F* > dxvsy2;
  vector<TH2F* > dxvsx2;
  vector<TH2F* > constdxvsy2;
  vector<TH2F* > constdxvsx2;
  vector<TH2F* > mupulldxvsy2;
  vector<TH2F* > mupulldxvsx2;
  vector<TH2F* > mupulldxvsxtx2;
  vector<TH2F* > tyvsy2;
  vector<TH2F* > tyvsx2;
  vector<TH2F* > txvsy2;
  vector<TH2F* > txvsx2;
  vector<TH2F* > dyvsty2;
  vector<TH2F* > dyvstx2;
  vector<TH2F* > dxvstx2;
  vector<TH2F* > dxvsty2;
  vector<TH1F* > fgdy2;
  vector<TH1F* > fgdx2;
  vector<TH1F* > tyTotDet2;

  for (Int_t dii=0; dii<5; dii++)
    {
      ny2.push_back(new TH1F (("ny2"+detname[dii]+detname[dii+1]).c_str(), ("ny ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 11,0.,10.));
      nx2.push_back(new TH1F (("nx2"+detname[dii]+detname[dii+1]).c_str(), ("nx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 11,0.,10.));
      dy2.push_back(new TH1F (("dy2"+detname[dii]+detname[dii+1]).c_str(), ("dy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 100,-50.,50.));
      dx2.push_back(new TH1F (("dx2"+detname[dii]+detname[dii+1]).c_str(), ("dx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 100,-50.,50.));
  
      dyvsy2.push_back(new TH2F (("dyvsy2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dyvsyM2.push_back(new TH2F (("dyvsyM2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y for -X("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dyvsyP2.push_back(new TH2F (("dyvsyP2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y for +X("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dyvsx2.push_back(new TH2F (("dyvsx2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdyvsy2.push_back(new TH2F (("constdyvsy2"+detname[dii]+detname[dii+1]).c_str(), ("dyvsy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdyvsx2.push_back(new TH2F (("constdyvsx2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldyvsy2.push_back(new TH2F (("mupulldyvsy2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldyvsyty2.push_back(new TH2F (("mupulldyvsyty2"+detname[dii]+detname[dii+1]).c_str(), ("dyty vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-700.,700.,20,-20.,20.));
      mupulldyvsx2.push_back(new TH2F (("mupulldyvsx2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      
      dxvsy2.push_back(new TH2F (("dxvsy2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      dxvsx2.push_back(new TH2F (("dxvsx2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdxvsy2.push_back(new TH2F (("constdxvsy2"+detname[dii]+detname[dii+1]).c_str(), ("dxvsy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      constdxvsx2.push_back(new TH2F (("constdxvsx2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldxvsy2.push_back(new TH2F (("mupulldxvsy2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldxvsx2.push_back(new TH2F (("mupulldxvsx2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));
      mupulldxvsxtx2.push_back(new TH2F (("mupulldxvsxtx2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs xtx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-20.,20.));

      tyvsy2.push_back(new TH2F (("tyvsy2"+detname[dii]+detname[dii+1]).c_str(), ("ty vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 35,-1200.,1200.,20,-1.,1.));
      tyvsx2.push_back(new TH2F (("tyvsx2"+detname[dii]+detname[dii+1]).c_str(), ("ty vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-1.,1.));
      txvsy2.push_back(new TH2F (("txvsy2"+detname[dii]+detname[dii+1]).c_str(), ("tx vs y ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-1.,1.));
      txvsx2.push_back(new TH2F (("txvsx2"+detname[dii]+detname[dii+1]).c_str(), ("tx vs x ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.,20,-1.,1.));

      dyvsty2.push_back(new TH2F (("dyvsty2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs ty ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));
      dxvstx2.push_back(new TH2F (("dxvstx2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs tx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));
      dxvsty2.push_back(new TH2F (("dxvsty2"+detname[dii]+detname[dii+1]).c_str(), ("dx vs ty ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));
      dyvstx2.push_back(new TH2F (("dyvstx2"+detname[dii]+detname[dii+1]).c_str(), ("dy vs tx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), BIN,-ANG,ANG,20,-20.,20.));

      fgdy2.push_back(new TH1F (("fgdy2"+detname[dii]+detname[dii+1]).c_str(), ("fgdy ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.));
      fgdx2.push_back(new TH1F (("fgdx2"+detname[dii]+detname[dii+1]).c_str(), ("fgdx ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 20,-1200.,1200.));

      tyTotDet2.push_back(new TH1F (("tyTotDet"+detname[dii]+detname[dii+1]).c_str(), ("total ty ("+detname[dii]+"-"+detname[dii+1]+")").c_str(), 25,-1.,1.));
    }

  TH1F *nyT2 = new TH1F("nyT2","ny",21, 0, 20);
  
  TH1F *nxT2 = new TH1F("nxT2","nx",21, 0, 20);
  
  TH1F *dyT2 = new TH1F("dyT2","dy",100.,-50.,50.);
  TH1F *SyT2f1 = new TH1F("SyT2f1","Sy FGD1",200,-961.,961.);
  TH1F *SxT2f1 = new TH1F("SxT2f1","Sx FGD1",200,-961.,961.);
  TH1F *SyT2f2 = new TH1F("SyT2f2","Sy FGD2",200,-961.,961.);
  TH1F *SxT2f2 = new TH1F("SxT2f2","Sx FGD2",200,-961.,961.);
  TH1F *SyzT2 = new TH1F("SyzT2","Syz",3600.,0.,3600.);
  TH1F *SxzT2 = new TH1F("SxzT2","Sxz",3600.,0.,3600.);

  TH1F *dxT2 = new TH1F("dxT2","dx",100.,-50.,50.);
  TH1F *dxN2 = new TH1F("dxN2","dx (x<0)",100.,-50.,50.);
  TH1F *dxP2 = new TH1F("dxP2","dx (x>0)",100.,-50.,50.);

  TH1F *dx3eN2 = new TH1F("dx3eN2","dx (TPC3-DSECAL) (x<0)",100.,-100.,100.);
  TH1F *dx3eP2 = new TH1F("dx3eP2","dx (TPC3-DSECAL) (x>0)",100.,-100.,100.);
  
  TH2F *dxvsxT2 = new TH2F("dxvsxT2","dx vs x",10,-900.,900.,20,-20.,20.); 
  TH2F *dxvsyT2 = new TH2F("dxvsyT2","dx vs y",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsxT2 = new TH2F("dyvsxT2","dy vs x",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsyT2 = new TH2F("dyvsyT2","dy vs y",10,-900.,900.,20,-20.,20.); 
  
  TH1F *tyT2 = new TH1F("tyT2","ty",BIN,-ANG,ANG);
  TH1F *txT2 = new TH1F("txT2","tx",BIN,-ANG,ANG);
  TH1F *txN2 = new TH1F("txN2","tx (x<0)",100.,-ANG,ANG);
  TH1F *txP2 = new TH1F("txP2","tx (x>0)",100.,-ANG,ANG);
  TH1F *tyeT2 = new TH1F("tyeT2","ty (DSECAL)",100.,-ANG,ANG);
  TH1F *txeT2 = new TH1F("txeT2","tx (DSECAL)",100.,-ANG,ANG);
  TH1F *txeN2 = new TH1F("txeN2","tx (DSECAL) (x<0)",100.,-ANG,ANG);
  TH1F *txeP2 = new TH1F("txeP2","tx (DSECAL) (x>0)",100.,-ANG,ANG);
  
  TH2F *txvsxT2 = new TH2F("txvsxT2","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *txvsyT2 = new TH2F("txvsyT2","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *tyvsxT2 = new TH2F("tyvsxT2","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *tyvsyT2 = new TH2F("tyvsyT2","  ",10,-900.,900.,20,-20.,20.); 
 
  TH2F *dyvstxT2 = new TH2F("dyvstxT2","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dyvstyT2 = new TH2F("dyvstyT2","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstxT2 = new TH2F("dxvstxT2","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstyT2 = new TH2F("dxvstyT2","  ",20,-ANG,ANG,20,-20.,20.);

  TH2F *dyvstx3e2 = new TH2F("dyvstx3e2","  ",BIN,-ANG,ANG,20,-20.,20.);
  
  TH2F *dxvsty3e2 = new TH2F("dxvsty3e2","  ",BIN,-ANG,ANG,20,-20.,20.);

  cout<<"here11"<<endl;

  cout<<"here1"<<endl;
   
  //open constants file to write out the constants
  ofstream oconstfl;
  oconstfl.open("megatrackalignconst.log", ios::out | ios::app);
  if(!oconstfl){cout<<"failed to open constants file to write"<<endl;}

  // Read in data files
  TChain *chain = new TChain("tpcfgdtree");
  //TChain *chain = new TChain("quickrecotree");

  ifstream infl("RunListB.log");
  if(!infl){cout<<"failed to open run list file"<<endl;}

  vector <std::string> filnamEString;
  std::string filnamEStringtmp;
  Int_t filcount = 0;
  
  while (true)
    {
      infl >> filnamEStringtmp;
      if(infl.eof()) break;
      if (filnamEStringtmp.find("tpcfgdoa") == 0)
       {
          filnamEString.push_back(filnamEStringtmp);
          filcount++;
        }
    }


  for(Int_t fii=0; fii<filcount; fii++)
    {
      chain->AddFile(filnamEString[fii].c_str());
      cout<<"files "<<filnamEString[fii]<<endl;
    }
 
  filnamEString.clear();
  infl.close();

  t2kalignvar *skt2kalignvar = new t2kalignvar;
  chain->SetBranchAddress("st2kalignvar", &skt2kalignvar);
 
  t2kalignvart *skt2kalignvart = new t2kalignvart;
  chain->SetBranchAddress("st2kalignvart", &skt2kalignvart);

  Long64_t nentries = chain->GetEntries();

  std::cout << " ENTRIES : " << nentries << std::endl; 

  //set some cut parameters
  //double chicut = 1., mupullcut = 2.;
  double chicut = 1.5, mupullcut = 2.5;

  Int_t avegtimeruns = 0, avegtimerunsE = 0;

  for (Long64_t iev = 0; iev < nentries; iev++)
    {
    
      chain->GetEntry(iev);

      if(iev == 0) avegtimeruns = skt2kalignvar->Stimestampev/100;
      if(iev == nentries-10) 
  {
    avegtimerunsE = (skt2kalignvar->Stimestampev/100 - avegtimeruns)/2;
    avegtimeruns = (avegtimeruns+skt2kalignvar->Stimestampev/100)/2;
  }
      hndofty->Fill(skt2kalignvar->Sty,skt2kalignvar->Sndoftrack);

      if(skt2kalignvar->Sndoftrack <= 36) continue;
      //if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack >= chicut) continue;
      //if(fabs(skt2kalignvar->Smupull)>mupullcut) continue; //!!!!!CAUTION!!! uncomment out for non horizontal planes
      if(skt2kalignvar->Stpctrkcount[0]> 1 ||
   skt2kalignvar->Stpctrkcount[1]> 1 ||
   skt2kalignvar->Stpctrkcount[2]> 1 ) continue;


      int tpc = TPC(skt2kalignvar->Sztrack);

      for (Int_t hii=0; hii<3; hii++)
  {
    if(skt2kalignvar->Stpcnum[hii] == 1 && skt2kalignvar->Sndoftrack > 0){redchitpctrk[hii]->Fill(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack);}
  }

      int nyhits11 = 0;
      int nyhits12 = 0;
      int nyhits22 = 0;
      int nyhits23 = 0;

 

      for(int i = 0; i < skt2kalignvar->Snyhits; i++ ){ // Loop over y hit

  int fgd =  FGD(skt2kalignvar->Syz[i]);

     //get some values of the track to use in the extrapolation
      Double_t vtxx = skt2kalignvar->Sxtrack;
      Double_t vtxy = skt2kalignvar->Sytrack;
      Double_t vtxz = skt2kalignvar->Sztrack;
      Double_t vtxangy = skt2kalignvar->Sty;
      Double_t vtxangx = skt2kalignvar->Stx;
      if(fabs(skt2kalignvar->Strkenddir.Z()) > 0. && ( fgd == -1 || fgd == -2) )
  {
    if( fgd == -1 || fgd == -2) {
      vtxx = skt2kalignvar->Strkendpos.X(); 
      vtxy = skt2kalignvar->Strkendpos.Y(); 
      vtxz = skt2kalignvar->Strkendpos.Z(); 
      vtxangy = (skt2kalignvar->Strkenddir.Y()/skt2kalignvar->Strkenddir.Z());
      vtxangx = (skt2kalignvar->Strkenddir.X()/skt2kalignvar->Strkenddir.Z());
   } 
  }

  //if(skt2kalignvar->Sqy[i] < 15.) continue;
  //try to calculate the extrapolation by hand
  Double_t newdy = 9999., newdx = 9999.;
  //if(skt2kalignvar->Srho != 0. && skt2kalignvar->Stx != 0)
  Double_t yoffgdpl = 9999.;
  if(skt2kalignvar->Srho != 0.)
    {

      Double_t centZP = 0.,centZM = 0.,centYP = 0.,centYM = 0.;
      centZP = vtxz + fabs(1./(skt2kalignvar->Srho * sqrt(1./(pow(vtxangy , 2))+1.)));
      centZM = vtxz - fabs(1./(skt2kalignvar->Srho * sqrt(1./(pow(vtxangy , 2))+1.)));
      centYP = vtxy + fabs(1./(skt2kalignvar->Srho * sqrt(pow(vtxangy , 2)+1.)));
      centYM = vtxy - fabs(1./(skt2kalignvar->Srho * sqrt(pow(vtxangy , 2)+1.)));

      //pick the proper center-of-circle pair of co-ordinates
      Double_t centZ = 0., centY = 0.;
      if(skt2kalignvar->Strkdir.fY > 0.)
        {
    if(skt2kalignvar->Srho < 0.){centZ = centZM;centY = centYP;}
    else{centZ = centZP;centY = centYM;}
        }
      else
        {
    if(skt2kalignvar->Srho < 0.){centZ = centZP;centY = centYP;}
    else{centZ = centZM;centY = centYM;}
        }


      //by now we have the center of the circle of curvature
      //find y value of the circle (curved from BField) when it intercepts the fgd hit's plane
      Double_t yoffgdplP = sqrt(pow((1./skt2kalignvar->Srho),2) - pow((skt2kalignvar->Syz[i]-centZ),2)) + centY;
      Double_t yoffgdplM = 0. - sqrt(pow((1./skt2kalignvar->Srho),2) - pow((skt2kalignvar->Syz[i]-centZ),2)) + centY;
   
      //take the min distance from the track vertex position
      //caution!! this assumes the track doesn't curve up on itself. I think....
      //double distPsq = pow((yoffgdplP - skt2kalignvar->Sytrack),2) + pow((skt2kalignvar->Syz[i] - skt2kalignvar->Sztrack),2);
      double distMsq = pow((yoffgdplM - vtxy),2) + pow((skt2kalignvar->Syz[i] - vtxz),2);
      double distPsq = pow((yoffgdplP - vtxy),2) + pow((skt2kalignvar->Syz[i] - vtxz),2);
      yoffgdpl = 0.;
      if(distPsq < distMsq){newdy = yoffgdplP - skt2kalignvar->Sy[i]; yoffgdpl = yoffgdplP; newstddytest9->Fill(yoffgdpl,skt2kalignvar->Syy[i]); }
      if(distPsq >= distMsq){newdy = yoffgdplM - skt2kalignvar->Sy[i]; yoffgdpl = yoffgdplM; newstddytest10->Fill(yoffgdpl,skt2kalignvar->Syy[i]);}
      if(fabs(yoffgdpl) > 850.) continue; //remove extrapolations outside the FGD boundary

      newstddy->Fill(newdy,skt2kalignvar->Sdy[i]);
      
      //to test prob
      if(skt2kalignvar->Strkdir.fY > 0.)
        {
    if(skt2kalignvar->Srho > 0.){newstddytest1->Fill(newdy,skt2kalignvar->Sdy[i]);}
    if(skt2kalignvar->Srho < 0.){newstddytest2->Fill(newdy,skt2kalignvar->Sdy[i]);}
    
        }
      else
        {
    if(skt2kalignvar->Srho > 0.){newstddytest3->Fill(newdy,skt2kalignvar->Sdy[i]);}
    else{newstddytest4->Fill(newdy,skt2kalignvar->Sdy[i]);}
  
        }
      if(skt2kalignvar->Strkdir.fZ > 0.)
        {
    if(skt2kalignvar->Srho > 0.){newstddytest5->Fill(newdy,skt2kalignvar->Sdy[i]);}
    else{newstddytest6->Fill(newdy,skt2kalignvar->Sdy[i]);}
  
        }
      else
        {
    if(skt2kalignvar->Srho > 0.){newstddytest7->Fill(newdy,skt2kalignvar->Sdy[i]);}
    else{newstddytest8->Fill(newdy,skt2kalignvar->Sdy[i]);}
  
        }

    }

  //check dx in a simillar way, but there is no BField effect in X
  double intrcpX =  vtxx - vtxangx * vtxz;
  double xoffgdpl = vtxangx * skt2kalignvar->Sxz[i] + intrcpX;
  newdx = xoffgdpl - skt2kalignvar->Sx[i];
  if(fabs(xoffgdpl) > 850.) continue; //remove extrapolations outside the FGD boundary

  if(newdy > 9990.) continue; //default
  if(newdx > 9990.) continue; //default

  Double_t deldy = newdy - skt2kalignvar->Sdy[i];
  Double_t deldx = newdx - skt2kalignvar->Sdx[i];
  //Double_t deldy = yoffgdpl - skt2kalignvar->Syy[i];
  deltady->Fill(deldy);
  deltadx->Fill(deldx);
  deltadyy->Fill(deldy,skt2kalignvar->Syy[i]);
  deltadyx->Fill(deldy,skt2kalignvar->Syx[i]);
  deltadxy->Fill(deldx,skt2kalignvar->Syy[i]);
  deltadxx->Fill(deldx,skt2kalignvar->Syx[i]);

  if(fabs(skt2kalignvar->Sdy[i]) > 50.) continue;
  if(fabs(newdy) > 50.) continue;


  if( fgd == -1 || fgd == 1 || fgd == -2 || fgd == 2)
    {
      hqy1->Fill(skt2kalignvar->Sqy[i]);
    } 
  if(skt2kalignvar->Sqy[i] < 15. || skt2kalignvar->Sqy[i] > 35.) continue;

  if (skt2kalignvar->Snyhits>0) nyT->Fill(skt2kalignvar->Snyhits);
  if (skt2kalignvar->Snxhits>0) nxT->Fill(skt2kalignvar->Snxhits);

  tyT->Fill(vtxangy);

  dyT->Fill(newdy); //skt2kalignvar->Sdy[i]
  if(skt2kalignvar->Syz[i] > 100 && skt2kalignvar->Syz[i] < 450) SyTf1->Fill(skt2kalignvar->Sy[i]);
  if(skt2kalignvar->Sxz[i] > 100 && skt2kalignvar->Sxz[i] < 450) SxTf1->Fill(skt2kalignvar->Sx[i]);
  if(skt2kalignvar->Syz[i] > 1450 && skt2kalignvar->Syz[i] < 1900) SyTf2->Fill(skt2kalignvar->Sy[i]);
  if(skt2kalignvar->Sxz[i] > 1450 && skt2kalignvar->Sxz[i] < 1900) SxTf2->Fill(skt2kalignvar->Sx[i]);
  SxzT->Fill(skt2kalignvar->Sxz[i]);
  SyzT->Fill(skt2kalignvar->Syz[i]);
  dyvsxT->Fill(skt2kalignvar->Syx[i],newdy);
  dyvsyT->Fill(skt2kalignvar->Sy[i],newdy);

  dyvstyT->Fill(vtxangy,newdy);
  dyvstxT->Fill(vtxangx,newdy);
  
  tyvsxT->Fill(skt2kalignvar->Syx[i],vtxangy);
  tyvsyT->Fill(skt2kalignvar->Sy[i],vtxangy);



  
  if( fgd == -1 && tpc == 1) {
    nyhits11++;
    dy[0]->Fill(newdy);
    tyTotDet[0]->Fill(vtxangy);
    dyvsy[0]->Fill(skt2kalignvar->Sy[i],newdy);
    fgdy[0]->Fill(skt2kalignvar->Sy[i]);
    fgdx[0]->Fill(skt2kalignvar->Sx[i]);
    tyvsyA[0]->Fill(skt2kalignvar->Syy[i],vtxangy);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP[0]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM[0]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx[0]->Fill(skt2kalignvar->Syx[i],newdy);
    dyvsty[0]->Fill(vtxangy,newdy);
    dyvstx[0]->Fill(vtxangx,newdy);



    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy[0]->Fill(skt2kalignvar->Sy[i],newdy);
      constdyvsx[0]->Fill(skt2kalignvar->Syx[i],newdy);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          if (vtxangy<1. && vtxangy>-1.)
          mupulldyvsy[0]->Fill(skt2kalignvar->Syy[i],newdy);
          //if (vtxangy<1. && vtxangy>-1. )
          //if (fabs(vtxangy) < 4.)

          mupulldyvsx[0]->Fill(skt2kalignvar->Syx[i],newdy);
         
          if (fabs(newdy) < 10.)mupulldyvsyty[0]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          tyvsy[0]->Fill(skt2kalignvar->Syy[i],vtxangy);
          if(fabs(skt2kalignvart->STruePrimMom.Z() > 0.))
       tyvsyPrim[0]->Fill(skt2kalignvart->STruePrimPos.Y(),skt2kalignvart->STruePrimMom.Y()/skt2kalignvart->STruePrimMom.Z());
        }
    }
      }
  }
  if( fgd == 1 && tpc == 2 ) {
    nyhits12++;
    dy[1]->Fill(newdy);
    tyTotDet[1]->Fill(vtxangy);
    dyvsy[1]->Fill(skt2kalignvar->Sy[i],newdy);
    fgdy[1]->Fill(skt2kalignvar->Sy[i]);
    fgdx[1]->Fill(skt2kalignvar->Sx[i]);
    tyvsyA[1]->Fill(skt2kalignvar->Syy[i],vtxangy);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP[1]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM[1]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx[1]->Fill(skt2kalignvar->Syx[i],newdy);
    dyvsty[1]->Fill(vtxangy,newdy);
    dyvstx[1]->Fill(vtxangx,newdy);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy[1]->Fill(skt2kalignvar->Sy[i],newdy);
      constdyvsx[1]->Fill(skt2kalignvar->Syx[i],newdy);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldyvsy[1]->Fill(skt2kalignvar->Syy[i],newdy);
          if (fabs(newdy) < 10.) mupulldyvsyty[1]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          mupulldyvsx[1]->Fill(skt2kalignvar->Syx[i],newdy);
          tyvsy[1]->Fill(skt2kalignvar->Syy[i],vtxangy);//skt2kalignvar->Sty);
          if(fabs(skt2kalignvart->STruePrimMom.Z() > 0.))
       tyvsyPrim[1]->Fill(skt2kalignvart->STruePrimPos.Y(),skt2kalignvart->STruePrimMom.Y()/skt2kalignvart->STruePrimMom.Z());
        }
    }
      }
  }
  if( fgd == -2 && tpc == 2) {
    nyhits22++;
    dy[2]->Fill(newdy);
    tyTotDet[2]->Fill(vtxangy);
    dyvsy[2]->Fill(skt2kalignvar->Sy[i],newdy);
    fgdy[2]->Fill(skt2kalignvar->Sy[i]);
    fgdx[2]->Fill(skt2kalignvar->Sx[i]);
    tyvsyA[2]->Fill(skt2kalignvar->Syy[i],vtxangy);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP[2]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM[2]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx[2]->Fill(skt2kalignvar->Syx[i],newdy);
    dyvsty[2]->Fill(vtxangy,newdy);
    dyvstx[2]->Fill(vtxangx,newdy);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy[2]->Fill(skt2kalignvar->Sy[i],newdy);
      constdyvsx[2]->Fill(skt2kalignvar->Syx[i],newdy);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldyvsy[2]->Fill(skt2kalignvar->Syy[i],newdy);
          if (fabs(newdy) < 10.) mupulldyvsyty[2]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          mupulldyvsx[2]->Fill(skt2kalignvar->Syx[i],newdy);
          tyvsy[2]->Fill(skt2kalignvar->Syy[i],vtxangy);//skt2kalignvar->Sty);
          if(fabs(skt2kalignvart->STruePrimMom.Z() > 0.))
      tyvsyPrim[2]->Fill(skt2kalignvart->STruePrimPos.Y(),skt2kalignvart->STruePrimMom.Y()/skt2kalignvart->STruePrimMom.Z());
        }
    }
      }
  }
  if( fgd == 2 && tpc == 3) {
    nyhits23++;
    dy[3]->Fill(newdy);
    tyTotDet[3]->Fill(vtxangy);
    dyvsy[3]->Fill(skt2kalignvar->Sy[i],newdy);
    fgdy[3]->Fill(skt2kalignvar->Sy[i]);
    fgdx[3]->Fill(skt2kalignvar->Sx[i]);
    tyvsyA[3]->Fill(skt2kalignvar->Syy[i],vtxangy);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP[3]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM[3]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx[3]->Fill(skt2kalignvar->Syx[i],newdy);
    dyvsty[3]->Fill(vtxangy,newdy);
    dyvstx[3]->Fill(vtxangx,newdy);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy[3]->Fill(skt2kalignvar->Sy[i],newdy);
      constdyvsx[3]->Fill(skt2kalignvar->Syx[i],newdy);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldyvsy[3]->Fill(skt2kalignvar->Syy[i],newdy);
          if (fabs(newdy) < 10.) mupulldyvsyty[3]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          mupulldyvsx[3]->Fill(skt2kalignvar->Syx[i],newdy);
          tyvsy[3]->Fill(skt2kalignvar->Syy[i],vtxangy);//skt2kalignvar->Sty);
          if(fabs(skt2kalignvart->STruePrimMom.Z() > 0.))
      tyvsyPrim[3]->Fill(skt2kalignvart->STruePrimPos.Y(),skt2kalignvart->STruePrimMom.Y()/skt2kalignvart->STruePrimMom.Z());
        }
    }
      }
  }

      }
      if (nyhits11>0) ny[0]->Fill(nyhits11);
      if (nyhits12>0) ny[1]->Fill(nyhits12);
      if (nyhits22>0) ny[2]->Fill(nyhits22);
      if (nyhits23>0) ny[3]->Fill(nyhits23);
    

      int nxhits11 = 0;
      int nxhits12 = 0;
      int nxhits22 = 0;
      int nxhits23 = 0;

      for(int i = 0; i < skt2kalignvar->Snxhits; i++ ){  // loop over x hits
 
   int fgd =  FGD(skt2kalignvar->Sxz[i]);
  
  //get some values of the track to use in the extrapolation
  Double_t vtxx = skt2kalignvar->Sxtrack;
  Double_t vtxy = skt2kalignvar->Sytrack;
  Double_t vtxz = skt2kalignvar->Sztrack;
  Double_t vtxangy = skt2kalignvar->Sty;
  Double_t vtxangx = skt2kalignvar->Stx;
  if(fabs(skt2kalignvar->Strkenddir.Z()) > 0. && ( fgd == -1 || fgd == -2) )
    {
      if( fgd == -1 || fgd == -2) {
        vtxx = skt2kalignvar->Strkendpos.X(); 
        vtxy = skt2kalignvar->Strkendpos.Y(); 
        vtxz = skt2kalignvar->Strkendpos.Z(); 
        vtxangy = (skt2kalignvar->Strkenddir.Y()/skt2kalignvar->Strkenddir.Z());
        vtxangx = (skt2kalignvar->Strkenddir.X()/skt2kalignvar->Strkenddir.Z());
      }
    }
  
  Double_t newdx = 9999.;
  //check dx in a simillar way, but there is no BField effect in X
  //so like the B=0 case
  double intrcpX =  vtxx - vtxangx * vtxz;
  double xoffgdpl = vtxangx * skt2kalignvar->Sxz[i] + intrcpX;
  newdx = xoffgdpl - skt2kalignvar->Sx[i];
  if(newdx > 9990.) continue; //default
  if(fabs(xoffgdpl) > 850.) continue; //remove extrapolations outside the FGD boundary

  newstddx->Fill(newdx,skt2kalignvar->Sdx[i]);


  if( fgd == -1 || fgd == 1 || fgd == -2 || fgd == 2)
    {
      hqx1->Fill(skt2kalignvar->Sqx[i]);
    } 

  if(skt2kalignvar->Sqx[i] < 15.) continue;

  hmom->Fill(skt2kalignvar->Smomentum);

  txT->Fill(vtxangx);
      
  dxT->Fill(newdx);
  if( skt2kalignvar->Sx[i] < 0 ) dxN->Fill(newdx);
  if( skt2kalignvar->Sx[i] > 0 ) dxP->Fill(newdx);
  dxvsxT->Fill(skt2kalignvar->Sx[i],newdx);
  dxvsyT->Fill(skt2kalignvar->Sxy[i],newdx);

  dxvstyT->Fill(vtxangy,newdx);
  dxvstxT->Fill(vtxangx,newdx);



  if( fgd == -1 && tpc == 1 ) {
    nxhits11++;
    dx[0]->Fill(newdx);
    dxvsx[0]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy[0]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx[0]->Fill(vtxangx,newdx);
    dxvsty[0]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx[0]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy[0]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx[0]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy[0]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx[0]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy[0]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }
  if( fgd == 1 && tpc == 2 ) {
    nxhits12++;
    dx[1]->Fill(newdx);
    dxvsx[1]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy[1]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx[1]->Fill(vtxangx,newdx);
    dxvsty[1]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx[1]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy[1]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx[1]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy[1]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx[1]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy[1]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }
  if( fgd == -2 && tpc == 2 ) {
    nxhits22++;
    dx[2]->Fill(newdx);
    dxvsx[2]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy[2]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx[2]->Fill(vtxangx,newdx);
    dxvsty[2]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx[2]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy[2]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx[2]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy[2]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx[2]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy[2]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }
  if( fgd == 2 && tpc == 3 ) {
    nxhits23++;
    dx[3]->Fill(newdx);
    dxvsx[3]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy[3]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx[3]->Fill(vtxangx,newdx);
    dxvsty[3]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx[3]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy[3]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx[3]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy[3]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx[3]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy[3]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }

      }
      if (nxhits11>0) nx[0]->Fill(nxhits11);
      if (nxhits12>0) nx[1]->Fill(nxhits12);
      if (nxhits22>0) nx[2]->Fill(nxhits22);
      if (nxhits23>0) nx[3]->Fill(nxhits23);

      // DSECAL
      //if (neyhits>0) ny[4]->Fill(neyhits);
      //if (nexhits>0) nx[4]->Fill(nexhits);

      // Loop over DSECAL y hits to get residuals for the TPC track
      //if (neyhits>0) tyeT->Fill(ty);
      int selEyhits = 0;
      for(int i = 0; i < skt2kalignvar->Sneyhits; i++ ){
  if (skt2kalignvar->SneyLayer[i] == 1) { // select first y layer
    double charge = skt2kalignvar->Seqy[i]; 
    if (charge>3.0 && charge<10.0) { // charge cut to remove noise hits
      dy[4]->Fill(skt2kalignvar->Sedy[i]);
      dyvsy[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sedy[i]);
      //fill a dyvsy for X +/-
      //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
      if(skt2kalignvar->Sex[i]+skt2kalignvar->Sedx[i] >= 0.)
        {
    dyvsyP[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sedy[i]);
        }
      else
        {
    dyvsyM[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sedy[i]);
        }
      dyvsx[4]->Fill(skt2kalignvar->Seyx[i],skt2kalignvar->Sedy[i]);
    
      dyvsty[4]->Fill(skt2kalignvar->Sty,skt2kalignvar->Sedy[i]);
      dyvstx3e->Fill(skt2kalignvar->Stx,skt2kalignvar->Sedy[i]);
    
      tyvsx[4]->Fill(skt2kalignvar->Seyx[i],skt2kalignvar->Sty);
      tyvsy[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sty);
      selEyhits++;
    }
  }
      }
      if (selEyhits>0) tyeT->Fill(skt2kalignvar->Sty);
      if (selEyhits>0) ny[4]->Fill(selEyhits);


      // Loop over DSECAL x hits to get residuals for the TPC track
      //if (nexhits>0) txeT->Fill(tx);
      int selExhits = 0;
      for(int i = 0; i < skt2kalignvar->Snexhits; i++ ){
  if (skt2kalignvar->SnexLayer[i] == 1) { // select first x layer
    double charge = skt2kalignvar->Seqx[i]; 
    if (charge>3.0 && charge<10.0) { // charge cut to remove noise hits
      dx[4]->Fill(skt2kalignvar->Sedx[i]);
      if( skt2kalignvar->Sex[i] < 0 ) dx3eN->Fill(skt2kalignvar->Sedx[i]);
      if( skt2kalignvar->Sex[i] > 0 ) dx3eP->Fill(skt2kalignvar->Sedx[i]);
      dxvsx[4]->Fill(skt2kalignvar->Sex[i],skt2kalignvar->Sedx[i]);
      dxvsy[4]->Fill(skt2kalignvar->Sexy[i],skt2kalignvar->Sedx[i]);
    
      dxvstx[4]->Fill(skt2kalignvar->Stx,skt2kalignvar->Sedx[i]);
      dxvsty3e->Fill(skt2kalignvar->Sty,skt2kalignvar->Sedx[i]);
    
      txvsx[4]->Fill(skt2kalignvar->Sex[i],skt2kalignvar->Stx);
      txvsy[4]->Fill(skt2kalignvar->Sexy[i],skt2kalignvar->Sty);
      selExhits++;
    }
  }
      }
      if (selExhits>0) txeT->Fill(skt2kalignvar->Stx);
      if (selExhits>0) nx[4]->Fill(selExhits);

    }
  
  chain->Delete();
  //inFile.close();

  //open new file to compare

  //TChain *chain2 = new TChain("tpcfgdtree");
  TChain *chain2 = new TChain("quickrecotree");

  ifstream infl2("RunListNoB.log");
  if(!infl2){cout<<"failed to open run list file"<<endl;}

  filcount = 0;
  
  while (true)
    {
      infl2 >> filnamEStringtmp;
      if(infl2.eof()) break;
      if (filnamEStringtmp.find("tpcfgdoa") == 0)
       {
          filnamEString.push_back(filnamEStringtmp);
          filcount++;
        }
    }



  for(Int_t fii=0; fii<filcount; fii++)
    {
      chain2->AddFile(filnamEString[fii].c_str());
      cout<<"files "<<filnamEString[fii]<<endl;
    }
 
  filnamEString.clear();
  infl2.close();

  chain2->SetBranchAddress("st2kalignvar", &skt2kalignvar);
  chain2->SetBranchAddress("st2kalignvart", &skt2kalignvart);

  Long64_t nentries2 = chain2->GetEntries();
  std::cout << " ENTRIES : " << nentries2 << std::endl; 

  for (Long64_t iev = 0; iev < nentries2; iev++)
    {
    
      chain2->GetEntry(iev);

      //if( TMath::Abs(skt2kalignvar->Stime) < 0.01 ) continue;

      if(skt2kalignvar->Sndoftrack <= 36) continue;
      //if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack >= chicut) continue;
      if(fabs(skt2kalignvar->Smupull)>mupullcut) continue;
      if(skt2kalignvar->Stpctrkcount[0]> 1 ||
   skt2kalignvar->Stpctrkcount[1]> 1 ||
   skt2kalignvar->Stpctrkcount[2]> 1 ) continue;
 //    if(skt2kalignvar->Stpctrkcount[0]< 1 ||
// 	 skt2kalignvar->Stpctrkcount[1]< 1 ||
// 	 skt2kalignvar->Stpctrkcount[2]< 1 ) continue;
//       if(fabs(skt2kalignvar->Stx) > 0.25) continue;


      int tpc = TPC(skt2kalignvar->Sztrack);

      for (Int_t hii=0; hii<3; hii++)
  {
    if(skt2kalignvar->Stpcnum[hii] == 1 && skt2kalignvar->Sndoftrack > 0){redchitpctrk2[hii]->Fill(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack);}
  }

      int nyhits11 = 0;
      int nyhits12 = 0;
      int nyhits22 = 0;
      int nyhits23 = 0;

      for(int i = 0; i < skt2kalignvar->Snyhits; i++ ){ // Loop over y hit
  //if(skt2kalignvar->Sqy[i] < 15. ) continue;

  int fgd =  FGD(skt2kalignvar->Syz[i]);

  //get some values of the track to use in the extrapolation
  Double_t vtxx = skt2kalignvar->Sxtrack;
  Double_t vtxy = skt2kalignvar->Sytrack;
  Double_t vtxz = skt2kalignvar->Sztrack;
  Double_t vtxangy = skt2kalignvar->Sty;
  Double_t vtxangx = skt2kalignvar->Stx;
  if(fabs(skt2kalignvar->Strkenddir.Z()) > 0. && ( fgd == -1 || fgd == -2) )
    {
      if( fgd == -1 || fgd == -2) {
        vtxx = skt2kalignvar->Strkendpos.X(); 
        vtxy = skt2kalignvar->Strkendpos.Y(); 
        vtxz = skt2kalignvar->Strkendpos.Z(); 
        vtxangy = (skt2kalignvar->Strkenddir.Y()/skt2kalignvar->Strkenddir.Z());
        vtxangx = (skt2kalignvar->Strkenddir.X()/skt2kalignvar->Strkenddir.Z());
      }
    }

  //try to calculate the extrapolation by hand
  Double_t newdy = 9999., newdx = 9999.;

  //considering B=0;
  //double intrcp =  skt2kalignvar->Sytrack - skt2kalignvar->Sty * skt2kalignvar->Sztrack;
  //double yoffgdpl = skt2kalignvar->Sty * skt2kalignvar->Syz[i] + intrcp;
  double intrcp = 99999.,yoffgdpl = 99999.;
   if(fabs(skt2kalignvar->Strkenddir.Z()) > 0.)
    {
      if( fgd == -1 || fgd == -2) {
        intrcp =  skt2kalignvar->Strkendpos.Y() - (skt2kalignvar->Strkenddir.Y()/skt2kalignvar->Strkenddir.Z()) * skt2kalignvar->Strkendpos.Z();
        yoffgdpl = (skt2kalignvar->Strkenddir.Y()/skt2kalignvar->Strkenddir.Z()) * skt2kalignvar->Syz[i] + intrcp;
      }
      else
        {
    intrcp =  skt2kalignvar->Sytrack - skt2kalignvar->Sty * skt2kalignvar->Sztrack;
    yoffgdpl = skt2kalignvar->Sty * skt2kalignvar->Syz[i] + intrcp;
        }
    }
  if(fabs(yoffgdpl) > 850.) continue; //remove extrapolations outside the FGD boundary
  newdy = yoffgdpl - skt2kalignvar->Sy[i];

  newstddy2->Fill(newdy,skt2kalignvar->Sdy[i]);

  double intrcpX = 99999.,xoffgdpl = 99999.;
   if(fabs(skt2kalignvar->Strkenddir.Z()) > 0.)
    {
      if( fgd == -1 || fgd == -2) {
        intrcpX =  skt2kalignvar->Strkendpos.X() - (skt2kalignvar->Strkenddir.X()/skt2kalignvar->Strkenddir.Z()) * skt2kalignvar->Strkendpos.Z();
        xoffgdpl = (skt2kalignvar->Strkenddir.X()/skt2kalignvar->Strkenddir.Z()) * skt2kalignvar->Sxz[i] + intrcpX;
      }
      else
        {
    intrcpX =  skt2kalignvar->Sxtrack - skt2kalignvar->Stx * skt2kalignvar->Sztrack;
    xoffgdpl = skt2kalignvar->Stx * skt2kalignvar->Sxz[i] + intrcpX;
        }
    }
  
  if(fabs(xoffgdpl) > 850.) continue; //remove extrapolations outside the FGD boundary
  newdx = xoffgdpl - skt2kalignvar->Sx[i];



  if(newdy > 9990.) continue; //default
  if(newdx > 9990.) continue; //default

  Double_t deldy = newdy - skt2kalignvar->Sdy[i];
  Double_t deldx = newdx - skt2kalignvar->Sdx[i];
  //Double_t deldy = yoffgdpl - skt2kalignvar->Syy[i];
  deltady2->Fill(deldy);

  deltadx2->Fill(deldx);
  //if (newdy > 100.) continue;
  if(fabs(newdy) > 50.) continue;



    if( fgd == -1 || fgd == 1 || fgd == -2 || fgd == 2)
    {
      hqy2->Fill(skt2kalignvar->Sqy[i]);
    } 
  
  if(skt2kalignvar->Sqy[i] < 15. || skt2kalignvar->Sqy[i] > 35.) continue;

  if (skt2kalignvar->Snyhits>0) nyT2->Fill(skt2kalignvar->Snyhits);
  if (skt2kalignvar->Snxhits>0) nxT2->Fill(skt2kalignvar->Snxhits);
     
  tyT2->Fill(skt2kalignvar->Sty);

  dyT2->Fill(newdy);
  if(skt2kalignvar->Syz[i] > 100 && skt2kalignvar->Syz[i] < 450) SyT2f1->Fill(skt2kalignvar->Sy[i]);
  if(skt2kalignvar->Sxz[i] > 100 && skt2kalignvar->Sxz[i] < 450) SxT2f1->Fill(skt2kalignvar->Sx[i]);
  if(skt2kalignvar->Syz[i] > 1450 && skt2kalignvar->Syz[i] < 1900) SyT2f2->Fill(skt2kalignvar->Sy[i]);
  if(skt2kalignvar->Sxz[i] > 1450 && skt2kalignvar->Sxz[i] < 1900) SxT2f2->Fill(skt2kalignvar->Sx[i]);
  SxzT2->Fill(skt2kalignvar->Sxz[i]);
  SyzT2->Fill(skt2kalignvar->Syz[i]);
  dyvsxT2->Fill(skt2kalignvar->Syx[i],newdy);
  dyvsyT2->Fill(skt2kalignvar->Sy[i],newdy);

  //dyvstyT2->Fill(skt2kalignvar->Sty,newdy);
  dyvstyT2->Fill(skt2kalignvar->Sty,newdy);
  dyvstxT2->Fill(skt2kalignvar->Stx,newdy);//newdy);
  
  tyvsxT2->Fill(skt2kalignvar->Syx[i],skt2kalignvar->Sty);
  tyvsyT2->Fill(skt2kalignvar->Sy[i],skt2kalignvar->Sty);


  
  if( fgd == -1 && tpc == 1) {
    nyhits11++;
    dy2[0]->Fill(newdy);
    tyTotDet2[0]->Fill(vtxangy);
    dyvsy2[0]->Fill(skt2kalignvar->Sy[i],newdy);//,skt2kalignvar->Sdy[i]);

    fgdy2[0]->Fill(skt2kalignvar->Sy[i]);
    fgdx2[0]->Fill(skt2kalignvar->Sx[i]);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP2[0]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM2[0]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx2[0]->Fill(skt2kalignvar->Syx[i],newdy);//skt2kalignvar->Sdy[i]);
    dyvsty2[0]->Fill(skt2kalignvar->Sty,newdy);//skt2kalignvar->Sdy[i]);
    dyvstx2[0]->Fill(vtxangx,newdy);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy2[0]->Fill(skt2kalignvar->Sy[i],newdy);//,skt2kalignvar->Sdy[i]);
      constdyvsx2[0]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldyvsy2[0]->Fill(skt2kalignvar->Syy[i],newdy);//,skt2kalignvar->Sdy[i]);
          //if (fabs(skt2kalignvar->Syy[i]*vtxangy) < 1.)
          //if (vtxangy != 0. && fabs(vtxangy)<3.)
          //if(fabs(skt2kalignvar->Syy[i])<600. && fabs(vtxangy)<1.5)
          if (fabs(newdy) < 10.) mupulldyvsyty2[0]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          mupulldyvsx2[0]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
          tyvsy2[0]->Fill(skt2kalignvar->Syy[i],vtxangy);//skt2kalignvar->Sty);
        }
    }
      }
  }
  if( fgd == 1 && tpc == 2 ) {
    nyhits12++;
    dy2[1]->Fill(newdy);
    tyTotDet2[1]->Fill(vtxangy);
    dyvsy2[1]->Fill(skt2kalignvar->Sy[i],newdy);//skt2kalignvar->Sdy[i]);
    fgdy2[1]->Fill(skt2kalignvar->Sy[i]);
    fgdx2[1]->Fill(skt2kalignvar->Sx[i]);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP2[1]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM2[1]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx2[1]->Fill(skt2kalignvar->Syx[i],newdy);//skt2kalignvar->Sdy[i]);
    dyvsty2[1]->Fill(skt2kalignvar->Sty,newdy);//skt2kalignvar->Sdy[i]);
    dyvstx2[1]->Fill(vtxangx,newdy);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy2[1]->Fill(skt2kalignvar->Sy[i],newdy);//,skt2kalignvar->Sdy[i]);
      constdyvsx2[1]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldyvsy2[1]->Fill(skt2kalignvar->Syy[i],newdy);//,skt2kalignvar->Sdy[i]);
          if (fabs(newdy) < 10.) mupulldyvsyty2[1]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          mupulldyvsx2[1]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
          tyvsy2[1]->Fill(skt2kalignvar->Syy[i],vtxangy);//skt2kalignvar->Sty);
        }
    }
      }
  }
  if( fgd == -2 && tpc == 2) {
    nyhits22++;
    dy2[2]->Fill(newdy);
    tyTotDet2[2]->Fill(vtxangy);
    dyvsy2[2]->Fill(skt2kalignvar->Sy[i],newdy);//,skt2kalignvar->Sdy[i]);
    fgdy2[2]->Fill(skt2kalignvar->Sy[i]);
    fgdx2[2]->Fill(skt2kalignvar->Sx[i]);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP2[2]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM2[2]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx2[2]->Fill(skt2kalignvar->Syx[i],newdy);//skt2kalignvar->Sdy[i]);
    dyvsty2[2]->Fill(skt2kalignvar->Sty,newdy);//skt2kalignvar->Sdy[i]);
    dyvstx2[2]->Fill(vtxangx,newdy);
    cout<<"events to look "<<skt2kalignvar->Snevt<<"\t"<<i<<"\t"<<skt2kalignvar->Sdy[i]<<"\t"<<skt2kalignvar->Stpctrkcount[0]<<"\t"<<skt2kalignvar->Stpctrkcount[1]<<"\t"<<skt2kalignvar->Stpctrkcount[2]<<endl;
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy2[2]->Fill(skt2kalignvar->Sy[i],newdy);//,skt2kalignvar->Sdy[i]);
      constdyvsx2[2]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldyvsy2[2]->Fill(skt2kalignvar->Syy[i],newdy);//,skt2kalignvar->Sdy[i]);
          if (fabs(newdy) < 10.) mupulldyvsyty2[2]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          mupulldyvsx2[2]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
          tyvsy2[2]->Fill(skt2kalignvar->Syy[i],vtxangy);//skt2kalignvar->Sty);
        }
    }
      }
  }
  if( fgd == 2 && tpc == 3) {
    nyhits23++;
    dy2[3]->Fill(newdy);
    tyTotDet2[3]->Fill(vtxangy);
    dyvsy2[3]->Fill(skt2kalignvar->Sy[i],newdy);//,skt2kalignvar->Sdy[i]);
    fgdy2[3]->Fill(skt2kalignvar->Sy[i]);
    fgdx2[3]->Fill(skt2kalignvar->Sx[i]);
    //fill a dyvsy for X +/-
    //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
    if(skt2kalignvar->Sx[i]+newdx >= 0.)
      {
        dyvsyP2[3]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    else
      {
        dyvsyM2[3]->Fill(skt2kalignvar->Sy[i],newdy);
      }
    dyvsx2[3]->Fill(skt2kalignvar->Syx[i],newdy);//skt2kalignvar->Sdy[i]);
    dyvsty2[3]->Fill(skt2kalignvar->Sty,newdy);//skt2kalignvar->Sdy[i]);
    dyvstx2[3]->Fill(vtxangx,newdy);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdyvsy2[3]->Fill(skt2kalignvar->Sy[i],newdy);//,newdy);
      constdyvsx2[3]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldyvsy2[3]->Fill(skt2kalignvar->Syy[i],newdy);//,skt2kalignvar->Sdy[i]);
          if (fabs(newdy) < 10.) mupulldyvsyty2[3]->Fill(skt2kalignvar->Syy[i]*vtxangy,newdy);
          mupulldyvsx2[3]->Fill(skt2kalignvar->Syx[i],newdy);//,skt2kalignvar->Sdy[i]);
          tyvsy2[3]->Fill(skt2kalignvar->Syy[i],vtxangy);//skt2kalignvar->Sty);
        }
    }
      }
  }

      }
      if (nyhits11>0) ny2[0]->Fill(nyhits11);
      if (nyhits12>0) ny2[1]->Fill(nyhits12);
      if (nyhits22>0) ny2[2]->Fill(nyhits22);
      if (nyhits23>0) ny2[3]->Fill(nyhits23);
    

      int nxhits11 = 0;
      int nxhits12 = 0;
      int nxhits22 = 0;
      int nxhits23 = 0;

      for(int i = 0; i < skt2kalignvar->Snxhits; i++ ){  // loop over x hits
  Double_t newdx = 9999.;
  //check dx in a simillar way, but there is no BField effect in X
  double intrcpX =  skt2kalignvar->Sxtrack - skt2kalignvar->Stx * skt2kalignvar->Sztrack;
  double xoffgdpl = skt2kalignvar->Stx * skt2kalignvar->Sxz[i] + intrcpX;
  newdx = xoffgdpl - skt2kalignvar->Sx[i];
  if(newdx > 9990.) continue; //default
  if(fabs(xoffgdpl) > 850.) continue; //remove extrapolations outside the FGD boundary

  newstddx2->Fill(newdx,skt2kalignvar->Sdx[i]);

  int fgd =  FGD(skt2kalignvar->Sxz[i]);

  //get some values of the track to use in the extrapolation
  Double_t vtxx = skt2kalignvar->Sxtrack;
  Double_t vtxy = skt2kalignvar->Sytrack;
  Double_t vtxz = skt2kalignvar->Sztrack;
  Double_t vtxangy = skt2kalignvar->Sty;
  Double_t vtxangx = skt2kalignvar->Stx;
  if(fabs(skt2kalignvar->Strkenddir.Z()) > 0. && ( fgd == -1 || fgd == -2) )
    {
      if( fgd == -1 || fgd == -2) {
        vtxx = skt2kalignvar->Strkendpos.X(); 
        vtxy = skt2kalignvar->Strkendpos.Y(); 
        vtxz = skt2kalignvar->Strkendpos.Z(); 
        vtxangy = (skt2kalignvar->Strkenddir.Y()/skt2kalignvar->Strkenddir.Z());
        vtxangx = (skt2kalignvar->Strkenddir.X()/skt2kalignvar->Strkenddir.Z());
      }
    }

  if( fgd == -1 || fgd == 1 || fgd == -2 || fgd == 2)
    {
      hqx2->Fill(skt2kalignvar->Sqx[i]);
    } 

  if(skt2kalignvar->Sqx[i] < 15.) continue;

  hmom2->Fill(skt2kalignvar->Smomentum);


  txT2->Fill(skt2kalignvar->Stx);
      
  dxT2->Fill(newdx);
  if( skt2kalignvar->Sx[i] < 0 ) dxN2->Fill(newdx);
  if( skt2kalignvar->Sx[i] > 0 ) dxP2->Fill(newdx);
  dxvsxT2->Fill(skt2kalignvar->Sx[i],newdx);
  dxvsyT2->Fill(skt2kalignvar->Sxy[i],newdx);

  dxvstyT2->Fill(skt2kalignvar->Sty,newdx);
  dxvstxT2->Fill(skt2kalignvar->Stx,newdx);



  if( fgd == -1 && tpc == 1 ) {
    nxhits11++;
    dx2[0]->Fill(newdx);
    dxvsx2[0]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy2[0]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx2[0]->Fill(skt2kalignvar->Stx,newdx);
    dxvsty2[0]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx2[0]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy2[0]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx2[0]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy2[0]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx2[0]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy2[0]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }
  if( fgd == 1 && tpc == 2 ) {
    nxhits12++;
    dx2[1]->Fill(newdx);
    dxvsx2[1]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy2[1]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx2[1]->Fill(skt2kalignvar->Stx,newdx);
    dxvsty2[1]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx2[1]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy2[1]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx2[1]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy2[1]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx2[1]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy2[1]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }
  if( fgd == -2 && tpc == 2 ) {
    nxhits22++;
    dx2[2]->Fill(newdx);
    dxvsx2[2]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy2[2]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx2[2]->Fill(skt2kalignvar->Stx,newdx);
    dxvsty2[2]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx2[2]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy2[2]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx2[2]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy2[2]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx2[2]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy2[2]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }
  if( fgd == 2 && tpc == 3 ) {
    nxhits23++;
    dx2[3]->Fill(newdx);
    dxvsx2[3]->Fill(skt2kalignvar->Sx[i],newdx);
    dxvsy2[3]->Fill(skt2kalignvar->Sxy[i],newdx);
    dxvstx2[3]->Fill(skt2kalignvar->Stx,newdx);
    dxvsty2[3]->Fill(vtxangy,newdx);
    if(skt2kalignvar->Sndoftrack > 0)
      {
        if(skt2kalignvar->Schioftrack/skt2kalignvar->Sndoftrack < chicut)
    {
      constdxvsx2[3]->Fill(skt2kalignvar->Sx[i],newdx);
      constdxvsy2[3]->Fill(skt2kalignvar->Syx[i],newdx);
      if(fabs(skt2kalignvar->Smupull)<mupullcut)
        {
          mupulldxvsx2[3]->Fill(skt2kalignvar->Sx[i],newdx);
          mupulldxvsy2[3]->Fill(skt2kalignvar->Sxy[i],newdx);
          if (fabs(newdx) < 10.) mupulldxvsxtx2[3]->Fill(skt2kalignvar->Sxx[i]*vtxangx,newdx);
          txvsy2[3]->Fill(skt2kalignvar->Sxy[i],skt2kalignvar->Stx);
        }
    }
      }
  }

      }
      if (nxhits11>0) nx2[0]->Fill(nxhits11);
      if (nxhits12>0) nx2[1]->Fill(nxhits12);
      if (nxhits22>0) nx2[2]->Fill(nxhits22);
      if (nxhits23>0) nx2[3]->Fill(nxhits23);

      // DSECAL
      //if (neyhits>0) ny2[4]->Fill(neyhits);
      //if (nexhits>0) nx2[4]->Fill(nexhits);

      // Loop over DSECAL y hits to get residuals for the TPC track
      //if (neyhits>0) tyeT2->Fill(ty);
      int selEyhits = 0;
      for(int i = 0; i < skt2kalignvar->Sneyhits; i++ ){
  if (skt2kalignvar->SneyLayer[i] == 1) { // select first y layer
    double charge = skt2kalignvar->Seqy[i]; 
    if (charge>3.0 && charge<10.0) { // charge cut to remove noise hits
      dy2[4]->Fill(skt2kalignvar->Sedy[i]);
      dyvsy2[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sedy[i]);
      //fill a dyvsy for X +/-
      //we could use Sxx for this, but there was a mistake in filling this for version9, so we use Sx+Sdx
      if(skt2kalignvar->Sex[i]+skt2kalignvar->Sedx[i] >= 0.)
        {
    dyvsyP2[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sedy[i]);
        }
      else
        {
    dyvsyM2[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sedy[i]);
        }
      dyvsx2[4]->Fill(skt2kalignvar->Seyx[i],skt2kalignvar->Sedy[i]);
    
      dyvsty2[4]->Fill(skt2kalignvar->Sty,skt2kalignvar->Sedy[i]);
      dyvstx3e->Fill(skt2kalignvar->Stx,skt2kalignvar->Sedy[i]);
    
      tyvsx2[4]->Fill(skt2kalignvar->Seyx[i],skt2kalignvar->Sty);
      tyvsy2[4]->Fill(skt2kalignvar->Sey[i],skt2kalignvar->Sty);
      selEyhits++;
    }
  }
      }
      if (selEyhits>0) tyeT2->Fill(skt2kalignvar->Sty);
      if (selEyhits>0) ny2[4]->Fill(selEyhits);


      // Loop over DSECAL x hits to get residuals for the TPC track
      //if (nexhits>0) txeT2->Fill(tx);
      int selExhits = 0;
      for(int i = 0; i < skt2kalignvar->Snexhits; i++ ){
  if (skt2kalignvar->SnexLayer[i] == 1) { // select first x layer
    double charge = skt2kalignvar->Seqx[i]; 
    if (charge>3.0 && charge<10.0) { // charge cut to remove noise hits
      dx2[4]->Fill(skt2kalignvar->Sedx[i]);
      if( skt2kalignvar->Sex[i] < 0 ) dx3eN2->Fill(skt2kalignvar->Sedx[i]);
      if( skt2kalignvar->Sex[i] > 0 ) dx3eP2->Fill(skt2kalignvar->Sedx[i]);
      dxvsx2[4]->Fill(skt2kalignvar->Sex[i],skt2kalignvar->Sedx[i]);
      dxvsy2[4]->Fill(skt2kalignvar->Sexy[i],skt2kalignvar->Sedx[i]);
    
      dxvstx2[4]->Fill(skt2kalignvar->Stx,skt2kalignvar->Sedx[i]);
      dxvsty3e2->Fill(skt2kalignvar->Sty,skt2kalignvar->Sedx[i]);
    
      txvsx2[4]->Fill(skt2kalignvar->Sex[i],skt2kalignvar->Stx);
      txvsy2[4]->Fill(skt2kalignvar->Sexy[i],skt2kalignvar->Stx);
      selExhits++;
    }
  }
      }
      if (selExhits>0) txeT2->Fill(skt2kalignvar->Stx);
      if (selExhits>0) nx2[4]->Fill(selExhits);

    }

  oconstfl<<avegtimeruns<<"\t";
  oconstfl<<avegtimerunsE<<"\t";

  //TFile * out;  // canvas output file
  //out = new TFile("tpcfgdalignhisto.root","recreate");

  // plot them ************************************************
  TCanvas *c0;
  TCanvas *c1;
  TCanvas *c1a;
  TCanvas *c1b;
  TCanvas *c1bx;
  TCanvas *c2;
  TCanvas *c3;
  TCanvas *c4;
  TCanvas *c5;
  TCanvas *c6;
  TCanvas *c7;
  TCanvas *c8;
  TCanvas *c9;
  TCanvas *c10;
  TCanvas *c11;
  TCanvas *c12;
  TCanvas *c13;
  TCanvas *c14;
  TCanvas *c15;
  TCanvas *c16;

  c0 = new TCanvas("c0","dy vs x");
  c0->Divide(2,2);
  vector <TLegend* > leg0;
  for(Int_t ci=0; ci<4; ci++)
    {
      c0->cd(ci+1);
 
      const int nhist = 2;//6;
 
      TH1F *hcol[nhist] = {fit1d(fitslices(mupulldyvsx[ci], 1), 1, "pol1", -700., 700.),
         fit1d(fitslices(mupulldyvsx2[ci], 8), 8, "pol1", -700., 700.)
      };
 
      std::string legname[nhist] = {

  "misaligned",
  "original"
  //"trkqual,mupull B>0",
  //"trkqual,mupull B=0"
      };

      hcol[0]->SetTitle("");
      hcol[0]->Draw();
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      //leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      //leg0[ci]->Draw("same");
      leg0.push_back((TLegend*)legshape3(nhist, hcol,legname));
      leg0[ci]->Draw("same");
    }
 
  c0->Update();
  c0->Print("tpcfgd_dyVSx_"+oFile+".eps");


  c1 = new TCanvas("c1","dy vs y");
  c1->Divide(2,2);
  vector <TLegend* > leg1;
  for(Int_t ci=0; ci<4; ci++)
    {
      c1->cd(ci+1);

      const int nhist = 2;//6;
 
      TH1F *hcol[nhist] = {fit1d(fitslices(mupulldyvsy[ci], 1), 1, "pol1", -700., 700.),
         fit1d(fitslices(mupulldyvsy2[ci], 8), 8, "pol1", -700., 700.)
      };

  
      std::string legname[nhist] = {

  "misaligned",
  "original"
  //"trkqual,mupull B>0",
  //"trkqual,mupull B=0"
      };

      hcol[0]->Draw();
      hcol[0]->SetTitle("");
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      //leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      //leg0[ci]->Draw("same");
      leg1.push_back((TLegend*)legshape3(nhist, hcol,legname));
      leg1[ci]->Draw("same");

    }

  c1->Update();
  c1->Print("tpcfgd_dyVSy_"+oFile+".eps");

  c1a = new TCanvas("c1a","dy vs y for XPM");
  vector <TLegend* > leg1a;
  c1a->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c1a->cd(ci+1);

      const int nhist = 4;//6;

      TH1F *hcol[nhist] = {fitslices(dyvsyP[ci], 1),
         fitslices(dyvsyP2[ci], 8),
         fitslices(dyvsyM[ci], 2),
         fitslices(dyvsyM2[ci], 5)
      };


      std::string legname[nhist] = {
  "misaligned, X>=0",
  "original, X>=0",
  "misaligned, X<0",
  "original, X<0"
      };

      hcol[0]->Draw();
      hcol[0]->SetTitle("");
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      leg1a.push_back((TLegend*)legshape1(nhist, hcol,legname));
      leg1a[ci]->Draw("same");


    }

  c1a->Update();
  c1a->Print("tpcfgd_dyVSyforXPM_"+oFile+".eps");

  c1b = new TCanvas("c1b","Sy");
  vector <TLegend* > leg1b;
  c1b->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c1b->cd(ci+1);

      const int nhist = 2;//6;
      TH1F *hcol42[nhist] = {colorhists(fgdy[ci], 1),colorhists(fgdy2[ci], 8)};
      hcol42[1]->Draw();
      hcol42[0]->Draw("same");


      std::string legname[nhist] = {
  "misaligned",
  "original"
      };

      leg1b.push_back((TLegend*)legshape2(nhist, hcol42,legname));
      leg1b[ci]->Draw("same");

    }

  c1b->Update();
  c1b->Print("tpcfgd_Sy_"+oFile+".eps");

  c1bx = new TCanvas("c1bx","Sx");
  vector <TLegend* > leg1bx;
  c1bx->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c1bx->cd(ci+1);

      const int nhist = 2;//6;
   
      TH1F *hcol42[nhist] = {colorhists(fgdx[ci], 1),colorhists(fgdx2[ci], 8)};
      hcol42[1]->Draw();
      hcol42[0]->Draw("same");

      std::string legname[nhist] = {
  "misal",
  "orig"
      };

 
 
      leg1bx.push_back((TLegend*)legshape1(nhist, hcol42,legname));
      leg1bx[ci]->Draw("same");

    }

  c1bx->Update();
  c1bx->Print("tpcfgd_Sx_"+oFile+".eps");

  TCanvas *c1tx = new TCanvas("c1tx","dx vs xtx", 1300, 800);
  c1tx->Divide(2,2);
  vector <TLegend* > leg1tx;
  for(Int_t ci=0; ci<4; ci++)
    {
      c1tx->cd(ci+1);
      c1tx->SetRightMargin(0.02);
      const int nhist = 2;//6;
 
      TH1F *hcol[nhist] = {fit1d(fitslices(mupulldxvsxtx[ci], 1), 1, "pol1", -300., 300.),
         fit1d(fitslices(mupulldxvsxtx2[ci], 8), 8, "pol1", -300., 300.)
      };

  
      std::string legname[nhist] = {

  "mis",
  "orig"
  //"trkqual,mupull B>0",
  //"trkqual,mupull B=0"
      };

      char fitname[256];
      strcpy (fitname,hcol[0]->GetName());
      strcat (fitname,"Fit");
      
      oconstfl<<hcol[0]->GetFunction(fitname)->GetParameter(1)<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetParError(1)<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetChisquare()<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetNDF()<<"\t";

      hcol[0]->Draw();
      hcol[0]->SetTitle("");
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      //leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      //leg0[ci]->Draw("same");
      leg1tx.push_back((TLegend*)legshape3(nhist, hcol,legname));
      leg1tx[ci]->Draw("same");

    }

  c1tx->Update();
  c1tx->Print("tpcfgd_dxVSxtx_"+oFile+".eps");

  TCanvas *c1ty = new TCanvas("c1ty","dy vs yty", 1300, 800);
  c1ty->Divide(2,2);
  vector <TLegend* > leg1ty;
  for(Int_t ci=0; ci<4; ci++)
    {
      c1ty->cd(ci+1);
      c1ty->SetRightMargin(0.02);
      const int nhist = 2;//6;
 
      TH1F *hcol[nhist] = {fit1d(fitslices(mupulldyvsyty[ci], 1), 1, "pol1", -300., 300.),
         fit1d(fitslices(mupulldyvsyty2[ci], 8), 8, "pol1", -300., 300.)
      };

  
      std::string legname[nhist] = {

  "mis",
  "orig"
  //"trkqual,mupull B>0",
  //"trkqual,mupull B=0"
      };

     char fitname[256];
     strcpy (fitname,hcol[0]->GetName());
     strcat (fitname,"Fit");
     
     oconstfl<<hcol[0]->GetFunction(fitname)->GetParameter(1)<<"\t";
     oconstfl<<hcol[0]->GetFunction(fitname)->GetParError(1)<<"\t";
     oconstfl<<hcol[0]->GetFunction(fitname)->GetChisquare()<<"\t";
     oconstfl<<hcol[0]->GetFunction(fitname)->GetNDF()<<"\t";

      hcol[0]->Draw();
      hcol[0]->SetTitle("");
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      //leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      //leg0[ci]->Draw("same");
      leg1ty.push_back((TLegend*)legshape3(nhist, hcol,legname));
      leg1ty[ci]->Draw("same");

    }

  c1ty->Update();
  c1ty->Print("tpcfgd_dyVSyty_"+oFile+".eps");

  TCanvas *c1tya = new TCanvas("c1tya","dy vs yty(2D)");
  c1tya->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c1tya->cd(ci+1);
      mupulldyvsyty[ci]->Draw("colz");
    }

  TCanvas *c1tyb = new TCanvas("c1tyb","dy vs yty(2D) nonmisaligned");
  c1tyb->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c1tyb->cd(ci+1);
      mupulldyvsyty2[ci]->Draw("colz");
    }


  c2 = new TCanvas("c2","dx vs x (FGD#:TPC#)");
  c2->Divide(2,2); 
  vector <TLegend* > leg2;
  for(Int_t ci=0; ci<4; ci++)
    {
      c2->cd(ci+1);
      const int nhist = 2;//6;
    
      TH1F *hcol[nhist] = {fit1d(fitslices(mupulldxvsx[ci], 1), 1, "pol1", -700., 700.),
         fit1d(fitslices(mupulldxvsx2[ci], 8), 8, "pol1", -700., 700.)
      };
   
      std::string legname[nhist] = {

  "misaligned",

  "original"
      };

      hcol[0]->Draw();
      hcol[0]->SetTitle("");
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      //leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      //leg0[ci]->Draw("same");
      leg2.push_back((TLegend*)legshape3(nhist, hcol,legname));
      leg2[ci]->Draw("same");

    }


  c2->Update();
  c2->Print("tpcfgd_dxVSx_"+oFile+".eps");


  c3 = new TCanvas("c3","dx vs y (FGD#:TPC#)");
  c3->Divide(2,2); 
  vector <TLegend* > leg3;
  for(Int_t ci=0; ci<4; ci++)
    {
      c3->cd(ci+1);
      const int nhist = 2;//6;
  
       TH1F *hcol[nhist] = {fit1d(fitslices(mupulldxvsy[ci], 1), 1, "pol1", -700., 1000.),
          fit1d(fitslices(mupulldxvsy2[ci], 8), 8, "pol1", -700., 1000.)

          //     TH1F *hcol[nhist] = {fit1d(fitslices(mupulldxvsy2[ci], 8), 8, "pol1", -700., 700.),
          //fit1d(fitslices(mupulldxvsy[ci], 1), 1, "pol1", -700., 700.)
      };
 
      std::string legname[nhist] = {

  "misaligned",

  "original"
      };

      char fitname[256];
      strcpy (fitname,hcol[0]->GetName());
      strcat (fitname,"Fit");
      
      oconstfl<<hcol[0]->GetFunction(fitname)->GetParameter(0)<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetParError(0)<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetParameter(1)<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetParError(1)<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetChisquare()<<"\t";
      oconstfl<<hcol[0]->GetFunction(fitname)->GetNDF()<<"\t";

      hcol[0]->Draw();
      hcol[0]->SetTitle("");
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }


 
      //leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      //leg0[ci]->Draw("same");
      leg3.push_back((TLegend*)legshape3(nhist, hcol,legname));
      leg3[ci]->Draw("same");
 
    }
 
  c3->Update();
  c3->Print("tpcfgd_dxVSy_"+oFile+".eps");



  c4 = new TCanvas("c4","dy and dx");
  c4->Divide(2,2); 
  c4->cd(1);

  const int nhist4 = 2;
  TH1F *hcol41[nhist4] = {fit1d(dyT, 1, "gaus", -40.,40.),fit1d(dyT2, 8, "gaus", -40.,40.)};
  std::string legname4[nhist4] = {"bothcuts misaligned","bothcuts original"};
  
  hscale(nhist4,hcol41);
  hcol41[0]->Draw();
  hcol41[1]->Draw("same");
  //gPad->SetLogy();
  
  TLegend *leg4 = (TLegend*)legshape2(nhist4, hcol41,legname4);
  leg4->Draw("same");

 
  c4->cd(2); 
  TH1F *hcol42[nhist4] = {fit1d(dxT, 1, "gaus", -40.,40.),fit1d(dxT2, 8, "gaus", -40.,40.)};
  hcol42[0]->Draw();
  hcol42[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol42,legname4);
  leg4->Draw("same");

  c4->cd(3); 
  TH1F *hcol43[nhist4] = {fit1d(dxN, 1, "gaus", -40.,40.),fit1d(dxN2, 8, "gaus", -40.,40.)};
  hcol43[0]->Draw();
  hcol43[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol43,legname4);
  leg4->Draw("same");

  c4->cd(4); 
  TH1F *hcol44[nhist4] = {fit1d(dxP, 1, "gaus", -40.,40.),fit1d(dxP2, 8, "gaus", -40.,40.)};
  hcol44[0]->Draw();
  hcol44[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol44,legname4);
  leg4->Draw("same");

  c4->Update();
  c4->Print("tpcfgd04_"+oFile+".eps");


  c5 = new TCanvas("c5","dy");
  c5->Divide(2,2); 
  for(Int_t ci=0; ci<4; ci++)
    {
      c5->cd(ci+1); 

      const int nhist = 2;
 
      TH1F *hcol[nhist] = {fit1d(dy[ci], 1, "gaus", -40.,40.),fit1d(dy2[ci], 8, "gaus", -40.,40.)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };
      
      oconstfl<<hcol[0]->GetMean()<<"\t";
      oconstfl<<hcol[0]->GetMeanError()<<"\t";
      //oconstfl<<hcol[0]->GetFunction(fitname)->GetChisquare()<<"\t";
      //oconstfl<<hcol[0]->GetFunction(fitname)->GetNDF()<<"\t";

      hscale(nhist,hcol);
      hcol[0]->Draw();
      hcol[1]->Draw("same");
      //gPad->SetLogy();
  
      TLegend *leg = ((TLegend*)legshape2(nhist, hcol,legname));
      leg->Draw("same");

 
    }


  c5->Update();
  c5->Print("tpcfgd05_"+oFile+".eps");


  c6 = new TCanvas("c6","dx");
  c6->Divide(2,2); 
  for(Int_t ci=0; ci<4; ci++)
    {
      c6->cd(ci+1); 
      const int nhist = 2;
  
      TH1F *hcol[nhist] = {fit1d(dx[ci], 1, "gaus", -40.,40.),fit1d(dx2[ci], 8, "gaus", -40.,40.)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };

      oconstfl<<hcol[0]->GetMean()<<"\t";
      oconstfl<<hcol[0]->GetMeanError()<<"\t";

      hscale(nhist,hcol);
      hcol[0]->Draw();
      hcol[1]->Draw("same");
      //gPad->SetLogy();
     
      TLegend *leg = ((TLegend*)legshape2(nhist, hcol,legname));
      leg->Draw("same");

 
    }

 
  c6->Update();
  c6->Print("tpcfgd06_"+oFile+".eps");


  c7 = new TCanvas("c7","dx vs tx");
  c7->Divide(2,2); 
  for(Int_t ci=0; ci<4; ci++)
    {
      c7->cd(ci+1); 
 
      const int nhist = 2;
      TH1F *hcol[nhist] = {fit1d(fitslices(dxvstx[ci], 1), 1, "pol1", -10., 10.),fit1d(fitslices(dxvstx2[ci], 8), 8, "pol1", -10., 10.)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };

      hcol[0]->Draw();
      hcol[1]->Draw("same");
       
      TLegend *leg = ((TLegend*)legshape1(nhist, hcol,legname));
      leg->Draw("same");

 
    }

  c7->Update(); 
  c7->Print("tpcfgd07_"+oFile+".eps");

  TCanvas *c7b = new TCanvas("c7b","dx vs ty");
  c7b->Divide(2,2); 
  for(Int_t ci=0; ci<4; ci++)
    {
      c7b->cd(ci+1); 
 
      const int nhist = 2;
      TH1F *hcol[nhist] = {fit1d(fitslices(dxvsty[ci], 1), 1, "pol1", -10., 10.),fit1d(fitslices(dxvsty2[ci], 8), 8, "pol1", -10., 10.)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };

      hcol[0]->Draw();
      hcol[1]->Draw("same");
       
      TLegend *leg = ((TLegend*)legshape1(nhist, hcol,legname));
      leg->Draw("same");

 
    }

  c7b->Update(); 
  c7b->Print("tpcfgd07b_"+oFile+".eps");


  c8 = new TCanvas("c8","dy vs ty");
  c8->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c8->cd(ci+1); 

      const int nhist = 2;
      TH1F *hcol[nhist] = {fit1d(fitslices(dyvsty[ci], 1), 1, "pol1", -10., 10.),fit1d(fitslices(dyvsty2[ci], 8), 8, "pol1", -10., 10.)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };

      hcol[0]->Draw();
      hcol[1]->Draw("same");
       
      TLegend *leg = ((TLegend*)legshape1(nhist, hcol,legname));
      leg->Draw("same");

  

    } 

  c8->Update(); 
  c8->Print("tpcfgd08_"+oFile+".eps");

  //draw y-projections of the dy vs ty histos
  //TCanvas *c8a = new TCanvas("c8a","dy vs ty slices", 2000.,2000.);
  TCanvas *c8a = new TCanvas("c8a","dy vs ty slices", 200.,200.);
  c8a->Divide(2,2);
  TH1 *hbins[4][dyvsty[0]->GetNbinsX()] ;
  TH1 *hbins2[4][dyvsty2[0]->GetNbinsX()] ;
  for(int di=0; di<4; di++)
    {
      TVirtualPad *pp = c8a->cd(di+1);
      pp->Divide(2,5);
      //Int_t hcnt = 3;
      for (int hi=0;hi<dyvsty[di]->GetNbinsX();hi++) {
  char name[256];
  sprintf(name,"%d_%d",di+1,hi+1); 
  pp->cd(hi+1);
  //hcnt++;
  //c8a->cd("3_1");
  hbins[di][hi] = dyvsty[di]->ProjectionY(Form("bin%d",hi+1),hi+1,hi+1);
  hbins[di][hi]->SetLineColor(2);
  hbins[di][hi]->Draw();
  hbins2[di][hi] = dyvsty2[di]->ProjectionY(Form("bin2_%d",hi+1),hi+1,hi+1);
  hbins2[di][hi]->SetLineColor(9);
  hbins2[di][hi]->Draw();
  hbins[di][hi]->Draw("same");
  //if(hi==0){hbins[di][hi]->Draw();}
  //else {hbins[di][hi]->Draw("same");}
      }
    }

  //draw y-projections of the dy vs ty histos
  //TCanvas *c8b = new TCanvas("c8b","dy vs ty slices 2", 2000.,2000.);
  TCanvas *c8b = new TCanvas("c8b","dy vs ty slices 2", 200.,200.);
  c8b->Divide(2,2);
  //TH1 *hbins2[4][dyvsty2[0]->GetNbinsX()] ;
  for(int di=0; di<4; di++)
    {
      TVirtualPad *pp2 = c8b->cd(di+1);
      pp2->Divide(2,5);
      //Int_t hcnt = 3;
      for (int hi=0;hi<dyvsty2[di]->GetNbinsX();hi++) {
  char name[256];
  sprintf(name,"%d_%d",di+1,hi+1); 
  pp2->cd(hi+1);
  //hcnt++;
  //c8a->cd("3_1");
  hbins2[di][hi] = dyvsty2[di]->ProjectionY(Form("bin2_%d",hi+1),hi+1,hi+1);
  hbins2[di][hi]->SetLineColor(9);
  hbins2[di][hi]->Draw();
  //if(hi==0){hbins[di][hi]->Draw();}
  //else {hbins[di][hi]->Draw("same");}
      }
    }

  c9 = new TCanvas("c9","dx,dy vs tx,ty");
  c9->Divide(2,2);
  c9->cd(1);
  TH1F *hcol91[nhist4] = {fit1d(fitslices(dyvstyT, 1), 1, "pol1", -10., 10.),
        fit1d(fitslices(dyvstyT2, 8), 8, "pol1", -10., 10.)};
  hcol91[0]->Draw();
  hcol91[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol91,legname4);
  leg4->Draw("same");
 
  c9->cd(2);
  TH1F *hcol92[nhist4] = {fit1d(fitslices(dyvstxT, 1), 1, "pol1", -10., 10.),
        fit1d(fitslices(dyvstxT2, 8), 8, "pol1", -10., 10.)};
  hcol92[0]->Draw();
  hcol92[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol92,legname4);
  leg4->Draw("same");
 
  c9->cd(3);
  TH1F *hcol93[nhist4] = {fit1d(fitslices(dxvstyT, 1), 1, "pol1", -10., 10.),
        fit1d(fitslices(dxvstyT2, 8), 8, "pol1", -10., 10.)};
  hcol93[0]->Draw();
  hcol93[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol93,legname4);
  leg4->Draw("same");

  c9->cd(4);
  TH1F *hcol94[nhist4] = {fit1d(fitslices(dxvstxT, 1), 1, "pol1", -10., 10.),
        fit1d(fitslices(dxvstxT2, 8), 8, "pol1", -10., 10.)};
  hcol94[0]->Draw();
  hcol94[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol94,legname4);
  leg4->Draw("same");

  c9->Update(); 
  c9->Print("tpcfgd09_"+oFile+".eps");


  c10 = new TCanvas("c10","nyhits");
  c10->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c10->cd(ci+1);
      const int nhist = 2;
      TH1F *hcol[nhist] = {colorhists(ny[ci],1),colorhists(ny2[ci],8)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };
      
      hcol[0]->Draw();
      hcol[1]->Draw("same");
       
      TLegend *leg = ((TLegend*)legshape2(nhist, hcol,legname));
      leg->Draw("same");

  
    }


  c10->Print("tpcfgd10_"+oFile+".eps");


  c11 = new TCanvas("c11","nxhits");
  c11->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c11->cd(ci+1);
      
      const int nhist = 2;
      TH1F *hcol[nhist] = {colorhists(nx[ci],1),colorhists(nx2[ci],8)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };
      
      hcol[0]->Draw();
      hcol[1]->Draw("same");
       
      TLegend *leg = ((TLegend*)legshape2(nhist, hcol,legname));
      leg->Draw("same");

 
    }  

 
  c11->Print("tpcfgd11_"+oFile+".eps");


  c12 = new TCanvas("c12","dy and dx");
  c12->Divide(2,2); 
  c12->cd(1); 
  TH1F *hcol121[nhist4] = {fit1d(dy[4], 1, "gaus", -40.,40.),fit1d(dy2[4], 8, "gaus", -40.,40.)};
  hcol121[0]->Draw();
  hcol121[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol121,legname4);
  leg4->Draw("same");
 
  c12->cd(2); 
  TH1F *hcol122[nhist4] = {fit1d(dx[4], 1, "gaus", -40.,40.),fit1d(dx2[4], 8, "gaus", -40.,40.)};
  hcol122[0]->Draw();
  hcol122[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol122,legname4);
  leg4->Draw("same");

  c12->cd(3); 
  TH1F *hcol123[nhist4] = {fit1d(dx3eN, 1, "gaus", -40.,40.),fit1d(dx3eN2, 8, "gaus", -40.,40.)};
  hcol123[0]->Draw();
  hcol123[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol123,legname4);
  leg4->Draw("same");

  c12->cd(4); 
  TH1F *hcol124[nhist4] = {fit1d(dx3eP, 1, "gaus", -40.,40.),fit1d(dx3eP2, 8, "gaus", -40.,40.)};
  hcol124[0]->Draw();
  hcol124[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol124,legname4);
  leg4->Draw("same");
 
  c12->Update();
  c12->Print("tpcfgd12_"+oFile+".eps");


  c13 = new TCanvas("c13","dy, dx vs y, x");
  c13->Divide(2,2); 
  c13->cd(1);
  TH1F *hcol131[nhist4] = {colorhists(fitslices2(dyvsy[4], 0),1),colorhists(fitslices2(dyvsy2[4], 0),8)};
  hcol131[0]->Draw();
  hcol131[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol131,legname4);
  leg4->Draw("same");

  c13->cd(2);
  TH1F *hcol132[nhist4] = {colorhists(fitslices2(dyvsx[4], 0),1),colorhists(fitslices2(dyvsx2[4], 0),8)};
  hcol132[0]->Draw();
  hcol132[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol132,legname4);
  leg4->Draw("same");

  c13->cd(3);
  TH1F *hcol133[nhist4] = {colorhists(fitslices2(dxvsx[4], 0),1),colorhists(fitslices2(dxvsx2[4], 0),8)};
  hcol133[0]->Draw();
  hcol133[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol133,legname4);
  leg4->Draw("same");

  c13->cd(4);
  TH1F *hcol134[nhist4] = {colorhists(fitslices2(dxvsy[4], 0),1),colorhists(fitslices2(dxvsy2[4], 0),8)};
  hcol134[0]->Draw();
  hcol134[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol134,legname4);
  leg4->Draw("same");

  c13->Update();
  c13->Print("tpcfgd13_"+oFile+".eps");


  c14 = new TCanvas("c14","dy, dx vs ty, tx");
  c14->Divide(2,2); 
  c14->cd(1);
  TH1F *hcol141[nhist4] = {colorhists(fitslices2(dyvsty[4], 0),1),colorhists(fitslices2(dyvsty2[4], 0),8)};
  hcol141[0]->Draw();
  hcol141[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol141,legname4);
  leg4->Draw("same");

  c14->cd(2);
  TH1F *hcol142[nhist4] = {colorhists(fitslices2(dyvstx3e, 0),1),colorhists(fitslices2(dyvstx3e2, 0),8)};
  hcol142[0]->Draw();
  hcol142[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol142,legname4);
  leg4->Draw("same");

  c14->cd(3);
  TH1F *hcol143[nhist4] = {colorhists(fitslices2(dxvstx[4], 0),1),colorhists(fitslices2(dxvstx2[4], 0),8)};
  hcol143[0]->Draw();
  hcol143[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol143,legname4);
  leg4->Draw("same");

  c14->cd(4);
  TH1F *hcol144[nhist4] = {colorhists(fitslices2(dxvsty3e, 0),1),colorhists(fitslices2(dxvsty3e2, 0),8)};
  hcol144[0]->Draw();
  hcol144[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol144,legname4);
  leg4->Draw("same");

  c14->Update();
  c14->Print("tpcfgd14_"+oFile+".eps");


  c15 = new TCanvas("c15","ny, nx, ty, tx");
  c15->Divide(2,2); 
  c15->cd(1);
  TH1F *hcol151[nhist4] = {colorhists(ny[4],1),colorhists(ny2[4],8)};
  hcol151[0]->Draw();
  hcol151[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol151,legname4);
  leg4->Draw("same");
 
  c15->cd(2);
  TH1F *hcol152[nhist4] = {colorhists(nx[4],1),colorhists(nx2[4],8)};
  hcol152[0]->Draw();
  hcol152[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol152,legname4);
  leg4->Draw("same");

  c15->cd(3);
  TH1F *hcol153[nhist4] = {colorhists(tyeT,1),colorhists(tyeT2,8)};
  hcol153[0]->Draw();
  hcol153[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol153,legname4);
  leg4->Draw("same");
   
  c15->cd(4);
  TH1F *hcol154[nhist4] = {colorhists(txeT,1),colorhists(txeT2,8)};
  hcol154[0]->Draw();
  hcol154[1]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol154,legname4);
  leg4->Draw("same");
   
  c15->Update();
  c15->Print("tpcfgd15_"+oFile+".eps");


  c16 = new TCanvas("c16","ny, nx, ty, tx");
  c16->Divide(2,2); 
  c16->cd(1);
  TH1F *hcol161[nhist4] = {colorhists(nyT,1),colorhists(nyT2,8)};
  hcol161[1]->Draw();
  hcol161[0]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol161,legname4);
  leg4->Draw("same");
 
  c16->cd(2);
  TH1F *hcol162[nhist4] = {colorhists(nxT,1),colorhists(nxT2,8)};
  hcol162[1]->Draw();
  hcol162[0]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol162,legname4);
  leg4->Draw("same");

  c16->cd(3);
  TH1F *hcol163[nhist4] = {colorhists(tyT,1),colorhists(tyT2,8)};
  hcol163[1]->Draw();
  hcol163[0]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol163,legname4);
  leg4->Draw("same");
  
  c16->cd(4);
  TH1F *hcol164[nhist4] = {colorhists(txT,1),colorhists(txT2,8)};
  hcol164[1]->Draw();
  hcol164[0]->Draw("same");
  leg4 = (TLegend*)legshape2(nhist4, hcol164,legname4);
  leg4->Draw("same");

  c16->Update();
  c16->Print("tpcfgd16_"+oFile+".eps");

  TCanvas *c17 = new TCanvas("c17","reduced #Chi^{2}");
  c17->Divide(3,1);
  for(Int_t hii=0; hii<3; hii++)
    {
      c17->cd(hii+1);
      const int nhist = 2;
      TH1F *hcol[nhist] = {colorhists(redchitpctrk[hii],1),colorhists(redchitpctrk2[hii],8)};
      //pois.push_back(new TF1(Form("poihisto%d",hj),"expo",201.,10000.)); //fitting from 0 to 10000
      //TH1F *hcol[nhist] = {fit1d(redchitpctrk[hii], 1, "expo", 0.1,20.),fit1d(redchitpctrk2[hii], 8, "expo", 0.1,20.)};
      std::string legname[nhist] = {
  "bothcuts misaligned",
  "bothcuts original"
      };


       
      TLegend *leg = ((TLegend*)legshape2(nhist, hcol,legname));
      leg->Draw("same");
 
      char fitname[256];
      strcpy (fitname,hcol[0]->GetName());
      strcat (fitname,"Fit");

      //TF1 *pois = new TF1(Form("poihisto%d",hj),"expo",201.,10000.));
      TF1 *poiso = new TF1("poihisto","expo",0.5,5.0);
      poiso->SetParameters(5.,1.);
      //poiso->SetLineColor(2);
      hcol[0]->Fit("poihisto","emr");
      //out->WriteTObject(poiso);
      oconstfl<<poiso->GetParameter(0)<<"\t";
      oconstfl<<poiso->GetParError(0)<<"\t";
      oconstfl<<poiso->GetParameter(1)<<"\t";
      oconstfl<<poiso->GetParError(1)<<"\t";
      hcol[0]->Draw();
      //poiso->Draw("lsame");
      //hcol[1]->Draw("same");
      poiso->Clear();
      //delete [] pois;
 

    }
  c17->Print("reducedchisq_"+oFile+".eps");

  TCanvas *c18 = new TCanvas("c18","dy vs y (2D)");
  c18->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c18->cd(ci+1);
      dyvsy[ci]->Draw("colz");
    }
  c18->Update();
  c18->Print("tpcfgd_dyVSy2d_"+oFile+".eps");

  TCanvas *c19 = new TCanvas("c19","dy vs x (2D)");
  c19->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c19->cd(ci+1);
      dyvsx[ci]->Draw("colz");
    }
  c19->Update();
  c19->Print("tpcfgd_dyVSx2d_"+oFile+".eps");


  TCanvas *c20 = new TCanvas("c20","dx vs x (2D)");
  c20->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c20->cd(ci+1);
      dxvsx[ci]->Draw("colz");
    }
  c20->Update();
  c20->Print("tpcfgd_dxVSx2d_"+oFile+".eps");

  TCanvas *c21 = new TCanvas("c21","dx vs y (2D)");
  c21->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c21->cd(ci+1);
      dxvsy[ci]->Draw("colz");
    }
  c21->Update();
  c21->Print("tpcfgd_dxVSy2d_"+oFile+".eps");


  TCanvas *c22 = new TCanvas("c22","Sx and Sy ");
  c22->Divide(2,3);
  
  c22->cd(1);
  TH1F *hcolSxf1[2] = {colorhists(SxTf1, 1),colorhists(SxT2f1, 8)};
  hcolSxf1[1]->Draw();
  hcolSxf1[0]->Draw("same");

  c22->cd(2);
  TH1F *hcolSyf1[2] = {colorhists(SyTf1, 1),colorhists(SyT2f1, 8)};
  hcolSyf1[1]->Draw();
  hcolSyf1[0]->Draw("same");

  c22->cd(3);
  TH1F *hcolSxf2[2] = {colorhists(SxTf2, 1),colorhists(SxT2f2, 8)};
  hcolSxf2[1]->Draw();
  hcolSxf2[0]->Draw("same");

  c22->cd(4);
  TH1F *hcolSyf2[2] = {colorhists(SyTf2, 1),colorhists(SyT2f2, 8)};
  hcolSyf2[1]->Draw();
  hcolSyf2[0]->Draw("same");

  TBox *yzt1 = new TBox(-763.581,-999.667*2.,-33.7323,1100.44*2.);
  TBox *yzt2 = new TBox(593.677,-998.687*2.,1323.99,1101.7*2.);
  TBox *yzt3 = new TBox(1954.07,-997.865*2.,2685.7,1102.13*2.);

  TBox *xzt1 = new TBox(-763.581,-905.887*2.,-33.7323,901.702*2.);
  TBox *xzt2 = new TBox(593.677,-904.871*2.,1323.99,901.079*2.);
  TBox *xzt3 = new TBox(1954.07,-906.725*2.,2685.7,900.526*2.);


  c22->cd(5);
  TH1F *hcolSxz[2] = {colorhists(SxzT, 1),colorhists(SxzT2, 8)};
  hcolSxz[1]->Draw();
  hcolSxz[0]->Draw("same");
  xzt1->Draw("same");
  xzt2->Draw("same");
  xzt3->Draw("same");

  c22->cd(6);
  TH1F *hcolSyz[2] = {colorhists(SyzT, 1),colorhists(SyzT2, 8)};
  hcolSyz[1]->Draw();
  hcolSyz[0]->Draw("same");
  yzt1->Draw("same");
  yzt2->Draw("same");
  yzt3->Draw("same");

  c22->Update();
  c22->Print("tpcfgd_SxSy_"+oFile+".eps");

  
  TCanvas *c23 = new TCanvas("c23","std vs new :dy ");
  c23->Divide(2,1);
  
  c23->cd(1);
  newstddy->SetMarkerSize(0.2);
  newstddy->SetMarkerColor(2);
  newstddy->SetMarkerStyle(24);
  newstddy->Draw("colz");
  newstddytest1->SetMarkerSize(0.4);
  newstddytest1->SetMarkerColor(3);
  newstddytest1->SetMarkerStyle(24);
  //newstddytest1->Draw("sames");
  newstddytest2->SetMarkerSize(0.4);
  newstddytest2->SetMarkerColor(4);
  newstddytest2->SetMarkerStyle(24);
  //newstddytest2->Draw("sames");
  newstddytest3->SetMarkerSize(0.4);
  newstddytest3->SetMarkerColor(5);
  newstddytest3->SetMarkerStyle(24);
  //newstddytest3->Draw("sames");
  newstddytest4->SetMarkerSize(0.4);
  newstddytest4->SetMarkerColor(6);
  newstddytest4->SetMarkerStyle(24);
  //newstddytest4->Draw("sames");
  newstddytest5->SetMarkerSize(0.4);
  newstddytest5->SetMarkerColor(7);
  newstddytest5->SetMarkerStyle(24);
  //newstddytest5->Draw("sames");
  newstddytest6->SetMarkerSize(0.4);
  newstddytest6->SetMarkerColor(8);
  newstddytest6->SetMarkerStyle(24);
  //newstddytest6->Draw("sames");
  newstddytest7->SetMarkerSize(0.4);
  newstddytest7->SetMarkerColor(9);
  newstddytest7->SetMarkerStyle(24);
  //newstddytest7->Draw("sames");
  newstddytest8->SetMarkerSize(0.4);
  newstddytest8->SetMarkerColor(11);
  newstddytest8->SetMarkerStyle(24);
  //newstddytest8->Draw("sames");
  newstddytest9->SetMarkerSize(0.4);
  newstddytest9->SetMarkerColor(41);
  newstddytest9->SetMarkerStyle(28);
  //newstddytest9->Draw("sames");
  newstddytest10->SetMarkerSize(0.4);
  newstddytest10->SetMarkerColor(39);
  newstddytest10->SetMarkerStyle(28);
  //newstddytest10->Draw("sames");
  c23->cd(2);
  newstddy2->Draw("colz");


  TCanvas *c23b = new TCanvas("c23b","std vs new :dy ");
  c23b->Divide(6,2);
  c23b->Divide(3,3);
 
  c23b->cd(1);
  newstddy->SetMarkerSize(0.2);
  newstddy->SetMarkerColor(2);
  newstddy->SetMarkerStyle(24);
  newstddy->Draw();
  c23b->cd(2);
  newstddytest1->SetMarkerSize(0.4);
  newstddytest1->SetMarkerColor(3);
  newstddytest1->SetMarkerStyle(24);
  newstddytest1->Draw();
  c23b->cd(3);
  newstddytest2->SetMarkerSize(0.4);
  newstddytest2->SetMarkerColor(4);
  newstddytest2->SetMarkerStyle(24);
  newstddytest2->Draw();
  c23b->cd(4);
  newstddytest3->SetMarkerSize(0.4);
  newstddytest3->SetMarkerColor(5);
  newstddytest3->SetMarkerStyle(24);
  newstddytest3->Draw();
  c23b->cd(5);
  newstddytest4->SetMarkerSize(0.4);
  newstddytest4->SetMarkerColor(6);
  newstddytest4->SetMarkerStyle(24);
  newstddytest4->Draw();
  c23b->cd(6);
  newstddytest5->SetMarkerSize(0.4);
  newstddytest5->SetMarkerColor(7);
  newstddytest5->SetMarkerStyle(24);
  newstddytest5->Draw();
  c23b->cd(7);
  newstddytest6->SetMarkerSize(0.4);
  newstddytest6->SetMarkerColor(8);
  newstddytest6->SetMarkerStyle(24);
  newstddytest6->Draw();

  c23b->cd(8);
  newstddytest9->SetMarkerSize(0.4);
  newstddytest9->SetMarkerColor(41);
  newstddytest9->SetMarkerStyle(28);
  newstddytest9->Draw();
  c23b->cd(9);
  newstddytest10->SetMarkerSize(0.4);
  newstddytest10->SetMarkerColor(39);
  newstddytest10->SetMarkerStyle(28);
  newstddytest10->Draw();
  
  TCanvas *c23x = new TCanvas("c23x","std vs new :dx ");
  c23x->Divide(2,1);
  
  c23x->cd(1);
  newstddx->SetMarkerSize(0.2);
  newstddx->SetMarkerColor(2);
  newstddx->SetMarkerStyle(24);
  newstddx->Draw("colz");
  c23x->cd(2);
  newstddx2->SetMarkerSize(0.2);
  newstddx2->SetMarkerColor(2);
  newstddx2->SetMarkerStyle(24);
  newstddx2->Draw("colz");

  TCanvas *c24 = new TCanvas("c24","charge ");
  c24->Divide(2,1);
  
  c24->cd(1);
  hqy2->SetLineColor(9);
  hqy2->Draw();
  hqy1->SetLineColor(2);
  hqy1->Draw("same");
  c24->cd(2);
  hqx2->SetLineColor(9);
  hqx2->Draw();
  hqx1->SetLineColor(2);
  hqx1->Draw("same");

  TCanvas *c25 = new TCanvas("c25","ndof vs ty ");
  c25->cd();
  hndofty->Draw("colz");

  TCanvas *c26 = new TCanvas("c26","delta(dy and dx) vs y and x (misaligned)");
  c26->Divide(2,2);
  c26->cd(1);
  deltadyy->Draw("colz");
  c26->cd(2);
  deltadyx->Draw("colz");
  c26->cd(3);
  deltadxy->Draw("colz");
  c26->cd(4);
  deltadxx->Draw("colz");



  TCanvas *c27 = new TCanvas("c27","delta(dy and dx) (misaligned and original)");
  c27->Divide(2,2);
  c27->cd(1);
  deltady->GetYaxis()->SetRangeUser(0.1,5e4);
  gPad->SetLogy();
  deltady->Draw();
  c27->cd(2);
  deltadx->GetYaxis()->SetRangeUser(0.1,5e4);
  gPad->SetLogy();
  deltadx->Draw();
  c27->cd(3);
  deltady2->GetYaxis()->SetRangeUser(0.1,5e4);
  gPad->SetLogy();
  deltady2->Draw();
  c27->cd(4);
  deltadx2->GetYaxis()->SetRangeUser(0.1,5e4);
  gPad->SetLogy();
  deltadx2->Draw();


  TCanvas *c28 = new TCanvas ("c28","momenta");
  c28->cd();
  TH1F *hcolmom[2] = {colorhists(hmom,1),colorhists(hmom2,8)};
  hcolmom[0]->Draw();
  hcolmom[1]->Draw("same");
std::string legnameM[2] = {

  "misaligned",

  " original"};
  TLegend *legmomM = (TLegend*)legshape2(2, hcolmom,legnameM);
  legmomM->Draw("same");

  TCanvas *c29 = new TCanvas ("c29","ty vs y");
  c29 = new TCanvas("c29","ty vs y");
  c29->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c29->cd(ci+1);

      const int nhist = 2;//6;
 
      c29->cd(1);
      //tyvsy[0]->Draw("colz");

      TH1F *hcol[nhist] = {fitslicesAng(tyvsy[ci], 1),
         fitslicesAng(tyvsy2[ci], 8)
      };

  
      std::string legname[nhist] = {

  "misaligned",
  "original"
  //"trkqual,mupull B>0",
  //"trkqual,mupull B=0"
      };

      hcol[0]->Draw();
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      leg0[ci]->Draw("same");

    }

  c29->Update();
  c29->Print("tpcfgd_tyVSy_"+oFile+".eps");

  TCanvas *c29a = new TCanvas ("c29a","ty vs y(2D)");
  c29a = new TCanvas("c29a","ty vs y");
  c29a->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c29a->cd(ci+1);
      tyvsy[ci]->Draw("colz");
    }
  c29a->Update();
  c29a->Print("tpcfgd_tyVSy_2D_"+oFile+".eps");

  TCanvas *c29p = new TCanvas ("c29p","Primary: ty vs y(2D)");
  c29a = new TCanvas("c29a","ty vs y");
  c29p->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c29p->cd(ci+1);
      tyvsyPrim[ci]->Draw("colz");
    }
  c29p->Update();
  c29p->Print("tpcfgd_tyVSy_2D_Primary"+oFile+".eps");


  TCanvas *c29b = new TCanvas ("c29b","All ty vs y(2D)");
  c29a = new TCanvas("c29a","ty vs y");
  c29b->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c29b->cd(ci+1);
      tyvsyA[ci]->Draw("colz");
    }  

  TCanvas *c30 = new TCanvas ("c30","tx vs y");
  c30 = new TCanvas("c30","tx vs y");
  c30->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c30->cd(ci+1);

      const int nhist = 2;//6;
 
      TH1F *hcol[nhist] = {fitslicesAng(txvsy[ci], 1),
         fitslicesAng(txvsy2[ci], 8)
      };

  
      std::string legname[nhist] = {

  "misaligned",
  "original"
  //"trkqual,mupull B>0",
  //"trkqual,mupull B=0"
      };

      hcol[0]->Draw();
      for(Int_t hi=0; hi<nhist; hi++)
  {
    hcol[hi]->Draw("same");
  }
 
      leg0.push_back((TLegend*)legshape1(nhist, hcol,legname));
      leg0[ci]->Draw("same");

    }

  c30->Update();
  c30->Print("tpcfgd_txVSy_"+oFile+".eps");


  TCanvas *c31 = new TCanvas("c31","ty per det boundary");
  vector <TLegend* > leg31;
 c31->Divide(2,2);
  for(Int_t ci=0; ci<4; ci++)
    {
      c31->cd(ci+1);
  
      const int nhist = 2;//6;
      TH1F *hcol42[nhist] = {colorhists(tyTotDet[ci], 1),colorhists(tyTotDet2[ci], 8)};
      hcol42[1]->Draw();
      hcol42[0]->Draw("same");


      std::string legname[nhist] = {
  "misaligned",
  "original"
      };

      leg31.push_back((TLegend*)legshape2(nhist, hcol42,legname));
      leg31[ci]->Draw("same");

    }
  c31->Update();
  c31->Print("tpcfgd_typerdet_"+oFile+".eps");

  out->Write();
  out->Close();

  oconstfl<<"\n";
  oconstfl.close();


}