////////////////////////////////////////////////////////////////////////
// $Id: TrackClusterSR.cxx,v 1.42 2006/12/01 18:57:23 rhatcher Exp $
//
// TrackClusterSR
//
// Author:  R. Lee 2001.02.26
////////////////////////////////////////////////////////////////////////

#include <cassert>
#include <iostream>
#include <string>
#include <cmath>

#include "TClass.h"

#include "CandTrackSR/TrackClusterSR.h"
#include "Validity/VldContext.h"
#include "Conventions/PlaneView.h"
#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "RecoBase/CandStripHandle.h"
#include "UgliGeometry/UgliGeometry.h"
#include "UgliGeometry/UgliGeomHandle.h"


ClassImp(TrackClusterSR)

//______________________________________________________________________
CVSID("$Id: TrackClusterSR.cxx,v 1.42 2006/12/01 18:57:23 rhatcher Exp $");

//______________________________________________________________________
TrackClusterSR::TrackClusterSR()
{
  fNStrip = 0;
  fMinStrip = 0;
  fMaxStrip = 0;
  fTPos = 0.;
  fTime3D = 0.;
  fTPosRMS = 0.;
  fZPos = 0.;
  fLPos = 0.;
  fCharge = 0.;
  fIsWide = false;
  fIsValid = true;
  fInShowerLikePlane = false;
  fParmMisalignmentError = 30*Munits::mm;
}

//______________________________________________________________________
TrackClusterSR::TrackClusterSR(Double_t misalignmentError)
{
  fNStrip = 0;
  fMinStrip = 0;
  fMaxStrip = 0;
  fTPos = 0.;
  fTime3D = 0.;
  fTPosRMS = 0.;
  fZPos = 0.;
  fLPos = 0.;
  fCharge = 0.;
  fIsWide = false;
  fIsValid = true;
  fInShowerLikePlane = false;
  fParmMisalignmentError = misalignmentError;          // 30*Munits::mm;
}

//______________________________________________________________________
TrackClusterSR::TrackClusterSR(CandStripHandle *strip, 
			       Double_t misalignmentError)
{
  fNStrip = 0;
  fMinStrip = 0;
  fMaxStrip = 0;
  fTPos = 0.;
  fTPosRMS = 0.;
  fZPos = 0.;
  fLPos = 0.;
  fCharge = 0.;
  fTime3D = 0.;
  fIsWide = false;
  fIsValid = true;
  fInShowerLikePlane = false;
  fParmMisalignmentError = misalignmentError;          // 30*Munits::mm;

  AddStrip(strip);
}

//______________________________________________________________________
TrackClusterSR::TrackClusterSR(const TrackClusterSR &rhs) : // Copy-ctor
  TObject(rhs)
, fNStrip(rhs.fNStrip)
, fMinStrip(rhs.fMinStrip)
, fMaxStrip(rhs.fMaxStrip)
, fZPos(rhs.fZPos)
, fCharge(rhs.fCharge)
, fTime3D(rhs.fTime3D)
, fTPos(rhs.fTPos)
, fTPosRMS(rhs.fTPosRMS)
, fLPos(rhs.fLPos)
, fIsWide(rhs.fIsWide)
, fIsValid(rhs.fIsValid)
, fInShowerLikePlane(rhs.fInShowerLikePlane)
, fParmMisalignmentError(rhs.fParmMisalignmentError)
{

// Copy TObjArray's contents
  TIter rhsIter(&rhs.fStripList);
  TObject *rhsObj;
  while ((rhsObj = rhsIter())) {
    CandStripHandle *cshObj;
    if ((cshObj = dynamic_cast<CandStripHandle*>(rhsObj))) {
      fStripList.Add(cshObj->DupHandle());
    }
  }
}

//______________________________________________________________________
TrackClusterSR::~TrackClusterSR()
{
  fStripList.Delete(); // Owned CandHandle*'s from TrackClusterSR vers 9
}

//______________________________________________________________________
void TrackClusterSR::AddStrip(CandStripHandle *strip)
{
  if (fStripList.GetLast()<0) {
    fZPos = strip->GetZPos();
  }

  PlexStripEndId strip_seid = strip->GetStripEndId();
  
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    CandStripHandle *strip1 = 
      dynamic_cast<CandStripHandle*>(fStripList.At(i));
    PlexStripEndId strip1_seid = strip1->GetStripEndId();
    if ( strip_seid == strip1_seid  &&  (strip1->GetTime()==strip->GetTime())  ) return;
  }

  Int_t stripnum = strip_seid.GetStrip();

  fNStrip++;
  if (fNStrip==1 || stripnum<fMinStrip) {
    fMinStrip = stripnum;
  }
  if (fNStrip==1 || stripnum>fMaxStrip) {
    fMaxStrip = stripnum;
  }
  fStripList.Add(strip->DupHandle());
  if (fNStrip==1) {
    fCharge  = strip->GetCharge();
    fTPos    = strip->GetTPos();
    fTPosRMS = 0.;
  } else {
    fCharge  = 0.;
    fTPos    = 0.;
    fTPosRMS = 0.;
    for (Int_t i=0; i<=fStripList.GetLast(); i++) {
      CandStripHandle *strip1 = 
        dynamic_cast<CandStripHandle*>(fStripList.At(i));
      Double_t charge = strip1->GetCharge();
      if (charge <= 0.0) {
        MSG("TrackSR", Msg::kWarning)
          << "AddStrip() strip1 charge "
          << charge << ", avoid corrupting avg/RMS of TPos." << endl;
        continue;  // skip polluting <TPos>,RMS(TPos) w/ funky charge
      }
      Double_t tpos   = strip1->GetTPos();
      //    fTPos    += charge*tpos;
      //  fTPosRMS += charge*tpos*tpos;
      //  fCharge  += charge;
      fTPos    += tpos;
      fTPosRMS += tpos*tpos;
      fCharge  += charge;
    }
    if (fCharge<=0.0) {
      MSG("TrackSR", Msg::kWarning)
        << "AddStrip() fCharge "
        << fCharge << ", avoid division." << endl;
      fCharge  = 1.0e-6;
      fTPos    = 0.0;
      fTPosRMS = 1.0e+6;
    }
    
    //    fTPos /= fCharge;
    //  fTPosRMS /= fCharge;
    // fTPosRMS -= fTPos*fTPos;
    fTPos /= fNStrip;
    fTPosRMS /= fNStrip;
    fTPosRMS -= fTPos*fTPos;
    if (fTPosRMS>=0.) {
      fTPosRMS = sqrt(fTPosRMS);
    } else {
      fTPosRMS = 0.;
    }
  } // fNStrip > 1

}

/// Assign "longitudinal" position of the cluster along the strip.
/// Triggers recalculation of the transverse position. 
void TrackClusterSR::SetLPos(Double_t lpos)
{
    MSG("TrackSR", Msg::kDebug)
          << "void TrackClusterSR::SetLPos(Double_t lpos)" << endl;
    MSG("TrackSR", Msg::kDebug)
          << "Setting fLPos = " << lpos << endl;
    fLPos = lpos;   
    RecalculateTPos();
} // void TrackClusterSR::SetLPos(Double_t)


/// Get rotation corrected tpos of the strip at position along the strip
/// corresponding to fLPos of this track cluster. Cache last tpos value
/// to avoid repeating calculations 
Double_t TrackClusterSR::GetRotationCorrectedTPos(const CandStripHandle* csh) const
{
    MSG("TrackSR", Msg::kDebug)
          << "Double_t TrackClusterSR::GetRotationCorrectedTPos(const CandStripHandle* csh) const" << endl;
          
    assert( csh && "Do not pass NULL pointers!!!");
    
    static CandStripHandle* lastCsh = 0;
    static Double_t lastLpos = 0.;
    static Double_t lastTpos = 0.;
        
    if ( csh != lastCsh || fLPos != lastLpos ) {
        UgliGeomHandle ugh(*csh->GetVldContext());    
        PlexStripEndId strip_seid = csh->GetStripEndId();
        UgliStripHandle ush = ugh.GetStripHandle(strip_seid);
        lastTpos = ush.GetTPos(fLPos);
    }        

    MSG("TrackSR", Msg::kDebug) << "TPos = " << csh->GetTPos() 
        << "; RotationCorrectedTPos = " << lastTpos << "; LPos = " << fLPos << endl;
    return lastTpos;
} // Double_t TrackClusterSR::GetRotationCorrectedTPos(const CandStripHandle*)


/// Recalculate transverse position - needs to be done when setting the
/// "longitudinal" postion along the strip fLPos.
void TrackClusterSR::RecalculateTPos()
{    
  // return when the cluster has no strips
  if ( fNStrip == 0 ) return;
  
  // running sums
  Double_t chargesum = 0.;
  Double_t tpossum = 0.;
  Double_t tposrmssum = 0.;
  
  // loop over strips
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    CandStripHandle *strip = 
      dynamic_cast<CandStripHandle*>(fStripList.At(i));
    
    // cache strip charge, tpos
    Double_t stripcharge = strip->GetCharge();
    Double_t striptpos = GetRotationCorrectedTPos(strip);
    
    // ignore strips with negative charge
    if (stripcharge <= 0.0) {
      MSG("TrackSR", Msg::kWarning)
        << "Strip charge is "
        << stripcharge << ", avoid corrupting avg/RMS of TPos." << endl;
      continue;  // skip polluting <TPos>,RMS(TPos) w/ funky charge
    }
    
    // increment running sums
    chargesum  += stripcharge;
    tpossum    += striptpos;
    tposrmssum += striptpos*striptpos;
  }    

  // set zero values and return if total charge is <=0.    
  if (chargesum <= 0.0) {
    MSG("TrackSR", Msg::kWarning)
        << "AddStrip() fCharge "
        << fCharge << ", avoid division." << endl;
    fCharge  = 1.0e-6;
    fTPos    = 0.0;
    fTPosRMS = 1.0e+6;
    return;
  }
  
  // update TrackCluster fTPos, fTPosRMS  
  fTPos = tpossum/fNStrip;
  fTPosRMS = tposrmssum/fNStrip - fTPos*fTPos;
  if ( fTPosRMS>=0.) {
      fTPosRMS = sqrt(fTPosRMS);
  } else {
      fTPosRMS = 0.;
  }
} // void TrackClusterSR::RecalculateTPos()

//______________________________________________________________________
Double_t TrackClusterSR::DTPos(Double_t tpos) const
{
  if (fNStrip==0) {
    return -99999.;
  }
  CandStripHandle *strip =
                     dynamic_cast<CandStripHandle*>(fStripList.First());
  if (GetNStrip()==1) {
    return fabs(strip->GetTPos()-tpos); 
  }
  Double_t dtpos = 9999999.;
  Bool_t above(0);
  Bool_t below(0);

// assume that trackcluster is contiguous set of strips
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    strip = dynamic_cast<CandStripHandle*>(fStripList.At(i));
    if (strip->GetTPos()>tpos) {
      above = 1;
    }
    if (strip->GetTPos()<tpos) {
      below = 1;
    }
    dtpos = min(dtpos,fabs(strip->GetTPos()-tpos));
  }
  if (above && below) {
    dtpos = 0.;
  }
  return dtpos;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetBegTime() const
{
  if (fNStrip==0) {
    return 0.;
  }
  Double_t begtime=0.;
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    CandStripHandle *strip =
                       dynamic_cast<CandStripHandle*>(fStripList.At(i));
    if (!i || strip->GetBegTime()<begtime) {
      begtime = strip->GetBegTime();
    }
  }
  return begtime;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetCharge() const
{
  if (fNStrip==0) {
    return 0.;
  }
  return fCharge;
}

//______________________________________________________________________
Int_t TrackClusterSR::GetMaxStrip() const
{
  return fMaxStrip;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetMaxTPos() const
{
  if (fNStrip==0) {
    return 0.;
  }
  Double_t tpos=-5.;
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    CandStripHandle *strip =
                       dynamic_cast<CandStripHandle*>(fStripList.At(i));
    if (!i || strip->GetTPos()>tpos) {
      tpos = strip->GetTPos();
    }
  }
  return tpos+0.02054;
}

//______________________________________________________________________
Int_t TrackClusterSR::GetMinStrip() const
{
  return fMinStrip;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetMinTPos() const
{
  if (fNStrip==0) {
    return 0.;
  }
  Double_t tpos=5.;
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    CandStripHandle *strip =
                       dynamic_cast<CandStripHandle*>(fStripList.At(i));
    if (!i || strip->GetTPos()<tpos) {
      tpos = strip->GetTPos();
    }
  }
  return tpos-0.02054;
}

//______________________________________________________________________
Int_t TrackClusterSR::GetNStrip() const
{
  return fNStrip;
}

//______________________________________________________________________
Int_t TrackClusterSR::GetPlane() const
{
  if (fStripList.GetLast()>=0) {
    CandStripHandle *strip =
                     dynamic_cast<CandStripHandle*>(fStripList.First());
    return strip->GetPlane();
  }
  else {
    return 0;
  }
}

//______________________________________________________________________
PlaneView::PlaneView_t TrackClusterSR::GetPlaneView() const
{
  if (fNStrip==0) {
    return PlaneView::kUnknown;
  }
  CandStripHandle *firststrip =
                     dynamic_cast<CandStripHandle*>(fStripList.First());
  return firststrip->GetPlaneView();
}

//______________________________________________________________________
const TObjArray* TrackClusterSR::GetStripList() const
{
  return &fStripList;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetTime3D() const
{
  return fTime3D;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetTPos() const
{
  if (fNStrip==0) {
    return 0.;
  }
  return fTPos;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetTPosError() const
{
  Double_t err = 0.;
  if (fNStrip==0) {
    err = 0.;
  } else if (fNStrip==1) {

// 4.108 cm / sqrt(12)
    err = 0.01185877;
  } else {
   
    err = max((GetMaxTPos()-GetMinTPos())*0.288,0.01185877);
	       //    err = max(fTPosRMS,0.01185877);
  }

  // add value in quadrature to reflect misalignments
  err = sqrt(err*err+fParmMisalignmentError*fParmMisalignmentError);
  return err;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetTPosRMS() const
{
  if (fNStrip==0) {
    return 0.;
  }
  return fTPosRMS;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetZPos() const
{
  return fZPos;
}

//______________________________________________________________________
Double_t TrackClusterSR::GetMisalignmentError() const
{
  return fParmMisalignmentError;
}

//______________________________________________________________________
Bool_t TrackClusterSR::InShowerLikePlane() const
{

// returns true if cluster is in a showerlike plane - ie one with a wide
// cluster
  return fInShowerLikePlane;
}

//______________________________________________________________________
Bool_t TrackClusterSR::IsContained(CandStripHandle *target) const
{ 
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    CandStripHandle *strip =
                       dynamic_cast<CandStripHandle*>(fStripList.At(i));

// check if object being pointed to is the same
    if (*strip == *target) {
      return 1;
    }
  }
  return 0;
}

//______________________________________________________________________
Bool_t TrackClusterSR::IsDoubleEnded() const
{
  for (Int_t i=0; i<=fStripList.GetLast(); i++) {
    CandStripHandle *strip =
                       dynamic_cast<CandStripHandle*>(fStripList.At(i));
    if (!i &&
           strip->GetVldContext()->GetDetector()==Detector::kNear) {
      return kFALSE;
    }
    if (strip->GetNDaughters()==2) {
      return kTRUE;
    }
  }
  return kFALSE;
}

//______________________________________________________________________
Bool_t TrackClusterSR::IsEquivalent(const TObject *rhs) const
{
  const TrackClusterSR *rcl = dynamic_cast<const TrackClusterSR*>(rhs);

// compare simple members
  if ( !(this->fNStrip            == rcl->fNStrip            &&
         this->fInShowerLikePlane == rcl->fInShowerLikePlane )) {
    return false;
  }

// compare TObjArray
  if (this->fStripList.GetEntries() != rcl->fStripList.GetEntries()) {
    return false;
  }
  else {
    TIter thisIter(&this->fStripList);
    TIter rhsIter(&rcl->fStripList);
    TObject *thisObj, *rhsObj;
    while ( ( (thisObj = thisIter.Next() ) != NULL) &&
            ( (rhsObj  = rhsIter.Next()  ) != NULL)    ) {
      CandStripHandle *thisStrip =
                                dynamic_cast<CandStripHandle*>(thisObj);
      CandStripHandle *rhsStrip =
                                 dynamic_cast<CandStripHandle*>(rhsObj);
      if (!thisStrip || !rhsStrip) return false; // not CandStripHandle*
      if (*thisStrip != *rhsStrip) return false;  // Deep Handle compare
    }
  }
  return true;
}

//______________________________________________________________________
Bool_t TrackClusterSR::IsValid() const
{
// returns true if cluster is a valid one to consider
  return fIsValid;
}

//______________________________________________________________________
Bool_t TrackClusterSR::IsWide() const
{

// returns true if cluster is a wide cluster
  return fIsWide;
}

//______________________________________________________________________
void TrackClusterSR::SetInShowerLikePlane(Bool_t showerLike)
{

// set the flag for this being in a shower like plane - ie 
// the plane contains at least one wide cluster
  fInShowerLikePlane = showerLike;
}

//______________________________________________________________________
void TrackClusterSR::SetIsValid(Bool_t isvalid)
{

// set the flag for this being a wide cluster - ie one with
// more strips in it than expected based on slope
  fIsValid = isvalid;
}

//______________________________________________________________________
void TrackClusterSR::SetIsWide(Bool_t iswide)
{

// set the flag for this being a wide cluster - ie one with
// more strips in it than expected based on slope
  fIsWide = iswide;
}

//______________________________________________________________________
void TrackClusterSR::SetTime3D(Double_t dvar)
{
  fTime3D = dvar;
}

//______________________________________________________________________
NavKey TrackClusterSR::KeyFromPlane(const TrackClusterSR *tcluster)
{
  return tcluster->GetPlane();
}

//______________________________________________________________________
NavKey TrackClusterSR::KeyFromPlaneTPosTime(
                                         const TrackClusterSR *tcluster)
{
   Int_t iplane = tcluster->GetPlane();
   Double_t tpos = tcluster->GetMinTPos();
   Int_t itpos = static_cast<Int_t>(tpos);
   Double_t time = tcluster->GetBegTime();
   Int_t itime = static_cast<Int_t>(time*.1/18.7);

   Int_t navkey = (iplane<<22) | (itpos<<11) | itime;

   return navkey;
}

//______________________________________________________________________
void TrackClusterSR::Streamer(TBuffer &R__b)
{

// Stream an object of class TrackClusterSR.
  if (R__b.IsReading()) {
    UInt_t R__s, R__c;
    Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
    TrackClusterSR::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
    if (R__v < 9) {

// fStripList contents owned from TrackClusterSR version 9
      CandHandle *ch;
      TIter cliter(&fStripList);
      while ((ch = dynamic_cast<CandHandle *>(cliter())))
        fStripList.AddAt(ch->DupHandle(), fStripList.IndexOf(ch));
    }
  }
  else {
    TrackClusterSR::Class()->WriteBuffer(R__b, this);
  }
}

XXXITRIMP(TrackClusterSR)