THcShower.cxx

    //
    for (Int_t i=fNclust-1; i>-1; i--) {

      THcShowerCluster* cluster = *(fClusterList->begin()+i);

      Double_t dx = TMath::Abs( clX(cluster) - XTrFront );

      if (dx <= (0.5*BlockThick[0] + fSlop)) {
	fNtracks++;  // lumber of shower tracks (Consistent with engine)
	if (dx < deltaX) {
	  mclust = i;
	  deltaX = dx;
	}
      }
    }
  }

  //Debug output.

  if (fdbg_tracks_cal) {
    cout << "---------------------------------------------------------------\n";
    cout << "Debug output from THcShower::MatchCluster for "
	 << GetApparatus()->GetName() << endl;
    cout << " event = " << fEvent << endl;

    cout << "  Track at DC:"
	 << "  X = " << Track->GetX()
	 << "  Y = " << Track->GetY()
	 << "  Theta = " << Track->GetTheta()
	 << "  Phi = " << Track->GetPhi()
	 << endl;
    cout << "  Track at the front of Calorimeter:"
	 << "  X = " << XTrFront
	 << "  Y = " << YTrFront
	 << "  Pathl = " << pathl
	 << endl;
    if (fvTest)
      cout << "  Fiducial volume test: inFidVol = " << inFidVol << endl;

    cout << "  Matched cluster #" << mclust << ",  delatX= " << deltaX
    	 << endl;
    cout << "---------------------------------------------------------------\n";
  }

  return mclust;
}

//_____________________________________________________________________________
Float_t THcShower::GetShEnergy(THaTrack* Track) {

  // Get total energy deposited in the cluster matched to the given
  // spectrometer Track.

  // Track coordinates at front of the calorimeter, initialize out of
  // acceptance.
  Double_t Xtr = -100.;
  Double_t Ytr = -100.;

  // Associate a cluster to the track.

  Int_t mclust = MatchCluster(Track, Xtr, Ytr);

  // Coordinate corrected total energy deposition in the cluster.

  Float_t Etrk = 0.;
  if (mclust >= 0) {         // if there is a matched cluster

    // Get matched cluster.

    THcShowerCluster* cluster = *(fClusterList->begin()+mclust);

    // Correct track energy depositions for the impact coordinate.

    for (UInt_t ip=0; ip<fNLayers; ip++) {

      // Coordinate correction factors for positive and negative sides,
      // different for double PMT counters in the 1-st two layes and for
      // single PMT counters in the rear two layers.
      Float_t corpos = 1.;
      Float_t corneg = 1.;

      if (ip < fNegCols) {
	corpos = Ycor(Ytr,0);
	corneg = Ycor(Ytr,1);
      }
      else {
	corpos = Ycor(Ytr);
	corneg = 0.;
      }

      // cout << ip << " clust energy pos = " <<  clEplane(cluster,ip,0)<< " clust energy pos = " <<  clEplane(cluster,ip,1) << endl;
      Etrk += clEplane(cluster,ip,0) * corpos;
      Etrk += clEplane(cluster,ip,1) * corneg;

    }   //over planes

  }   //mclust>=0

  //Debug output.

  if (fdbg_tracks_cal) {
    cout << "---------------------------------------------------------------\n";
    cout << "Debug output from THcShower::GetShEnergy for "
	 << GetApparatus()->GetName() << endl;
    cout << " event = " << fEvent << endl;

    cout << "  Track at the calorimeter: X = " << Xtr << "  Y = " << Ytr;
    if (mclust >= 0)
      cout << ", matched cluster #" << mclust << "." << endl;
    else
      cout << ", no matched cluster found." << endl;

    cout << "  Coordinate corrected track energy = " << Etrk << "." << endl;
    cout << "---------------------------------------------------------------\n";
  }

  return Etrk;
}

//_____________________________________________________________________________
Int_t THcShower::FineProcess( TClonesArray& tracks )
{

  // Shower energy assignment to the spectrometer tracks.
  //

  Int_t Ntracks = tracks.GetLast()+1;   // Number of reconstructed tracks
      Double_t Xtr = -100.;
      Double_t Ytr = -100.;
  for (Int_t itrk=0; itrk<Ntracks; itrk++) {
    THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
    if (theTrack->GetIndex()==0) {
      fEtrack=theTrack->GetEnergy();
      fEtrackNorm=fEtrack/theTrack->GetP();
      Xtr = -100.;
      Ytr = -100.;               
      fNclustTrack = MatchCluster(theTrack, Xtr, Ytr);
      fXTrack=Xtr;
      fYTrack=Ytr;
      if (fNclustTrack>=0) {
	THcShowerCluster* cluster = *(fClusterList->begin()+fNclustTrack);
	fXclustTrack=clX(cluster);
        fYclustTrack=clY(cluster);
      }
      if (fHasArray) {
	fNclustArrayTrack = fArray->MatchCluster(theTrack,Xtr,Ytr);      
        if (fNclustArrayTrack>=0) {
          fXclustArrayTrack=fArray->GetClX();
          fYclustArrayTrack=fArray->GetClY();
          fSizeClustArray=fArray->GetClSize();
          fNblockHighEnergy=fArray->GetClMaxEnergyBlock();
          fXTrackArray=Xtr;
          fYTrackArray=Ytr;
        }
      }
    }
  }       //over tracks

  //Debug output.

  if (fdbg_tracks_cal) {
    cout << "---------------------------------------------------------------\n";
    cout << "Debug output from THcShower::FineProcess for "
	 << GetApparatus()->GetName() << endl;
    cout << " event = " << fEvent << endl;
    cout << "  Number of tracks = " << Ntracks << endl;
    if (fNclustTrack>=0) {
      cout << " matching cluster info " << endl;
      cout << " X info = " << fXclustTrack << " " << Xtr << " " << fXclustTrack-Xtr << endl;
      cout << " Y info = " << fYclustTrack << " " << Ytr << " " << fYclustTrack-Ytr << endl;
    }

    for (Int_t itrk=0; itrk<Ntracks; itrk++) {
      THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
      cout << "  Track " << itrk << ": "
	   << "  X = " << theTrack->GetX()
	   << "  Y = " << theTrack->GetY()
	   << "  Theta = " << theTrack->GetTheta()
	   << "  Phi = " << theTrack->GetPhi()
	   << "  Energy = " << theTrack->GetEnergy() 
	   << "  Energy/Ptrack = " <<  fEtrackNorm << endl;
    }

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

  return 0;
}

Double_t THcShower::GetNormETot( ){
  return fEtotNorm;
}

ClassImp(THcShower)
////////////////////////////////////////////////////////////////////////////////