THcShower.cxx

    XTrBack += GetOrigin().X();   // from local coord. system
    YTrBack += GetOrigin().Y();   // to the spectrometer system

    inFidVol = (XTrFront <= fvXmax) && (XTrFront >= fvXmin) &&
               (YTrFront <= fvYmax) && (YTrFront >= fvYmin) &&
               (XTrBack <= fvXmax) && (XTrBack >= fvXmin) &&
               (YTrBack <= fvYmax) && (YTrBack >= fvYmin);

  }

  // Match a cluster to the track.

  Int_t mclust = -1;    // The match cluster #, initialize with a bogus value.
  Double_t deltaX = kBig;   // Track to cluster distance

  if (inFidVol) {

    // Since hits and clusters are in reverse order (with respect to Engine),
    // search backwards to be consistent with Engine.
    //
    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, UInt_t NLayers, UInt_t L0) {

  // Get part of energy deposited in the cluster matched to the given
  // spectrometer Track, limited by range of layers from L0 to L0+NLayers-1.

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

    char prefix = GetApparatus()->GetName()[0];

    // Correct track energy depositions for the impact coordinate.

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

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

      if (ip < fNegCols) {
	if (prefix == 'H') {              //HMS 1-st 2 layers
	  corpos = Ycor(Ytr,0);
	  corneg = Ycor(Ytr,1);
	}
	else {                            //SHMS Preshower
	  corpos = YcorPr(Ytr,0);
	  corneg = YcorPr(Ytr,1);
	}
      }
      else {
	corpos = Ycor(Ytr);
	corneg = 0.;
      }

      // cout << ip << " clust energy pos = " <<  clEplane(cluster,ip,0)<< " clust energy neg = " <<  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 << "  Layers from " << L0+1 << " to " << L0+NLayers << ".\n";
    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();
      fEPRtrack=GetShEnergy(theTrack,1);
      fEPRtrackNorm=fEPRtrack/theTrack->GetP();
      fETotTrackNorm=fEtot/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

  if (Ntracks>0) {
    for(UInt_t ip=0;ip<fNLayers;ip++) {
      fPlanes[ip]-> AccumulateStat(tracks);
    }
    if(fHasArray) fArray->AccumulateStat(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)
////////////////////////////////////////////////////////////////////////////////