////////////////////////////////////////////////////////////////////////
// $Id: AlgMSTClusterList.cxx,v 1.24 2002/08/24 23:45:53 miyagawa Exp $
//
// AlgMSTClusterList
//
// Begin_Html<img src="../../pedestrians.gif" align=center>
// <a href="../source_warning.html">Warning for beginners</a>.<br> 
//
// This is an Algorithm class to create a list of CandMSTClusters from a
// CandDigiPairs contained in a CandSliceList.  It first uses a minimum
// spanning tree algorithm to cluster digit pairs together by proximity.
// It then connects clusters on adjacent planes to account for the gap
// between super-modules.  It then connects clusters across the central
// coil hole.
//
// Author:  P.S. Miyagawa 9/2000
//
// Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
// <a href="../ALG_Classes.html"> Algorithm Classes</a> (part of
// <a href="../index.html">The MINOS Class User Guide</a>)End_Html
////////////////////////////////////////////////////////////////////////

#include <cassert>

#include "Algorithm/AlgConfig.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Candidate/CandContext.h"
#include "MessageService/MsgService.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "BubbleSpeak/AlgMSTClusterList.h"
#include "BubbleSpeak/CandMSTCluster.h"
#include "BubbleSpeak/CandMSTClusterHandle.h"
#include "BubbleSpeak/CandDigiPairHandle.h"
#include "BubbleSpeak/ClusterBox.h"

ClassImp(AlgMSTClusterList)

//......................................................................

CVSID("$Id: AlgMSTClusterList.cxx,v 1.24 2002/08/24 23:45:53 miyagawa Exp $");

//......................................................................

AlgMSTClusterList::AlgMSTClusterList()
{
//
//  Purpose:    Default constructor.
//
//  Arguments:  n/a
//
//  Return:     n/a
//
}

//......................................................................

AlgMSTClusterList::~AlgMSTClusterList()
{
//
//  Purpose:    Default destructor.
//
//  Arguments:  n/a
//
//  Return:     n/a
//
}

//......................................................................

Bool_t AlgMSTClusterList::IsCloseZ(const ClusterBox &b1,
                               const ClusterBox &b2, Int_t sep) const
{
//
//  Purpose:  Determine whether two boxes are close in plane number.
//            Used primarily to bridge the gap between super-modules.
//
//  Arguments:
//    b1,b2     in    The boxes to compare.
//    sep       in    Max separation in plane number between box ends.
//
//  Return:   kTRUE   if box ends are within sep
//            kFALSE  otherwise
//

   Int_t diff1 = b1.GetPlo() - b2.GetPhi();
   Int_t diff2 = b2.GetPlo() - b1.GetPhi();
   return (diff1 <= sep && diff1 > 0) || (diff2 <= sep && diff2 > 0);
}

//......................................................................

Bool_t AlgMSTClusterList::IsInCross(const ClusterBox &b1,
          const ClusterBox &b2, Float_t hol, Float_t dist) const
{
//
//  Purpose:  Determine whether the two boxes are an inward crossing.
//
//  Arguments:
//    b1,b2     in    The boxes to compare.
//    hol       in    Max radius at inward crossing end.
//    dist      in    Max separation between inward ends of boxes
//
//  Return:   kTRUE   if the clusters form an inward crossing
//            kFALSE  if not
//
//  Notes:    Assumes that the digit pair clusters are filled with
//            strips of opposite orientation.
//

// Check whether inward end of clusters overlap.
   if (TMath::Abs(b1.GetYhi() - b2.GetYhi()) < dist) {

// Get transverse positions of inward ends of each cluster.
      Float_t xmin = b1.GetXin();
      Float_t ymin = b2.GetXin();

// Check radius.
      Float_t rad2 = xmin * xmin + ymin * ymin;
      Float_t hol2 = hol * hol;
      return (rad2 < hol2);
   }

// No overlap.
   else return kFALSE;
}

//......................................................................

Bool_t AlgMSTClusterList::IsOutCross(const ClusterBox &b1,
          const ClusterBox &b2, Float_t hol, Float_t dist) const
{
//
//  Purpose:  Determine whether the two boxes are an outward crossing.
//
//  Arguments:
//    b1,b2     in    The boxes to compare.
//    hol       in    Max radius at outward crossing end.
//    dist      in    Max separation between outward ends of boxes
//
//  Return:   kTRUE   if the clusters form an outward crossing
//            kFALSE  if not
//
//  Notes:    Assumes that the digit pair clusters are filled with
//            strips of opposite orientation.
//

// Check whether outward end of clusters overlap.
   if (TMath::Abs(b1.GetYlo() - b2.GetYlo()) < dist) {

// Get transverse positions of outward ends of each cluster.
      Float_t xmin = b1.GetXout();
      Float_t ymin = b2.GetXout();

// Check radius.
      Float_t rad2 = xmin * xmin + ymin * ymin;
      Float_t hol2 = hol * hol;
      return (rad2 < hol2);
   }

// No overlap.
   else return kFALSE;
}

//......................................................................

Bool_t AlgMSTClusterList::IsOverlap(const ClusterBox &b1,
          const ClusterBox &b2, Float_t dist) const
{
//
//  Purpose:  Determine whether the digit pairs in two boxes overlap.
//
//  Arguments:
//    b1,b2     in    The boxes to compare.
//    dist      in    Max separation between digit pairs for boxes to
//                    overlap.
//
//  Return:   kTRUE   if any digit pairs of this box are within dist of
//                    any digit pairs of the other box
//            kFALSE  otherwise
//

// Check whether envelopes overlap.
   if (((b1.GetXlo() - b2.GetXhi()) > dist)
       || ((b2.GetXlo() - b1.GetXhi()) > dist)
       || ((b1.GetYlo() - b2.GetYhi()) > dist)
       || ((b2.GetYlo() - b1.GetYhi()) > dist)) {
      return kFALSE;   // No envelope overlap.
   }

// Check whether any digit pairs overlap.
   TIter b1Itr(&b1);
   TIter b2Itr(&b2);
   while (CandDigiPairHandle *b1Chh =
             dynamic_cast<CandDigiPairHandle *>(b1Itr())) {
      Float_t dist2 = dist * dist;
      while (CandDigiPairHandle *b2Chh =
                dynamic_cast<CandDigiPairHandle *>(b2Itr())) {
         Float_t sepX = b1Chh->GetTPos() - b2Chh->GetTPos();
         Float_t sepY = b1Chh->GetZPos() - b2Chh->GetZPos();
         Float_t sep2  = sepX * sepX + sepY * sepY;
         if (sep2 <= dist2) return kTRUE;   // Overlap found.
      }
      b2Itr.Reset();
   }
   return kFALSE;   // No overlap.
}

//......................................................................

void AlgMSTClusterList::RunAlg(AlgConfig &ac, CandHandle &ch,
                               CandContext &cx)
{
//
//  Purpose:  Group together CandDigiPairs into MSTClusters that are
//            stored in a new CandMSTClusterList.
//
//  Argument:
//    ac        in    AlgConfig containing cluster parameters.
//    ch        in    Handle to the new CandMSTClusterList to fill.
//    cx        in    CandContext containing the CandSliceList to
//                    cluster.
//
//  Return:   n/a
//

   MSG("BubAlg", Msg::kVerbose)
      << "Starting AlgMSTClusterList::RunAlg()" << endl;

// Check for CandSliceListHandle input.
   assert(cx.GetDataIn());
   assert(cx.GetDataIn()->InheritsFrom("CandSliceListHandle"));
   const CandSliceListHandle *cslh =
               dynamic_cast<const CandSliceListHandle*>(cx.GetDataIn());

// Iterate over CandSlices.
   TIter cshItr(cslh->GetDaughterIterator());
   while (CandSliceHandle *csh =
             dynamic_cast<CandSliceHandle *>(cshItr())) {
      RunAlgOnSlice(ac, ch, cx, csh);
   }
}

//......................................................................

void AlgMSTClusterList::RunAlgOnSlice(AlgConfig &ac, CandHandle &ch,
                                  CandContext &cx, CandSliceHandle *csh)
{
//
//  Purpose:  Group together CandDigiPairs into MSTClusters that are
//            stored in a new CandMSTClusterList.
//
//  Argument:
//    ac        in    AlgConfig containing cluster parameters.
//    ch        in    Handle to the new CandMSTClusterList to fill.
//    cx        in    CandContext containing information needed for call
//                    to AlgMSTCluster.
//    csh       in    CandSlice containing CandDigiPairs to group into a
//                    CandMSTCluster.
//
//  Return:   n/a
//

   MSG("BubAlg", Msg::kVerbose)
      << "Starting AlgMSTClusterList::RunAlgOnSlice()" << endl;

// Save config parameters.
   Double_t distmax = ac.GetDouble("DistMax");
   Double_t holemax = ac.GetDouble("HoleMax");
   Int_t planemax = ac.GetInt("PlaneMax");

// Initialize cluster boxes.
   TObjArray *cbayU = new TObjArray();
   TObjArray *cbayV = new TObjArray();
   TObjArray *arrayA = 0;
   TObjArray *arrayB = 0;
   TObjArray *shield = 0;
   assert(csh);
   TIter chhItr(csh->GetDaughterIterator());
   while (CandDigiPairHandle *chh =
             dynamic_cast<CandDigiPairHandle *>(chhItr())) {

      // Check for digits from veto shield.
      if (chh->GetStripEndId().IsVetoShield()) {
         if (!shield) shield = new TObjArray();
         shield->Add(chh);
         continue;
      }

      // Separate by strip orientation.
      switch (chh->GetPlaneView()) {
      case PlaneView::kU:
         cbayU->Add(new ClusterBox(chh));
         break;
      case PlaneView::kV:
         cbayV->Add(new ClusterBox(chh));
         break;
      case PlaneView::kA:
         if (!arrayA) arrayA = new TObjArray();
         arrayA->Add(chh);
         break;
      case PlaneView::kB:
         if (!arrayB) arrayB = new TObjArray();
         arrayB->Add(chh);
         break;
      default:
         MSG("BubAlg", Msg::kWarning)
            << "Invalid plane orientation." << endl;
      }
   }

// MSTCluster using MST algorithm.
   RunMSTAlg(distmax, cbayU);
   RunMSTAlg(distmax, cbayV);

// Join clusters across super-modules.
   RunCloseZAlg(planemax, cbayU);
   RunCloseZAlg(planemax, cbayV);

// Join clusters across the coil hole.
   RunHoleCrossAlg(holemax, distmax, cbayU, cbayV);

// General setup for creating new CandMSTClusters.
   MSG("BubAlg", Msg::kVerbose)
      << "AlgFactory &af = AlgFactory::GetInstance();" << endl;
   AlgFactory &af = AlgFactory::GetInstance();

   MSG("BubAlg", Msg::kVerbose)
      << "AlgHandle ah = af.GetAlgHandle(\"AlgMSTCluster\","
      << "\"default\");" << endl;
   AlgHandle ah = af.GetAlgHandle("AlgMSTCluster", "default");

   MSG("BubAlg", Msg::kVerbose)
      << "Create CandContext instance." << endl;
   CandContext cxx(this, cx.GetMom());
   cxx.SetCandRecord(cx.GetCandRecord());

// Iterate over ClusterBoxes to create CandMSTClusters.
   // U-view clusters
   TIter cbUItr(cbayU);
   while (ClusterBox *cb = dynamic_cast<ClusterBox *>(cbUItr())) {
      cxx.SetDataIn(cb);
      CandMSTClusterHandle clh = CandMSTCluster::MakeCandidate(ah, cxx);
      clh.SetCandSlice(csh);
      ch.AddDaughterLink(clh);
   }
   cbayU->Delete();
   delete cbayU;
   cbayU = 0;

   // V-view clusters
   TIter cbVItr(cbayV);
   while (ClusterBox *cb = dynamic_cast<ClusterBox *>(cbVItr())) {
      cxx.SetDataIn(cb);
      CandMSTClusterHandle clh = CandMSTCluster::MakeCandidate(ah, cxx);
      clh.SetCandSlice(csh);
      ch.AddDaughterLink(clh);
   }
   cbayV->Delete();
   delete cbayV;
   cbayV = 0;

   // A-view clusters
   if (arrayA) {
      cxx.SetDataIn(arrayA);
      CandMSTClusterHandle clh = CandMSTCluster::MakeCandidate(ah, cxx);
      clh.SetCandSlice(csh);
      ch.AddDaughterLink(clh);
      delete arrayA;
      arrayA = 0;
   }

   // B-view clusters
   if (arrayB) {
      cxx.SetDataIn(arrayB);
      CandMSTClusterHandle clh = CandMSTCluster::MakeCandidate(ah, cxx);
      clh.SetCandSlice(csh);
      ch.AddDaughterLink(clh);
      delete arrayB;
      arrayB = 0;
   }

   // Veto shield clusters
   if (shield) {
      cxx.SetDataIn(shield);
      CandMSTClusterHandle clh = CandMSTCluster::MakeCandidate(ah, cxx);
      clh.SetCandSlice(csh);
      ch.AddDaughterLink(clh);
      delete shield;
      shield = 0;
   }
}

//......................................................................

void AlgMSTClusterList::RunCloseZAlg(Int_t sep, TObjArray *&toay)
{
//
//  Purpose:  Group clusters together across the gap between super-
//            modules using plane number.
//
//  Arguments:
//    sep     in      Maximum plane number separation between clusters.
//    toay    in/out  Array with initial cluster boxes. It is replaced
//                    by an array with final cluster boxes.
//
//  Return:   n/a
//

   MSG("BubAlg", Msg::kVerbose)
      << "Starting AlgMSTClusterList::RunCloseZAlg()" << endl;

// Check for single cluster.
   Int_t nboxnew = toay->GetEntries();
   if (nboxnew <= 1) return;

   Int_t nboxold;

// Loop until no clusters are joined.
   do {

// Save start number of boxes for later check.
      nboxold = nboxnew;

// Create new array to store ClusterBoxes.
      TObjArray *cbay = new TObjArray();

// Store first ClusterBox in new array.
      TIter tyItr(toay);
      cbay->Add(toay->Remove(tyItr()));

// Iterate over remaining ClusterBoxes
      while (ClusterBox *ty = dynamic_cast<ClusterBox *>(tyItr())) {
         toay->Remove(ty);

// Iterate over new ClusterBoxes to search for overlap.
         TIter cbItr(cbay);
         while (ClusterBox *cb = dynamic_cast<ClusterBox *>(cbItr())) {
            if (IsCloseZ(*cb, *ty, sep)) {
               cb->Join(*ty);
               delete ty;
               ty = 0;
               break;
            }
         }

// No overlap found, so add as new ClusterBox.
         cbay->Add(ty);
      }

// Clean-up for next loop.
      toay->Delete();
      delete toay;
      toay = cbay;
      nboxnew = toay->GetEntries();
   } while (nboxnew != nboxold);
}

//......................................................................

void AlgMSTClusterList::RunHoleCrossAlg(Float_t hole, Float_t dist,
                                   TObjArray *&toayU, TObjArray *&toayV)
{
//
//  Purpose:  Join clusters across the coil hole.
//
//  Arguments:
//    hole    in      Maximum hole crossing radius in std length units.
//    dist    in      Maximum separation between ends of cluster boxes.
//    toayU   in/out  Array with digit pair clusters of U-view strips.
//    toayV   in/out  Array with digit pair clusters of V-view strips.
//
//  Return:   n/a
//
//  Notes:    Current implementation allows only clusters in particular
//            direction to be joined.
//

   MSG("BubAlg", Msg::kVerbose)
      << "Starting AlgMSTClusterList::RunHoleCrossAlg()" << endl;

   TObjArray outay;
   TObjArray inay;

// Iterate over clusters in both sets.
   TIter cbUItr(toayU);
   TIter cbVItr(toayV);
   while (ClusterBox *cbU = dynamic_cast<ClusterBox *>(cbUItr())) {
      while (ClusterBox *cbV = dynamic_cast<ClusterBox *>(cbVItr())) {

// Check whether outward crossing.
         if (IsOutCross(*cbU, *cbV, hole, dist)) {
            TObjArray *outpair = new TObjArray(2);
            outpair->AddAt(cbU, 0);
            outpair->AddAt(cbV, 1);
            outay.Add(outpair);
         }

// Check whether inward crossing.
         if (IsInCross(*cbU, *cbV, hole, dist)) {
            TObjArray *inpair = new TObjArray(2);
            inpair->AddAt(cbU, 0);
            inpair->AddAt(cbV, 1);
            inay.Add(inpair);
         }
      }
      cbVItr.Reset();
   }

// Iterate over all outward crossing clusters.
   TIter outItr(&outay);
   while (TObjArray *outpair = dynamic_cast<TObjArray*>(outItr())) {
      ClusterBox *outU = dynamic_cast<ClusterBox*>(outpair->At(0));
      if (!outU) continue;
      ClusterBox *outV = dynamic_cast<ClusterBox*>(outpair->At(1));
      if (!outV) continue;

// Iterate over all inward crossing clusters.
      TIter inItr(&inay);
      while (TObjArray *inpair = dynamic_cast<TObjArray*>(inItr())) {
         ClusterBox *inU = dynamic_cast<ClusterBox*>(inpair->At(0));
         if (!inU) continue;
         ClusterBox *inV = dynamic_cast<ClusterBox*>(inpair->At(1));
         if (!inV) continue;

// Join matching clusters.
         if ((inU->GetYhi() < outU->GetYlo())
             && (inV->GetYhi() < outV->GetYlo())) {
            inU->Join(*outU);
            inV->Join(*outV);
            delete toayU->Remove(outU);
            outU = 0;
            delete toayV->Remove(outV);
            outV = 0;
            delete outay.Remove(outpair);
            outpair = 0;
            delete inay.Remove(inpair);
            inpair = 0;
            break;
         }
      }
   }
}

//......................................................................

void AlgMSTClusterList::RunMSTAlg(Float_t dist, TObjArray *&toay)
{
//
//  Purpose:  Group CandDigiPairs into cluster using a minimum spanning
//            tree algorithm.
//
//  Arguments:
//    dist    in      Maximum separation between digit pairs within a
//                    cluster
//    toay    in/out  Array with initial cluster boxes. It is replaced
//                    by an array with final cluster boxes.
//
//  Return:   n/a
//

   MSG("BubAlg", Msg::kVerbose)
      << "Starting AlgMSTClusterList::RunMSTAlg()" << endl;

// Check for empty array or single cluster.
   Int_t nboxnew = toay->GetEntries();
   if (nboxnew <= 1) return;

   Int_t nboxold;

// Loop until no clusters are joined.
   do {

// Save start number of boxes for later check.
      nboxold = nboxnew;

// Create new array to store ClusterBoxes.
      TObjArray *cbay = new TObjArray();

// Store first ClusterBox in new array.
      TIter tyItr(toay);
      cbay->Add(toay->Remove(tyItr()));

// Iterate over remaining ClusterBoxes
      while (ClusterBox *ty = dynamic_cast<ClusterBox *>(tyItr())) {
         toay->Remove(ty);

// Iterate over new ClusterBoxes to search for overlap.
         TIter cbItr(cbay);
         while (ClusterBox *cb = dynamic_cast<ClusterBox *>(cbItr())) {
            if (IsOverlap(*cb, *ty, dist)) {
               cb->Join(*ty);
               delete ty;
               ty = 0;
               break;
            }
         }

// No overlap found, so add as new ClusterBox.
         cbay->Add(ty);
      }

// Clean-up for next loop.
      toay->Delete();
      delete toay;
      toay = cbay;
      nboxnew = toay->GetEntries();
   } while (nboxnew != nboxold);
}

//......................................................................

void AlgMSTClusterList::Trace(const char *c) const
{
//
//  Purpose:  Trace the AlgMSTClusterList.
//
//  Arguments:
//    c          in    String tag for the trace.
//
//  Return:   n/a
//

  MSG("BubCand", Msg::kDebug)
     << "**********Begin AlgMSTClusterList::Trace(\"" << c << "\")"
     << endl
     << "**********End AlgMSTClusterList::Trace(\"" << c << "\")"
     << endl;
}