////////////////////////////////////////////////////////////////////////
// $Id: AlgBandClusterList.cxx,v 1.5 2003/08/22 18:35:43 miyagawa Exp $
//
// AlgBandClusterList
//
// 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
// CandDigiPairs contained in a CandSliceList.  It first defines a
// central track from strips that have digits on both ends.  It then
// places any strips that are in a band about this track into the main
// cluster.  Any remaining strips are placed in separate clusters.
//
// Author:  P.S. Miyagawa 2/2002
//
// 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 "UgliGeometry/UgliGeomHandle.h"
#include "BubbleSpeak/AlgBandClusterList.h"
#include "BubbleSpeak/CandMSTCluster.h"
#include "BubbleSpeak/CandMSTClusterHandle.h"
#include "BubbleSpeak/CandDigiPairHandle.h"

ClassImp(AlgBandClusterList)

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

CVSID("$Id: AlgBandClusterList.cxx,v 1.5 2003/08/22 18:35:43 miyagawa Exp $");

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

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

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

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

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

void AlgBandClusterList::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 CandDigiPairList to
//                    cluster.
//
//  Return:   n/a
//

   MSG("BubAlg", Msg::kVerbose)
      << "Starting AlgBandClusterList::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 AlgBandClusterList::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 AlgBandClusterList::RunAlgOnSlice()" << endl;

// Save config parameters.
   Double_t bandspan = ac.GetDouble("BandSpan");

// Define central track.
   Double_t utsum  = 0;
   Double_t ut2sum = 0;
   Double_t utzsum = 0;
   Double_t uzsum  = 0;
   Double_t uz2sum = 0;
   Double_t uwtsum = 0;
   Double_t vtsum  = 0;
   Double_t vt2sum = 0;
   Double_t vtzsum = 0;
   Double_t vzsum  = 0;
   Double_t vz2sum = 0;
   Double_t vwtsum = 0;
   TIter chhItr(csh->GetDaughterIterator());
   while (const CandDigiPairHandle *chh =
             dynamic_cast<const CandDigiPairHandle *>(chhItr())) {

      // Select strips with double-sided readout.
      Bool_t east = kFALSE;
      Bool_t west = kFALSE;
      TIter cdhItr(chh->GetDaughterIterator());
      while (const CandDigitHandle *cdh =
                dynamic_cast<const CandDigitHandle *>(cdhItr())) {
         StripEnd::StripEnd_t ende = cdh->GetPlexSEIdAltL().GetEnd();
         if (ende == StripEnd::kEast) east = kTRUE;
         else if (ende == StripEnd::kWest) west = kTRUE;

         // Update stats for double-sided.
         if (east && west) {
            Double_t tpos = chh->GetTPos();
            Double_t zpos = chh->GetZPos();
            Double_t wt = chh->GetCharge();
            switch (chh->GetPlaneView()) {
            case PlaneView::kU:
               utsum  += wt * tpos;
               ut2sum += wt * tpos * tpos;
               utzsum += wt * tpos * zpos;
               uzsum  += wt * zpos;
               uz2sum += wt * zpos * zpos;
               uwtsum += wt;
               break;
            case PlaneView::kV:
               vtsum  += wt * tpos;
               vt2sum += wt * tpos * tpos;
               vtzsum += wt * tpos * zpos;
               vzsum  += wt * zpos;
               vz2sum += wt * zpos * zpos;
               vwtsum += wt;
               break;
            default:
               break;
            }
            continue;
         }
      }
   }

   // Calculate central track parameters.
   Double_t udet = uwtsum * uz2sum - uzsum * uzsum;
   Double_t uinter, uslope;
   if (udet != 0.) {
      uinter = (utsum * uz2sum - utzsum * uzsum) / udet;
      uslope = (uwtsum * utzsum - uzsum * utsum) / udet;
   }
   else {
      uinter = (uwtsum!=0. ? utsum / uwtsum : -1e12);
      uslope = 1e12;
   }
   Double_t uslope2_1 = uslope * uslope + 1;

   Double_t vdet = vwtsum * vz2sum - vzsum * vzsum;
   Double_t vinter, vslope;
   if (vdet != 0.) {
      vinter = (vtsum * vz2sum - vtzsum * vzsum) / vdet;
      vslope = (vwtsum * vtzsum - vzsum * vtsum) / vdet;
   }
   else {
      vinter = (vwtsum!=0. ? vtsum / vwtsum : -1e12);
      vslope = 1e12;
   }
   Double_t vslope2_1 = vslope * vslope + 1;

// Find strip width using first digitized strip.
   Float_t cellwidth = 0;
   chhItr.Reset();
   while (const CandDigiPairHandle *chht =
             dynamic_cast<const CandDigiPairHandle *>(chhItr())) {
      UgliGeomHandle ugh(*chht->GetVldContext());
      UgliStripHandle ush = ugh.GetStripHandle(chht->GetStripEndId());
      if (ush.IsValid()) {
         cellwidth = ush.GetHalfWidth();
         cellwidth *= 2.0;
         break;
      }
   }

// Set difference check limit.
   Double_t d2check = bandspan * cellwidth;
   d2check *= d2check;

// Separate strips into clusters.
   TObjArray *mainU = 0;
   TObjArray *mainV = 0;
   TObjArray *others = 0;
   TObjArray *arrayA = 0;
   TObjArray *arrayB = 0;
   TObjArray *shield = 0;
   chhItr.Reset();
   while (CandDigiPairHandle *chh =
             dynamic_cast<CandDigiPairHandle *>(chhItr())) {
      Double_t xdat = chh->GetZPos();
      Double_t ydat = chh->GetTPos();
      Double_t d2diff = 0;

      // 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:
         d2diff = ydat - uinter - uslope * xdat;
         d2diff *= d2diff;
         d2diff /= uslope2_1;
         if (d2diff <= d2check) {
            if (!mainU) mainU = new TObjArray();
            mainU->Add(chh);
         }
         else {
            if (!others) others = new TObjArray();
            others->Add(chh);
         }
         break;

      case PlaneView::kV:
         d2diff = ydat - vinter - vslope * xdat;
         d2diff *= d2diff;
         d2diff /= vslope2_1;
         if (d2diff <= d2check) {
            if (!mainV) mainV = new TObjArray();
            mainV->Add(chh);
         }
         else {
            if (!others) others = new TObjArray();
            others->Add(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;
      }
   }

// 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 arrays to create CandMSTClusters.
   // Main U-view cluster.
   if (mainU) {
      cxx.SetDataIn(mainU);
      CandMSTClusterHandle clh = CandMSTCluster::MakeCandidate(ah, cxx);
      clh.SetCandSlice(csh);
      ch.AddDaughterLink(clh);
      delete mainU;
      mainU = 0;
   }

   // Main V-view cluster.
   if (mainV) {
      cxx.SetDataIn(mainV);
      CandMSTClusterHandle clh = CandMSTCluster::MakeCandidate(ah, cxx);
      clh.SetCandSlice(csh);
      ch.AddDaughterLink(clh);
      delete mainV;
      mainV = 0;
   }

   // Other clusters.
   if (others) {
      TIter chhItr(others);
      while (CandDigiPairHandle *chh =
                dynamic_cast<CandDigiPairHandle *>(chhItr())) {
         TObjArray *toay = new TObjArray();
         toay->Add(chh);
         cxx.SetDataIn(toay);
         CandMSTClusterHandle clh
                               = CandMSTCluster::MakeCandidate(ah, cxx);
         clh.SetCandSlice(csh);
         ch.AddDaughterLink(clh);
         delete toay;
         toay = 0;
      }
      delete others;
      others = 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 AlgBandClusterList::Trace(const char *c) const
{
//
//  Purpose:  Trace the AlgBandClusterList.
//
//  Arguments:
//    c          in    String tag for the trace.
//
//  Return:   n/a
//

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