//==================================================================
//
//	TEleCand.cc
//
// PURPOSE: Example of how to read in new top electron blocks
//
// DESCRIPTION:
// Make histograms of electron variables for central and plug electrons
// Make N-1 plots (all TEleCand cuts applied except on variable plotted) for central electrons
// Calculate invariant mass of tight-tight and tight-loose electron pairs.
//
// INPUT is your Stntuple file - please specify name and location in Stntuple/top/ele_cand.C
// OUTPUT is called ele_cand.hbk and ia a file containing lots of histograms!
//
// INSTRUCTIONS:
// You can analyse your Stntuple file with TEleCand like this
// .x Stntuple/top/ele_cand.C
//
// You can display a histogram from the output file like this
// .x Stntuple/top/ele_example.C("ele_cand.hbk");
//
// You can make a postscript file with all 1D histograms from the output file  like this
// .x Stntuple/top/ele_display.C("ele_cand.hbk");
// .! gs electron.ps
// 
//	Version 0.0:      P. Koehn     5/5/01
//	Version 0.9:      R. Hughes    7/23/01
//	Version 1.0:      E. Thomson   2/6/02
//
//==================================================================
#if !defined (__CINT__) || defined (__MAKECINT__)
#include "TF1.h"
#include <iostream>
#include <iomanip>
#include <algorithm>

#include "Stntuple/loop/TStnAna.hh"

#endif
#include "TEleCand.hh"

//Add an electron directory to an existing histogram file
//_____________________________________________________________________________
TEleCand::TEleCand(TFile* file, const char* name, const char* title): TStnModule(name,title){
  //
  // create a subdirectory in existing file
  histFile = file;
  histDir = name;
  dir = histFile->mkdir(histDir);
  //
  // Initialize cuts
  init();
}
//Make a new histogram file, add an electron directory to it
//_____________________________________________________________________________
TEleCand::TEleCand(const char* name, const char* title): TStnModule(name,title){
  //
  // create a subdirectory in new file
  histFile = new TFile("TEleCand.hbk","RECREATE","test file");
  histDir = name;
  dir = histFile->mkdir(histDir);
  //
  // Initialize cuts
  init();
}

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

//Initialise electron cuts
//_____________________________________________________________________________
void TEleCand::init() {
  //
  printf("----- Init: ---- %s\n",GetName());
  NElectronMax = 20;
  totEvents = 0;
  eleEvents = 0;

  for(int i=0;i<20;i++) nPassCut[i] = 0;
  nameCut[0] = "All enabled cuts";
  nameCut[1] = "Et";
  nameCut[2] = "TrackPt";
  nameCut[3] = "EOverP";
  nameCut[4] = "DelX";
  nameCut[5] = "DelZ";
  nameCut[6] = "Lshr";
  nameCut[7] = "HadEm";
  nameCut[8] = "Isol";
  nameCut[9] = "Chi2Strip";

  nCem = 0;
  nZees = 0;
  //
  // Default is to make histograms
  enableEleHist(true);
  
  //Use TopEventModule tight and loose collections?
  enableEleTopEvent(false);

  //Use pre TopEventModule code?
  enableOldCode(false);

  //
  // Define the default cut values for high PT electrons
  enableEleEtCut(false);
  setEleEtCut(20.0);

  enableElePtCut(false);
  setElePtCut(15.0);

  enableEleEOverPCut(false);
  setEleEOverPCut(2.0);

  enableEleDelXCut(false);
  setEleDelXCut(3.0);

  enableEleDelZCut(false);
  setEleDelZCut(3.0);

  enableEleLshrCut(false);
  setEleLshrCut(0.2);

  enableEleHadEmCut(false);
  setEleHadEmCut(0.1);

  enableEleIsolCut(false);
  setEleIsolCut(0.1);

  enableEleChi2StripCut(false);
  setEleChi2StripCut(10.0);
//
// Turn off this cut - not defined yet 11/27/01
  enableEleTisolCut(false);
  setEleTisolCut(0.1);
//
// Turn off this cut - not defined yet 11/27/01
  enableEleDelZvrtCut(false);
  setEleDelZvrtCut(5.0);
//
// Turn off this cut - not defined yet 11/27/01
  enableEleZvrtCut(false);
  setEleZvrtCut(60.0);
//
// Remove conversion electrons
  enableEleRemoveConv(false);

}
//_____________________________________________________________________________
int TEleCand::BeginJob() {

  printf("----- BeginJob: ---- %s\n",GetName());
  printf("  Electron Et cut enabled       %d value  %f\n",isEnableEtCut(),getEleEtCut());
  printf("  Electron Pt cut enabled       %d value  %f\n",isEnablePtCut(),getElePtCut());
  printf("  Electron EOverP cut enabled   %d value  %f\n",isEnableEOverPCut(),getEleEOverPCut());
  printf("  Electron DexX cut enabled     %d value  %f\n",isEnableDelXCut(),getEleDelXCut());
  printf("  Electron DelZ cut enabled     %d value  %f\n",isEnableDelZCut(),getEleDelZCut());
  printf("  Electron Lshr cut enabled     %d value  %f\n",isEnableLshrCut(),getEleLshrCut());
  printf("  Electron Had/Em cut enabled   %d value  %f\n",isEnableHadEmCut(),getEleHadEmCut());
  printf("  Electron Isol cut enabled     %d value  %f\n",isEnableIsolCut(),getEleIsolCut());
  printf("  Electron Trk Isol cut enabled %d value  %f\n",isEnableTisolCut(),getEleTisolCut());
  printf("  Electron Chi2Strip cut enabled     %d value  %f\n",isEnableChi2StripCut(),getEleChi2StripCut());
  printf("  Electron DelZvrt cut enabled  %d value  %f\n",isEnableDelZvrtCut(),getEleDelZvrtCut());
  printf("  Electron Zvrt cut enabled     %d value  %f\n",isEnableZvrtCut(),getEleZvrtCut());


  // register the Track and Jet data blocks
  
  pHeaderBlock = (TStnHeaderBlock*) fAna->RegisterDataBlock("HeaderBlock","TStnHeaderBlock");
  if (! pHeaderBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","HeaderBlock");
    fEnabled = 0;
  }
  pTrackBlock = (TStnTrackBlock*) fAna->RegisterDataBlock("TrackBlock","TStnTrackBlock");
  if (! pTrackBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","TrackBlock");
    fEnabled = 0;
  }
  pJetBlock = (TStnJetBlock*) fAna->RegisterDataBlock("JetBlock","TStnJetBlock");
  if (! pJetBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","JetBlock");
  }

  // register the standard electron block
  pEleBlock = (TStnElectronBlock*) fAna->RegisterDataBlock("ElectronBlock","TStnElectronBlock");
  if (! pEleBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","EleBlock");
    fEnabled = 0;
  }

  // register blocks for tight/loose electrons defined by TopMods/src/TopEventModule.cc
  if(isEnableTopEvent()){
    
    pTopTEleBlock = (TStnElectronBlock*) fAna->RegisterDataBlock("TopTElectronBlock","TStnElectronBlock");
    if (! pTopTEleBlock) {
      printf(" >>> branch *** %s *** doesn't exist \n","TopTElectronBlock");
      fEnabled = 0;
    }  
    
    pTopLEleBlock = (TStnElectronBlock*) fAna->RegisterDataBlock("TopLElectronBlock","TStnElectronBlock");
    if (! pTopLEleBlock) {
      printf(" >>> branch *** %s *** doesn't exist \n","TopLElectronBlock");
      fEnabled = 0;
    }  

    //The top summary block
    pTopSummaryBlock = (TStnTopSummaryBlock*) fAna->RegisterDataBlock("TopSummaryBlock","TStnTopSummaryBlock");
    if (! pTopSummaryBlock) {
      printf(" >>> branch *** %s *** doesn't exist \n","TopSummaryBlock");
      fEnabled = 0;
    }

  }

  // Book the histograms
  if (isEnableHist()) {
    BookHistograms("default");
    if (isEnableTopEvent()){
      BookHistograms("tight");
      BookHistograms("loose");
      BookZBosonHistograms("ZBoson");
    }
  }

  return 0;
}

//_____________________________________________________________________________
int TEleCand::BeginRun() {
  printf("----- BeginRun: ---- %s\n",GetName());
  return 0;
}

//_____________________________________________________________________________
void TEleCand::BookHistograms(const char* subdirName) {
  printf("----- BookHistograms: ---- %s\n",GetName());
  
  histFile->cd(histDir); 
  //make a new subdirectory
  TDirectory* dirdef = gDirectory->mkdir(subdirName);
  gDirectory->cd(subdirName);
  
  new TH1F("DetCode","Detector code",3,-0.5,2.5);
  // CENTRAL

  //Best estimate of 4-momentum
  new TH1F("ncem_E","Energy",50, 0.0, 100.0);
  new TH1F("ncem_Et","Electron Transverse energy",50, 0.0, 100.0);
  new TH1F("ncem_Phi","Electron Phi",100, 0.0, 6.2832);  
  new TH1F("ncem_Eta","Electron Eta",80,-4.0, 4.0);  

  //Seed tower location of em cluster
  new TH1F("ncem_SeedPhi","EM Seed Phi",24, -0.5, 23.5);
  new TH1F("ncem_SeedEta","EM Seed Eta",52, -0.5, 51.5);
  new TH1F("ncem_Wedge","EM Seed Wedge (accessor for central)",48, -0.5, 47.5);  
  new TH1F("ncem_Side","EM Seed Side",2,-0.5, 1.5);  

  //Em cluster angles
  new TH1F("ncem_EmPhi","EM Cluster Phi",100, 0.0, 6.2832);  
  new TH1F("ncem_EmDetEta","EM Cluster Detector Eta",80,-4.0, 4.0);  
  new TH1F("ncem_EmZv","EM Cluster z vertex",100, -200.0, 200.0);
  new TH1F("ncem_EmEvEta","EM Cluster Eta using z vertex",80,-4.0, 4.0);  

  //Em cluster energy
  new TH1F("ncem_EmEt","Em Transverse Energy",50, 0.0, 100.0);
  new TH1F("ncem_HadEt","Had Transverse Energy",100, 0.0, 20.0);
  new TH1F("ncem_EmE","Em Energy",50, 0.0, 100.0);
  new TH1F("ncem_HadE","Had Energy",100, 0.0, 20.0);

  //ID variables and isolation
  new TH1F("ncem_HadEm","HAD energy / EM energy",50, 0.0, 0.5);
  new TH1F("ncem_Lshr","Lshr",120, -3.0, 3.0);
  new TH1F("ncem_Isol","Isolation",200, 0.0, 2.0);

  //Track quantities
  new TH1F("ncem_EOverP","E/P",100, 0.0, 5.0);  
  new TH1F("ncem_Pt","Transverse momentum",50, 0.0, 100.0);  
  new TH1F("ncem_Charge","Track charge",3, -1.5, 1.5);
  new TH1F("ncem_Z0","Track Z0",100, -200.0, 200.0);
  new TH1F("ncem_D0","Track D0",120, -3.0, 3.0);
  new TH1F("ncem_Phi0","Track Phi0",100, 0.0, 6.2832);  
  new TH1F("ncem_Eta0","Track Eta",80,-4.0, 4.0);  

  new TH1F("ncem_Chi2Track","Track Chi2",50,0.0, 100.0);  
  new TH1F("ncem_Chi2Cot","Track COT Chi2",50,0.0, 100.0);  
  new TH1F("ncem_Chi2Svx","Track SVX Chi2",50,0.0, 100.0);  
  new TH1F("ncem_nDof","Track Fit number of degrees of freedom (assumes npar=5)",110,-0.5, 109.5);  
  new TH1F("ncem_Chi2perDof","Track Chi2 per degree of freedom (assumes npar=5)",50,0.0, 10.0);  

  new TH1F("ncem_NAxHits","Track Axial Hits",50,-0.5, 49.5);  
  new TH1F("ncem_NStHits","Track Stereo Hits",50,-0.5, 49.5);  
  new TH1F("ncem_NTotHits","Track COT Hits",100,-0.5, 99.5);  

  new TH1F("ncem_NSvxHits","Track SVX Hits",20,-0.5, 19.5);  
  new TH1F("ncem_NSvxRPhiHits","Track SVX RPhi Hits",20,-0.5, 19.5);  
  new TH1F("ncem_NSvxSASHits","Track SVX Small Angle Stereo Hits",20,-0.5, 19.5);  

  new TH1F("ncem_ClusterSize","Cluster Size",20,-0.5, 19.5);  
  new TH1F("ncem_NTracks","Number of matched tracks",20,-0.5, 19.5); 
 
  new TH1F("ncem_NCesStrip","Number of matched CES strips",20,-0.5, 19.5);  
  new TH1F("ncem_ECesStrip","CES strip energy",50,0.0, 100.0);  
  new TH1F("ncem_ZCes","Z(CES,strip) ",65,0.0, 260.0); //4cm
  new TH1F("ncem_DelZ","DeltaZ",50,-10.0, 10.0);
  new TH1F("ncem_Chi2Strip","Strip Chi2",50, 0.0, 50.0);

  new TH1F("ncem_NCesWire","Number of matched CES wires",20,-0.5, 19.5);  
  new TH1F("ncem_ECesWire","CES wire energy",50,0.0, 100.0);  
  new TH1F("ncem_XCes","X(CES,wire) ",30,-30.0, 30.0); //2cm
  new TH1F("ncem_DelX","DeltaX",50,-10.0, 10.0);
  new TH1F("ncem_Chi2Wire","Wire Chi2",50, 0.0, 50.0);

  //2D plots
  new TH2F("ncem_Isol_Et","Isolation vs Transverse Energy",50, 0.0, 100.0, 50, 0.0, 0.5);
  new TH2F("ncem_HadEm_E","HAD/EM Energy vs Energy",50, 0.0, 150.0, 50, 0.0, 0.125);
  new TH2F("ncem_EOverP_E","E/P vs Energy",50, 0.0, 100.0, 50, 0.0, 5.);
  new TH2F("ncem_DelX_Et","Delta X vs Transverse Energy",50, 0.0, 100.0, 50, -10.0, 10.0);
  new TH2F("ncem_DelZ_Et","Delta Z vs Transverse Energy",50, 0.0, 100.0, 50, -10.0, 10.0);
  new TH2F("ncem_DelX_Wedge","Delta X vs Wedge (W,E)",48, -0.5, 47.5, 50, -10.0, 10.0);
  new TH2F("ncem_DelZ_Wedge","Delta Z vs Wedge (W,E)",48, -0.5, 47.5, 50, -10.0, 10.0);
  new TH2F("ncem_Phi0_EmPhi","Track Phi0  vs EmCluster Phi",50, 0.0, 6.2832, 50, 0.0, 6.2832);    
  new TH2F("ncem_DelEStrip_E","(E_Strip - E_EM) vs Energy",50, 0.0, 100.0, 50, -25.0, 25.0);
  new TH2F("ncem_DelEWire_E","(E_Wire - E_EM) vs Energy",50, 0.0, 100.0, 50, -25.0, 25.0);

  // PLUG

  //Best estimate of 4-momentum
  new TH1F("npem_E","Energy",100, 0.0, 200.0);
  new TH1F("npem_Et","Electron Transverse energy",50, 0.0, 100.0);
  new TH1F("npem_Phi","Electron Phi",100, 0.0, 6.2832);  
  new TH1F("npem_Eta","Electron Eta",80,-4.0, 4.0);  

  //Seed tower location of em cluster
  new TH1F("npem_SeedPhi","EM Seed Phi",48, -0.5, 47.5);
  new TH1F("npem_SeedEta","EM Seed Eta",52, -0.5, 51.5);
  new TH1F("npem_Wedge","EM Seed Wedge ",96, -0.5, 95.5);  
  new TH1F("npem_Side","EM Seed Side",2,-0.5, 1.5);  

  //Em cluster angles
  new TH1F("npem_EmPhi","EM Cluster Phi",100, 0.0, 6.2832);  
  new TH1F("npem_EmDetEta","EM Cluster Detector Eta",80,-4.0, 4.0);  
  new TH1F("npem_EmZv","EM Cluster z vertex",100, -200.0, 200.0);
  new TH1F("npem_EmEvEta","EM Cluster Eta using z vertex",80,-4.0, 4.0);  

  //Em cluster energy
  new TH1F("npem_EmEt","Em Transverse Energy",50, 0.0, 100.0);
  new TH1F("npem_HadEt","Had Transverse Energy",100, 0.0, 20.0);
  new TH1F("npem_EmE","Em Energy",100, 0.0, 200.0);
  new TH1F("npem_HadE","Had Energy",100, 0.0, 20.0);

  //ID variables and isolation
  new TH1F("npem_HadEm","HAD energy / EM energy",50, 0.0, 0.5);
  new TH1F("npem_Chi2Three","3x3 Chi2",50, 0.0, 50.0);
  new TH1F("npem_Isol","Isolation",200, 0.0, 2.0);

  //Track quantities
  new TH1F("npem_EOverP","E/P",100, 0.0, 5.0);  
  new TH1F("npem_Pt","Transverse momentum",50, 0.0, 100.0);  
  new TH1F("npem_Charge","Track charge",3, -1.5, 1.5);
  new TH1F("npem_Z0","Track Z0",100, -200.0, 200.0);
  new TH1F("npem_D0","Track D0",120, -3.0, 3.0);
  new TH1F("npem_Phi0","Track Phi0",100, 0.0, 6.2832);  
  new TH1F("npem_Eta0","Track Eta",80,-4.0, 4.0);  

  new TH1F("npem_Chi2Track","Track Chi2",50,0.0, 100.0);  
  new TH1F("npem_Chi2Cot","Track COT Chi2",50,0.0, 100.0);  
  new TH1F("npem_Chi2Svx","Track SVX Chi2",50,0.0, 100.0);  
  new TH1F("npem_nDof","Track Fit number of degrees of freedom (assumes npar=5)",110,-0.5, 109.5);  
  new TH1F("npem_Chi2perDof","Track Chi2 per degree of freedom (assumes npar=5)",50,0.0, 10.0);  

  new TH1F("npem_NAxHits","Track Axial Hits",50,-0.5, 49.5);  
  new TH1F("npem_NStHits","Track Stereo Hits",50,-0.5, 49.5);  
  new TH1F("npem_NTotHits","Track COT Hits",100,-0.5, 99.5);  

  new TH1F("npem_NSvxHits","Track SVX Hits",20,-0.5, 19.5);  
  new TH1F("npem_NSvxRPhiHits","Track SVX RPhi Hits",20,-0.5, 19.5);  
  new TH1F("npem_NSvxSASHits","Track SVX Small Angle Stereo Hits",20,-0.5, 19.5);  

  new TH1F("npem_ClusterSize","Cluster Size",20,-0.5, 19.5);  
  new TH1F("npem_NTracks","Number of matched tracks",20,-0.5, 19.5);  

  new TH1F("npem_NPesClusters","Number of matched PES clusters",20,-0.5, 19.5);  
  new TH1F("npem_EPesU","PES U energy",100,0.0, 200.0);  
  new TH1F("npem_EPesV","PES V energy",100,0.0, 200.0);  
  new TH1F("npem_PesPhi","PES Phi",50, 0.0, 6.2832);  
  new TH1F("npem_PesEta","PES Eta",80,-4.0, 4.0);  
  new TH1F("npem_DelR","Delta R PES-PEM ",50,0.0, 0.5);  

  new TH2F("npem_Isol_Et","Isolation vs Transverse Energy",50, 0.0, 100.0, 50, 0.0, 0.5);
  new TH2F("npem_HadEm_E","HAD/EM Energy vs Energy",50, 0.0, 200.0, 50, 0.0, 0.125);
  new TH2F("npem_EOverP_E","E/P vs Energy",50, 0.0, 200.0, 50, 0.0, 5.);
  new TH2F("npem_Chi2Three_Et","3x3 Chi2 vs Transverse Energy",50, 0.0, 100.0, 50, 0.0, 50.0);
  new TH2F("npem_DelR_Et","Delta R vs Transverse Energy",50, 0.0, 100.0, 50, 0.0, 1.0);
  new TH2F("npem_DelR_Wedge","Delta R vs Wedge (W,E)",48, -0.5, 47.5, 50, 0.0, 1.0);
  new TH2F("npem_Phi0_EmPhi","Track Phi0  vs EmCluster Phi",50, 0.0, 6.2832, 50, 0.0, 6.2832);    
  new TH2F("npem_PesPhi_EmPhi","Pes Phi  vs EmCluster Phi",50, 0.0, 6.2832, 50, 0.0, 6.2832);    
  new TH2F("npem_PesEta_EmEta","Pes Eta  vs EmCluster Eta",50, -4.0, 4.0, 50, -4.0, 4.0);    
  new TH2F("npem_DelEU_E","(E_U - E_EM) vs Energy",50, 0.0, 200.0, 50, -25.0, 25.0);
  new TH2F("npem_DelEV_E","(E_V - E_EM) vs Energy",50, 0.0, 200.0, 50, -25.0, 25.0);


  gDirectory->pwd();
  gDirectory->ls();

}
//_____________________________________________________________________________
void TEleCand::BookZBosonHistograms(const char* subdirName) {
  printf("----- BookZBosonHistograms: ---- %s\n",GetName());
  
  //Book standard detector histograms
  BookHistograms(subdirName);
  
  //Book extra-special Z boson histograms
  pairMass = new TH1F("pairMass","Pair Mass",100, 0.0, 200.0);
    
  gDirectory->pwd();
  gDirectory->ls();

}

//==================================
//       Event(int ientry)
//==================================
int TEleCand::Event(int ientry) {
  
  if (getDebugLevel() > 0) {
    printf(" ---- TEleCand ientry = %7i\n",ientry);
  }

  // Counter for number of events processed so far
  totEvents++;
  
  // Get the header info
  pHeaderBlock->GetEntry(ientry);
  int iev = pHeaderBlock->EventNumber();
  int irun = pHeaderBlock->RunNumber();
  if (ientry%100 == 0)  printf("  Events processed %d current event is %d run %d \n",totEvents,iev,irun);


  // Get the top summary info
  pTopSummaryBlock->GetEntry(ientry);
  //pTopSummaryBlock->Print();

  // Get the jet info
  //pJetBlock->GetEntry(ientry);
  //int njets = pJetBlock->NJets();
  
  // Get the tracking info
  pTrackBlock->GetEntry(ientry);
  int ntracks = pTrackBlock->NTracks();

  // Get access to the conversion electron module
  //TStnAna * TestAna = TStnModule::GetAna();
  //TConvEleCand* eleConvCand = (TConvEleCand*) TestAna->TStnAna::GetModule("EleConv");
 
  // Get the electron blocks
  pEleBlock->GetEntry(ientry);
  int nelectrons = pEleBlock->NElectrons();
  
  int nTopTelectrons = 0;
  int nTopLelectrons = 0;
  if(isEnableTopEvent()){
    pTopTEleBlock->GetEntry(ientry);
    pTopLEleBlock->GetEntry(ientry);
    nTopTelectrons = pTopTEleBlock->NElectrons();
    nTopLelectrons = pTopLEleBlock->NElectrons();
  }
  
  //
  // Fill the histograms 
  if (isEnableHist()) {
    
    // Fill histograms for default electrons 
    if(nelectrons >0){
      GetHistogramPointers("default");      
      for (int i=0; i<nelectrons; i++) {
	TStnElectron* electron = pEleBlock->Electron(i);
	FillHistograms(electron);
      }
    }
    
    //
    // Fill histograms for tight electrons
    if(nTopTelectrons >0){      
      GetHistogramPointers("tight");
      for (int i=0; i<nTopTelectrons; i++) {
	TStnElectron* electron = pTopTEleBlock->Electron(i);
	FillHistograms(electron);
      }
    }
    //
    // Fill histograms for loose electrons
    if(nTopLelectrons >0){      
      GetHistogramPointers("loose");
      for (int i=0; i<nTopLelectrons; i++) {
	TStnElectron* electron = pTopLEleBlock->Electron(i);
	FillHistograms(electron);
      }
    }
  }

  
  //Z boson efficiency study 
  if(isEnableTopEvent()){
    GetHistogramPointers("ZBoson");
    
    int nZCandidates = 0;
    // Now loop over the tight electrons
    //For each tight electron, calculate invariant mass with each of the other tight/loose electrons
    for (int i=0; i<nTopTelectrons; i++) {
      TStnElectron* tight = pTopTEleBlock->Electron(i); 
      
      //Loop over other Tight electrons
      for(int j=i+1; j<nTopTelectrons;j++){
	TStnElectron* ele2 = pTopTEleBlock->Electron(j);
	if(isZBoson(tight,ele2)) nZCandidates++;
      }	  
      
      //Loop over loose electrons
      for(int j=0; j<nTopLelectrons;j++){
	TStnElectron* ele2 = pTopLEleBlock->Electron(j);
	
	//Require that loose electron passes isolation cut too
	if(getDebugLevel()>1) cout << " Isolation " << ele2->Iso()/ele2->Et() 
				   << " Cut" << getEleIsolCut() << endl;
	
	if((ele2->Iso()/ele2->Et()) < getEleIsolCut() ){
	  if(isZBoson(tight,ele2)) nZCandidates++;
	}
      }	  
    }
    
    if(nZCandidates > 1) cout <<  "WARNING More than 1 Z candidate  " << nZCandidates << " in this event "<< iev << endl;
  }//Endif for isEnableTopEvent

  //
  //Old Code before TopEventModule existed
  //Make TEleCand cuts on Stntuple standard electrons, set flags for good electrons and link to closest jet
  if(isEnableOldCode())oldCodeElectron();
  
  return 0;		       
}
//==================================
//       GetHistogramPointers()
//==================================
bool TEleCand::GetHistogramPointers(const char* subdirName) {

  histFile->cd(histDir);
  gDirectory->cd(subdirName);   
  if(getDebugLevel()>2) gDirectory->pwd();
  if(getDebugLevel()>3) gDirectory->ls();

  DetCode= (TH1F*) gDirectory->Get("DetCode");
  ncem_E= (TH1F*) gDirectory->Get("ncem_E");
  ncem_Et= (TH1F*) gDirectory->Get("ncem_Et");
  ncem_Phi= (TH1F*) gDirectory->Get("ncem_Phi");
  ncem_Eta= (TH1F*) gDirectory->Get("ncem_Eta");

  ncem_SeedPhi= (TH1F*) gDirectory->Get("ncem_SeedPhi");
  ncem_SeedEta= (TH1F*) gDirectory->Get("ncem_SeedEta");
  ncem_Wedge= (TH1F*) gDirectory->Get("ncem_Wedge");
  ncem_Side= (TH1F*) gDirectory->Get("ncem_Side");

  ncem_EmPhi= (TH1F*) gDirectory->Get("ncem_EmPhi");
  ncem_EmDetEta= (TH1F*) gDirectory->Get("ncem_EmDetEta");
  ncem_EmZv= (TH1F*) gDirectory->Get("ncem_EmZv");
  ncem_EmEvEta= (TH1F*) gDirectory->Get("ncem_EmEvEta");

  ncem_EmEt= (TH1F*) gDirectory->Get("ncem_EmEt");
  ncem_HadEt= (TH1F*) gDirectory->Get("ncem_HadEt");
  ncem_EmE= (TH1F*) gDirectory->Get("ncem_EmE");
  ncem_HadE= (TH1F*) gDirectory->Get("ncem_HadE");

  ncem_HadEm= (TH1F*) gDirectory->Get("ncem_HadEm");
  ncem_Lshr= (TH1F*) gDirectory->Get("ncem_Lshr");
  ncem_Isol= (TH1F*) gDirectory->Get("ncem_Isol");

  ncem_EOverP= (TH1F*) gDirectory->Get("ncem_EOverP");
  ncem_Pt= (TH1F*) gDirectory->Get("ncem_Pt");
  ncem_Charge= (TH1F*) gDirectory->Get("ncem_Charge");
  ncem_Z0= (TH1F*) gDirectory->Get("ncem_Z0");  
  ncem_D0= (TH1F*) gDirectory->Get("ncem_D0");  
  ncem_Phi0= (TH1F*) gDirectory->Get("ncem_Phi0");  
  ncem_Eta0= (TH1F*) gDirectory->Get("ncem_Eta0");  

  ncem_Chi2Track= (TH1F*) gDirectory->Get("ncem_Chi2Track");
  ncem_Chi2Cot= (TH1F*) gDirectory->Get("ncem_Chi2Cot");
  ncem_Chi2Svx= (TH1F*) gDirectory->Get("ncem_Chi2Svx");
  ncem_nDof= (TH1F*) gDirectory->Get("ncem_nDof");
  ncem_Chi2perDof= (TH1F*) gDirectory->Get("ncem_Chi2perDof");

  ncem_NAxHits= (TH1F*) gDirectory->Get("ncem_NAxHits");
  ncem_NStHits= (TH1F*) gDirectory->Get("ncem_NStHits");
  ncem_NTotHits= (TH1F*) gDirectory->Get("ncem_NTotHits");
  ncem_NSvxHits= (TH1F*) gDirectory->Get("ncem_NSvxHits");
  ncem_NSvxRPhiHits= (TH1F*) gDirectory->Get("ncem_NSvxRPhiHits");
  ncem_NSvxSASHits= (TH1F*) gDirectory->Get("ncem_NSvxSASHits");

  ncem_ClusterSize= (TH1F*) gDirectory->Get("ncem_ClusterSize");
  ncem_NTracks= (TH1F*) gDirectory->Get("ncem_NTracks");

  ncem_NCesStrip= (TH1F*) gDirectory->Get("ncem_NCesStrip");
  ncem_ECesStrip= (TH1F*) gDirectory->Get("ncem_ECesStrip");
  ncem_ZCes= (TH1F*) gDirectory->Get("ncem_ZCes");
  ncem_DelZ= (TH1F*) gDirectory->Get("ncem_DelZ");
  ncem_Chi2Strip= (TH1F*) gDirectory->Get("ncem_Chi2Strip");

  ncem_NCesWire= (TH1F*) gDirectory->Get("ncem_NCesWire");
  ncem_ECesWire= (TH1F*) gDirectory->Get("ncem_ECesWire");
  ncem_XCes= (TH1F*) gDirectory->Get("ncem_XCes");
  ncem_DelX= (TH1F*) gDirectory->Get("ncem_DelX");
  ncem_Chi2Wire= (TH1F*) gDirectory->Get("ncem_Chi2Wire");

  ncem_Isol_Et    = (TH2F*) gDirectory->Get("ncem_Isol_Et");
  ncem_HadEm_E    = (TH2F*) gDirectory->Get("ncem_HadEm_E");
  ncem_EOverP_E   = (TH2F*) gDirectory->Get("ncem_EOverP_E");
  ncem_Phi0_EmPhi = (TH2F*) gDirectory->Get("ncem_Phi0_EmPhi");
  ncem_DelX_Et    = (TH2F*) gDirectory->Get("ncem_DelX_Et");
  ncem_DelZ_Et    = (TH2F*) gDirectory->Get("ncem_DelZ_Et");
  ncem_DelX_Wedge = (TH2F*) gDirectory->Get("ncem_DelX_Wedge");
  ncem_DelZ_Wedge = (TH2F*) gDirectory->Get("ncem_DelZ_Wedge");
  ncem_DelEStrip_E = (TH2F*) gDirectory->Get("ncem_DelEStrip_E");
  ncem_DelEWire_E = (TH2F*) gDirectory->Get("ncem_DelEWire_E");

  npem_E= (TH1F*) gDirectory->Get("npem_E");
  npem_Et= (TH1F*) gDirectory->Get("npem_Et");
  npem_Phi= (TH1F*) gDirectory->Get("npem_Phi");
  npem_Eta= (TH1F*) gDirectory->Get("npem_Eta");

  npem_SeedPhi= (TH1F*) gDirectory->Get("npem_SeedPhi");
  npem_SeedEta= (TH1F*) gDirectory->Get("npem_SeedEta");
  npem_Wedge= (TH1F*) gDirectory->Get("npem_Wedge");
  npem_Side= (TH1F*) gDirectory->Get("npem_Side");

  npem_EmPhi= (TH1F*) gDirectory->Get("npem_EmPhi");
  npem_EmDetEta= (TH1F*) gDirectory->Get("npem_EmDetEta");
  npem_EmZv= (TH1F*) gDirectory->Get("npem_EmZv");
  npem_EmEvEta= (TH1F*) gDirectory->Get("npem_EmEvEta");

  npem_EmEt= (TH1F*) gDirectory->Get("npem_EmEt");
  npem_HadEt= (TH1F*) gDirectory->Get("npem_HadEt");
  npem_EmE= (TH1F*) gDirectory->Get("npem_EmE");
  npem_HadE= (TH1F*) gDirectory->Get("npem_HadE");

  npem_HadEm= (TH1F*) gDirectory->Get("npem_HadEm");
  npem_Isol= (TH1F*) gDirectory->Get("npem_Isol");
  npem_Chi2Three= (TH1F*) gDirectory->Get("npem_Chi2Three");

  npem_EOverP= (TH1F*) gDirectory->Get("npem_EOverP");
  npem_Pt= (TH1F*) gDirectory->Get("npem_Pt");
  npem_Charge= (TH1F*) gDirectory->Get("npem_Charge");
  npem_Z0= (TH1F*) gDirectory->Get("npem_Z0");  
  npem_D0= (TH1F*) gDirectory->Get("npem_D0");  
  npem_Phi0= (TH1F*) gDirectory->Get("npem_Phi0");  
  npem_Eta0= (TH1F*) gDirectory->Get("npem_Eta0");  

  npem_Chi2Track= (TH1F*) gDirectory->Get("npem_Chi2Track");
  npem_Chi2Cot= (TH1F*) gDirectory->Get("npem_Chi2Cot");
  npem_Chi2Svx= (TH1F*) gDirectory->Get("npem_Chi2Svx");
  npem_nDof= (TH1F*) gDirectory->Get("npem_nDof");
  npem_Chi2perDof= (TH1F*) gDirectory->Get("npem_Chi2perDof");

  npem_NAxHits= (TH1F*) gDirectory->Get("npem_NAxHits");
  npem_NStHits= (TH1F*) gDirectory->Get("npem_NStHits");
  npem_NTotHits= (TH1F*) gDirectory->Get("npem_NTotHits");
  npem_NSvxHits= (TH1F*) gDirectory->Get("npem_NSvxHits");
  npem_NSvxRPhiHits= (TH1F*) gDirectory->Get("npem_NSvxRPhiHits");
  npem_NSvxSASHits= (TH1F*) gDirectory->Get("npem_NSvxSASHits");

  npem_ClusterSize= (TH1F*) gDirectory->Get("npem_ClusterSize");
  npem_NTracks= (TH1F*) gDirectory->Get("npem_NTracks");

  npem_NPesClusters= (TH1F*) gDirectory->Get("npem_NPesClusters");
  npem_PesEta= (TH1F*) gDirectory->Get("npem_PesEta");
  npem_PesPhi= (TH1F*) gDirectory->Get("npem_PesPhi");
  npem_EPesU= (TH1F*) gDirectory->Get("npem_EPesU");
  npem_EPesV= (TH1F*) gDirectory->Get("npem_EPesV");
  npem_DelR= (TH1F*) gDirectory->Get("npem_DelR");

  npem_Isol_Et    = (TH2F*) gDirectory->Get("npem_Isol_Et");
  npem_HadEm_E    = (TH2F*) gDirectory->Get("npem_HadEm_E");
  npem_EOverP_E   = (TH2F*) gDirectory->Get("npem_EOverP_E");
  npem_Chi2Three_Et   = (TH2F*) gDirectory->Get("npem_Chi2Three_Et");
  npem_DelR_Et   = (TH2F*) gDirectory->Get("npem_DelR_Et");
  npem_DelR_Wedge   = (TH2F*) gDirectory->Get("npem_DelR_Wedge");
  npem_Phi0_EmPhi = (TH2F*) gDirectory->Get("npem_Phi0_EmPhi");
  npem_PesPhi_EmPhi = (TH2F*) gDirectory->Get("npem_PesPhi_EmPhi");
  npem_PesEta_EmEta = (TH2F*) gDirectory->Get("npem_PesEta_EmEta");
  npem_DelEU_E = (TH2F*) gDirectory->Get("npem_DelEU_E");
  npem_DelEV_E = (TH2F*) gDirectory->Get("npem_DelEV_E");

  return true;		       
}


//==================================
//       FillHistograms()
//==================================
void TEleCand::FillHistograms(TStnElectron* electron) {
  
  // Make these plots for all electron candidates, with no TEleCand specified cuts
  DetCode->Fill(electron->DetectorCode());

  if(electron->IsCentral()){
    //Best estimate of 4 momentum
    ncem_E->Fill(electron->Momentum()->E());
    ncem_Et->Fill(electron->Momentum()->Perp());
    ncem_Phi->Fill(cdfphi(electron->Momentum()->Phi()));
    ncem_Eta->Fill(electron->Momentum()->Eta());

    //Seed tower location of em cluster
    ncem_SeedPhi->Fill(electron->SeedIPhi());
    ncem_SeedEta->Fill(electron->SeedIEta());
    ncem_Wedge->Fill(electron->WedgeNumber());
    ncem_Side->Fill(electron->Side());

    //Em cluster angles
    ncem_EmPhi->Fill(electron->EmClusPhi());
    ncem_EmDetEta->Fill(electron->DetEta());
    ncem_EmZv->Fill(electron->EmClusZv());
    ncem_EmEvEta->Fill(electron->EmClusEvEta());

    //Em cluster energy
    ncem_EmEt->Fill(electron->EmEt());
    ncem_HadEt->Fill(electron->HadEt());
    ncem_EmE->Fill(electron->EmE());
    ncem_HadE->Fill(electron->HadE());
    ncem_HadEm->Fill(electron->HadEm());

    //Topology and id variables
    ncem_Isol->Fill(electron->Iso()/electron->Momentum()->Perp());
    ncem_Lshr->Fill(electron->Lshr());
    

    //Track quantities
    int eleTrack = electron->TrackNumber();  
    if(eleTrack != -1 && electron->NTracks() !=0 && eleTrack < pTrackBlock->NTracks()){
      TStnTrack* track = pTrackBlock->Track(eleTrack);

      ncem_EOverP->Fill(electron->EOverP());
      ncem_Pt->Fill(electron->TrackPt());
      ncem_Charge->Fill(electron->Charge());

      ncem_Z0->Fill(track->Z0());
      ncem_D0->Fill(track->D0());
      ncem_Phi0->Fill(cdfphi(track->Phi0()));
      ncem_Eta0->Fill(track->Eta());


      ncem_NAxHits->Fill(track->NCotHitsAx());
      ncem_NStHits->Fill(track->NCotHitsSt());
      ncem_NTotHits->Fill(track->NCotHitsTot());
      ncem_NSvxHits->Fill(track->NSvxHits());
      ncem_NSvxRPhiHits->Fill(track->NSvxRPhiHits());
      ncem_NSvxSASHits->Fill(track->NSvxSASHits());

      //Chi2
      ncem_Chi2Track->Fill(track->Chi2());
      if(track->NCotHitsTot() > 0) ncem_Chi2Cot->Fill(track->Chi2Cot());
      if(track->NSvxHits() > 0) ncem_Chi2Svx->Fill(track->Chi2Svx());

      //way to calculate number of degrees of freedom???
      int nDof = track->NCotHitsTot() + track->NSvxHits() - 5;
      ncem_nDof->Fill(nDof);
      if(nDof > 0) ncem_Chi2perDof->Fill(track->Chi2()/nDof);     	

      ncem_DelX->Fill(electron->DelX());
      ncem_DelZ->Fill(electron->DelZ());

      ncem_EOverP_E->Fill(electron->Momentum()->E(),electron->EOverP() );
      ncem_Phi0_EmPhi->Fill(electron->EmClusPhi(),track->Phi0() );

    }
    
    //CES
    ncem_ClusterSize->Fill(electron->ClusterSize());
    ncem_NTracks->Fill(electron->NTracks());

    ncem_NCesStrip->Fill(electron->NCesClusters(0));
    ncem_ECesStrip->Fill(electron->CesEnergy(0));
    ncem_ZCes->Fill(electron->ZCes());
    ncem_Chi2Strip->Fill(electron->Chi2Strip());      

    ncem_NCesWire->Fill(electron->NCesClusters(1));
    ncem_ECesWire->Fill(electron->CesEnergy(1));
    ncem_XCes->Fill(electron->XCes());
    ncem_Chi2Wire->Fill(electron->Chi2Wire());      

    //2D plots
    ncem_Isol_Et->Fill(electron->Et(),electron->Iso()/electron->Et() );
    ncem_HadEm_E->Fill(electron->Momentum()->E(),electron->HadEm() );
    ncem_DelX_Et->Fill(electron->Et(),electron->DelX() );
    ncem_DelZ_Et->Fill(electron->Et(),electron->DelZ() );
    ncem_DelX_Wedge->Fill(electron->WedgeNumber(),electron->DelX() );
    ncem_DelZ_Wedge->Fill(electron->WedgeNumber(),electron->DelZ() );
    ncem_DelEStrip_E->Fill(electron->Momentum()->E(), electron->CesEnergy(0)-electron->Momentum()->E() );
    ncem_DelEWire_E->Fill(electron->Momentum()->E(), electron->CesEnergy(1)-electron->Momentum()->E() );


  }
  if(electron->IsPlug()){
    
    npem_E->Fill(electron->Momentum()->E());
    npem_Et->Fill(electron->Momentum()->Perp());
    npem_Phi->Fill(cdfphi(electron->Momentum()->Phi()));
    npem_Eta->Fill(electron->Momentum()->Eta());

    //Seed tower location of em cluster
    npem_SeedPhi->Fill(electron->SeedIPhi());
    npem_SeedEta->Fill(electron->SeedIEta());
    npem_Wedge->Fill((electron->SeedIPhi() + 48*electron->Side()));
    npem_Side->Fill(electron->Side());

    //Em cluster angles
    npem_EmPhi->Fill(electron->EmClusPhi());
    npem_EmDetEta->Fill(electron->DetEta());
    npem_EmZv->Fill(electron->EmClusZv());
    npem_EmEvEta->Fill(electron->EmClusEvEta());

    //Em cluster energy
    npem_EmEt->Fill(electron->EmEt());
    npem_HadEt->Fill(electron->HadEt());
    npem_EmE->Fill(electron->EmE());
    npem_HadE->Fill(electron->HadE());
    npem_HadEm->Fill(electron->HadEm());

    //Topology and id variables
    npem_Isol->Fill(electron->Iso()/electron->Momentum()->Perp());
    

    //Track quantities
    int eleTrack = electron->TrackNumber();  
    if(eleTrack != -1 && electron->NTracks() !=0 && eleTrack < pTrackBlock->NTracks()){
      TStnTrack* track = pTrackBlock->Track(eleTrack);

      npem_EOverP->Fill(electron->EOverP());
      npem_Pt->Fill(electron->TrackPt());
      npem_Charge->Fill(electron->Charge());

      npem_Z0->Fill(track->Z0());
      npem_D0->Fill(track->D0());
      npem_Phi0->Fill(track->Phi0());
      npem_Eta0->Fill(track->Eta());

      npem_NAxHits->Fill(track->NCotHitsAx());
      npem_NStHits->Fill(track->NCotHitsSt());
      npem_NTotHits->Fill(track->NCotHitsTot());
      npem_NSvxHits->Fill(track->NSvxHits());
      npem_NSvxRPhiHits->Fill(track->NSvxRPhiHits());
      npem_NSvxSASHits->Fill(track->NSvxSASHits());

      //Chi2
      npem_Chi2Track->Fill(track->Chi2());
      if(track->NCotHitsTot() > 0) npem_Chi2Cot->Fill(track->Chi2Cot());
      if(track->NSvxHits() > 0) npem_Chi2Svx->Fill(track->Chi2Svx());

      //way to calculate number of degrees of freedom???
      int nDof = track->NCotHitsTot() + track->NSvxHits() - 5;
      npem_nDof->Fill(nDof);
      if(nDof > 0) npem_Chi2perDof->Fill(track->Chi2()/nDof);     	

      npem_EOverP_E->Fill(electron->Momentum()->E(),electron->EOverP() );
      npem_Phi0_EmPhi->Fill(electron->EmClusPhi(),track->Phi0() );

    }
    
    npem_ClusterSize->Fill(electron->ClusterSize());
    npem_NTracks->Fill(electron->NTracks());

    npem_Chi2Three->Fill(electron->Chi2Three());
    npem_NPesClusters->Fill(electron->NPesClusters());
    if(electron->NPesClusters() >0){
      npem_PesEta->Fill(electron->PesEta());
      npem_PesPhi->Fill(electron->PesPhi());
      npem_EPesU->Fill(electron->PesEnergy(0));
      npem_EPesV->Fill(electron->PesEnergy(1));
      
      float deltaR = sqrt( (electron->PesEta() - electron->DetEta()) * 
			   (electron->PesEta() - electron->DetEta()) +
			   (electron->PesPhi() - electron->EmClusPhi())  * 
			   (electron->PesPhi() - electron->EmClusPhi())  );      
      npem_DelR->Fill(deltaR);
      
      npem_DelR_Et->Fill(electron->Et(),deltaR );
      npem_DelR_Wedge->Fill(electron->SeedIPhi() + 48*electron->Side(),deltaR );
      npem_PesPhi_EmPhi->Fill(electron->EmClusPhi(),electron->PesPhi() );
      npem_PesEta_EmEta->Fill(electron->DetEta(),electron->PesEta() );
      npem_DelEU_E->Fill(electron->Momentum()->E(), electron->PesEnergy(0)-electron->Momentum()->E() );
      npem_DelEV_E->Fill(electron->Momentum()->E(), electron->PesEnergy(1)-electron->Momentum()->E() );
    }

    //2D plots
    npem_Isol_Et->Fill(electron->Et(),electron->Iso()/electron->Et() );
    npem_HadEm_E->Fill(electron->Momentum()->E(),electron->HadEm() );
    npem_Chi2Three_Et->Fill(electron->Et(),electron->Chi2Three() );


  }
    // Was this electron found as a conversion?
    bool eleConv = false; //eleConvCand->IsaConvElectron(i);

  return;		       
}

//------------------------------------------------------------------------
// This function returns the phi value in 0 to 2*PI range.
//------------------------------------------------------------------------
float TEleCand::cdfphi(const float &phi) {
	Float_t pi = 3.141592653;
	Float_t twopi = 2*pi;
	Float_t myphi = phi;
	if (myphi < 0) {
		while (myphi < 0) {
			myphi = myphi + twopi;
		}
	} else {
		while (myphi > twopi) {
			myphi = myphi - twopi;
		}
	}
	return myphi;
}
//==================================
//       EndJob()
//==================================
int TEleCand::EndJob() {

  printf("----- end job: ---- %s\n",GetName());
  printf("  Events processed %d with good electrons %d \n",
	 totEvents,eleEvents);

  Efficiency(isEnableTopEvent());

  return 0;
}

//==========================================================
bool TEleCand::isZBoson(TStnElectron* tight, TStnElectron* ele2) {

  bool isZ = false;
  
  
  // Reconstruct the Z boson...take the first lep. and add the second.
  const TLorentzVector* p1 = (TLorentzVector*)tight->Momentum();
  const TLorentzVector* p2 = (TLorentzVector*)ele2->Momentum();
  TLorentzVector *pair = new TLorentzVector();
  pair->operator+=(*p1);
  pair->operator+=(*p2);
  pairMass->Fill(pair->M());
  
  if (getDebugLevel() > 4) cout << endl << " Tight electron E     " << tight->Momentum()->E() 
			    << " Px " << tight->Momentum()->Px() 
			    << " Py " << tight->Momentum()->Py() 
			    << " Pz " << tight->Momentum()->Pz()
			    << " M  " << tight->Momentum()->M() << endl;
  
  if (getDebugLevel() > 4) cout << " Candidate electron E " << ele2->Momentum()->E() 
			    << " Px " << ele2->Momentum()->Px() 
			    << " Py " << ele2->Momentum()->Py() 
			    << " Pz " << ele2->Momentum()->Pz()
			    << " M  " << ele2->Momentum()->M() << endl;

  if (getDebugLevel() > 4) cout << " Pair                 " << pair->E() 
			    << " Px " << pair->Px() 
			    << " Py " << pair->Py() 
			    << " Pz " << pair->Pz() 
			    << " M  " << pair->M() << endl;

  //Is the pair mass within a window around the Z mass?
  double upperZ = 110.0;
  double lowerZ = 70.0;
  if(pair->M() >= lowerZ && pair->M() <= upperZ){
    isZ = true;

    //Fill histograms for the second leg
    FillHistograms(ele2);

    //
    //Count electrons that pass various cuts
    nZees++;
    if(ele2->DetectorCode() == 0){
      nCem++;
      if(ele2->Et()         > getEleEtCut() ) nPassCut[1]++;
      if(ele2->TrackPt()    > getElePtCut() ) nPassCut[2]++;
      if(ele2->EOverP()     < getEleEOverPCut() ) nPassCut[3]++;
      if(fabs(ele2->DelX()) < getEleDelXCut() )nPassCut[4]++;
      if(fabs(ele2->DelZ()) < getEleDelZCut() )nPassCut[5]++;
      if(ele2->Lshr()       < getEleLshrCut() )nPassCut[6]++;
      if(ele2->HadEm()      < getEleHadEmCut() )nPassCut[7]++;
      if((ele2->Iso()/ele2->Et()) < getEleIsolCut() ) nPassCut[8]++;
      if(ele2->Chi2Strip()  < getEleChi2StripCut() )nPassCut[9]++;
      

      if( ( (!isEnableEtCut())       || (ele2->Et()         > getEleEtCut()) ) &&
	  ((!isEnablePtCut())       || (ele2->TrackPt()> getElePtCut()) ) &&
	  ((!isEnableEOverPCut())   || (ele2->EOverP()     < getEleEOverPCut()) ) &&
	  ((!isEnableDelXCut())     || (fabs(ele2->DelX()) < getEleDelXCut()) ) &&
	  ((!isEnableDelZCut())     || (fabs(ele2->DelZ()) < getEleDelZCut()) ) &&
	  ((!isEnableLshrCut())     || (ele2->Lshr()       < getEleLshrCut()) ) &&
	  ((!isEnableHadEmCut())    || (ele2->HadEm()      < getEleHadEmCut()) ) &&
	  ((!isEnableIsolCut())     || ((ele2->Iso()/ele2->Et()) < getEleIsolCut()) ) &&
	  ((!isEnableTisolCut())    || (ele2->Tiso()       < getEleTisolCut()) ) &&
	  ((!isEnableChi2StripCut())|| (ele2->Chi2Strip()  < getEleChi2StripCut()) )
	  ){
	nPassCut[0]++;
      }
    }
  }
  return isZ;
}
//==========================================================
void TEleCand::Efficiency(bool ENABLE) {

  if(ENABLE){
    cout << "Number of Z candidates is "<< nZees << " and in CEM is " << nCem << endl;
    for(int i=0;i<20;i++){
      cout << "Number passing cuts " << nameCut[i] << " = " << nPassCut[i] << endl;
    }
  }


}
//==================================
//       myFindTrackZ0()
//==================================
double TEleCand::myFindTrackZ0(TStnElectron* electron) {
 
  double elez0 = -999.0;
  //
  //Use the Stntuple track link to find the electron track
  int eleTrack = electron->TrackNumber();  
  if(eleTrack != -1 && electron->NTracks() !=0 && eleTrack < pTrackBlock->NTracks()){
    TStnTrack* track = pTrackBlock->Track(eleTrack);
    elez0 = track->Z0();

    if(getDebugLevel() > 1) printf(" eleTrack pt %f phi %f eta %f z0 %f d0 %f index %d NTracks % d \n",
				   track->Pt(),track->Phi0(),track->Eta(),track->Z0(), track->D0(),
				   eleTrack, electron->NTracks());  
  }
  return elez0;
}

//==================================
//       myFindClosestJet()
//==================================
int TEleCand::myFindClosestJet(TStnElectron* electron) {

  int jetIndex = -1;  
  // Find the jet (from original, not reclustered, collection) that is closest to the electron
  float minDeltaR = 1000.0;
  for (int j=0; j<pJetBlock->NJets(); j++) {
    TStnJet* jet = pJetBlock->Jet(j);
    float deltaR = sqrt( (jet->Momentum()->Eta() - electron->Momentum()->Eta()) * 
			 (jet->Momentum()->Eta() - electron->Momentum()->Eta()) +
			 (jet->Momentum()->Phi() - electron->Momentum()->Phi()) *
			 (jet->Momentum()->Phi() - electron->Momentum()->Phi()) );
    if (deltaR < minDeltaR) {
      minDeltaR = deltaR;
      jetIndex = j;
    }
  }

  if (jetIndex != -1){
    TStnJet* elejet = pJetBlock->Jet(jetIndex);
    //
    // Some printout
    if (getDebugLevel() > 0) {
      printf("\n Electron Et %f phi %f eta %f cluster pt %f track pt %f \n closejet Et %f phi %f eta %f emf %f DelR %f\n",
	     electron->Et(),cdfphi(electron->Momentum()->Phi()),electron->Momentum()->Eta(),
	     electron->Momentum()->Pt(),electron->TrackPt(),
	     elejet->Et(),cdfphi(elejet->Momentum()->Phi()),elejet->Momentum()->Eta(),elejet->EmFract(),minDeltaR);
    }
  }
  
  return jetIndex;
}

//old code, before TopEventModule!
//Make TEleCand cuts on Stntuple electrons, flag good electrons
//==================================
//       oldCodeElectron()
//==================================
void TEleCand::oldCodeElectron() {
  
  // Loop over the standard electrons
  bool foundEle = false;
  NgoodElectron = 0;
  
  for (int i=0; i<pEleBlock->NElectrons(); i++) {
    TStnElectron* electron = pEleBlock->Electron(i);
    
    // Was this electron found as a conversion?
    bool eleConv = false; //eleConvCand->IsaConvElectron(i);
    //
    //What is the closest jet to this electron?
    int jetIndex = myFindClosestJet(electron);
    ElectronJetIndex[i] = jetIndex;
    //
    // Make the cuts.  Note that if a cut is not enabled, then NO REQUIREMENT
    // on that variable is made (it is as if the variable is not cut on).
    if((electron->DetectorCode() == 0) &&
       ((!isEnableEtCut())       || (electron->Et()         > getEleEtCut()) ) &&
       ((!isEnablePtCut())       || (electron->TrackPt()> getElePtCut()) ) &&
       ((!isEnableRemoveConv())  || (!eleConv)) &&
       ((!isEnableEOverPCut())   || (electron->EOverP()     < getEleEOverPCut()) ) &&
       ((!isEnableDelXCut())     || (fabs(electron->DelX()) < getEleDelXCut()) ) &&
       ((!isEnableDelZCut())     || (fabs(electron->DelZ()) < getEleDelZCut()) ) &&
       ((!isEnableLshrCut())     || (electron->Lshr()       < getEleLshrCut()) ) &&
       ((!isEnableHadEmCut())    || (electron->HadEm()      < getEleHadEmCut()) ) &&
       ((!isEnableIsolCut())     || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ((!isEnableTisolCut())    || (electron->Tiso()       < getEleTisolCut()) ) &&
       ((!isEnableChi2StripCut())|| (electron->Chi2Strip()  < getEleChi2StripCut()) )
       ){
      NgoodElectronList[i] = true;
      foundEle = true;
      if (NgoodElectron < NElectronMax) {
	goodElectron[NgoodElectron] = electron;
	NgoodElectron++;
      }else{
	printf("Too many good Electrons (%d)! NElectronMax=%d must be increased! \n",NgoodElectron,NElectronMax);
      }
    } else {
      NgoodElectronList[i] = false;
    }
  }
  //
  // If at least one good electron was found....
  if (foundEle) {
    eleEvents++;
    if (getDebugLevel() > 0) {
      printf(" Good electron run %d event %d\n",pHeaderBlock->RunNumber(),pHeaderBlock->EventNumber());
    }
  }

  
  return;
}
/*
//Temporary holding place for this code 
  //TEleCand specified cut plots
  // Make plots of each cut variable, requiring all of the OTHER, ENABLED cuts
  if(  (electron->DetectorCode() == 0) &&
       ( true                      || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_Et->Fill(electron->Et());
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( true                      || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_Pt->Fill(electron->TrackPt());
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( true                      || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_EOverP->Fill(electron->EOverP());
  
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( true                      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_DelX->Fill(electron->DelX());
  
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( true                      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_DelZ->Fill(electron->DelZ());
  
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( true                      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_Lshr->Fill(electron->Lshr());
  
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( true                      || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_HadEm->Fill(electron->HadEm());
  
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( true                      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( (!isEnableChi2StripCut()) || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_Isol->Fill(electron->Iso()/electron->Et());
  
  
  if(  (electron->DetectorCode() == 0) &&
       ( (!isEnableEtCut())        || (electron->Et()     > getEleEtCut()) ) &&
       ( (!isEnablePtCut())        || (electron->TrackPt()> getElePtCut()) ) &&
       ( (!isEnableEOverPCut())    || (electron->EOverP() < getEleEOverPCut()) ) &&
       ( (!isEnableDelXCut())      || (electron->DelX()   < getEleDelXCut()) ) &&
       ( (!isEnableDelZCut())      || (electron->DelZ()   < getEleDelZCut()) ) &&
       ( (!isEnableLshrCut())      || (electron->Lshr()   < getEleLshrCut()) ) &&
       ( (!isEnableHadEmCut())     || (electron->HadEm()  < getEleHadEmCut()) ) &&
       ( (!isEnableIsolCut())      || ((electron->Iso()/electron->Et()) < getEleIsolCut()) ) &&
       ( true                      || (electron->Chi2Strip()            < getEleChi2StripCut()) ) &&
       ( (!isEnableRemoveConv())   || (!eleConv)  ) 
       ) cem_Chi2Strip->Fill(electron->Chi2Strip());
//Temporary holding place for this code */