//_____________________________________________________________________________
// XFT Analysis Module
//_____________________________________________________________________________


#include <iostream>
#include "TF1.h"
#include "Helix.hh"
#include "TCanvas.h"
#include "Stntuple/obj/TStnHeaderBlock.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"

#include "TXftAnaModule.hh"

float konst = 0.00014989*14.116;
float halfpi = 0.5*3.14159;
float twopi = 2.0*3.14159;
float speedLight = 29.98;     // Speed of light in cm/nsec

int minCell[] = {999,999,999,999};
int maxCell[] = {-1,-1,-1,-1};
float cotRadius[8][12];
float cotPhi[8][12];
// wire variables...
int nsn = 12;
//	float dwlr = 0.5*0.762;// half distance between  wires
float dwlr = 0.5*0.7112;
// Number of pixels
int npix_cell[] = {6,12,6,12,6,6,6,6};  
// Number of cells
/* 168, 192 , 240,  288  , 336,  384,  432,   480*/  
int ncl[] = {42*4, 48*4, 60*4, 72*4, 84*4, 96*4, 108*4, 120*4};
float rsns[] = {46.0, 58.534, 70.0, 82.055, 94.0, 105.575, 119.0, 129.096};
//	float[] rsns = {46.0, 58.0, 70.0, 82.0, 94.0, 106.0, 119.0, 131.0};
/*  .0374,             .0327        .0262          .0218*/
float dclw[] = {twopi/ncl[0], 
		twopi/ncl[1],  
		twopi/ncl[2], 
		twopi/ncl[3],  
		twopi/ncl[4], 
		twopi/ncl[5],  
		twopi/ncl[6], 
		twopi/ncl[7]};
float pixPhi[] = {twopi/(ncl[0]*npix_cell[0]),  twopi/(ncl[1]*npix_cell[1]),  
		  twopi/(ncl[2]*npix_cell[2]),  twopi/(ncl[3]*npix_cell[3]),
		  twopi/(ncl[4]*npix_cell[4]),  twopi/(ncl[5]*npix_cell[5]),  
		  twopi/(ncl[6]*npix_cell[6]),  twopi/(ncl[7]*npix_cell[7])};


float ang = 35.0 * twopi / 360.0; /* .61 */
float drift = 0.0088;	    // drift = 100 microns/nsec
float phz0 = 0.0;
// These variables set constants for the XFT
float wireEfficiency = 1.00;
float promptTime = 44.0;
float delayedTime = 88.0;
float cutoffTime = 132.0;
int prtlvl = 0;

int xft_Total = 0;
int xft_unmatched = 0;
int index_match = 0;
int nfid = 0;
int nfid_matched = 0;
int ntrigger =0;
int ntrigger_matched =0;

// ClassImp(TXftAnaModule)
//_____________________________________________________________________________
TXftAnaModule::TXftAnaModule(const char* name, const char* title):
  TStnModule(name,title)
{
}

//_____________________________________________________________________________
TXftAnaModule::~TXftAnaModule() {
}

//_____________________________________________________________________________
void TXftAnaModule::SetJobName(const char* jobname) {
  jobName = jobname;
}
//_____________________________________________________________________________
void TXftAnaModule::BookHistograms() {

  printf("Job Name: %s\n", jobName); 
  //
  char filename[120];
  strcat(strcpy(filename,jobName),".hbk");
  printf("Histogram output file: %s\n", filename); 
  //
  // Make the histogram output file name
  char mess[120];
  strcat(strcpy(mess,jobName)," file analyzed with xftRoot.C");
  printf("Histogram output file root message: %s\n", mess);
  outfile = new TFile(filename,"RECREATE",mess);
  //
  // Delete the old histos
  GetListOfHistograms()->Delete();
  //
  // Now define histos
  
  // XFT Pixels 
  xftNHits = new TH1F("xftNHits","XFT Number of Hits",1000,0.,10000.);
  xftNPixels = new TH1F("xftNPixels","XFT Number of Pixels",1000,0,5000);
  xftNTracks = new TH1F("xftNTracks","XFT Number of 4-layer tracks",300,0.,300.);
  offNTracks = new TH1F("offNTracks","Offline Number of tracks",1000,0.,1000.);
  
  
  SL1PixelChip = new TH1F("SL1PixelChip","SL1 Pixels per Chip",48,0.,6.2832);
  SL2PixelChip = new TH1F("SL2PixelChip","SL2 Pixels per Chip",72,0.,6.2832);
  SL3PixelChip = new TH1F("SL3PixelChip","SL3 Pixels per Chip",96,0.,6.2832);
  SL4PixelChip = new TH1F("SL4PixelChip","SL4 Pixels per Chip",120,0.,6.2832);
  
  SL1PixelCell = new TH1F("SL1PixelCell","SL1 Pixels per Cell",192,0.,6.2832);
  SL2PixelCell = new TH1F("SL2PixelCell","SL2 Pixels per Cell",288,0.,6.2832);
  SL3PixelCell = new TH1F("SL3PixelCell","SL3 Pixels per Cell",384,0.,6.2832);
  SL4PixelCell = new TH1F("SL4PixelCell","SL4 Pixels per Cell",480,0.,6.2832);
  
  // All XFT tracks
  xftPtBinShort = new TH1F("xftPtBinShort","XFT Pt Bin Short Tracks",125,0.,124.);
  xftPtBin = new TH1F("xftPtBin","XFT Pt Bin",125,0.,124.);
  xftPtVal = new TH1F("xftPtVal","XFT Pt Val",1000,-200.,200.);
  xftShort = new TH1F("xftShort","XFT Short Bit",2,0.,2.);
  xftIso = new TH1F("xftIso","XFT Iso",2,0.,2.);
  xftAbsPhi = new TH1F("xftAbsPhi","XFT AbsolutePhi",288,0.0,6.2832);
  xftCrate = new TH1F("xftCrate","XFT Crate",3,0.,3.);
  xftBoard = new TH1F("xftBoard","XFT Board",24,0.,24.);
  xftChip = new TH1F("xftChip","XFT Chip",288,0.,288.);
  xftPhiBin = new TH1F("xftPhiBin","XFT Phi Bin",2304,0.,288.);
  //
  //Trigger studies
  //
  xftAbsPhi4 = new TH1F("xftAbsPhi4","XFT AbsolutePhi above 4 GeV/c",288,0.0,6.2832);
  xftAbsPhi8 = new TH1F("xftAbsPhi8","XFT AbsolutePhi above 8 GeV/c",288,0.0,6.2832);
  xftpt_cem_match = new TH1F("xftpt_cem_match","XFT Pt - cem 8 pt 8",125,0.,124.);
  xftphi_cem_match = new TH1F("xftphi_cem_match","XFT AbsolutePhi - cem 8 pt 8",288,0.0,6.2832);
  xftpt_cem_nomatch4 = new TH1F("xftpt_cem_nomatch4","XFT Pt - cem 0 pt 4",125,0.,124.);
  xftphi_cem_nomatch4 = new TH1F("xftphi_cem_nomatch4","XFT AbsolutePhi - cem 0 pt 4",288,0.0,6.2832);
  //
  xftpt_cmu_match = new TH1F("xftpt_cmu_match","XFT Pt - cmu 6 pt 8",125,0.,124.);
  xftphi_cmu_match = new TH1F("xftphi_cmu_match","XFT AbsolutePhi - cmu 6 pt 8",288,0.0,6.2832);
  xftpt_cmu_nomatch8 = new TH1F("xftpt_cmu_nomatch8","XFT Pt - cmu 0 pt 8",125,0.,124.);
  xftphi_cmu_nomatch8 = new TH1F("xftphi_cmu_nomatch8","XFT AbsolutePhi - cmu 0 pt 8",288,0.0,6.2832);
  //
  //Want one with jet trigger
  xftpt_jet = new TH1F("xftpt_jet","XFT Pt - jet trigger",125,0.,124.);
  xftphi_jet = new TH1F("xftphi_jet","XFT AbsolutePhi - jet trigger",288,0.0,6.2832);
  
  
  //Fake XFT tracks
  xftFakePtBin = new TH1F("xftFakePtBin","XFT UnmatchedPt Bin",125,0.,124.);
  xftFakePtVal = new TH1F("xftFakePtVal","XFT UnmatchedPt Val",1000,-200.,200.);
  xftFakeShort = new TH1F("xftFakeShort","XFT UnmatchedShort Bit",2,0.,2.);
  xftFakeIso = new TH1F("xftFakeIso","XFT UnmatchedIso",2,0.,2.);
  xftFakeAbsPhi = new TH1F("xftFakeAbsPhi","XFT UnmatchedAbsolutePhi",288,0.0,6.2832);
  xftFakeCrate = new TH1F("xftFakeCrate","XFT UnmatchedCrate",3,0.,3.);
  xftFakeBoard = new TH1F("xftFakeBoard","XFT UnmatchedBoard",24,0.,24.);
  xftFakeChip = new TH1F("xftFakeChip","XFT UnmatchedChip",288,0,288);
  xftFakePhiBin = new TH1F("xftFakePhiBin","XFT FakePhiBin",2304,0.,288);
  
  // All OFFLINE tracks
  offPhi = new TH1F("offPhi","Offline Phi",288,0.0,6.2832);
  offd0 = new TH1F("offd0","Offline d0",50,-2.0,2.0);
  offphid0 = new TH2F("offphid0","Offline phi vs d0",50,0.0,6.2832,50,-2.0,2.0);
  offz0 = new TH1F("offz0","Offline z0",50,-50.0,50.0);
  offPt = new TH1F("offPt","Offline Pt",2000,0.0,200.0);
  offCotan = new TH1F("offCotan","Offline cotan",50,-2.0, 2.0);
  offEta = new TH1F("offEta","Offline eta",50,-2.0, 2.0);
  CotHitsTot = new TH1F("CotHitsTot","Offline Total COT Hits",100,0,100);
  CotHitsAx = new TH1F("CotHitsAx","Offline Axial COT Hits",50,0,50);
  CotHitsSt = new TH1F("CotHitsSt","Offline Stereo COT Hits",50,0,50);
  offZSl3 = new TH1F("offZSl3","Offline Z(ax sl3)",150,-150,150);
  offZSl2 = new TH1F("offZSl2","Offline Z(ax sl2)",150,-150,150);
  offZSl1 = new TH1F("offZSl1","Offline Z(ax sl1)",150,-150,150);
  offZSl0 = new TH1F("offZSl0","Offline Z(ax sl0)",150,-150,150);
  CotHitsTotZSl3 = new TH2F("CotHitsTotZSl3","Offline Total COT Hits versus ZSl3",150,-150,150,100,0,100);
  CotHitsAxZSl3 = new TH2F("CotHitsAxZSl3","Offline Axial COT Hits versus  ZSl3",150,-150,150,50,0,50);
  CotHitsStZSl3 = new TH2F("CotHitsStZSl3","Offline Stereo COT Hits versus ZSl3",150,-150,150,50,0,50);
  
  //No fiducial requirement
  allPhi = new TH1F("allPhi","all Phi",288,0.0,6.2832);
  alld0 = new TH1F("alld0","all d0",50,-2.0,2.0);
  allz0 = new TH1F("allz0","all z0",50,-50.0,50.0);
  allPt = new TH1F("allPt","all Pt",2000,0.0,200.0);
  allCotan = new TH1F("allCotan","all cotan",50,-2.0, 2.0);
  allEta = new TH1F("allEta","all eta",50,-2.0, 2.0);
  allZSl3 = new TH1F("allZSl3","all Z(ax sl3)",60,-150,150);
  allZSl2 = new TH1F("allZSl2","all Z(ax sl2)",60,-150,150);
  allZSl1 = new TH1F("allZSl1","all Z(ax sl1)",60,-150,150);
  allZSl0 = new TH1F("allZSl0","all Z(ax sl0)",60,-150,150);
  qallPhi = new TH1F("qallPhi","qall Phi",288,0.0,6.2832);
  qalld0 = new TH1F("qalld0","qall d0",50,-2.0,2.0);
  qallz0 = new TH1F("qallz0","qall z0",50,-50.0,50.0);
  qallPt = new TH1F("qallPt","qall Pt",2000,0.0,200.0);
  qallCotan = new TH1F("qallCotan","qall cotan",50,-2.0, 2.0);
  qallEta = new TH1F("qallEta","qall eta",50,-2.0, 2.0);
  qallZSl3 = new TH1F("qallZSl3","qall Z(ax sl3)",150,-150,150);
  qallZSl2 = new TH1F("qallZSl2","qall Z(ax sl2)",150,-150,150);
  qallZSl1 = new TH1F("qallZSl1","qall Z(ax sl1)",150,-150,150);
  qallZSl0 = new TH1F("qallZSl0","qall Z(ax sl0)",150,-150,150);
  
  // OFFLINE tracks found by XFT 
  offPhiFound = new TH1F("offPhiFound","Offline XFT Matched Phi",288,0.0,6.2832);
  offd0Found = new TH1F("offd0Found","Offline XFT Matched d0",50,-2.0,2.0);
  offz0Found = new TH1F("offz0Found","Offline XFT Matched z0",50,-50.0,50.0);
  offPtFound = new TH1F("offPtFound","Offline XFT Matched Pt",2000,0.0,200.0);
  offCotanFound = new TH1F("offCotanFound","Offline XFT Matched cotan",50,-2.0, 2.0);
  offEtaFound = new TH1F("offEtaFound","Offline XFT Matched eta",50,-2.0, 2.0);
  CotHitsTotFound = new TH1F("CotHitsTotFound","Offline Total COT Hits, Found XFT",100,0,100);
  CotHitsAxFound = new TH1F("CotHitsAxFound","Offline Axial COT Hits, Found XFT",50,0,50);
  CotHitsStFound = new TH1F("CotHitsStFound","Offline Stereo COT Hits, Found XFT",50,0,50);
  offZSl3Found = new TH1F("offZSl3Found","Offline XFT Matched Z(ax sl3)",150,-150,150);
  offZSl2Found = new TH1F("offZSl2Found","Offline XFT Matched Z(ax sl2)",150,-150,150);
  offZSl1Found = new TH1F("offZSl1Found","Offline XFT Matched Z(ax sl1)",150,-150,150);
  offZSl0Found = new TH1F("offZSl0Found","Offline XFT Matched Z(ax sl0)",150,-150,150);
  CotHitsTotZSl3Found = new TH2F("CotHitsTotZSl3Found","Offline Matched Total COT Hits versus ZSl3",150,-150,150,100,0,100);
  CotHitsAxZSl3Found = new TH2F("CotHitsAxZSl3Found","Offline Matched Axial COT Hits versus ZSl3",150,-150,150,50,0,50);
  CotHitsStZSl3Found = new TH2F("CotHitsStZSl3Found","Offline Matched Stereo COT Hits versus ZSl3",150,-150,150,50,0,50);
  
  // OFFLINE tracks missed by XFT
  offPhiMissed = new TH1F("offPhiMissed","Offline XFT Missed Phi",288,0.0,6.2832);
  offd0Missed = new TH1F("offd0Missed","Offline XFT Missed d0",50,-2.0,2.0);
  offz0Missed = new TH1F("offz0Missed","Offline XFT Missed z0",50,-50.0,50.0);
  offPtMissed = new TH1F("offPtMissed","Offline XFT Missed Pt",2000,0.0,200.0);
  offCotanMissed = new TH1F("offCotanMissed","Offline XFT Missed cotan",50,-2.0, 2.0);
  offEtaMissed = new TH1F("offEtaMissed","Offline XFT Missed eta",50,-2.0, 2.0);
  CotHitsTotMissed = new TH1F("CotHitsTotMissed","Offline Total COT Hits, Missed XFT",100,0,100);
  CotHitsAxMissed = new TH1F("CotHitsAxMissed","Offline Axial COT Hits, Missed XFT",50,0,50);
  CotHitsStMissed = new TH1F("CotHitsStMissed","Offline Stereo COT Hits, Missed XFT",50,0,50);
  offZSl3Missed = new TH1F("offZSl3Missed","Offline XFT Missed Z(ax sl3)",150,-150,150);
  offZSl2Missed = new TH1F("offZSl2Missed","Offline XFT Missed Z(ax sl2)",150,-150,150);
  offZSl1Missed = new TH1F("offZSl1Missed","Offline XFT Missed Z(ax sl1)",150,-150,150);
  offZSl0Missed = new TH1F("offZSl0Missed","Offline XFT Missed Z(ax sl0)",150,-150,150);
  CotHitsTotZSl3Missed = new TH2F("CotHitsTotZSl3Missed","Offline XFT Missed Total COT Hits versus ZSl3",150,-150,150,100,0,100);
  CotHitsAxZSl3Missed = new TH2F("CotHitsAxZSl3Missed","Offline XFT Missed Axial COT Hits versus ZSl3",150,-150,150,50,0,50);
  CotHitsStZSl3Missed = new TH2F("CotHitsStZSl3Missed","Offline XFT Missed Stereo COT Hits versus ZSl3",150,-150,150,50,0,50);
  
  //Efficiency = (Offline tracks found by XFT)/ (All Offline tracks)
  offPhiEff = new TH1F("offPhiEff","Offline XFT Eff Phi",288,0.0,6.2832);
  offd0Eff = new TH1F("offd0Eff","Offline XFT Eff d0",50,-2.0,2.0);
  offz0Eff = new TH1F("offz0Eff","Offline XFT Eff z0",50,-50.0,50.0);
  offPtEff = new TH1F("offPtEff","Offline XFT Eff Pt",2000,0.0,200.0);
  offCotanEff = new TH1F("offCotanEff","Offline XFT Eff cotan",50,-2.0, 2.0);
  offEtaEff = new TH1F("offEtaEff","Offline XFT Eff eta",50,-2.0, 2.0);
  CotHitsTotEff = new TH1F("CotHitsTotEff","Offline Total COT Hits",100,0,100);
  CotHitsAxEff = new TH1F("CotHitsAxEff","Offline Axial COT Hits",50,0,50);
  CotHitsStEff = new TH1F("CotHitsStEff","Offline Stereo COT Hits",50,0,50);
  offZSl3Eff = new TH1F("offZSl3Eff","Offline XFT Eff Z(ax sl3)",150,-150,150);
  offZSl2Eff = new TH1F("offZSl2Eff","Offline XFT Eff Z(ax sl2)",150,-150,150);
  offZSl1Eff = new TH1F("offZSl1Eff","Offline XFT Eff Z(ax sl1)",150,-150,150);
  offZSl0Eff = new TH1F("offZSl0Eff","Offline XFT Eff Z(ax sl0)",150,-150,150);
	
  // Offline pt spectra 
  offPtTriggerAll1p5 = new TH1F("offPtTriggerAll1p5","Offline Pt - all",1000,0.0,20.0);
  offPtTriggerAll4 = new TH1F("offPtTriggerAll4","Offline Pt - all",200,0.0,20.0);
  offPtTriggerAll8 = new TH1F("offPtTriggerAll8","Offline Pt - all",200,0.0,20.0);
  
  offPtTrigger1p5 = new TH1F("offPtTrigger1p5","Offline Pt - XFT matched ",1000,0.0,20.0);
  offPtTrigger4 = new TH1F("offPtTrigger4","Offline Pt above 4 GeV/c - XFT matched ",200,0.0,20.0);
  offPtTrigger8 = new TH1F("offPtTrigger8","Offline Pt above 8 GeV/c- XFT matched ",200,0.0,20.0);
  
  offInvPtTriggerAll = new TH1F("offInvPtTriggerAll","Offline 1/Pt - all",100,0.0,1.0);
  offInvPtTrigger1p5 = new TH1F("offInvPtTrigger1p5","Offline 1/Pt - XFT matched ",100,0.0,1.0);
  offInvPtTrigger4 = new TH1F("offInvPtTrigger4","Offline 1/Pt above 4 GeV/c - XFT matched ",100,0.0,1.0);
  offInvPtTrigger8 = new TH1F("offInvPtTrigger8","Offline 1/Pt above  GeV/c - XFT matched ",100,0.0,1.0);
  
  // Offline pt spectra for negatively charged tracks
  moffPtTriggerAll1p5 = new TH1F("moffPtTriggerAll1p5","-ve Offline Pt - all",1000,0.0,20.0);
  moffPtTriggerAll4 = new TH1F("moffPtTriggerAll4","-ve Offline Pt above 4 GeV/c - all",200,0.0,20.0);
  moffPtTriggerAll8 = new TH1F("moffPtTriggerAll8","-ve Offline Pt above 4 GeV/c - all",200,0.0,20.0);
  moffPtTrigger1p5 = new TH1F("moffPtTrigger1p5","-ve Offline Pt - XFT matched",1000,0.0,20.0);
  moffPtTrigger4 = new TH1F("moffPtTrigger4","-ve Offline Pt above 4 GeV/c - XFT matched",200,0.0,20.0);
  moffPtTrigger8 = new TH1F("moffPtTrigger8","-ve Offline Pt above 8 GeV/c- XFT matched ",200,0.0,20.0);
  
  // Offline pt spectra for positively charged tracks
  poffPtTriggerAll1p5 = new TH1F("poffPtTriggerAll1p5","+ve Offline Pt - all",1000,0.0,20.0);
  poffPtTriggerAll4 = new TH1F("poffPtTriggerAll4","+ve Offline Pt above 4 GeV/c - all",200,0.0,20.0);
  poffPtTriggerAll8 = new TH1F("poffPtTriggerAll8","+ve Offline Pt above 4 GeV/c - all",200,0.0,20.0);
  poffPtTrigger1p5 = new TH1F("poffPtTrigger1p5","+ve Offline Pt - XFT matched",1000,0.0,20.0);
  poffPtTrigger4 = new TH1F("poffPtTrigger4","+ve Offline Pt above 4 GeV/c - XFT matched",200,0.0,20.0);
  poffPtTrigger8 = new TH1F("poffPtTrigger8","+ve Offline Pt above 8 GeV/c- XFT matched ",200,0.0,20.0);
  
  // Turn-on curves
  turnon1p5 = new TH1F("turnon1p5","Turn on at 1.5 GeV/c",1000,0.0,20.0);
  turnon4 = new TH1F("turnon4","Turn on at 4 GeV/c ",200,0.0,20.0);
  turnon8 = new TH1F("turnon8","Turn on at 8 GeV/c ",200,0.0,20.0);
  
  iturnon1p5 = new TH1F("iturnon1p5","Turn on (1/pT) at 1.5 GeV/c",100,0.0,1.0);
  iturnon4 = new TH1F("iturnon4","Turn on (1/pT) at 8 GeV/c ",100,0.0,1.0);
  iturnon8 = new TH1F("iturnon8","Turn on (1/pT) at 8 GeV/c ",100,0.0,1.0);
  
  // Turn-on curves for negatively charged tracks
  mturnon1p5 = new TH1F("mturnon1p5","-ve Turn on at 1.5 GeV/c",1000,0.0,20.0);
  mturnon4 = new TH1F("mturnon4","-ve Turn on above at 4 GeV/c ",200,0.0,20.0);
  mturnon8 = new TH1F("mturnon8","-ve Turn on above at 8 GeV/c ",200,0.0,20.0);
  
  // Turn-on curves for positively charged tracks
  pturnon1p5 = new TH1F("pturnon1p5","+ve Turn on at 1.5 GeV/c",1000,0.0,20.0);
  pturnon4 = new TH1F("pturnon4","+ve Turn on at 4 GeV/c",200,0.0,20.0);
  pturnon8 = new TH1F("pturnon8","+ve Turn on at 8 Gev/c",200,0.0,20.0);
  
  
  // Pt, phi resolutions for all, positive, negative tracks
  dPhiHist = new TH1F("dPhiHist","DeltaPhi(XFT-Offline)",60,-0.06,0.06);
  dPtHist = new TH1F("dPtHist","DeltaPt(XFT-Offline)/Pt**2",50,-0.2,0.2);
  
  PdPhiHist = new TH1F("PdPhiHist","DeltaPhi(XFT-Offline) positives",100,-0.1,0.1);
  PdPtHist = new TH1F("PdPtHist","DeltaPt(XFT-Offline)/Pt**2  positives",50,-0.2,0.2);
  
  MdPhiHist = new TH1F("MdPhiHist","DeltaPhi(XFT-Offline) negatives",100,-0.1,0.1);
  MdPtHist = new TH1F("MdPtHist","DeltaPt(XFT-Offline)/Pt**2 negatives",50,-0.2,0.2);
  //
  // Pt, phi resolutions vs Phi
  dPhiPhiHist = new TH2F("dPhiPhiHist","DeltaPhi(XFT-Offline) vs Phi",25,0.0,6.2832,25,-0.1,0.1);
  dPtPhiHist = new TH2F("dPtPhiHist","DeltaPt(XFT-Offline)/Pt**2 vs Phi",25,0.0,6.2832,25,-0.2,0.2);
  
  // Pt bin offline vs Pt bin XFT
  PtBinPtBinHist = new TH2F("PtBinPtBinHist","Pt bin Off vs Pt Bin XFT ",97,0.0,96.0,97,0.0,96.0);
  
  // Pt bin difference vs chip, phi for all, +ve, -ve tracks
  dPtBinChipHist = new TH2F("dPtBinChipHist","DeltaPtBin vs Chip",288,0,288,25,-5.0,5.0);
  dPtBinPhiHist = new TH2F("dPtBinPhiHist","DeltaPtBin vs Phi",25,0.0,6.2832,25,-5.0,5.0);
  dPtBinPhiPos = new TH2F("dPtBinPhiPos","DeltaPtBin Pos vs Phi",25,0.0,6.2832,25,-5.0,5.0);
  dPtBinPhiNeg = new TH2F("dPtBinPhiNeg","DeltaPtBin Neg vs Phi",25,0.0,6.2832,25,-5.0,5.0);
  //
  // Corrected pt resolution and pt resolution vs phi for +ve and -ve tracks
  dPtNewHist = new TH1F("dPtNewHist","DeltaPt(XFT-Offline)/Pt**2 new",25,-0.2,0.2);
  dPtBinPhiPosCorr = new TH2F("dPtBinPhiPosCorr","Corr DeltaPt(bin) Pos vs Phi",25,0.0,6.2832,25,-5.0,5.0);
  dPtBinPhiNegCorr = new TH2F("dPtBinPhiNegCorr","Corr DeltaPt(bin) Neg vs Phi",25,0.0,6.2832,25,-5.0,5.0);
  
  // Angle between pairs of tracks for all, opposite, same sign pairs and those removed as duplicates a la XTRP
  PairPhi = new TH1F("PairPhi","Pair Phi",180,0.0,180.);
  PairPhiOppo = new TH1F("PairPhiOppo","Pair Phi opposite sign",180,0.0,180.);
  PairPhiSame = new TH1F("PairPhiSame","Pair Phi same sign",180,0.0,180.);
  
  duplPairPhi = new TH1F("duplPairPhi","Duplicate Pair Phi",1800,0.0,180.);
  duplPairPhiOppo = new TH1F("duplPairPhiOppo","Duplicate Pair Phi opposite sign",1800,0.0,180.);
  duplPairPhiSame = new TH1F("duplPairPhiSame","Duplicate Pair Phi same sign",1800,0.0,180.);
  
  //XTC hits
  cellsPrompt0 = new TH1F("cellsPrompt0","Cells with Prompt Hits, Ax0",192,0,192);
  cellsPrompt1 = new TH1F("cellsPrompt1","Cells with Prompt Hits, Ax1",288,0,288);
  cellsPrompt2 = new TH1F("cellsPrompt2","Cells with Prompt Hits, Ax2",384,0,384);
  cellsPrompt3 = new TH1F("cellsPrompt3","Cells with Prompt Hits, Ax3",480,0,480);
  
  cellsDelayed0 = new TH1F("cellsDelayed0","Cells with Delayed Hits, Ax0",192,0,192);
  cellsDelayed1 = new TH1F("cellsDelayed1","Cells with Delayed Hits, Ax1",288,0,288);
  cellsDelayed2 = new TH1F("cellsDelayed2","Cells with Delayed Hits, Ax2",384,0,384);
  cellsDelayed3 = new TH1F("cellsDelayed3","Cells with Delayed Hits, Ax3",480,0,480);
  
  cellsPD0 = new TH1F("cellsPD0","Cells with PD Hits, Ax0",192,0,192);
  cellsPD1 = new TH1F("cellsPD1","Cells with PD Hits, Ax1",288,0,288);
  cellsPD2 = new TH1F("cellsPD2","Cells with PD Hits, Ax2",384,0,384);
  cellsPD3 = new TH1F("cellsPD3","Cells with PD Hits, Ax3",480,0,480);
  
  wiresDelayed0 = new TH1F("wiresDelayed0","wires with Delayed Hits, Ax0",192*12,0,192*12);
  wiresDelayed1 = new TH1F("wiresDelayed1","wires with Delayed Hits, Ax1",288*12,0,288*12);
  wiresDelayed2 = new TH1F("wiresDelayed2","wires with Delayed Hits, Ax2",384*12,0,384*12);
  wiresDelayed3 = new TH1F("wiresDelayed3","wires with Delayed Hits, Ax3",480*12,0,480*12);
  
  wiresPrompt0 = new TH1F("wiresPrompt0","wires with Prompt Hits, Ax0",192*12,0,192*12);
  wiresPrompt1 = new TH1F("wiresPrompt1","wires with Prompt Hits, Ax1",288*12,0,288*12);
  wiresPrompt2 = new TH1F("wiresPrompt2","wires with Prompt Hits, Ax2",384*12,0,384*12);
  wiresPrompt3 = new TH1F("wiresPrompt3","wires with Prompt Hits, Ax3",480*12,0,480*12);
  
  wiresPD0 = new TH1F("wiresPD0","wires with PD Hits, Ax0",192*12,0,192*12);
  wiresPD1 = new TH1F("wiresPD1","wires with PD Hits, Ax1",288*12,0,288*12);
  wiresPD2 = new TH1F("wiresPD2","wires with PD Hits, Ax2",384*12,0,384*12);
  wiresPD3 = new TH1F("wiresPD3","wires with PD Hits, Ax3",480*12,0,480*12);
  
  //L3 tracks, extrapolate and note hits used in track fit.
  wiresMissedTDC0 = new TH1F("wiresMissedTDC0","wires with Missed TDC Hits, Ax0",192*12,0,192*12);
  wiresMissedTDC1 = new TH1F("wiresMissedTDC1","wires with Missed TDC Hits, Ax1",288*12,0,288*12);
  wiresMissedTDC2 = new TH1F("wiresMissedTDC2","wires with Missed TDC Hits, Ax2",384*12,0,384*12);
  wiresMissedTDC3 = new TH1F("wiresMissedTDC3","wires with Missed TDC Hits, Ax3",480*12,0,480*12);
  
  wiresFoundTDC0 = new TH1F("wiresFoundTDC0","wires with Found TDC Hits, Ax0",192*12,0,192*12);
  wiresFoundTDC1 = new TH1F("wiresFoundTDC1","wires with Found TDC Hits, Ax1",288*12,0,288*12);
  wiresFoundTDC2 = new TH1F("wiresFoundTDC2","wires with Found TDC Hits, Ax2",384*12,0,384*12);
  wiresFoundTDC3 = new TH1F("wiresFoundTDC3","wires with Found TDC Hits, Ax3",480*12,0,480*12);
  
  //L3 tracks, extrapolate and try to match to XTC hits
  cellsFound0 = new TH1F("cellsFound0","Cells with Found Hits, Ax0",192,0,192);
  cellsFound1 = new TH1F("cellsFound1","Cells with Found Hits, Ax1",288,0,288);
  cellsFound2 = new TH1F("cellsFound2","Cells with Found Hits, Ax2",384,0,384);
  cellsFound3 = new TH1F("cellsFound3","Cells with Found Hits, Ax3",480,0,480);
  
  cellsMissed0 = new TH1F("cellsMissed0","Cells with Missed Hits, Ax0",192,0,192);
  cellsMissed1 = new TH1F("cellsMissed1","Cells with Missed Hits, Ax1",288,0,288);
  cellsMissed2 = new TH1F("cellsMissed2","Cells with Missed Hits, Ax2",384,0,384);
  cellsMissed3 = new TH1F("cellsMissed3","Cells with Missed Hits, Ax3",480,0,480);
  
  
  wiresMissed0 = new TH1F("wiresMissed0","wires with Missed Hits, Ax0",192*12,0,192*12);
  wiresMissed1 = new TH1F("wiresMissed1","wires with Missed Hits, Ax1",288*12,0,288*12);
  wiresMissed2 = new TH1F("wiresMissed2","wires with Missed Hits, Ax2",384*12,0,384*12);
  wiresMissed3 = new TH1F("wiresMissed3","wires with Missed Hits, Ax3",480*12,0,480*12);
  
  
  HitDcell0 = new TH1F("HitDcell0","Offline-Hit Delta cell, Ax0",20,-10,10);
  HitDcell1 = new TH1F("HitDcell1","Offline-Hit Delta cell, Ax1",20,-10,10);
  HitDcell2 = new TH1F("HitDcell2","Offline-Hit Delta cell, Ax2",20,-10,10);
  HitDcell3 = new TH1F("HitDcell3","Offline-Hit Delta cell, Ax3",20,-10,10);
  
  foundHits0 = new TH2F("foundHits0","Hits near track, Ax0",192,0,192,13,0,13);
  foundHits1 = new TH2F("foundHits1","Hits near track, Ax1",288,0,288,13,0,13);
  foundHits2 = new TH2F("foundHits2","Hits near track, Ax2",384,0,384,13,0,13);
  foundHits3 = new TH2F("foundHits3","Hits near track, Ax3",480,0,480,13,0,13);
  
  foundHitsIso0 = new TH2F("foundHitsIso0","HitsIso near track, Ax0",192,0,192,13,0,13);
  foundHitsIso1 = new TH2F("foundHitsIso1","HitsIso near track, Ax1",288,0,288,13,0,13);
  foundHitsIso2 = new TH2F("foundHitsIso2","HitsIso near track, Ax2",384,0,384,13,0,13);
  foundHitsIso3 = new TH2F("foundHitsIso3","HitsIso near track, Ax3",480,0,480,13,0,13);
  
  foundHitsIsoUM0 = new TH2F("foundHitsIsoUM0","HitsIsoUM near track, Ax0",192,0,192,13,0,13);
  foundHitsIsoUM1 = new TH2F("foundHitsIsoUM1","HitsIsoUM near track, Ax1",288,0,288,13,0,13);
  foundHitsIsoUM2 = new TH2F("foundHitsIsoUM2","HitsIsoUM near track, Ax2",384,0,384,13,0,13);
  foundHitsIsoUM3 = new TH2F("foundHitsIsoUM3","HitsIsoUM near track, Ax3",480,0,480,13,0,13);
  
  foundHitsZSl0 = new TH2F("foundHitsZSl0","Hits near track versus z, Ax0",150,-150,150,13,0,13);
  foundHitsZSl1 = new TH2F("foundHitsZSl1","Hits near track versus z, Ax1",150,-150,150,13,0,13);
  foundHitsZSl2 = new TH2F("foundHitsZSl2","Hits near track versus z, Ax2",150,-150,150,13,0,13);
  foundHitsZSl3 = new TH2F("foundHitsZSl3","Hits near track versus z, Ax3",150,-150,150,13,0,13);
  
  foundHitsIsoZSl0 = new TH2F("foundHitsIsoZSl0","HitsIso near track versus z, Ax0",150,-150,150,13,0,13);
  foundHitsIsoZSl1 = new TH2F("foundHitsIsoZSl1","HitsIso near track versus z, Ax1",150,-150,150,13,0,13);
  foundHitsIsoZSl2 = new TH2F("foundHitsIsoZSl2","HitsIso near track versus z, Ax2",150,-150,150,13,0,13);
  foundHitsIsoZSl3 = new TH2F("foundHitsIsoZSl3","HitsIso near track versus z, Ax3",150,-150,150,13,0,13);
  
  foundHitsIsoUMZSl0 = new TH2F("foundHitsIsoUMZSl0","HitsIsoUM near track versus z, Ax0",150,-150,150,13,0,13);
  foundHitsIsoUMZSl1 = new TH2F("foundHitsIsoUMZSl1","HitsIsoUM near track versus z, Ax1",150,-150,150,13,0,13);
  foundHitsIsoUMZSl2 = new TH2F("foundHitsIsoUMZSl2","HitsIsoUM near track versus z, Ax2",150,-150,150,13,0,13);
  foundHitsIsoUMZSl3 = new TH2F("foundHitsIsoUMZSl3","HitsIsoUM near track versus z, Ax3",150,-150,150,13,0,13);
  
}


//_____________________________________________________________________________
int TXftAnaModule::BeginJob() {

  // register the data block(s)
  //
  // First do the Trigger
  fTriggerBlock = (TStnTriggerBlock*) RegisterDataBlock("TriggerBlock",
							"TStnTriggerBlock");
  if (! fTriggerBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","triggerBlock");
    fEnabled = 0;
  }
  //
  // Second do the Header
  fHeaderBlock = (TStnHeaderBlock*) RegisterDataBlock("HeaderBlock",
						      "TStnHeaderBlock");
  if (! fHeaderBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","headerBlock");
    fEnabled = 0;
  }
  //
  // Third do the XFT
  fXftBlock = (TXftBlock*) RegisterDataBlock("XftBlock",
					     "TXftBlock");
  if (! fXftBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","XftBlock");
    fEnabled = 0;
  }
//
// The XFT info is in the following format
//
// finder hit info (in Stntuple/Stntuple/data/TXftHit.hh)
//   int      hPrompt; //
//   int      hDelayed; //
//   int      hSl;           //
//   int      hAbsCell;           //  
//   int      hLocalCell;           //
//   int      hCrate;           //
//   int      hBoard;           //
//   int      hChip;           //
//
// To loop over and access the hits:
//   int nHits = fXftBlock->fNxftHits;
//   for (Int_t jx=0; jx<nHits; jx++) {
//    	TXftHit& xfthit = *(fXftBlock->XftHit(jx));
//       ...
//       The hit info is xfthit.hPrompt
//
// finder pixel info (in Stntuple/Stntuple/data/TXftPixel.hh)
//   int      pSl;           //
//   int      pAbsPix;           //
//   int      pCrate;           //
//   int      pBoard;           //
//   int      pChip;           //
//   int      pLocalPix;           //
//   int      pLocalCell;           //
//
// To loop over and access the pixels:
//   int nPixels = fXftBlock->fNxftPixels;
//   for (Int_t jx=0; jx<nPixels; jx++) {
//    	TXftPixel& xftpixel = *(fXftBlock->XftPixel(jx));
//       ...
//       The pixel info is xftpixel.pAbsPix
//
//	linker tracks (in Stntuple/Stntuple/data/TXftTrack.hh)
//   int      ptBin;         // bin number
//   int      miniPhi;       // phi within Linker chip
//   int      shortBit;       // short track bit
//   int      isolationBit;  // isolation
//   int      linkerChip;     // Linker on board
//   int      linkerBoard;   // Linker board for this track
//   int      linkerCrate;   // Linker crate for this track
//   float    absPhi;        // absolute phi in radians
//   float    ptVal;            // pt value (after bin)
//
// To loop over and access the tracks:
//   int nTracks = fXftBlock->fNxftTracks;
//   for (Int_t jx=0; jx<nTracks; jx++) {
//    	TXftTrack& xfttrack = *(fXftBlock->XftTrack(jx));
//       ...
//       The pt info is xfttrack.ptBin


  fTrackBlock = (TStnTrackBlock*) RegisterDataBlock("TrackBlock",
						    "TStnTrackBlock");
  if (! fTrackBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","TrackBlock");
    fEnabled = 0;
  }
  
  // book histograms
  BookHistograms();
  //
  // Read in the road database for 4 layer roads
  char filename[200];
  sprintf(filename,"%s/%s",
	  getenv("PROJECT_DIR"),
	  "/XFTSim/XFTDatabase/Linker4LayerRoads_090099.dat");
  roadFile(filename);
  //
  // Read in the road database for 3 layer roads
  char filename3[200];
  sprintf(filename3,"%s/%s",
	  getenv("PROJECT_DIR"),
	  "/XFTSim/XFTDatabase/Linker3LayerRoads_050100.dat");
  road3File(filename3);
  
  // Read in the cot database
  char filenameCOT[200];
  sprintf(filenameCOT,"%s",
	  "Stntuple/ana/xft_cot_geom.dat");
  cotGeom(filenameCOT);
  
  return 0;
}


//_____________________________________________________________________________
int TXftAnaModule::BeginRun() {
  return 0;
}

//_____________________________________________________________________________
int TXftAnaModule::Event(int ientry) {
  
  //  Efficiency: numerator is number of offline tracks with pt above this matched to XFT
  //  Efficiency: denominator is number of offline tracks with pt above this
  //  Fakes: number of XFT tracks not matched to an offline track with pt above this
  float OFFPT_THRESHOLD = 1.5;
  //
  // Minimum pt for entries into pt spectra plots for turn-on curves	
  float OFFPT_MINIMUM = 1.2;
  //
  // each of these is a TBranch, move the branch to this event
  fXftBlock->GetEntry(ientry);
  fTrackBlock->GetEntry(ientry);
  fTriggerBlock->GetEntry(ientry);
  //
  // Header info
  int runNumber = fHeaderBlock->RunNumber();
  int eventNumber = fHeaderBlock->EventNumber();
  //
  // get number of xft objects
  int nHits = fXftBlock->fNxftHits;
  int nPixels = fXftBlock->fNxftPixels;
  int nTracks = fXftBlock->fNxftTracks;
  int oTracks = fTrackBlock->fNTracks;
  xftNHits->Fill(nHits);
  xftNPixels->Fill(nPixels);

  
// Get the trigger bits for 127032 [28,176
  
  int pTriggerBit_cem_nomatch4 = fTriggerBlock->Tfrd()->PrescaledL1Bit(53);
  int pTriggerBit_cmu_nomatch8 = fTriggerBlock->Tfrd()->PrescaledL1Bit(39);
  int TriggerBit_em8 = fTriggerBlock->Tfrd()->PrescaledL1Bit(3);
  int TriggerBit_jet10 = fTriggerBlock->Tfrd()->PrescaledL1Bit(7);
  int TriggerBit_cem_match8 = fTriggerBlock->Tfrd()->PrescaledL1Bit(17);
  
  int TriggerBit_cem_nomatch4 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(53);
  int TriggerBit_cmu_nomatch8 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(39);
  int TriggerBit_cem_match4 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(15);
  int TriggerBit_cmu_match8 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(47);
  
  // Get the trigger bits for 123499 [15,169]
  if(runNumber == 123499){
    pTriggerBit_cem_nomatch4 = fTriggerBlock->Tfrd()->PrescaledL1Bit(53);
    pTriggerBit_cmu_nomatch8 = fTriggerBlock->Tfrd()->PrescaledL1Bit(35);
    TriggerBit_em8 = fTriggerBlock->Tfrd()->PrescaledL1Bit(1);
    TriggerBit_jet10 = fTriggerBlock->Tfrd()->PrescaledL1Bit(3);
    
    TriggerBit_cem_nomatch4 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(53);
    TriggerBit_cmu_nomatch8 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(35);
    TriggerBit_cem_match4 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(13);
    TriggerBit_cmu_match8 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(41);
  }
  //Get the trigger bits for 123917 [22,170]
  if(runNumber == 123917) {
    pTriggerBit_cem_nomatch4 = fTriggerBlock->Tfrd()->PrescaledL1Bit(53);
    pTriggerBit_cmu_nomatch8 = fTriggerBlock->Tfrd()->PrescaledL1Bit(39);
    TriggerBit_em8 = fTriggerBlock->Tfrd()->PrescaledL1Bit(3);
    TriggerBit_jet10 = fTriggerBlock->Tfrd()->PrescaledL1Bit(7);
    
    TriggerBit_cem_nomatch4 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(53);
    TriggerBit_cmu_nomatch8 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(39);
    TriggerBit_cem_match4 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(15);
    TriggerBit_cmu_match8 = fTriggerBlock->Tfrd()->UnprescaledL1Bit(47);
  }
  
  //	if(pTriggerBit_cem_nomatch4 == 0 && pTriggerBit_cmu_nomatch8 == 0) return;
  
  //	int TriggerBit_jet5 = fTriggerBlock->Tfrd()->PrescaledL1Bit(5);
  //	int TriggerBit_met = fTriggerBlock->Tfrd()->PrescaledL1Bit(7);
  //	int TriggerBit_sumet = fTriggerBlock->Tfrd()->PrescaledL1Bit(9);
  //	int TriggerBit_mb = fTriggerBlock->Tfrd()->PrescaledL1Bit(21);
  //	int TriggerBit_zb = fTriggerBlock->Tfrd()->PrescaledL1Bit(23);
  
  int unprescaledWord0 = fTriggerBlock->Tfrd()->Word(3);
  int unprescaledWord1 = fTriggerBlock->Tfrd()->Word(4);
  int prescaledWord0 = fTriggerBlock->Tfrd()->Word(5);
  int prescaledWord1 = fTriggerBlock->Tfrd()->Word(6);
  int myron = fTriggerBlock->Tsid()->MyronFlag();
  
  //	if(TriggerBit_cem_match8 == 0) return;	
  // Some printout
  //   printf("Number of XFT Hits %d\n",nHits);
  //   printf("Number of XFT Pixels %d\n",nPixels);
  //   printf("Number of XFT Tracks %d\n",nTracks);
  //  printf("Number of Offline Tracks %d\n",oTracks);
  
  //	if(nTracks>0) printf("\n Event %d NXFT tracks %d NCOT tracks%d",eventNumber,nTracks,oTracks);
  //-------------------------------------------------------------------
  // LOOP OVER XFT HITS
  //-------------------------------------------------------------------
  for (Int_t hx=0; hx<nHits; hx++) {
    TXftHit& xfthit = *(fXftBlock->XftHit(hx));
    if (xfthit.hAbsCell < minCell[xfthit.hSl-1]) minCell[xfthit.hSl-1] = xfthit.hAbsCell;
    if (xfthit.hAbsCell > maxCell[xfthit.hSl-1]) maxCell[xfthit.hSl-1] = xfthit.hAbsCell;
    for (Int_t iwire=0; iwire<12; iwire++) {
      int iprompt = (xfthit.hPrompt >> iwire) & 1;
      int idelayed = (xfthit.hDelayed >> iwire) & 1;
      if (iprompt==1 && idelayed==0) {
	if (xfthit.hSl==1) cellsPrompt0->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==2) cellsPrompt1->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==3) cellsPrompt2->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==4) cellsPrompt3->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==1) wiresPrompt0->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==2) wiresPrompt1->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==3) wiresPrompt2->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==4) wiresPrompt3->Fill(xfthit.hAbsCell*12 + iwire);
      }
      if (iprompt==0 && idelayed==1) {
	if (xfthit.hSl==1) cellsDelayed0->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==2) cellsDelayed1->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==3) cellsDelayed2->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==4) cellsDelayed3->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==1) wiresDelayed0->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==2) wiresDelayed1->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==3) wiresDelayed2->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==4) wiresDelayed3->Fill(xfthit.hAbsCell*12 + iwire);
      }
      if (iprompt==1 && idelayed==1) {
	if (xfthit.hSl==1) cellsPD0->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==2) cellsPD1->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==3) cellsPD2->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==4) cellsPD3->Fill(xfthit.hAbsCell);
	if (xfthit.hSl==1) wiresPD0->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==2) wiresPD1->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==3) wiresPD2->Fill(xfthit.hAbsCell*12 + iwire);
	if (xfthit.hSl==4) wiresPD3->Fill(xfthit.hAbsCell*12 + iwire);
      }
    }
  }
  
  

  //-------------------------------------------------------------------
  // LOOP OVER XFT PIXELS
  // Will get two entries for high pt track segments in SL3 and SL4
  //-------------------------------------------------------------------
  for (Int_t jx=0; jx<nPixels; jx++) {
    TXftPixel& xftpixel = *(fXftBlock->XftPixel(jx));
    //
    // SL1 Offset: 6 pixels + (4 cells * 12 pixels); number of pixels 192 cells *12 pixels/cell
    // SL2 Offset: 6 pixels + (4 cells * 12 pixels); number of pixels 288 cells *12 pixels/cell
    // SL3 Offset: 3 pixels + (8 cells *  6 pixels); number of pixels 384 cells * 6 pixels/cell
    // SL4 Offset: 3 pixels + (4 cells *  6 pixels); number of pixels 480 cells * 6 pixels/cell
    //
    // Bin histograms per chip and per cell
    //
    float cellPhi = 0.;
    float chipPhi = 0.;
    if(xftpixel.pSl == 1){
      int zeroPix = abs(xftpixel.pAbsPix) +6;
      cellPhi =(zeroPix / 12.) *(twopi/192.);
      chipPhi =(zeroPix / 48.) *(twopi/48.);
      SL1PixelCell->Fill( cellPhi );
      SL1PixelChip->Fill( chipPhi );
      
    }else if(xftpixel.pSl == 2){	
      int zeroPix = abs(xftpixel.pAbsPix) +6;
      cellPhi =(zeroPix / 12.) *(twopi/288.);
      chipPhi =(zeroPix / 48.) *(twopi/72.);
      SL2PixelCell->Fill( cellPhi );
      SL2PixelChip->Fill( chipPhi );
      
    }else if(xftpixel.pSl == 3){	
      int zeroPix = abs(xftpixel.pAbsPix) +3;
      cellPhi =(zeroPix / 6.) *(twopi/384.);
      chipPhi =(zeroPix / 24.) *(twopi/96.);
      SL3PixelCell->Fill( cellPhi );
      SL3PixelChip->Fill( chipPhi );
      
    }else if(xftpixel.pSl == 4){	
      int zeroPix = abs(xftpixel.pAbsPix)+3;
      cellPhi =(zeroPix / 6.) *(twopi/480.);
      chipPhi =(zeroPix / 24.) *(twopi/120.);
      SL4PixelCell->Fill( cellPhi );
      SL4PixelChip->Fill( chipPhi );
    }								
    //			printf("\n Phi %f, SL %d, Crate %d, Board %d, Chip %d, LocalCell %d, LocalPix %d, ABsPix %d",
    //				cellPhi,xftpixel.pSl,xftpixel.pCrate,xftpixel.pBoard, xftpixel.pChip,
    //				xftpixel.pLocalCell,xftpixel.pLocalPix,xftpixel.pAbsPix);
    
  }
  //---------------------------------------------------------------------------
  // LOOP OVER XFT TRACKS
  //---------------------------------------------------------------------------
  //	Fill pt, phi, crate, board, chip histograms
  //  Try to match xft track to an offline track
  //	Unmatched tracks are fakes - fill pt,phi,crate,board,chip histograms
  //---------------------------------------------------------------------------
  
  int xftN4layer = 0;
  for (Int_t jx=0; jx<nTracks; jx++) {
    TXftTrack& xfttrack = *(fXftBlock->XftTrack(jx));
    if(xfttrack.shortBit     == 1) xftPtBinShort->Fill(xfttrack.ptBin);
    //
    // ANALYSE 4-LAYER TRACKS ONLY
    if(xfttrack.isolationBit == 1){	
      
      xftN4layer++;
      float ptval = getPtValue(xfttrack.ptBin,xfttrack.isolationBit,
			       xfttrack.shortBit,xfttrack.absPhi);
      
      float phiBin = (288./2304.)*(xfttrack.miniPhi
				   +8*(xfttrack.linkerChip
				       +12*(xfttrack.linkerBoard+8*xfttrack.linkerCrate)));
      
      float phiZero = getPhiZero(ptval,xfttrack.absPhi);
      
      if (prtlvl > 0) {
	printf("\nEvent %d, myron %d, XFT pT bin %d, pT %f, phi %f, phi0 %f, short %d ,iso %d, crate %d, board %d, chip %d, minphi %d",
	       eventNumber, myron, xfttrack.ptBin, ptval, xfttrack.absPhi, phiZero,
	       xfttrack.shortBit, xfttrack.isolationBit,
	       xfttrack.linkerCrate, xfttrack.linkerBoard, xfttrack.linkerChip, xfttrack.miniPhi);
	printf("\nPrescaled trigger words %x %x, Unprescaled trigger words %x %x",	
	       prescaledWord0, prescaledWord1, unprescaledWord0, unprescaledWord1);
      }
      //
      // Fill pt, phi, crate, board, chip histograms for all XFT tracks
      xftPtBin->Fill(xfttrack.ptBin);
      xftPtVal->Fill(ptval);
      
      xftShort->Fill(xfttrack.shortBit);
      xftIso->Fill(xfttrack.isolationBit);
      
      xftAbsPhi->Fill(xfttrack.absPhi);
      xftCrate->Fill(xfttrack.linkerCrate);
      xftBoard->Fill(xfttrack.linkerBoard+ 8*xfttrack.linkerCrate);
      xftChip->Fill(xfttrack.linkerChip
		    +12*(xfttrack.linkerBoard+8*xfttrack.linkerCrate));
      xftPhiBin->Fill(phiBin);			
      //
      //Fill ptBin and Phi plots for various trigger bits			
      if(TriggerBit_jet10 == 1){
	xftpt_jet->Fill(xfttrack.ptBin);
	xftphi_jet->Fill(xfttrack.absPhi);
      }
      
      if(xfttrack.ptBin >35 && xfttrack.ptBin < 60){ 	
	xftAbsPhi8->Fill(xfttrack.absPhi);				
	if(TriggerBit_cmu_nomatch8 == 1){				 
	  xftpt_cmu_nomatch8->Fill(xfttrack.ptBin);
	  xftphi_cmu_nomatch8->Fill(xfttrack.absPhi);
	}
	if(TriggerBit_cmu_match8 == 1){
	  xftpt_cmu_match->Fill(xfttrack.ptBin);
	  xftphi_cmu_match->Fill(xfttrack.absPhi);
	}
      }
      
      if(xfttrack.ptBin >23 && xfttrack.ptBin < 72){ 	
	xftAbsPhi4->Fill(xfttrack.absPhi);				
	if(TriggerBit_cem_nomatch4 == 1){				 
	  xftpt_cem_nomatch4->Fill(xfttrack.ptBin);
	  xftphi_cem_nomatch4->Fill(xfttrack.absPhi);
	}
	if(TriggerBit_cem_match4 == 1){
	  xftpt_cem_match->Fill(xfttrack.ptBin);
	  xftphi_cem_match->Fill(xfttrack.absPhi);
	}
      }
      //
      //TRY TO MATCH THIS XFT TRACK TO AN OFFLINE TRACK
      //Require ofline track pt to be above some minimum value
      //Match if within X pixels in 3 out of 4 superlayers
      bool matched = false;
      if (ptval != 0.0) { //indicates pt bin is out of array bounds!
	xft_Total++;
	//
	// Loop over L3 tracks and try to match to the XFT track
	for (Int_t jg=0; jg<oTracks; jg++) {
	  TStnTrack& track = *(fTrackBlock->Track(jg));
	  if (fabs(track.Pt()) > OFFPT_MINIMUM) {
	    int prt = 0;
	    int nmatch = matchTracks(prt,ptval,phiZero, 0.0, 0.0, 0.0,
				     track.Charge()*track.Pt(),track.Phi0(),track.Lam0(),
				     track.D0(),track.Z0());
	    if (nmatch >= 3) matched = true;
	    
	    if (prtlvl > 0) printf("   Gen Pt %f, phi %f, pz %f \n",
				   track.Pt()*track.Charge(),track.Phi0(),track.Pt()*track.Lam0());
	    if (prtlvl > 0) printf("   match = %d \n",nmatch);
	  }
	}
	//
	//No match => FAKE XFT track 
	//Fill pt, phi, crate, board, chip histograms for fake xft tracks				
	if (!matched) {
	  xft_unmatched++;
	  xftFakePtBin->Fill(xfttrack.ptBin);
	  xftFakePtVal->Fill(ptval);
	  
	  xftFakeShort->Fill(xfttrack.shortBit);
	  xftFakeIso->Fill(xfttrack.isolationBit);
	  
	  xftFakeAbsPhi->Fill(xfttrack.absPhi);
	  xftFakeCrate->Fill(xfttrack.linkerCrate);
	  xftFakeBoard->Fill(xfttrack.linkerBoard+ 8*xfttrack.linkerCrate);
	  xftFakeChip->Fill(xfttrack.linkerChip
			    +12*(xfttrack.linkerBoard+8*xfttrack.linkerCrate));
	  xftFakePhiBin->Fill(phiBin);			
	}	
      } //ptval check
    }//4-layer track selection
  }//loop over xft tracks
  
  //Only fill if we have at least 1 XFT 4-layer track
  //Interested to see if high fake rate coming from L3 skipping track reco if photon trigger
  
  if(xftNTracks > 0 ){
    xftNTracks->Fill(xftN4layer);
    offNTracks->Fill(oTracks);
  }	
  //---------------------------------------------------------------------------
  // LOOP OVER OFFLINE TRACKS TWICE
  //---------------------------------------------------------------------------
  // First Flag fiducial, offline track with the highest pt in a linker chip
  // Second Fill histograms for these flagged tracks at least X stereo hits
  //
  // Initialise arrays that will store information on highest pt track in each linker chip
  int trackIndex[288];
  float maxPt[288];
  for (Int_t ilink = 0; ilink<288; ilink++) {
    trackIndex[ilink] = -1;
    maxPt[ilink]      = 0.;
  }
  
  // FIRST LOOP OVER ALL OFFLINE TRACKS
  bool goodTrack[2000];
  if(oTracks > 2000) oTracks = 2000; //avoid array overflows
  for(Int_t jg=0; jg<oTracks; jg++){
    TStnTrack& track = *(fTrackBlock->Track(jg));
    goodTrack[jg] = false;
    //
    //speed things up by skipping very many, very low pt tracks
    if(fabs(track.Pt()) > OFFPT_MINIMUM){ 
      if(fabs(track.Pt()) > OFFPT_THRESHOLD){
	//
	// Find the z-postion of the track at each SL
	float zpos[4];
	for (int isl=1; isl<8; isl=isl+2) {
	  float xyz[3];
	  track_Prop(rsns[isl], track.C0(),cdfphi(track.Phi0()),
		     track.Lam0(),track.D0(),track.Z0(),xyz);
	  zpos[(isl-1)/2] = xyz[2];
	}
	allPhi->Fill(track.Phi0());
	alld0->Fill(track.D0());
	allz0->Fill(track.Z0());
	allCotan->Fill(track.Lam0());
	allEta->Fill(track.Eta());
	allPt->Fill(track.Pt());
	allZSl3->Fill(zpos[3]);
	allZSl2->Fill(zpos[2]);
	allZSl1->Fill(zpos[1]);
	allZSl0->Fill(zpos[0]);
	if(track.NCotHitsSt() > 20 && track.NCotHitsAx() > 20){
	  qallPhi->Fill(track.Phi0());
	  qalld0->Fill(track.D0());
	  qallz0->Fill(track.Z0());
	  qallCotan->Fill(track.Lam0());
	  qallEta->Fill(track.Eta());
	  qallPt->Fill(track.Pt());
	  qallZSl3->Fill(zpos[3]);
	  qallZSl2->Fill(zpos[2]);
	  qallZSl1->Fill(zpos[1]);
	  qallZSl0->Fill(zpos[0]);
	}
      }
      //
      // Calculate linker chip number for fiducial, offline tracks			
      if( fiducial(track.Charge()*track.Pt(),track.Phi0(),track.Lam0(),track.D0(),track.Z0())){
	float rpix = returnHits(5,track.C0(),cdfphi(track.Phi0()),
				track.Lam0(),track.D0(),track.Z0());
        int jlink = (int) (rpix / 8.0);
	//									
	// Highest pt, fiducial, offline track in this linker chip so far => update track index[] and maxPt[]
	if( fabs(track.Pt()) > maxPt[jlink] ){
	  trackIndex[jlink] = jg;
	  maxPt[jlink]      = fabs(track.Pt());
	}
      }
    }
  }
  //
  // Now set flag for highest pt, fiducial, offline track in each linker chip	 
  for (Int_t ilink = 0; ilink<288; ilink++) {
    if( trackIndex[ilink] != -1 ) goodTrack[ trackIndex[ilink] ] = true;		
  }	
  //
  // SECOND LOOP OVER ALL OFFLINE TRACKS
  float zpos[4];
  for (Int_t jg=0; jg<oTracks; jg++) {
    TStnTrack& track = *(fTrackBlock->Track(jg));      
    if(goodTrack[jg] && track.NCotHitsSt() > 20 && track.NCotHitsAx() > 20){
      //
      // Is this track isolated?
      int tNear=0;
      for (Int_t kg=0; kg<oTracks; kg++) {
	if (kg != jg) {
	  TStnTrack& ktrack = *(fTrackBlock->Track(kg));
	  int closeLayers = nearTracks(track.C0(),track.Phi0(),
				       track.Lam0(),track.D0(),track.Z0(),
				       ktrack.C0(),ktrack.Phi0(),
				       ktrack.Lam0(),ktrack.D0(),ktrack.Z0());
	  if (closeLayers >= 1) tNear++;
	}
      }
      //
      // Find the z-postion of the track at each SL
      for (int isl=1; isl<8; isl=isl+2) {
	float xyz[3];
	track_Prop(rsns[isl], track.C0(),cdfphi(track.Phi0()),
		   track.Lam0(),track.D0(),track.Z0(),xyz);
	zpos[(isl-1)/2] = xyz[2];
	//               std::cout << "zpos " << zpos[(isl-1)/2] 
	//               		<< " isl " << (isl-1)/2 << std::endl;
      }
      //
      //Fill pt spectra histograms for turn-on curves
      float inversePt = 0.;
      if(track.Pt() > 0.1) inversePt = 1./track.Pt();
      
      offPtTriggerAll1p5->Fill(track.Pt());
      offPtTriggerAll4->Fill(track.Pt());
      offPtTriggerAll8->Fill(track.Pt());
      offInvPtTriggerAll->Fill(inversePt);
      if(track.Charge()>0){ //positive charge
	poffPtTriggerAll1p5->Fill(track.Pt());
	poffPtTriggerAll4->Fill(track.Pt());
	poffPtTriggerAll8->Fill(track.Pt());
      }else{		//negative charge
	moffPtTriggerAll1p5->Fill(track.Pt());
	moffPtTriggerAll4->Fill(track.Pt());
	moffPtTriggerAll8->Fill(track.Pt());
      }
      //
      //Fill histograms for efficiencies (offline pt above 1.5 GeV/c)
      if(fabs(track.Pt()) > OFFPT_THRESHOLD){
	nfid++;
	offPhi->Fill(track.Phi0());
	offd0->Fill(track.D0());
	offz0->Fill(track.Z0());
	offCotan->Fill(track.Lam0());
	offEta->Fill(track.Eta());
	offPt->Fill(track.Pt());
	offphid0->Fill(cdfphi(track.Phi0()),track.D0());
	CotHitsTot->Fill(track.NCotHitsTot());
	CotHitsAx->Fill(track.NCotHitsAx());
	CotHitsSt->Fill(track.NCotHitsSt());
	offZSl3->Fill(zpos[3]);
	offZSl2->Fill(zpos[2]);
	offZSl1->Fill(zpos[1]);
	offZSl0->Fill(zpos[0]);
	CotHitsTotZSl3->Fill(zpos[3],track.NCotHitsTot());
	CotHitsAxZSl3->Fill(zpos[3],track.NCotHitsAx());
	CotHitsStZSl3->Fill(zpos[3],track.NCotHitsSt());
      }
      //
      // Now see if this matches an XFT track
      bool matched = false;
      for (Int_t jx=0; jx<nTracks; jx++) {
	TXftTrack& xfttrack = *(fXftBlock->XftTrack(jx));
	
	float ptval = getPtValue(xfttrack.ptBin,xfttrack.isolationBit,
				 xfttrack.shortBit,xfttrack.absPhi);
	if (ptval != 0.0 && (xfttrack.isolationBit == 1)) {
	  float phiZero = getPhiZero(ptval,xfttrack.absPhi);
	  int prt = 0;
	  int nmatch = 0;
	  nmatch = matchTracks(prt,ptval,phiZero, 0.0, 0.0, 0.0,
			       track.Charge()*track.Pt(),track.Phi0(),track.Lam0(),
			       track.D0(),track.Z0());
	  //
	  // If match is good enough, remember track
	  if (nmatch >= 3) {
	    matched = true;
	    index_match = jx;
	  }
	}//ptval and 4-layer checks
      }//loop over xft tracks
      //
      //HISTOGRAMS FOR FOUND AND MISSED OFFLINE TRACKS
      if (matched) {
	
	TXftTrack& xfttrack = *(fXftBlock->XftTrack(index_match));
	float ptval = getPtValue(xfttrack.ptBin,xfttrack.isolationBit,
				 xfttrack.shortBit,xfttrack.absPhi);
	//
	//Fill pt spectra histograms for found tracks for turn-on curves
	offPtTrigger1p5->Fill(track.Pt());
	offInvPtTrigger1p5->Fill(inversePt);
	if(track.Charge()>0)	poffPtTrigger1p5->Fill(track.Pt());
	if(track.Charge()<0)	moffPtTrigger1p5->Fill(track.Pt());
	
	if(xfttrack.ptBin >23 && xfttrack.ptBin < 72){ 	
	  offPtTrigger4->Fill(track.Pt());
	  offInvPtTrigger4->Fill(inversePt);
	  if(track.Charge()>0)	poffPtTrigger4->Fill(track.Pt());
	  if(track.Charge()<0)	moffPtTrigger4->Fill(track.Pt());
	}
	if(xfttrack.ptBin >35 && xfttrack.ptBin < 60){
	  offPtTrigger8->Fill(track.Pt());
	  offInvPtTrigger8->Fill(inversePt);
	  if(track.Charge()>0)	poffPtTrigger8->Fill(track.Pt());
	  if(track.Charge()<0)	moffPtTrigger8->Fill(track.Pt());
	}
	
	//
	//Fill histograms for found tracks for efficiencies
	if(fabs(track.Pt()) > OFFPT_THRESHOLD){
	  nfid_matched++;
	  offPhiFound->Fill(track.Phi0());
	  offd0Found->Fill(track.D0());
	  offz0Found->Fill(track.Z0());
	  offCotanFound->Fill(track.Lam0());
	  offPtFound->Fill(track.Pt());
	  offEtaFound->Fill(track.Eta());
	  CotHitsTotFound->Fill(track.NCotHitsTot());
	  CotHitsAxFound->Fill(track.NCotHitsAx());
	  CotHitsStFound->Fill(track.NCotHitsSt());
	  offZSl3Found->Fill(zpos[3]);
	  offZSl2Found->Fill(zpos[2]);
	  offZSl1Found->Fill(zpos[1]);
	  offZSl0Found->Fill(zpos[0]);
	  CotHitsTotZSl3Found->Fill(zpos[3],track.NCotHitsTot());
	  CotHitsAxZSl3Found->Fill(zpos[3],track.NCotHitsAx());
	  CotHitsStZSl3Found->Fill(zpos[3],track.NCotHitsSt());
	  //
	  //FILL RESOLUTION HISTOGRAMS
	  //
	  // Pt, Phi resolution for all, positive and negative tracks 
	  //					float phiZero = getPhiZero(ptval,xfttrack.absPhi);
	  float rtest = -10.0;
	  float phiZero = getPhiAtRadius(ptval,xfttrack.absPhi,rtest);
	  
	  dPhiHist->Fill(phiZero - track.Phi0());
	  dPtHist->Fill((ptval-track.Charge()*track.Pt())/ (track.Pt()*track.Pt()));
	  if (ptval > 0.0) {
	    PdPhiHist->Fill(phiZero - track.Phi0());
	    PdPtHist->Fill((ptval-track.Charge()*track.Pt())/ (track.Pt()*track.Pt()));
	  } else {
	    MdPhiHist->Fill(phiZero - track.Phi0());
	    MdPtHist->Fill((ptval-track.Charge()*track.Pt())/ (track.Pt()*track.Pt()));
	  }
	  //
	  //Apply pt(phi) correction for beam-spot offset if necessary
	  float ptvalCorr = getPtValue(xfttrack.ptBin,xfttrack.isolationBit,
				       xfttrack.shortBit,xfttrack.absPhi);
	  int ptBinCorr = nearestPtBin(ptvalCorr);
	  
	  if (ptvalCorr != 0.0) {
	    dPtNewHist->Fill((ptvalCorr-track.Charge()*track.Pt()) / (track.Pt()*track.Pt()));
	  }
	  
	  //Phi, pt resolution vs phi						
	  dPhiPhiHist->Fill(track.Phi0(),phiZero - track.Phi0());
	  dPtPhiHist->Fill(cdfphi(track.Phi0()),(ptval-track.Charge()*track.Pt()) / (track.Pt()*track.Pt()));
	  
	  //
	  // Find the nearest true pt bin for the offline track
	  int ptBinOff = nearestPtBin(track.Charge()*track.Pt());
	  if (ptBinOff < 128) {
	    PtBinPtBinHist->Fill(ptBinOff,xfttrack.ptBin);
	    
	    // Plot (PtBin Offline - PtBin XFT) vs chip, phi for all, positive and negative tracks
	    float diffBin = ptBinOff - xfttrack.ptBin;
	    float diffBinCorr = ptBinOff - ptBinCorr;
	    
	    dPtBinChipHist->Fill(xfttrack.linkerChip+12*xfttrack.linkerBoard 
				 +12*8*xfttrack.linkerCrate,diffBin);
	    dPtBinPhiHist->Fill(track.Phi0(),diffBin);
            if (track.Charge() > 0.0) {
	      dPtBinPhiPos->Fill(track.Phi0(),diffBin);
	      dPtBinPhiPosCorr->Fill(track.Phi0(),diffBinCorr);
            } else {
	      dPtBinPhiNeg->Fill(track.Phi0(),diffBin);
	      dPtBinPhiNegCorr->Fill(track.Phi0(),diffBinCorr);
            }
	  }
	  if (prtlvl > 0) printf("\nOffline Pt %f, bin %d, xft pt %f bin %d corrected pt %f bin %d\n",
				 track.Charge()*track.Pt(), ptBinOff,
				 ptval, xfttrack.ptBin, ptvalCorr, ptBinCorr);						
	  
	}//Min pt cut for efficiencies
	//
	//FILL MISSED HISTOGRAMS - XFT INEFFICIENCY 											
      } else {
	if(fabs(track.Pt()) > OFFPT_THRESHOLD){
	  offPhiMissed->Fill(track.Phi0());
	  offd0Missed->Fill(track.D0());
	  offz0Missed->Fill(track.Z0());
	  offPtMissed->Fill(track.Pt());
	  offCotanMissed->Fill(track.Lam0());
	  offEtaMissed->Fill(track.Eta());
	  CotHitsTotMissed->Fill(track.NCotHitsTot());
	  CotHitsAxMissed->Fill(track.NCotHitsAx());
	  CotHitsStMissed->Fill(track.NCotHitsSt());
	  offZSl3Missed->Fill(zpos[3]);
	  offZSl2Missed->Fill(zpos[2]);
	  offZSl1Missed->Fill(zpos[1]);
	  offZSl0Missed->Fill(zpos[0]);
	  CotHitsTotZSl3Missed->Fill(zpos[3],track.NCotHitsTot());
	  CotHitsAxZSl3Missed->Fill(zpos[3],track.NCotHitsAx());
	  CotHitsStZSl3Missed->Fill(zpos[3],track.NCotHitsSt());
	  
	  if (prtlvl > 0) {
	    printf("\nUnmatched Run %d Event %d \n",
		   runNumber, eventNumber);
	    printf("\nUnmatched Gen Pt %f, phi %f, cotan %f Z0 %f d0 %f \n",
		   track.Charge()*track.Pt(), track.Phi0(), track.Lam0(),
		   track.Z0(), track.D0());
	    printf("\n  Alg %d COT total hits %d ax hits %d ax segs %d\n",
		   track.Algorithm(), track.NCotHitsTot(),
		   track.NCotHitsAx(), track.NAxSeg());
	  }	
	}//Min pt cut for efficiencies				
      } //matched check
      //
      //
      // TRY TO MATCH PROMPT/DELAYED HITS
      if(fabs(track.Pt()) > OFFPT_THRESHOLD){
	for (int isl=1; isl<8; isl=isl+2) {
	  int nHit = 0;
	  for (Int_t iwire=0; iwire<12; iwire++) {
	    //
	    // For this superlayer and wire, which cell does the track go through?
	    int icell = returnCell(isl,iwire,
				   track.C0(),cdfphi(track.Phi0()),
				   track.Lam0(),track.D0(),track.Z0());
	    
	    //
	    // Is there a hit attached to this track?
	    int tdcHit = 0;
	    int absWire = 12*isl + iwire;
	    if (absWire < 32) {
	      int wirePos = absWire;
	      tdcHit = (track.CotHitMask(0) >> wirePos) & 1;
	    } else if (absWire < 64) {
	      int wirePos = absWire - 32;
	      tdcHit = (track.CotHitMask(1) >> wirePos) & 1;
	    } else {
	      int wirePos = absWire - 64;
	      tdcHit = (track.CotHitMask(2) >> wirePos) & 1;
	    }
	    if (tdcHit == 0) {
	      if ((isl+1)/2==1) wiresMissedTDC0->Fill(icell*12 + iwire);
	      if ((isl+1)/2==2) wiresMissedTDC1->Fill(icell*12 + iwire);
	      if ((isl+1)/2==3) wiresMissedTDC2->Fill(icell*12 + iwire);
	      if ((isl+1)/2==4) wiresMissedTDC3->Fill(icell*12 + iwire);
	    }else{
	      if ((isl+1)/2==1) wiresFoundTDC0->Fill(icell*12 + iwire);
	      if ((isl+1)/2==2) wiresFoundTDC1->Fill(icell*12 + iwire);
	      if ((isl+1)/2==3) wiresFoundTDC2->Fill(icell*12 + iwire);
	      if ((isl+1)/2==4) wiresFoundTDC3->Fill(icell*12 + iwire);
	    }
	    
	    //						printf("isl %d iwire %d icell %d pt %f \n",
	    //                  		isl,iwire,icell,konst/track.C0());
	    //
	    // Now loop over the hits in the XFT hit block, and look for one within 3 cells
	    bool foundHit = false;
	    for (Int_t hx=0; hx<nHits; hx++) {
	      TXftHit& xfthit = *(fXftBlock->XftHit(hx));
	      //
	      // Make sure the SL's match 
	      if (xfthit.hSl == (isl+1)/2) {
		//
		// Look for a hit on this wire              
		int iprompt = (xfthit.hPrompt >> iwire) & 1;
		int idelayed = (xfthit.hDelayed >> iwire) & 1;
		int dCell = icell - xfthit.hAbsCell;
		if (dCell > 0.5*ncl[isl]) dCell = dCell - ncl[isl];
		if (dCell < -0.5*ncl[isl]) dCell = dCell + ncl[isl];
		//                        printf("iprompt %d idelayed %d cell %d \n",
		//                        		iprompt,idelayed,xfthit.hAbsCell);
		if ((abs(dCell) < 10) && ((iprompt == 1) || (idelayed == 1)) ) {
		  if ((dCell == -1) || (dCell == 0)) {
		    foundHit = true;
		  }
		  if (xfthit.hSl == 1) {
		    HitDcell0->Fill(icell - xfthit.hAbsCell);
		  } else if (xfthit.hSl == 2) {
		    HitDcell1->Fill(icell - xfthit.hAbsCell);
		  } else if (xfthit.hSl == 3) {
		    HitDcell2->Fill(icell - xfthit.hAbsCell);
		  } else {
		    HitDcell3->Fill(icell - xfthit.hAbsCell);
		  }
		}  // end check on dCell and hit present
	      }  // end check on SL
	    } // end loop over hits
	    if (foundHit) { 
	      nHit++;
	      if ((isl+1)/2==1) cellsFound0->Fill(icell);
	      if ((isl+1)/2==2) cellsFound1->Fill(icell);
	      if ((isl+1)/2==3) cellsFound2->Fill(icell);
	      if ((isl+1)/2==4) cellsFound3->Fill(icell);
	    } else {
	      if ((isl+1)/2==1) cellsMissed0->Fill(icell);
	      if ((isl+1)/2==2) cellsMissed1->Fill(icell);
	      if ((isl+1)/2==3) cellsMissed2->Fill(icell);
	      if ((isl+1)/2==4) cellsMissed3->Fill(icell);
	      if ((isl+1)/2==1) wiresMissed0->Fill(icell*12 + iwire);
	      if ((isl+1)/2==2) wiresMissed1->Fill(icell*12 + iwire);
	      if ((isl+1)/2==3) wiresMissed2->Fill(icell*12 + iwire);
	      if ((isl+1)/2==4) wiresMissed3->Fill(icell*12 + iwire);
	    }               
	  } // end loop over wires
	  //
	  // Now accumulate stats
	  //Assign to cell in middle of SL
	  int mcell = returnCell(isl,6,
				 track.C0(),cdfphi(track.Phi0()),
				 track.Lam0(),track.D0(),track.Z0());
	  
	  if (isl == 1) {						
	    foundHits0->Fill(mcell,nHit);
	    foundHitsZSl0->Fill(zpos[0],nHit);
	    if (tNear == 0) {
	      foundHitsIso0->Fill(mcell,nHit);
	      foundHitsIsoZSl0->Fill(zpos[0],nHit);
	    }
	    if ((tNear == 0)&&!matched) {
	      foundHitsIsoUM0->Fill(mcell,nHit);
	      foundHitsIsoUMZSl0->Fill(zpos[0],nHit);
	    }
	  } else if (isl == 3) {
	    foundHits1->Fill(mcell,nHit);
	    foundHitsZSl1->Fill(zpos[1],nHit);
	    if (tNear == 0) {
	      foundHitsIso1->Fill(mcell,nHit);
	      foundHitsIsoZSl1->Fill(zpos[1],nHit);
	    }
	    if ((tNear == 0)&&!matched){
	      foundHitsIsoUM1->Fill(mcell,nHit);
	      foundHitsIsoUMZSl1->Fill(zpos[1],nHit);
	    }
	  } else if (isl == 5) {
	    foundHits2->Fill(mcell,nHit);
	    foundHitsZSl2->Fill(zpos[2],nHit);
	    if (tNear == 0){
	      foundHitsIso2->Fill(mcell,nHit);
	      foundHitsIsoZSl2->Fill(zpos[2],nHit);
	    }
	    if ((tNear == 0)&&!matched){
	      foundHitsIsoUM2->Fill(mcell,nHit);
	      foundHitsIsoUMZSl2->Fill(zpos[2],nHit);
	    }
	  } else if (isl == 7) {
	    foundHits3->Fill(mcell,nHit);
	    foundHitsZSl3->Fill(zpos[3],nHit);
	    if (tNear == 0){
	      foundHitsIso3->Fill(mcell,nHit);
	      foundHitsIsoZSl3->Fill(zpos[3],nHit);
	    }
	    if ((tNear == 0)&&!matched){
	      foundHitsIsoUM3->Fill(mcell,nHit);
	      foundHitsIsoUMZSl3->Fill(zpos[3],nHit);
	    }
	  }
	} // end loop over sl's
      } //pt threshold check
    } //fiducial etc check
  } //loop over Offline tracks
  
  //
  //Try to make 2-trigger plots
  bool really =false;
  if(really){
    //
    //Loop over xft tracks, plot phi between track pairs over with pt>2GeV/c
    for (Int_t jx=0; jx<nTracks; jx++) {
      TXftTrack& xfttrack1 = *(fXftBlock->XftTrack(jx));
      if(xfttrack1.isolationBit ==1){
	float ptval1 =getPtValue(xfttrack1.ptBin,xfttrack1.isolationBit,
				 xfttrack1.shortBit,xfttrack1.absPhi);
	if(fabs(ptval1) >2.0){
	  //Track over 2 GeV and 4-layer
	  for (Int_t kx=jx+1; kx<nTracks; kx++) {
	    TXftTrack& xfttrack2 = *(fXftBlock->XftTrack(kx));
	    if(xfttrack2.isolationBit ==1){
	      float ptval2 =getPtValue(xfttrack2.ptBin,xfttrack2.isolationBit,
				       xfttrack2.shortBit,xfttrack2.absPhi);
	      if(fabs(ptval2) > 2.0){
		
		float pairAngle = fabs(xfttrack1.absPhi - xfttrack2.absPhi) * 360./twopi;
		if(pairAngle >180.) pairAngle = 360. - pairAngle;
		
		//Check for duplicates (would be removed by xtrp)
		bool dupl = xtrp_duplicate_removal(xfttrack1.linkerChip,xfttrack1.miniPhi,
						   xfttrack1.linkerBoard,xfttrack1.linkerCrate,
						   xfttrack2.linkerChip,xfttrack2.miniPhi,
						   xfttrack2.linkerBoard,xfttrack2.linkerCrate);
		if(dupl){
		  printf("\n\n Duplicate tracks - not putting this into the pair plot");
		  printf("\nEvent %d XFT Ptbin %d, pt %f, absphi %f miniphi %d short %d iso %d crate %d board %d chip %d",
			 eventNumber,xfttrack1.ptBin,ptval1,
			 xfttrack1.absPhi,xfttrack1.miniPhi,
			 xfttrack1.shortBit,xfttrack1.isolationBit,
			 xfttrack1.linkerCrate,xfttrack1.linkerBoard,xfttrack1.linkerChip);
		  printf("\nEvent %d XFT Ptbin %d, pt %f, absphi %f miniphi %d short %d iso %d crate %d board %d chip %d",
			 eventNumber,xfttrack2.ptBin,
			 ptval2,xfttrack2.absPhi,xfttrack2.miniPhi,
			 xfttrack2.shortBit,xfttrack2.isolationBit,
			 xfttrack2.linkerCrate,xfttrack2.linkerBoard,xfttrack2.linkerChip);
		  //
		  //both signs
		  duplPairPhi->Fill(pairAngle);
		  if(ptval1*ptval2 >0.0){ //same sign histo
		    duplPairPhiSame->Fill(pairAngle);
		  }else{						//opposite sign histo
		    duplPairPhiOppo->Fill(pairAngle);
		  }
		}else{					
		  //both signs
		  PairPhi->Fill(pairAngle);
		  if(ptval1*ptval2 >0.0){ //same sign histo
		    PairPhiSame->Fill(pairAngle);
		  }else{						//opposite sign histo
		    PairPhiOppo->Fill(pairAngle);
		  }
		}//duplicate check
	      }//pt track 2 check
	    }//4layer check
	  }//for loop over second track
	  
	}//pt track 1 check
      }//4 layer check					
    }//for loop over first track	
  }
  return 0;
}


//_____________________________________________________________________________
int TXftAnaModule::EndJob() {
  //
  //
  // Prinout
  printf("\n\nTotal fiducial offline =%d, found (matched to XFT) =%d \n",
	 nfid,nfid_matched);
  printf("Total XFT = %d, fakes (unmatched to offline) =%d \n",xft_Total,xft_unmatched);
  // Cal the efficiency histos
  //
  // Efficiency histograms
  effCalc(offPhiFound,offPhi,offPhiEff);
  effCalc(offd0Found,offd0,offd0Eff);
  effCalc(offz0Found,offz0,offz0Eff);
  effCalc(offPtFound,offPt,offPtEff);
  effCalc(offCotanFound,offCotan,offCotanEff);
  effCalc(offEtaFound,offEta,offEtaEff);
  effCalc(CotHitsTotFound,CotHitsTot,CotHitsTotEff);
  effCalc(CotHitsAxFound,CotHitsAx,CotHitsAxEff);
  effCalc(CotHitsStFound,CotHitsSt,CotHitsStEff);
  effCalc(offZSl0Found,offZSl0,offZSl0Eff);
  effCalc(offZSl1Found,offZSl1,offZSl1Eff);
  effCalc(offZSl2Found,offZSl2,offZSl2Eff);
  effCalc(offZSl3Found,offZSl3,offZSl3Eff);
  
  effCalc(offPtTrigger1p5,offPtTriggerAll1p5,turnon1p5);
  effCalc(offPtTrigger4,offPtTriggerAll4,turnon4);
  effCalc(offPtTrigger8,offPtTriggerAll8,turnon8);
  
  effCalc(offInvPtTrigger1p5,offInvPtTriggerAll,iturnon1p5);
  effCalc(offInvPtTrigger4,offInvPtTriggerAll,iturnon4);
  effCalc(offInvPtTrigger8,offInvPtTriggerAll,iturnon8);
  
  effCalc(moffPtTrigger1p5,moffPtTriggerAll1p5,mturnon1p5);
  effCalc(moffPtTrigger4,moffPtTriggerAll4,mturnon4);
  effCalc(moffPtTrigger8,moffPtTriggerAll8,mturnon8);
  
  effCalc(poffPtTrigger1p5,poffPtTriggerAll1p5,pturnon1p5);
  effCalc(poffPtTrigger4,poffPtTriggerAll4,pturnon4);
  effCalc(poffPtTrigger8,poffPtTriggerAll8,pturnon8);
  //
  // Write the histogram file out
  outfile->Write();
  outfile->Close();
  //
  // A final message from our sponsors
  printf("----- end job: ---- %s\n",GetName());
  return 0;
}

//-----------------------------------------------
// 4 layer road files
//-----------------------------------------------
void TXftAnaModule::roadFile(char* pathName)
{
  int index = 0;
  int ipix_cell[4];
  int mcpix[4];
  int nPtbins = 0;
  
  // Setup the Input file Stream
  // ---------------------------
  std::ifstream f1(pathName);
  
  
  // Read in the header pixels per cell for each layer
  // -------------------------------------------------
  for(index=0; index < 4; index++)
    {
      f1 >> ipix_cell[index];
    }
  // Read in the core pixels per layer
  // ---------------------------------
  index = 0;
	for(index=0; index < 4; index++)
	  {
	    f1 >> mcpix[index];
	    numPix[index][1] = mcpix[index];
	  }
	
	
	// Read in the number of linkers
	// -----------------------------
	f1 >> nlink;
	
	// Read in the number of Pt bins
	// -----------------------------
	f1 >> nPtbins;
	
	// Read in the Pt bins
	// -------------------
	int temp;
	for(index=0; index < nPtbins; index++)
	  {
	    f1 >> temp;
	    f1 >> ptbins[index];
	    printf("4pt bin %d pt value %g \n",index,ptbins[index]);
	  }
	f1.close();
}//End method

// ------------------------------
// 3 layer Road Files
// ------------------------------
void TXftAnaModule::road3File(char* pathName)
{
  int index = 0;
  int ipix_cell[4];
  int mcpix[4];
  int nPtbins = 0;
  
  std::ifstream f1(pathName);
  
  // Read in the header pixels per cell for each layer
  // -------------------------------------------------
  for(index=0; index < 4; index++)
    {
      f1>>ipix_cell[index];
    }
  // Read in the core pixels per layer
  // ---------------------------------
  for(index=0; index < 4; index++)
    {
      f1 >> mcpix[index];
    }
  // Read in the number of linkers
  // -----------------------------
  f1>>nlink3;
  if(nlink3 != nlink) 
    printf("**** WARNING **** 4-layer and 3-layer roads databases have different number of linkers\n");
  // Read in the number of Pt bins
  // -----------------------------
  f1 >> nPtbins;
  // Read in the Pt bins
  // -------------------
  int temp;
  for(index=0; index < nPtbins; index++)
    {
      f1 >> temp;
      f1 >> ptbins3[index];
      printf("3pt bin %d pt value %g \n",index,ptbins3[index]);
    }
  f1.close();
}

//---------------------------------------------------------------------------------------------
// CONVERT XFT PT BIN TO A VALUE IN GeV/c.
//  Apply corrections for shifted beamspot here if necessary.
//---------------------------------------------------------------------------------------------
float TXftAnaModule::getPtValue(int ptBin, int tfnd4, int tfnd3, float phi) {
  int binind;
  float chg;
  
  if(ptBin < 48){ //Negative tracks
    // Old corrections for shifted beamspot in April 2001 data
    //		 float binOff = 0.336 + 0.953*sin(phi + 1.304);
    //		 float binOff = 0.5 + 0.221 + 0.953*sin(phi + 1.311);
    //	   ptBin = ptBin + binOff;
    binind = abs(ptBin - 47);
    chg = -1.0;
  } else {				//Positive tracks
    // Old corrections for shifted beamspot in April 2001 data	
    //		 float binOff = 0.701 + 1.06*sin(phi + 1.304);
    //		 float binOff = 0.5 + 0.727 + 1.02*sin(phi + 1.213);
    //	   ptBin = ptBin + binOff;
    binind = ptBin - 48;
    chg = 1.0;
  }
  float ptValue = 0.0;
  if (tfnd4 == 1) {
    if (binind <48) {
      ptValue = chg*ptbins[binind];
    }
  } else if (tfnd3 == 1) {
    if (binind < 16) {
      ptValue = chg*ptbins3[binind];
    }
  }
  return ptValue;
}

//------------------------------------------------------------------------
//CONVERT XFT PHI AT AXIAL SL 3 TO PHI 0 (DIRECTION AT INTERACTION REGION)
//------------------------------------------------------------------------
float TXftAnaModule::getPhiZero(float ptXFT, float phiXFT)
{
  float cval = konst / ptXFT;
  float phiprop = phiXFT - asin(105.575*cval);
  if (ptXFT > 0.0) {
    phiprop = phiprop + 0.0025;
  } else {
    phiprop = phiprop + 0.0008;
  }
  return phiprop;
}


//-------------------------------------------------------------------------------
// MATCH 2 TRACKS - USE POSITION AT EACH SUPERLAYER
//-------------------------------------------------------------------------------
int TXftAnaModule::matchTracks(int prt, 
			       float pt1, float phi1, float cotan1, float d1, float z1, 
			       float pt2, float phi2, float cotan2, float d2, float z2)
{
  int matchPix = 0;
  float c1 = konst / pt1;
  float c2 = konst / pt2;
  
  // This loops over all four of the axial superlayers
  for (int i=1; i<8; i=i+2)
    {
      float rpix1 = returnHits(i,c1,cdfphi(phi1),cotan1,d1,z1);
      float rpix2 = returnHits(i,c2,cdfphi(phi2),cotan2,d2,z2);
      float pixDiff = rpix1 - rpix2;
      if (pixDiff < (-0.5*npix_cell[i]*ncl[i]))
	{
	  pixDiff = pixDiff + (float)(npix_cell[i]*ncl[i]);
	}
      if (pixDiff > (0.5*npix_cell[i]*ncl[i]))
	{
	  pixDiff = (float)(npix_cell[i]*ncl[i]) - pixDiff;
	}
      if (pixDiff < 0) pixDiff = - pixDiff;
      
      if (pixDiff < 10.0)
	{
	  matchPix++;
	}
      if (prt > 0) printf("      pixD %f rpix1 %f rpix2 %f match %d \n",
			  pixDiff,rpix1,rpix2,matchPix);
    }
  return matchPix;
}

//-----------------------------------------------------------------------------
// FIND THE CELL THIS TRACK GOES THROUGH IN A GIVEN SUPERLAYER
//-----------------------------------------------------------------------------
float TXftAnaModule::returnHits (int layer, float cot_curv, 
				 float cot_phi0, float cot_cotan, 
				 float cot_d0, float cot_z0)
{
  
  //
  // Local vars
  float phi_int;
  float pix_hit;
  float xyz[3];
  //
  // Find the cell this track goes through
  track_Prop(rsns[layer], 
	     cot_curv, cot_phi0, cot_cotan, 
	     cot_d0, cot_z0, xyz);
  if (fabs(xyz[2]) < 150)
    {
      phi_int = cdfphi(atan2(xyz[1],xyz[0]));
      pix_hit = phi_int / pixPhi[layer];
    }
  else
    {
      pix_hit = -999;
    }
  return pix_hit;
}

//---------------------------------------------------------------------------------------
// APPLY FIDUCIAL CUTS - THESE TRACKS SHOULD BE FOUND BY XFT
// Fiducial if track goes through all 4 AXIAL layers
// Removed pt cut
//-----------------------------------------------------------------------------------------
bool TXftAnaModule::fiducial (float pt, 
			      float cot_phi0, float cot_cotan, 
			      float cot_d0, float cot_z0)
{
  bool fid = false;
  
  float cot_curv = konst / pt;
  //
  // Check SL1
  float xyz_sl1[3];
  float xyz_sl7[3];
  track_Prop(rsns[1], 
	     cot_curv, cot_phi0, cot_cotan, 
	     cot_d0, cot_z0, xyz_sl1);
  track_Prop(rsns[7], 
	     cot_curv, cot_phi0, cot_cotan, 
	     cot_d0, cot_z0, xyz_sl7);
  
  if ((fabs(xyz_sl1[2])<150.0)
      && (fabs(xyz_sl7[2])<150.0) ) fid = true;
  return fid;
}

/**
* Calculates the x,y,z position of a charged track passing through the
* COT, at a given radius from the origin.
* @param radius is the radius from x=0,y=0 to determine the x,y,z
* @param curvature is the half curvature for the charged track
* @param PhiZeroIn is the phi at the origin for the charged track
* @param cotan is the  cot(theta) for the charged track
* @param D_O is the impact parameter for the charged track
* @param Z_O is the z position at the distance of closest approach
*            to x=0,y=0
* @return Returns xyz position in a 3-element array.
*/
void TXftAnaModule::track_Prop(float radius, float curvature, float PhiZeroIn, 
			       float cotan, float D_0, float Z_0, float* xyz) {
  float epsilon = 1.0;        // epsilon is outgoing so is 1
  track_Prop(epsilon,radius,curvature,PhiZeroIn,cotan,D_0,Z_0,xyz);
}

/**
* Calculates the x,y,z position of a charged track passing through the
* COT, at a given radius from the origin.
* @param radius is the radius from x=0,y=0 to determine the x,y,z
* @param curvature is the half curvature for the charged track
* @param PhiZeroIn is the phi at the origin for the charged track
* @param cotan is the  cot(theta) for the charged track
* @param D_O is the impact parameter for the charged track
* @param Z_O is the z position at the distance of closest approach
*            to x=0,y=0
* @return Returns xyz position in a 3-element array.
*/
void TXftAnaModule::track_Prop(float epsilon, float radius, float curvature, 
			       float PhiZeroIn, float cotan, float D_O, float Z_0, float* xyz) {
  //
  //		float xyz[] = new float[3];
  //
  // calculate some constants on entry
  //		float epsilon = 1.0;        // epsilon is outgoing so is 1
  //	printf("in track_prop rad %f d0 %f \n",radius,D_O);
  float B = curvature* sqrt((radius*radius-D_O*D_O)/(1.0+2.0*curvature*D_O));
  float U_O = cos(PhiZeroIn);
  float V_O = sin(PhiZeroIn);
  float rho = 2.0*curvature;
  float XO = -D_O*V_O;
  float YO = D_O*U_O;			
  
  float cosPhiZero = cos(PhiZeroIn);
  float sinPhiZero = sin(PhiZeroIn);
  
  float Radius_of_Arc = fabs(1.0/(2*curvature)); 	
  //  float X_C = -(1/(2*curvature)+D_O)*sinPhiZero;  // X_C calculated here
  //float Y_C =  (1/(2*curvature)+D_O)*cosPhiZero;  // Y_C calculated here
  //			
  xyz[0] = XO + (U_O*2.0*epsilon*B*sqrt(1.0-B*B)-V_O*2.0*B*B)/rho;
  xyz[1] = YO + (V_O*2.0*epsilon*B*sqrt(1.0-B*B)+U_O*2.0*B*B)/rho;
  xyz[2] = Z_0 + ((cotan)*asin(B))/curvature;
  
}

//------------------------------------------------------------------------
// This function returns the phi value in 0 to 2*PI range.
//------------------------------------------------------------------------
float TXftAnaModule::cdfphi(float phi) {
  if (phi < 0) {
    while (phi < 0) {
      phi = phi + twopi;
    }
  } else {
    while (phi > twopi) {
      phi = phi - twopi;
    }
  }
  return phi;
}

//----------------------------------------------------------------
// CONVERT A PT VALUE IN GeV/c TO THE NEAREST XFT PT BIN
//----------------------------------------------------------------
int TXftAnaModule::nearestPtBin(float ptVal) {
  
  float pt = fabs(ptVal);
  int binind = 128;
  float mindelpt = 1000000.0;
  float delpt;
  if (pt > 1.50) {
    for (int j = 0; j<48; j++) {
      delpt = fabs(pt-ptbins[j]);
      //			printf("delpt %f mindelpt %f bin %d \n",delpt,mindelpt,binind);
      if (delpt < mindelpt) {
	mindelpt = delpt;
	binind = j;
      }
    }
  }
  if (ptVal < 0.0) {
    binind = 47 - binind;
  } else {
    binind = 48 + binind;
  }
  
  return binind;
  
}

//-----------------------------------------------------------------------------------------------
//DIVIDE VARIOUS HISTOGRAMS FOR EFFICIENCY PLOTS
//-----------------------------------------------------------------------------------------------
void TXftAnaModule::effCalc (TH1F* numHist, TH1F* denomHist, TH1F* effHist) 
{
  
  // Loop over the Pt Fnd and Gen plots to determine efficiency
  for(Int_t i=0; i<(denomHist->GetNbinsX()+1); i++) 
    {
      float effFnd=0.0;
      float efferr=0.0;
      float gennum = denomHist->GetBinContent(i);
      float numobs = numHist->GetBinContent(i);
      
      
      if(gennum>0.5) {
	effFnd = numobs/gennum;
	efferr = sqrt(effFnd*(1.0-effFnd)/gennum);
      }
      
      effHist->SetBinContent(i,effFnd);
      effHist->SetBinError(i,efferr);
    }
  
}

//-----------------------------------------------------------------------------------------
//TRY TO REMOVE DUPLICATE TRACKS A LA XTRP
//-----------------------------------------------------------------------------------------	
bool TXftAnaModule::xtrp_duplicate_removal(int chip1,int miniphi1,int board1,int crate1,
					   int chip2,int miniphi2,int board2,int crate2)
{
  bool dupl = false;
  //
  //SAME CRATE
  //Check if on same board - chips 5 &6
  if(crate1 == crate2){
    if(board1 == board2){
      if(chip1 == 5 && chip2 == 6){
	if(miniphi1 >5 && miniphi2 <2) dupl=true;
	
      }else if(chip2 == 5 && chip1 == 6){
	if(miniphi2 >5 && miniphi1 <2) dupl=true;
      }
    }
    //						
    //Check if on neighbouring boards
    
    if( (board1 == (board2-1)) &&
	(chip1 == 11)  && (chip2 == 0) &&
	(miniphi1 > 5) && (miniphi2 < 2))	dupl=true;
    
    if( (board2 == (board1-1)) &&	
	(chip2 == 11) && (chip1 == 0)  &&
	(miniphi2 >5) && (miniphi1 <2) ) dupl=true;
    
  }else{
    //DIFFERENT CRATES
    //Check if on neighbouring boards in neigbouring crates
    
    if( (board1 == 7) && (board2 == 0) ){
      if( crate2 == ((crate1 + 1)%3) ){
	if(chip1==11 && chip2==0){
	  if(miniphi1 >5 && miniphi2 <2) dupl=true;
	}
      }
    }else if( (board2 == 7) && (board1==0) ){
      if( crate1 == ((crate2 + 1)%3) ){
	if(chip2==11 && chip1==0){
	  if(miniphi2 >5 && miniphi1 <2) dupl=true;
	}
      }
    }
  }
  return dupl;
}


void TXftAnaModule::cotGeom(char* pathName)
{
  // Setup the Input file Stream
  // ---------------------------
  std::ifstream f1(pathName);


  // Read in the header pixels per cell for each layer
  // -------------------------------------------------
  float rtemp1,rtemp2,rtemp3,rtemp4,rtemp5;
  int itemp1,itemp2,itemp3;
  for(int sl=0; sl < 8; sl++)
    {
      for(int wire=0; wire<12; wire++) {
      	f1 >> itemp1
	   >> itemp2
	   >> itemp3
	   >> rtemp1
	   >> rtemp2
	   >> rtemp3
	   >> cotRadius[sl][wire]
	   >> cotPhi[sl][wire]
	   >> rtemp4
	   >> rtemp5;
	printf("sl %i wire %i radius %f phi %f\n",sl,wire,
	       cotRadius[sl][wire],
	       cotPhi[sl][wire]);
	//	      printf(" %i %i %i %f %f %f %f %f %f %f \n",
	//         sl,cell,wire,x,y,r,lrad,phi,dlwr,delphi);
      	
        }
    }
  f1.close();
}//End method
// ********************************************
// **             nearTracks                **
// ********************************************
int TXftAnaModule::nearTracks( 
			      float curv1, float phi1, float cotan1, float d1, float z1, 
			      float curv2, float phi2, float cotan2, float d2, float z2)
{
  int near = 0;
  float xyz[3];
  for (int isl=1; isl<8; isl=isl+2) {
    //
    // Find the where track goes in this layer
    track_Prop(rsns[isl], 
	       curv1, phi1, cotan1, d1, z1, xyz);
    if (fabs(xyz[2]) < 150.0)
      {
      	
	float phi_int1 = cdfphi(atan2(xyz[1],xyz[0]));
	//
	// Find the where track goes in this layer
	track_Prop(rsns[isl], 
		   curv2, phi2, cotan2, d2, z2, xyz);
	if (fabs(xyz[2]) < 150.0)
	  {
	    
	    float phi_int2 = cdfphi(atan2(xyz[1],xyz[0]));
	    //
	    // How close are these?
	    float phidiff = cdfphi(phi_int1 - phi_int2);
            if (phidiff < 5*dclw[isl]) near++;
	  }
      }
  }
  
  return near;
}



int TXftAnaModule::returnCell (int layer, int wire, 
			       float cot_curv, 
			       float cot_phi0, float cot_cotan, 
			       float cot_d0, float cot_z0)
{
  
//
  // Local vars
  float phi_int;
  int cell_hit;
  float xyz[3];
  //
  // Find the cell this track goes through
  track_Prop(cotRadius[layer][wire], 
	     cot_curv, cot_phi0, cot_cotan, 
	     cot_d0, cot_z0, xyz);
  //	printf("returncell rad, curv, phi, cot, d0, z0, z %f %f %f %f %f %f %f\n",
  //   		cotRadius[layer][wire], 
  //						cot_curv, cot_phi0, cot_cotan, 
  //						cot_d0, cot_z0, xyz[2]);
  if (fabs(xyz[2]) < 150.0)
    {
      phi_int = cdfphi(atan2(xyz[1],xyz[0]));
      phi_int = phi_int - cotPhi[layer][wire];
      cell_hit = cdfphi(phi_int) / dclw[layer];
    } else {
      cell_hit = -999;
   }
  //   printf(" phi_int %f dclw %f cell %i \n",phi_int,dclw[layer],cell_hit);
  
  return cell_hit;
}
double TXftAnaModule::getPhiAtRadius(float ptXFT, float phiXFT, float Radius)
{
  double cval = konst / ptXFT;
  double beam_x = -0.1035;
  double beam_y = 0.386;
  //	double beam_x = 0.0;
  //	double beam_y = 0.0;
  
  // Get Distance from offset beam position to ASL3 point.
  
  double x = 105.575*cos(phiXFT) - beam_x;
  double y = 105.575*sin(phiXFT) - beam_y;
  double rdis = sqrt(x*x+y*y);
  
  double phiBeam = 3.14159+atan2(-y,-x) - asin(rdis*cval);
  
  // Now create a Helix
  double ch = 1.0;
  if(cval<0.0) ch = -1.0;
  //	printf(" ptXFT %f phiXFT %f Rdis %f phiBeam %f Chg %f \n",ptXFT,phiXFT,rdis,phiBeam,ch);
  Helix *xHel = new Helix(fabs(ptXFT),phiBeam,0.0,beam_x,beam_y,0.0,ch,1.4116);	
  
  // Now get the Phi at the outer silicon layer
  double phiprop = xHel->getPhi0();
  if(Radius>0.) phiprop = xHel->getPhiAtR(Radius);	
  
  //	phiprop = phiBeam;

  return phiprop;
}