THcHodoscope.cxx

/** \class THcHodoscope
    \ingroup Detectors

\brief Generic hodoscope consisting of multiple
planes with multiple paddles with phototubes on both ends.

This differs from Hall A scintillator class in that it is the whole
hodoscope array, not just one plane.

*/

#include "THcSignalHit.h"
#include "THcShower.h"
#include "THcCherenkov.h"
#include "THcHallCSpectrometer.h"
#include "THcHitList.h"
#include "THcRawShowerHit.h"
#include "TClass.h"
#include "math.h"
#include "THaSubDetector.h"

#include "THcHodoscope.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcDetectorMap.h"
#include "THaGlobals.h"
#include "THaCutList.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "VarDef.h"
#include "VarType.h"
#include "THaTrack.h"
#include "TClonesArray.h"
#include "TMath.h"

#include "THaTrackProj.h"
#include <vector>

#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cassert>

using namespace std;
using std::vector;

//_____________________________________________________________________________
THcHodoscope::THcHodoscope( const char* name, const char* description,
				  THaApparatus* apparatus ) :
  THaNonTrackingDetector(name,description,apparatus)
{
  // Constructor

  //fTrackProj = new TClonesArray( "THaTrackProj", 5 );
  // Construct the planes
  fNPlanes = 0;			// No planes until we make them
  fStartTime=-1e5;
  fGoodStartTime=kFALSE;
}

//_____________________________________________________________________________
THcHodoscope::THcHodoscope( ) :
  THaNonTrackingDetector()
{
  // Constructor
}

//_____________________________________________________________________________
void THcHodoscope::Setup(const char* name, const char* description)
{
  /**
     Create the scintillator plane objects for the hodoscope.
     
     Uses the Xhodo_num_planes and Xhodo_plane_names to get the number of
     planes and their names.

     Gets a pointer to the Cherenkov named "cer" ("hgcer" in the case of the SHMS.)

  */     
  if( IsZombie()) return;

  // fDebug = 1;  // Keep this at one while we're working on the code

  char prefix[2];

  prefix[0]=tolower(GetApparatus()->GetName()[0]);
  prefix[1]='\0';

  TString temp(prefix[0]);
  fSHMS=kFALSE;
  if (temp == "p" ) fSHMS=kTRUE;
  TString histname=temp+"_timehist";
  hTime = new TH1F(histname,"",400,0,200);
  // cout << " fSHMS = " << fSHMS << endl;
  string planenamelist;
  DBRequest listextra[]={
    {"hodo_num_planes", &fNPlanes, kInt},
    {"hodo_plane_names",&planenamelist, kString},
    {"hodo_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
    {"hodo_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
    {0}
  };
  //fNPlanes = 4; 		// Default if not defined
  fTDC_RefTimeCut = 0;		// Minimum allowed reference times
  fADC_RefTimeCut = 0;
  gHcParms->LoadParmValues((DBRequest*)&listextra,prefix);

  cout << "Plane Name List : " << planenamelist << endl;

  vector<string> plane_names = vsplit(planenamelist);
  // Plane names
  if(plane_names.size() != (UInt_t) fNPlanes) {
    cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
    // Should quit.  Is there an official way to quit?
  }
  fPlaneNames = new char* [fNPlanes];
  for(Int_t i=0;i<fNPlanes;i++) {
    fPlaneNames[i] = new char[plane_names[i].length()+1];
    strcpy(fPlaneNames[i], plane_names[i].c_str());
  }

  // Probably shouldn't assume that description is defined
  char* desc = new char[strlen(description)+100];
  fPlanes = new THcScintillatorPlane* [fNPlanes];
  for(Int_t i=0;i < fNPlanes;i++) {
    strcpy(desc, description);
    strcat(desc, " Plane ");
    strcat(desc, fPlaneNames[i]);
    fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i+1, this); // Number planes starting from zero!!
    cout << "Created Scintillator Plane " << fPlaneNames[i] << ", " << desc << endl;
  }

  // Save the nominal particle mass
  THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
  fPartMass = app->GetParticleMass();
  fBetaNominal = app->GetBetaAtPcentral();



  if (fSHMS) {
    fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("hgcer"));
  } else {
    fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("cer"));
  }

  delete [] desc;
}

//_____________________________________________________________________________
THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
{
  // cout << "In THcHodoscope::Init()" << endl;
  Setup(GetName(), GetTitle());

  char EngineDID[] = "xSCIN";
  EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
  if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
    static const char* const here = "Init()";
    Error( Here(here), "Error filling detectormap for %s.", EngineDID );
    return kInitError;
  }

  // Should probably put this in ReadDatabase as we will know the
  // maximum number of hits after setting up the detector map
  // But it needs to happen before the sub detectors are initialized
  // so that they can get the pointer to the hitlist.
  cout << " Hodo tdc ref time cut = " << fTDC_RefTimeCut << " " << fADC_RefTimeCut << endl;

  InitHitList(fDetMap, "THcRawHodoHit", fDetMap->GetTotNumChan()+1,
	      fTDC_RefTimeCut, fADC_RefTimeCut);

  EStatus status;
  // This triggers call of ReadDatabase and DefineVariables
  if( (status = THaNonTrackingDetector::Init( date )) )
    return fStatus=status;

  for(Int_t ip=0;ip<fNPlanes;ip++) {
    if((status = fPlanes[ip]->Init( date ))) {
      return fStatus=status;
    }
  }

  fNScinHits     = new Int_t [fNPlanes];
  fGoodPlaneTime = new Bool_t [fNPlanes];
  fNPlaneTime    = new Int_t [fNPlanes];
  fSumPlaneTime  = new Double_t [fNPlanes];

  //  Double_t  fHitCnt4 = 0., fHitCnt3 = 0.;

  // Int_t m = 0;
  // fScinHit = new Double_t*[fNPlanes];
  // for ( m = 0; m < fNPlanes; m++ ){
  //   fScinHit[m] = new Double_t[fNPaddle[0]];
  // }

  for (int ip=0; ip<fNPlanes; ++ip) {
    fScinHitPaddle.push_back(std::vector<Int_t>(fNPaddle[ip], 0));
  }

  fPresentP = 0;
  THaVar* vpresent = gHaVars->Find(Form("%s.present",GetApparatus()->GetName()));
  if(vpresent) {
    fPresentP = (Bool_t *) vpresent->GetValuePointer();
  }

  return fStatus = kOK;
}
//_____________________________________________________________________________
Int_t THcHodoscope::ReadDatabase( const TDatime& date )
{
  /**
     Read this detector's parameters from the ThcParmlist.

     This function is called by THaDetectorBase::Init() once at the
     beginning of the analysis.
  */
  //  static const char* const here = "ReadDatabase()";
  char prefix[2];
  char parname[100];

  // Determine which spectrometer in order to construct
  // the parameter names (e.g. hscin_1x_nr vs. sscin_1x_nr)

  prefix[0]=tolower(GetApparatus()->GetName()[0]);
  //
  prefix[1]='\0';
  strcpy(parname,prefix);
  strcat(parname,"scin_");
  //  Int_t plen=strlen(parname);
  // cout << " readdatabse hodo fnplanes = " << fNPlanes << endl;

  CreateMissReportParms(Form("%sscin",prefix));

  fBetaNoTrk = 0.;
  fBetaNoTrkChiSq = 0.;

  fNPaddle = new UInt_t [fNPlanes];
  fFPTime = new Double_t [fNPlanes];
  fPlaneCenter = new Double_t[fNPlanes];
  fPlaneSpacing = new Double_t[fNPlanes];

  prefix[0]=tolower(GetApparatus()->GetName()[0]);
  //
  prefix[1]='\0';

  for(Int_t i=0;i<fNPlanes;i++) {

    DBRequest list[]={
      {Form("scin_%s_nr",fPlaneNames[i]), &fNPaddle[i], kInt},
      {0}
    };


    gHcParms->LoadParmValues((DBRequest*)&list,prefix);
  }

  // GN added
  // reading variables from *hodo.param
  fMaxScinPerPlane=fNPaddle[0];
  for (Int_t i=1;i<fNPlanes;i++) {
    fMaxScinPerPlane=(fMaxScinPerPlane > fNPaddle[i])? fMaxScinPerPlane : fNPaddle[i];
  }
  // need this for "padded arrays" i.e. 4x16 lists of parameters (GN)
  fMaxHodoScin=fMaxScinPerPlane*fNPlanes;
  if (fDebug>=1)  cout <<"fMaxScinPerPlane = "<<fMaxScinPerPlane<<" fMaxHodoScin = "<<fMaxHodoScin<<endl;

  fHodoVelLight=new Double_t [fMaxHodoScin];
  fHodoPosSigma=new Double_t [fMaxHodoScin];
  fHodoNegSigma=new Double_t [fMaxHodoScin];
  fHodoPosMinPh=new Double_t [fMaxHodoScin];
  fHodoNegMinPh=new Double_t [fMaxHodoScin];
  fHodoPosPhcCoeff=new Double_t [fMaxHodoScin];
  fHodoNegPhcCoeff=new Double_t [fMaxHodoScin];
  fHodoPosTimeOffset=new Double_t [fMaxHodoScin];
  fHodoNegTimeOffset=new Double_t [fMaxHodoScin];
  fHodoPosPedLimit=new Int_t [fMaxHodoScin];
  fHodoNegPedLimit=new Int_t [fMaxHodoScin];
  fHodoPosInvAdcOffset=new Double_t [fMaxHodoScin];
  fHodoNegInvAdcOffset=new Double_t [fMaxHodoScin];
  fHodoPosInvAdcLinear=new Double_t [fMaxHodoScin];
  fHodoNegInvAdcLinear=new Double_t [fMaxHodoScin];
  fHodoPosInvAdcAdc=new Double_t [fMaxHodoScin];
  fHodoNegInvAdcAdc=new Double_t [fMaxHodoScin];

  //New Time-Walk Calibration Parameters
  fHodoVelFit=new Double_t [fMaxHodoScin];
  fHodoCableFit=new Double_t [fMaxHodoScin];
  fHodo_LCoeff=new Double_t [fMaxHodoScin];
  fHodoPos_c1=new Double_t [fMaxHodoScin];
  fHodoNeg_c1=new Double_t [fMaxHodoScin];
  fHodoPos_c2=new Double_t [fMaxHodoScin];
  fHodoNeg_c2=new Double_t [fMaxHodoScin];
  fHodoSigmaPos=new Double_t [fMaxHodoScin];
  fHodoSigmaNeg=new Double_t [fMaxHodoScin];

  fNHodoscopes = 2;
  fxLoScin = new Int_t [fNHodoscopes];
  fxHiScin = new Int_t [fNHodoscopes];
  fyLoScin = new Int_t [fNHodoscopes];
  fyHiScin = new Int_t [fNHodoscopes];
  fHodoSlop = new Double_t [fNPlanes];
  fTdcOffset = new Int_t [fNPlanes];
  fAdcTdcOffset = new Double_t [fNPlanes];
  fHodoPosAdcTimeWindowMin = new Double_t [fMaxHodoScin];
  fHodoPosAdcTimeWindowMax = new Double_t [fMaxHodoScin];
  fHodoNegAdcTimeWindowMin = new Double_t [fMaxHodoScin];
  fHodoNegAdcTimeWindowMax = new Double_t [fMaxHodoScin];


  for(Int_t ip=0;ip<fNPlanes;ip++) { // Set a large default window
   fTdcOffset[ip] = 0 ;
   fAdcTdcOffset[ip] = 0.0 ;
  }

  DBRequest list[]={
    {"cosmicflag",                       &fCosmicFlag,            kInt,            0,  1},
    {"NumPlanesBetaCalc",                       &fNumPlanesBetaCalc,            kInt,            0,  1},
    {"start_time_center",                &fStartTimeCenter,                      kDouble},
    {"start_time_slop",                  &fStartTimeSlop,                        kDouble},
    {"scin_tdc_to_time",                 &fScinTdcToTime,                        kDouble},
    {"scin_tdc_min",                     &fScinTdcMin,                           kDouble},
    {"scin_tdc_max",                     &fScinTdcMax,                           kDouble},
    {"tof_tolerance",                    &fTofTolerance,          kDouble,         0,  1},
    {"pathlength_central",               &fPathLengthCentral,                    kDouble},
    {"hodo_pos_sigma",                   &fHodoPosSigma[0],       kDouble,  fMaxHodoScin, 1},
    {"hodo_neg_sigma",                   &fHodoNegSigma[0],       kDouble,  fMaxHodoScin, 1},
    {"hodo_pos_ped_limit",               &fHodoPosPedLimit[0],    kInt,     fMaxHodoScin, 1},
    {"hodo_neg_ped_limit",               &fHodoNegPedLimit[0],    kInt,     fMaxHodoScin, 1},
    {"tofusinginvadc",                   &fTofUsingInvAdc,        kInt,            0,  1},
    {"xloscin",                          &fxLoScin[0],            kInt,     (UInt_t) fNHodoscopes},
    {"xhiscin",                          &fxHiScin[0],            kInt,     (UInt_t) fNHodoscopes},
    {"yloscin",                          &fyLoScin[0],            kInt,     (UInt_t) fNHodoscopes},
    {"yhiscin",                          &fyHiScin[0],            kInt,     (UInt_t) fNHodoscopes},
    {"track_eff_test_num_scin_planes",   &fTrackEffTestNScinPlanes,                 kInt},
    {"cer_npe",                          &fNCerNPE,               kDouble,         0,  1},
    {"normalized_energy_tot",            &fNormETot,              kDouble,         0,  1},
    {"hodo_slop",                        fHodoSlop,               kDouble,  (UInt_t) fNPlanes},
    {"debugprintscinraw",                &fdebugprintscinraw,               kInt,  0,1},
    {"hodo_tdc_offset",                  fTdcOffset,              kInt,     (UInt_t) fNPlanes, 1},
    {"hodo_adc_tdc_offset",              fAdcTdcOffset,           kDouble,  (UInt_t) fNPlanes, 1},
    {"hodo_PosAdcTimeWindowMin",         fHodoPosAdcTimeWindowMin, kDouble,  (UInt_t) fMaxHodoScin, 1},
    {"hodo_PosAdcTimeWindowMax",         fHodoPosAdcTimeWindowMax, kDouble,  (UInt_t) fMaxHodoScin, 1},
    {"hodo_NegAdcTimeWindowMin",         fHodoNegAdcTimeWindowMin, kDouble,  (UInt_t) fMaxHodoScin, 1},
    {"hodo_NegAdcTimeWindowMax",         fHodoNegAdcTimeWindowMax, kDouble,  (UInt_t) fMaxHodoScin, 1},
    {"dumptof",                          &fDumpTOF,               kInt,    0, 1},
    {"TOFCalib_shtrk_lo",                &fTOFCalib_shtrk_lo,               kDouble,    0, 1},
    {"TOFCalib_shtrk_hi",                &fTOFCalib_shtrk_hi,               kDouble,    0, 1},
    {"TOFCalib_cer_lo",                &fTOFCalib_cer_lo,               kDouble,    0, 1},
    {"TOFCalib_beta_lo",                &fTOFCalib_beta_lo,               kDouble,    0, 1},
    {"TOFCalib_beta_hi",                &fTOFCalib_beta_hi,               kDouble,    0, 1},
    {"dumptof_filename",                 &fTOFDumpFile,           kString, 0, 1},
    {0}
  };

  // Defaults if not defined in parameter file

  for(UInt_t ip=0;ip<fMaxHodoScin;ip++) {
    fHodoPosAdcTimeWindowMin[ip] = -1000.;
    fHodoPosAdcTimeWindowMax[ip] = 1000.;
    fHodoNegAdcTimeWindowMin[ip] = -1000.;
    fHodoNegAdcTimeWindowMax[ip] = 1000.;

    fHodoPosPedLimit[ip] = 0.0;
    fHodoNegPedLimit[ip] = 0.0;
    fHodoPosSigma[ip] = 0.2;
    fHodoNegSigma[ip] = 0.2;
  }
  fTOFCalib_shtrk_lo=-kBig;
  fTOFCalib_shtrk_hi= kBig;
  fTOFCalib_cer_lo=-kBig;
  fTOFCalib_beta_lo=-kBig;
  fTOFCalib_beta_hi= kBig;
  fdebugprintscinraw=0;
  fDumpTOF = 0;
  fTOFDumpFile="";
  fTofUsingInvAdc = 1;
  fTofTolerance = 3.0;
  fNCerNPE = 2.0;
  fNormETot = 0.7;
  fCosmicFlag=0;
  fNumPlanesBetaCalc=4;
  // Gets added to each reference time corrected raw TDC value
  // to make sure valid range is all positive.

  gHcParms->LoadParmValues((DBRequest*)&list,prefix);
  if (fCosmicFlag==1) cout << "Setup for cosmics in TOF"<< endl;
  // cout << " cosmic flag = " << fCosmicFlag << endl;
  if(fDumpTOF) {
    fDumpOut.open(fTOFDumpFile.c_str());
    if(fDumpOut.is_open()) {
      //fDumpOut << "Hodoscope Time of Flight calibration data" << endl;
    } else {
      fDumpTOF = 0;
      cout << "WARNING: Unable to open TOF Dump file " << fTOFDumpFile << endl;
      cout << "Data for TOF calibration not being written." << endl;
    }
  }

  // cout << " x1 lo = " << fxLoScin[0]
  //      << " x2 lo = " << fxLoScin[1]
  //      << " x1 hi = " << fxHiScin[0]
  //      << " x2 hi = " << fxHiScin[1]
  //      << endl;

  // cout << " y1 lo = " << fyLoScin[0]
  //      << " y2 lo = " << fyLoScin[1]
  //      << " y1 hi = " << fyHiScin[0]
  //      << " y2 hi = " << fyHiScin[1]
  //      << endl;

  // cout << "Hdososcope planes hits for trigger = " << fTrackEffTestNScinPlanes
  //      << " normalized energy min = " << fNormETot
  //      << " number of photo electrons = " << fNCerNPE
  //      << endl;

  //Set scin Velocity/Cable to default
  for (UInt_t i=0; i<fMaxHodoScin; i++) {
    fHodoVelLight[i] = 15.0;
  }

  if (fTofUsingInvAdc) {
    DBRequest list2[]={
      {"hodo_vel_light",                   &fHodoVelLight[0],       kDouble,  fMaxHodoScin, 1},
      {"hodo_pos_invadc_offset",&fHodoPosInvAdcOffset[0],kDouble,fMaxHodoScin},
      {"hodo_neg_invadc_offset",&fHodoNegInvAdcOffset[0],kDouble,fMaxHodoScin},
      {"hodo_pos_invadc_linear",&fHodoPosInvAdcLinear[0],kDouble,fMaxHodoScin},
      {"hodo_neg_invadc_linear",&fHodoNegInvAdcLinear[0],kDouble,fMaxHodoScin},
      {"hodo_pos_invadc_adc",&fHodoPosInvAdcAdc[0],kDouble,fMaxHodoScin},
      {"hodo_neg_invadc_adc",&fHodoNegInvAdcAdc[0],kDouble,fMaxHodoScin},
      {0}
    };
    
    gHcParms->LoadParmValues((DBRequest*)&list2,prefix);
  };
  /* if (!fTofUsingInvAdc) {
      DBRequest list3[]={
    {"hodo_vel_light",                   &fHodoVelLight[0],       kDouble,  fMaxHodoScin},
    {"hodo_pos_minph",                   &fHodoPosMinPh[0],       kDouble,  fMaxHodoScin},
    {"hodo_neg_minph",                   &fHodoNegMinPh[0],       kDouble,  fMaxHodoScin},
    {"hodo_pos_phc_coeff",               &fHodoPosPhcCoeff[0],    kDouble,  fMaxHodoScin},
    {"hodo_neg_phc_coeff",               &fHodoNegPhcCoeff[0],    kDouble,  fMaxHodoScin},
    {"hodo_pos_time_offset",             &fHodoPosTimeOffset[0],  kDouble,  fMaxHodoScin},
    {"hodo_neg_time_offset",             &fHodoNegTimeOffset[0],  kDouble,  fMaxHodoScin},
    {0}   
    };
  
    
    gHcParms->LoadParmValues((DBRequest*)&list3,prefix);
  */
     DBRequest list4[]={
    {"hodo_velFit",                      &fHodoVelFit[0],   kDouble,  fMaxHodoScin, 1},
    {"hodo_cableFit",                    &fHodoCableFit[0], kDouble,  fMaxHodoScin, 1},
    {"hodo_LCoeff",                      &fHodo_LCoeff[0],  kDouble,  fMaxHodoScin, 1},
    {"c1_Pos",                           &fHodoPos_c1[0],   kDouble,  fMaxHodoScin, 1},
    {"c1_Neg",                           &fHodoNeg_c1[0],   kDouble,  fMaxHodoScin, 1},
    {"c2_Pos",                           &fHodoPos_c2[0],   kDouble,  fMaxHodoScin, 1},
    {"c2_Neg",                           &fHodoNeg_c2[0],   kDouble,  fMaxHodoScin, 1},
    {"TDC_threshold",                    &fTdc_Thrs,        kDouble, 0, 1},
    {"hodo_PosSigma",                   &fHodoSigmaPos[0], kDouble,  fMaxHodoScin, 1},
    {"hodo_NegSigma",                   &fHodoSigmaNeg[0], kDouble,  fMaxHodoScin, 1},
    {0}   
     };
     
     fTdc_Thrs = 1.0;
     //Set Default Values if NOT defined in param file
     for (UInt_t i=0; i<fMaxHodoScin; i++)
       {
	 
	 //Turn OFF Time-Walk Correction if param file NOT found
	 fHodoPos_c1[i] = 0.0;
	 fHodoPos_c2[i] = 0.0;
	 fHodoNeg_c1[i] = 0.0;
	 fHodoNeg_c2[i] = 0.0;
       }
     for (UInt_t i=0; i<fMaxHodoScin; i++)                                                                    
       {  
	 //Set scin Velocity/Cable to default
	 fHodoCableFit[i] = 0.0;
	 fHodoVelFit[i] = 15.0;
	 //set time coeff between paddles to default
	 fHodo_LCoeff[i] = 0.0;
	 
       }
    
     gHcParms->LoadParmValues((DBRequest*)&list4,prefix);
  
  if (fDebug >=1) {
    cout <<"******* Testing Hodoscope Parameter Reading ***\n";
    cout<<"StarTimeCenter = "<<fStartTimeCenter<<endl;
    cout<<"StartTimeSlop = "<<fStartTimeSlop<<endl;
    cout <<"ScintTdcToTime = "<<fScinTdcToTime<<endl;
    cout <<"TdcMin = "<<fScinTdcMin<<" TdcMax = "<<fScinTdcMax<<endl;
    cout <<"TofTolerance = "<<fTofTolerance<<endl;
    cout <<"*** VelLight ***\n";
    for (Int_t i1=0;i1<fNPlanes;i1++) {
      cout<<"Plane "<<i1<<endl;
      for (UInt_t i2=0;i2<fMaxScinPerPlane;i2++) {
	cout<<fHodoVelLight[GetScinIndex(i1,i2)]<<" ";
      }
      cout <<endl;
    }
    cout <<endl<<endl;
    // check fHodoPosPhcCoeff
    /*
    cout <<"fHodoPosPhcCoeff = ";
    for (int i1=0;i1<fMaxHodoScin;i1++) {
      cout<<this->GetHodoPosPhcCoeff(i1)<<" ";
    }
    cout<<endl;
    */
  }
  //
  if ((fTofTolerance > 0.5) && (fTofTolerance < 10000.)) {
    cout << "USING "<<fTofTolerance<<" NSEC WINDOW FOR FP NO_TRACK CALCULATIONS.\n";
  }
  else {
    fTofTolerance= 3.0;
    cout << "*** USING DEFAULT 3 NSEC WINDOW FOR FP NO_TRACK CALCULATIONS!! ***\n";
  }
  fIsInit = true;
  return kOK;
}

//_____________________________________________________________________________
Int_t THcHodoscope::DefineVariables( EMode mode )
{
  /**
    Initialize global variables for histograms and Root tree
  */
  // cout << "THcHodoscope::DefineVariables called " << GetName() << endl;
  if( mode == kDefine && fIsSetup ) return kOK;
  fIsSetup = ( mode == kDefine );

  // Register variables in global list

  RVarDef vars[] = {
    // Move these into THcHallCSpectrometer using track fTracks
    {"beta",       "Beta including track info",                "fBeta"},
    {"betanotrack",       "Beta from scintillator hits",                "fBetaNoTrk"},
    {"betachisqnotrack",  "Chi square of beta from scintillator hits",  "fBetaNoTrkChiSq"},
    {"fpHitsTime",        "Time at focal plane from all hits",            "fFPTimeAll"},
    {"starttime",         "Hodoscope Start Time",                         "fStartTime"},
    {"goodstarttime",     "Hodoscope Good Start Time (logical flag)",                    "fGoodStartTime"},
    {"goodscinhit",       "Hit in fid area",                              "fGoodScinHits"},
    {"TimeHist_Sigma",       "",                              "fTimeHist_Sigma"},
    {"TimeHist_Peak",       "",                              "fTimeHist_Peak"},
    {"TimeHist_Hits",       "",                              "fTimeHist_Hits"},
     { 0 }
  };
  return DefineVarsFromList( vars, mode );
  //  return kOK;
}

//_____________________________________________________________________________
THcHodoscope::~THcHodoscope()
{
  // Destructor. Remove variables from global list.

  delete [] fFPTime;
  delete [] fPlaneCenter;
  delete [] fPlaneSpacing;

  if( fIsSetup )
    RemoveVariables();
  if( fIsInit )
    DeleteArrays();

  for( int i = 0; i < fNPlanes; ++i ) {
    delete fPlanes[i];
    delete [] fPlaneNames[i];
  }
  delete [] fPlanes;
  delete [] fPlaneNames;
}

//_____________________________________________________________________________
void THcHodoscope::DeleteArrays()
{
  // Delete member arrays. Used by destructor.
  // Int_t k;
  // for( k = 0; k < fNPlanes; k++){
  //   delete [] fScinHit[k];
  // }
  // delete [] fScinHit;

  delete [] fxLoScin;             fxLoScin = NULL;
  delete [] fxHiScin;             fxHiScin = NULL;
  delete [] fyLoScin;             fyLoScin = NULL;
  delete [] fyHiScin;             fyHiScin = NULL;
  delete [] fHodoSlop;            fHodoSlop = NULL;

  delete [] fNPaddle;             fNPaddle = NULL;
  delete [] fHodoVelLight;        fHodoVelLight = NULL;
  delete [] fHodoPosSigma;        fHodoPosSigma = NULL;
  delete [] fHodoNegSigma;        fHodoNegSigma = NULL;
  delete [] fHodoPosMinPh;        fHodoPosMinPh = NULL;
  delete [] fHodoNegMinPh;        fHodoNegMinPh = NULL;
  delete [] fHodoPosPhcCoeff;     fHodoPosPhcCoeff = NULL;
  delete [] fHodoNegPhcCoeff;     fHodoNegPhcCoeff = NULL;
  delete [] fHodoPosTimeOffset;   fHodoPosTimeOffset = NULL;
  delete [] fHodoNegTimeOffset;   fHodoNegTimeOffset = NULL;
  delete [] fHodoPosPedLimit;     fHodoPosPedLimit = NULL;
  delete [] fHodoNegPedLimit;     fHodoNegPedLimit = NULL;
  delete [] fHodoPosInvAdcOffset; fHodoPosInvAdcOffset = NULL;
  delete [] fHodoNegInvAdcOffset; fHodoNegInvAdcOffset = NULL;
  delete [] fHodoPosInvAdcLinear; fHodoPosInvAdcLinear = NULL;
  delete [] fHodoNegInvAdcLinear; fHodoNegInvAdcLinear = NULL;
  delete [] fHodoPosInvAdcAdc;    fHodoPosInvAdcAdc = NULL;
  delete [] fHodoNegInvAdcAdc;    fHodoNegInvAdcAdc = NULL;
  delete [] fGoodPlaneTime;       fGoodPlaneTime = NULL;
  delete [] fNPlaneTime;          fNPlaneTime = NULL;
  delete [] fSumPlaneTime;        fSumPlaneTime = NULL;
  delete [] fNScinHits;           fNScinHits = NULL;
  delete [] fTdcOffset;           fTdcOffset = NULL;
  delete [] fAdcTdcOffset;        fAdcTdcOffset = NULL;
  delete [] fHodoNegAdcTimeWindowMin;    fHodoNegAdcTimeWindowMin = NULL;
  delete [] fHodoNegAdcTimeWindowMax;    fHodoNegAdcTimeWindowMax = NULL;
  delete [] fHodoPosAdcTimeWindowMin;    fHodoPosAdcTimeWindowMin = NULL;
  delete [] fHodoPosAdcTimeWindowMax;    fHodoPosAdcTimeWindowMax = NULL;
   
  delete [] fHodoVelFit;                 fHodoVelFit = NULL;
  delete [] fHodoCableFit;               fHodoCableFit = NULL;
  delete [] fHodo_LCoeff;                fHodo_LCoeff = NULL;
  delete [] fHodoPos_c1;                 fHodoPos_c1 = NULL;
  delete [] fHodoNeg_c1;                 fHodoNeg_c1 = NULL;
  delete [] fHodoPos_c2;                 fHodoPos_c2 = NULL;
  delete [] fHodoNeg_c2;                 fHodoNeg_c2 = NULL;
  delete [] fHodoSigmaPos;               fHodoSigmaPos = NULL;
  delete [] fHodoSigmaNeg;               fHodoSigmaNeg = NULL;
}

//_____________________________________________________________________________
void THcHodoscope::Clear( Option_t* opt )
{
  /*! \brief Clears variables
   *
   *  Called by  THcHodoscope::Decode
   *
   */
  fTimeHist_Sigma=  kBig;
  fTimeHist_Peak=  kBig;
  fTimeHist_Hits=  kBig;

  fBeta = 0.0;
  fBetaNoTrk = 0.0;
  fBetaNoTrkChiSq = 0.0;
  fStartTime  = -1000.;
  fFPTimeAll= -1000.;
  fGoodStartTime = kFALSE;
  fGoodScinHits = 0;
 
  if( *opt != 'I' ) {
    for(Int_t ip=0;ip<fNPlanes;ip++) {
      fPlanes[ip]->Clear();
      fFPTime[ip]=0.;
      fPlaneCenter[ip]=0.;
      fPlaneSpacing[ip]=0.;
      for(UInt_t iPaddle=0;iPaddle<fNPaddle[ip]; ++iPaddle) {
	fScinHitPaddle[ip][iPaddle]=0;
      }
    }
  }
  fdEdX.clear();
  fNScinHit.clear();
  fNClust.clear();
  fClustSize.clear();
  fClustPos.clear();
  fThreeScin.clear();
  fGoodScinHitsX.clear();
  fGoodFlags.clear();
}

//_____________________________________________________________________________
Int_t THcHodoscope::Decode( const THaEvData& evdata )
{
  /*! \brief Decodes raw data and processes raw data into hits for each instance of  THcScintillatorPlane
   *
   *  - Reads raw data using THcHitList::DecodeToHitList
   *  - If one wants to subtract pedestals (assumed to be a set of data at beginning of run)
   *    + Must define "Pedestal_event" cut in the cuts definition file
   *    + For each "Pedestal_event" calls THcScintillatorPlane::AccumulatePedestals and returns
   *    + After First event which is not a  "Pedestal_event" calls THcScintillatorPlane::CalculatePedestals
   *  - For each scintillator plane THcScintillatorPlane::ProcessHits
   *  - Calls THcHodoscope::EstimateFocalPlaneTime
   *
   *
   */
  // Get the Hall C style hitlist (fRawHitList) for this event
  Bool_t present = kTRUE;	// Suppress reference time warnings
  if(fPresentP) {		// if this spectrometer not part of trigger
    present = *fPresentP;
  }
  fNHits = DecodeToHitList(evdata, !present);

  //
  // GN: print event number so we can cross-check with engine
  // if (evdata.GetEvNum()>1000)
  //   cout <<"\nhcana_event " << evdata.GetEvNum()<<endl;

  fCheckEvent = evdata.GetEvNum();
  fEventType =  evdata.GetEvType();

  if(gHaCuts->Result("Pedestal_event")) {
    Int_t nexthit = 0;
    for(Int_t ip=0;ip<fNPlanes;ip++) {

      nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
    }
    fAnalyzePedestals = 1;	// Analyze pedestals first normal events
    return(0);
  }
  if(fAnalyzePedestals) {
    for(Int_t ip=0;ip<fNPlanes;ip++) {

      fPlanes[ip]->CalculatePedestals();
    }
    fAnalyzePedestals = 0;	// Don't analyze pedestals next event
  }

  // Let each plane get its hits
  Int_t nexthit = 0;

  fNfptimes=0;
  Int_t thits = 0;
  for(Int_t ip=0;ip<fNPlanes;ip++) {

    fPlaneCenter[ip] = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
    fPlaneSpacing[ip] = fPlanes[ip]->GetSpacing();

    //    nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
    // GN: select only events that have reasonable TDC values to start with
    // as per the Engine h_strip_scin.f
    nexthit = fPlanes[ip]->ProcessHits(fRawHitList,nexthit);
    thits+=fPlanes[ip]->GetNScinHits();
  }
  fStartTime=-1000;
  if (thits>0 ) EstimateFocalPlaneTime();

  if (fdebugprintscinraw == 1) {
    for(UInt_t ihit = 0; ihit < fNRawHits ; ihit++) {
//    THcRawHodoHit* hit = (THcRawHodoHit *) fRawHitList->At(ihit);
//    cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
//	 << hit->fADC_pos << " " << hit->fADC_neg << " "  <<  hit->fTDC_pos
//	 << " " <<  hit->fTDC_neg << endl;
    }
    cout << endl;
  }


  return fNHits;
}

//_____________________________________________________________________________
void THcHodoscope::EstimateFocalPlaneTime()
{
  /*! \brief Calculates the Drift Chamber start time and fBetaNoTrk (velocity determined without track info)
   *
   *  - Called by  THcHodoscope::Decode
   *  - selects good scintillator paddle hits
   *     + loops through hits in each scintillator plane and fills histogram array, "timehist", with corrected times for positive
   *        and negative ends of each paddle
   *     + Determines the peak of "timehist"
   *
   */
  Int_t ihit=0;
  Int_t nscinhits=0;		// Total # hits with at least one good tdc
  hTime->Reset();
  //
  for(Int_t ip=0;ip<fNPlanes;ip++) {
    Int_t nphits=fPlanes[ip]->GetNScinHits();
    nscinhits += nphits;
    TClonesArray* hodoHits = fPlanes[ip]->GetHits();
    for(Int_t i=0;i<nphits;i++) {
      THcHodoHit *hit = (THcHodoHit*)hodoHits->At(i);
      if(hit->GetHasCorrectedTimes()) {
	Double_t postime=hit->GetPosTOFCorrectedTime();
	Double_t negtime=hit->GetNegTOFCorrectedTime();
	hTime->Fill(postime);
	hTime->Fill(negtime);
      }
    }
  }
  //
  //
  ihit = 0;
  Double_t fpTimeSum = 0.0;
  fNfptimes=0;
  Int_t  Ngood_hits_plane=0;
  Double_t Plane_fptime_sum=0.0;
  Bool_t goodplanetime[fNPlanes];
  Bool_t twogoodtimes[nscinhits];
  Double_t tmin = 0.5*hTime->GetMaximumBin();
  fTimeHist_Peak=  tmin;
  fTimeHist_Sigma=  hTime->GetRMS();
  fTimeHist_Hits=  hTime->Integral();
  for(Int_t ip=0;ip<fNumPlanesBetaCalc;ip++) {
    goodplanetime[ip] = kFALSE;
    Int_t nphits=fPlanes[ip]->GetNScinHits();
    TClonesArray* hodoHits = fPlanes[ip]->GetHits();
    Ngood_hits_plane=0;
    Plane_fptime_sum=0.0;
    for(Int_t i=0;i<nphits;i++) {
      THcHodoHit *hit = (THcHodoHit*)hodoHits->At(i);
      twogoodtimes[ihit] = kFALSE;
      if(hit->GetHasCorrectedTimes()) {
	Double_t postime=hit->GetPosTOFCorrectedTime();
	Double_t negtime=hit->GetNegTOFCorrectedTime();
	if ((postime>(tmin-fTofTolerance)) && (postime<(tmin+fTofTolerance)) &&
	    (negtime>(tmin-fTofTolerance)) && (negtime<(tmin+fTofTolerance)) ) {
	  hit->SetTwoGoodTimes(kTRUE);
	  twogoodtimes[ihit] = kTRUE;	// Both tubes fired
	  Int_t index=hit->GetPaddleNumber()-1;	 //
	  Double_t fptime;
	  if(fCosmicFlag==1) {
	    fptime = hit->GetScinCorrectedTime()
	      + (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
	      / (29.979 * fBetaNominal);
	  }else{
	    fptime = hit->GetScinCorrectedTime()
	      - (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
	      / (29.979 * fBetaNominal);
	  }
          Ngood_hits_plane++;
	  Plane_fptime_sum+=fptime;
	  fpTimeSum += fptime;
	  fNfptimes++;
	  goodplanetime[ip] = kTRUE;
	} else {
	  hit->SetTwoGoodTimes(kFALSE);
	}
      }
      ihit++;
    }
    if (Ngood_hits_plane) fPlanes[ip]->SetFpTime(Plane_fptime_sum/float(Ngood_hits_plane));
    fPlanes[ip]->SetNGoodHits(Ngood_hits_plane);
  }

  if(fNfptimes>0) {
    fStartTime = fpTimeSum/fNfptimes;
    fGoodStartTime=kTRUE;
    fFPTimeAll = fStartTime ;
  } else {
    fStartTime = fStartTimeCenter;
    fGoodStartTime=kFALSE;
    fFPTimeAll = fStartTime ;
  }
    //
   //
  hTime->Reset();
  //
  if((goodplanetime[0]||goodplanetime[1]) &&(goodplanetime[2]||goodplanetime[3])) {

    Double_t sumW = 0.;
    Double_t sumT = 0.;
    Double_t sumZ = 0.;
    Double_t sumZZ = 0.;
    Double_t sumTZ = 0.;
    Int_t ihhit = 0;

    for(Int_t ip=0;ip<fNumPlanesBetaCalc;ip++) {
      Int_t nphits=fPlanes[ip]->GetNScinHits();
      TClonesArray* hodoHits = fPlanes[ip]->GetHits();
      
      for(Int_t i=0;i<nphits;i++) {
	Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;

	if(twogoodtimes[ihhit]){

	  Double_t sigma = 0.0;
	  if(fTofUsingInvAdc)
	    {
	      sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
					   TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) );
	    }
	  else{
	    sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoSigmaPos[GetScinIndex(ip,index)],2) +
					 TMath::Power( fHodoSigmaNeg[GetScinIndex(ip,index)],2) ) );
	  }

	  Double_t scinWeight = 1 / TMath::Power(sigma,2);
	  Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
 	  //	  cout << "hit = " << ihhit + 1 << "   zpos = " << zPosition << "   sigma = " << sigma << endl;
	  //cout << "fHodoSigma+ = " << fHodoSigmaPos[GetScinIndex(ip,index)] << endl;
	  sumW  += scinWeight;
	  sumT  += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
	  sumZ  += scinWeight * zPosition;
	  sumZZ += scinWeight * ( zPosition * zPosition );
	  sumTZ += scinWeight * zPosition * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();

	} // condition of good scin time
	ihhit ++;
      } // loop over hits of plane
    } // loop over planes

    Double_t tmp = sumW * sumZZ - sumZ * sumZ ;
    Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ;
    Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;

    if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {

      fBetaNoTrk = tmp / tmpDenom;
      fBetaNoTrkChiSq = 0.;
      ihhit = 0;

      for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){                           // Loop over planes
	Int_t nphits=fPlanes[ip]->GetNScinHits();
	TClonesArray* hodoHits = fPlanes[ip]->GetHits();

	for(Int_t i=0;i<nphits;i++) {
	  Int_t index=((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber()-1;

	  if(twogoodtimes[ihhit]) {

	    Double_t zPosition = fPlanes[ip]->GetZpos() + (index%2)*fPlanes[ip]->GetDzpos();
	    Double_t timeDif = ( ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() - t0 );
	   
	    Double_t sigma = 0.0;
	    if(fTofUsingInvAdc){
	      sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoPosSigma[GetScinIndex(ip,index)],2) +
					   TMath::Power( fHodoNegSigma[GetScinIndex(ip,index)],2) ) );
	    }
	    else {
	      sigma = 0.5 * ( TMath::Sqrt( TMath::Power( fHodoSigmaPos[GetScinIndex(ip,index)],2) +
					   TMath::Power( fHodoSigmaNeg[GetScinIndex(ip,index)],2) ) );
	    }

	    fBetaNoTrkChiSq += ( ( zPosition / fBetaNoTrk - timeDif ) *
				 ( zPosition / fBetaNoTrk - timeDif ) ) / ( sigma * sigma );


	  } // condition for good scin time
	  ihhit++;
	} // loop over hits of a plane
      } // loop over planes

      Double_t pathNorm = 1.0;

      fBetaNoTrk = fBetaNoTrk * pathNorm;
      fBetaNoTrk = fBetaNoTrk / 29.979;    // velocity / c

    }  // condition for fTmpDenom
    else {
      fBetaNoTrk = 0.;
      fBetaNoTrkChiSq = -2.;
    } // else condition for fTmpDenom
    //
    fGoodEventTOFCalib=kFALSE;
    if ((fNumPlanesBetaCalc==4)&&goodplanetime[0]&&goodplanetime[1]&&goodplanetime[2]&&goodplanetime[3]&&fPlanes[0]->GetNGoodHits()==1&&fPlanes[1]->GetNGoodHits()==1&&fPlanes[2]->GetNGoodHits()==1&&fPlanes[3]->GetNGoodHits()==1) fGoodEventTOFCalib=kTRUE;
    if ((fNumPlanesBetaCalc==3)&&goodplanetime[0]&&goodplanetime[1]&&goodplanetime[2]&&fPlanes[0]->GetNGoodHits()==1&&fPlanes[1]->GetNGoodHits()==1&&fPlanes[2]->GetNGoodHits()==1) fGoodEventTOFCalib=kTRUE;
    //
    //
  }
}

//_____________________________________________________________________________
Int_t THcHodoscope::ApplyCorrections( void )
{
  return(0);
}
//_____________________________________________________________________________
Double_t THcHodoscope::TimeWalkCorrection(const Int_t& paddle,
					     const ESide side)
{
  return(0.0);
}


//_____________________________________________________________________________
Int_t THcHodoscope::CoarseProcess( TClonesArray& tracks )
{


  Int_t ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
  // -------------------------------------------------

  //  fDumpOut << " ntrack =  " << ntracks  << endl;

  if (ntracks > 0 ) {

    // **MAIN LOOP: Loop over all tracks and get corrected time, tof, beta...
    vector<Double_t> nPmtHit(ntracks);
    vector<Double_t> timeAtFP(ntracks);
    fdEdX.reserve(ntracks);
    fGoodFlags.reserve(ntracks);
    for ( Int_t itrack = 0; itrack < ntracks; itrack++ ) { // Line 133
      nPmtHit[itrack]=0;
      timeAtFP[itrack]=0;

      THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(itrack) );
      if (!theTrack) return -1;

      for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
	fGoodPlaneTime[ip] = kFALSE;
	fNScinHits[ip] = 0;
	fNPlaneTime[ip] = 0;
	fSumPlaneTime[ip] = 0.;
      }