AlgDeMuxCosmics.cxx

Go to the documentation of this file.

// $Id: AlgDeMuxCosmics.cxx,v 1.101 2007/02/04 05:55:06 rhatcher Exp $
//
// AlgDeMuxCosmics
//
// An Algorithm class that DeMultiplexes a cosmic ray event and adjusts the
//weights in the PlexSEIdAltLists of the CandDigits from that event
//
// Author:  B. Rebel 11/2000

#include <cassert>

#include "Candidate/CandContext.h"
#include "Candidate/CandHandle.h"
#include "CandData/CandHeader.h"
#include "CandData/CandRecord.h"
#include "CandDigit/CandDigitHandle.h"
#include "CandDigit/CandDigitListHandle.h"
#include "CandDigit/CandDeMuxDigitHandle.h"
#include "CandDigit/CandDeMuxDigitListHandle.h"
#include "MessageService/MsgService.h"
#include "Algorithm/AlgConfig.h"
#include "DeMux/DmxStatus.h"
#include "DeMux/DmxHypothesis.h"
#include "DeMux/DmxPlane.h"
#include "DeMux/DmxShowerPlane.h"
#include "DeMux/DmxMuonPlane.h"
#include "DeMux/DmxUtilities.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "DeMux/AlgDeMuxCosmics.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "TObjArray.h"

#include "TFolder.h"
#include "TObjArray.h"
#include "TROOT.h"
#include "TMath.h"

//define the NavKey to sort the planes
static NavKey KeyOnPlane( const DmxPlane *plane){
   return plane->GetPlaneNumber();
}

//define the NavKeys to select U and V planes
NavKey KeyOnUView( const DmxPlane *plane){
   return plane->GetPlaneView() == PlaneView::kU;
}

NavKey KeyOnVView( const DmxPlane *plane){
   return plane->GetPlaneView() == PlaneView::kV;
}

//define the NavKeys to select only valid U and V planes
NavKey KeyOnValidU( const DmxPlane *plane){
   return plane->GetPlaneView() == PlaneView::kU && (Int_t)plane->IsValid() == 1;
}

NavKey KeyOnValidV( const DmxPlane *plane){
   return plane->GetPlaneView() == PlaneView::kV && (Int_t)plane->IsValid() == 1;
}

ClassImp(AlgDeMuxCosmics)
  

//-----------------------------------------------------------------------

CVSID("$Id: AlgDeMuxCosmics.cxx,v 1.101 2007/02/04 05:55:06 rhatcher Exp $");

//-----------------------------------------------------------------------

AlgDeMuxCosmics::AlgDeMuxCosmics() :
  fRequiredMatedSignalRatio(0.33),
  fSlopeRMS(0.),
  fStrayPlanes(0),
  fPlanesInSet(6),
  fStrayCut(6)
{
  //default constructor
}

//-----------------------------------------------------------------------

AlgDeMuxCosmics::~AlgDeMuxCosmics()
{

}

//-----------------------------------------------------------------------
//run the demuxing on an event

void AlgDeMuxCosmics::RunAlg(AlgConfig &acd, CandHandle &ch, CandContext & /* cx */)
{

  assert( ch.InheritsFrom("CandDigitListHandle") );
  CandDeMuxDigitListHandle &cdlh = dynamic_cast<CandDeMuxDigitListHandle&>(ch);
  
  MSG("AlgDeMuxCosmics", Msg::kDebug) << cdlh.GetNDaughters() << " digits in AlgDeMuxCosmics" << endl;
  
  VldContext vldc = ch.GetCandRecord()->GetCandHeader()->GetVldContext();
  if (vldc.GetDetector() != Detector::kFar) {
    static int msglimit = 20; // no more than 20 warnings about the detector
    if (msglimit) {
      MSG("DMX",Msg::kDebug) 
        << "AlgDeMuxCosmics::RunAlg() can not DeMux "
        << Detector::AsString(vldc.GetDetector())
        << " detector.  Skip futher DeMux processing."
        << endl;
      if (--msglimit == 0) 
        MSG("DMX",Msg::kDebug) 
          << " ... last message of this type." << endl;
    }
    return;
  }

  //get the AlgConfigDeMux object
  fPlanesInSet = acd.GetInt("PlanesInSet");
  fRequiredMatedSignalRatio = acd.GetDouble("RatioMatedSignalForValid");
  fStrayCut = acd.GetInt("StrayDeltaStripLimit");
  //MSG("DMX", Msg::kDebug) << "starting cosmics algorithm" << endl;

  //get the DmxStatus object
  // Find the DmxStatus object or create one for needed scratch space
  DmxStatus *status = dynamic_cast<DmxStatus *>(gROOT->GetRootFolder()->FindObject("Loon/DeMux/DmxStatus"));
  bool tempStatus = false;
  if (status == 0) {
    MSG("DMXX", Msg::kDebug) << "//root/Loon/DeMux/DmxStatus not found."
                              << " Create a temporary DmxStatus." << endl;
    status = new DmxStatus;      // Make a temporary DmxStatus if needed
    tempStatus = true;
  }
  else status->ResetStatus(); 

  status->SetEventNumber(ch.GetCandRecord()->GetCandHeader()->GetSnarl());
  
  //instantiate the DmxInitialize object and then fill DmxStatus with the event information     
  DmxUtilities util;
  
  util.FillPlaneArray(status, cdlh, acd);
  const TObjArray *planeArray = status->GetPlaneArray();
  
  //create the TObjectItr over planes and program it to sort the planes
  DmxPlaneItr planeItr(planeArray);
  DmxPlaneKeyFunc *pnKF = planeItr.CreateKeyFunc();
  pnKF->SetFun(KeyOnPlane);
  planeItr.GetSet()->AdoptSortKeyFunc(pnKF);
  pnKF = 0;
  
  planeItr.ResetFirst();

  // while(planeItr.IsValid()){
//     cout << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber() << endl;
//     planeItr.Next();
//   }
//   planeItr.ResetFirst();

  status->SetNumberOfPlanes(planeItr.SizeSelect());
  status->SetVertexPlaneNumber(util.FindVertexPlane(planeItr));
  
  //demux the event if there is an identified vertex
  if( status->GetVertexPlaneNumber() != -1){
    
    planeItr.ResetFirst();
    planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), 500);
    planeItr.GetSet()->Update();
    //set the end plane number since there was a vertex plane
    status->SetEndPlaneNumber(util.FindEndPlane(planeItr));
    planeItr.GetSet()->ClearSlice(false);
    planeItr.ResetFirst();

    //set the z position of the vertex plane
    planeItr.GetSet()->Slice(status->GetVertexPlaneNumber());
    planeItr.GetSet()->Update();
    planeItr.ResetFirst();
    status->SetVertexPlaneZPosition(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetZPosition());

    planeItr.GetSet()->ClearSlice(false);
    planeItr.ResetFirst();

    //make the iterators for the window, valid planes, and vertex
    DmxPlaneItr vertexPlaneItr(planeArray);
    DmxPlaneItr validPlaneItr(planeArray);
    DmxPlaneItr windowPlaneItr(planeArray);
       
    //create a KeyFunc to sort planes by number
    DmxPlaneKeyFunc *vertexPlaneKF = vertexPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *validPlaneKF = validPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *validWindowPlaneKF = windowPlaneItr.CreateKeyFunc();
       
    //program the KeyFunc with the sort function
    vertexPlaneKF->SetFun(KeyOnPlane);
    validPlaneKF->SetFun(KeyOnPlane);
    validWindowPlaneKF->SetFun(KeyOnPlane);
       
    //get the NavSet from the iterator and pass the KeyFunc to it
    vertexPlaneItr.GetSet()->AdoptSortKeyFunc(vertexPlaneKF);
    validPlaneItr.GetSet()->AdoptSortKeyFunc(validPlaneKF);
    windowPlaneItr.GetSet()->AdoptSortKeyFunc(validWindowPlaneKF);
       
    //clear the KF pointer because we no longer own the KeyFunc
    vertexPlaneKF = 0;
    validPlaneKF = 0;
    validWindowPlaneKF = 0;
    planeItr.ResetFirst();
    validPlaneItr.ResetFirst();
    windowPlaneItr.ResetFirst();
        
    //now program a key function to select on plane orientation - U first
    DmxPlaneKeyFunc *orientValidUKF = validPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientValidWindowUKF = windowPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientUKF = planeItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientValidVKF = validPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientValidWindowVKF = windowPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientVKF = planeItr.CreateKeyFunc();
    
    Int_t viewCtr = 0;
    while( viewCtr < 2 ){
      fSlopeRMS = 0.;
      fStrayPlanes = 0;
      
      if( viewCtr == 0){
        //program this function with the select function
        orientUKF->SetFun(KeyOnUView);
        orientValidUKF->SetFun(KeyOnValidU);
        orientValidWindowUKF->SetFun(KeyOnValidU);
                
        //adopt it as a selection function
        planeItr.GetSet()->AdoptSelectKeyFunc(orientUKF);
        orientUKF = 0;
        planeItr.Reset();
        validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidUKF);
        orientValidUKF = 0;
        validPlaneItr.Reset();
        windowPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidWindowUKF);
        orientValidWindowUKF = 0;
        windowPlaneItr.Reset();
      }
      else if(viewCtr == 1){
        //program this function with the select function
        orientVKF->SetFun(KeyOnVView);
        orientValidVKF->SetFun(KeyOnValidV);
        orientValidWindowVKF->SetFun(KeyOnValidV);
                
        //adopt it as a selection function
        planeItr.GetSet()->AdoptSelectKeyFunc(orientVKF);
        orientVKF = 0;
        planeItr.Reset();
        validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidVKF);
        orientValidVKF = 0;
        validPlaneItr.Reset();
        windowPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidWindowVKF);
        orientValidWindowVKF = 0;
        windowPlaneItr.Reset();
      }

      //  validPlaneItr.ResetLast();
//      MSG("AlgDeMuxCosmics", Msg::kDebug)<<"event " << status->GetEventNumber() << endl;
//        while(validPlaneItr.IsValid()){
//      MSG("DMX1", Msg::kDebug)<<"\tvalid plane " 
//                             << dynamic_cast<DmxPlane* >(validPlaneItr.Ptr())->GetPlaneNumber() << endl; 
//      validPlaneItr.Prev();
//        }
      
//        validPlaneItr.Reset();
      
      //  DmxPlane *plane = 0;
      
//        while( (plane = dynamic_cast<DmxPlane *>(validPlaneItr())) ){
//      MSG("DMX1", Msg::kDebug) << "\t" << plane->GetPlaneNumber() <<endl;
//        }
//        validPlaneItr.Reset();

      //slice to those planes from the vertex on
      planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
      planeItr.GetSet()->Update();

      validPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
      validPlaneItr.GetSet()->Update();

      windowPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
      windowPlaneItr.GetSet()->Update();
      
      //if only 2 planes set them to their best hyps
      if(validPlaneItr.SizeSelect() == 2){
        validPlaneItr.ResetFirst();

        //variables for straight line fit
        Float_t cogBeg = -1., zBeg = -1., cogEnd = -1., zEnd = -1.;

        while( validPlaneItr.IsValid() ){
          
          //only set it if it is a shower plane, muon planes are set in the constructor
          if(validPlaneItr.Ptr()->GetPlaneType() == DmxPlaneTypes::kShower
             && !validPlaneItr.Ptr()->IsGolden())
            dynamic_cast<DmxShowerPlane *>(validPlaneItr.Ptr())->SetStrips("best");

          //cout << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetCoG() << "\t" 
          //   << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber() 
          //   << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetZPosition() << endl;
          
          //set the variables for the straighe line fit between the two planes
          if(cogBeg == -1.){
            cogBeg = validPlaneItr.Ptr()->GetCoG();
            zBeg = validPlaneItr.Ptr()->GetZPosition();
          }
          else{
            cogEnd = validPlaneItr.Ptr()->GetCoG();
            zEnd = validPlaneItr.Ptr()->GetZPosition();
          }
          
          validPlaneItr.Next();
        }//end loop over planes
        
        //set the nonvalid planes if they exist
        if(planeItr.SizeSelect()>2 && cogEnd-cogBeg != 0. && cogBeg != -1. && cogEnd != -1.){
          Double_t slope = (cogEnd-cogBeg)/(zEnd - zBeg);
          Double_t intercept = cogBeg - slope*zBeg;
          //cout << slope << "\t" << intercept << endl;
          SetPlanesToDeterminedFit(planeItr, intercept, slope);
          fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
        }

        //fStrayPlanes -= util.CheckFit(validPlaneItr);
        validPlaneItr.ResetFirst();

      }//end if only 2 valid planes
      else if(validPlaneItr.SizeSelect() >2 ){

        //if there are only enough planes for 1 instances of the sliding window + initial fit,
        //change fPlanesInSet to give you 3 + the initial fit, but at least 4 always
        if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==6) fPlanesInSet -= 2;
        else if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) 
                && fPlanesInSet==5) fPlanesInSet -= 1;

        if(validPlaneItr.SizeSelect() > fPlanesInSet ){
          UseSlidingWindow(planeItr, validPlaneItr, windowPlaneItr, status, 
                           status->GetVertexPlaneNumber(),
                           status->GetEndPlaneNumber());        
        }
        else if(validPlaneItr.SizeSelect() <= fPlanesInSet){
          UseSingleFit(planeItr, validPlaneItr, status);
        }
        
        //check the end points of the event - send the validPlaneItr to the function
        //fStrayPlanes -= util.CheckFit(validPlaneItr);
        
        //what about multiple muons?  look for another vertex after the end plane
        planeItr.GetSet()->ClearSlice(false);
        validPlaneItr.GetSet()->ClearSlice(false);
        windowPlaneItr.GetSet()->ClearSlice(false);
        vertexPlaneItr.GetSet()->Slice(status->GetEndPlaneNumber()+1, 500);
        vertexPlaneItr.GetSet()->Update();

        //MSG("DMX1", Msg::kDebug)<< "event =  " << status->GetEventNumber() 
        //                   << "\tend plane = " << status->GetEndPlaneNumber() << endl;
        
        Int_t nextVertex = -1;
        if(vertexPlaneItr.SizeSelect() > 3){ nextVertex = util.FindVertexPlane(vertexPlaneItr);}
        
        while( nextVertex != -1){
          
          status->SetMultipleMuon(true); 
          cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent);
          vertexPlaneItr.GetSet()->ClearSlice(false);
          vertexPlaneItr.GetSet()->Slice(nextVertex, 500);
          vertexPlaneItr.GetSet()->Update();
          
          //get the next end plane
          Int_t endPlane = util.FindEndPlane(vertexPlaneItr);

          //MSG("DMX1", Msg::kDebug)<< "\tin double muon loop, vertex at plane " << nextVertex 
          //             << "\tend plane = " << endPlane << "\tplanes in slice = " 
          //             << vertexPlaneItr.SizeSelect() << endl;
          
          planeItr.GetSet()->Slice(nextVertex, endPlane);
          planeItr.GetSet()->Update();
          
          validPlaneItr.GetSet()->Slice(nextVertex, endPlane);
          validPlaneItr.GetSet()->Update();
          
          windowPlaneItr.GetSet()->Slice(nextVertex, endPlane);
          windowPlaneItr.GetSet()->Update();
                    
          //if only 2 planes set them to their best hyps
          if(validPlaneItr.SizeSelect() == 2){
            validPlaneItr.ResetFirst();
            
            //variables for straight line fit
            Float_t cogBeg = -1., zBeg = -1., cogEnd = -1., zEnd = -1.;
            
            while( validPlaneItr.IsValid() ){
              
              //only set it if it is a shower plane, muon planes are set in the constructor
              if(validPlaneItr.Ptr()->GetPlaneType() == DmxPlaneTypes::kShower
                 && !validPlaneItr.Ptr()->IsGolden())
                dynamic_cast<DmxShowerPlane *>(validPlaneItr.Ptr())->SetStrips("best");
              
              //set the variables for the straighe line fit between the two planes
              if(cogBeg == -1.){
                cogBeg = validPlaneItr.Ptr()->GetCoG();
                zBeg = validPlaneItr.Ptr()->GetZPosition();
              }
              else{
                cogEnd = validPlaneItr.Ptr()->GetCoG();
                zEnd = validPlaneItr.Ptr()->GetZPosition();
              }
              
              validPlaneItr.Next();
            }//end loop over planes
            
            if(planeItr.SizeSelect()>2 && cogEnd-cogBeg != 0. && cogBeg != -1. && cogEnd != -1.){
              Double_t slope = (cogEnd-cogBeg)/(zEnd - zBeg);
              Double_t intercept = cogBeg - slope*zBeg;
              SetPlanesToDeterminedFit(planeItr, intercept, slope);
              fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
            }
            
            //fStrayPlanes -= util.CheckFit(validPlaneItr);
            validPlaneItr.ResetFirst();
            
          }//end if only 2 valid planes
          else if(validPlaneItr.SizeSelect() > 2){
            //if there are only enough planes for 1 instances of the sliding window + initial fit,
            //change fPlanesInSet to give you 3 + the initial fit
            if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==6) fPlanesInSet -= 2;
            else if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) 
                    && fPlanesInSet==5) fPlanesInSet -= 1;
                    
            if(validPlaneItr.SizeSelect() > fPlanesInSet ){
              UseSlidingWindow(planeItr, validPlaneItr, windowPlaneItr, status, 
                               nextVertex, endPlane);   
            }
            else if(validPlaneItr.SizeSelect() <= fPlanesInSet ){
              UseSingleFit(planeItr, validPlaneItr, status);
            }
          }//end if more than 2 valid planes in event

          //check the end points of the event - send the validPlaneItr to the function
          //fStrayPlanes -= util.CheckFit(validPlaneItr);

          //look for another vertex after the end plane
          planeItr.GetSet()->ClearSlice(false);
          validPlaneItr.GetSet()->ClearSlice(false);
          windowPlaneItr.GetSet()->ClearSlice(false);
          vertexPlaneItr.GetSet()->ClearSlice(false);

          vertexPlaneItr.GetSet()->Slice(endPlane+1, 500);
          vertexPlaneItr.GetSet()->Update();

          //MSG("DMX1", Msg::kDebug)<< "\tlook for new vertex, slice from " << endPlane+1 
          //             << "\tto plane 500, planes in slice = " 
          //             << vertexPlaneItr.SizeSelect() << endl;
          
          vertexPlaneItr.Reset();
          
          if(vertexPlaneItr.SizeSelect() > 0){
            nextVertex = util.FindVertexPlane(vertexPlaneItr);
          }
          else{ nextVertex = -1; }
        }//end loop to find next vertex for spread out multiples
      }//end if more than 2 valid planes in set
      else if(validPlaneItr.SizeSelect() < 2){ 
        MSG("DMX", Msg::kDebug)<< "Event " << status->GetEventNumber() 
                               << " not demuxed; only " 
                               << validPlaneItr.SizeSelect() 
                               << " valid planes in view" << endl;

        status->SetValidPlanesFailure(true); 
        cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kTooFewValidPlanes);
      } //not enough planes to demux

      //done using the first view, so clear the selection function
      //MSG("Dmx", Msg::kDebug) <<"Done with View " << viewCtr << endl;
      
      status->SetFigureOfMerit(validPlaneItr.SizeSelect(), fStrayPlanes);
      if(viewCtr == 0 && status->GetEventDeMuxed() ){ 
        status->SetUStrayPlanes(fStrayPlanes);
        status->SetUValidPlanes(validPlaneItr.SizeSelect());
        cdlh.SetNumValidPlanesU(validPlaneItr.SizeSelect());
        cdlh.SetNumStrayPlanesU(fStrayPlanes);
        if( status->GetFigureOfMeritFailure() ){
          //see if the event failed the figure of Merit test, if so, see if it is an 
          //overlapping muon
          cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kEventFailedFilter);
          if(util.IsOverlappingMultiple(validPlaneItr, 1.*status->GetVertexPlaneNumber(), 
                                        acd.GetDouble("HoughInterceptRMSLimit"), 
                                        acd.GetDouble("HoughSlopeRMSLimit"),
                                        acd.GetDouble("HoughPeakLimit"))){
            status->SetUOverlappingMultiple(true);
            cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent);
          }
        }
      }
      if(viewCtr == 1 && status->GetEventDeMuxed() ){ 
        status->SetVStrayPlanes(fStrayPlanes);
        status->SetVValidPlanes(validPlaneItr.SizeSelect());
        cdlh.SetNumValidPlanesV(validPlaneItr.SizeSelect());
        cdlh.SetNumStrayPlanesV(fStrayPlanes);
        //see if the event failed the figure of Merit test, if so, see if it is an 
        //overlapping muon
        if( status->GetFigureOfMeritFailure() ){
          cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kEventFailedFilter);
          if(util.IsOverlappingMultiple(validPlaneItr, 1.*status->GetVertexPlaneNumber(), 
                                        acd.GetDouble("HoughInterceptRMSLimit"), 
                                        acd.GetDouble("HoughSlopeRMSLimit"),
                                        acd.GetDouble("HoughPeakLimit"))){
            status->SetVOverlappingMultiple(true);
            cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent);
          }
        }
      }
      vertexPlaneItr.GetSet()->ClearSlice(false);
      planeItr.GetSet()->ClearSlice(false);
      planeItr.GetSet()->AdoptSelectKeyFunc(0);
      planeItr.Reset();
      validPlaneItr.GetSet()->ClearSlice(false);
      validPlaneItr.GetSet()->AdoptSelectKeyFunc(0);
      validPlaneItr.Reset();
      windowPlaneItr.GetSet()->ClearSlice(false);
      windowPlaneItr.GetSet()->AdoptSelectKeyFunc(0);
      windowPlaneItr.Reset();

      ++viewCtr;
    } //end loop over views

    //if the event was demuxed, check to see how the demuxing and timing jive
    //otherwise see if you need to set the nonphysical solution flag
    if( status->GetEventDeMuxed() ){
      status->SetAverageTimingOffset(util.CheckFitWithTiming(validPlaneItr));
      cdlh.SetAvgTimeOffset(status->GetAverageTimingOffset());
    }
    else if(status->GetNonPhysicalFailure() )  
      cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNonPhysicalStripSolution);
         
    
  }//end if there was a vertex
  else{ 
    status->SetNoVertexFailure(true);
    cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNoVertex);
    MSG("DMX", Msg::kDebug) << "no vertex found for event " << status->GetEventNumber() << endl;
  }

  //get rid of the temporary status object if you made it.
  if(tempStatus) delete status;

  return;
}

//______________________________________________________________________
void AlgDeMuxCosmics::Trace(const char * /* c */) const
{
}


//----------------------------------------------------------------------------------
//this is the method used in the current algorithm
void AlgDeMuxCosmics::FindCosmicSolution(DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr,
                                         DmxStatus *status, Double_t &a1, Double_t &a2,
                                         Double_t &a3, Double_t &a4, Double_t &chiSq, 
                                         Float_t offset)
{
  //variables to keep track of the best reconstruction
  Double_t fitA1 = 0.;
  Double_t fitA2 = 0.;
  Double_t fitA3 = 0.;
  Double_t fitA4 = 0.;
  Double_t useA1 = -10000.;
  Double_t useA2 = -10000.;
  Double_t useA3 = -10000.;
  Double_t useA4 = -10000.;
  Double_t bestChi2 = 1.e20;
  Double_t chi2 = 1.000001e20;
  Int_t direction = status->GetEventDirection();

  DmxUtilities util;

  //make the lever arm for the demuxer the # of valid planes or 6, whichever
  //is smaller
  UInt_t leverArm = validPlaneItr.SizeSelect();
    
  if(leverArm > fPlanesInSet){ leverArm = fPlanesInSet;}
  
  //MSG("DMXX", Msg::kDebug) << "in FindCosmicSolution" << endl;

  //MSG("DMX", Msg::kDebug) << "\tlever arm = " << leverArm << endl;

  //declare an array to tell the fitter which hypotheses to use
  Int_t hypsToUse[] = {0,0,0,0,0,0};
  
  for(Int_t plane0Hyp = 1; plane0Hyp < 4; plane0Hyp++){
    hypsToUse[0] = plane0Hyp;
    for(Int_t plane1Hyp = 1; plane1Hyp < 4; plane1Hyp++){
      hypsToUse[1] = plane1Hyp;
      if(leverArm == 2){
        
        FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
        
        if(chi2 < bestChi2){
          bestChi2 = chi2;
          useA1 = fitA1;
          useA2 = fitA2;
          useA3 = fitA3;
          useA4 = fitA4;
        }//end if its a better straight line fit
      }
      else if(leverArm >= 3){
        for(Int_t plane2Hyp = 1; plane2Hyp < 4; plane2Hyp++){
          hypsToUse[2] = plane2Hyp;
          if(leverArm == 3){
            
            FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
            
            if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
              bestChi2 = chi2;
              useA1 = fitA1;
              useA2 = fitA2;
              useA3 = fitA3;
              useA4 = fitA4;

              //see if you can drop some planes and improve the fit
              if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
                FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
                
                if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
                  bestChi2 = chi2;
                  useA1 = fitA1;
                  useA2 = fitA2;
                  useA3 = fitA3;
                  useA4 = fitA4;
                }
                validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
              }//end if planes should be dropped because they are way off
            }//end if its a better straight line fit
          }
          else if(leverArm >= 4){
            for(Int_t plane3Hyp = 1; plane3Hyp < 4; plane3Hyp++){
              hypsToUse[3] = plane3Hyp;
              if(leverArm == 4){
                
                FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
                
                if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
                  bestChi2 = chi2;
                  useA1 = fitA1;
                  useA2 = fitA2;
                  useA3 = fitA3;
                  useA4 = fitA4;

                  //see if you can drop some planes and improve the fit
                  if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
                    FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
                    
                    if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
                      bestChi2 = chi2;
                      useA1 = fitA1;
                      useA2 = fitA2;
                      useA3 = fitA3;
                      useA4 = fitA4;
                    }
                    validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
                  }//end if planes should be dropped because they are way off
                }//end if its a better straight line fit
              }
              else if(leverArm >= 5){
                for(Int_t plane4Hyp = 1; plane4Hyp < 4; plane4Hyp++){
                  hypsToUse[4] = plane4Hyp;
                  
                  if(leverArm == 5){
                    
                    FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
                    
                    if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
                      bestChi2 = chi2;
                      useA1 = fitA1;
                      useA2 = fitA2;
                      useA3 = fitA3;
                      useA4 = fitA4;

                      //see if you can drop some planes and improve the fit
                      if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
                        FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
                        
                        if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
                          bestChi2 = chi2;
                          useA1 = fitA1;
                          useA2 = fitA2;
                          useA3 = fitA3;
                          useA4 = fitA4;
                        }
                        validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
                      }//end if planes should be dropped because they are way off
                    }//end if its a better straight line fit
                  }
                  else if(leverArm == 6){
                    for(Int_t plane5Hyp = 1; plane5Hyp < 4; plane5Hyp++){
                      
                      hypsToUse[5] = plane5Hyp;
                      FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
                      if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
                        bestChi2 = chi2;
                        useA1 = fitA1;
                        useA2 = fitA2;  
                        useA3 = fitA3;
                        useA4 = fitA4;
                     
                        //see if you can drop some planes and improve the fit
                        if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
                          FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
                          if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2,  fitA3, fitA4, offset) ){
                            bestChi2 = chi2;
                            useA1 = fitA1;
                            useA2 = fitA2;
                            useA3 = fitA3;
                            useA4 = fitA4;
                      
                          }
                          validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
                        } //end if planes should be dropped because they are way off
                      }//end if its a better straight line fit
                    }//end for loop over 6th plane's hyps
                  }//end if using 6 planes
                }//end for loop over 5th plane's hyps
              }//end if at least 5 planes
            } //loop over 4th plane
          } //end if at least 4 planes
        } //loop over 3rd plane
      } //end if at least 3 planes
    } //loop over 2nd plane
  } //loop over 1st plane
  
  //force the fit to the one found
  if(useA1 != -10000. && useA2 != -10000. && useA3 != 10000. && useA4 != 10000.){
    a1 = useA1;
    a2 = useA2;
    a3 = useA3;
    a4 = useA4;
    chiSq = bestChi2;
  }
  
  //MSG("DMX", Msg::kDebug) << "finished in FindCosmicsSolution" << endl;
  return;
}

//----------------------------------------------------------------------------------
void AlgDeMuxCosmics::FindWindowCosmicSolution(DmxPlaneItr &validPlaneItr, 
                                               DmxStatus *status, Double_t &a1, 
                                               Double_t &a2, Double_t &a3, 
                                               Double_t &a4,Double_t &chiSq, 
                                               Float_t offset)
                                               
{
  //MSG("DMXX", Msg::kDebug) << "in FindWindowCosmicSolution" << endl;

  //variables to keep track of the best reconstruction
  Double_t fitA1L = 0.;
  Double_t fitA2L = 0.;
  Double_t chi2L = 1.0e20;
//   Double_t fitA1Q = 0.;
//   Double_t fitA2Q = 0.;
//   Double_t fitA3Q = 0.;
//   Double_t chi2Q = 1.0e20;
  Double_t fitA1 = 0.;
  Double_t fitA2 = 0.;
  Double_t fitA3 = 0.;
  Double_t fitA4 = 0.;
  Double_t chi2 = 1.0e20;
  Double_t chi2Fit = 1.0e20;

  Int_t direction = status->GetEventDirection();

  //make the lever arm for the demuxer the # of valid planes or 6, whichever
  //is smaller
  Int_t leverArm = fPlanesInSet;

  //offset is the lowest z pos in the window. subtract it from each x array value
  //- this has the effect of changing your y axis location 
  //relative to the z = 0 location so that when doing higher order fits, you 
  //dont get really stupid numbers
 
  //declare the variables to do a least squares fit to a straight line.  
  //use three chiSq values to keep track of the 3 different hyps in the
  //unset plane - also have 3 different y values for those hyps and 3
  //different weights
  Double_t x[] = {0.,0.,0.,0.,0.,0.};
  Double_t y[] = {0.,0.,0.,0.,0.,0.};
  Double_t weight[] = {1.,1.,1.,1.,1.,1.};
  Double_t weight2 = 1.;
  Double_t weight3 = 1.;
  Double_t y2 = 0.;
  Double_t y3 = 0.;

  //variable to let you know if the unset plane is a golden plane or not
  bool lastGolden = false;

  Double_t a1L = 0.;
  Double_t a2L = 0.;
  Double_t a1LB = 0.;
  Double_t a2LB = 0.;
  Double_t a1LS = 0.;
  Double_t a2LS = 0.;
  Double_t a1LT = 0.;
  Double_t a2LT = 0.;
 //  Double_t a1QB = 0.;
//   Double_t a2QB = 0.;
//   Double_t a3QB = 0.;
//   Double_t a1QS = 0.;
//   Double_t a2QS = 0.;
//   Double_t a3QS = 0.;
//   Double_t a1QT = 0.;
//   Double_t a2QT = 0.;
//   Double_t a3QT = 0.;
//   Double_t a1Q = 0.;
//   Double_t a2Q = 0.;
//   Double_t a3Q = 0.;
  
  Double_t chiSqL = 1.e20;
  Double_t chiSqLB = 1.e20;
  Double_t chiSqLS = 1.e20;
  Double_t chiSqLT = 1.e20;
//   Double_t chiSqQB = 1.e20;
//   Double_t chiSqQS = 1.e20;
//   Double_t chiSqQT = 1.e20;
//   Double_t chiSqQ = 1.e20;

  Double_t a1Fit = 0.;
  Double_t a2Fit = 0.;
  Double_t a3Fit = 0.;
  Double_t a4Fit = 0.;

  
  DmxUtilities util;

  if(direction == -1) validPlaneItr.ResetLast();
  else if( direction == 1) validPlaneItr.ResetFirst();
  
  Int_t planeCtr = 0;
  DmxPlane *plane = 0;
      
  while( (plane = validPlaneItr()) && planeCtr < leverArm){

    //cout << "plane = " << plane->GetPlaneNumber() << endl;

    x[planeCtr] = plane->GetZPosition() - offset;
    y[planeCtr] = plane->GetCoG();
    weight[planeCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge());

    if(plane->IsGolden()){

      //golden plane, so make the weight smaller to give it more
      //clout in the fit
      weight[planeCtr] /= TMath::Sqrt(10.);
      if(planeCtr == leverArm-1){
        lastGolden = true;
        y2 = y[planeCtr];
        y3 = y[planeCtr];
        weight2 = weight[planeCtr];
        weight3 = weight[planeCtr];
      }
    }

    //only make changes if not a golden plane
    if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden() ){

      //make sure you dont divide by zero when getting the mated signal ratio of
      //the set hypothesis - could be that you set the plane to a hypothesis
      //with no mated signal, even though it is a valid plane
      if(planeCtr < leverArm-1 
         && dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio() > 0.){
        weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()); 
//      MSG("Dmx", Msg::kDebug) << dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio() << endl;
      }
      else if( planeCtr == leverArm-1 ){          
        
        y2 = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
        y3 = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();

        //dont have to worry about dividing by zero here because you havent
        //set the plane yet and you know that the best hypotheses in the planes
        //must have at least 0.3 mated signal
        weight2 = weight[planeCtr]/TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); 
        weight3 = weight[planeCtr]/TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());         

        weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
      }
    }

    ++planeCtr;
    
  } //end loop over planes to find variables for fit
      
  //do fit using function in DmxUtilities

  //the best hyp fit first
  util.FindLinearFit(x, y, weight, leverArm, a1LB, a2LB, chiSqLB);
  //util.FindQuadraticFit(x, y, weight, leverArm, a1QB, a2QB, a3QB, chiSqQB);

  //find fit for using the last plane's second best hyp
  y[leverArm-1] = y2;
  weight[leverArm-1] = weight2;
  util.FindLinearFit(x, y, weight, leverArm, a1LS, a2LS, chiSqLS);
  //util.FindQuadraticFit(x, y, weight, leverArm, a1QS, a2QS, a3QS, chiSqQS);

  //find fit using the last plane's third best hyp
  y[leverArm-1] = y3;
  weight[leverArm-1] = weight3;
  util.FindLinearFit(x, y, weight, leverArm, a1LT, a2LT, chiSqLT);
  //util.FindQuadraticFit(x, y, weight, leverArm, a1QT, a2QT, a3QT, chiSqQT);

  //util.FindCubicFit(x, y, weight, leverArm, a1, a2, a3, a4, chiSq);

 //  MSG("DMXX", Msg::kDebug) << "\tlinear fit1 " << chiSqLB << "\t" << a1LB << "\t" << a2LB << endl;
//   MSG("DMXX", Msg::kDebug) << "\tlinear fit2 " << chiSqLS << "\t" << a1LS << "\t" << a2LS << endl;
//   MSG("DMXX", Msg::kDebug) << "\tlinear fit3 " << chiSqLT << "\t" << a1LT << "\t" << a2LT << endl;
  //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
  //              << "\t" << a2Q << "\t" << a3Q << endl;
   
  //pick the better linear fit
  if(chiSqLB<=chiSqLS && chiSqLB<=chiSqLT){
    chiSqL = chiSqLB;
    a1L = a1LB;
    a2L = a2LB;
  }
  else if(chiSqLS<chiSqLB && chiSqLS<=chiSqLT){
    chiSqL = chiSqLS;
    a1L = a1LS;
    a2L = a2LS;
  }
  else if(chiSqLT<chiSqLB && chiSqLT<chiSqLS){
    chiSqL = chiSqLT;
    a1L = a1LT;
    a2L = a2LT;
  }

  //pick the better quadratic fit
  // if(chiSqQB<=chiSqQS && chiSqQB<=chiSqQT){
//     chiSqQ = chiSqQB;
//     a1Q = a1QB;
//     a2Q = a2QB;
//     a3Q = a3QB;
//   }
//   else if(chiSqQS<chiSqQB && chiSqQS<=chiSqQT){
//     chiSqQ = chiSqQS;
//     a1Q = a1QS;
//     a2Q = a2QS;
//     a3Q = a3QS;
//   }
//   else if(chiSqQT<chiSqQB && chiSqQT<chiSqQS){
//     chiSqQ = chiSqQT;
//     a1Q = a1QT;
//     a2Q = a2QT;
//     a3Q = a3QT;
//   }

  //pick the better fit - linear or quadratic
  //if(chiSqL<=chiSqQ){
    chi2Fit = chiSqL;
    a1Fit = a1L;
    a2Fit = a2L;
    a3Fit = 0.;
    a4Fit = 0.;
    //}
 //  else{
//     chi2Fit = chiSqQ;
//     a1Fit = a1Q;
//     a2Fit = a2Q;
//     a3Fit = a3Q;
//     a4Fit = 0.;
//   }

    //MSG("DMXX", Msg::kDebug) << "\tinitial fit " << chi2Fit << "\t" << a1Fit << "\t" << a2Fit 
    //            << "\t" << a3Fit << "\t" << a4Fit << endl;
  
  //redo the fit using only the previously set planes and only if you have more than 3 planes in the set
  if(leverArm > 3){
    
    //do fit using function in DmxUtilities
    util.FindLinearFit(x, y, weight, leverArm-1, fitA1L, fitA2L, chi2L);
    //util.FindQuadraticFit(x, y, weight, planeCtr, fitA1Q, fitA2Q, fitA3Q, chi2Q);
    //util.FindCubicFit(x, y, weight, planeCtr, fitA1, fitA2, fitA3, fitA4, chiSq);
    
    //MSG("DMXX", Msg::kDebug) << "\trefined linear fit " << chi2L << "\t" << fitA1L << "\t" << fitA2L << endl;
    //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
    //              << "\t" << a2Q << "\t" << a3Q << endl;
  
    //pick the better fit - linear or quadratic
    //if(chi2L<=chi2Q){
    chi2 = chi2L;
    fitA1 = fitA1L;
    fitA2 = fitA2L;
    fitA3 = 0.;
    fitA4 = 0.;
    //}
    //  else{
    //       chi2 = chi2Q;
    //       fitA1 = fitA1Q;
    //       fitA2 = fitA2Q;
    //       fitA3 = fitA3Q;
    //       fitA4 = 0.;
    //     }

    
    if( lastGolden ){

      if(direction == 1) validPlaneItr.ResetLast();
      else if( direction == -1) validPlaneItr.ResetFirst();
      
      plane = validPlaneItr();
      
      MSG("AlgDmx", Msg::kDebug) << "last plane number = " << plane->GetPlaneNumber() 
                                 << " this is a golden plane " 
                                 << plane->IsGolden() << endl;
     
      
      //unset the IsGolden flag for the last plane if the difference between the
      //fit center of gravity and the plane's center of gravity is > 0.5m and
      //any of the following
      //1. the golden CoG is different from the previous set valid plane's
      //   CoG by more than 0.25m/plane (plane to plane spacing is 0.1188m)
      //2. chi^2 of fit without the golden plane is an order of magnitude
      //   smaller than the chi^2 of the fit with it
      //3. the average difference between the set CoG's and the golden CoG and
      //   the fit CoG's for all planes in the window is greater than 0.5m

      //cout << x[leverArm-1] << " " << x[leverArm-2] << " " << a1Fit << " " << a2Fit << endl;

      Double_t fitCoG = (a1Fit + a2Fit*(x[leverArm-1])
                         + a3Fit*TMath::Power(x[leverArm-1],2) 
                         + a4Fit*TMath::Power(x[leverArm-1],3));

      Double_t fitDiff = 0.;
      for(Int_t ii = 0; ii<leverArm; ii++){
        fitDiff += TMath::Abs((a1Fit + a2Fit*(x[ii]) + a3Fit*TMath::Power(x[ii],2) 
                               + a4Fit*TMath::Power(x[ii],3) - y[ii]));
      }
      Double_t planeDiff = TMath::Abs(x[leverArm-1]-x[leverArm-2])/0.1188;
      
      MSG("DmxAlg", Msg::kDebug) << "fit is " << fitCoG << " golden CoG = " 
                                << plane->GetCoG() << endl
                                <<" y diff = " << TMath::Abs(y[leverArm-2]-plane->GetCoG()) 
                                << "/" << planeDiff*.25 << endl
                                <<" fitDiff = " << fitDiff <<"/" << 0.5*leverArm << endl
                                << "chi2Fit = " << chi2Fit << "/" << 10.*chi2 << endl;
      
      if(TMath::Abs(plane->GetCoG()-fitCoG)>0.5 && (fitDiff>0.5*leverArm
         || TMath::Abs(y[leverArm-2]-plane->GetCoG())>(0.25*planeDiff)
         || chi2Fit>10.*chi2)){
        
        plane->SetGolden(false);
        MSG("DmxAlg", Msg::kDebug) << plane->GetPlaneNumber() << " was golden " 
                                   << plane->IsGolden() << endl;
        lastGolden = false;
      }
    }//end check of last plane is golden

    //if the fit found by dropping the last plane has a better chi^2 than the 
    //fit using all planes, and that last plane was not golden, then pass the
    //fit from the first planes in the window out.  if the last plane was initially
    //golden, this only happens for the last time through.
    //
    if(chi2 < chi2Fit && !lastGolden){
    
      chi2Fit = chi2;
      a1Fit = fitA1;
      a2Fit = fitA2;
      a3Fit = fitA3;
      a4Fit = fitA4;
    }
   
  }//end for loop to improve fit by dropping a plane

  chiSq = chi2Fit;  
  a1 = a1Fit;
  a2 = a2Fit;
  a3 = a3Fit;
  a4 = a4Fit;
  
  //MSG("DMXX", Msg::kDebug) << "\tchosen linear fit " << chiSq << "\t" << a1 << "\t" << a2 << endl;

  validPlaneItr.Reset();

  return;
}

//------------------------------------------------------------------------------------------
//this method is used in the algorithm
void AlgDeMuxCosmics::UseSingleFit(DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, 
                                   DmxStatus *status)
{
  //MSG("DMXX", Msg::kDebug) << "in UseSingleFit" << endl;

  DmxUtilities util = DmxUtilities();
  
  Double_t a1F = -10000.;
  Double_t a2F = -10000.;
  Double_t a3F = -10000.;
  Double_t a4F = -10000.;
  Double_t chiSqF = 1.1e20;
  Double_t a1B = -10000.;
  Double_t a2B = -10000.;
  Double_t a3B = -10000.;
  Double_t a4B = -10000.;
  Double_t chiSqB = 1.1e20;
 
  validPlaneItr.ResetFirst();
  Float_t offset = validPlaneItr.Ptr()->GetZPosition();
  validPlaneItr.ResetLast();
  //set the event direction forwards
  status->SetEventDirection(1);
  
  //cout << "calling FindCosmicSolution" << endl;
  FindCosmicSolution(planeItr, validPlaneItr, status, a1F, a2F, a3F, a4F, chiSqF, offset);
  
  //set the event direction to backwards
  status->SetEventDirection(-1);
  
  FindCosmicSolution(planeItr, validPlaneItr, status, a1B, a2B, a3B, a4B, chiSqB, offset);
  
  DmxPlane *plane = 0;
  if( chiSqF <= chiSqB && a1F != -10000. && a2F != -10000. 
      && a3F != -10000. && a4F != 10000.){
    
    //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
    //set it as non-golden.
    validPlaneItr.Reset();

    while( (plane = planeItr()) ){
      Double_t fitCoG = (a1F + a2F*plane->GetZPosition() 
                         + a3F*TMath::Power(plane->GetZPosition(),2) 
                         + a4F*TMath::Power(plane->GetZPosition(),3));
      if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
    }

    SetPlanesToDeterminedFit(planeItr, a1F, a2F, a3F, a4F, offset); 
   
    fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
  }
  else if(chiSqB < chiSqF && a1B != -10000. && a2B != -10000. 
          && a3B != 10000. && a4B != 10000.){
   
    //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
    //set it as non-golden.
    validPlaneItr.Reset();
  
    while( (plane = planeItr()) ){
      Double_t fitCoG = (a1B + a2B*plane->GetZPosition() 
                         + a3B*TMath::Power(plane->GetZPosition(),2) 
                         + a4B*TMath::Power(plane->GetZPosition(),3));
      if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
    }

    SetPlanesToDeterminedFit(planeItr, a1B, a2B, a3B, a4B, offset); 
    
    fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
}
  else if( a1B == -10000. && a2B == -10000. && a3B != 10000. && a4B != 10000.
           && a1F == -10000. && a2F == -10000. && a3F != -10000. && a4F != 10000.){
    
    MSG("DMX", Msg::kDebug)<< "Event " << status->GetEventNumber() 
                           << " not demuxed; can't find a physical solution" 
                           << endl;
    status->SetNonPhysicalFailure(true); 
  }
  
  
  return;
}

//-------------------------------------------------------------------------------------------
//this method is the one used in the current algorithm
void AlgDeMuxCosmics::UseSlidingWindow(DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr,
                                       DmxPlaneItr &windowPlaneItr, DmxStatus *status, 
                                       Int_t vertex, Int_t endPlane)
{
  //MSG("DMXX", Msg::kDebug) << "in UseSlidingWindow" << endl;

  //start the sliding window - keep the window to 6 valid planes 
  //and use the fit parameters you found for the previous set.  if the 
  //new plane is a shower plane, set to the hypothesis nearest the fit.  
  //if none of the 3 best hypotheses is within 12 strips of the fit, then
  //just use the same fit parameters.  redo the fit, and then set all 
  //non-set planes in the window. move the the window down another plane 
  //and do it again.

  DmxUtilities util;

  UInt_t validPlaneCnt = validPlaneItr.SizeSelect();
  UInt_t loopCtr = 0;
  UInt_t ctr = 0;
  Int_t firstPlane = 0;
  Int_t almostLastPlane = 0;
  Int_t lastPlane = 0;
  Double_t a1 = -10000.;
  Double_t a2 = -10000.;
  Double_t a3 = -10000.;
  Double_t a4 = -10000.;
  Double_t chiSq = 1.1e20;
  Float_t offset = 0.;
  DmxPlane *plane = 0;
  Int_t direction = 1;  //status->GetEventDirection();

  // MSG("dmx", Msg::kDebug) << "Event = " << status->GetEventNumber() << "\tvertex = " << vertex 
//               << "\tend = " << endPlane << endl;
//   while(validPlaneItr.IsValid() ){
//     MSG("dmx", Msg::kDebug) << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber()
//                         << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->IsGolden()
//                         <<endl;
//     validPlaneItr.Next();
//   }

  validPlaneItr.ResetFirst();
  windowPlaneItr.ResetFirst();

  //now demux the event the <= makes sure that you find an initial fit for the first
  //set of planes, then go from there
  while(loopCtr <= (validPlaneCnt - fPlanesInSet) ){

      ctr = 0;
      
      //MSG("dmx", Msg::kDebug) << "find new window" << endl;
      //MSG("dmx", Msg::kDebug) << "\tdirection = " << direction;

      while( validPlaneItr.IsValid() && ctr < fPlanesInSet){
        
        plane = validPlaneItr.Ptr();
        //MSG("dmx", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
        //             << endl;
        
        if(direction == 1){ 
          if( ctr == 0){ 
            firstPlane = plane->GetPlaneNumber(); 
            offset = plane->GetZPosition();
          }
          if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); }
          if( ctr == fPlanesInSet-1){ lastPlane = plane->GetPlaneNumber(); }
          validPlaneItr.Next();
        }
        else if(direction == -1){ 
          if( ctr == 0){ lastPlane = plane->GetPlaneNumber(); }
          if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); }
          if( ctr == fPlanesInSet-1){ 
            firstPlane = plane->GetPlaneNumber(); 
            offset = plane->GetZPosition();
          }
          validPlaneItr.Prev();
        }
        ++ctr;
        
      }

      //MSG("dmx", Msg::kDebug) << "end find new window" << endl;
      ctr = 0;
            
      //MSG("dmx", Msg::kDebug) << "\tfirst plane = " << firstPlane << "\talmost last plane = " 
      //                     << almostLastPlane << "\tlast plane = " << lastPlane << endl;
      
      //move the window
      windowPlaneItr.GetSet()->ClearSlice(false);
      planeItr.GetSet()->ClearSlice(false);

      windowPlaneItr.GetSet()->Slice(firstPlane, lastPlane);
      windowPlaneItr.GetSet()->Update();
      planeItr.GetSet()->Slice(firstPlane, lastPlane);
      planeItr.GetSet()->Update();

      if( loopCtr == 0){
        
        //set the initial window - try it both ways and take the best one as
        //determined by the chiSq
        Double_t a1B = -10000.;
        Double_t a2B = -10000.;
        Double_t a3B = -10000.;
        Double_t a4B = -10000.;
        Double_t chiSqB = 1.1e20;
        
        //set the event direction forwards
        status->SetEventDirection(1);

        //you are looking at the first planes of the event, so your offset should be the
        //z position of the vertex
        Float_t offsetF = status->GetVertexZPosition();

        FindCosmicSolution(planeItr, windowPlaneItr, status, a1, a2, a3, a4, chiSq, offsetF);

        //MSG("dmx", Msg::kDebug) << "\tforward direction" << a1 << "\t" << a2 
        //             << "\t" << chiSq << "\t" << offset << endl;
        
        //set the event direction to backwards
        Int_t forwardLastPlane = lastPlane;

        status->SetEventDirection(-1);
        direction = -1;
        validPlaneItr.ResetLast();
        windowPlaneItr.GetSet()->ClearSlice(false);
        planeItr.GetSet()->ClearSlice(false);
        windowPlaneItr.ResetLast();
        ctr = 0;
        
        //MSG("dmx", Msg::kDebug) << "check reverse direction" << endl;
        while( validPlaneItr.IsValid() && ctr < fPlanesInSet){

          plane = validPlaneItr.Ptr();
          // MSG("dmx", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
          //              << endl;
          
          if( ctr == 0){ lastPlane = plane->GetPlaneNumber(); }
          if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); }
          if( ctr == fPlanesInSet-1){ 
            firstPlane = plane->GetPlaneNumber(); 
            offset = plane->GetZPosition();
          }
          ++ctr;
          validPlaneItr.Prev();
        }
        //MSG("dmx", Msg::kDebug) << "end check reverse direction" << endl;
        ctr = 0;
        
        //MSG("dmx", Msg::kDebug) << "first plane = " << firstPlane << "\talmost last plane = " 
        //              << almostLastPlane << "\tlast plane = " << lastPlane << endl;
                
        windowPlaneItr.GetSet()->Slice(firstPlane, lastPlane);
        windowPlaneItr.GetSet()->Update();
        windowPlaneItr.ResetLast();
        planeItr.GetSet()->Slice(firstPlane, lastPlane);
        planeItr.GetSet()->Update();

        FindCosmicSolution(planeItr, windowPlaneItr, status, a1B, a2B, a3B, a4B, chiSqB, offset);
        
        //MSG("dmx", Msg::kDebug) << "\tbackward direction" << a1B << "\t" << a2B 
        //             << "\t" << chiSqB << "\t" << offset << endl;
        
        if( chiSq <= chiSqB && a1 != -10000. && a2 != -10000. && a3 != 10000. && a4 != 10000.){
          planeItr.GetSet()->ClearSlice(false);

          //slice from the vertex to the last plane in the forward direction to get them
          //all
          planeItr.GetSet()->Slice(vertex, forwardLastPlane);
          planeItr.GetSet()->Update();

          //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
          //set it as non-golden.
          windowPlaneItr.GetSet()->ClearSlice(false);
          windowPlaneItr.GetSet()->Slice(vertex, forwardLastPlane);
          windowPlaneItr.GetSet()->Update();

          while( (plane = windowPlaneItr()) ){
            Double_t fitCoG = (a1 + a2*plane->GetZPosition() 
                               + a3*TMath::Power(plane->GetZPosition(),2) 
                               + a4*TMath::Power(plane->GetZPosition(),3));
            if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
          }

          //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << vertex 
          //     << "\tto " << forwardLastPlane << endl;

          SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offsetF); 

          fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
            
          status->SetEventDirection(1);
          direction = 1;
          validPlaneItr.ResetFirst();
          ctr = 0;
          
          //get the validPlaneItr back to where it should be for this direction
          while(ctr<fPlanesInSet){
            validPlaneItr.Next();
            ++ctr;
          }
          
          ctr = 0;
          
        }           
        else if(chiSqB < chiSq && a1B != -10000. && a2B != -10000. && a3B != 10000. && a4B != 10000.){
          
          //set a, b, and chiSq to the backwards values
          a1 = a1B;
          a2 = a2B;
          a3 = a3B;
          a4 = a4B;
          chiSq = chiSqB;

          //slice from the end plane to the first plane
          planeItr.GetSet()->ClearSlice(false);
          planeItr.GetSet()->Slice(firstPlane, endPlane);
          planeItr.GetSet()->Update();

          //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
          //set it as non-golden.
          windowPlaneItr.ResetFirst();
          while( (plane = windowPlaneItr()) ){
            Double_t fitCoG = (a1 + a2*plane->GetZPosition() 
                               + a3*TMath::Power(plane->GetZPosition(),2) 
                               + a4*TMath::Power(plane->GetZPosition(),3));
            if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
          }
          //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << firstPlane 
          //             << "\tto " << endPlane << endl;

          SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offset);

          fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);

          //event direction is already backwards and the validPlaneItr is where it should be
        }
        else if( a1B == -10000. && a2B == -10000. && a1 == -10000. 
                 && a2 == -10000. && a3 == -10000. && a4 == -10000. 
                 && a3B == -10000. && a4B == -10000.){
           
          MSG("dmx", Msg::kDebug)<< "Event " << status->GetEventNumber() 
                                 << " not demuxed;" 
                                 << " can't find a physical solution" << endl;
          status->SetNonPhysicalFailure(true);
        }
        //MSG("dmx", Msg::kDebug) << "\tloop 0 fit - " << a1 << "\t" << a2 << "\t" << a3 << "\tchi^2 = " 
        //             << chiSq << "\tstray planes = " << fStrayPlanes << endl;
      }//end if initial set
      else{
        //set the window
        
        FindWindowCosmicSolution(windowPlaneItr, status, a1, a2, a3, a4, chiSq, offset);

        //MSG("dmx", Msg::kDebug) << "\tloop " << loopCtr << " " << firstPlane << " " << lastPlane 
        //             << "\t" << a1 << "\t" << a2 << "\t" << a3 << "\tchi^2 = " 
        //             << chiSq << "\tstray planes = " << fStrayPlanes << endl;
        
        planeItr.GetSet()->ClearSlice(false);
        if(direction == 1 ){
          firstPlane = almostLastPlane + 1;
          //if its the last loop, select the rest of the unset planes to set, otherwise
          //stay in the loop bounds
          if(loopCtr == (validPlaneItr.SizeSelect()-fPlanesInSet)) lastPlane = endPlane;
        }
        else if(direction == -1 ){
          lastPlane = almostLastPlane-1;
          if(loopCtr == (validPlaneItr.SizeSelect()-fPlanesInSet)) firstPlane = vertex;
        }
        planeItr.GetSet()->Slice(firstPlane, lastPlane);
        planeItr.GetSet()->Update();

        //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << firstPlane 
        //             << "\tto " << lastPlane << endl;
        
        SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offset);
        fStrayPlanes += util.FindPlanesOffFit(planeItr,  fStrayCut);
        //status->SetHistogramEntries(b, a, plane->GetPlaneView());
      }


      ctr = 0;

      //set the validPlaneItr to the opposite direction and back up 5 planes
      //MSG("dmx", Msg::kDebug) << "move iterator to first plane in new window" << endl;
      while(ctr<fPlanesInSet-1 && validPlaneItr.IsValid()){
        
        //MSG("dmx", Msg::kDebug) << "plane number = " 
        //             << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber()
        //             << endl;
        ++ctr;
        if(direction == 1){ validPlaneItr.Prev(); }
        else if(direction == -1){ validPlaneItr.Next(); }
      }
      //MSG("dmx", Msg::kDebug) << "end move iterator to first plane in new window" << endl;
      ctr = 0;
 
      //should go through the loop and get the new window
      ++loopCtr;
  }

  return;
}

//-----------------------------------------------------------------------------
//this method is used in the current algorithm
//find the best line for a set of planes
void AlgDeMuxCosmics::FindFit(DmxPlaneItr &planeItr, Double_t &prevChiSq, 
                              Double_t &a1, Double_t &a2, Double_t &a3, 
                              Double_t &a4, Int_t* hypotheses, 
                              Int_t direction, Int_t leverArm, Float_t offset)
{

  //MSG("DMXX", Msg::kDebug) << "in FindFit" << endl;

  //cout << "in find fit for first n planes" << endl;
  Double_t a1L = 0.;
  Double_t a2L = 0.;
  Double_t a1Q = 0.;
  Double_t a2Q = 0.;
  Double_t a3Q = 0.;
  
  Double_t chiSqL = 1.e20;
  Double_t chiSqQ = 1.e20;

  //declare the variables to do a least squares fit to a straight line.  
  //weight is the inverse fraction of the total charge on the plane that 
  //is on opposite ends of common strips. multiply the inverse by a 
  //sensible number to take account for digits with 1000+ adc's

  //the arrays are the number of planes in the set, ie the lever arm
  Double_t x[] = {0.,0.,0.,0.,0.,0.};
  Double_t y[] = {0.,0.,0.,0.,0.,0.};
  Double_t weight[] = {1.,1.,1.,1.,1.,1.};

  DmxUtilities util;

  //cout << "got util object" << endl;
  if(direction == -1){ planeItr.ResetLast(); }
  else if( direction == 1){ planeItr.ResetFirst(); }

  //planeItr.ResetFirst();

  Int_t planeCtr = 0;
  Int_t arrayCtr = 0;
  DmxPlane *plane = 0;

  
  //MSG("DMX", Msg::kDebug) << "\t" << hypotheses[0] << "\t" << hypotheses[1] 
  //             << "\t" << hypotheses[2] << "\t" << hypotheses[3] 
  //             << "\t" << hypotheses[4] << "\t" << hypotheses[5] << endl;
  
  //loop over the plane iterator to fill the variables

  while( (plane = planeItr()) && planeCtr < leverArm){

    if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
      x[arrayCtr] = plane->GetZPosition() - offset;
     
      //get the maximum possible weight first, then for each hypothesis multiply by the square of the fraction
      weight[arrayCtr] = 1./ TMath::Sqrt(plane->GetPlaneCharge());
      if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.);
      
      y[arrayCtr] = plane->GetCoG();
      
      //only do this if the shower plane is not golden - if it is, you know where it should be
      //reconstructed to
      if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){

        if(hypotheses[planeCtr] == 1){
          //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetCoG();
          weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
        }
        //plane->GetCoG() returns the best hypothesis CoG so only change the value of y if you are wanting
        //the 2nd or 3rd best hypotheses
        else if(hypotheses[planeCtr] == 2){ 
          //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetCoG();
          y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
          weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
        }
        else if(hypotheses[planeCtr] == 3){ 
          //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetCoG();
          y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
          weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
        }
      }
      ++arrayCtr;
      
    }//end if the plane isnt marked as sucking && it is in the window bounds
    
    ++planeCtr;
  } //end loop over planes to find variables for fit

  //cout << "filled arrays" << endl;

  Double_t chiSq = 0.;
  
  //use the DmxUtilities function to find a linear fit
  util.FindLinearFit(x, y, weight, arrayCtr, a1L, a2L, chiSqL);
  //util.FindQuadraticFit(x, y, weight, arrayCtr, a1Q, a2Q, a3Q, chiSqQ);
  //util.FindCubicFit(x, y, weight, arrayCtr, a1, a2, a3, a4, chiSq);

  //MSG("DMXX", Msg::kDebug) << "\tlinear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl;
  //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
  //                      << "\t" << a2Q << "\t" << a3Q << endl;

  //pick the better fit
  if(chiSqL<=chiSqQ){
    prevChiSq = chiSqL;
    a1 = a1L;
    a2 = a2L;
    a3 = 0.;
    a4 = 0.;
  }
  else{
    prevChiSq = chiSqQ;
    a1 = a1Q;
    a2 = a2Q;
    a3 = a3Q;
    a4 = 0.;
  }
  //MSG("DMX", Msg::kDebug) << "\tinitial fit " << a1 << "\t" << a2 << "\t" << a3 << "\t" << a4 << chiSq << endl; 
 
  planeItr.Reset();
  planeCtr = 0;
  arrayCtr = 0;
  //cout << "reset the iterator" << endl;
  
  //redo the fit as many times as there are planes, dropping each plane in 
  //turn to see if it produces a better chi^2 value - only do it if you have 
  //more than 3 planes - ie you can drop one and still do a reasonable fit
  //2 planes = great straight line every time.
  if(leverArm > 3){
    a1L = 0.;
    a2L = 0.;
    a1Q = 0.;
    a2Q = 0.;
    a3Q = 0.;

    Double_t fitA1 = 0.;
    Double_t fitA2 = 0.;
    Double_t fitA3 = 0.;
    Double_t fitA4 = 0.;

    for(Int_t i = 0; i < leverArm; i++){
      arrayCtr = 0;

      //cout <<"in loop " << i << endl;
      //loop over the plane iterator to fill the variables
      while( (plane = planeItr()) && planeCtr < leverArm){
        
        //use those planes not marked as kTRUE - a mark of kTRUE means to drop the plane
        //from the set when finding the fit. ie if its true, the plane truely sucks
        if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){

          if(i != planeCtr){
            x[arrayCtr] = plane->GetZPosition() - offset;
            y[arrayCtr] = plane->GetCoG();

            //cout << "filling arrays " << arrayCtr << " " << x[arrayCtr] << " " << y[arrayCtr] 
            // << " " << weightSq[arrayCtr] << endl;

            weight[arrayCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge());
            if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.);
      
            //only make changes if not a golden plane
            if( plane->GetPlaneType() == DmxPlaneTypes::kShower  && !plane->IsGolden()){
              //only change y if the hypothesis isnt the best one - the unset plane returns the best
              //cog from the GetCoG() method
              
              if(hypotheses[planeCtr] == 1){
                weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio());
              }
              else if(hypotheses[planeCtr] == 2){ 
                y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
                weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
              }
              else if(hypotheses[planeCtr] == 3){ 
                y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
                weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
              }
            }

            ++arrayCtr;
            //else{ MSG("DMX", Msg::kDebug) << "\tdrop plane " << plane->GetPlaneNumber() << endl; }
          } //end if plane is not the dropped one
        }//end if plane is in window bounds
        ++planeCtr;
      } //end loop over planes to find variables for fit
            
      //use the DmxUtilities function to find a linear fit
      
      util.FindLinearFit(x, y, weight, arrayCtr, a1L, a2L, chiSqL);
      //util.FindQuadraticFit(x, y, weight, arrayCtr, a1Q, a2Q, a3Q, chiSqQ);
      //util.FindCubicFit(x, y, weight, arrayCtr, fitA1, fitA2, fitA3, fitA4, chiSq);
            
      //MSG("DMXX", Msg::kDebug) << "\tlinear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl;
      //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
      //                              << "\t" << a2Q << "\t" << a3Q << endl;
      //pick the better fit
      if(chiSqL<=chiSqQ){
        chiSq = chiSqL;
        fitA1 = a1L;
        fitA2 = a2L;
      }
      else{
        chiSq = chiSqQ;
        fitA1 = a1Q;
        fitA2 = a2Q;
        fitA3 = a3Q;
      }
    
      planeItr.Reset();

      planeCtr = 0;
      if(chiSq < prevChiSq){
        prevChiSq = chiSq;
        a1 = fitA1;
        a2 = fitA2;
        a3 = fitA3;
        a4 = fitA4;
        //MSG("DMX", Msg::kDebug) << "\trefined fit, drop plane " << i 
        //             << "\t" << prevChiSq << "\t" << a << "\t" << b << endl;
      }
      
    }//end for loop to improve fit by dropping a plane
  }//end if enough planes to improve fit

  planeItr.ResetFirst();
  //MSG("DMX", Msg::kDebug) << "\tfinal fit" << "\t" << a << "\t" << b << "\t" << prevChiSq << endl; 

  //cout << "finished find fit for first n planes" << endl;
  
  return;
}

//------------------------------------------------------------------------------------------
//function to set the planes in the event to the straight line fit
//this one is used in the current algorithm
void AlgDeMuxCosmics::SetPlanesToDeterminedFit(DmxPlaneItr &planeItr, Double_t a1, Double_t a2,
                                               Double_t a3, Double_t a4, Float_t offset)
{
  //MSG("DMXX", Msg::kDebug) << "in SetPlanesToDeterminedFit" << endl;

  planeItr.ResetFirst();

  DmxPlane *plane = 0;

  while( planeItr.IsValid() ){

    plane = planeItr.Ptr();
    
    //MSG("DMXX", Msg::kDebug)<<"\tin SetPlanesToDeterminedFit " << plane->GetPlaneNumber() 
    //             <<" has golden flag " << (Int_t)plane->IsGolden() << endl;
    
    Float_t fitCoG = (a1 + (plane->GetZPosition()-offset)*a2 
                      + TMath::Power(plane->GetZPosition()-offset,2)*a3
                      + TMath::Power(plane->GetZPosition()-offset,3)*a4);

    //MSG("DMXX", Msg::kDebug)<< a1 << "\t" << a2 << "\t" << a3 
    //             << "\t" << a4 << "\t" << fitCoG 
    //             << "\t" << plane->GetZPosition() << endl;
    
    //golden planes have already been set in their constructors
    if( !plane->IsGolden() ){
      //MSG("DMXX", Msg::kDebug)<<"\tsetting plane " << plane->GetPlaneNumber() 
      //                     << " to " << fitCoG << endl;
      plane->SetStrips(fitCoG); 
    }
        
    planeItr.Next();
  }

  planeItr.Reset();
  
  //cout << "finished set planes for fit" << endl;
  return;
}

//------------------------------------------------------------------------------------------
//function to set the planes in the event to the straight line fit
//used for cases with only 2 valid planes in a view
void AlgDeMuxCosmics::SetPlanesToDeterminedFit(DmxPlaneItr &planeItr, Double_t a, Double_t b)
{

  //MSG("DMXX", Msg::kDebug) << "in SetPlanesToDeterminedFit" << endl;

  planeItr.ResetFirst();

  DmxPlane *plane = 0;

  while( planeItr.IsValid() ){

    plane = planeItr.Ptr();
      
    //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
    //MSG("DMXX", Msg::kDebug)<<"\tin SetPlanesToDeterminedFit " << plane->GetPlaneNumber() 
    //             <<" has golden flag " << (Int_t)plane->IsGolden() << endl;
    Float_t fitCoG = a + (plane->GetZPosition() * b);
    //MSG("DMXX", Msg::kDebug)<< a << "\t" << b << "\t" << fitCoG << endl;
    
    //golden planes have already been set in their constructors
    if( !plane->IsGolden() ){
      //MSG("DMXX", Msg::kDebug)<<"\tsetting plane " << plane->GetPlaneNumber() << endl;
      plane->SetStrips(fitCoG); 
    }
    //}//end if plane is in the window bounds
    
    planeItr.Next();
  }

  planeItr.Reset();
  
  return;
}

//------------------------------------------------------------------------------
//function to identify planes whose 3 best hypotheses are all way off from the fit line
//make a cut for them to be at least 12 strips off. only pass this an iterator over valid 
//planes
Bool_t AlgDeMuxCosmics::FindPlanesToDropFromFit(DmxPlaneItr &planeItr, Double_t a1, Double_t a2,
                                                Double_t a3, Double_t a4,
                                                Int_t direction, Int_t leverArm, Float_t offset)
{

  //MSG("DMXX", Msg::kDebug) << "in FindPlanesToDropFromFit" << endl;

  Double_t diff1 = 0.;
  Double_t diff2 = 0.;
  Double_t diff3 = 0.;
  bool planesDropped = false;

  Int_t dropPlanes = 0;
  Int_t planeCtr = 0;
  if( direction == 1 ){ planeItr.ResetFirst(); }
  else if( direction == -1 ){ planeItr.ResetLast(); }
 
  //MSG("DMX", Msg::kDebug)<<"\tin FindPlanesToDropFromFit"; 
  
  DmxPlane *plane = 0;
  //find the plane farthest off from the fit, keep going until you have just 3 planes left
  while( (plane = planeItr()) && dropPlanes<(Int_t)(0.1*leverArm) && planeCtr<leverArm){

    Double_t fitCoG = (a1 + (plane->GetZPosition()-offset)*a2 
                       + TMath::Power(plane->GetZPosition()-offset,2)*a3
                       + TMath::Power(plane->GetZPosition()-offset,3)*a4);
  
      //MSG("DMX", Msg::kDebug)<<"\t" << plane->GetPlaneNumber() << endl;
      
      if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
        if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){
          diff1 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG());
          diff2 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG());
          diff3 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG());
        }
        else if( plane->GetPlaneType() == DmxPlaneTypes::kMuon ){
          diff1 = TMath::Abs(fitCoG - plane->GetCoG());
          diff2 = TMath::Abs(fitCoG - plane->GetCoG());
          diff3 = TMath::Abs(fitCoG - plane->GetCoG());
        }
        
        //if all the 3 best hypotheses are way off from the fit, that plane is most 
        //likely messing up the fit.  so mark it as kTRUE - ie it, it truely sucks
        if(diff1 > .504 && diff2 > .504 && diff3 > .504 && !plane->IsGolden() ){
          planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr());
          planesDropped = true;
          //MSG("DMX", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl;
          //decrement the number of planes now used for a fit.
          ++dropPlanes;
        }
        else if(diff1 > 0.504 && plane->IsGolden()){
          
          //if this was a golden plane but it just doesnt work, make a non-golden plane
          planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr());
          plane->SetGolden(false);
          planesDropped = true;
          //MSG("DMX", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl;
          //decrement the number of planes now used for a fit.
          ++dropPlanes;
        }//end if dropping a golden plane
      } //end if the mask hasnt already been set for this plane
      
      ++planeCtr;
  }
  planeItr.ResetFirst();
  
  //cout << "finished drop planes from fit" << endl;

  return planesDropped;
}

//----------------------------------------------------------------------------------
// void AlgDeMuxCosmics::FindWinCosmicSolution(DmxPlaneItr &validPlaneItr, DmxStatus *status, 
//                                          Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4,
//                                          Double_t &chiSq, Float_t offset)
                                               
// {
//   MSG("DMXX", Msg::kDebug) << "in FindWinCosmicSolution" << endl;

//   //variables to keep track of the best reconstruction
//   Double_t fitA1 = 0.;
//   Double_t fitA2 = 0.;
//   Double_t fitA3 = 0.;
//   Double_t fitA4 = 0.;
//   Double_t useA1 = -10000.;
//   Double_t useA2 = -10000.;
//   Double_t useA3 = -10000.;
//   Double_t useA4 = -10000.;
//   Double_t bestChi2 = 1.0e20;
//   Double_t chi2 = 1.0000001e20;
//   Int_t direction = status->GetEventDirection();

//   //make the lever arm for the demuxer the # of valid planes or 6, whichever
//   //is smaller
//   Int_t leverArm = fPlanesInSet;

//   for(Int_t unSetPlaneHyp = 1; unSetPlaneHyp < 4; unSetPlaneHyp++){
        
//     FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, 
//          unSetPlaneHyp, direction, leverArm, offset);
    
//     if(chi2 < bestChi2 ){
//       bestChi2 = chi2;
//       useA1 = fitA1;
//       useA2 = fitA2;
//       useA3 = fitA3;
//       useA4 = fitA4;
//     }//end if its a better straight line fit
//   }
  
//   //check to make sure that the fit parameters aren't the defaults and 
//   //that the bestChi2 fit for the window isn't horrible compared to the previous
//   //fit
//   if(useA1 != -10000. && useA2 != -10000. && useA3 != -10000. 
//      && useA4 != -10000. ){//&& bestChi2<10.*chiSq){
//     a1 = useA1;
//     a2 = useA2;
//     a3 = useA3;
//     a4 = useA4;
//     chiSq = bestChi2;
//   }

//   //MSG("DMXX", Msg::kDebug) << "\tchosen linear fit " << chiSq << "\t" << a1 << "\t" << a2 << endl;

//   return;
// }

//---------------------------------------------------------------------------------------
//find the best line for a set of planes
// void AlgDeMuxCosmics::FindFit(DmxPlaneItr &planeItr, Double_t &prevChiSq, 
//                            Double_t &a1, Double_t &a2, Double_t &a3,
//                            Double_t &a4, Int_t unSetPlaneHyp,
//                            Int_t direction, Int_t leverArm, Float_t offset)
// {

//   MSG("DMXX", Msg::kDebug) << "in FindFit" << endl;

//   //offset is the lowest z pos in the window. subtract it from each x array value
//   //- this has the effect of changing your y axis location 
//   //relative to the z = 0 location so that when doing higher order fits, you 
//   //dont get really stupid numbers
 
//   //cout << "in find fit for window" << endl;

//   //declare the variables to do a least squares fit to a straight line.  
  
//   Double_t x[] = {0.,0.,0.,0.,0.,0.};
//   Double_t y[] = {0.,0.,0.,0.,0.,0.};
//   Double_t weight[] = {1.,1.,1.,1.,1.,1.};
//   Double_t chiSq = 0.;

//   Double_t a1L = 0.;
//   Double_t a2L = 0.;
//   Double_t a1Q = 0.;
//   Double_t a2Q = 0.;
//   Double_t a3Q = 0.;
  
//   Double_t chiSqL = 1.e20;
//   Double_t chiSqQ = 1.e20;

//   DmxUtilities util;

//   if(direction == -1){ planeItr.ResetLast(); }
//   else if( direction == 1){ planeItr.ResetFirst(); }
  
//   Int_t planeCtr = 0;
//   DmxPlane *plane = 0;
      
//   while( (plane = planeItr()) && planeCtr < leverArm){

//     //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//       x[planeCtr] = plane->GetZPosition() - offset;
//       y[planeCtr] = plane->GetCoG();
//       weight[planeCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge());
//       if(plane->IsGolden()) weight[planeCtr] /= TMath::Sqrt(10.);
      
//       //only make changes if not a golden plane
//       if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
//      if( planeCtr < leverArm-1 ){
//        weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()); 
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 1 ){          
//        weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 2 ){
//        y[planeCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
//        weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); 
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 3 ){
//        y[planeCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
//        weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());      
//      }
//       }
    
     
//       ++planeCtr;

//       //}//end if plane is in bounds for the window.
//   } //end loop over planes to find variables for fit
  
//   planeItr.Reset();
//   planeCtr = 0;
    
//   //do fit using function in DmxUtilities
//   util.FindLinearFit(x, y, weight, leverArm, a1L, a2L, chiSqL);
//   //util.FindQuadraticFit(x, y, weight, leverArm, a1Q, a2Q, a3Q, chiSqQ);
//   //util.FindCubicFit(x, y, weight, leverArm, a1, a2, a3, a4, chiSq);

//   MSG("DMXX", Msg::kDebug) << "\tlinear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl;
//   //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
//   //           << "\t" << a2Q << "\t" << a3Q << endl;
   
//   //pick the better fit
//   if(chiSqL<=chiSqQ){
//     prevChiSq = chiSqL;
//     a1 = a1L;
//     a2 = a2L;
//     a3 = 0.;
//     a4 = 0.;
//   }
//   else{
//     prevChiSq = chiSqQ;
//     a1 = a1Q;
//     a2 = a2Q;
//     a3 = a3Q;
//     a4 = 0.;
//   }

//   //MSG("DMXX", Msg::kDebug) << "\tinitial fit " << prevChiSq << "\t" << a1 << "\t" << a2 
//   //           << "\t" << a3 << "\t" << a4 << endl;
  
//   //redo the fit using only the previously set planes and only if you have more than 3 planes in the set
//   if(leverArm > 3){
//     a1L = 0.;
//     a2L = 0.;
//     a1Q = 0.;
//     a2Q = 0.;
//     a3Q = 0.;
 
//     Double_t fitA1 = 0.;
//     Double_t fitA2 = 0.;
//     Double_t fitA3 = 0.;
//     Double_t fitA4 = 0.;
   
//     //loop over the plane iterator to fill the variables
//     //use all but the last plane in the set to see what your fit is without it
//     while( (plane = planeItr()) && planeCtr < (leverArm - 1)){
      
//       x[planeCtr] = plane->GetZPosition() - offset;
//       y[planeCtr] = plane->GetCoG();
//       //weight[planeCtr] = 100. / plane->GetPlaneCharge();
//       //if(plane->IsGolden() ) weight[planeCtr]  = 0.001; // /= TMath::Sqrt(10.);
//       weight[planeCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge());
//       if(plane->IsGolden() ) weight[planeCtr] /= TMath::Sqrt(10.);
      
//       if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
//      weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio());
//       }
      
//       ++planeCtr;
      
//     } //end loop over planes to find variables for fit
    
    
//     planeItr.Reset();
       
//     //do fit using function in DmxUtilities
//     util.FindLinearFit(x, y, weight, planeCtr, a1L, a2L, chiSqL);
//     //util.FindQuadraticFit(x, y, weight, planeCtr, a1Q, a2Q, a3Q, chiSqQ);
//     //util.FindCubicFit(x, y, weight, planeCtr, fitA1, fitA2, fitA3, fitA4, chiSq);
    
//     MSG("DMXX", Msg::kDebug) << "\trefined linear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl;
//     //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
//     //                   << "\t" << a2Q << "\t" << a3Q << endl;
//     //pick the better fit
//     if(chiSqL<=chiSqQ){
//       chiSq = chiSqL;
//       fitA1 = a1L;
//       fitA2 = a2L;
//     }
//     else{
//       chiSq = chiSqQ;
//       fitA1 = a1Q;
//       fitA2 = a2Q;
//       fitA3 = a3Q;
//     }
    
//     //MSG("DMXX", Msg::kDebug) << "\tdrop fit " << chiSq << "\t" << fitA1 << "\t" << fitA2 
//     //                   << "\t" << fitA3 << "\t" << fitA4 << endl << "\t\t" << x[0] << "\t"
//     //                   << x[planeCtr] << "\t" << planeCtr << endl;
//     planeCtr = 0;

//     bool droppedGolden = false;

//     while( (plane = planeItr()) ){
      
//       if( planeCtr == (leverArm - 1) && plane->IsGolden() && unSetPlaneHyp == 3){
//      droppedGolden = true;

//      //if this is the 3rd time through the loop and the difference in the golden CoG and the
//      //fitCoG without the golden plane is greater than 0.25m, and the fit is inside the detector
//      //then set the plane as non-Golden it just fooled the golden cut beforee
//      //or, go with the better fit if it is way better - ie the chi^2 is a factor of 10
//      //or more less than the fit with all planes, golden or not. again, only do that 
//      //the last time through

//      Double_t fitCoG = (fitA1 + fitA2*(plane->GetZPosition()-offset)
//                         + fitA3*TMath::Power(plane->GetZPosition()-offset,2) 
//                         + fitA4*TMath::Power(plane->GetZPosition()-offset,3));

//      if(TMath::Abs(fitCoG)<=4.0
//         && (TMath::Abs(plane->GetCoG()-fitCoG)>0.25 || prevChiSq>10.*chiSq)){
//        plane->SetGolden(false);
//        //MSG("DmxAlg", Msg::kDebug) << " is golden " << (Int_t)plane->IsGolden() << endl;
//        droppedGolden = false;
//      }
//       }//end check of last plane is golden
//       ++planeCtr;
        
//     } //end loop to find chiSq
    
//     planeCtr = 0;
//     planeItr.Reset();

//     //if the fit found by dropping the last plane has a better chi^2 than the 
//     //fit using all planes, and that last plane was not golden, then pass the
//     //fit from the first planes in the window out.  if the last plane was initially
//     //golden, this only happens for the last time through.
//     //
//     if(chiSq < prevChiSq && !droppedGolden){
    
//       prevChiSq = chiSq;
//       a1 = fitA1;
//       a2 = fitA2;
//       a3 = fitA3;
//       a4 = fitA4;
//     }
   
//     planeItr.Reset();
    
//   }//end for loop to improve fit by dropping a plane

//   planeItr.ResetFirst();

//   //cout << "finished find fit for window" << endl;

//   return;
// }

//-------------------------------------------------------------------------------------
//function to identify planes whose 3 best hypotheses are all way off from the fit line
//make a cut for them to be at least 12 strips off. only pass this an iterator over valid 
//planes
// Bool_t AlgDeMuxCosmics::FindPlanesToDropFromFit(DmxPlaneItr &planeItr, Double_t a, Double_t b, 
//                                              Int_t direction, Int_t leverArm, 
//                                                 Int_t /* firstPlane */, Int_t /* lastPlane */)
// {
//   MSG("DMXX", Msg::kDebug) << "in FindPlanesToDropFromFit" << endl;

//   Double_t diff1 = 0.;
//   Double_t diff2 = 0.;
//   Double_t diff3 = 0.;
//   bool planesDropped = false;

//   Int_t dropPlanes = 0;
//   Int_t planeCtr = 0;
//   if( direction == 1 ){ planeItr.ResetFirst(); }
//   else if( direction == -1 ){ planeItr.ResetLast(); }
 
//   //MSG("DMX", Msg::kDebug)<<"\tin FindPlanesToDropFromFit"; 
  
//   DmxPlane *plane = 0;
//   //find the plane farthest off from the fit, keep going until you have just 3 planes left
//   while( (plane = planeItr()) && dropPlanes<(Int_t)(0.1*leverArm) && planeCtr<leverArm){

//     //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//       Double_t fitCoG = a + ((plane->GetZPosition() * 1.0) * b);
//       //MSG("DMX", Msg::kDebug)<<"\t" << plane->GetPlaneNumber() << endl;
      
//       if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
//      if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){
//        diff1 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG());
//        diff2 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG());
//        diff3 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG());
//      }
//      else if( plane->GetPlaneType() == DmxPlaneTypes::kMuon ){
//        diff1 = TMath::Abs(fitCoG - plane->GetCoG());
//        diff2 = TMath::Abs(fitCoG - plane->GetCoG());
//        diff3 = TMath::Abs(fitCoG - plane->GetCoG());
//      }
        
//      //if all the 3 best hypotheses are way off from the fit, that plane is most 
//      //likely messing up the fit.  so mark it as kTRUE - ie it, it truely sucks
//      if(diff1 > .504 && diff2 > .504 && diff3 > .504 && !plane->IsGolden() ){
//        planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr());
//        planesDropped = true;
//        //MSG("DMX", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl;
//        //decrement the number of planes now used for a fit.
//        ++dropPlanes;
//      }
//      else if(diff1 > 0.504 && plane->IsGolden()){
          
//        //if this was a golden plane but it just doesnt work, make a non-golden plane
//        planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr());
//        plane->SetGolden(false);
//        planesDropped = true;
//        //MSG("DMX", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl;
//        //decrement the number of planes now used for a fit.
//        ++dropPlanes;
//      }//end if dropping a golden plane
//       } //end if the mask hasnt already been set for this plane
      
//       ++planeCtr;
//       //}//end if plane is in window bounds
//   }
//   planeItr.ResetFirst();
  
//   return planesDropped;
// }

//----------------------------------------------------------------------------------
// void AlgDeMuxCosmics::FindCosmicSolution(DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr,
//                                       DmxStatus *status, Double_t &a, Double_t &b,
//                                       Double_t &chiSq, Int_t firstPlane, Int_t lastPlane)
// {
//   MSG("DMXX", Msg::kDebug) << "in FindCosmicSolution" << endl;

//   //variables to keep track of the best reconstruction
//   Double_t fitA = 0.;
//   Double_t fitB = 0.;
//   Double_t useA = -10000.;
//   Double_t useB = -10000.;
//   Double_t bestChi2 = 1.e20;
//   Double_t chi2 = 1.0000001e20;
//   Int_t direction = status->GetEventDirection();

//   DmxUtilities util;
//   //make the lever arm for the demuxer the # of valid planes or 6, whichever
//   //is smaller
//   UInt_t leverArm = validPlaneItr.SizeSelect();
    
//   if(leverArm > fPlanesInSet){ leverArm = fPlanesInSet;}
  
//   //MSG("DMX", Msg::kDebug) << "\tlever arm = " << leverArm << endl;

//   //declare an array to tell the fitter which hypotheses to use
//   Int_t hypsToUse[] = {0,0,0,0,0,0};
  
//   for(Int_t plane0Hyp = 1; plane0Hyp < 4; plane0Hyp++){
//     hypsToUse[0] = plane0Hyp;
//     for(Int_t plane1Hyp = 1; plane1Hyp < 4; plane1Hyp++){
//       hypsToUse[1] = plane1Hyp;
//       if(leverArm == 2){
        
//      FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
        
//      if(chi2 < bestChi2){
//        bestChi2 = chi2;
//        useA = fitA;
//        useB = fitB;
//      }//end if its a better straight line fit
//       }
//       else if(leverArm >= 3){
//      for(Int_t plane2Hyp = 1; plane2Hyp < 4; plane2Hyp++){
//        hypsToUse[2] = plane2Hyp;
//        if(leverArm == 3){
            
//          FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
            
//          if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//            bestChi2 = chi2;
//            useA = fitA;
//            useB = fitB;
              
//            //see if you can drop some planes and improve the fit
//            if( FindPlanesToDropFromFit(validPlaneItr, fitA, fitB, direction, leverArm, firstPlane, lastPlane) ){
//              FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
                
//              if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//                bestChi2 = chi2;
//                useA = fitA;
//                useB = fitB;
//              }
//              validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
//            }//end if planes should be dropped because they are way off
//          }//end if its a better straight line fit
//        }
//        else if(leverArm >= 4){
//          for(Int_t plane3Hyp = 1; plane3Hyp < 4; plane3Hyp++){
//            hypsToUse[3] = plane3Hyp;
//            if(leverArm == 4){
                
//              FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
                
//              if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//                bestChi2 = chi2;
//                useA = fitA;
//                useB = fitB;
                  
//                //see if you can drop some planes and improve the fit
//                if( FindPlanesToDropFromFit(validPlaneItr, fitA, fitB, direction, leverArm, firstPlane, lastPlane) ){
//                  FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
                    
//                  if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//                    bestChi2 = chi2;
//                    useA = fitA;
//                    useB = fitB;
//                  }
//                  validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
//                }//end if planes should be dropped because they are way off
//              }//end if its a better straight line fit
//            }
//            else if(leverArm >= 5){
//              for(Int_t plane4Hyp = 1; plane4Hyp < 4; plane4Hyp++){
//                hypsToUse[4] = plane4Hyp;
                  
//                if(leverArm == 5){
                    
//                  FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
                    
//                  if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//                    bestChi2 = chi2;
//                    useA = fitA;
//                    useB = fitB;
                      
//                    //see if you can drop some planes and improve the fit
//                    if( FindPlanesToDropFromFit(validPlaneItr, fitA, fitB, direction, leverArm, firstPlane, lastPlane) ){
//                      FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
                        
//                      if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//                        bestChi2 = chi2;
//                        useA = fitA;
//                        useB = fitB;
//                      }
//                      validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
//                    }//end if planes should be dropped because they are way off
//                  }//end if its a better straight line fit
//                }
//                else if(leverArm == 6){
//                  for(Int_t plane5Hyp = 1; plane5Hyp < 4; plane5Hyp++){
                      
//                    hypsToUse[5] = plane5Hyp;
//                    FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
//                    if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//                      bestChi2 = chi2;
//                      useA = fitA;
//                      useB = fitB;
//                      //MSG("DMX", Msg::kDebug) << "\tfindstraightlinefit" << useA << "\t" << useB << "\t" << chi2 << endl;
                     
//                      //see if you can drop some planes and improve the fit
//                      if( FindPlanesToDropFromFit(validPlaneItr, fitA, fitB, direction, leverArm, firstPlane, lastPlane) ){
//                        FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, hypsToUse, direction, leverArm, firstPlane, lastPlane);
//                        if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA, fitB, firstPlane, lastPlane) ){
//                          bestChi2 = chi2;
//                          useA = fitA;
//                          useB = fitB;
//                          //MSG("DMX", Msg::kDebug) << "\tfindstraightlinedropfit" << useA << "\t" << useB << "\t" << chi2 << endl;
                      
//                        }
//                        validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
//                      } //end if planes should be dropped because they are way off
//                    }//end if its a better straight line fit
//                  }//end for loop over 6th plane's hyps
//                }//end if using 6 planes
//              }//end for loop over 5th plane's hyps
//            }//end if at least 5 planes
//          } //loop over 4th plane
//        } //end if at least 4 planes
//      } //loop over 3rd plane
//       } //end if at least 3 planes
//     } //loop over 2nd plane
//   } //loop over 1st plane
  
//   //force the fit to the one found
//   if(useA != -10000. && useB != -10000.){
//     a = useA;
//     b = useB;
//     chiSq = bestChi2;

//     //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
//     //set it as non-golden.
//     validPlaneItr.Reset();
//     DmxPlane *plane = 0;

//     while( (plane = dynamic_cast<DmxPlane *>(planeItr())) ){
//       if(plane->IsGolden() && TMath::Abs(a+(b*plane->GetZPosition()))>0.25) plane->SetGolden(false);
//     }
//     validPlaneItr.Reset();

//   }
  
//   //MSG("DMX", Msg::kDebug) << "\t" << a << "\t" << b << "\t" << chiSq << endl;
//   return;
// }

//----------------------------------------------------------------------------------
// void AlgDeMuxCosmics::FindWindowCosmicSolution(DmxPlaneItr &validPlaneItr, DmxStatus *status, 
//                                             Double_t &a, Double_t &b, Double_t &chiSq, Int_t firstPlane, 
//                                             Int_t lastPlane)
                                               
// {
//   MSG("DMXX", Msg::kDebug) << "in FindWindowCosmicSolution" << endl;

//   //variables to keep track of the best reconstruction
//   Double_t fitA = 0.;
//   Double_t fitB = 0.;
//   Double_t useA = -10000.;
//   Double_t useB = -10000.;
//   Double_t bestChi2 = 1.e20;
//   Double_t chi2 = 1.000001e20;
//   Int_t direction = status->GetEventDirection();

//   //make the lever arm for the demuxer the # of valid planes or 6, whichever
//   //is smaller
//   Int_t leverArm = fPlanesInSet;

//   for(Int_t unSetPlaneHyp = 1; unSetPlaneHyp < 4; unSetPlaneHyp++){
        
//      FindStraightLineFit(validPlaneItr, chi2, fitA, fitB, unSetPlaneHyp, direction, leverArm, firstPlane, lastPlane);
        
//      if(chi2 < bestChi2 ){
//        bestChi2 = chi2;
//        useA = fitA;
//        useB = fitB;
//      }//end if its a better straight line fit
//   }
  
//   a = useA;
//   b = useB;
//   chiSq = bestChi2;

//   return;
// }

//------------------------------------------------------------------------------------
//find the best line for a set of planes
// void AlgDeMuxCosmics::FindStraightLineFit(DmxPlaneItr &planeItr, Double_t &prevChiSq, 
//                                        Double_t &a, Double_t &b, Int_t* hypotheses, 
//                                        Int_t direction, Int_t leverArm,
//                                           Int_t /* firstPlane */, Int_t /* lastPlane */)
// {
//    MSG("DMXX", Msg::kDebug) << "in FindStraightLineFit" << endl;


//   //declare the variables to do a least squares fit to a straight line.  all
//   //equations taken from Bevington, second edition.  fit to the line
//   //CoG = a + b*plane
//   //where CoG is the center of gravity for plane number plane.
//   //the variables used are
//   // pSqDivSSqSum = sum of (plane number)^2 / (weight)^2
//   // cogDivSSqSum = sum of (center of gravity)^2 / (weight)^2
//   // pDivSSqSum = sum of plane number / (weight)^2
//   // pCoGDivSSqSum = sum of (plane number * center of gravity) / (weight)^2
//   // invSSqSum = sum of weight^-2
//   // pDivSSqSum = sum of plane number / weight^2
//   //weight is the inverse fraction of the total charge on the plane that is on opposite ends of common strips
//   //multiply the inverse by a sensible number to take account for digits with 1000+ adc's

//   Double_t pSqDivSSqSum = 0.;
//   Double_t cogDivSSqSum = 0.;
//   Double_t pDivSSqSum = 0.;
//   Double_t pCoGDivSSqSum = 0.;
//   Double_t invSSqSum = 0.;

//   if(direction == -1){ planeItr.ResetLast(); }
//   else if( direction == 1){ planeItr.ResetFirst(); }

//   //planeItr.ResetFirst();

//   Int_t planeCtr = 0;
//   DmxPlane *plane = 0;
      
//   //MSG("DMX", Msg::kDebug) << "\t" << hypotheses[0] << "\t" << hypotheses[1] 
//   //          << "\t" << hypotheses[2] << "\t" << hypotheses[3] 
//   //          << "\t" << hypotheses[4] << "\t" << hypotheses[5] << endl;
//   //loop over the plane iterator to fill the variables

//   while( (plane = planeItr()) && planeCtr < leverArm){

//     if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
//       //&& plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//       Double_t x = plane->GetZPosition();

//       //get the maximum possible weight first, then for each hypothesis multiply 
//       //by the square of the fraction
//       Double_t weightSq = 10000. / (plane->GetPlaneCharge() * plane->GetPlaneCharge());
//       if(plane->IsGolden()) weightSq = 0.000001;
      
//       Double_t y = plane->GetCoG();
      
//       //only do this if the shower plane is not golden - if it is, you know where it should be
//       //reconstructed to
//       if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){

//      if(hypotheses[planeCtr] == 1){
//        //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetCoG();
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
//      }
//      //plane->GetCoG() returns the best hypothesis CoG so only change the value of y if you are wanting
//      //the 2nd or 3rd best hypotheses
//      else if(hypotheses[planeCtr] == 2){ 
//        //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetCoG();
//        y = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
//      }
//      else if(hypotheses[planeCtr] == 3){ 
//        //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetCoG();
//        y = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
//      }
//       }
      
//       pSqDivSSqSum += (x * x) / weightSq;
//       cogDivSSqSum += y / weightSq;
//       pDivSSqSum += x / weightSq;
//       pCoGDivSSqSum += (x * y) / weightSq;
//       invSSqSum += 1. / weightSq;
//       //MSG("DMX", Msg::kDebug) << "\t" << x << "\t" << y << "\t" << hypotheses[planeCtr] << "\t" << weightSq 
//       //                  << "\t" << pSqDivSSqSum << "\t" << cogDivSSqSum << "\t" 
//       //                  << pDivSSqSum << "\t" << pCoGDivSSqSum << "\t" << invSSqSum << endl;
  
//     }//end if the plane isnt marked as sucking && it is in the window bounds
    
//     ++planeCtr;
//   } //end loop over planes to find variables for fit
  
//   Int_t degreesOfFreedom = planeCtr-1;
//   planeItr.Reset();
//   planeCtr = 0;
  
//   //find the values of a and b
//   Double_t delta = (invSSqSum * pSqDivSSqSum) - (pDivSSqSum * pDivSSqSum);
//   a = (1. / delta) * ((pSqDivSSqSum * cogDivSSqSum) - (pDivSSqSum * pCoGDivSSqSum));
//   b = (1. / delta) * ((invSSqSum * pCoGDivSSqSum) - (pDivSSqSum * cogDivSSqSum));
//   //MSG("DMX", Msg::kDebug) << "\t" << delta << "\t" << a << "\t" << b << endl; 
      
//   //find the chi^2 value using the typical equation for chi^2 -
//   //chi^2 = SUM[(y_i - y(x_i)) / sigma_i]^2

//   Double_t chiSq = 0.;
  
//   while( (plane = planeItr()) && planeCtr < leverArm){
    
//     if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
//      //&& plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//       Double_t xi = plane->GetZPosition();
//       Double_t weight = 100. / plane->GetPlaneCharge();
//       Double_t yi = plane->GetCoG();
//       if(plane->IsGolden()) weight = 0.001;
      
//       //only change the value if the plane isnot golden
//       if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){

//      if(hypotheses[planeCtr] == 1){
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio(); 
//      }
//      //plane->GetCoG() returns the best hypothesis CoG so only change the value of y if you are wanting
//      //the 2nd or 3rd best hypotheses
//      else if(hypotheses[planeCtr] == 2){ 
//        yi = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
//      }
//      else if(hypotheses[planeCtr] == 3){ 
//        yi = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
//      }
//       }      

//       Double_t yxi = a + (b * xi);
//       Double_t quotient = (yi - yxi) / weight;
      
//       chiSq += quotient * quotient;
    
//     }//end if plane is in window bounds
  
//     ++planeCtr;
//   } //end loop to find chiSq
  
  
//   prevChiSq = chiSq/(1.*degreesOfFreedom);
//   //MSG("DMX", Msg::kDebug) << "\tinitial fit " << prevChiSq << "\t" << a << "\t" << b << endl;
    
//   planeItr.Reset();
//   planeCtr = 0;

//   //redo the fit as many times as there are planes, dropping each plane in 
//   //turn to see if it produces a better chi^2 value - only do it if you have 
//   //more than 3 planes - ie you can drop one and still do a reasonable fit
//   //2 planes = great straight line every time.
//   if(leverArm > 3){
//     Double_t fitA = 0.;
//     Double_t fitB = 0.;
    
//     for(Int_t i = 0; i < leverArm; i++){
//       planeCtr = 0;
//       pSqDivSSqSum = 0;
//       cogDivSSqSum = 0;
//       pDivSSqSum = 0;
//       pCoGDivSSqSum = 0;
//       invSSqSum = 0;
      
//       //loop over the plane iterator to fill the variables
//       while( (plane = planeItr()) && planeCtr < leverArm){
        
//      //use those planes not marked as kTRUE - a mark of kTRUE means to drop the plane
//      //from the set when finding the fit. ie if its true, the plane truely sucks
//      if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
//          //&& plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//        if(i != planeCtr){
//          Double_t x = plane->GetZPosition();
//          Double_t weightSq = 10000. / (plane->GetPlaneCharge() * plane->GetPlaneCharge());
//          Double_t y = plane->GetCoG();
//          if(plane->IsGolden()) weightSq = 0.000001;
      
//          //only make changes if not a golden plane
//          if( plane->GetPlaneType() == DmxPlaneTypes::kShower  && !plane->IsGolden()){
//            //only change y if the hypothesis isnt the best one - the unset plane returns the best
//            //cog from the GetCoG() method
              
//            if(hypotheses[planeCtr] == 1){
//              weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio() * 
//                                dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio());
//            }
//            else if(hypotheses[planeCtr] == 2){ 
//              y = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
//              weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()
//                                * dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
//            }
//            else if(hypotheses[planeCtr] == 3){ 
//              y = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
//              weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio()
//                                * dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
//            }
//          }
          
//          pSqDivSSqSum += (x * x) / weightSq;
//          cogDivSSqSum += y / weightSq;
//          pDivSSqSum += x / weightSq;
//          pCoGDivSSqSum += (x * y) / weightSq;
//          invSSqSum += 1. / weightSq;
//          //MSG("DMX", Msg::kDebug) << "\t" << x << "\t" << y << "\t" << hypotheses[planeCtr] << endl; //"\t" << weightSq << "\t" 
//          //               << pSqDivSSqSum << "\t" << cogDivSSqSum << "\t" 
//          //               << pDivSSqSum << "\t" << pCoGDivSSqSum << "\t" << invSSqSum << endl;
        
//        //else{ MSG("DMX", Msg::kDebug) << "\tdrop plane " << plane->GetPlaneNumber() << endl; }
//        } //end if plane is not the dropped one
//      }//end if plane is in window bounds
//         ++planeCtr;
//       } //end loop over planes to find variables for fit
      
//       degreesOfFreedom = planeCtr-2;
//       //find the values of a and b
//       delta = (invSSqSum * pSqDivSSqSum) - (pDivSSqSum * pDivSSqSum);
//       fitA = (1. / delta) * ((pSqDivSSqSum * cogDivSSqSum) - (pDivSSqSum * pCoGDivSSqSum));
//       fitB = (1. / delta) * ((invSSqSum * pCoGDivSSqSum) - (pDivSSqSum * cogDivSSqSum));
//       //MSG("DMX", Msg::kDebug) << "\tfit " << "\t" << fitA << "\t" << fitB << endl;
  
//       //find the chi^2 value using the typical equation for chi^2 -
//       //chi^2 = SUM[(y_i - y(x_i)) / sigma_i]^2
    
//       planeItr.Reset();
      
//       chiSq = 0.;
//       planeCtr = 0;
      
//       while( (plane = planeItr()) && planeCtr < leverArm){
//      //MSG("DmxAlg" , Msg::kDebug) << "plane " << plane->GetPlaneNumber() << " mask " 
//      //                 << (Int_t)planeItr.GetSet()->GetMasks().GetMask(plane)
//      //                 << endl;
//      if(!planeItr.GetSet()->GetMasks().GetMask(plane)){
//        //&& plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//        if( i != planeCtr ){
//          Double_t xi = plane->GetZPosition();
//          Double_t weight = 100. / plane->GetPlaneCharge();
//          Double_t yi = plane->GetCoG();
//          if(plane->IsGolden()) weight = 0.001;
          
//          //only make changes if not a golden plane
//          if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
//            if(hypotheses[planeCtr] == 1){
//              weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
//            }
//            //plane->GetCoG() returns the best hypothesis CoG so only change the value of y if you are wanting
//            //the 2nd or 3rd best hypotheses
//            else if(hypotheses[planeCtr] == 2){ 
//              yi = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
//              weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
//            }
//            else if(hypotheses[planeCtr] == 3){ 
//              yi = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
//              weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
//            }
//          }      
          
//          Double_t yxi = fitA + (fitB * xi);
//          Double_t quotient = (yi - yxi) / weight;
//          //add the contribution to the chiSq if its not the plane you are dropping
//          chiSq += quotient * quotient;
//        } //end if plane is not the one to be dropped
//       }//end if plane is in bounds of the window
//         ++planeCtr;
//       } //end loop to find chiSq
                     
//       planeItr.Reset();

//       planeCtr = 0;
//       chiSq /= 1.*degreesOfFreedom;
//       if(chiSq < prevChiSq){
//      prevChiSq = chiSq;
//      a = fitA;
//      b = fitB;
//      //MSG("DMX", Msg::kDebug) << "\trefined fit, drop plane " << i 
//      //             << "\t" << prevChiSq << "\t" << a << "\t" << b << endl;
//       }
      
//     }//end for loop to improve fit by dropping a plane
//   }//end if enough planes to improve fit

//   planeItr.ResetFirst();
//   //MSG("DMX", Msg::kDebug) << "\tfinal fit" << "\t" << a << "\t" << b << "\t" << prevChiSq << endl; 
  
//   return;
// }

//---------------------------------------------------------------------------------------
//find the best line for a set of planes
// void AlgDeMuxCosmics::FindStraightLineFit(DmxPlaneItr &planeItr, Double_t &prevChiSq, 
//                                        Double_t &a, Double_t &b, Int_t unSetPlaneHyp,
//                                        Int_t direction, Int_t leverArm, 
//                                           Int_t /* firstPlane */, Int_t /* lastPlane */)
// {

//   MSG("DMXX", Msg::kDebug) << "in FindStraightLineFit" << endl;
  
//   //declare the variables to do a least squares fit to a straight line.  all
//   //equations taken from Bevington, second edition.  fit to the line
//   //CoG = a + b*plane
//   //where CoG is the center of gravity for plane number plane.
//   //the variables used are
//   // pSqDivSSqSum = sum of (plane number)^2 / (weight)^2
//   // cogDivSSqSum = sum of (center of gravity)^2 / (weight)^2
//   // pDivSSqSum = sum of plane number / (weight)^2
//   // pCoGDivSSqSum = sum of (plane number * center of gravity) / (weight)^2
//   // invSSqSum = sum of weight^-2
//   // pDivSSqSum = sum of plane number / weight^2
//   //weight is the total charge on the plane
  
//   Double_t pSqDivSSqSum = 0.;
//   Double_t cogDivSSqSum = 0.;
//   Double_t pDivSSqSum = 0.;
//   Double_t pCoGDivSSqSum = 0.;
//   Double_t invSSqSum = 0.;

//   if(direction == -1){ planeItr.ResetLast(); }
//   else if( direction == 1){ planeItr.ResetFirst(); }
  
//   Int_t planeCtr = 0;
//   DmxPlane *plane = 0;
      
//   while( (plane = planeItr()) && planeCtr < leverArm){

//     //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//       Double_t x = plane->GetZPosition();
//       Double_t weightSq = 10000. / (plane->GetPlaneCharge() * plane->GetPlaneCharge());
//       Double_t y = plane->GetCoG();
//       if(plane->IsGolden()) weightSq = 0.000001;
      
//       //only make changes if not a golden plane
//       if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
//      if( planeCtr < leverArm-1 ){
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()); 
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 1 ){          
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 2 ){
//        y = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); 
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 3 ){
//        y = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());  
//      }
//       }
    
//       pSqDivSSqSum += (x * x) / weightSq;
//       cogDivSSqSum += y / weightSq;
//       pDivSSqSum += x / weightSq;
//       pCoGDivSSqSum += (x * y) / weightSq;
//       invSSqSum += 1. / weightSq;
//       //MSG("DMX", Msg::kDebug) << "\t" << x << "\t" << y << endl; //<< "\t" << weightSq << "\t" 
//       //<< pSqDivSSqSum << "\t" << cogDivSSqSum << "\t" 
//       //<< pDivSSqSum << "\t" << pCoGDivSSqSum << "\t" << invSSqSum << endl;
    
//       ++planeCtr;

//       //}//end if plane is in bounds for the window.
//   } //end loop over planes to find variables for fit
  
//   Int_t degreesOfFreedom = planeCtr-1;
//   planeItr.Reset();
//   planeCtr = 0;
  
//   //find the values of a and b
//   Double_t delta = (invSSqSum * pSqDivSSqSum) - (pDivSSqSum * pDivSSqSum);
//   a = (1. / delta) * ((pSqDivSSqSum * cogDivSSqSum) - (pDivSSqSum * pCoGDivSSqSum));
//   b = (1. / delta) * ((invSSqSum * pCoGDivSSqSum) - (pDivSSqSum * cogDivSSqSum));
//   //MSG("DMX", Msg::kDebug) << "\t" << delta << "\t" << a << "\t" << b << endl; 
      
//   //find the chi^2 value using the typical equation for chi^2 -
//   //chi^2 = SUM[(y_i - y(x_i)) / sigma_i]^2

//   Double_t chiSq = 0.;
  
//   while( (plane = planeItr()) && planeCtr < leverArm){

//     //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//       Double_t xi = plane->GetZPosition();
//       Double_t weight = 100. / plane->GetPlaneCharge();
//       Double_t yi = plane->GetCoG();
//       if(plane->IsGolden()) weight = 0.001;
      
//       if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
//      if( planeCtr < leverArm-1 ){
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 1 ){          
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 2 ){
//        yi = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
//      }
//      else if( planeCtr == leverArm-1 && unSetPlaneHyp == 3 ){
//        yi = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
//      }
//       }

//       Double_t yxi = a + (b * xi);
//       Double_t quotient = (yi - yxi) / weight;
      
//       chiSq += quotient * quotient;
      
//       ++planeCtr;
//       //}//end if plane is in window bounds
//   } //end loop to find chiSq
  
  
//   prevChiSq = chiSq/(1.*degreesOfFreedom);

//   //MSG("DMXX", Msg::kDebug) << "\tinitial fit " << prevChiSq << "\t" << a << "\t" << b << endl;
  
//   planeItr.Reset();
//   planeCtr = 0;

//   //redo the fit using only the previously set planes and only if you have more than 3 planes in the set
//   if(leverArm > 3){
//     Double_t fitA = 0.;
//     Double_t fitB = 0.;
    
//     planeCtr = 0;
//     pSqDivSSqSum = 0;
//     cogDivSSqSum = 0;
//     pDivSSqSum = 0;
//     pCoGDivSSqSum = 0;
//     invSSqSum = 0;
    
//     //loop over the plane iterator to fill the variables
//     while( (plane = planeItr()) && planeCtr < (leverArm - 1)){
        
//       //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//      //use all but the last plane in the set to see what your fit is without it
//      Double_t x = plane->GetZPosition();
//      Double_t weightSq = 10000. / (plane->GetPlaneCharge() * plane->GetPlaneCharge());
//      Double_t y = plane->GetCoG();
//      if(plane->IsGolden()) weightSq = 0.000001;
      
//      if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
//        weightSq *= 1. / (dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()
//                          * dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio());
//      }
        
//      pSqDivSSqSum += (x * x) / weightSq;
//      cogDivSSqSum += y / weightSq;
//      pDivSSqSum += x / weightSq;
//      pCoGDivSSqSum += (x * y) / weightSq;
//      invSSqSum += 1. / weightSq;
//      //MSG("DMX", Msg::kDebug) << "\t" << x << "\t" << y << endl; //"\t" << weightSq << "\t" 
//      //                   << pSqDivSSqSum << "\t" << cogDivSSqSum << "\t" 
//      //                   << pDivSSqSum << "\t" << pCoGDivSSqSum << "\t" << invSSqSum << endl;
        
//      //else{ MSG("DMX", Msg::kDebug) << "\tdrop plane " << plane->GetPlaneNumber() << endl; }
        
//      ++planeCtr;
//      //}//end if plane is in the window bounds
//     } //end loop over planes to find variables for fit
    
//     degreesOfFreedom = planeCtr-2;
//     //find the values of a and b
//     delta = (invSSqSum * pSqDivSSqSum) - (pDivSSqSum * pDivSSqSum);
//     fitA = (1. / delta) * ((pSqDivSSqSum * cogDivSSqSum) - (pDivSSqSum * pCoGDivSSqSum));
//     fitB = (1. / delta) * ((invSSqSum * pCoGDivSSqSum) - (pDivSSqSum * cogDivSSqSum));
    
//     //find the chi^2 value using the typical equation for chi^2 -
//     //chi^2 = SUM[(y_i - y(x_i)) / sigma_i]^2
    
//     planeItr.Reset();
    
//     chiSq = 0.;
//     planeCtr = 0;
    
//     bool droppedGolden = false;

//     while( (plane = planeItr()) ){
      
//       //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
//       if( planeCtr < (leverArm - 1)){
//      Double_t xi = plane->GetZPosition();
//      Double_t weight = 100. / plane->GetPlaneCharge();
//      Double_t yi = plane->GetCoG();
//      if(plane->IsGolden()) weight = 0.001;
      
//      if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
//        weight *= 1. / dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
//      }
        
//      Double_t yxi = fitA + (fitB * xi);
//      Double_t quotient = (yi - yxi) / weight;
//      //add the contribution to the chiSq if its not the plane you are dropping
//      chiSq += quotient * quotient;
//       }
//       //see if the dropped plane was a golden plane
//       else if( planeCtr == (leverArm - 1) && plane->IsGolden() && unSetPlaneHyp == 3){
//      droppedGolden = true;

//      //MSG("DmxAlg", Msg::kDebug) << "dropped plane " << plane->GetPlaneNumber() << " is golden " 
//      //                << (Int_t)plane->IsGolden() << "\tplane CoG = " 
//      //                << plane->GetCoG() << "fit CoG = " 
//      //                << fitA+(fitB*plane->GetZPosition()) 
//      //                << " a = " << fitA << " b = " << fitB 
//      //                << " chi^2 = " << chiSq << endl;
        
//      //if this is the 3rd time through the loop and the difference in the golden CoG and the
//      //fitCoG without the golden plane is greater than 0.25m, and the fit is inside the detector
//      //then set the plane as non-Golden it just fooled the golden cut beforee
//      //or, go with the better fit if it is way better - ie the chi^2 is a factor of 10
//      //or more less than the fit with all planes, golden or not. again, only do that 
//      //the last time through
//      if(TMath::Abs(fitA+(fitB*plane->GetZPosition()))<4.0
//         && (TMath::Abs(plane->GetCoG()-(fitA+(fitB*plane->GetZPosition())))>0.25 || prevChiSq>10.*chiSq/degreesOfFreedom)){
//        plane->SetGolden(false);
//        //MSG("DmxAlg", Msg::kDebug) << " is golden " << (Int_t)plane->IsGolden() << endl;
//        droppedGolden = false;
//      }
//       }//end check of last plane is golden
//       ++planeCtr;
//      //}//end if plane is in window bounds
//     } //end loop to find chiSq
    
//     planeCtr = 0;
//     chiSq /= 1.*degreesOfFreedom;
//     //MSG("DMXX", Msg::kDebug) << "\trefined fit "  
//     //                  << "\t" << chiSq << "\t" << fitA << "\t" << fitB << endl;

//     planeItr.Reset();

//     //if the fit found by dropping the last plane has a better chi^2 than the 
//     //fit using all planes, and that last plane was not golden, then pass the
//     //fit from the first planes in the window out.  if the last plane was initially
//     //golden, this only happens for the last time through.
//     //
//     if(chiSq < prevChiSq && !droppedGolden){
    
//       prevChiSq = chiSq;
//       a = fitA;
//       b = fitB;
//     }
   
//     planeItr.Reset();
    
//   }//end for loop to improve fit by dropping a plane

//   planeItr.ResetFirst();
//   return;
// }

---