THcShowerArray.cxx

      piccolumn = 0;
    }
  }
#endif

  //Debug output.


  return(ihit);
}
//_____________________________________________________________________________
Int_t THcShowerArray::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
{
  // Extract data for this plane from hit list and accumulate in
  // arrays for subsequent pedestal calculations.

  Int_t nrawhits = rawhits->GetLast()+1;

  Int_t ihit = nexthit;

  while(ihit < nrawhits) {

    THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);

    if(hit->fPlane != fLayerNum) {
      break;
    }

    Int_t element = hit->fCounter - 1; // Should check if in range

		// Should check that counter # is in range
		Int_t adc = 0;
		if (fUsingFADC) {
			adc = hit->GetRawAdcHitPos().GetData(
				fPedSampLow, fPedSampHigh, fDataSampLow, fDataSampHigh
			);
		}
		else {
			adc = hit->GetData(0);
		}

    if(adc <= fPedLimit[element]) {
      fPedSum[element] += adc;
      fPedSum2[element] += adc*adc;
      fPedCount[element]++;
      if(fPedCount[element] == fMinPeds/5) {
	fPedLimit[element] = 100 + fPedSum[element]/fPedCount[element];
      }
    }
    ihit++;
  }
  fNPedestalEvents++;

  // Debug output.

  if ( ((THcShower*) GetParent())->fdbg_raw_cal ) {

    cout << "---------------------------------------------------------------\n";
    cout << "Debug output from THcShowerArray::AcculatePedestals for "
    	 << GetParent()->GetPrefix() << ":" << endl;

    cout << "Processed hit list for plane " << GetName() << ":\n";

    for (Int_t ih=nexthit; ih<nrawhits; ih++) {

      THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ih);

      if(hit->fPlane != fLayerNum) {
	break;
      }

			Int_t adc = 0;
			if (fUsingFADC) {
				adc = hit->GetRawAdcHitPos().GetData(
					fPedSampLow, fPedSampHigh, fDataSampLow, fDataSampHigh
				);
			}
			else {
				adc = hit->GetData(0);
			}

      cout << "  hit " << ih << ":"
	   << "  plane =  " << hit->fPlane
	   << "  counter = " << hit->fCounter
	   << "  ADC = " << adc
	   << endl;
    }

    cout << "---------------------------------------------------------------\n";

  }

  return(ihit);
}

//_____________________________________________________________________________
void THcShowerArray::CalculatePedestals( )
{
  // Use the accumulated pedestal data to calculate pedestals.

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

    fPed[i] = ((Float_t) fPedSum[i]) / TMath::Max(1, fPedCount[i]);
    fSig[i] = sqrt(((Float_t)fPedSum2[i])/TMath::Max(1, fPedCount[i])
		   - fPed[i]*fPed[i]);
    fThresh[i] = fPed[i] + TMath::Min(50., TMath::Max(10., 3.*fSig[i]));

  }

  // Debug output.

  if ( ((THcShower*) GetParent())->fdbg_raw_cal ) {

    cout << "---------------------------------------------------------------\n";
    cout << "Debug output from THcShowerArray::CalculatePedestals for "
    	 << GetParent()->GetPrefix() << ":" << endl;

    cout << "  ADC pedestals and thresholds for calorimeter plane "
	 << GetName() << endl;
    for(Int_t i=0; i<fNelem;i++) {
      cout << "  element " << i << ": "
	   << "  Pedestal = " << fPed[i]
	   << "  threshold = " << fThresh[i]
	   << endl;
    }
    cout << "---------------------------------------------------------------\n";

  }

}
//_____________________________________________________________________________
void THcShowerArray::InitializePedestals( )
{
  fNPedestalEvents = 0;

  fPedSum = new Int_t [fNelem];
  fPedSum2 = new Int_t [fNelem];
  fPedCount = new Int_t [fNelem];

  fSig = new Float_t [fNelem];
  fPed = new Float_t [fNelem];
  fThresh = new Float_t [fNelem];

  for(Int_t i=0;i<fNelem;i++) {
    fPedSum[i] = 0;
    fPedSum2[i] = 0;
    fPedCount[i] = 0;
  }
}

//------------------------------------------------------------------------------

// Fiducial volume limits.

Double_t THcShowerArray::fvXmin() {
  THcShower* fParent;
  fParent = (THcShower*) GetParent();
  return fXPos[0][0] - fXStep/2 + fParent->fvDelta;
}

Double_t THcShowerArray::fvYmax() {
  THcShower* fParent;
  fParent = (THcShower*) GetParent();
  return fYPos[0][0] + fYStep/2 - fParent->fvDelta;
}

Double_t THcShowerArray::fvXmax() {
  THcShower* fParent;
  fParent = (THcShower*) GetParent();
  return fXPos[fNRows-1][fNColumns-1] + fXStep/2 - fParent->fvDelta;
}

Double_t THcShowerArray::fvYmin() {
  THcShower* fParent;
  fParent = (THcShower*) GetParent();
  return fYPos[fNRows-1][fNColumns-1] - fYStep/2 + fParent->fvDelta;
}
Double_t THcShowerArray::clMaxEnergyBlock(THcShowerCluster* cluster) {
  Double_t max_energy=-1.;
  Double_t max_block=-1.;
  for (THcShowerClusterIt it=(*cluster).begin(); it!=(*cluster).end(); ++it) {
    if ( (**it).hitE() > max_energy )   {
       max_energy = (**it).hitE();
       max_block = ((**it).hitColumn())*fNRows + (**it).hitRow()+1;
    }
  }
  return max_block;
}

//_____________________________________________________________________________
Int_t THcShowerArray::AccumulateStat(TClonesArray& tracks )
{
  // Accumumate statistics for efficiency calculations.
  //
  // Choose electron events in gas Cherenkov with good Chisq of the best track.
  // Project best track to Array,
  // calculate module number for the track,
  // accrue number of tracks for the module,
  // accrue number of hits for the module, if it is hit.
  // Accrue total numbers of tracks and hits for Array.

  THaTrack* BestTrack = static_cast<THaTrack*>( tracks[0]);
  if (BestTrack->GetChi2()/BestTrack->GetNDoF() > fStatMaxChi2) return 0;

  THcHallCSpectrometer *app=dynamic_cast<THcHallCSpectrometer*>(GetApparatus());

  THaDetector* detc;
  if (GetParent()->GetPrefix()[0] == 'P')
    detc = app->GetDetector("hgcer");
  else
    detc = app->GetDetector("cer");
  
  THcCherenkov* hgcer = dynamic_cast<THcCherenkov*>(detc);
  if (!hgcer) {
    cout << "****** THcShowerArray::AccumulateStat: HGCer not found! ******"
	 << endl;
    return 0;
  }

  if (hgcer->GetCerNPE() < fStatCerMin) return 0;
  
  Double_t XTrk = kBig;
  Double_t YTrk = kBig;
  Double_t pathl = kBig;

  // Track interception with Array. The coordinates are in the calorimeter's
  // local system.

  fOrigin = GetOrigin();
  THcShower* fParent = (THcShower*) GetParent();
  fParent->CalcTrackIntercept(BestTrack, pathl, XTrk, YTrk);

  // Transform coordiantes to the spectrometer's coordinate system.
  XTrk += GetOrigin().X();
  YTrk += GetOrigin().Y();
						     
  for (Int_t i=0; i<fNelem; i++) {

    Int_t row = i%fNRows;
    Int_t col =i/fNRows;

    if (TMath::Abs(XTrk - fXPos[row][col]) < fStatSlop &&
	TMath::Abs(YTrk - fYPos[row][col]) < fStatSlop) {

      fStatNumTrk.at(i)++;
      fTotStatNumTrk++;
      
      if (fGoodAdcPulseInt.at(i) > 0.) {
	fStatNumHit.at(i)++;
	fTotStatNumHit++;
      }
      
    }
    
  }

  if ( ((THcShower*) GetParent())->fdbg_tracks_cal ) {
    cout << "---------------------------------------------------------------\n";
    cout << "THcShowerArray::AccumulateStat:" << endl;
    cout << "   Chi2/NDF = " <<BestTrack->GetChi2()/BestTrack->GetNDoF() <<endl;
    cout << "   HGCER Npe = " << hgcer->GetCerNPE() << endl;
    cout << "   XTrk, YTrk = " << XTrk << "  " << YTrk << endl;						     
    for (Int_t i=0; i<fNelem; i++) {
      
      Int_t row = i%fNRows;
      Int_t col =i/fNRows;

      if (TMath::Abs(XTrk - fXPos[row][col]) < fStatSlop &&
	  TMath::Abs(YTrk - fYPos[row][col]) < fStatSlop) {

	cout << "   Module " << i << ", X=" << fXPos[i/fNRows][i%fNRows]
	     << ", Y=" << fYPos[i/fNRows][i%fNRows] << " matches track" << endl;

	if (fGoodAdcPulseInt.at(i) > 0.)
	  cout << "   PulseIntegral = " << fGoodAdcPulseInt.at(i) << endl;
      }
    }
      
    cout << "---------------------------------------------------------------\n";
    //    getchar();
  }
  
  return 1;
}