THcHodoscope.cxx

	    }
	    
	    Double_t tdc_neg = hit->GetNegTDC();
	    if(tdc_neg >=fScinTdcMin && tdc_neg <= fScinTdcMax ) {
	      Double_t adc_neg = hit->GetNegADC();
	      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,adc_neg));
	      } else {
		timen -= fHodoNegPhcCoeff[fPIndex]*
		  TMath::Sqrt(TMath::Max(0.0,adc_neg/fHodoNegMinPh[fPIndex]-1.0))
		  + pathn/fHodoVelLight[fPIndex]
		  + fHodoNegTimeOffset[fPIndex];
	      }
	      fTOFPInfo[ihhit].scin_neg_time = timen;
	      timen -=  zcor;
	      fTOFPInfo[ihhit].time_neg = timen;
	      
	      for ( Int_t k = 0; k < 200; k++ ){ // Line 230
		Double_t tmin = 0.5 * ( k + 1 );
		if ( ( timen > tmin ) && ( timen < ( tmin + fTofTolerance ) ) )
		  timehist[k] ++;
	      }
	    }
	  } // 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;
	
      // Find bin with the most hits
      Int_t jmax = 0; // Line 240
      Int_t maxhit = 0;
      for ( Int_t k = 0; k < 200; k++ ){
	if ( timehist[k] > maxhit ){
	  jmax = k+1;
	  maxhit = timehist[k];
	}
      }
	

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

	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) {
	  if ( fTOFPInfo[ih].keep_pos ) { // 301
	    fTOFCalc[ih].good_tdc_pos = kTRUE;
	    fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
	    if(fDumpTOF && ntracks==1) {
	      fDumpOut << "1 " << ip+1 << " " << paddle+1 << " " <<
		hit->GetPosTDC()*fScinTdcToTime << " " << fTOFPInfo[ih].pathp << 
		" "  << fTOFPInfo[ih].zcor << " " << fTOFPInfo[ih].time_pos <<
		" " << hit->GetPosADC() << endl;
	    }
	  }	    
	  if ( fTOFPInfo[ih].keep_neg ) { //
	    fTOFCalc[ih].good_tdc_neg = kTRUE;
	    fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
	    if(fDumpTOF && ntracks==1) {
	      fDumpOut << "2 " << ip+1 << " " << paddle+1 << " " <<
		hit->GetNegTDC()*fScinTdcToTime << " " << fTOFPInfo[ih].pathn << 
		" "  << fTOFPInfo[ih].zcor << " " << fTOFPInfo[ih].time_neg <<
		" " << hit->GetNegADC() << endl;
	    }
	  }	    

	  // ** 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.;
	      fTOFCalc[ih].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] + 
							fHodoNegSigma[fPIndex] * fHodoNegSigma[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;
	      fTOFCalc[ih].scin_sigma = fHodoPosSigma[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;
	      fTOFCalc[ih].scin_sigma = fHodoNegSigma[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 );
      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 < fNPlanes; 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 = FPTimeSum/nFPTimeSum;
      
      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

  //-----------------------------------------------------------------------
  //
  //   Trnslation 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.

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

    std::vector<Double_t> scin_temp;
    fScinHitPaddle.push_back(scin_temp); // Create array of hits per plane

    for (UInt_t ipaddle = 0; ipaddle < fNPaddle[0]; ipaddle++ ){
	  fScinHitPaddle[ip].push_back(0.0);
	  fScinHitPaddle[ip][ipaddle] = 0.0;	  
    }
  }

  for(Int_t ip = 0; ip < fNPlanes; ip++ ) {
    if (!fPlanes[ip])
      return -1;

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

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

  //  *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 < fNPlanes; ip++ ) {
    // Planes ip = 0 = 1X
    // Planes ip = 2 = 2X
    if (!fPlanes[ip]) return -1;
    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.
  Int_t icount;
  for (Int_t ip = 0; ip < 3; ip +=2 ) {
    icount = 0;
    if ( fScinHitPaddle[ip][0] == 1 )
      icount ++;

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

      if ( ( fScinHitPaddle[ip][ipaddle] == 0 ) &&
	   ( fScinHitPaddle[ip][ipaddle + 1] == 1 ) )
	icount ++;
      
    } // Loop over  paddles

    fNClust[ip] = icount;
    icount = 0;

    for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[0] - 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

    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.

  for (Int_t ip = 1; ip < 4; ip +=2 ) {
    // Planes ip = 1 = 1Y
    // Planes ip = 3 = 2Y
    if (!fPlanes[ip]) return -1;

    icount = 0;
    if ( fScinHitPaddle[ip][0] == 1 )
      icount ++;

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

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

    } // Loop over Y paddles

    fNClust[ip] = icount;
    icount = 0;

    for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[1] - 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

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

  }// Loop over Y planes

  // *now put some "tracking" like cuts on the hslopes, based only on scins...
  // *by "slope" here, I mean the difference in the position of scin hits in two
  // *like-planes.  For example, a track that those great straight through will 
  // *have a slope of zero.  If it moves one scin over from s1x to s2x it has an
  // *x-slope of 1...  I pick the minimum slope if there are multiple scin hits.

  Double_t bestXpScin = 100.0;
  Double_t bestYpScin = 100.0;

  for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[0]; ipaddle++ ){
    for (Int_t ipaddle2 = 0; ipaddle2 < (Int_t) fNPaddle[0]; ipaddle2++ ){

      if ( ( fScinHitPaddle[0][ipaddle] == 1 ) &&
	   ( fScinHitPaddle[2][ipaddle2] == 1 ) ){

	Double_t slope = TMath::Abs(ipaddle - ipaddle2);

	if ( slope < bestXpScin ) {
	  bestXpScin = slope;

	}
      }
    }  // Second loop over X paddles
  } // First loop over X paddles


  for (Int_t ipaddle = 0; ipaddle < (Int_t) fNPaddle[1]; ipaddle++ ){
    for (Int_t ipaddle2 = 0; ipaddle2 < (Int_t) fNPaddle[1]; ipaddle2++ ){

      if ( ( fScinHitPaddle[1][ipaddle] == 1 ) &&
	   ( fScinHitPaddle[3][ipaddle2] == 1 ) ){

	Double_t slope = TMath::Abs(ipaddle - ipaddle2);

	if ( slope < bestYpScin ) {
	  bestYpScin = slope;	
	}
      }
    }  // Second loop over Y paddles
  } // First loop over Y paddles

  // *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);
  }

  // *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.param...
  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;
  }


  if ( !fChern || !fShower ) { 
    return 0;    
  }

  
  if ( ( fGoodScinHits == 1 ) && ( fShower->GetNormETot() > fNormETot ) &&
       ( fChern->GetCerNPE() > fNCerNPE ) )
    fScinShould = 1;
  
  if ( ( fGoodScinHits == 1 ) && ( fShower->GetNormETot() > fNormETot ) &&
       ( fChern->GetCerNPE() > fNCerNPE ) && ( tracks.GetLast() + 1 > 0 ) ) {
      fScinDid = 1;
  }
	    if(fDumpTOF && ntracks==1) {
	      fDumpOut << "0 "  << endl;
	    }
  
  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() {
  return fPathLengthCentral;
}
ClassImp(THcHodoscope)
////////////////////////////////////////////////////////////////////////////////