//--------------------------------------------------------------------------
// File and Version Information:
// 
//
// Description:
//      Module used to validate Ces clustering
//      in the data.
//
// Environment:
//	Software developed for the CDFII Detector
//
// Author List:
//   02-03-2001      Michael Riveline		Original Author
//
//
//------------------------------------------------------------------------

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

//-------------
// C Headers --
//-------------
#include <assert.h>

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

//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "BaBar/Cdf.hh"
#include "BankTools/ToFromYbos.hh"
#include "BaBar/Cdf.hh"

#include "Edm/EventRecord.hh"
#include "Edm/ConstHandle.hh"
#include "Edm/Handle.hh"
#include "Trybos/TRY_Bank_Number.hh"
#include "ElectronObjects/EmSeed.hh"
#include "ElectronObjects/EmCluster.hh"
#include "StripChamberGeometry/HWScDetectorNode.hh"
#include "TrackingObjects/Storable/CdfTrackColl.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"


#include "AbsEnv/AbsEnv.hh"
#include "CalorObjects/CesClusterColl.hh"
#include "ShowerMax/StpDefs.hh"

//-------------------------------
// Hbook Classes 
//-------------------------------
#include "HepTuple/HepHBookFileManager.h"
#include "HepTuple/HepHist1D.h"
#include "HepTuple/HepHist2D.h"
#include <HepTuple/HepHBookNtuple.h>
#include <HepTuple/HepNtuple.h>
#include <HepTuple/HepHist1D.h>
#include <HepTuple/HepHist2D.h>
#include <HepTuple/Column.h>
#include <HepTuple/Block.h>

//-----------------------------------------------------------------------
// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
//-----------------------------------------------------------------------

// static const char rcsid[] = "$Id: CentralStripClusterTest.cc,v 1.10 2001/03/01 21:38:10 riveline Exp $";

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

CentralStripClusterTest::CentralStripClusterTest(
    const char* const theName,
    const char* const theDescription )
  : AppModule( theName, theDescription ),
    _stripPar( this ),
    _clustPar( this ),
    _paramSet( this ),
    //Other talk-to stuff
    _rcut("Rcut",this,10.),
    _histfile("histfile",this,"ces_test.hbook"),
    _calhistos("Cal_histograms",this,true),
    _trackhistos("Track_histograms",this,true)
  
{

  
  commands( )->append( &_stripPar._seedThrsParam );
  commands( )->append( &_stripPar._stripThrsParam );
  commands( )->append( &_stripPar._printStripsParam );
  commands( )->append( &_stripPar._numberStripsParam );
  commands( )->append( &_stripPar._numberWiresParam );
  commands( )->append( &_stripPar._collectionStrategyParam );
  commands( )->append( &_stripPar._ptcutParam );
  commands( )->append( &_stripPar._etaminParam );
  commands( )->append( &_stripPar._etamaxParam );

  commands( )->append( &_clustPar._seedClusterThrsParam );
  commands( )->append( &_clustPar._stripClusterThrsParam );
  commands( )->append( &_clustPar._printClusterParam );
  commands( )->append( &_clustPar._numberStripsClusterParam );
  commands( )->append( &_clustPar._numberWiresClusterParam );
  commands( )->append( &_clustPar._chi2cutParam );
  commands( )->append( &_clustPar._energycutParam );



  commands( )->append( &_paramSet._algName );
  commands( )->append( &_paramSet._towerEtaMin );
  commands( )->append( &_paramSet._towerEtaMax );
  commands( )->append( &_paramSet._towerEtMin );
  commands( )->append( &_paramSet._seedEMEtMin );
  commands( )->append( &_paramSet._seedHad2EMMax );
  commands( )->append( &_paramSet._daughterERatioMax );
  commands( )->append( &_paramSet._clusterEMEtMin );
  commands( )->append( &_paramSet._useTrigTwrHad );
  commands( )->append( &_paramSet._clusterHad2EMMax );
  commands( )->append( &_paramSet._overrideEMEtMin );
  commands( )->append( &_paramSet._deltaEtaCEMMax );
  commands( )->append( &_paramSet._deltaPhiCEMMax );
  commands( )->append( &_paramSet._deltaEtaPEMMax );
  commands( )->append( &_paramSet._deltaPhiPEMMax );  

  _rcut.addDescription("\tEta-Phi cut for track-CES matching");
  commands()->append(&_rcut);

  _histfile.addDescription("\tOutput root file name.");
  commands()->append(&_histfile);

  _calhistos.addDescription("\tDo you want CAL-CES comparison?");
  commands()->append(&_calhistos);
  
  _trackhistos.addDescription("\tDo you want Track-CES comparison?");
  commands()->append(&_trackhistos);



}


//--------------
// Destructor --
//--------------

CentralStripClusterTest::~CentralStripClusterTest( )
{
  cout << "Goodbye from CentralStripClusterTest" << endl;
  // Delete the PhysicsTowerDataMaker we created at beginJob time
  delete _ptdMaker; _ptdMaker = 0;
  delete _manager;
}

//--------------
// Operations --
//--------------

AppResult CentralStripClusterTest::beginJob( AbsEvent* aJob )
{


  _stripPar.print();
  _paramSet.print();
  _clustPar.print();
  cout << "Value of rcut: " << _rcut.value() << endl;
  cout << "Name of the histogram: " << _histfile.value() .c_str()<< endl;
  cout << "CAL comparison? " << _calhistos.value() << endl;
  cout << "Track comparison? " << _trackhistos.value() << endl;


  // Initialize the strip clustering
  doCluster.get_ces_const();

  // establish database connection
  //  _StripCol.initJob();



  // Add the monitoring histograms


  //Initialize HBOOK
  _manager = new HepHBookFileManager( _histfile.value().c_str() ); 
  assert( 0 != _manager );
  
  // Track-based histograms

  // Strip Energy
  _EsHisto = &_manager->hist1D( "E(CES,strips) (track based)", 50 , 0., 25.0, 10 );
  assert( 0 != _EsHisto );  

  // Strip Z position (local)
  _ZlsHisto = &_manager->hist1D( "Zl(CES,strips) (track based)", 50 , 0., 270.0, 11 );
  assert( 0 != _ZlsHisto );  

  // Strip Z position (global)
  _ZgsHisto = &_manager->hist1D( "Zg(CES,strips) (track based)", 70 , -270., 270.0, 12 );
  assert( 0 != _ZgsHisto );  

  // Chi2 of the fit (strip)
  _Chi2sHisto = &_manager->hist1D( "Chi2(CES,strips) (track based)", 40 , 0., 40.0, 13 );
  assert( 0 != _Chi2sHisto );  

  // Difference between strip Z position and track extrapolated Z position
  _DZsHisto = &_manager->hist1D( "DZ(CES,strips) (track based)", 70 , -10., 10.0, 14 );
  assert( 0 != _DZsHisto );  

  // 2D plot comparing Z position and Z track
  _ZtZsHisto = &_manager->hist2D( "Z(Track) vs Z(CES)", 100 , -240., 240.0, 100, -240., 240., 17 );
  assert( 0 != _ZtZsHisto );  

  // Ratio E(strips)/P(Track)
  _EsPHisto = &_manager->hist1D( "E(CES,strips)/P  (track based)", 50 , 0., 5.0, 15 );
  assert( 0 != _EsPHisto );  
  
  // Ratio E(strips)/E(Em Cluster)
  _EsEmHisto = &_manager->hist1D( "E(CES,strips)/E(Em)  (track based)", 50 , 0., 5.0, 16 );
  assert( 0 != _EsPHisto );  



  // Wire Energy
  _EwHisto = &_manager->hist1D( "E(CES,wires) (track based)", 50 , 0., 25.0, -10 );
  assert( 0 != _EwHisto );  

  // Wire X position (local)
  _XwHisto = &_manager->hist1D( "X(CES,wires) (track based)", 60 , -30., 30.0, -11 );
  assert( 0 != _XwHisto );  

  // Wire Phi position (global)
  _PhiwHisto = &_manager->hist1D( "Phi(CES,wires) (track based)", 70 , 0., 7.0, -12 );
  assert( 0 != _PhiwHisto );  

  // Chi2 of the fit (wire)
  _Chi2wHisto = &_manager->hist1D( "Chi2(CES,wires) (track based)", 40 , 0., 40.0, -13 );
  assert( 0 != _Chi2wHisto );  

  // Difference between strip Z position and track extrapolated Z position
  _DXwHisto = &_manager->hist1D( "DX(CES,wires) (track based)", 70 , -10., 10.0, -14 );
  assert( 0 != _DXwHisto );  

  // 2D plot comparing Phi position and Phi track
  _PhitPhiwHisto = &_manager->hist2D( "Phi(Track) vs Phi(CES) (track based)", 100 , 0., 7.0, 100, 0, 7, -17 );
  assert( 0 != _PhitPhiwHisto );  


  // Ratio E(wires)/P(Track)
  _EwPHisto = &_manager->hist1D( "E(CES,wires)/P (track based)", 50 , 0., 5.0, -15 );
  assert( 0 != _EwPHisto );  
  
  // Ratio E(wires)/E(Em Cluster)
  _EwEmHisto = &_manager->hist1D( "E(CES,wires)/E(Em) (track based)", 50 , 0., 5.0, -16 );
  assert( 0 != _EwPHisto );  


  // 1 cal histogram, just to check if everything is all right.


  // Ratio E(CAL)/P(track)
  _EcalPHisto = &_manager->hist1D( "E(cal)/P(track) (track based)", 50 , 0., 3.0, 20 );
  assert( 0 != _EcalPHisto );  


  // Histograms for unbiased CES clusters



  // Strip Energy
  _EsHisto_1 = &_manager->hist1D( "E(CES,strip) unbiased", 50 , 0., 25.0, 100 );
  assert( 0 != _EsHisto_1 );  

  // Strip Z position (local)
  _ZlsHisto_1 = &_manager->hist1D( "Zl(CES,strip) unbiased", 70 , 0., 270.0, 110 );
  assert( 0 != _ZlsHisto_1 );  

  // Strip Z position (global)
  _ZgsHisto_1 = &_manager->hist1D( "Zg(CES,strip) unbiased", 70 , -270., 270.0, 120 );
  assert( 0 != _ZgsHisto_1 );  

  // Chi2 of the fit (strip)
  _Chi2sHisto_1 = &_manager->hist1D( "Chi2(CES,strip) unbiased", 40 , 0., 40.0, 130 );
  assert( 0 != _Chi2sHisto_1 );  

  // Difference between strip Z position and track extrapolated Z position
  _DZsHisto_1 = &_manager->hist1D( "DZ(CES,strip) unbiased", 70 , -10., 10.0, 140 );
  assert( 0 != _DZsHisto_1 );  

  // Ratio E(strip)/P(Track)
  _EsPHisto_1 = &_manager->hist1D( "E(CES,strip)/P unbiased", 50 , 0., 5.0, 150 );
  assert( 0 != _EsPHisto_1 );  
  
  // Ratio E(strips)/E(Em Cluster)
  _EsEmHisto_1 = &_manager->hist1D( "E(CES,strips)/E(Em) unbiased", 50 , 0., 5.0, 160 );
  assert( 0 != _EsPHisto_1 );  

  // 2D plot comparing Z position and Z track
  _ZtZsHisto_1 = &_manager->hist2D( "Z(Track) vs Z(CES) unbiased", 100 , -240., 240.0, 100, -240., 240., 170 );
  assert( 0 != _ZtZsHisto_1 );  



  // Wire Energy
  _EwHisto_1 = &_manager->hist1D( "E(CES,wires) unbiased", 50 , 0., 25.0, -100 );
  assert( 0 != _EwHisto_1 );  

  // Wire X position (local)
  _XwHisto_1 = &_manager->hist1D( "X(CES,wires) unbiased", 60 , -30., 30.0, -110 );
  assert( 0 != _XwHisto_1 );  

  // Wire Phi position (global)
  _PhiwHisto_1 = &_manager->hist1D( "Phi(CES,wires) unbiased", 70 , 0., 7.0, -120 );
  assert( 0 != _PhiwHisto_1 );  

  // Chi2 of the fit (wire)
  _Chi2wHisto_1 = &_manager->hist1D( "Chi2(CES,wires) unbiased", 40 , 0., 40.0, -130 );
  assert( 0 != _Chi2wHisto_1 );  

  // Difference between strip Z position and track extrapolated Z position
  _DXwHisto_1 = &_manager->hist1D( "DX(CES,wires) unbiased", 100 , -10., 10.0, -140 );
  assert( 0 != _DXwHisto_1 );  


  // Ratio E(wires)/P(Track)
  _EwPHisto_1 = &_manager->hist1D( "E(CES,wires)/P unbiased", 50 , 0., 5.0, -150 );
  assert( 0 != _EwPHisto_1 );  
  
  // Ratio E(wires)/E(Em Cluster)
  _EwEmHisto_1 = &_manager->hist1D( "E(CES,wires)/E(Em) unbiased", 50 , 0., 5.0, -160 );
  assert( 0 != _EwPHisto_1 );  


  // 2D plot comparing Phi position and Phi track
  _PhitPhiwHisto_1 = &_manager->hist2D( "Phi(Track) vs Phi(CES) unbiased", 100 , 0., 7.0, 100, 0, 7.0, -170 );
  assert( 0 != _PhitPhiwHisto_1 );  





  // Ratio E(CAL)/P(track)
  _EcalPHisto_1 = &_manager->hist1D( "E(cal)/P(track) (unbiased CES)", 50 , 0., 3.0, 200 );
  assert( 0 != _EcalPHisto_1 );  


  // Declare our PhysicsTowerDataMaker which handles filling our
  // PhysicsTowerData collection for each event.
  _ptdMaker = new PhysicsTowerDataMaker();
  _clusterStrategy = new SeededEmClustering();


  return AppResult::OK;
}


AppResult CentralStripClusterTest::beginRun( AbsEvent* aRun )
{
  //  _StripCol.initRun(_verbose.value());
  return AppResult::OK;
}


AppResult CentralStripClusterTest::event( AbsEvent* anEvent )
{
  // Get Run and Event Number from AbsEnv
  
  AbsEnv* pEnv = AbsEnv::instance( );
  _nrun = pEnv->runNumber( );
  _nevt = pEnv->trigNumber( );

  // Is it MC or data?

  bool mcflag = pEnv->monteFlag();
  
  

  _debug = 0;
  


  //
  // Initialize cal clustering
  //


  
  std::string     collectionSelector = "DefaultPTD";
  if (PhysicsTowerData::find(_chTowerCollection, collectionSelector, anEvent) 
      ==  PhysicsTowerData::ERROR) {
    PhysicsTowerParams params;  
    _chTowerCollection = _ptdMaker->create(anEvent, params, collectionSelector);

  

  calorMakeView( _chTowerCollection, _clusterableTowers,
                 TotEtGreater(_paramSet.towerEtMin()) );
  }


  //
  // Initialize strip clustering
  //


  _StripCol.initCollection(anEvent);


  //
  // Store tracks which cross the CES in the view matchingTracks
  //

  CdfTrackView * MatchingList = new CdfTrackView;
  CdfTrackView::CollType & matchingTracks = MatchingList->contents();


  CdfTrackView_h  theTracks;
  CdfTrackView::defTracks( theTracks );  // get default selected tracks
  
  
  for (CdfTrackView::const_iterator i = theTracks->contents().begin(); i != theTracks->contents().end(); ++i) {
    
    const CdfTrack_lnk & theTrack = *i;
    
    if (theTrack->pt() > _stripPar.ptcut() && theTrack->pseudoRapidity() > _stripPar.etamin_track() && theTrack->pseudoRapidity() < _stripPar.etamax_track()) {
      int match = _StripCol.extrapolate_track( theTrack, mcflag );
      if (match) matchingTracks.push_back(theTrack);
    }
  }



  
  //
  // Check what collection strategy is used: 0 --> order strips by decreasing
  //                                               energy and use strips with
  //                                               highest energy as seeds
  //                                         1 --> use track based seed 
  //                                               (extrapolate track to CES
  //                                                and use strip corresponding
  //                                                to extrapolated track).
  //                                         2 --> use both
  //
  //





  if (_stripPar.collectionStrategy() == 1 || _stripPar.collectionStrategy()==2){
    
    
    // Loop over all tracks (need cut on track momentum/eta before doing the
    // clustering... not yet implemented).
    

    for (CdfTrackView::const_iterator i = matchingTracks.begin(); i != matchingTracks.end(); ++i) {
      
      const CdfTrack_lnk & theTrack = *i;
      float p_track = (theTrack->momentum()).mag();

      float zl_strip;
      float zg_strip;
      float e_strip;
      float chi2_strip;
      
      float x_wire;
      float phi_wire;
      float e_wire;
      float chi2_wire;

      
      _StripCol.extrapolate_track( theTrack, mcflag );
      
      
      
      _module = _StripCol.getModule();
      _barrel = _StripCol.getSide();
      _eta_global = _StripCol.Eta_global();
      _phi_global = _StripCol.Phi_global();
      _xloc = _StripCol.get_Xlocal();
      _zloc = _StripCol.get_Zlocal();

      
      float theta = 2 * atan (exp(-_eta_global));
      
      float z_track=0.;
      if (theta){
	z_track = HWScDetectorNode::STR_RADIUS/tan(theta);	
      }
      //
      // Create the strip clustering in the vector MyStrips
      //
      
      bool cess = false;
      bool cesw = false;
      
      const vector<ScStrip> MyStrips = _StripCol.get_track_based_strips( &_stripPar);
      if (MyStrips.size()){	
	// Get the extrapolated track position at the radius of the CES
	
	CesCluster myCluster = doCluster.do_cluster_strips( &MyStrips, &_clustPar, theTrack, _zloc );
	
	if (doCluster.is_valid()){
	  cess = true;
	  zl_strip=myCluster.fitted_position();
	  zg_strip=myCluster.global_position();
	  e_strip=myCluster.raw_energy();
	  chi2_strip=myCluster.chi2();

	  //	  if (e_strip > esmax[_module][_barrel]){
	  //	    esmax[_module][_barrel] = e_strip;
	  //	    etamax[_module][_barrel] = eta_strip;
	  //	  }


	  _ZlsHisto->accumulate(zl_strip);
	  _ZgsHisto->accumulate(zg_strip);
	  _EsHisto->accumulate(e_strip);
	  _Chi2sHisto->accumulate(chi2_strip);
	  _DZsHisto->accumulate(_zloc-zl_strip);
	  _EsPHisto->accumulate(e_strip/p_track);	  
	  _ZtZsHisto->accumulate(zg_strip,z_track);


	  if (_debug || _verbose.value()){
	    myCluster.print();
	  }	
	}	  
      }
	
      //
      // Wires clustering
      //
      
      const vector<ScWire> MyWires = _StripCol.get_track_based_wires( &_stripPar);
      if (MyWires.size()){	
	
	CesCluster myCluster = doCluster.do_cluster_wires( &MyWires, &_clustPar, theTrack, _xloc );	
	if (doCluster.is_valid()){
	  cesw = true;
	  x_wire=myCluster.fitted_position();
	  phi_wire=myCluster.global_position();
	  e_wire=myCluster.raw_energy();
	  chi2_wire=myCluster.chi2();
	  
	  //	  if (e_wire > ewmax[_module][_barrel]){
	  //	    ewmax[_module][_barrel] = e_wire;
	  //	    phimax[_module][_barrel] = phi_wire;	      
	  //	  }

	  _XwHisto->accumulate(x_wire);
	  _PhiwHisto->accumulate(phi_wire);
	  _EwHisto->accumulate(e_wire);
	  _Chi2wHisto->accumulate(chi2_wire);
	  _DXwHisto->accumulate(_xloc-x_wire);
	  _EwPHisto->accumulate(e_wire/p_track);
	  _PhitPhiwHisto->accumulate(phi_wire,_phi_global);	
	  
	  if (_debug || _verbose.value()){
	    myCluster.print();
	  }	
	}	  
      }
      
      
      //
      // Get track based CEM cluster
      //

      if (_calhistos.value()){
  
	if (cess && cesw){
	  EmCluster em_cluster = get_CEM_cluster( _eta_global , _phi_global);
	  
	  float e_cal=em_cluster.totalEnergy();
	  float et_em=em_cluster.emEt();
	  float hadem_cal=em_cluster.hadEm();
	  
	  if (et_em > _paramSet._clusterEMEtMin.value() && hadem_cal < _paramSet._clusterHad2EMMax.value()){
	    _EcalPHisto->accumulate(e_cal/p_track);
	    _EsEmHisto->accumulate(e_strip/e_cal);
	    _EwEmHisto->accumulate(e_wire/e_cal);
	  }
	}
      }
    }  
  } 
  


  if (_stripPar.collectionStrategy() == 0 || _stripPar.collectionStrategy()==2){
    
    //
    // Strategy #2  --> find all seed strips (with energy > seedThrs).
    //

    int ces_s=0;
    int ces_w=0;
    int  module_strip[100];
    int side_strip[100];
    float zl_strip[100];
    float zg_strip[100];
    float e_strip[100];
    float chi2_strip[100];

    int  module_wire[100];
    int side_wire[100];
    float x_wire[100];
    float phi_wire[100];
    float e_wire[100];
    float chi2_wire[100];

	
    bool cess = false;
    bool cesw = false;



    
    const vector<ScStrip> StripSeed = _StripCol.findStripSeed( &_stripPar );    
    if (StripSeed.size()){
      for(vector <ScStrip>::const_iterator it_strips = StripSeed.begin(); it_strips!= StripSeed.end(); it_strips++){
	
	const vector<ScStrip> MyStrips = _StripCol.get_seed_based_strips( &_stripPar, it_strips);
	if (MyStrips.size()){	
	  CesCluster myCluster = doCluster.do_cluster_strips( &MyStrips, &_clustPar);
	  if (doCluster.is_valid()){
	    cess = true;

	    module_strip[ces_s]=myCluster.module();
	    side_strip[ces_s]=myCluster.side();
	    zl_strip[ces_s]=myCluster.fitted_position();
	    zg_strip[ces_s]=myCluster.global_position();
	    e_strip[ces_s]=myCluster.raw_energy();
	    chi2_strip[ces_s]=myCluster.chi2();	

	    _ZlsHisto_1->accumulate(zl_strip[ces_s]);
	    _ZgsHisto_1->accumulate(zg_strip[ces_s]);
	    _EsHisto_1->accumulate(e_strip[ces_s]);
	    _Chi2sHisto_1->accumulate(chi2_strip[ces_s]);


	    ces_s++;
	    
	    if (_debug || _verbose.value()){
	      myCluster.print();
	    }
	  }
	}
      }
    }
    
    //
    // Wires
    //
    
    const vector<ScWire> WireSeed = _StripCol.findWireSeed( &_stripPar );
    
    if (WireSeed.size()){
      for(vector <ScWire>::const_iterator it_wires = WireSeed.begin(); it_wires!= WireSeed.end(); it_wires++){
	
	
	const vector<ScWire> MyWires = _StripCol.get_seed_based_wires( &_stripPar, it_wires);
	if (MyWires.size()){	
	  CesCluster myCluster = doCluster.do_cluster_wires( &MyWires, &_clustPar );
	  if (doCluster.is_valid()){
	    cesw = true;

	    module_wire[ces_w]=myCluster.module();
	    side_wire[ces_w]=myCluster.side();
	    x_wire[ces_w]=myCluster.fitted_position();
	    phi_wire[ces_w]=myCluster.global_position();
	    e_wire[ces_w]=myCluster.raw_energy();
	    chi2_wire[ces_w]=myCluster.chi2();

	    _XwHisto_1->accumulate(x_wire[ces_w]);
	    _PhiwHisto_1->accumulate(phi_wire[ces_w]);
	    _EwHisto_1->accumulate(e_wire[ces_w]);
	    _Chi2wHisto_1->accumulate(chi2_wire[ces_w]);

	    ces_w++;	   
	    
	    if (_debug || _verbose.value()){
	      myCluster.print();
	    }	  
	  }
	}
      }
    }


    //
    // Loop over all clusters. Find tracks and EM cluster close by
    // Select the highest energy cluster in each wedge and compare it to
    // the corresponding EM cluster.
    //


    float esmax[24][2];
    float ewmax[24][2];
    float etamax[24][2];
    float phimax[24][2];
    float pmax[24][2];

    for (int mod=0;mod<24;mod++){
      for(int sid=0;sid<2;sid++){
	esmax[mod][sid]=0.;
	ewmax[mod][sid]=0.;
	pmax[mod][sid]=0.;
	etamax[mod][sid]=0.;
	phimax[mod][sid]=0.;
      }
    }

    if (cess && cesw){
      for(int is=0; is<ces_s; is++){
	for(int iw=0; iw<ces_w; iw++){	


	  // Check that both clusters are in the same module/side
	  if (module_wire[iw] == module_strip[is]){
	    if (side_wire[iw] == side_strip[is]){



	      // Calculate eta and phi
	      float phi = phi_wire[iw];
	      float theta=0;
	      float eta =0;
	      float rces = float(HWScDetectorNode::STR_RADIUS);
	      theta = atan2(rces,zg_strip[is]);
	      //
	      // This doesn't work
	      // theta = atan(rces/zg_strip[is]);
	      //
	      if (theta < 0) theta += 2* M_PI;
	      eta = -log(tan(theta/2));


	      int mod = module_strip[is];
	      int sid = side_strip[is];


	      // Get the maximum energy 

	      if (e_strip[is] > esmax[mod][sid]){
		esmax[mod][sid] = e_strip[is];
	        etamax[mod][sid] = eta;
	      }

	      if (e_wire[iw] > ewmax[mod][sid]){
		ewmax[mod][sid] = e_wire[iw];
	        phimax[mod][sid] = phi;
	      }


	      //
	      // Track comparison
	      //
	      
	      
	      if (_trackhistos.value()){
		
		// Distance track-cluster in eta-phi space
		float delta_track_ces=1000;
		float delta_z = 1000;
		float delta_x = 1000;
		float zt1=1000;
		float phit1=1000;
		float p_track=1000.;
		
		//
		// Loop over all tracks crossing CES and select the track
		// closest to the CES cluster in eta-phi
		//
		
		for (CdfTrackView::const_iterator i = matchingTracks.begin(); i != matchingTracks.end(); ++i) {
		  
		  const CdfTrack_lnk & theTrack = *i;
		  float p_int=(theTrack->momentum()).mag();
		  
		  // Get the maximum p in the module/side (for cal comparison)
		  if (_calhistos.value()){
		    if (p_int > pmax[mod][sid]){
		      pmax[mod][sid]=p_int;
		    }
		  }
		  
		  _StripCol.extrapolate_track( theTrack, mcflag );
		  
		  if (_StripCol.getModule() == module_strip[is] && 
		      _StripCol.getSide() == side_strip[is]){
		    float eta_track = _StripCol.Eta_global();
		    float theta_track = 2 * atan (exp(-eta_track));
		    
		    float xloc = _StripCol.get_Xlocal();
		    float zloc = _StripCol.get_Zlocal();
		    
		    
		    float z_track=0.;
		    if (theta_track){
		      z_track = HWScDetectorNode::STR_RADIUS/tan(theta_track);
		    }
		    
		    
		    float phi_track = _StripCol.Phi_global();
		    
		    float delta_int = sqrt( (eta_track-eta)*(eta_track-eta) + (phi_track -phi)*(phi_track-phi));
		    if (delta_int < delta_track_ces){
		      delta_track_ces = delta_int;
		      delta_z = zloc - zl_strip[is];
		      delta_x = xloc - x_wire[iw];
		      zt1 = z_track;
		      phit1 = phi_track;		      
		      p_track = p_int;
		    }
		  }   // if module(track) ==  module(ces)
		} // loop over matching tracks
		
		//
		// Fill track_CES histograms
		//
		
		if (delta_track_ces < _rcut.value()){
		  _DZsHisto_1 ->accumulate(delta_z);
		  _DXwHisto_1 ->accumulate(delta_x);
		  _ZtZsHisto_1 ->accumulate(zg_strip[is],zt1);
		  _PhitPhiwHisto_1 ->accumulate(phi_wire[iw],phit1);
		  _EsPHisto_1->accumulate(e_strip[is]/p_track);
		  _EwPHisto_1->accumulate(e_wire[iw]/p_track);
		}
	      } // if statement for track histos
	    }   // side_wire=side_strip
	  }     // module_wire=module_strip
	}       //loop on strips
      }         // loop on wires

      // Do we really want to get the cal histograms?
      
      if (_calhistos.value()){
	for (int mod=0; mod<24; mod++){
	  for(int sid=0; sid<2; sid++){
	    if (etamax[mod][sid]){
	      EmCluster em_cluster = get_CEM_cluster(etamax[mod][sid],phimax[mod][sid]);
	      
	      float e_cal=em_cluster.totalEnergy();
	      float et_em=em_cluster.emEt();
	      float hadem_cal=em_cluster.hadEm();
	      
	      if (et_em > _paramSet._clusterEMEtMin.value()   && hadem_cal < _paramSet._clusterHad2EMMax.value()){
		_EsEmHisto_1->accumulate(esmax[mod][sid]/e_cal);
		_EwEmHisto_1->accumulate(ewmax[mod][sid]/e_cal);
		if (_trackhistos.value()){
		  if (pmax[mod][sid])
		    _EcalPHisto_1->accumulate(e_cal/pmax[mod][sid]);
		}
	      }
	    }
	  }
	}
      }
    }
  }
  
  // Destroy the tracks matching CES
  MatchingList->destroy();
  return AppResult::OK;
}

 

AppResult CentralStripClusterTest::endRun( AbsEvent* aRun )
{
  return AppResult::OK;
}

AppResult CentralStripClusterTest::endJob( AbsEvent* aJob )
{
  return AppResult::OK;
}

AppResult CentralStripClusterTest::abortJob( AbsEvent* aJob )
{
   return AppResult::OK;
}

EmCluster CentralStripClusterTest::get_CEM_cluster(float eta_seed, float phi_seed)
{
  //
  // First find out which tower has been hit (atower)
  //
  
  size_t ieta = _geom.iEta(eta_seed);
  size_t iphi = _geom.iPhi(ieta,phi_seed);
  const PhysicsTower*      atower;

  // Initialize EM cluster to 0
  EmCluster em_cluster;
  
  
  TowerKey seedKey( ieta, iphi );
  Towers::const_iterator  ciTowers =
    std::find_if( _clusterableTowers.begin(), _clusterableTowers.end(),
		  TowerKey_eq( seedKey ));
  if (ciTowers!= _clusterableTowers.end()){
    atower = *ciTowers;
    
    // Calculate the Hadron/EM energy ratio
    // Guard against divide by zero error if EM energy is zero
    double had2em;
    if (atower->emEnergy( ) > 0.0) {
      had2em = atower->hadEnergy( ) / atower->emEnergy( );
    } else {
      had2em = 0.0;
    }

    
    // Store this tower as an EM seed if it is within the eta limits for
    // clustering and if it passes EM energy cut and Had/EM cut
    if ( atower->iEta() >= _paramSet.towerEtaMin()  &&
	 atower->iEta() <= _paramSet.towerEtaMax()  &&
	 atower->emEt( ) > _paramSet.seedEMEtMin()  &&
	 had2em < _paramSet.seedHad2EMMax() ) {
      EmSeed     new_seed( atower );
      //new_seed.setVertexZ( towers.vertexZ() );
      
      EmCluster new_cluster( new_seed );
      PhysicsTowerView                usedTowers;
      
      _clusterStrategy->build(new_cluster, _clusterableTowers,
			      usedTowers, _paramSet);
      em_cluster = EmCluster(new_cluster);
    }
  }
  return(em_cluster);
    
}