//--------------------------------------------------------------------------
// Description:
// ------------
//
//	Class SimpleNtuple : ROOT ntuples made easy(?)
//
//--------------------------------------------------------------------------

//-----------------------
// This Class's Header --
//-----------------------
#include "HighLevelObjects/SimpleNtuple.hh"

//---------------
// C++ Headers --
//---------------
#include <assert.h>
#include <iostream>
#include <vector>

//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "BaBar/Cdf.hh"
#include "Edm/EventRecord.hh"
#include "Edm/ConstHandle.hh"
#include "Edm/Handle.hh"
#include "TrackingObjects/Storable/CdfTrackColl.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Tracks/track_sort.hh"
#include "ElectronObjects/CdfEmObjectColl.hh"
#include "HighLevelObjects/ElectronVariables.hh"
#include "AbsEnv/AbsEnv.hh"
#include "ElectronObjects/CdfEmObject.hh"
#include "ElectronObjects/CdfEmObjectColl.hh"
#include "CalorObjects/Pes2dCluster.hh"
#include "CalorObjects/Pes2dClusterColl.hh"
#include "CalorObjects/PesCluster.hh"
#include "CalorObjects/PesClusterColl.hh"
#include "CalorObjects/PlugStrip.hh"
#include "CalorObjects/PlugStripColl.hh"

#ifdef __GNUG__
#pragma implementation
#endif

#include <string>
#include <cstdio>
#include <assert.h>
#include <iostream>
#include <iomanip>
#include "cern_i/packlib.h"
#include "Experiment/Experiment.hh"
#include "AbsEnv/AbsEnv.hh"
#include "Calor/Cal2.hh"
#include "Bbos/bank_interface.hh"
#include "TH1.h"
#include "TH2.h"
#include "TVector2.h"
#include "MetObjects/CdfMet.hh"
#include "Met/CdfMetPuffer.hh"
#include "JetObjects/CdfJetColl.hh"
#include "JetObjects/JetAlgParams.hh"
#include "Jet/JetUtils.hh"
#include "SimulationObjects/HEPG_StorableBank.hh"
#include "VertexObjects/VertexColl.hh"
#include "VertexObjects/Vertex.hh"

//----------------
// Constructors --/////////////////////////////////////////
//----------------

SimpleNtuple::SimpleNtuple(
			   const char* const theName,
			   const char* const theDescription )
  : AppModule( theName, theDescription ),
    SimpleNtuplefile( "SimpleNtuplefile", this, "SimpleNtuple.root")
{
  std::cout << " " << std::endl;
  std::cout << " " << std::endl;
  std::cout << "   Hello from module SimpleNtuple.   " << std::endl;
  std::cout << " " << std::endl;
  commands()->append( &SimpleNtuplefile );
}

//////////////////////////////////////////////////////////

//--------------
// Destructor --\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//--------------

SimpleNtuple::~SimpleNtuple( )
{
  std::cout << " " << std::endl;
  std::cout << " " << std::endl;
  std::cout << "   Goodbye from module SimpleNtuple.   " << std::endl;
  std::cout << " " << std::endl;
}

//\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

AppResult SimpleNtuple::beginRun( AbsEvent *aJob) {
  
  std::cout << "Begin Run SimpleNtuple" << std::endl;
  
  return AppResult::OK;
}

AppResult SimpleNtuple::beginJob( AbsEvent *aJob) {
  
  std::cout << "Begin Job SimpleNtuple" << std::endl;
  
  hfile = new TFile(SimpleNtuplefile.value().c_str(),"RECREATE","Simple ROOT Ntuple");
  
  std::cout << "Opening File " << SimpleNtuplefile.value().c_str() << std::endl ;
  
  std::cout << "Creating Branches" << std::endl ;
  
  // comment out creation of "ntuple".  this is more limited than using 
  // trees, so opt for creation of TTree with multiple branches instead...
  
  // ntuple = new TNtuple("ntuple","Electron Quantities","Run:Event:HadOverEm");
  
  SimpleNtupleTree = new TTree("SimpleNtupleTree","SimpleNtupleTree");

  SimpleNtupleTree->Branch("numEventObjs",&numEventObjs,"numEventObjs/I");
  SimpleNtupleTree->Branch("numHepgObjs",&numHepgObjs,"numHepgObjs/I");
  SimpleNtupleTree->Branch("numVertexObjs",&numVertexObjs,"numVertexObjs/I");
  SimpleNtupleTree->Branch("numEmObjs",&numEmObjs,"numEmObjs/I");
  SimpleNtupleTree->Branch("numMetObjs",&numMetObjs,"numMetObjs/I");
  SimpleNtupleTree->Branch("numPes2dObjs",&numPes2dObjs,"numPes2dObjs/I");
  SimpleNtupleTree->Branch("numPesObjs",&numPesObjs,"numPesObjs/I");
  SimpleNtupleTree->Branch("numJetObjs",&numJetObjs,"numJetObjs/I");  

  SimpleNtupleTree->Branch("runNumber",&runNumber,"runNumber[numEventObjs]/I");
  SimpleNtupleTree->Branch("eventNumber",&eventNumber,"eventNumber[numEventObjs]/I");

  SimpleNtupleTree->Branch("hepgElectronEta",&hepgElectronEta,"hepgElectronEta[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronPhi",&hepgElectronPhi,"hepgElectronPhi[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronEt",&hepgElectronEt,"hepgElectronEt[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronPx",&hepgElectronPx,"hepgElectronPx[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronPy",&hepgElectronPy,"hepgElectronPy[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronPz",&hepgElectronPz,"hepgElectronPz[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronE",&hepgElectronE,"hepgElectronE[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronVx",&hepgElectronVx,"hepgElectronVx[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronVy",&hepgElectronVy,"hepgElectronVy[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronVz",&hepgElectronVz,"hepgElectronVz[numHepgObjs]/D");
  SimpleNtupleTree->Branch("hepgElectronVt",&hepgElectronVt,"hepgElectronVt[numHepgObjs]/D");

  SimpleNtupleTree->Branch("vertexPrimX",&vertexPrimX,"vertexPrimX[numVertexObjs]/D");
  SimpleNtupleTree->Branch("vertexPrimY",&vertexPrimY,"vertexPrimY[numVertexObjs]/D");
  SimpleNtupleTree->Branch("vertexPrimZ",&vertexPrimZ,"vertexPrimZ[numVertexObjs]/D");

  SimpleNtupleTree->Branch("seedEnergy",&seedEnergy,"seedEnergy[numEmObjs]/D");
  SimpleNtupleTree->Branch("seedEt",&seedEt,"seedEt[numEmObjs]/D");
  SimpleNtupleTree->Branch("region",&region,"region[numEmObjs]/I");
  SimpleNtupleTree->Branch("side",&side,"side[numEmObjs]/I");
  SimpleNtupleTree->Branch("etaDetector",&etaDetector,"etaDetector[numEmObjs]/D");
  SimpleNtupleTree->Branch("etaEvent",&etaEvent,"etaEvent[numEmObjs]/D");
  SimpleNtupleTree->Branch("etaRMS",&etaRMS,"etaRMS[numEmObjs]/D");
  SimpleNtupleTree->Branch("phi",&phi,"phi[numEmObjs]/D");
  SimpleNtupleTree->Branch("phiRMS",&phiRMS,"phiRMS[numEmObjs]/D");
  SimpleNtupleTree->Branch("emCalEx",&emCalEx,"emCalEx[numEmObjs]/D");
  SimpleNtupleTree->Branch("emCalEy",&emCalEy,"emCalEy[numEmObjs]/D");
  SimpleNtupleTree->Branch("emCalEz",&emCalEz,"emCalEz[numEmObjs]/D");
  SimpleNtupleTree->Branch("pemPesMatch",&pemPesMatch,"pemPesMatch[numEmObjs]/D");
  SimpleNtupleTree->Branch("pemPesDeta",&pemPesDeta,"pemPesDeta[numEmObjs]/D");
  SimpleNtupleTree->Branch("pemPesDphi",&pemPesDphi,"pemPesDphi[numEmObjs]/D");
  SimpleNtupleTree->Branch("zVertex",&zVertex,"zVertex[numEmObjs]/D");
  SimpleNtupleTree->Branch("zCentroid",&zCentroid,"zCentroid[numEmObjs]/D");
  SimpleNtupleTree->Branch("emEnergy",&emEnergy,"emEnergy[numEmObjs]/D");
  SimpleNtupleTree->Branch("emEt",&emEt,"emEt[numEmObjs]/D");
  //  SimpleNtupleTree->Branch("correctedEt",&correctedEt,"correctedEt[numEmObjs]/D");
  SimpleNtupleTree->Branch("hadEnergy",&hadEnergy,"hadEnergy[numEmObjs]/D");
  SimpleNtupleTree->Branch("hadEt",&hadEt,"hadEt[numEmObjs]/D");
  SimpleNtupleTree->Branch("totalEnergy",&totalEnergy,"totalEnergy[numEmObjs]/D");
  SimpleNtupleTree->Branch("totalEt",&totalEt,"totalEt[numEmObjs]/D");
  SimpleNtupleTree->Branch("hadEm",&hadEm,"hadEm[numEmObjs]/D");
  SimpleNtupleTree->Branch("totalIso4",&totalIso4,"totalIso4[numEmObjs]/D");
  SimpleNtupleTree->Branch("totalIso7",&totalIso7,"totalIso7[numEmObjs]/D");
  SimpleNtupleTree->Branch("emIso4",&emIso4,"emIso4[numEmObjs]/D");
  SimpleNtupleTree->Branch("emIso7",&emIso7,"emIso7[numEmObjs]/D");
  SimpleNtupleTree->Branch("hadIso4",&hadIso4,"hadIso4[numEmObjs]/D");
  SimpleNtupleTree->Branch("hadIso7",&hadIso7,"hadIso7[numEmObjs]/D");
  SimpleNtupleTree->Branch("lshr2",&lshr2,"lshr2[numEmObjs]/D");
  SimpleNtupleTree->Branch("lshr",&lshr,"lshr[numEmObjs]/D");
  SimpleNtupleTree->Branch("trackPt",&trackPt,"trackPt[numEmObjs]/D");
  SimpleNtupleTree->Branch("charge",&charge,"charge[numEmObjs]/D");
  SimpleNtupleTree->Branch("eOverP",&eOverP,"eOverP[numEmObjs]/D");
  SimpleNtupleTree->Branch("xCes",&xCes,"xCes[numEmObjs]/D");
  SimpleNtupleTree->Branch("chiW",&chiW,"chiW[numEmObjs]/D");
  SimpleNtupleTree->Branch("deltaX",&deltaX,"deltaX[numEmObjs]/D");
  SimpleNtupleTree->Branch("zCes",&zCes,"zCes[numEmObjs]/D");
  SimpleNtupleTree->Branch("chiS",&chiS,"chiS[numEmObjs]/D");
  SimpleNtupleTree->Branch("deltaZ",&deltaZ,"deltaZ[numEmObjs]/D");
  SimpleNtupleTree->Branch("numPes2dClusters",&numPes2dClusters,"numPes2dClusters[numEmObjs]/I");
  SimpleNtupleTree->Branch("pesEmObjEta",&pesEmObjEta,"pesEmObjEta[numEmObjs]/D");
  SimpleNtupleTree->Branch("pesEmObjPhi",&pesEmObjPhi,"pesEmObjPhi[numEmObjs]/D");
  SimpleNtupleTree->Branch("pesTestEmObjMatching",&pesTestEmObjMatching,"pesTestEmObjMatching[numEmObjs]/D");
  SimpleNtupleTree->Branch("pesEmObjWidth",&pesEmObjWidth,"pesEmObjWidth[numEmObjs]/D");
  SimpleNtupleTree->Branch("trackIso",&trackIso,"trackIso[numEmObjs]/D");
  SimpleNtupleTree->Branch("nTowers",&nTowers,"nTowers[numEmObjs]/D");
  SimpleNtupleTree->Branch("nTracks",&nTracks,"nTracks[numEmObjs]/D");
  SimpleNtupleTree->Branch("pem3x3FitTowers",&pem3x3FitTowers,"pem3x3FitTowers[numEmObjs]/I");
  SimpleNtupleTree->Branch("pem3x3FitDetEta",&pem3x3FitDetEta,"pem3x3FitDetEta[numEmObjs]/D");
  SimpleNtupleTree->Branch("pem3x3FitPhi",&pem3x3FitPhi,"pem3x3FitPhi[numEmObjs]/D");
  SimpleNtupleTree->Branch("pem3x3EmEnergy",&pem3x3EmEnergy,"pem3x3EmEnergy[numEmObjs]/D");
  SimpleNtupleTree->Branch("pem3x3Chisq",&pem3x3Chisq,"pem3x3Chisq[numEmObjs]/D");

  SimpleNtupleTree->Branch("met",&met,"met[numMetObjs]/D");
  SimpleNtupleTree->Branch("metPhi",&metPhi,"metPhi[numMetObjs]/D");
  SimpleNtupleTree->Branch("metZV",&metZV,"metZV[numMetObjs]/D");
  SimpleNtupleTree->Branch("metESum",&metESum,"metESum[numMetObjs]/D");
  SimpleNtupleTree->Branch("metEtSum",&metEtSum,"metEtSum[numMetObjs]/D");
  SimpleNtupleTree->Branch("metExSum",&metExSum,"metExSum[numMetObjs]/D");
  SimpleNtupleTree->Branch("metEySum",&metEySum,"metEySum[numMetObjs]/D");

  SimpleNtupleTree->Branch("pes2dWidth",&pes2dWidth,"pes2dWidth[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dEnergy",&pes2dEnergy,"pes2dEnergy[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dX",&pes2dX,"pes2dX[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dY",&pes2dY,"pes2dY[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dZ",&pes2dZ,"pes2dZ[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dGlblEta",&pes2dGlblEta,"pes2dGlblEta[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dGlblPhi",&pes2dGlblPhi,"pes2dGlblPhi[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dUEnergy",&pes2dUEnergy,"pes2dUEnergy[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dVEnergy",&pes2dVEnergy,"pes2dVEnergy[numPes2dObjs]/D");
  SimpleNtupleTree->Branch("pes2dNumUStrips",&pes2dNumUStrips,"pes2dNumUStrips[numPes2dObjs]/I");
  SimpleNtupleTree->Branch("pes2dNumVStrips",&pes2dNumVStrips,"pes2dNumVStrips[numPes2dObjs]/I");

  SimpleNtupleTree->Branch("pesWidthStrips",&pesWidthStrips,"pesWidthStrips[numPesObjs]/D");
  SimpleNtupleTree->Branch("pesEnergy",&pesEnergy,"pesEnergy[numPesObjs]/D");
  SimpleNtupleTree->Branch("pesSide",&pesSide,"pesSide[numPesObjs]/I");
  SimpleNtupleTree->Branch("pesView",&pesView,"pesView[numPesObjs]/I");
  SimpleNtupleTree->Branch("pesOctant",&pesOctant,"pesOctant[numPesObjs]/I");
  SimpleNtupleTree->Branch("pesCentroidStrip",&pesCentroidStrip,"pesCentroidStrip[numPesObjs]/D");
  SimpleNtupleTree->Branch("pesNumStrips",&pesNumStrips,"pesNumStrips[numPesObjs]/I");
  SimpleNtupleTree->Branch("pesLastStrip",&pesLastStrip,"pesLastStrip[numPesObjs]/I");
  SimpleNtupleTree->Branch("pesSeedStrip",&pesSeedStrip,"pesSeedStrip[numPesObjs]/I");
  SimpleNtupleTree->Branch("pesCentroidCM",&pesCentroidCM,"pesCentroidCM[numPesObjs]/D");
  SimpleNtupleTree->Branch("pesWidthCM",&pesWidthCM,"pesWidthCM[numPesObjs]/D");
  SimpleNtupleTree->Branch("pesClusterZ",&pesClusterZ,"pesClusterZ[numPesObjs]/D");

  SimpleNtupleTree->Branch("jetEtaDetector",&jetEtaDetector,"jetEtaDetector[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetGuardEnergy",&jetGuardEnergy,"jetGuardEnergy[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetEmFraction",&jetEmFraction,"jetEmFraction[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetCentroidEta",&jetCentroidEta,"jetCentroidEta[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetCentroidPhi",&jetCentroidPhi,"jetCentroidPhi[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetEtaMoment",&jetEtaMoment,"jetEtaMoment[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetPhiMoment",&jetPhiMoment,"jetPhiMoment[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetEtaPhiMoment",&jetEtaPhiMoment,"jetEtaPhiMoment[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetPx",&jetPx,"jetPx[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetPy",&jetPy,"jetPy[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetPz",&jetPz,"jetPz[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetTotEnergy",&jetTotEnergy,"jetTotEnergy[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetCentroidEt",&jetCentroidEt,"jetCentroidEt[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetTotP",&jetTotP,"jetTotP[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetTotPt",&jetTotPt,"jetTotPt[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetPt",&jetPt,"jetPt[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetPtSquared",&jetPtSquared,"jetPtSquared[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetTotEt",&jetTotEt,"jetTotEt[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetEt",&jetEt,"jetEt[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetEta",&jetEta,"jetEta[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetTheta",&jetTheta,"jetTheta[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetPhi",&jetPhi,"jetPhi[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetMassSquared",&jetMassSquared,"jetMassSquared[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetMass",&jetMass,"jetMass[numJetObjs]/D");
  SimpleNtupleTree->Branch("jetRapidity",&jetRapidity,"jetRapidity[numJetObjs]/D");

  std::cout << "Ntuple framework created." << std::endl;
  
  hfile->ls(); // just to verify what we have constructed...
  
  return AppResult::OK;
}

//----------------------------------------------------------------------------
// Begin filling quantities here...
//----------------------------------------------------------------------------
AppResult SimpleNtuple::event( AbsEvent* anEvent ){
  
  //  hfile->cd();

  numEventObjs = 0; numHepgObjs = 0; numVertexObjs = 0; numEmObjs = 0; 
  numMetObjs = 0; numPes2dObjs = 0; numPesObjs = 0; numJetObjs = 0;

  CdfTrackView_h theTracks;
  CdfTrackView::defTracks( theTracks );

  CdfEmObjectColl::const_handle emOb_coll;
  Pes2dClusterColl::const_handle pes2d_coll;
  Pes2dClusterColl::const_handle second_pes2d_coll;
  CdfMet_ch met_ch;
  PesClusterColl_ch pes_coll;
  CdfJetColl_ch jet_ch;
  PlugStripColl::const_handle plugStrip_coll;

  // first fill run and event numbers:

  runNumber[numEventObjs]   = AbsEnv::instance()->runNumber();
  eventNumber[numEventObjs] = AbsEnv::instance()->trigNumber();
  numEventObjs++;

  /////// Start with HEPG information, if there is any.../////////////////////

  EventRecord::ConstIterator hepgIter(anEvent, "HEPG_StorableBank");
  if(hepgIter.is_valid()) {
    const HEPG_StorableBank& hepg= *(ConstHandle<HEPG_StorableBank>(hepgIter));
    int nhep = hepg.n_particles();
    //    std::cout << "number of hepg particles = " << nhep << std::endl;
    for(HEPG_StorableBank::particleIter it(hepg); it.is_valid() ; ++it){
      int idpart = hepg.idhep(it);
      //      std::cout << "particle id = " << idpart << std::endl;
      int idmother = hepg.jmo1hep(it);
      //      std::cout << "mother particle number = " << idmother << std::endl;
      int idpart_mother = hepg.idhep((hepg.jmo1hep(it))-1);
      //      std::cout << "mother particle id = " << idpart_mother << std::endl;

      // NOTE:  The following is written specifically for W electrons.  You
      // must edit the code here and recompile if you are interested in 
      // something different.
      if (abs(idpart)==11 && abs(idpart_mother)==24) {
	HepLorentzVector v(hepg.Px(it),hepg.Py(it),hepg.Pz(it),hepg.E(it));
	hepgElectronEta[numHepgObjs] = v.eta();
	hepgElectronPhi[numHepgObjs] = v.phi();
	hepgElectronEt[numHepgObjs]  = v.et();
	hepgElectronPx[numHepgObjs]  = hepg.Px(it);
	hepgElectronPy[numHepgObjs]  = hepg.Py(it);
	hepgElectronPz[numHepgObjs]  = hepg.Pz(it);
	hepgElectronE[numHepgObjs]   = hepg.E(it);
	hepgElectronVx[numHepgObjs]  = hepg.Vx(it);
	hepgElectronVy[numHepgObjs]  = hepg.Vy(it);
	hepgElectronVz[numHepgObjs]  = hepg.Vz(it);
	hepgElectronVt[numHepgObjs]  = hepg.Vt(it);
	numHepgObjs++;
      }
    }
  }

  ///////////// Now get vertex information from vxprim if present...//////////

  for (EventRecord::ConstIterator iter(anEvent,"VertexColl");iter.is_valid();
       ++iter) {
    ConstHandle<VertexColl> vertexset(iter);
    Link<Vertex> primary = VertexColl::findPrimary(anEvent);
    if (primary.is_nonnull()) {
      HepPoint3D coord = primary->vertex();
      vertexPrimX[numVertexObjs] = coord.x();
      vertexPrimY[numVertexObjs] = coord.y();
      vertexPrimZ[numVertexObjs] = coord.z();
      numVertexObjs++;
    }
  }

  /////////////////////// Now start EmObject section...///////////////////////

//   double primary_lepton_Z0 = -1000.0;

//   if (CdfEmObjectColl::find(emOb_coll)) {
//     CdfEmObjectColl::const_iterator i;
//     for (i = emOb_coll->contents().begin();
// 	 i != emOb_coll->contents().end() ; ++i ) { 
//       const EmCluster* cluster = (*i).getEmCluster();
//       CdfEmObject::Track_link  seed_track = (*i).maxPtTrack();
//       if (seed_track.is_nonnull() && (*i).lshrIsValid(3)) {
// 	double lshr = (*i).lshr(3);
// 	double tiso = ElectronVariables::eleTiso(seed_track,theTracks,0.4,10.0,0.2);
// 	const Hep3Vector p  = seed_track->momentum();
// 	const HepDouble mag_of_p = p.Hep3Vector::mag();
// 	double eOp = (cluster->emEnergy())/mag_of_p;
// 	if ( (cluster->region()) == 0 &&
// 	     (cluster->emEt()) > 20.0 &&
// 	     (lshr < 0.2) &&
// 	     (cluster->hadEm()) < (.05 + (0.00045*(cluster->emEnergy()))) &&
// 	     ((cluster->emIsolationEt4())/(cluster->emEt())) < 0.1 &&
// 	     (0.5 < eOp) && (eOp < 2.0) &&
// 	     (tiso < 5.0) ) {
// 	  primary_lepton_Z0 = seed_track->z0() ;
// 	  //	  std::cout << "primary lepton Z0 = " << primary_lepton_Z0 << std::endl;
// 	}
//       }
//     }
//     //    std::cout << "got here..." << std::endl;
//   }

  if (CdfEmObjectColl::find(emOb_coll)) {
    CdfEmObjectColl::const_iterator i;
    for (i = emOb_coll->contents().begin();
	 i != emOb_coll->contents().end() ; ++i ) {

      const EmCluster* cluster = (*i).getEmCluster();
      CdfEmObject::Track_link  seed_track = (*i).maxPtTrack();
      const EmSeed& theSeed = cluster->seed();
      //      if (primary_lepton_Z0 != -1000.0 ) {
      //      cluster->setVertexZ(primary_lepton_Z0);
      //      }
      seedEnergy[numEmObjs]     = theSeed.emE();
      seedEt[numEmObjs]         = theSeed.emEt();
      region[numEmObjs]         = cluster->region();
      side[numEmObjs]           = cluster->side();
      etaDetector[numEmObjs]    = cluster->emEtaDetector();
      etaEvent[numEmObjs]       = cluster->emEtaEvent();
      etaRMS[numEmObjs]         = cluster->emEtaRMS();
      phi[numEmObjs]            = cluster->emPhi();
      phiRMS[numEmObjs]         = cluster->emPhiRMS();
      HepLorentzVector cal4v    = cluster->fourMomentum();
      emCalEx[numEmObjs]        = cal4v.px();
      emCalEy[numEmObjs]        = cal4v.py();
      emCalEz[numEmObjs]        = cal4v.pz();
      pemPesMatch[numEmObjs] = -9999.0;
      pemPesDeta[numEmObjs]  = -9999.0;
      pemPesDphi[numEmObjs]  = -9999.0;
      if (cluster->region()==1){
	if (Pes2dClusterColl::find(pes2d_coll)) {
	  Pes2dClusterColl::const_iterator j;
	  vector<double> distances;
	  distances.clear();
	  vector<double> deltaeta;
	  deltaeta.clear();
	  vector<double> deltaphi;
	  deltaphi.clear();
	  for (j = pes2d_coll->contents().begin();
	       j != pes2d_coll->contents().end() ; ++j ) {
	    double pesgeta = (*j).plug2dClusterEta();
	    double pesgphi = (*j).plug2dClusterPhi();
	    double deta = (cluster->pem3x3FitDetEta())-(pesgeta);
	    deltaeta.push_back(deta);
	    //std::cout << "delta eta: " << deta << std::endl;
	    double dphi = (cluster->pem3x3FitPhi())-(pesgphi);
	    if (dphi>3.1415927) {dphi = 6.2831854-dphi;}
	    if (dphi<-3.1415927) {dphi = -6.2831854-dphi;}
	    deltaphi.push_back(dphi);
	    //std::cout << "delta phi: " << dphi << std::endl;
	    double dist = sqrt( (deta*deta) + (dphi*dphi) );
	    distances.push_back(dist);
	    //std::cout << "match distance: " << dist << std::endl;
	  }
	  double mindist = *(min_element(distances.begin(),distances.end()));
	  pemPesMatch[numEmObjs] = mindist;
	  double mindeta = *(min_element(deltaeta.begin(),deltaeta.end()));
	  pemPesDeta[numEmObjs] = mindeta;
	  double mindphi = *(min_element(deltaphi.begin(),deltaphi.end()));
	  pemPesDphi[numEmObjs] = mindphi;
	}
      }
      zVertex[numEmObjs]        = cluster->vertexZ();
      zCentroid[numEmObjs]      = cluster->zCentroid();
      emEnergy[numEmObjs]       = cluster->emEnergy();
      emEt[numEmObjs]           = cluster->emEt();
      // CorrectedEt[numEmObjs]    = ElectronVariables::eleCorrectedEt(*i);
      hadEnergy[numEmObjs]      = cluster->hadEnergy();
      hadEt[numEmObjs]          = cluster->hadEt();
      totalEnergy[numEmObjs]    = cluster->totalEnergy();
      totalEt[numEmObjs]        = cluster->totalEt();
      hadEm[numEmObjs]          = cluster->hadEm();
      totalIso4[numEmObjs]      = cluster->totalIsolationEt4();
      totalIso7[numEmObjs]      = cluster->totalIsolationEt7();
      emIso4[numEmObjs]         = cluster->emIsolationEt4();
      emIso7[numEmObjs]         = cluster->emIsolationEt7();
      hadIso4[numEmObjs]        = cluster->hadIsolationEt4();
      hadIso7[numEmObjs]        = cluster->hadIsolationEt7();
      if (((*i).lshrIsValid(2))) {lshr2[numEmObjs] = (*i).lshr(2);}
      else {lshr2[numEmObjs] = -9999.0;}
      if (((*i).lshrIsValid(3))) {lshr[numEmObjs] = (*i).lshr(3);}
      else {lshr[numEmObjs] = -9999.0;}
      // eOverP[numEmObjs]         = ElectronVariables::eleEOverP(*i);
      xCes[numEmObjs]                 = -9999.0;
      chiW[numEmObjs]                 = -9999.0;
      deltaX[numEmObjs]               = -9999.0;
      zCes[numEmObjs]                 = -9999.0;
      chiS[numEmObjs]                 = -9999.0;
      deltaZ[numEmObjs]               = -9999.0;
      pesEmObjEta[numEmObjs]          = -9999.0;
      pesEmObjPhi[numEmObjs]          = -9999.0;
      pesTestEmObjMatching[numEmObjs] = -9999.0;
      pesEmObjWidth[numEmObjs]        = -9999.0;
      if (region[numEmObjs]==0) {
	ValueVectorLink<CesClusterColl> strip_cluster;
	ValueVectorLink<CesClusterColl> wire_cluster;
	wire_cluster = (*i).bestMatchingCesCluster(CesCluster::WIRE);
	if (wire_cluster.is_nonnull()) {
	  xCes[numEmObjs] = wire_cluster->fitted_z_strip();
	  chiW[numEmObjs] = wire_cluster->chi2();
	  if (seed_track) {
	    deltaX[numEmObjs] = (*i).maxPtTrackLocalCoord().x()
	      -wire_cluster->localCoord();
	  }
	}
	strip_cluster = (*i).bestMatchingCesCluster(CesCluster::STRIP);
	if (strip_cluster.is_nonnull()) {
	  zCes[numEmObjs] = strip_cluster->fitted_z_strip();
	  chiS[numEmObjs] = strip_cluster->chi2();
	  if (seed_track) {
	    deltaZ[numEmObjs] = (*i).maxPtTrackLocalCoord().z()
	      -strip_cluster->globalCoord();
	  }
	}
      }
      else if (region[numEmObjs]==1) {
	ValueVectorLink<Pes2dClusterColl> pes_cluster;
	numPes2dClusters[numEmObjs] = (*i).matchingPes2dClusters().size();
	pes_cluster = (*i).bestMatchingPes2dCluster();
	if (pes_cluster.is_nonnull()) {
	  pesEmObjEta[numEmObjs]   = pes_cluster->plug2dClusterEta();
	  pesEmObjPhi[numEmObjs]   = pes_cluster->plug2dClusterPhi();
	  double testDeltaEta = ((cluster->pem3x3FitDetEta())-(pes_cluster->plug2dClusterEta()));
	  double testDeltaPhi = ((cluster->pem3x3FitPhi())-(pes_cluster->plug2dClusterPhi()));
	  if (testDeltaPhi>3.1415927) {testDeltaPhi = 6.2831854-testDeltaPhi;}
	  if (testDeltaPhi<-3.1415927) {testDeltaPhi = -6.2831854-testDeltaPhi;}
	  double sqrtdist = sqrt(testDeltaEta*testDeltaEta + testDeltaPhi*testDeltaPhi);
	  pesTestEmObjMatching[numEmObjs] = sqrtdist;
	  pesEmObjWidth[numEmObjs] = pes_cluster->clusterWidth();
	}
      }
      // xCes[numEmObjs] = ElectronVariables::eleXces(*i);
      // deltaZ[numEmObjs] = ElectronVariables::eleDeltaZ(*i);
      if (seed_track){
	trackIso[numEmObjs] = ElectronVariables::eleTiso(seed_track,theTracks,0.4,10.0,0.2);
	trackPt[numEmObjs]  = seed_track->pt();
	charge[numEmObjs]   = seed_track->charge();
	const Hep3Vector p  = seed_track->momentum();
	const HepDouble mag_of_p = p.Hep3Vector::mag();
	eOverP[numEmObjs]   = (cluster->emEnergy())/mag_of_p;
      }
      else { 
	trackIso[numEmObjs] = -9999.0;
	trackPt[numEmObjs]  = -9999.0;
	charge[numEmObjs]   = -9999.0;
	eOverP[numEmObjs]   = -9999.0;
      }
      nTowers[numEmObjs] = cluster->nTowers();
      nTracks[numEmObjs] = (*i).numberTracks();
      pem3x3FitTowers[numEmObjs] = cluster->pem3x3FitTowers();
      pem3x3FitDetEta[numEmObjs] = cluster->pem3x3FitDetEta();
      pem3x3FitPhi[numEmObjs]    = cluster->pem3x3FitPhi();
      pem3x3EmEnergy[numEmObjs]  = cluster->pem3x3EmEnergy();
      pem3x3Chisq[numEmObjs]     = cluster->pem3x3Chisq();
      numEmObjs++ ;
    }
  }

  ////////////////////// Now start Met section...///////////////////////

  CdfMet::Error metStatus = CdfMet::find(met_ch);
  if(metStatus == CdfMet::ERROR){
    ERRLOG(ELerror,"SimpleNtuple.exe: ") << "Error: Could not find CdfMet" << endmsg;
    //return AppResult::ERROR;
    met[numMetObjs]      = -9999.0;
    metPhi[numMetObjs]   = -9999.0;
    metZV[numMetObjs]    = -9999.0;
    metESum[numMetObjs]  = -9999.0;
    metEtSum[numMetObjs] = -9999.0;
    metExSum[numMetObjs] = -9999.0;
    metEySum[numMetObjs] = -9999.0;
  }
  else {
    met[numMetObjs]      = met_ch->met();
    // why does Pasha do the following:
    // Metphi[numMetObjs] = TVector2::Phi_0_2pi(atan2(-met_ch->eySum(),-met_ch->exSum()));
    metPhi[numMetObjs]   = met_ch->phi();
    metZV[numMetObjs]    = met_ch->zVertex();
    metESum[numMetObjs]  = met_ch->eSum();
    metEtSum[numMetObjs] = met_ch->etSum();
    metExSum[numMetObjs] = -met_ch->exSum();
    metEySum[numMetObjs] = -met_ch->eySum();
  }
  numMetObjs++;
  
  ///////////////////// Now start pes2d section...////////////////////////////
  
  if (Pes2dClusterColl::find(second_pes2d_coll)) {
    Pes2dClusterColl::const_iterator iter;
    for (iter = second_pes2d_coll->contents().begin();
	 iter != second_pes2d_coll->contents().end() ; ++iter ) {
      pes2dWidth[numPes2dObjs]      = (*iter).clusterWidth();
      pes2dEnergy[numPes2dObjs]     = (*iter).energy();
      pes2dX[numPes2dObjs]          = (*iter).plug2dClusterX();
      pes2dY[numPes2dObjs]          = (*iter).plug2dClusterY();
      pes2dZ[numPes2dObjs]          = (*iter).plug2dClusterZ();
      pes2dGlblEta[numPes2dObjs]    = (*iter).plug2dClusterEta();
      pes2dGlblPhi[numPes2dObjs]    = (*iter).plug2dClusterPhi();
      pes2dUEnergy[numPes2dObjs]    = (*iter).plug2dUenergy();
      pes2dVEnergy[numPes2dObjs]    = (*iter).plug2dVenergy();
      pes2dNumUStrips[numPes2dObjs] = (*iter).plug2dUstrips();
      pes2dNumVStrips[numPes2dObjs] = (*iter).plug2dVstrips();
      numPes2dObjs++;
    }
  }
  else {
    pes2dWidth[numPes2dObjs]      = -9999.0;
    pes2dEnergy[numPes2dObjs]     = -9999.0;
    pes2dX[numPes2dObjs]          = -9999.0;
    pes2dY[numPes2dObjs]          = -9999.0;
    pes2dZ[numPes2dObjs]          = -9999.0;
    pes2dGlblEta[numPes2dObjs]    = -9999.0;
    pes2dGlblPhi[numPes2dObjs]    = -9999.0;
    pes2dEnergy[numPes2dObjs]     = -9999.0;
    pes2dUEnergy[numPes2dObjs]    = -9999.0;
    pes2dVEnergy[numPes2dObjs]    = -9999.0;
    pes2dNumUStrips[numPes2dObjs] = -9999.0;
    pes2dNumVStrips[numPes2dObjs] = -9999.0;
    numPes2dObjs++;
  }

  //////////////////////// Now start pes section...///////////////////////////

  if (PesClusterColl::find(pes_coll)) {
    PesClusterColl::const_iterator it;
    for (it = pes_coll->contents().begin();
	 it != pes_coll->contents().end() ; ++it ) {
      pesWidthStrips[numPesObjs]   = (*it).clusterWidth();
      pesEnergy[numPesObjs]        = (*it).energy();
      pesSide[numPesObjs]          = (*it).side();
      pesView[numPesObjs]          = (*it).view();
      pesOctant[numPesObjs]        = (*it).octant();
      pesCentroidStrip[numPesObjs] = (*it).centroid();
      pesNumStrips[numPesObjs]     = (*it).numStrips();
      pesLastStrip[numPesObjs]     = (*it).lastStrip();
      pesSeedStrip[numPesObjs]     = (*it).seedStrip();
      pesCentroidCM[numPesObjs]    = (*it).centroidCM();
      pesWidthCM[numPesObjs]       = (*it).widthCM();
      pesClusterZ[numPesObjs]      = (*it).clusterZ();
      numPesObjs++;
    }
  }
  else {
    pesWidthStrips[numPesObjs]   = -9999.0;
    pesEnergy[numPesObjs]        = -9999.0;
    pesSide[numPesObjs]          = -9999.0;
    pesView[numPesObjs]          = -9999.0;
    pesOctant[numPesObjs]        = -9999.0;
    pesCentroidStrip[numPesObjs] = -9999.0;
    pesNumStrips[numPesObjs]     = -9999.0;
    pesLastStrip[numPesObjs]     = -9999.0;
    pesSeedStrip[numPesObjs]     = -9999.0;
    pesCentroidCM[numPesObjs]    = -9999.0;
    pesWidthCM[numPesObjs]       = -9999.0;
    pesClusterZ[numPesObjs]      = -9999.0;
    numPesObjs++;
  }

  /////////////////////// Now start jet section...///////////////////////////

  std::string description = "JetCluModule-cone0.4";
  CdfJetColl::Error jetStatus = CdfJetColl::find(anEvent, jet_ch, description);
  if(jetStatus == CdfJetColl::ERROR){ 
    jetEtaDetector[numJetObjs]  = -9999.0;
    jetGuardEnergy[numJetObjs]  = -9999.0;
    jetEmFraction[numJetObjs]   = -9999.0;
    jetCentroidEta[numJetObjs]  = -9999.0;
    jetCentroidPhi[numJetObjs]  = -9999.0;
    // HepLorentzVector fourMomentum();
    jetEtaMoment[numJetObjs]    = -9999.0;
    jetPhiMoment[numJetObjs]    = -9999.0;
    jetEtaPhiMoment[numJetObjs] = -9999.0;
    jetPx[numJetObjs]           = -9999.0;
    jetPy[numJetObjs]           = -9999.0;
    jetPz[numJetObjs]           = -9999.0;
    jetTotEnergy[numJetObjs]    = -9999.0;
    jetCentroidEt[numJetObjs]   = -9999.0;
    jetTotP[numJetObjs]         = -9999.0;
    jetTotPt[numJetObjs]        = -9999.0;
    jetPt[numJetObjs]           = -9999.0;
    jetPtSquared[numJetObjs]    = -9999.0;
    jetTotEt[numJetObjs]        = -9999.0;
    jetEt[numJetObjs]           = -9999.0;
    jetEta[numJetObjs]          = -9999.0;
    jetTheta[numJetObjs]        = -9999.0;
    jetPhi[numJetObjs]          = -9999.0;
    jetMassSquared[numJetObjs]  = -9999.0;
    jetMass[numJetObjs]         = -9999.0;
    jetRapidity[numJetObjs]     = -9999.0;
    numJetObjs++;
  } 
  else {
    const JetList& jets = jet_ch->contents();
    ConstJetIter iter;
    for( iter = jets.begin() ; iter != jets.end() ; ++iter){
      jetEtaDetector[numJetObjs]  = iter->etaDetector();
      jetGuardEnergy[numJetObjs]  = iter->guardEnergy(); //sum of energy in twrs near edge of plug & twrs brdrng 30 & 90 degree cracks.
      jetEmFraction[numJetObjs]   = iter->emFraction(); //sum of emEnergy of the towers of the jet divided by totE of the jet.
      jetCentroidEta[numJetObjs]  = iter->centroidEta(); 
      jetCentroidPhi[numJetObjs]  = iter->centroidPhi();
      // HepLorentzVector fourMomentum();
      jetEtaMoment[numJetObjs]    = iter->etaEtaMoment(); //Measure the spreading out of energy distribution of the jets in eta.
      jetPhiMoment[numJetObjs]    = iter->phiPhiMoment(); //Same as above in phi.
      jetEtaPhiMoment[numJetObjs] = iter->etaPhiMoment(); //Same as above at 45 degrees in eta-phi.
      jetPx[numJetObjs]           = iter->px();
      jetPy[numJetObjs]           = iter->py();
      jetPz[numJetObjs]           = iter->pz();
      jetTotEnergy[numJetObjs]    = iter->totEnergy();
      jetCentroidEt[numJetObjs]   = iter->centroidEt();
      jetTotP[numJetObjs]         = iter->totP(); //sqrt(px^2 + py^2 + pz^2)
      jetTotPt[numJetObjs]        = iter->totPt(); //qrt(px^2 + py^2)
      jetPt[numJetObjs]           = iter->pt(); //Same as above.
      jetPtSquared[numJetObjs]    = iter->ptSquared();
      jetTotEt[numJetObjs]        = iter->totEt(); //totEnergy * (Pt / P)
      jetEt[numJetObjs]           = iter->et(); //Same as above.
      jetEta[numJetObjs]          = iter->eta(); //1/2 ln( totP + pz / totP - pz)
      jetTheta[numJetObjs]        = iter->theta(); //arctan( pt / pz)
      jetPhi[numJetObjs]          = iter->phi(); //arctan( py / px)
      jetMassSquared[numJetObjs]  = iter->massSquared(); //totE^2 - ( px^2 + py^2 + pz^2)
      jetMass[numJetObjs]         = iter->mass();
      jetRapidity[numJetObjs]     = iter->rapidity(); //1/2 ln( totE + pz / totE - pz)
      numJetObjs++;
    }
  }

  SimpleNtupleTree->Fill();
  return AppResult::OK;
}

//----------------------------------------------------------------------------
// End filling.
//----------------------------------------------------------------------------

AppResult SimpleNtuple::endRun( AbsEvent* aRun )
{
  std::cout << name( ) << " end Run" << std::endl;
  return AppResult::OK;
}

AppResult SimpleNtuple::endJob( AbsEvent* aJob )
{
  std::cout << name() << " end Job" << std::endl;
  hfile->Write();
  hfile->Close();
  return AppResult::OK;
}

AppResult SimpleNtuple::abortJob( AbsEvent* aJob )
{
  std::cout << name( ) << " abort Job" << std::endl;
  return AppResult::OK;
}

AppResult SimpleNtuple::abortRun( AbsEvent* aJob )
{
  std::cout << name( ) << " abort Run" << std::endl;
  return AppResult::OK;
}

SimpleNtuple::Command::Command() {
}

SimpleNtuple::Command::~Command() {
}

SimpleNtuple::Command::Command(const char* name, AppModule* module) : 
  APPCommand(name,module) {
}

void SimpleNtuple::Command::show() const {
}

bool SimpleNtuple::Command::isShowable() const {
  return true;
}

string SimpleNtuple::Command::description() const {
  string retval;
  return retval;
}

// 
// end
//