////////////////////////////////////////////////////////////////////////
// $Id: AlgStraightCluster.cxx,v 1.13 2007/01/15 19:52:00 rhatcher Exp $
//
// AlgStraightCluster
//
// Begin_Html<img src="../../pedestrians.gif" align=center>
// <a href="../source_warning.html">Warning for beginners</a>.<br> 
//
// This is an Algorithm class to fill a CandStraightCluster with
// CandDigiPairHandles. The CandDigiPairHandles point to CandDigiPairs
// that lie close to a straight line.
//
// Author:  P.S. Miyagawa 10/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 "Candidate/CandContext.h"
#include "MessageService/MsgService.h"
#include "BubbleSpeak/AlgStraightCluster.h"
#include "BubbleSpeak/CandStraightClusterHandle.h"
#include "BubbleSpeak/CandDigiPairHandle.h"
#include "BubbleSpeak/CandDigiPairListHandle.h"

#include "TMath.h"

ClassImp(AlgStraightCluster)

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

CVSID("$Id: AlgStraightCluster.cxx,v 1.13 2007/01/15 19:52:00 rhatcher Exp $");

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

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

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

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

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

void AlgStraightCluster::RunAlg(AlgConfig & /* ac */, CandHandle &ch,
                                               CandContext &cx)
{
//
//  Purpose:  Fills the CandStraightCluster with CandDigiPairs
//            constructed using the supplied CandDigiPairList.
//
//  Argument:
//    ac        in    AlgConfig (not used).
//    ch        in    Handle to the CandStraightCluster to fill.
//    cx        in    CandContext containing a CandDigiPairList.
//
//  Return:   n/a
//

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

// Check for class of handle.
   assert(ch.InheritsFrom("CandStraightClusterHandle"));
   CandStraightClusterHandle &clh =
                           dynamic_cast<CandStraightClusterHandle&>(ch);

   assert(cx.GetDataIn());

// Check for class of input and create iterator.
   TIter *chhItr = 0;
   if (cx.GetDataIn()->InheritsFrom("TObjArray")) {
      const TObjArray *toay =
                       dynamic_cast<const TObjArray*>(cx.GetDataIn());
      chhItr = new TIter(toay);
   }
   else if (cx.GetDataIn()->InheritsFrom("CandDigiPairListHandle")) {
      const CandDigiPairListHandle *chlh =
            dynamic_cast<const CandDigiPairListHandle*>(cx.GetDataIn());
      chhItr = new TIter(chlh->GetDaughterIterator());
   }
   else if (cx.GetDataIn()->InheritsFrom("CandMSTClusterHandle")) {
      const CandMSTClusterHandle *cch =
            dynamic_cast<const CandMSTClusterHandle*>(cx.GetDataIn());
      chhItr = new TIter(cch->GetDaughterIterator());
      clh.SetCandSlice(cch->GetCandSlice());
   }

// Initialize running totals.
   Float_t wtsum = 0;
   Float_t tsum  = 0;
   Float_t t2sum = 0;
   Float_t zsum  = 0;
   Float_t z2sum = 0;
   Float_t tzsum = 0;

   Float_t tmin  =  1e7;
   Float_t tmax  = -1e7;
   Float_t zmin  =  1e7;
   Float_t zmax  = -1e7;

// Iterate over CandDigiPair daughters.
   while (CandDigiPairHandle *chh =
             dynamic_cast<CandDigiPairHandle *>((*chhItr)())) {
      if (!chh->GetCandRecord()) chh->SetCandRecord(cx.GetCandRecord());
      clh.AddDaughterLink(*chh);

// Update running totals.
      Float_t tpos = chh->GetTPos();
      Float_t zpos = chh->GetZPos();
      Float_t wt   = chh->GetCharge();

      wtsum += wt;
      tsum  += wt * tpos;
      t2sum += wt * tpos * tpos;
      zsum  += wt * zpos;
      z2sum += wt * zpos * zpos;
      tzsum += wt * tpos * zpos;

      if (tpos < tmin) tmin = tpos;
      if (tpos > tmax) tmax = tpos;
      if (zpos < zmin) zmin = zpos;
      if (zpos > zmax) zmax = zpos;
   }

// Calculate cluster statistics and fit values.
   MSG("BubAlg", Msg::kVerbose) << "Calculating fit parameters." << endl;
   Int_t fitmode;
   Float_t det, inter, intererr, slope, slopeerr;
   if (zmax == zmin) {
      MSG("BubAlg", Msg::kVerbose) << "zmax == zmin" << endl;
      fitmode = 1;
      det = 0;
      inter = zmin;
      intererr = 0;
      slope = 0;
      slopeerr = 0;
   }
   else if (tmax == tmin) {
      MSG("BubAlg", Msg::kVerbose) << "tmax == tmin" << endl;
      fitmode = 0;
      det = 0;
      inter = tmin;
      intererr = 0;
      slope = 0;
      slopeerr = 0;
   }
   else if ((zmax - zmin) < (tmax - tmin)) {
      MSG("BubAlg", Msg::kVerbose) << "zdiff < tdiff" << endl;
      fitmode = 1;
      det = wtsum * t2sum - tsum * tsum;
      if (!det) det = 1e-9;
      inter = (-tzsum * tsum + zsum * t2sum) / det;
      intererr = 1.18 * TMath::Sqrt(TMath::Abs(tsum / det));
      slope = (-tsum * zsum + wtsum * tzsum) / det;
      slopeerr = 1.18 * TMath::Sqrt(TMath::Abs(wtsum / det));
   }
   else {
      MSG("BubAlg", Msg::kVerbose) << "tdiff < zdiff" << endl;
      fitmode = 0;
      det = wtsum * z2sum - zsum * zsum;
      if (!det) det = 1e-9;
      inter = (-tzsum * zsum + tsum * z2sum) / det;
      intererr = 1.18 * TMath::Sqrt(TMath::Abs(zsum / det));
      slope = (-zsum * tsum + wtsum * tzsum) / det;
      slopeerr = 1.18 * TMath::Sqrt(TMath::Abs(wtsum / det));
   }

// Set cluster statistics and fit values.
   clh.SetWtSum(wtsum);
   clh.SetTsum(tsum);
   clh.SetT2sum(t2sum);
   clh.SetTZsum(tzsum);
   clh.SetZsum(zsum);
   clh.SetZ2sum(z2sum);

   clh.SetFitMode(fitmode);
   clh.SetFitDet(det);
   clh.SetFitInter(inter);
   clh.SetFitInterErr(intererr);
   clh.SetFitSlope(slope);
   clh.SetFitSlopeErr(slopeerr);

// Clean-up.
   delete chhItr;
   chhItr = 0;
   MSG("BubAlg", Msg::kVerbose)
      << "Completed AlgStraightCluster::RunAlg" << endl;
}

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

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

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