// $Header: TrackRecSmear.cxx $

#include "TrackRecSmear.h"
#include "SmearTrack.h"
#include "math.h"

ClassImp(TrackRecSmear)

//----------------------------------------------------------------------
//        constructor: create everything
 TrackRecSmear::TrackRecSmear(GetParameters* gp) {

  m_recon = new SmearTrack(gp);
  m_parameters = gp;
};

//----------------------------------------------------------------------
//        execute: 'reconstruct' the tracks from McPart
void
 TrackRecSmear::doit(Event* event)
{

  m_event = event;

  TObjArray* McPList = m_event->MCparticles();
  TObjArray* TrackerList = m_event->Tracker();
  TObjArray* VXDList = m_event->VXD();

  m_event->Tracks()->Delete();

 // create tracks for particles within the acceptance
	
  Int_t nparts = McPList->GetEntries();

  McPart* mcP;

  for( int part=0; part< nparts; part++) {
    mcP = (McPart*)McPList->At(part);
    int statusCode = mcP->GetStatus();

    // these are Gismo status codes: 8 = PRIMARY; 1 = DECAYED.

    int decayed = ( statusCode==8 || statusCode==1 ); 

    if ( decayed ) continue;

    float* p = mcP->GetMomentum();
    if ( mcP->GetCharge() == 0 ) continue;

    float pTot = 0.;
    for ( int i=0; i < 3; i++) {
      pTot += p[i]*p[i];

    }
    pTot = sqrt(pTot);
    float cos_theta = p[2]/pTot;
    float abscos = fabs(cos_theta);

    // find out how many hits in the VXD and Tracker come from this mcP 

    int nHits = 0;
    McPart* mcPtrkr;
    Int_t index = -1;

    int nTrkHits = TrackerList->GetEntries();

    for (int trkr=0; trkr < nTrkHits ; trkr++) {
      Tracker_Hit* tkHit = (Tracker_Hit*)TrackerList->At(trkr);
      index = tkHit->GetParticle();
      mcPtrkr = (McPart*)McPList->At(index);

      if ( mcPtrkr == mcP ) nHits += 1;
    }

    int nVXDHits = VXDList->GetEntries();
    for (int vxd=0; vxd < nVXDHits ; vxd++) {
      VXD_Hit* vxdHit = (VXD_Hit*)VXDList->At(vxd);
      index = vxdHit->GetParticle();
      mcPtrkr = (McPart*)McPList->At(index);

      if ( mcPtrkr == mcP ) nHits += 1;
    }

    if ( (nHits >= m_parameters->getMinHit()) && 
	    ( pTot > m_parameters->getPTMin()) ) {
      m_event->Tracks()->Add(m_recon->makeTrack(mcP,index));

    }
  }

}


//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void 
 TrackRecSmear::spew(FILE* ofile)const
//-------------------------------
{

// now write out the tracks

  fprintf(ofile,"TrackSystem n %i n", m_event->Tracks()->GetEntries());
  TIter next(m_event->Tracks());
  Track* tk;
  while(tk = (Track*)next() ) {
    fprintf(ofile,"%f %f %f %f %f %f %f %in",
          tk->GetMomentum()[0],
          tk->GetMomentum()[1],
          tk->GetMomentum()[2],
          tk->GetPosition()[0],
          tk->GetPosition()[1],
          tk->GetPosition()[2],
          tk->GetCharge(),
          (int)(tk->GetParticle()));

  }
  fprintf(ofile,"endn");


}