//---------------------------------------------//
//  Class:  XFT_Efficiency                       //
//---------------------------------------------//

#include "Electroweak/XFT_Efficiency.hh"

#include "HepTuple/HepFileManager.h"
#include "HepTuple/HepRootFileManager.h"
#include "HepTuple/HepRootNtuple.h"

#include "Extrapolator/AbsExtrapolator.hh"
#include "Extrapolator/MExtrapolator.hh"

#include "TriggerObjects/TL2D_StorableBank.hh"
#include "VertexObjects/Beamline.hh"

#ifdef USE_CDFEDM2
#include "TrackingObjects/Tracks/CdfTrack.hh"
#include "TrackingObjects/Storable/CdfTrackView.hh"
#include "TrackingObjects/Storable/CdfTrackColl.hh"
#else
#include "Banks/TRCK_Bank.hh"
#include "Trybos/TRY_Record_Iter_Same.hh"
#endif // USE_CDFEDM2

#ifdef DEFECT_OLD_IOSTREAM_HEADERS
#include <iostream>
#else
#include <iostream>
#  ifndef DEFECT_NO_STDLIB_NAMESPACES
using std::dec;
using std::cout;
using std::endl;
#  endif
#endif


XFT_Efficiency::XFT_Efficiency( const char* const theName,
                                  const char* const theDescription )
  : HepHistModule( theName, theDescription )
  , _event_number(0)
{
}


XFT_Efficiency::~XFT_Efficiency() {}

AppResult XFT_Efficiency::beginJob( AbsEvent* aJob ) 
{

  _d0 = new TH1F("d0","",80,-0.2,0.2);
  _z0 = new TH1F("z0","",200,-100.,100.);
  _p0 = new TH1F("p0","",64,0.,6.4);
  _ct = new TH1F("ct","",60,-1.5,1.5);
  _cu = new TH1F("cu","",100,-0.0002116,0.0002116);
  _pt = new TH1F("pt","",150,0.,50.);
  _c2 = new TH1F("c2","",100,0.,5.);
  _sl = new TH1F("sl","",8,0.,8.);
  _psl6 = new TH1F("psl6","",64,0.,6.4);
  _xftpt = new TH1F("xftpt","",150,0.,50.);
  _xftcu = new TH1F("xftcu","",100,-0.0002116,0.0002116);
  _xftptvpt = new TH2F("xftptvpt","",150,0.,50.,150,0.,50.);
  _xftcuvcu = new TH2F("xftcuvcu","",100,-0.0002116,0.0002116,
		       100,-0.002116,0.002116);
  _xftpsl6 = new TH1F("xftpsl6","",64,0.,6.4);
  _dpsl6 = new TH1F("dpsl6","",200,-0.2,0.2);
  _dpmatchsl6 = new TH1F("dpmatchsl6","",200,-0.1,0.1);

  // find beamline
  _beam = new CotBeam();
  return AppResult::OK;
}

AppResult XFT_Efficiency::endJob( AbsEvent* aJob ) 
{
  if (_beam)
    delete _beam;

  // Print out the counters.
  cout << "-------------------------------------------------------------" << endl;
  cout << "XFT_Efficiency statistics" << endl;
  cout << "Number of events: " << _event_number << endl;
  cout << "-------------------------------------------------------------" << endl;
  return AppResult::OK;
}

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

AppResult XFT_Efficiency::beginRun( AbsEvent* aRun ) 
{
  _beam->loadRun();
  if ( _beam->isValid() ) {
    Beamline userBeam = _beam->getBeamline();
  }
  return AppResult::OK;

}

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


AppResult XFT_Efficiency::event( AbsEvent* anEvent ) {

  // Initialize for this event
  ++_event_number;

  Beamline userBeam;
  if ( _beam->isValid() )
    userBeam = _beam->getBeamline();

  // Loop over tracks. 
  CdfTrackView_h hView;
  CdfTrackView::Error trackStatus ;
  trackStatus = CdfTrackView::defTracks(hView,"");
  CdfTrackView::CollType& tracks = hView->contents();

  for(CdfTrackView::const_iterator iter = tracks.begin(); 
      iter != tracks.end(); ++iter){
    const CdfTrack * trkPtr = iter->operator->();
    const CdfTrack& trkRef = **iter;

    // Look for COT stand-alone parent track
    CdfTrack_clnk cotParentTrack = 0; //initialize to null    
    bool foundCotParent = false;
    if(trkPtr->algorithm().value() == CdfTrack::CotStandAloneAlg){
      cotParentTrack=trkPtr;
      foundCotParent=true;
    }
    else if(trkPtr->algorithm().value()>CdfTrack::CotStandAloneAlg){
      cotParentTrack=trkPtr;
      while(cotParentTrack->algorithm().value() > CdfTrack::CotStandAloneAlg
	    && cotParentTrack->parent()) {
	cotParentTrack = cotParentTrack->parent();
      }
      if(cotParentTrack->algorithm().value()==CdfTrack::CotStandAloneAlg)
	foundCotParent=true;
    } // find COT parent

    if(foundCotParent){

      int nhtot = 0;
      int nhits[8];
      int nax = 0;
      int nst = 0;
      for ( int sl = 0 ; sl < 8 ; ++sl ){
	nhits[sl] = trkRef.numCTHits(sl);
	nhtot += nhits[sl];
	if (nhits[sl]>4) {
	  if (sl%2 == 0) nst++;
	  else nax++;
	}
      }

      float p0 = cotParentTrack->phi0();
      float z0 = cotParentTrack->z0();

      float beamx = userBeam.BeamX()+z0*userBeam.SlopeX();
      float beamy = userBeam.BeamY()+z0*userBeam.SlopeY();

      float d0 = cotParentTrack->d0() + beamx*sin(p0) - beamy*cos(p0);
      float cu = cotParentTrack->curvature();
      if (nax > 2 && nst > 2 && 
	  fabs(d0) < 0.2 && 
	  fabs(cu) < 0.0002116) {

	_cu->Fill(cu);
	_pt->Fill(fabs(0.002115939/cu));
	_p0->Fill(p0);
	_ct->Fill(cotParentTrack->lambda());
	_d0->Fill(d0);
	_z0->Fill(z0);
	_c2->Fill(cotParentTrack->chi2()/(nhtot-5));
	_sl->Fill(cotParentTrack->algorithm().CotSeedSl());

	// now find associated XFT track
	AbsExtrapolator* extrap;
	extrap = new MExtrapolator(CdfDetector::instance());
	AbsExtrapolator::Error status;

	//const CdfTrack& cotTrack = **cotParentTrack;
	Helix helix((*cotParentTrack).getTrajectory());

	HepPoint3D position =
	  helix.getPosition(helix.getPathLengthAtRhoEquals(d0));
	
	Hep3Vector momentum = (*cotParentTrack).momentum();

	// Initialize the extrapolator
	extrap->setInitialPosition(position);
	extrap->setInitialMomentum(momentum);
	if (cu < 0) {
	  extrap->setParticleCharge(-1);
	  extrap->setParticleType(AbsExtrapolator::MUM);
	}
	else {
	  extrap->setParticleCharge(1);
	  extrap->setParticleType(AbsExtrapolator::MUP);
	}

	HepPoint3D currentLocation;
	HepVector3D currentDirection;

	// Extrapolate
	float extrapolationRadius = 106.0; // from TDR, COT section, SL 6 radius
	float phi_sl6 = -999.;
	status = extrap->extrapolateToRadius(extrapolationRadius);
	if (status == AbsExtrapolator::OK) {

	  currentLocation = extrap->getCurrentPosition();
	  currentDirection = extrap->getCurrentMomentum();

	  phi_sl6 = currentLocation.phi();
	  if (phi_sl6 < 0.) phi_sl6 += 2.*PI;
	  _psl6->Fill(phi_sl6);
	}
	else {
	  std::cout << "Bad extrapolation!  d0: " << d0 << ", cu: " 
		    << cu << std::endl;
	  return AppResult::OK;
	}
	delete extrap;

	// Loop through the XFT tracks to find the closest
	// and next to closest in phi

	double deltaPhiBest = 999.;
	double deltaPhiNextBest = 999.;

	double phi1 = -999.;
	double phi2 = -999.;
	double curv1 = 999.;
	double curv2 = 999.;

	// 4-layer tracks
	double ptBinValue[48] = {
	  192.00, 64.00, 38.40, 27.42, 21.33, 17.45, 14.77,
	  12.80, 11.29, 10.10,  9.14,  8.34,  7.68,  7.11,
	  6.62,  6.19,  5.81,  5.48,  5.18,  4.92,  4.68,
	  4.46,  4.26,  4.09,  3.92,  3.76,  3.62,  3.49,
	  3.37,  3.25,  3.14,  3.04,  2.91,  2.74,  2.59,
	  2.46,  2.34,  2.23,  2.13,  2.04,  1.95,  1.88,
	  1.81,  1.74,  1.68,  1.63,  1.57,  1.52 }; 

	// Loop over TL2D banks
	for ( EventRecord::ConstIterator iTL2D(anEvent, "TL2D_StorableBank");
	      iTL2D.is_valid(); ++iTL2D) {

	  ConstHandle<TL2D_StorableBank> TL2D(iTL2D);
	  
	  int validTracks = TL2D->get_num_XTRP_tracks();

	  for( int i = 0; i < validTracks; i++ ) {

	    double phi = 2.0 * PI * (TL2D->get_XTRP_phi(i) / 2304.0);
	    int ptbin = TL2D->get_XTRP_track_curvature(i);

	    double ptValue = 0.0;
 
	    if(ptbin<48) {
	      ptValue = (-1.0) * ptBinValue[47-ptbin];
	    } else if(ptbin<96) {
	      ptValue = ptBinValue[ptbin-48];
	    } else {
	      std::cout << "XFT Bin out of range!  Bin #: " 
			<< ptbin << std::endl;
	      ptValue = 0.01;
	    }
	    _xftcu->Fill(0.00211593/ptValue);
	    _xftpsl6->Fill(phi);
	    _xftpt->Fill(fabs(ptValue));

	    double deltaPhi = phi_sl6 - phi;
	    while (deltaPhi > PI) deltaPhi -= 2.0 * PI;
	    while (deltaPhi < -PI) deltaPhi += 2.0 * PI;

	    _dpsl6->Fill(deltaPhi);

	    if (fabs(deltaPhi) < fabs(deltaPhiBest)) {
	      deltaPhiNextBest = deltaPhiBest;
	      phi2 = phi1;
	      curv2 = curv1;
	      deltaPhiBest = deltaPhi;
	      phi1 = phi;
	      curv1 = 0.00211593/ptValue;
	    }
	    else if (fabs(deltaPhi) < fabs(deltaPhiNextBest)) {
	      deltaPhiNextBest = deltaPhi;
	      phi2 = phi;
	      curv2 = 0.00211593/ptValue;
	    } // find best match
	  } // loop over XFT tracks

	  _dpmatchsl6->Fill(deltaPhiBest);

	  if (fabs(deltaPhiBest) < 0.01) {
	    _xftcuvcu->Fill(cu,curv1);
	    _xftptvpt->Fill(fabs(0.002115939/cu),
			    fabs(0.002115939/curv1));
	  } else {
	    // set failures to pt = 20 MeV
	    _xftcuvcu->Fill(cu,0.0021159539/0.1);
	    _xftptvpt->Fill(fabs(0.002115939/cu),
			    0.02115939);
	  }

	} // loop over TL2D	  
      } // good track requirements
    } // have COT parent
  } // loop over tracks

  return AppResult::OK;
}