Skip to content
Snippets Groups Projects
THcHodoscope.cxx 54.3 KiB
Newer Older
  • Learn to ignore specific revisions
  • 	    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;
    
    	      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
    
    	} // First loop over hits in a plane <---------
    
    	//-----------------------------------------------------------------------------------------------
    
    	//------------- First large loop over scintillator hits ends here --------------------
    
    	//-----------------------------------------------------------------------------------------------
    
          // 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];
    
    	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;
    
    	//         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) {
    	  if ( fTOFPInfo[ih].keep_pos ) { // 301
    	    fTOFCalc[ih].good_tdc_pos = kTRUE;
    
    	    fGoodFlags[itrack][ip][iphit].goodTdcPos = kTRUE;
    
    	    if(fDumpTOF && ntracks==1 && fGoodEventTOFCalib) {
                  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;
    
    	  if ( fTOFPInfo[ih].keep_neg ) { //
    	    fTOFCalc[ih].good_tdc_neg = kTRUE;
    
    	    fGoodFlags[itrack][ip][iphit].goodTdcNeg = kTRUE;
    
    	    if(fDumpTOF && ntracks==1 && fGoodEventTOFCalib) {
                 fDumpOut << fixed << setprecision(2);
                 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;
    
    
    	  // ** 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;
    
    	    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() );
    
    Zafar's avatar
    Zafar committed
    	      }
    
    	      if ( fTOFCalc[ih].good_tdc_neg ){
    
    		fdEdX[itrack][fNScinHit[itrack]-1]=
    		  TMath::Max( 0., hit->GetNegADC() );
    	      } else{
    		fdEdX[itrack][fNScinHit[itrack]-1]=0.0;
    
    Zafar's avatar
    Zafar committed
    	      }
    
    	    }
    	    // --------------------------------------------------------------------------------------------
    
    	  } // 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];
    
    	  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 );
    
         //
           //------------------------------------------------------------------------------
    
    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] ) ) ){
    
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	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() );
    
    	    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
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	Double_t tmp = sumW * sumZZ - sumZ * sumZ ;
    	Double_t t0 = ( sumT * sumZZ - sumZ * sumTZ ) / tmp ;
    	Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	if ( TMath::Abs( tmpDenom ) > ( 1 / 10000000000.0 ) ) {
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  beta = tmp / tmpDenom;
    
    	  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
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  beta = beta / pathNorm;
    
    	  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;
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          if ( nFPTime != 0 ){
    
          	timeAtFP[itrack] = ( sumFPTime / nFPTime );
    
    Zafar's avatar
    Zafar committed
          }
    
          //
          // ---------------------------------------------------------------------------
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          Double_t FPTimeSum=0.0;
          Int_t nFPTimeSum=0;
    
          for (Int_t ip = 0; ip < temp_planes; ip++ ){
    
    	if ( fNPlaneTime[ip] != 0 ){
    	  fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
    	  FPTimeSum += fSumPlaneTime[ip];
    	  nFPTimeSum += fNPlaneTime[ip];
    
    	  fFPTime[ip] = 1000. * ( ip + 1 );
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          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);
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          theTrack->SetFPTime(fptime);
          theTrack->SetBeta(beta);
          theTrack->SetBetaChi2( betaChiSq );
    
    Stephen A. Wood's avatar
    Stephen A. Wood committed
          theTrack->SetNPMT(nPmtHit[itrack]);
          theTrack->SetFPTime( timeAtFP[itrack]);
    
    Zafar's avatar
    Zafar committed
    
    
        } // Main loop over tracks ends here.
    
      } // If condition for at least one track
    
    
       if(fDumpTOF && ntracks==1 && fGoodEventTOFCalib) {
             fDumpOut << "0 "  << endl;
    
      //-----------------------------------------------------------------------
      //
      //   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 < temp_planes; 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;
    
        }
      }
    
      //  *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 < temp_planes; 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 ++;
    
    
    hallc-online's avatar
    hallc-online committed
        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 ++;
    
        } // Loop over  paddles
    
        fNClust[ip] = icount;
        icount = 0;
    
    
    hallc-online's avatar
    hallc-online committed
        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
    
        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.
    
    
    hallc-online's avatar
    hallc-online committed
      for (Int_t ip = 1; ip < temp_planes; ip +=2 ) {
    
        // Planes ip = 1 = 1Y
        // Planes ip = 3 = 2Y
        if (!fPlanes[ip]) return -1;
    
        icount = 0;
        if ( fScinHitPaddle[ip][0] == 1 )
          icount ++;
    
    
    hallc-online's avatar
    hallc-online committed
        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 ++;
    
        } // Loop over Y paddles
    
        fNClust[ip] = icount;
        icount = 0;
    
    
    hallc-online's avatar
    hallc-online committed
        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
    
        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);
      }
    
    
    hallc-online's avatar
    hallc-online committed
      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
    
    hallc-online's avatar
    hallc-online committed
      // * 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
      for (Int_t itrk=0; itrk<Ntracks; itrk++) {
        THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
        if (theTrack->GetIndex()==0) {
        fBeta=theTrack->GetBeta();
        }
      }       //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() {
      return fPathLengthCentral;
    }
    
    ClassImp(THcHodoscope)
    ////////////////////////////////////////////////////////////////////////////////