// $Header: TrackFullSmear.cxx $

#include "TrackFullSmear.h"
#include "SmearTrack.h"
#include "math.h"

//______________________________________________________________________
//
// TrackFullSmear
//
// oversees creation of tracks by applying acceptance cuts on McParts

ClassImp(TrackFullSmear)
 TrackFullSmear::TrackFullSmear(GetParameters* gp) {
  // constructor: sets up track smearer

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

//----------------------------------------------------------------------
//        execute: 'reconstruct' the tracks from McPart
void
 TrackFullSmear::doit(Event* event)
{
  // apply acceptance cuts on appropriate McParts and get SmearTrack t
  // make the tracks

  m_event = event;

  TObjArray* McPList = event->MCparticles();
  TObjArray* TkList = event->Tracks();

  TkList->Delete();

 // create tracks for particles within the acceptance
	
  TIter next(McPList);
  McPart* mcP;

  while( mcP = (McPart*)next() ) {

    int statusCode = mcP->GetStatus();

    int decayed = ( 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);
    int index = McPList->IndexOf(mcP);
    if ( (abscos < m_parameters->getPolarOuter()) && 
	    ( pTot > m_parameters->getPTMin() ) ) {
      TkList->Add(m_recon->makeTrack(mcP,index));
    }
  }


}


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

// writes out the tracks to output file

  fprintf(ofile,"Tracks 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");


}