THcHodoscope.cxx

      }
      std::vector<Double_t> dedx_temp;
      fdEdX.push_back(dedx_temp); // Create array of dedx per hit
      std::vector<std::vector<GoodFlags> > goodflagstmp1;
      fGoodFlags.push_back(goodflagstmp1);

      Int_t nFPTime = 0;
      Double_t betaChiSq = -3;
      Double_t beta = 0;
      //      timeAtFP[itrack] = 0.;
      Double_t sumFPTime = 0.; // Line 138
      fNScinHit.push_back(0);

      //! Calculate all corrected hit times and histogram
      //! This uses a copy of code below. Results are save in time_pos,neg
      //! including the z-pos. correction assuming nominal value of betap
      //! Code is currently hard-wired to look for a peak in the
      //! range of 0 to 100 nsec, with a group of times that all
      //! agree withing a time_tolerance of time_tolerance nsec. The normal
      //! peak position appears to be around 35 nsec.
      //! NOTE: if want to find particles with beta different than
      //! reference particle, need to make sure this is big enough
      //! to accomodate difference in TOF for other particles
      //! Default value in case user hasnt defined something reasonable


      // Loop over scintillator planes.
      // In ENGINE, its loop over good scintillator hits.
      hTime->Reset();
      fTOFCalc.clear();   // SAW - Can we
      fTOFPInfo.clear();  // SAW - combine these two?
      Int_t ihhit = 0;		// Hit # overall

      for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {

	std::vector<GoodFlags> goodflagstmp2;
	fGoodFlags[itrack].push_back(goodflagstmp2);

	fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
	TClonesArray* hodoHits = fPlanes[ip]->GetHits();

	Double_t zPos = fPlanes[ip]->GetZpos();
	Double_t dzPos = fPlanes[ip]->GetDzpos();

	// first loop over hits with in a single plane
	for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
	  // iphit is hit # within a plane
	  THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);

	  fTOFPInfo.push_back(TOFPInfo());
	  // Can remove these as we will initialize in the constructor
	  //	  fTOFPInfo[ihhit].time_pos = -99.0;
	  //	  fTOFPInfo[ihhit].time_neg = -99.0;
	  //	  fTOFPInfo[ihhit].keep_pos = kFALSE;
	  //	  fTOFPInfo[ihhit].keep_neg = kFALSE;
	  fTOFPInfo[ihhit].scin_pos_time = 0.0;
	  fTOFPInfo[ihhit].scin_neg_time = 0.0;
	  fTOFPInfo[ihhit].hit = hit;
	  fTOFPInfo[ihhit].planeIndex = ip;
	  fTOFPInfo[ihhit].hitNumInPlane = iphit;
	  fTOFPInfo[ihhit].onTrack = kFALSE;

	  Int_t paddle = hit->GetPaddleNumber()-1;
	  Double_t zposition = zPos + (paddle%2)*dzPos;

	  Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() *
	    ( zposition ); // Line 183

	  Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() *
	    ( zposition ); // Line 184

	  Double_t scinTrnsCoord, scinLongCoord;
	  if ( ( ip == 0 ) || ( ip == 2 ) ){ // !x plane. Line 185
	    scinTrnsCoord = xHitCoord;
	    scinLongCoord = yHitCoord;
	  } else if ( ( ip == 1 ) || ( ip == 3 ) ){ // !y plane. Line 188
	    scinTrnsCoord = yHitCoord;
	    scinLongCoord = xHitCoord;
	  } else { return -1; } // Line 195

	  fTOFPInfo[ihhit].scinTrnsCoord = scinTrnsCoord;
	  fTOFPInfo[ihhit].scinLongCoord = scinLongCoord;

	  Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();

	  // Index to access the 2d arrays of paddle/scintillator properties
	  Int_t fPIndex = GetScinIndex(ip,paddle);
          
	  if ( TMath::Abs( scinCenter - scinTrnsCoord ) <
	       ( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293

	    fTOFPInfo[ihhit].onTrack = kTRUE;
	    Double_t zcor = zposition/(29.979*fBetaNominal)*
	      TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
			  + theTrack->GetPhi()*theTrack->GetPhi());
	    fTOFPInfo[ihhit].zcor = zcor;
	    if (fCosmicFlag) {
	      Double_t zcor = -zposition/(29.979*1.0)*
		TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
			    + theTrack->GetPhi()*theTrack->GetPhi());
	      fTOFPInfo[ihhit].zcor = zcor;
	    }
	    Double_t tdc_pos = hit->GetPosTDC();
	    if(tdc_pos >=fScinTdcMin && tdc_pos <= fScinTdcMax ) {
	      Double_t adc_pos = hit->GetPosADC();
	      Double_t adcamp_pos = hit->GetPosADCpeak();
	      Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord;
	      fTOFPInfo[ihhit].pathp = pathp;
	      Double_t timep = tdc_pos*fScinTdcToTime;
	      if(fTofUsingInvAdc) {
		timep -= fHodoPosInvAdcOffset[fPIndex]
		  + pathp/fHodoPosInvAdcLinear[fPIndex]
		  + fHodoPosInvAdcAdc[fPIndex]
		  /TMath::Sqrt(TMath::Max(20.0*.020,adc_pos));
	      } else {
	        //Double_t tw_corr_pos = fHodoPos_c1[fPIndex]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) -  fHodoPos_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
		Double_t tw_corr_pos=0.;
		pathp=scinLongCoord;
		if (adcamp_pos>0) tw_corr_pos = 1./pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) -  1./pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);            
		timep += -tw_corr_pos + fHodo_LCoeff[fPIndex]+ pathp/fHodoVelFit[fPIndex];
	      }
	      fTOFPInfo[ihhit].scin_pos_time = timep;
 	      timep -= zcor;
	      fTOFPInfo[ihhit].time_pos = timep;

              hTime->Fill(timep);
	    }
	    Double_t tdc_neg = hit->GetNegTDC();
	    if(tdc_neg >=fScinTdcMin && tdc_neg <= fScinTdcMax ) {
	      Double_t adc_neg = hit->GetNegADC();
	      Double_t adcamp_neg = hit->GetNegADCpeak();
	      Double_t pathn =  scinLongCoord - fPlanes[ip]->GetPosRight();
	      fTOFPInfo[ihhit].pathn = pathn;
	      Double_t timen = tdc_neg*fScinTdcToTime;
	      if(fTofUsingInvAdc) {
		timen -= fHodoNegInvAdcOffset[fPIndex]
		  + pathn/fHodoNegInvAdcLinear[fPIndex]
		  + fHodoNegInvAdcAdc[fPIndex]
		  /TMath::Sqrt(TMath::Max(20.0*.020,adc_neg));
	      } else {
		pathn=scinLongCoord ;
		Double_t tw_corr_neg =0 ;
		if (adcamp_neg >0) tw_corr_neg= 1./pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[fPIndex]) -  1./pow(200./fTdc_Thrs, fHodoNeg_c2[fPIndex]);              
		timen += -tw_corr_neg- 2*fHodoCableFit[fPIndex] + fHodo_LCoeff[fPIndex]- pathn/fHodoVelFit[fPIndex];

	      }
	      fTOFPInfo[ihhit].scin_neg_time = timen;
	      timen -=  zcor;
	      fTOFPInfo[ihhit].time_neg = timen;
              hTime->Fill(timen);
	    }
	  } // condition for cenetr on a paddle
	  ihhit++;
	} // First loop over hits in a plane <---------

	//-----------------------------------------------------------------------------------------------
	//------------- First large loop over scintillator hits ends here --------------------
	//-----------------------------------------------------------------------------------------------
      }
      Int_t nhits=ihhit;


      if(0.5*hTime->GetMaximumBin() > 0) {
	Double_t tmin = 0.5*hTime->GetMaximumBin();
	
	for(Int_t ih = 0; ih < nhits; ih++) { // loop over all scintillator hits
	  if ( ( fTOFPInfo[ih].time_pos > (tmin-fTofTolerance) ) && ( fTOFPInfo[ih].time_pos < ( tmin + fTofTolerance ) ) ) {
	    fTOFPInfo[ih].keep_pos=kTRUE;
	  }
	  if ( ( fTOFPInfo[ih].time_neg > (tmin-fTofTolerance) ) && ( fTOFPInfo[ih].time_neg < ( tmin + fTofTolerance ) ) ){
	    fTOFPInfo[ih].keep_neg=kTRUE;
	  }
	}
      }

      //---------------------------------------------------------------------------------------------
      // ---------------------- Scond loop over scint. hits in a plane ------------------------------
      //---------------------------------------------------------------------------------------------
      for(Int_t ih=0; ih < nhits; ih++) {
	THcHodoHit *hit = fTOFPInfo[ih].hit;
	Int_t iphit = fTOFPInfo[ih].hitNumInPlane;
	Int_t ip = fTOFPInfo[ih].planeIndex;
	//         fDumpOut << " looping over hits = " << ih << " plane = " << ip+1 << endl;
	GoodFlags flags;
	// Flags are used by THcHodoEff
	fGoodFlags[itrack][ip].push_back(flags);
	fGoodFlags[itrack][ip][iphit].onTrack = kFALSE;
	fGoodFlags[itrack][ip][iphit].goodScinTime = kFALSE;
	fGoodFlags[itrack][ip][iphit].goodTdcNeg = kFALSE;
	fGoodFlags[itrack][ip][iphit].goodTdcPos = kFALSE;

	fTOFCalc.push_back(TOFCalc());
	// Do we set back to false for each track, or just once per event?
	fTOFCalc[ih].good_scin_time = kFALSE;
	// These need a track index too to calculate efficiencies
	fTOFCalc[ih].good_tdc_pos = kFALSE;
	fTOFCalc[ih].good_tdc_neg = kFALSE;
	fTOFCalc[ih].pindex = ip;

	Int_t paddle = hit->GetPaddleNumber()-1;
	fTOFCalc[ih].hit_paddle = paddle;
	fTOFCalc[ih].good_raw_pad = paddle;

	//	Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
	//	Double_t scinTrnsCoord = fTOFPInfo[ih].scinTrnsCoord;
	//	Double_t scinLongCoord = fTOFPInfo[ih].scinLongCoord;

	Int_t fPIndex = GetScinIndex(ip,paddle);

	if (fTOFPInfo[ih].onTrack) {
	  fGoodFlags[itrack][ip][iphit].onTrack = kTRUE;
	  if ( fTOFPInfo[ih].keep_pos ) { // 301
	    fTOFCalc[ih].good_tdc_pos = kTRUE;
	    fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
	  }
	  if ( fTOFPInfo[ih].keep_neg ) { //
	    fTOFCalc[ih].good_tdc_neg = kTRUE;
	    fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
	  }
	  // ** Calculate ave time for scin and error.
	  if ( fTOFCalc[ih].good_tdc_pos ){
	    if ( fTOFCalc[ih].good_tdc_neg ){
	      fTOFCalc[ih].scin_time  = ( fTOFPInfo[ih].scin_pos_time +
					  fTOFPInfo[ih].scin_neg_time ) / 2.;
	      fTOFCalc[ih].scin_time_fp  = ( fTOFPInfo[ih].time_pos +
					     fTOFPInfo[ih].time_neg ) / 2.;
	      
	      if (fTofUsingInvAdc){
		fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] +
						       fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex] )/2.;
	      }
	      else {
		fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoSigmaPos[fPIndex] * fHodoSigmaPos[fPIndex] +
						       fHodoSigmaNeg[fPIndex] * fHodoSigmaNeg[fPIndex] )/2.;
	      }

	      fTOFCalc[ih].good_scin_time = kTRUE;
	      fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
	    } else{
	      fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_pos_time;
	      fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_pos;
	      
	      if (fTofUsingInvAdc){
		fTOFCalc[ih].scin_sigma = fHodoPosSigma[fPIndex];
	      }
	      else{
		fTOFCalc[ih].scin_sigma = fHodoSigmaPos[fPIndex];
	      }
	      
	      fTOFCalc[ih].good_scin_time = kTRUE;
	      fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
	    }
	  } else {
	    if ( fTOFCalc[ih].good_tdc_neg ){
	      fTOFCalc[ih].scin_time = fTOFPInfo[ih].scin_neg_time;
	      fTOFCalc[ih].scin_time_fp = fTOFPInfo[ih].time_neg;
	      if (fTofUsingInvAdc){
		fTOFCalc[ih].scin_sigma = fHodoNegSigma[fPIndex];
	      }
	      else{
		fTOFCalc[ih].scin_sigma = fHodoSigmaNeg[fPIndex];
	      }
	      fTOFCalc[ih].good_scin_time = kTRUE;
	      fGoodFlags[itrack][ip][iphit].goodScinTime = kTRUE;
	    }
	  } // In h_tof.f this includes the following if condition for time at focal plane
	    // // because it is written in FORTRAN code

	    // c     Get time at focal plane
	  if ( fTOFCalc[ih].good_scin_time ){

	    // scin_time_fp doesn't need to be an array
	    // Is this any different than the average of time_pos and time_neg?
	    //	    Double_t scin_time_fp = ( fTOFPInfo[ih].time_pos +
	    //				      fTOFPInfo[ih].time_neg ) / 2.;
	    Double_t scin_time_fp = fTOFCalc[ih].scin_time_fp;

	    sumFPTime = sumFPTime + scin_time_fp;
	    nFPTime ++;

	    fSumPlaneTime[ip] = fSumPlaneTime[ip] + scin_time_fp;
	    fNPlaneTime[ip] ++;
	    fNScinHit[itrack] ++;

	    if ( ( fTOFCalc[ih].good_tdc_pos ) && ( fTOFCalc[ih].good_tdc_neg ) ){
	      nPmtHit[itrack] = nPmtHit[itrack] + 2;
	    } else {
	      nPmtHit[itrack] = nPmtHit[itrack] + 1;
	    }

	    fdEdX[itrack].push_back(0.0);

	    // --------------------------------------------------------------------------------------------
	    if ( fTOFCalc[ih].good_tdc_pos ){
	      if ( fTOFCalc[ih].good_tdc_neg ){
		fdEdX[itrack][fNScinHit[itrack]-1]=
		  TMath::Sqrt( TMath::Max( 0., hit->GetPosADC() * hit->GetNegADC() ) );
	      } else{
		fdEdX[itrack][fNScinHit[itrack]-1]=
		  TMath::Max( 0., hit->GetPosADC() );
	      }
	    } else{
	      if ( fTOFCalc[ih].good_tdc_neg ){
		fdEdX[itrack][fNScinHit[itrack]-1]=
		  TMath::Max( 0., hit->GetNegADC() );
	      } else{
		fdEdX[itrack][fNScinHit[itrack]-1]=0.0;
	      }
	    }
	    // --------------------------------------------------------------------------------------------


	  } // time at focal plane condition
	} // on track condition

	  // ** See if there are any good time measurements in the plane.
	if ( fTOFCalc[ih].good_scin_time ){
	  fGoodPlaneTime[ip] = kTRUE;
	  fTOFCalc[ih].dedx = fdEdX[itrack][fNScinHit[itrack]-1];
	} else {
	  fTOFCalc[ih].dedx = 0.0;
	}
	
      } // Second loop over hits of a scintillator plane ends here
      theTrack->SetGoodPlane3( fGoodPlaneTime[2] ? 1 : 0 );
      if (fNumPlanesBetaCalc==4) theTrack->SetGoodPlane4( fGoodPlaneTime[3] ? 1 : 0 );
      //
      //------------------------------------------------------------------------------
      //------------------------------------------------------------------------------
      //------------------------------------------------------------------------------
      //------------------------------------------------------------------------------
      //------------------------------------------------------------------------------
      //------------------------------------------------------------------------------
      //------------------------------------------------------------------------------
      //------------------------------------------------------------------------------

      // * * Fit beta if there are enough time measurements (one upper, one lower)
      // From h_tof_fit
      if ( ( ( fGoodPlaneTime[0] ) || ( fGoodPlaneTime[1] ) ) &&
	   ( ( fGoodPlaneTime[2] ) || ( fGoodPlaneTime[3] ) ) ){

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

	for(Int_t ih=0; ih < nhits; ih++) {
	  Int_t ip = fTOFPInfo[ih].planeIndex;

	  if ( fTOFCalc[ih].good_scin_time ) {

	    Double_t scinWeight = 1 / ( fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma );
	    Double_t zPosition = ( fPlanes[ip]->GetZpos()
				   +( fTOFCalc[ih].hit_paddle % 2 ) *
				   fPlanes[ip]->GetDzpos() );

	    sumW  += scinWeight;
	    sumT  += scinWeight * fTOFCalc[ih].scin_time;
	    sumZ  += scinWeight * zPosition;
	    sumZZ += scinWeight * ( zPosition * zPosition );
	    sumTZ += scinWeight * zPosition * fTOFCalc[ih].scin_time;

	  } // condition of good scin time
	} // loop over hits

	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 ) ) {

	  beta = tmp / tmpDenom;
	  betaChiSq = 0.;

	  for(Int_t ih=0; ih < nhits; ih++) {
	    Int_t ip = fTOFPInfo[ih].planeIndex;

	    if ( fTOFCalc[ih].good_scin_time ){

	      Double_t zPosition = ( fPlanes[ip]->GetZpos() + ( fTOFCalc[ih].hit_paddle % 2 ) *
				     fPlanes[ip]->GetDzpos() );
	      Double_t timeDif = ( fTOFCalc[ih].scin_time - t0 );
	      betaChiSq += ( ( zPosition / beta - timeDif ) *
			     ( zPosition / beta - timeDif ) )  /
		( fTOFCalc[ih].scin_sigma * fTOFCalc[ih].scin_sigma );

	    } // condition for good scin time
	  } // loop over hits

	  Double_t pathNorm = TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
					   theTrack->GetPhi()   * theTrack->GetPhi() );
	  // Take angle into account
	  beta = beta / pathNorm;
	  beta = beta / 29.979;    // velocity / c

	} // condition for fTmpDenom
	else {
	  beta = 0.;
	  betaChiSq = -2.;
	} // else condition for fTmpDenom
      } else {
	beta = 0.;
	betaChiSq = -1;
      }

      if ( nFPTime != 0 ){
      	timeAtFP[itrack] = ( sumFPTime / nFPTime );
      }
      //
      // ---------------------------------------------------------------------------

      Double_t FPTimeSum=0.0;
      Int_t nFPTimeSum=0;
      for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
	if ( fNPlaneTime[ip] != 0 ){
	  fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
	  FPTimeSum += fSumPlaneTime[ip];
	  nFPTimeSum += fNPlaneTime[ip];
	} else {
	  fFPTime[ip] = 1000. * ( ip + 1 );
	}
      }
      Double_t fptime=-1000;
      if (nFPTimeSum>0) fptime = FPTimeSum/nFPTimeSum;
      fFPTimeAll = fptime;
      Double_t dedx=0.0;
      for(UInt_t ih=0;ih<fTOFCalc.size();ih++) {
	if(fTOFCalc[ih].good_scin_time) {
	  dedx = fTOFCalc[ih].dedx;
	  break;
	}
      }
      theTrack->SetDedx(dedx);
      theTrack->SetFPTime(fptime);
      theTrack->SetBeta(beta);
      theTrack->SetBetaChi2( betaChiSq );
      theTrack->SetNPMT(nPmtHit[itrack]);
      theTrack->SetFPTime( timeAtFP[itrack]);


    } // Main loop over tracks ends here.

  } // If condition for at least one track


  //OriginalTrackEffTest();
  TrackEffTest();


  return 0;

}
//
void THcHodoscope::TrackEffTest(void)
{
  //Double_t PadLow[4];
  //Double_t PadHigh[4];
  // assume X planes are 0,2 and Y planes are 1,3
  std::array<int,4> PadLow = {fxLoScin[0], fyLoScin[0], fxLoScin[1], fyLoScin[1]};
  std::array<int,4> PadHigh = {fxHiScin[0], fyHiScin[0], fxHiScin[1], fyHiScin[1]};
  
  std::array<double,4> PadPosLo;
  std::array<double,4> PadPosHi;

  for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
    Double_t lowtemp = fPlanes[ip]->GetPosCenter(PadLow[ip] - 1) + fPlanes[ip]->GetPosOffset();
    Double_t hitemp  = fPlanes[ip]->GetPosCenter(PadHigh[ip] - 1) + fPlanes[ip]->GetPosOffset();
    if (lowtemp < hitemp) {
      PadPosLo[ip]=lowtemp;
      PadPosHi[ip]=hitemp;
    } else {
      PadPosLo[ip]=hitemp;
      PadPosHi[ip]=lowtemp;
    }
  }  
  //
  const Int_t           MaxNClus = 5;
  std::vector<Int_t>    iw(MaxNClus, 0);
  std::vector<Double_t> dw(MaxNClus, 0);
  for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
    fNClust.push_back(0);
    fClustSize.push_back(iw);
    fClustPos.push_back(dw);
  }

  //{
  //  std::vector<int>  test_vector = {1,2,5,6,7, 9,10, 20};
  //  auto test_res = hcana::find_discontinuity(test_vector);
  //  std::cout << " find_discontinuity    test: " << *test_res << "\n";
  //  std::cout << " count_discontinuities test: " << hcana::count_discontinuities(test_vector)
  //            << "\n";
  //  auto cont_ranges  = hcana::get_discontinuities(test_vector);
  //  int ir = 0;
  //    std::cout << "ranges: \n";
  //  for(auto r : cont_ranges) {
  //    for(auto v = r.first; v< r.second; v++) {
  //      std::cout << *v << " ";
  //    }
  //    std::cout << "\n";
  //  }
  //}
  //{
  //  std::vector<int>  test_vector = {0,2,4,8,10,12,13, 16,18,20,23};
  //  auto              test_res    = hcana::find_discontinuity(
  //      test_vector, [](const int& x1, const int& x2) { return x1 + 2 == x2; });
  //  std::cout << " find_discontinuity test2: " << *test_res << "\n";
  //  std::cout << " count_discontinuities test2: "
  //            << hcana::count_discontinuities(
  //                   test_vector, [](const int& x1, const int& x2) { return x1 + 2 == x2; })
  //            << "\n";
  //}

  using hit_iter = std::vector<THcHodoHit*>::iterator;
  using HitRangeVector = typename std::vector<std::pair<hit_iter, hit_iter>>;
  using ClusterPositions = typename std::vector<double>;
  std::vector<std::vector<THcHodoHit*>> hit_vecs;
  std::vector<HitRangeVector>           hit_clusters;
  std::vector<ClusterPositions>         pos_clusters;
  std::vector<std::vector<int>>         good_clusters;
  std::array<int,4>                     n_good_clusters;
  // Loop over each plane
  for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){

    // Get the hit array
    TClonesArray* hodoHits = fPlanes[ip]->GetHits();

    // Create vector for good hits.
    hit_vecs.push_back(std::vector<THcHodoHit*>());
    for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
      THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
      // keep only those with "two good times"
      if ( hit->GetTwoGoodTimes() ) {
        hit_vecs[ip].push_back(hit);
      }
    }

    // Cluster (group adjacent) hits 
    if( hit_vecs[ip].size() > 0 ) {
      auto& hv = hit_vecs[ip];
      hcana::sort(hv, [](const THcHodoHit* h1, const THcHodoHit* h2) {
            return h1->GetPaddleNumber() < h2->GetPaddleNumber();
      });

      // get the clusters for this layer
      hit_clusters.push_back(
          hcana::get_discontinuities(hv, [](const THcHodoHit* h1, const THcHodoHit* h2) {
            return h1->GetPaddleNumber() + 1 == h2->GetPaddleNumber();
          }));
    } else {
      hit_clusters.push_back( HitRangeVector{} );
    } 

    // Calculate each cluster's mean position (in one coordinate)
    ClusterPositions clust_positions;
    std::vector<int> clust_bound;
    for(auto& r : hit_clusters.back()) {
      double a_pos  = 0.0;
      int    n_hit = 0;
      for(auto chit = r.first; chit < r.second; chit++) {
        a_pos += fPlanes[ip]->GetPosCenter((*chit)->GetPaddleNumber()-1) + fPlanes[ip]->GetPosOffset();
        n_hit++;
      }
      // take the average position
      double avg_pos = a_pos / double(n_hit);
      // Calculate if it is the bound by the upper and lower limits where we expect
      // full tracks to be reconstructed.
      int is_bound = int((avg_pos >= PadPosLo[ip]) &&  (avg_pos <= PadPosHi[ip]));
      clust_bound.push_back(is_bound);
      if(is_bound) {
      // we are only interested in clusters positioned inthe "sweet spot"
        clust_positions.push_back(avg_pos);
        n_good_clusters[ip]++;
        if (n_hit > 10) {
          _logger->warn("cluster in hodoscope track efficiency has {} hits", n_hit);
        }
      }
    }
    pos_clusters.push_back(clust_positions);
    good_clusters.push_back(clust_bound);
    //_logger->info("{} clusters in Hod plane {}", hit_clusters.back().size(), ip );
  }

  bool has_both_x_clusters          = (n_good_clusters[0] > 0) && (n_good_clusters[2] > 0);
  bool has_both_y_clusters          = (n_good_clusters[1] > 0) && (n_good_clusters[3] > 0);
  bool at_least_one_good_x_clusters = (n_good_clusters[0] > 0) || (n_good_clusters[2] > 0);
  bool at_least_one_good_y_clusters = (n_good_clusters[1] > 0) || (n_good_clusters[3] > 0);
  bool good_x_match = false;
  bool good_y_match = false;

  // Superficial matching. Find out if clusters in front and back could possibly belong to the same
  // track. This done by checking to see if the position differences  are less than 10 cm;
  // ie, |x1-x2|<10  or |y1-y2| <10
  // TODO do actual matching with tracks elsewhere
  if(has_both_x_clusters) {
    for(auto x1_pos : pos_clusters[0] ) {
      for(auto x2_pos : pos_clusters[2] ) {
        if( std::abs(x1_pos-x2_pos) < 10.0) {
          good_x_match = true;
          break;
        }
      }
    }
  }
  if(has_both_y_clusters) {
    for(auto y1_pos : pos_clusters[1] ) {
      for(auto y2_pos : pos_clusters[3] ) {
        if( std::abs(y1_pos-y2_pos) < 10.0) {
          good_y_match = true;
          break;
        }
      }
    }
  }

  fGoodScinHits = 0;
  // (2+2) 4 good hits 
  if( good_y_match &&  good_x_match )  {
    fGoodScinHits = 1;
  }
  // (2+1) 3 good hits
  if( at_least_one_good_y_clusters &&  good_x_match )  {
    fGoodScinHits = 1;
  }
  // (1+2) 3 good hits
  if( at_least_one_good_x_clusters &&  good_y_match )  {
    fGoodScinHits = 1;
  }


  //for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
  //  //std::cout << ip << " " << fPlanes[ip]->GetHits()->GetEntries() << " hits\n" ;
  //  TClonesArray* hodoHits    = fPlanes[ip]->GetHits();
  //  Int_t         prev_padnum = -100;

  //  for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
  //    THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
  //    Int_t padnum  = hit->GetPaddleNumber();
  //    if ( hit->GetTwoGoodTimes() ) {
  //      if ( padnum==prev_padnum+1 ) {
  //        fClustSize[ip][fNClust[ip]-1]=fClustSize[ip][fNClust[ip]-1]+1;
  //        fClustPos[ip][fNClust[ip]-1]=fClustPos[ip][fNClust[ip]-1]+fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset();
  //        //  cout << "Add to cluster  pl = " << ip+1 << " hit = " << iphit << " pad = " << padnum << " clus =  " << fNClust[ip] << " cl size = " << fClustSize[ip][fNClust[ip]-1] << " pos " << fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset() << endl;
  //      } else {
  //        if (fNClust[ip]<MaxNClus) fNClust[ip]++;
  //        fClustSize[ip][fNClust[ip]-1]=1;
  //        fClustPos[ip][fNClust[ip]-1]=fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset();
  //        //  cout << " New clus pl = " << ip+1 << " hit = " << iphit << " pad = " << padnum << " clus = " << fNClust[ip] << " cl size = " << fClustSize[ip][fNClust[ip]] << " pos " << fPlanes[ip]->GetPosCenter(padnum-1)+ fPlanes[ip]->GetPosOffset() << endl;
  //      }
  //      prev_padnum=padnum;
  //    }
  //  }
  //}
  //
  //Bool_t inside_bound[4]={kFALSE,kFALSE,kFALSE,kFALSE};
  //for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {	 
  //  for(Int_t ic = 0; ic <fNClust[ip] ; ic++ ) {
  //    fClustPos[ip][ic]=fClustPos[ip][ic]/fClustSize[ip][ic];
  //    inside_bound[ip] = fClustPos[ip][ic]>=PadPosLo[ip] &&  fClustPos[ip][ic]<=PadPosHi[ip];
  //    //cout << "plane = " << ip+1 << " Cluster = " << ic+1 << " size = " << fClustSize[ip][ic]<< " pos = " << fClustPos[ip][ic] << " inside = " << inside_bound[ip] << " lo = " << PadPosLo[ip]<< " hi = " << PadPosHi[ip]<< endl;
  //  }
  //}
  //
  //std::array<int,4> good_for_track_test = {0,0,0,0};
  //Int_t sum_good_track_test=0;
  //for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
  //  if (fNClust[ip]==1 && inside_bound[ip] && fClustSize[ip][0]<=2) good_for_track_test[ip]=1;
  //  //cout << " good for track = " << good_for_track_test[ip] << endl;
  //  sum_good_track_test+=good_for_track_test[ip];
  //}	 

  // these should be input parameters
  //Double_t trackeff_scint_ydiff_max= 10. ;
  //Double_t trackeff_scint_xdiff_max= 10. ;

  //Bool_t xdiffTest=kFALSE;
  //Bool_t ydiffTest=kFALSE;

  //fGoodScinHits = 0;
  //if (fTrackEffTestNScinPlanes == 4) {
  //  if (fTrackEffTestNScinPlanes==sum_good_track_test) {
  //    xdiffTest= TMath::Abs(fClustPos[0][0]-fClustPos[2][0])<trackeff_scint_xdiff_max;
  //    ydiffTest= TMath::Abs(fClustPos[1][0]-fClustPos[3][0])<trackeff_scint_ydiff_max;
  //    if (xdiffTest && ydiffTest) fGoodScinHits = 1;
  //  }
  //}
  ////
  //if (fTrackEffTestNScinPlanes == 3) {
  //  if (fTrackEffTestNScinPlanes==sum_good_track_test) {
  //    if(good_for_track_test[0]==1&&good_for_track_test[2]==1) {
  //      xdiffTest= TMath::Abs(fClustPos[0][0]-fClustPos[2][0])<trackeff_scint_xdiff_max;
  //      ydiffTest=kTRUE;
  //    }
  //    if (good_for_track_test[1]==1&&good_for_track_test[3]==1) {
  //      xdiffTest=kTRUE;
  //      ydiffTest= TMath::Abs(fClustPos[1][0]-fClustPos[3][0])<trackeff_scint_ydiff_max;
  //    }
  //    if (xdiffTest && ydiffTest) fGoodScinHits = 1;
  //  }
  //  if (sum_good_track_test==4) {
  //    xdiffTest= TMath::Abs(fClustPos[0][0]-fClustPos[2][0])<trackeff_scint_xdiff_max;
  //    ydiffTest= TMath::Abs(fClustPos[1][0]-fClustPos[3][0])<trackeff_scint_ydiff_max;
  //    if (xdiffTest && ydiffTest) fGoodScinHits = 1;
  //  }
  //}
  //       
  //	cout << " good scin = " << fGoodScinHits << " " << sum_good_track_test << " " << xdiffTest  << " " << ydiffTest<< endl;
  //cout << " ************" << endl;
  //
}
//
//
//
void THcHodoscope::OriginalTrackEffTest(void)
{
  /**
      Translation of h_track_tests.f file for tracking efficiency
  */

  //************************now look at some hodoscope tests
  //  *second, we move the scintillators.  here we use scintillator cuts to see
  //  *if a track should have been found.
  cout << " enter track eff" <<  fNumPlanesBetaCalc << endl;
  for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
    cout << " loop over planes " <<  ip+1 << endl;

    TClonesArray* hodoHits = fPlanes[ip]->GetHits();
    //    TClonesArray* scinPosTDC = fPlanes[ip]->GetPosTDC();
    //    TClonesArray* scinNegTDC = fPlanes[ip]->GetNegTDC();

    fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
    cout << " hits =  " << fNScinHits[ip]  << endl;
    for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
      Int_t paddle = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;

      fScinHitPaddle[ip][paddle] = 1;
      cout << " hit =  " << iphit+1 << " " << paddle+1   << endl;
    }
  }

  //  *next, look for clusters of hits in a scin plane.  a cluster is a group of
  //  *adjacent scintillator hits separated by a non-firing scintillator.
  //  *Wwe count the number of three adjacent scintillators too.  (A signle track
  //  *shouldn't fire three adjacent scintillators.

  for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
    // Planes ip = 0 = 1X
    // Planes ip = 2 = 2X
    fNClust.push_back(0);
    fThreeScin.push_back(0);
  }

  // *look for clusters in x planes... (16 scins)  !this assume both x planes have same
  // *number of scintillators.
  cout << " looking for cluster in x planes" << endl;
  Int_t icount;
  for (Int_t ip = 0; ip < 3; ip +=2 ) {
    icount = 0;
    if ( fScinHitPaddle[ip][0] == 1 ) icount ++;
    cout << "plane =" << ip <<  "check if paddle 1 hit " << icount << endl;

    for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 1; ipaddle++ ){
      // !look for number of clusters of 1 or more hits
      if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
	   ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
	icount ++;
      cout << " paddle =  " << ipaddle+1 << " " << icount << endl;
 
    } // Loop over  paddles
    cout << "Two  cluster in plane =  " << ip+1 << " " << icount << endl;
    fNClust[ip] = icount;
    icount = 0;

    for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 2; ipaddle++ ){
      // !look for three or more adjacent hits

      if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
	   ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
	   ( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
	icount ++;
    } // Second loop over paddles
    cout << "Three  clusters in plane =  " << ip+1 << " " << icount << endl;

    if ( icount > 0 )
      fThreeScin[ip] = 1;

  } // Loop over X plane

  // *look for clusters in y planes... (10 scins)  !this assume both y planes have same
  // *number of scintillators.
  cout << " looking for cluster in y planes" << endl;

  for (Int_t ip = 1; ip < fNumPlanesBetaCalc; ip +=2 ) {
    // Planes ip = 1 = 1Y
    // Planes ip = 3 = 2Y

    icount = 0;
    if ( fScinHitPaddle[ip][0] == 1 ) icount ++;
    cout << "plane =" << ip <<  "check if paddle 1 hit " << icount << endl;

    for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 1; ipaddle++ ){
      //  !look for number of clusters of 1 or more hits

      if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
	   ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
	icount ++;
      cout << " paddle =  " << ipaddle+1 << " " << icount << endl;

    } // Loop over Y paddles
    cout << "Two  cluster in plane =  " << ip+1 << " " << icount << endl;

    fNClust[ip] = icount;
    icount = 0;

    for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[ip] - 2; ipaddle++ ){
      // !look for three or more adjacent hits

      if ( ( fScinHitPaddle[ip][ipaddle] == 1 ) &&
	   ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) &&
	   ( fScinHitPaddle[ip][ipaddle + 2] == 1 ) )
	icount ++;

    } // Second loop over Y paddles
    cout << "Three  clusters in plane =  " << ip+1 << " " << icount << endl;

    if ( icount > 0 )
      fThreeScin[ip] = 1;

  }// Loop over Y planes


  // *next we mask out the edge scintillators, and look at triggers that happened
  // *at the center of the acceptance.  To change which scins are in the mask
  // *change the values of h*loscin and h*hiscin in htracking.param

  //      fGoodScinHits = 0;
  for (Int_t ifidx = fxLoScin[0]; ifidx < (Int_t) fxHiScin[0]; ifidx ++ ){
    fGoodScinHitsX.push_back(0);
  }

  fHitSweet1X=0;
  fHitSweet2X=0;
  fHitSweet1Y=0;
  fHitSweet2Y=0;
  // *first x plane.  first see if there are hits inside the scin region
  for (Int_t ifidx = fxLoScin[0]-1; ifidx < fxHiScin[0]; ifidx ++ ){
    if ( fScinHitPaddle[0][ifidx] == 1 ){
      fHitSweet1X = 1;
      fSweet1XScin = ifidx + 1;
    }
  }

  // *  next make sure nothing fired outside the good region
  for (Int_t ifidx = 0; ifidx < fxLoScin[0]-1; ifidx ++ ){
    if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
  }
  for (Int_t ifidx = fxHiScin[0]; ifidx < (Int_t) fNPaddle[0]; ifidx ++ ){
    if ( fScinHitPaddle[0][ifidx] == 1 ){ fHitSweet1X = -1; }
  }

  // *second x plane.  first see if there are hits inside the scin region
  for (Int_t ifidx = fxLoScin[1]-1; ifidx < fxHiScin[1]; ifidx ++ ){
    if ( fScinHitPaddle[2][ifidx] == 1 ){
      fHitSweet2X = 1;
      fSweet2XScin = ifidx + 1;
    }
  }
  // *  next make sure nothing fired outside the good region
  for (Int_t ifidx = 0; ifidx < fxLoScin[1]-1; ifidx ++ ){
    if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
  }
  for (Int_t ifidx = fxHiScin[1]; ifidx < (Int_t) fNPaddle[2]; ifidx ++ ){
    if ( fScinHitPaddle[2][ifidx] == 1 ){ fHitSweet2X = -1; }
  }

  // *first y plane.  first see if there are hits inside the scin region
  for (Int_t ifidx = fyLoScin[0]-1; ifidx < fyHiScin[0]; ifidx ++ ){
    if ( fScinHitPaddle[1][ifidx] == 1 ){
      fHitSweet1Y = 1;
      fSweet1YScin = ifidx + 1;
    }
  }
  // *  next make sure nothing fired outside the good region
  for (Int_t ifidx = 0; ifidx < fyLoScin[0]-1; ifidx ++ ){
    if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
  }
  for (Int_t ifidx = fyHiScin[0]; ifidx < (Int_t) fNPaddle[1]; ifidx ++ ){
    if ( fScinHitPaddle[1][ifidx] == 1 ){ fHitSweet1Y = -1; }
  }

  // *second y plane.  first see if there are hits inside the scin region
  for (Int_t ifidx = fyLoScin[1]-1; ifidx < fyHiScin[1]; ifidx ++ ){
    if ( fScinHitPaddle[3][ifidx] == 1 ){
      fHitSweet2Y = 1;
      fSweet2YScin = ifidx + 1;
    }
  }

  // *  next make sure nothing fired outside the good region
  for (Int_t ifidx = 0; ifidx < fyLoScin[1]-1; ifidx ++ ){
    if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
  }
  for (Int_t ifidx = fyHiScin[1]; ifidx < (Int_t) fNPaddle[3]; ifidx ++ ){
    if ( fScinHitPaddle[3][ifidx] == 1 ){ fHitSweet2Y = -1; }
  }

  fTestSum = fHitSweet1X + fHitSweet2X + fHitSweet1Y + fHitSweet2Y;

  // * now define a 3/4 or 4/4 trigger of only good scintillators the value
  // * is specified in htracking
  if ( fTestSum >= fTrackEffTestNScinPlanes ){
    fGoodScinHits = 1;
    for (Int_t ifidx = fxLoScin[0]; ifidx < fxHiScin[0]; ifidx ++ ){
      if ( fSweet1XScin == ifidx )
	fGoodScinHitsX[ifidx] = 1;
    }
  }

  // * require front/back hodoscopes be close to each other
  if ( ( fGoodScinHits == 1 ) && ( fTrackEffTestNScinPlanes == 4 ) ){
    if ( TMath::Abs( fSweet1XScin - fSweet2XScin ) > 3 )
      fGoodScinHits = 0;
    if ( TMath::Abs( fSweet1YScin - fSweet2YScin ) > 2 )
      fGoodScinHits = 0;
  }
//
}
//_____________________________________________________________________________
Int_t THcHodoscope::FineProcess( TClonesArray&  tracks  )
{
  Int_t Ntracks = tracks.GetLast()+1;   // Number of reconstructed tracks
  Double_t hitPos;
  Double_t hitDistance;
  Int_t ih=0;
  for (Int_t itrk=0; itrk<Ntracks; itrk++) {
    THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
    if (theTrack->GetIndex()==0) {
      fBeta=theTrack->GetBeta();
      fFPTimeAll=theTrack->GetFPTime();
      Double_t shower_track_enorm = theTrack->GetEnergy()/theTrack->GetP();
      for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
	Double_t pl_xypos=0;
	Double_t pl_zpos=0;
	Int_t   num_good_pad=0;
	TClonesArray* hodoHits = fPlanes[ip]->GetHits();
	for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
	  THcHodoHit *hit = fTOFPInfo[ih].hit;
	  if ( fTOFCalc[ih].good_tdc_pos && fTOFCalc[ih].good_tdc_neg ) {
	    Bool_t sh_pid=(shower_track_enorm > fTOFCalib_shtrk_lo && shower_track_enorm < fTOFCalib_shtrk_hi);
	    Bool_t beta_pid=( fBeta > fTOFCalib_beta_lo &&  fBeta < fTOFCalib_beta_hi);
	    // cer_pid is true if there is no Cherenkov detector
	    Bool_t cer_pid=(fCherenkov?(fCherenkov->GetCerNPE()>fTOFCalib_cer_lo):kTRUE);
	    if(fDumpTOF && Ntracks==1 && fGoodEventTOFCalib && sh_pid && beta_pid && cer_pid) {
	      fDumpOut << fixed << setprecision(2);
	      fDumpOut  << showpoint << " 1" << setw(3) << ip+1 << setw(3) << hit->GetPaddleNumber()  << setw(10) << hit->GetPosTDC()*fScinTdcToTime  << setw(10) << fTOFPInfo[ih].pathp << setw(10) << fTOFPInfo[ih].zcor  << setw(10) << fTOFPInfo[ih].time_pos << setw(10) << hit->GetPosADC() << endl;
	      fDumpOut  << showpoint << " 2" << setw(3) << ip+1 << setw(3) << hit->GetPaddleNumber() << setw(10) << hit->GetNegTDC()*fScinTdcToTime  << setw(10) << fTOFPInfo[ih].pathn << setw(10) << fTOFPInfo[ih].zcor  << setw(10) << fTOFPInfo[ih].time_neg << setw(10) << hit->GetNegADC() << endl;
	    }
	    Int_t padind = ((THcHodoHit*)hodoHits->At(iphit))->GetPaddleNumber()-1;
	    pl_xypos+=fPlanes[ip]->GetPosCenter(padind)+ fPlanes[ip]->GetPosOffset();
	    pl_zpos+=fPlanes[ip]->GetZpos()+ (padind%2)*fPlanes[ip]->GetDzpos();
	    num_good_pad++; 
	  } 
	  ih++;
	  //	  cout << ip << " " << iphit << " " <<  fGoodFlags[itrk][ip][iphit].goodScinTime << " " <<   fGoodFlags[itrk][ip][iphit].goodTdcPos << " " << fGoodFlags[itrk][ip][iphit].goodTdcNeg << endl;
	}
	hitDistance=kBig;
	if (num_good_pad !=0 ) {
	  pl_xypos=pl_xypos/num_good_pad;
	  pl_zpos=pl_zpos/num_good_pad;
	  hitPos = theTrack->GetY() + theTrack->GetPhi()*pl_zpos ;
	  if (ip%2 == 0)  hitPos = theTrack->GetX() + theTrack->GetTheta()*pl_zpos ;
	  hitDistance =  hitPos - pl_xypos;
	  fPlanes[ip]->SetTrackXPosition(theTrack->GetX() + theTrack->GetTheta()*pl_zpos );
	  fPlanes[ip]->SetTrackYPosition(theTrack->GetY() + theTrack->GetPhi()*pl_zpos );
	}
	//      cout << " ip " << ip << " " << hitPos << " " << pl_xypos << " " << pl_zpos << " " <<  hitDistance << endl;
	fPlanes[ip]->SetHitDistance(hitDistance);
      }
      if(ih>0&&fDumpTOF && Ntracks==1 && fGoodEventTOFCalib && shower_track_enorm > fTOFCalib_shtrk_lo && shower_track_enorm < fTOFCalib_shtrk_hi ) {
	fDumpOut << "0 "  << endl;
      }
    }
  }     //over tracks
  //
  return 0;
}
//_____________________________________________________________________________
Int_t THcHodoscope::GetScinIndex( Int_t nPlane, Int_t nPaddle ) {
  // GN: Return the index of a scintillator given the plane # and the paddle #
  // This assumes that both planes and
  // paddles start counting from 0!
  // Result also counts from 0.
  return fNPlanes*nPaddle+nPlane;
}
//_____________________________________________________________________________
Int_t THcHodoscope::GetScinIndex( Int_t nSide, Int_t nPlane, Int_t nPaddle ) {
  return nSide*fMaxHodoScin+fNPlanes*nPaddle+nPlane-1;
}
//_____________________________________________________________________________
Double_t THcHodoscope::GetPathLengthCentral() {