Skip to content
Snippets Groups Projects
THcHodoscope.cxx 56.2 KiB
Newer Older
  • Learn to ignore specific revisions
  • 	  
    	  // ** Check if scin is on track
    
    	  if ( TMath::Abs( fScinCenter - fScinTrnsCoord ) >
    	       ( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
    
    	    // * * Check for good TDC
    
    	    if ( ( ((THcSignalHit*)scinPosTDC->At(iphit))->GetData() > fScinTdcMin ) &&
    		 ( ((THcSignalHit*)scinPosTDC->At(iphit))->GetData() < fScinTdcMax ) &&
    		 ( fTOFPInfo[iphit].keep_pos ) ) { // 301
    
    	      
    	      // ** Calculate time for each tube with a good tdc. 'pos' side first.
    
    	      fTOFCalc[ihhit].good_tdc_pos = kTRUE;
    
    	      Double_t fADCph = ((THcSignalHit*)scinPosADC->At(iphit))->GetData();
    	      fTOFPInfo[iphit].adcPh = fADCph;
    	      Double_t fPath = fPlanes[ip]->GetPosLeft() - fScinLongCoord;
    	      fTOFPInfo[iphit].path = fPath;
    
    	      // * Convert TDC value to time, do pulse height correction, correction for
    
    	      // * propogation of light thru scintillator, and offset.	      
    	      Double_t fTime = ((THcSignalHit*)scinPosTDC->At(iphit))->GetData() * fScinTdcToTime;
    	      fTime = fTime - ( fHodoPosPhcCoeff[fPIndex] * TMath::Sqrt( TMath::Max( 0. , 
    					( ( fADCph / fHodoPosMinPh[fPIndex] ) - 1 ) ) ) );
    	      fTime = fTime - ( fPath / fHodoVelLight[fPIndex] );
    	      fTOFPInfo[iphit].time = fTime;
    	      fTOFPInfo[iphit].scin_pos_time = fTime - fHodoPosTimeOffset[fPIndex];
    
    	      
    	    } // check for good pos TDC condition
    	    
    	    // ** Repeat for pmts on 'negative' side
    
    	    if ( ( ((THcSignalHit*)scinNegTDC->At(iphit))->GetData() > fScinTdcMin ) &&
    		 ( ((THcSignalHit*)scinNegTDC->At(iphit))->GetData() < fScinTdcMax ) &&
    		 ( fTOFPInfo[iphit].keep_neg ) ) { //
    
    	      
    	      // ** Calculate time for each tube with a good tdc. 'pos' side first.
    
    	      fTOFCalc[ihhit].good_tdc_neg = kTRUE;
    
    	      Double_t fADCph = ((THcSignalHit*)scinNegADC->At(iphit))->GetData();
    	      fTOFPInfo[iphit].adcPh = fADCph;
    	      //	      Double_t fPath = fPlanes[ip]->GetPosRight() - fScinLongCoord;
    	      Double_t fPath = fScinLongCoord - fPlanes[ip]->GetPosRight();
    	      fTOFPInfo[iphit].path = fPath;
    
    	      // * Convert TDC value to time, do pulse height correction, correction for
    	      // * propogation of light thru scintillator, and offset.
    
    	      Double_t fTime = ((THcSignalHit*)scinNegTDC->At(iphit))->GetData() * fScinTdcToTime;
    	      fTime = fTime - ( fHodoNegPhcCoeff[fPIndex] *
    			   TMath::Sqrt( TMath::Max( 0. , ( ( fADCph / fHodoNegMinPh[fPIndex] ) - 1 ) ) ) );
    	      fTime = fTime - ( fPath / fHodoVelLight[fPIndex] );
    	      fTOFPInfo[iphit].time = fTime;
    	      fTOFPInfo[iphit].scin_neg_time = fTime - fHodoNegTimeOffset[fPIndex];
          
    
    	    } // check for good neg TDC condition
    	    
    	    // ** Calculate ave time for scin and error.
    
    	    if ( fTOFCalc[ihhit].good_tdc_pos ){
    	      if ( fTOFCalc[ihhit].good_tdc_neg ){	
    
    		fTOFCalc[ihhit].scin_time  = ( fTOFPInfo[iphit].scin_pos_time + 
    					       fTOFPInfo[iphit].scin_neg_time ) / 2.;
    		fTOFCalc[ihhit].scin_sigma = TMath::Sqrt( fHodoPosSigma[fPIndex] * fHodoPosSigma[fPIndex] + 
    							  fHodoNegSigma[fPIndex] * fHodoNegSigma[fPIndex] )/2.;
    
    		fTOFCalc[ihhit].good_scin_time = kTRUE;
    
    		fTOFCalc[ihhit].scin_time = fTOFPInfo[iphit].scin_pos_time;
    
    		fTOFCalc[ihhit].scin_sigma = fHodoPosSigma[fPIndex];
    
    		fTOFCalc[ihhit].good_scin_time = kTRUE;
    
    	      if ( fTOFCalc[ihhit].good_tdc_neg ){
    		fTOFCalc[ihhit].scin_time = fTOFPInfo[iphit].scin_neg_time;
    
    		fTOFCalc[ihhit].scin_sigma = fHodoNegSigma[fPIndex];
    
    		fTOFCalc[ihhit].good_scin_time = 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[ihhit].good_scin_time ){
    
    	      // scin_time_fp doesn't need to be an array
    	      Double_t scin_time_fp = fTOFCalc[ihhit].scin_time -
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	       	( fPlanes[ip]->GetZpos() + ( paddle % 2 ) * fPlanes[ip]->GetDzpos() ) /
    	       	( 29.979 * betaP ) *
    
    	       	TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() +
    	       		     theTrack->GetPhi() * theTrack->GetPhi() );
    
    	      fSumfpTime = fSumfpTime + scin_time_fp;
    
    Zafar's avatar
    Zafar committed
    	      fNfpTime ++;
    
    	      fSumPlaneTime[ip] = fSumPlaneTime[ip] + scin_time_fp;
    
    	      fNPlaneTime[ip] ++;
    	      fNScinHit[itrack] ++;
    
    	      if ( ( fTOFCalc[ihhit].good_tdc_pos ) && ( fTOFCalc[ihhit].good_tdc_neg ) ){
    	      	fNPmtHit[itrack] = fNPmtHit[itrack] + 2;
    
    Zafar's avatar
    Zafar committed
    	      }
    	      else {
    
    	      	fNPmtHit[itrack] = fNPmtHit[itrack] + 1;
    	      }
    
    Zafar's avatar
    Zafar committed
    
    
    	      // --------------------------------------------------------------------------------------------
    
    	      if ( fTOFCalc[ihhit].good_tdc_pos ){
    		if ( fTOFCalc[ihhit].good_tdc_neg ){
    		  fdEdX[itrack][fNScinHit[itrack]-1]=
    		    TMath::Sqrt( TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(iphit))->GetData() *
                                                     ((THcSignalHit*)scinNegADC->At(iphit))->GetData() ) );
    
    Zafar's avatar
    Zafar committed
    		}
    		else{
    
    		  fdEdX[itrack][fNScinHit[itrack]-1]=
    		    TMath::Max( 0., ((THcSignalHit*)scinPosADC->At(iphit))->GetData() );
    
    Zafar's avatar
    Zafar committed
    	       	}
    	      }
    	      else{
    
    		if ( fTOFCalc[ihhit].good_tdc_neg ){
    		  fdEdX[itrack][fNScinHit[itrack]-1]=
    		    TMath::Max( 0., ((THcSignalHit*)scinNegADC->At(iphit))->GetData() );
    
    Zafar's avatar
    Zafar committed
    		}
    		else{
    
    		  fdEdX[itrack][fNScinHit[itrack]-1]=0.0;
    
    Zafar's avatar
    Zafar committed
    		}
    	      }
    
    	      // --------------------------------------------------------------------------------------------
    
    
    
    	    } // time at focal plane condition
    	  } // on track else condition
    	  
    	  // ** See if there are any good time measurements in the plane.
    
    	  if ( fTOFCalc[ihhit].good_scin_time ){
    
    	    fGoodPlaneTime[ip] = kTRUE;
    
    
    	  // Can this be done after looping over hits and planes?
    
    	  if ( fGoodPlaneTime[2] )	theTrack->SetGoodPlane3( 1 );
    
    Zafar's avatar
    Zafar committed
    	  if ( !fGoodPlaneTime[2] )	theTrack->SetGoodPlane3( 0 );
    
    	  if ( fGoodPlaneTime[3] )	theTrack->SetGoodPlane4( 1 );
    
    Zafar's avatar
    Zafar committed
    	  if ( !fGoodPlaneTime[3] )	theTrack->SetGoodPlane4( 0 );
    
    	} // Second loop over hits of a scintillator plane ends here
          } // Loop over scintillator planes ends here
    
    
    Zafar's avatar
    Zafar committed
          //------------------------------------------------------------------------------
          //------------------------------------------------------------------------------
          //------------------------------------------------------------------------------
          //------------------------------------------------------------------------------
          //------------------------------------------------------------------------------
    
          //------------------------------------------------------------------------------
          //------------------------------------------------------------------------------
          //------------------------------------------------------------------------------
    
    
          // * * Fit beta if there are enough time measurements (one upper, one lower)
    
          if ( ( ( fGoodPlaneTime[0] ) || ( fGoodPlaneTime[1] ) ) && 
    	   ( ( fGoodPlaneTime[2] ) || ( fGoodPlaneTime[3] ) ) ){	
    	
    
    	Double_t fSumW, fSumT, fSumZ, fSumZZ, fSumTZ;
    	Double_t fScinWeight, fTmp, fT0, fTmpDenom, fPathNorm, fZPosition, fTimeDif;
    
    	fSumW = 0.;	  fSumT = 0.;       fSumZ = 0.;     fSumZZ = 0.;	 fSumTZ = 0.;
    	fScinWeight = 0.;  fTmp = 0.;        fT0 = 0.;       fTmpDenom = 0.; 
    	fPathNorm = 0.;	  fZPosition = 0.;  fTimeDif = 0.;
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	for (Int_t ip = 0; ip < fNPlanes; ip++ ){
    
    Zafar's avatar
    Zafar committed
    	  fNScinHits[ip] = fPlanes[ip]->GetNScinHits();	  
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
    
    	    if ( fTOFCalc[ihhit].good_scin_time ) {
    
    	      fScinWeight = 1 / ( fTOFCalc[ihhit].scin_sigma * fTOFCalc[ihhit].scin_sigma );
    	      fZPosition = ( fPlanes[ip]->GetZpos() + ( fTOFCalc[ihhit].hit_paddle % 2 ) * 
    			     fPlanes[ip]->GetDzpos() );
    
    	      fSumW  += fScinWeight;
    	      fSumT  += fScinWeight * fTOFCalc[ihhit].scin_time;
    	      fSumZ  += fScinWeight * fZPosition;
    	      fSumZZ += fScinWeight * ( fZPosition * fZPosition );
    	      fSumTZ += fScinWeight * fZPosition * fTOFCalc[ihhit].scin_time;
    
    	      	      
    	    } // condition of good scin time
    
    	  } // loop over hits of plane
    	} // loop over planes
    
    	
    	fTmp = fSumW * fSumZZ - fSumZ * fSumZ ;
    	fT0 = ( fSumT * fSumZZ - fSumZ * fSumTZ ) / fTmp ;
    	fTmpDenom = fSumW * fSumTZ - fSumZ * fSumT;
    	
    	if ( TMath::Abs( fTmpDenom ) > ( 1 / 10000000000.0 ) ) {
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  beta = fTmp / fTmpDenom;
    	  betaChiSq = 0.;	  
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  for (Int_t ip = 0; ip < fNPlanes; ip++ ){                           // Loop over planes
    
    Zafar's avatar
    Zafar committed
    	    fNScinHits[ip] = fPlanes[ip]->GetNScinHits();	  
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	    for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){                    // Loop over hits of a plane
    
    	      if ( fTOFCalc[ihhit].good_scin_time ){
    
    		fZPosition = ( fPlanes[ip]->GetZpos() + ( fTOFCalc[ihhit].hit_paddle % 2 ) * 
    			       fPlanes[ip]->GetDzpos() );
    		fTimeDif = ( fTOFCalc[ihhit].scin_time - fT0 );		
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    		betaChiSq += ( ( fZPosition / beta - fTimeDif ) * 
    				( fZPosition / beta - fTimeDif ) )  / 
    
    		              ( fTOFCalc[ihhit].scin_sigma * fTOFCalc[ihhit].scin_sigma );
    
    		
    	      } // condition for good scin time
    
    	    } // loop over hits of a plane
    	  } // loop over planes
    	  
    
    	  fPathNorm = TMath::Sqrt( 1. + theTrack->GetTheta() * theTrack->GetTheta() + 
    				       theTrack->GetPhi()   * theTrack->GetPhi() );
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  beta = beta / fPathNorm;
    	  beta = beta / 29.979;    // velocity / c	  
    
    	  
    	}  // condition for fTmpDenom	
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  beta = 0.;
    	  betaChiSq = -2.;
    
    	} // else condition for fTmpDenom
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	beta = 0.;
    	betaChiSq = -1;
    
    Zafar's avatar
    Zafar committed
          if ( fNfpTime != 0 ){
    
          	fTimeAtFP[itrack] = ( fSumfpTime / fNfpTime ); 
    
    Zafar's avatar
    Zafar committed
          }
    
          //
          // ---------------------------------------------------------------------------
    
                
          Double_t fFPTimeSum=0.0;
          Int_t fNfpTimeSum=0;
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          for (Int_t ip = 0; ip < fNPlanes; ip++ ){
    
    	if ( fNPlaneTime[ip] != 0 ){
    	  fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
    
    	  fFPTimeSum += fSumPlaneTime[ip];
    	  fNfpTimeSum += fNPlaneTime[ip];
    
    	  fFPTime[ip] = 1000. * ( ip + 1 );
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          Double_t fptime = fFPTimeSum/fNfpTimeSum;
    
          // This can't be right.  Plus if there are no hits, then
          // it is undefined.
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          theTrack->SetFPTime(fptime);
          theTrack->SetBeta(beta);
          theTrack->SetBetaChi2( betaChiSq );
    
          theTrack->SetNPMT(fNPmtHit[itrack]);
          theTrack->SetFPTime( fTimeAtFP[itrack]);
    
    Zafar's avatar
    Zafar committed
    
    
        } // 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;
    
        scinPosTDC = fPlanes[ip]->GetPosTDC();
        scinNegTDC = fPlanes[ip]->GetNegTDC();
    
        fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
        for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
          Int_t paddlePos = ((THcSignalHit*)scinPosTDC->At(iphit))->GetPaddleNumber()-1;
          Int_t paddleNeg = ((THcSignalHit*)scinPosTDC->At(iphit))->GetPaddleNumber()-1;
          if ( paddlePos != paddleNeg )
    	return -1;
    
          fScinHitPaddle[ip][paddlePos] = 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.
    
      fBestXpScin = 100.0;
      fBestYpScin = 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 ) ){
    
    	fSlope = TMath::Abs(ipaddle - ipaddle2);
    
    	if ( fSlope < fBestXpScin ) {
    	  fBestXpScin = fSlope;
    
    
        }  // 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 ) ){
    
    	fSlope = TMath::Abs(ipaddle - ipaddle2);
    
    	if ( fSlope < fBestYpScin ) {
    	  fBestYpScin = fSlope;	
    
        }  // 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;
      }
    
      const char* apparatus_name = "H";
      THcHallCSpectrometer *app = static_cast<THcHallCSpectrometer*>(GetApparatus());
    
      if( ( strcmp(app->GetName(), apparatus_name) == 0 ) ) {
    
        if ( ( fGoodScinHits == 1 ) && ( fShower->GetNormETot() > 0.7 ) &&
    	 ( fChern->GetCerNPE() > 2.0 ) )
          fScinShould ++;
    
        if ( ( fGoodScinHits == 1 ) && ( fShower->GetNormETot() > 0.7 ) &&
         	 ( fChern->GetCerNPE() > 2.0 ) && ( tracks.GetLast() + 1 > 0 ) ) {
          fScinDid ++;
          fGood_hits ++;
        }
      }
    
    //_____________________________________________________________________________
    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)
    ////////////////////////////////////////////////////////////////////////////////