THcHallCSpectrometer.cxx

/** \class THcHallCSpectrometer
    \ingroup Base

\brief A standard Hall C spectrometer apparatus

 Uses the standard optics polynomials to trace back to target

 Contains no standard detectors.  All detectors must be added.

 The usual name of this object is either "H", "S", "P"
 for HMS, SOS, or suPerHMS respectively

 In method ReadDatabase  calls the THaSpectrometer::SetCentralAngles(th_geo,ph_geo,bend_down)
  bend_down is flag that should equal kFALSE for Hall C spectrometers
  th_geo and ph_geo in degrees.
  Define the LAB coordinate system so the +Z points to beam dump, +Y is up, +X is beam left
  th_geo is rotation about the Y axis in XZ plane to coordinates X',Y=Y',Z'
  After th_geo rotation, ph_geo rotation is about the X' axis in the  Y'Z' plane.
  Calls THaAnalysisObject::GeoToSph to calculate the spherical angles. th_sph and ph_sph
  th_sph is rotation about the Y axis in XZ plane to coordinates X',Y=Y',Z
  After th_sph rotation, ph_sph rotation is about the original Z axis
  In Lab coordinate system:
       X = r*sin(th_geo)*cos(ph_geo)          X = r*sin(th_sph)*cos(ph_sph)
       Y = r*sin(ph_geo)                      Y = r*sin(th_sph)*sin(ph_sph)
       Z = r*cos(th_geo)*cos(ph_geo)          Z = r*cos(th_sph)

      cos(th_sph) = cos(th_geo)*cos(ph_geo)
      cos(ph_sph) = sin(th_geo)*cos(ph_geo)/sqrt(1-cos^2(th_geo)*cos^2(ph_geo))

      GeoToSph is coded so that 
       1) negative th_geo and ph_geo = 0 returns th_sph=abs(th_geo) and ph_sph =180
       2) positive th_geo and ph_geo = 0 returns th_sph=th_geo and ph_sph =0

   Using the spherical angles, the TRotation fToLabRot and inverse fToTraRot are calculated
    fToLabRot is rotation matrix from the spectrometer TRANSPORT system to Lab system
     TRANSPORT coordinates are +X_tra points vertically down, +Z_tra is along the central ray and +Y_tra = ZxX
     For ph_sph = 0   X_lab = Y_tra*cos(th_sph) + Z_tra*sin(th_sph)
                      Y_lab = -X_tra
                      Z_lab = -Y_tra*sin(th_sph) + Z_tra*cos(th_sph)
     For ph_sph = 180 X_lab = Y_tra*cos(th_sph) - Z_tra*sin(th_sph)
                      Y_lab = -X_tra
                      Z_lab = Y_tra*sin(th_sph) + Z_tra*cos(th_sph)


                      

  

\author S. A. Wood

*/
//////////////////////////////////////////////////////////////////////////

#include "THcHallCSpectrometer.h"
#include "THaTrackingDetector.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "THaTrack.h"
#include "THaTrackProj.h"
#include "THaTriggerTime.h"
#include "TMath.h"
#include "TList.h"

#include "THcRawShowerHit.h"
#include "THcSignalHit.h"
#include "THcShower.h"
#include "THcHitList.h"
#include "THcHodoscope.h"

#include <vector>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>

using std::vector;
using namespace std;

//_____________________________________________________________________________
THcHallCSpectrometer::THcHallCSpectrometer( const char* name, const char* description ) :
  THaSpectrometer( name, description )
{
  // Constructor. Defines the standard detectors for the HRS.
  //  AddDetector( new THaTriggerTime("trg","Trigger-based time offset"));

  //sc_ref = static_cast<THaScintillator*>(GetDetector("s1"));

  SetTrSorting(kTRUE);
}

//_____________________________________________________________________________
THcHallCSpectrometer::~THcHallCSpectrometer()
{
  // Destructor

  DefineVariables( kDelete );
}

//_____________________________________________________________________________
Int_t THcHallCSpectrometer::DefineVariables( EMode mode )
{
  // Define/delete standard variables for a spectrometer (tracks etc.)
  // Can be overridden or extended by derived (actual) apparatuses
  if( mode == kDefine && fIsSetup ) return kOK;
  THaSpectrometer::DefineVariables( mode );
  fIsSetup = ( mode == kDefine );
  RVarDef vars[] = {
    { "tr.betachisq", "Chi2 of beta", "fTracks.THaTrack.GetBetaChi2()"},
    { 0 }
  };


  return DefineVarsFromList( vars, mode );
}

//_____________________________________________________________________________
Bool_t THcHallCSpectrometer::SetTrSorting( Bool_t set )
{
  if( set )
    fProperties |= kSortTracks;
  else
    fProperties &= ~kSortTracks;

  return set;
}

//_____________________________________________________________________________
Bool_t THcHallCSpectrometer::GetTrSorting() const
{
  return ((fProperties & kSortTracks) != 0);
}

//_____________________________________________________________________________
void THcHallCSpectrometer::InitializeReconstruction()
{
  fNReconTerms = 0;
  fReconTerms.clear();
  fAngSlope_x = 0.0;
  fAngSlope_y = 0.0;
  fAngOffset_x = 0.0;
  fAngOffset_y = 0.0;
  fDetOffset_x = 0.0;
  fDetOffset_y = 0.0;
  fZTrueFocus = 0.0;
}
//_____________________________________________________________________________
Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
{

  static const char* const here = "THcHallCSpectrometer::ReadDatabase";

#ifdef WITH_DEBUG
  cout << "In THcHallCSpectrometer::ReadDatabase()" << endl;
#endif

  const char* detector_name = "hod";
  //THaApparatus* app = GetDetector();
  THaDetector* det = GetDetector("hod");
  // THaDetector* det = app->GetDetector( shower_detector_name );

  if( dynamic_cast<THcHodoscope*>(det) ) {
    fHodo = static_cast<THcHodoscope*>(det);     // fHodo is a membervariable
  } else {
    Error("THcHallCSpectrometer", "Cannot find hodoscope detector %s",
    	  detector_name );
    fHodo = NULL;
  }


  THaDetector* detdc = GetDetector("dc");
  if( dynamic_cast<THcDC*>(detdc) ) {
    fDC = static_cast<THcDC*>(detdc);     // fHodo is a membervariable
  } else {
    Error("THcHallCSpectrometer", "Cannot find detector DC");
    fDC = NULL;
  }


  // Get the matrix element filename from the variable store
  // Read in the matrix

  InitializeReconstruction();

  char prefix[2];

#ifdef WITH_DEBUG
  cout << " GetName() " << GetName() << endl;
#endif
  
  prefix[0]=tolower(GetName()[0]);
  prefix[1]='\0';

  string reconCoeffFilename;
  DBRequest list[]={
    {"_recon_coeff_filename", &reconCoeffFilename,     kString               },
    {"theta_offset",          &fThetaOffset,           kDouble               },
    {"phi_offset",            &fPhiOffset,             kDouble               },
    {"delta_offset",          &fDeltaOffset,           kDouble               },
    {"thetacentral_offset",   &fThetaCentralOffset,    kDouble               },
    {"_oopcentral_offset",    &fOopCentralOffset,      kDouble               },
    {"pcentral_offset",       &fPCentralOffset,        kDouble               },
    {"pcentral",              &fPcentral,              kDouble               },
    {"theta_lab",             &fTheta_lab,             kDouble               },
    {"partmass",              &fPartMass,              kDouble               },
    {"phi_lab",               &fPhi_lab,               kDouble,         0,  1},
    {"sel_using_scin",        &fSelUsingScin,          kInt,            0,  1},
    {"sel_using_prune",       &fSelUsingPrune,         kInt,            0,  1},
    {"sel_ndegreesmin",       &fSelNDegreesMin,        kDouble,         0,  1},
    {"sel_dedx1min",          &fSeldEdX1Min,           kDouble,         0,  1},
    {"sel_dedx1max",          &fSeldEdX1Max,           kDouble,         0,  1},
    {"sel_betamin",           &fSelBetaMin,            kDouble,         0,  1},
    {"sel_betamax",           &fSelBetaMax,            kDouble,         0,  1},
    {"sel_etmin",             &fSelEtMin,              kDouble,         0,  1},
    {"sel_etmax",             &fSelEtMax,              kDouble,         0,  1},
    {"hodo_num_planes",       &fNPlanes,               kInt                  },
    {"scin_2x_zpos",          &fScin2XZpos,            kDouble,         0,  1},
    {"scin_2x_dzpos",         &fScin2XdZpos,           kDouble,         0,  1},
    {"scin_2y_zpos",          &fScin2YZpos,            kDouble,         0,  1},
    {"scin_2y_dzpos",         &fScin2YdZpos,           kDouble,         0,  1},
    {"prune_xp",              &fPruneXp,               kDouble,         0,  1},
    {"prune_yp",              &fPruneYp,               kDouble,         0,  1},
    {"prune_ytar",            &fPruneYtar,             kDouble,         0,  1},
    {"prune_delta",           &fPruneDelta,            kDouble,         0,  1},
    {"prune_beta",            &fPruneBeta,             kDouble,         0,  1},
    {"prune_df",              &fPruneDf,               kDouble,         0,  1},
    {"prune_chibeta",         &fPruneChiBeta,          kDouble,         0,  1},
    {"prune_npmt",            &fPruneNPMT,           kDouble,         0,  1},
    {"prune_fptime",          &fPruneFpTime,             kDouble,         0,  1},
    {0}
  };

  // Default values
  fSelUsingScin = 0;
  fSelUsingPrune = 0;
  fPhi_lab = 0.;
  gHcParms->LoadParmValues((DBRequest*)&list,prefix);

  EnforcePruneLimits();

#ifdef WITH_DEBUG
  cout <<  "\n\n\nhodo planes = " <<  fNPlanes << endl;
  cout <<  "sel using scin = "    <<  fSelUsingScin << endl;
  cout <<  "fPruneXp = "          <<  fPruneXp << endl;
  cout <<  "fPruneYp = "          <<  fPruneYp << endl;
  cout <<  "fPruneYtar = "        <<  fPruneYtar << endl;
  cout <<  "fPruneDelta = "       <<  fPruneDelta << endl;
  cout <<  "fPruneBeta = "        <<  fPruneBeta << endl;
  cout <<  "fPruneDf = "          <<  fPruneDf << endl;
  cout <<  "fPruneChiBeta = "     <<  fPruneChiBeta << endl;
  cout <<  "fPruneFpTime = "      <<  fPruneFpTime << endl;
  cout <<  "fPruneNPMT = "        <<  fPruneNPMT << endl;
  cout <<  "sel using prune = "   <<  fSelUsingPrune << endl;
#endif
  cout <<  "fPartMass = "         <<  fPartMass << endl;
  cout <<  "fPcentral = "         <<  fPcentral << " " <<fPCentralOffset << endl;
  cout <<  "fThetalab = "         <<  fTheta_lab << " " <<fThetaCentralOffset << endl;
  fPcentral= fPcentral*(1.+fPCentralOffset/100.);
  // Check that these offsets are in radians
  fTheta_lab=fTheta_lab + fThetaCentralOffset*TMath::RadToDeg();
  Double_t ph = fPhi_lab+fPhiOffset*TMath::RadToDeg();
  cout << "Central angles = " << fTheta_lab << endl;
  // SetCentralAngles method in podd THaSpectrometer
  // fTheta_lab and ph are geographical angles, converts to spherical coordinates
  // Need to set fTheta_lab to negative for spectrometer like HMS on beam right
  //   This gives phi_sph = 180 th_sph=abs(th_geo)
  // Computes TRotation fToLabRot and fToTraRot
  Bool_t bend_down = kFALSE;
  SetCentralAngles(fTheta_lab, ph, bend_down);
  Double_t off_x = 0.0, off_y = 0.0, off_z = 0.0;
  fPointingOffset.SetXYZ( off_x, off_y, off_z );

  //
  ifstream ifile;
  ifile.open(reconCoeffFilename.c_str());
  if(!ifile.is_open()) {
    Error(here, "error opening reconstruction coefficient file %s",reconCoeffFilename.c_str());
    //    return kInitError; // Is this the right return code?
    return kOK;
  }

  string line="!";
  int good=1;
  while(good && line[0]=='!') {
    good = getline(ifile,line).good();
    //    cout << line << endl;
  }
  // Read in focal plane rotation coefficients
  // Probably not used, so for now, just paste in fortran code as a comment
  while(good && line.compare(0,4," ---")!=0) {
    //  if(line(1:13).eq.'h_ang_slope_x')read(line,1201,err=94)h_ang_slope_x
    //  if(line(1:13).eq.'h_ang_slope_y')read(line,1201,err=94)h_ang_slope_y
    //  if(line(1:14).eq.'h_ang_offset_x')read(line,1201,err=94)h_ang_offset_x
    //  if(line(1:14).eq.'h_ang_offset_y')read(line,1201,err=94)h_ang_offset_y
    //  if(line(1:14).eq.'h_det_offset_x')read(line,1201,err=94)h_det_offset_x
    //  if(line(1:14).eq.'h_det_offset_y')read(line,1201,err=94)h_det_offset_y
    //  if(line(1:14).eq.'h_z_true_focus')read(line,1201,err=94)h_z_true_focus
    good = getline(ifile,line).good();
  }
  // Read in reconstruction coefficients and exponents
  line=" ";
  good = getline(ifile,line).good();
  //  cout << line << endl;
  fNReconTerms = 0;
  fReconTerms.clear();
  fReconTerms.reserve(500);
  //cout << "Reading matrix elements" << endl;
  while(good && line.compare(0,4," ---")!=0) {
    fReconTerms.push_back(reconTerm());
    sscanf(line.c_str()," %le %le %le %le %1d%1d%1d%1d%1d"
	   ,&fReconTerms[fNReconTerms].Coeff[0],&fReconTerms[fNReconTerms].Coeff[1]
	   ,&fReconTerms[fNReconTerms].Coeff[2],&fReconTerms[fNReconTerms].Coeff[3]
	   ,&fReconTerms[fNReconTerms].Exp[0]
	   ,&fReconTerms[fNReconTerms].Exp[1]
	   ,&fReconTerms[fNReconTerms].Exp[2]
	   ,&fReconTerms[fNReconTerms].Exp[3]
	   ,&fReconTerms[fNReconTerms].Exp[4]);
    fNReconTerms++;
    good = getline(ifile,line).good();
  }
  cout << "Read " << fNReconTerms << " matrix element terms"  << endl;
  if(!good) {
    Error(here, "Error processing reconstruction coefficient file %s",reconCoeffFilename.c_str());
    return kInitError; // Is this the right return code?
  }
  return kOK;
}

//_____________________________________________________________________________
void THcHallCSpectrometer::EnforcePruneLimits()
{
  // Enforce minimum values for the prune cuts

  fPruneXp      = TMath::Max( 0.08, fPruneXp);
  fPruneYp      = TMath::Max( 0.04, fPruneYp);
  fPruneYtar    = TMath::Max( 4.0,  fPruneYtar);
  fPruneDelta   = TMath::Max( 13.0, fPruneDelta);
  fPruneBeta    = TMath::Max( 0.1,  fPruneBeta);
  fPruneDf      = TMath::Max( 1.0,  fPruneDf);
  fPruneChiBeta = TMath::Max( 2.0,  fPruneChiBeta);
  fPruneFpTime  = TMath::Max( 5.0,  fPruneFpTime);
  fPruneNPMT    = TMath::Max( 6.0,  fPruneNPMT);
}

//_____________________________________________________________________________
Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
{
  // Reconstruct target coordinates for all tracks found in the focal plane.

  // In Hall A, this is passed off to the tracking detectors.
  // In Hall C, we do the target traceback here since the traceback should
  // not depend on which tracking detectors are used.

  fNtracks = tracks.GetLast()+1;

  for (Int_t it=0;it<tracks.GetLast()+1;it++) {
    THaTrack* track = static_cast<THaTrack*>( tracks[it] );
    Double_t xptar=kBig,yptar=kBig,ytar=kBig,delta=kBig;
    Double_t xtar=0;
    CalculateTargetQuantities(track,xtar,xptar,ytar,yptar,delta); 
    // Transfer results to track
    // No beam raster yet
    //; In transport coordinates phi = hyptar = dy/dz and theta = hxptar = dx/dz
    //;    but for unknown reasons the yp offset is named  htheta_offset
    //;    and  the xp offset is named  hphi_offset

    track->SetTarget(0.0, ytar*100.0, xptar+fPhiOffset, yptar+fThetaOffset);
    track->SetDp(delta*100.0+fDeltaOffset);	// Percent.  (Don't think podd cares if it is % or fraction)
    // There is an hpcentral_offset that needs to be applied somewhere.
    // (happly_offs)
    Double_t ptemp = fPcentral*(1+track->GetDp()/100.0);
    track->SetMomentum(ptemp);
    TVector3 pvect_temp;
    TransportToLab(track->GetP(),track->GetTTheta(),track->GetTPhi(),pvect_temp);
    track->SetPvect(pvect_temp);
  }

  if (fHodo==0 || (( fSelUsingScin == 0 ) && ( fSelUsingPrune == 0 )) ) {
    BestTrackSimple();
  } else if (fHodo!=0 && fSelUsingPrune !=0) {
    BestTrackUsingPrune();
  } else if (fHodo!=0){
    BestTrackUsingScin();
  }


  return 0;
}
//
void THcHallCSpectrometer::CalculateTargetQuantities(THaTrack* track,Double_t gbeam_y,Double_t  xptar,Double_t ytar,Double_t yptar,Double_t delta) 
{
    Double_t hut[5];
    Double_t hut_rot[5];

    hut[0] = track->GetX()/100.0 + fZTrueFocus*track->GetTheta() + fDetOffset_x;//m
    hut[1] = track->GetTheta() + fAngOffset_x;//radians
    hut[2] = track->GetY()/100.0 + fZTrueFocus*track->GetPhi() + fDetOffset_y;//m
    hut[3] = track->GetPhi() + fAngOffset_y;//radians

    hut[4] = -gbeam_y/100.0;

    // Retrieve the focal plane coordnates
    // Do the transpormation
    // Stuff results into track
    hut_rot[0] = hut[0];
    hut_rot[1] = hut[1] + hut[0]*fAngSlope_x;
    hut_rot[2] = hut[2];
    hut_rot[3] = hut[3] + hut[2]*fAngSlope_y;
    hut_rot[4] = hut[4];

    // Compute COSY sums
    Double_t sum[4];
    for(Int_t k=0;k<4;k++) {
      sum[k] = 0.0;
    }
    for(Int_t iterm=0;iterm<fNReconTerms;iterm++) {
      Double_t term=1.0;
      for(Int_t j=0;j<5;j++) {
	if(fReconTerms[iterm].Exp[j]!=0) {
	  term *= pow(hut_rot[j],fReconTerms[iterm].Exp[j]);
	}
      }
      for(Int_t k=0;k<4;k++) {
	sum[k] += term*fReconTerms[iterm].Coeff[k];
      }
    }
    xptar=sum[0];
    ytar=sum[1];
    yptar=sum[2];
    delta=sum[3];
}
//
//_____________________________________________________________________________
Int_t THcHallCSpectrometer::TrackCalc()
{

  if( fNtracks > 0 ) {
    Int_t hit_gold_track=0; // find track with index =0 which is best track
    for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){
      THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
      if (aTrack->GetIndex()==0) hit_gold_track=itrack;  
    }
    
    fDC->SetFocalPlaneBestTrack(hit_gold_track);
    fGoldenTrack = static_cast<THaTrack*>( fTracks->At(hit_gold_track) );
    fTrkIfo      = *fGoldenTrack;
    fTrk         = fGoldenTrack;
  } else
    fGoldenTrack = NULL;


  return TrackTimes( fTracks );
}

//_____________________________________________________________________________
Int_t THcHallCSpectrometer::BestTrackSimple()
{

  if( GetTrSorting() )   fTracks->Sort();

  // Assign index=0 to the best track, 
    for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){
      THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
      aTrack->SetIndex(1);  
      if (itrack==0) aTrack->SetIndex(0);  
    }

  return(0);
}

//_____________________________________________________________________________
Int_t THcHallCSpectrometer::BestTrackUsingScin()
{
  Double_t chi2Min;

  if( fNtracks > 0 ) {
    fGoodTrack = -1;
    chi2Min = 10000000000.0;
    Int_t y2Dmin = 100;
    Int_t x2Dmin = 100;

    Bool_t* x2Hits        = new Bool_t [fHodo->GetNPaddles(2)];
    Bool_t* y2Hits        = new Bool_t [fHodo->GetNPaddles(3)];
    for (UInt_t j = 0; j < fHodo->GetNPaddles(2); j++ ){
      x2Hits[j] = kFALSE;
    }
    for (UInt_t j = 0; j < fHodo->GetNPaddles(3); j++ ){
      y2Hits[j] = kFALSE;
    }

    for (Int_t rawindex=0; rawindex<fHodo->GetTotHits(); rawindex++) {
      Int_t ip = fHodo->GetGoodRawPlane(rawindex);
      if(ip==2) {	// X2
	Int_t goodRawPad = fHodo->GetGoodRawPad(rawindex);
	x2Hits[goodRawPad] = kTRUE;
      } else if (ip==3) { // Y2
	Int_t goodRawPad = fHodo->GetGoodRawPad(rawindex);
	y2Hits[goodRawPad] = kTRUE;
      }
    }

    for (Int_t itrack = 0; itrack < fNtracks; itrack++ ){
      Double_t chi2PerDeg;

      THaTrack* aTrack = static_cast<THaTrack*>( fTracks->At(itrack) );
      if (!aTrack) return -1;

      if ( aTrack->GetNDoF() > fSelNDegreesMin ){
	chi2PerDeg =  aTrack->GetChi2() / aTrack->GetNDoF();

	if( ( aTrack->GetDedx()    > fSeldEdX1Min )   &&
	    ( aTrack->GetDedx()    < fSeldEdX1Max )   &&
	    ( aTrack->GetBeta()    > fSelBetaMin  )   &&
	    ( aTrack->GetBeta()    < fSelBetaMax  )   &&
	    ( aTrack->GetEnergy()  > fSelEtMin    )   &&
	    ( aTrack->GetEnergy()  < fSelEtMax    ) )
	  {
	    Int_t x2D, y2D;

	    if ( fNtracks > 1 ){
	      Double_t hitpos3  = aTrack->GetX() + aTrack->GetTheta() * ( fScin2XZpos + 0.5 * fScin2XdZpos );
	      Int_t icounter3  = TMath::Nint( ( hitpos3 - fHodo->GetPlaneCenter(2) ) / fHodo->GetPlaneSpacing(2) ) + 1;
	      Int_t hitCnt3  = TMath::Max( TMath::Min(icounter3, (Int_t) fHodo->GetNPaddles(2) ) , 1); // scin_2x_nr = 16
	      //	      fHitDist3 = fHitPos3 - ( fHodo->GetPlaneSpacing(2) * ( hitCnt3 - 1 ) + fHodo->GetPlaneCenter(2) );
	      Double_t hitpos4 = aTrack->GetY() + aTrack->GetPhi() * ( fScin2YZpos + 0.5 * fScin2YdZpos );
	      Int_t icounter4  = TMath::Nint( ( fHodo->GetPlaneCenter(3) - hitpos4 ) / fHodo->GetPlaneSpacing(3) ) + 1;
	      Int_t hitCnt4  = TMath::Max( TMath::Min(icounter4, (Int_t) fHodo->GetNPaddles(3) ) , 1); // scin_2y_nr = 10
	      //	      fHitDist4 = fHitPos4 - ( fHodo->GetPlaneCenter(3) - fHodo->GetPlaneSpacing(3) * ( hitCnt4 - 1 ) );
	      // Plane 3
	      Int_t mindiff=1000;
	      for (UInt_t i = 0; i <  fHodo->GetNPaddles(2); i++ ){
		if ( x2Hits[i] ) {
		  Int_t diff = TMath::Abs((Int_t)hitCnt3-(Int_t)i-1);
		  if (diff < mindiff) mindiff = diff;
		}
	      }
	      if(mindiff < 1000) {
		x2D = mindiff;
	      } else {
		x2D = 0;	// Is this what we really want if there were no hits on this plane?
	      }

	      // Plane 4
	      mindiff = 1000;
	      for (UInt_t i = 0; i < fHodo->GetNPaddles(3); i++ ){
		if ( y2Hits[i] ) {
		  Int_t diff = TMath::Abs((Int_t)hitCnt4-(Int_t)i-1);
		  if (diff < mindiff) mindiff = diff;
		}
	      }
	      if(mindiff < 1000) {
		y2D = mindiff;
	      } else {
		y2D = 0;	// Is this what we really want if there were no hits on this plane?
	      }
	    } else { // Only a single track
	      x2D = 0.;
	      y2D = 0.;
	    }

	    if ( y2D <= y2Dmin ) {
	      if ( y2D < y2Dmin ) {
		x2Dmin = 100;
		chi2Min = 10000000000.;
	      } // y2D min

	      if ( x2D <= x2Dmin ){
		if ( x2D < x2Dmin ){
		  chi2Min = 10000000000.0;
		} // condition x2D
		if ( chi2PerDeg < chi2Min ){

		  fGoodTrack = itrack; // fGoodTrack = itrack
		  y2Dmin = y2D;
		  x2Dmin = x2D;
		  chi2Min = chi2PerDeg;

		  fGoldenTrack = static_cast<THaTrack*>( fTracks->At( fGoodTrack ) );
		  fTrkIfo      = *fGoldenTrack;
		  fTrk         = fGoldenTrack;
		}
	      } // condition x2D
	    } // condition y2D
	  } // conditions for dedx, beta and trac energy
      } // confition for fNFreeFP greater than fSelNDegreesMin
    } // loop over tracks

    // Fall back to using track with minimum chi2
    if ( fGoodTrack == -1 ){

      chi2Min = 10000000000.0;
      for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
	Double_t chi2PerDeg;
	THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
	if (!aTrack) return -1;

	if ( aTrack->GetNDoF() > fSelNDegreesMin ){
	  chi2PerDeg =  aTrack->GetChi2() / aTrack->GetNDoF();

	  if ( chi2PerDeg < chi2Min ){

	    fGoodTrack = iitrack;
	    chi2Min = chi2PerDeg;

	  }
	}
      } // loop over trakcs
      // Set index for fGoodTrack = 0
      for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
 	THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
        aTrack->SetIndex(1);
	if (iitrack==fGoodTrack) aTrack->SetIndex(0);
      }
     //
    }

  }

  return(0);
}

//_____________________________________________________________________________
Int_t THcHallCSpectrometer::BestTrackUsingPrune()
{


  Int_t nGood;
  Double_t chi2Min;

  if ( fNtracks > 0 ) {
    chi2Min   = 10000000000.0;
    fGoodTrack = 0;
    Bool_t* keep      = new Bool_t [fNtracks];
    Int_t* reject    = new Int_t  [fNtracks];

    THaTrack *testTracks[fNtracks];

    // ! Initialize all tracks to be good
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      keep[ptrack] = kTRUE;
      reject[ptrack] = 0;
      testTracks[ptrack] = static_cast<THaTrack*>( fTracks->At(ptrack) );
      if (!testTracks[ptrack]) return -1;
    }

    // ! Prune on xptar
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( TMath::Abs( testTracks[ptrack]->GetTTheta() ) < fPruneXp ) && ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if ( nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( TMath::Abs( testTracks[ptrack]->GetTTheta() ) >= fPruneXp ){
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 1;
	}
      }
    }

    // ! Prune on yptar
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( TMath::Abs( testTracks[ptrack]->GetTPhi() ) < fPruneYp ) && ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( TMath::Abs( testTracks[ptrack]->GetTPhi() ) >= fPruneYp ){
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 2;

	}
      }
    }

    // !     Prune on ytar
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( TMath::Abs( testTracks[ptrack]->GetTY() ) < fPruneYtar ) && ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( TMath::Abs( testTracks[ptrack]->GetTY() ) >= fPruneYtar ){
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 10;
	}
      }
    }

    // !     Prune on delta
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( TMath::Abs( testTracks[ptrack]->GetDp() ) < fPruneDelta ) && ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( TMath::Abs( testTracks[ptrack]->GetDp() ) >= fPruneDelta ){
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 20;
	}
      }
    }

    // !     Prune on beta
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      Double_t p = testTracks[ptrack]->GetP();
      Double_t betaP = p / TMath::Sqrt( p * p + fPartMass * fPartMass );
      if ( ( TMath::Abs( testTracks[ptrack]->GetBeta() - betaP ) < fPruneBeta ) && ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	Double_t p = testTracks[ptrack]->GetP();
	Double_t betaP = p / TMath::Sqrt( p * p + fPartMass * fPartMass );
	if ( TMath::Abs( testTracks[ptrack]->GetBeta() - betaP ) >= fPruneBeta ) {
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 100;
	}
      }
    }

    // !     Prune on deg. freedom for track chisq
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( testTracks[ptrack]->GetNDoF() >= fPruneDf ) && ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( testTracks[ptrack]->GetNDoF() < fPruneDf ){
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 200;
	}
      }
    }

    //!     Prune on num pmt hits
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( testTracks[ptrack]->GetNPMT() >= fPruneNPMT ) && ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( testTracks[ptrack]->GetNPMT() < fPruneNPMT ){
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 100000;
	}
      }
    }

    // !     Prune on beta chisqr
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( testTracks[ptrack]->GetBetaChi2() < fPruneChiBeta ) &&
	   ( testTracks[ptrack]->GetBetaChi2() > 0.01 )          &&   ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( ( testTracks[ptrack]->GetBetaChi2() >= fPruneChiBeta ) ||
	     ( testTracks[ptrack]->GetBetaChi2() <= 0.01          ) ){
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 1000;
	}
      }
    }

    // !     Prune on fptime
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( TMath::Abs( testTracks[ptrack]->GetFPTime() - fHodo->GetStartTimeCenter() ) < fPruneFpTime )  &&
	   ( keep[ptrack] ) ){
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( TMath::Abs( testTracks[ptrack]->GetFPTime() - fHodo->GetStartTimeCenter() ) >= fPruneFpTime ) {
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 2000;
	}
      }
    }

    // !     Prune on Y2 being hit
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( testTracks[ptrack]->GetGoodPlane4() == 1 )  &&  keep[ptrack] ) {
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( testTracks[ptrack]->GetGoodPlane4() != 1 ) {
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 10000;
	}
      }
    }

    // !     Prune on X2 being hit
    nGood = 0;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
      if ( ( testTracks[ptrack]->GetGoodPlane3() == 1 )  &&  keep[ptrack] ) {
	nGood ++;
      }
    }
    if (nGood > 0 ) {
      for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){
	if ( testTracks[ptrack]->GetGoodPlane3() != 1 ) {
	  keep[ptrack] = kFALSE;
	  reject[ptrack] = reject[ptrack] + 20000;
	}
      }
    }

    // !     Pick track with best chisq if more than one track passed prune tests
    Double_t chi2PerDeg = 0.;
    for (Int_t ptrack = 0; ptrack < fNtracks; ptrack++ ){

      chi2PerDeg =  testTracks[ptrack]->GetChi2() / testTracks[ptrack]->GetNDoF();

      if ( ( chi2PerDeg < chi2Min ) && ( keep[ptrack] ) ){
	fGoodTrack = ptrack;
	chi2Min = chi2PerDeg;
      }
    }
      // Set index=0 for fGoodTrack 
      for (Int_t iitrack = 0; iitrack < fNtracks; iitrack++ ){
 	THaTrack* aTrack = dynamic_cast<THaTrack*>( fTracks->At(iitrack) );
        aTrack->SetIndex(1);
	if (iitrack==fGoodTrack) aTrack->SetIndex(0);
      }
     //


  } 

  return(0);
}

//_____________________________________________________________________________
Int_t THcHallCSpectrometer::TrackTimes( TClonesArray* Tracks ) {
  // Do the actual track-timing (beta) calculation.
  // Use multiple scintillators to average together and get "best" time at S1.
  //
  // To be useful, a meaningful timing resolution should be assigned
  // to each Scintillator object (part of the database).


  return 0;
}

//_____________________________________________________________________________
Int_t THcHallCSpectrometer::ReadRunDatabase( const TDatime& date )
{
  // Override THaSpectrometer with nop method.  All needed kinamatics
  // read in ReadDatabase.

  return kOK;
}

//_____________________________________________________________________________
ClassImp(THcHallCSpectrometer)