DmxUtilities.cxx

Go to the documentation of this file.

//$Id: DmxUtilities.cxx,v 1.36 2004/11/03 22:22:49 brebel Exp $
//
//DmxUtilities.cxx
//
//DmxUtilities is the class that sorts digits into plane objects and 
//keeps track of all the info about an event you want to stick
//around - ie TObjArray of planes, vertex and end planes, etc
//
//Author:  B. Rebel 10/2001

#include "DeMux/DmxUtilities.h"
#include "MessageService/MsgService.h"
#include "MessageService/MsgFormat.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "DeMux/DmxHypothesis.h"
#include "DeMux/DmxShowerPlane.h"
#include "DeMux/DmxMuonPlane.h"
#include "DeMux/DmxStatus.h"
#include "TMath.h"
#include "TH2.h"
#include "TMatrix.h"
#include "Conventions/Munits.h"
#include "Conventions/PlaneView.h"
#include "Conventions/StripEnd.h"
#include "Algorithm/AlgConfig.h"
#include "CandDigit/CandDeMuxDigit.h"
#include "CandDigit/CandDeMuxDigitHandle.h"

//define the NavKey to sort the digits by plane
static NavKey KeyFromPlane( const CandDeMuxDigitHandle *cdh){
  return cdh->GetPlexSEIdAltL().GetPlane();}


ClassImp(DmxUtilities)
//______________________________________________________________________
CVSID("$Id: DmxUtilities.cxx,v 1.36 2004/11/03 22:22:49 brebel Exp $");

//---------------------------------------------------------------------
DmxUtilities::DmxUtilities()
{
}

//---------------------------------------------------------------------
void DmxUtilities::FillPlaneArray(DmxStatus *status,  
                                                                  CandDeMuxDigitListHandle &cdlh, AlgConfig &acd) 
{
                
        TObjArray *planeArray = new TObjArray();
        MSG("DmxUtil", Msg::kDebug) << cdlh.GetNDaughters() << " digits in DmxUtil" << endl;

        //make two iterators - one over the whole set and one to iterate over each plane
        //plane in turn
        CandDeMuxDigitHandleItr fCdhit = cdlh.GetDaughterIterator();
        
        //create the KeyFunc
        CandDeMuxDigitHandleKeyFunc *planeKF = fCdhit.CreateKeyFunc();
        
        //program the KeyFunc with the sort function
        planeKF->SetFun(KeyFromPlane);
        
        //get the NavSet from the iterator and pass the KeyFunc to it
        fCdhit.GetSet()->AdoptSortKeyFunc(planeKF);
        
        //clear the KF pointer because we no longer own the KeyFunc
        planeKF = 0;

        //initialize variable to keep track of the current plane number
        Int_t planeNum = 0;

        //MSG("DmxUtil", Msg::kDebug) << "current plane = " << planeNum << endl;

        CandDeMuxDigitHandle *currentDigit = 0;
        CandDeMuxDigitHandle *digit = 0;

        //go through the sorted Digits
        while( (currentDigit = fCdhit()) ){

                Int_t plane = 0;

                //find the plane number corresponding to that Digit
                plane = currentDigit->GetPlexSEIdAltL().GetPlane();
    
                //MSG("DmxUtil", Msg::kDebug) << "plane number of current digit = " << plane 
                //                  << endl;
    
                //if the plane number of the current digit is different from the
                //current plane number, then add another DmxPlaneDebug object to 
                //the TObjArray container for that plane number and increment the
                //fNumberOfPlanes variable
    
                if(plane != planeNum && !currentDigit->GetPlexSEIdAltL().IsVetoShield()){
      
                        planeNum = plane;

                        //make a new iterator for this plane's digit.  dont move this outside the
                        //loop because then you have one iterator for all planes versus 1 iterator
                        //per plane.  makes a big difference.
                        CandDeMuxDigitHandleItr planeCdhit(cdlh.GetDaughterIterator());
                        CandDeMuxDigitHandleKeyFunc *plane1KF = planeCdhit.CreateKeyFunc();
                        plane1KF->SetFun(KeyFromPlane);
                        planeCdhit.GetSet()->AdoptSortKeyFunc(plane1KF,true,false);
                        plane1KF = 0;

                        //get the slice for this plane
                        planeCdhit.GetSet()->Slice(planeNum);
                        planeCdhit.GetSet()->Update();
                        
                        Int_t numDigits = 0;
                        if( acd.GetInt("UseCandDigitMask")==1 ){
                                //loop over the digits in this plane and figure out which ones might be 
                                //from cross talk.  set the cutoff value for signal fraction to be in 
                                //AlgConfigDeMux
                                Float_t eastSet0[] = {0., 0., 0., 0., 0., 0., 0., 0., 
                                                                          0., 0., 0., 0., 0., 0., 0., 0.};
                                Float_t westSet0[] = {0., 0., 0., 0., 0., 0., 0., 0., 
                                                  0., 0., 0., 0., 0., 0., 0., 0.};
                                Float_t eastSet1[] = {0., 0., 0., 0., 0., 0., 0., 0., 
                                                                          0., 0., 0., 0., 0., 0., 0., 0.};
                                Float_t westSet1[] = {0., 0., 0., 0., 0., 0., 0., 0., 
                                                                          0., 0., 0., 0., 0., 0., 0., 0.};
                                Float_t eastSet2[] = {0., 0., 0., 0., 0., 0., 0., 0., 
                                                                          0., 0., 0., 0., 0., 0., 0., 0.};
                                Float_t westSet2[] = {0., 0., 0., 0., 0., 0., 0., 0., 
                                                                          0., 0., 0., 0., 0., 0., 0., 0.};
        
                                FillHitPixels(planeCdhit, 0, eastSet0, westSet0);
                                FillHitPixels(planeCdhit, 1, eastSet1, westSet1);
                                FillHitPixels(planeCdhit, 2, eastSet2, westSet2);
                                Int_t pixel = -1;
                                Float_t crossTalk = 0.; 

                                numDigits = planeCdhit.SizeSelect();

                                //loop over the digits in this plane to identify low signal and xtalk
                                //digits.  also set the DeMuxVetoFlag for all digits in this plane - if the
                                //digits have a strip end in the demuxed hypothesis, the DeMuxVetoFlag will
                                //be unset when the digit is demuxed.
                                while( (digit = planeCdhit()) ){

                                        pixel =  GetDigitPixel(digit->GetChannelId().GetVaChannel());
                                        crossTalk = 0.;
                
                                        currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(1);
          
                                        if(digit->GetCharge()
                                           <acd.GetDouble("AveragePEGainConversion")*acd.GetDouble("IgnoreSignalLimit")){
                                                digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kLowSignal);
                                                //MSG("DmxUtil", Msg::kDebug) << "digit with charge " << digit->GetCharge() 
                                                //                     << " has too little charge" << endl;
                                                --numDigits;
                                        }
            
                                        if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast 
                                                && digit->GetDeMuxDigitFlagWord() == 0){

                                                if(digit->GetChannelId().GetVaChip() == 0)
                                                        crossTalk = CheckForXTalk(pixel, eastSet0);
                                                else if(digit->GetChannelId().GetVaChip() == 1)
                                                        crossTalk = CheckForXTalk(pixel, eastSet1);
                                                else if(digit->GetChannelId().GetVaChip() == 2)
                                                        crossTalk = CheckForXTalk(pixel, eastSet2);

                                                //mark it if its cross talk and subtract its charge to the total for the side
                                                if(digit->GetCharge()<(acd.GetDouble("XTalkFractionLimit")*crossTalk)){
                                                        digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kXTalk);
                                                        --numDigits;
                                                        //MSG("DmxUtil", Msg::kDebug) << "digit with charge = " << digit->GetCharge() 
                                                        //                << " is xtalk" << endl;
                                                }
                                        }//end if digit is from the east side and not already flagged
                                        else if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest &&
                                                         digit->GetDeMuxDigitFlagWord() == 0){

                                                if(digit->GetChannelId().GetVaChip() == 0)
                                                        crossTalk = CheckForXTalk(pixel, westSet0);
                                                else if(digit->GetChannelId().GetVaChip() == 1)
                                                        crossTalk = CheckForXTalk(pixel, westSet1);
                                                else if(digit->GetChannelId().GetVaChip() == 2)
                                                        crossTalk = CheckForXTalk(pixel, westSet2);
          
                                                 //mark it if its not cross talk and add its charge to the total for the side
                                                if(digit->GetCharge()<(acd.GetDouble("XTalkFractionLimit")*crossTalk)){
                                                        digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kXTalk);
                                                        --numDigits;
                                                        //MSG("DmxUtil", Msg::kDebug) << "digit with charge = " << digit->GetCharge() 
                                                        //                << " is xtalk" << endl;
                                                }
                                        }//end if digit is from the west side and not already marked
          
                                }//end loop over digits to mark xtalk digits
                        }
                        //if not marking xtalk digits, 
                        //then set numDigits to size of the whole set
                        else numDigits = planeCdhit.SizeSelect(); 

                        planeCdhit.ResetFirst();
      
                        //determine the type of plane object by the number of digits it
                        //has.  if > 2, its a shower plane, if  <= 2 its 
                        //a muon plane.  dont do a thing if the plane is a veto shield plane
     
                        if(numDigits > 2 && planeNum>-1){
                                planeArray->AddLast(new DmxShowerPlane(acd, planeCdhit, planeNum));
                                //MSG("DmxUtil", Msg::kDebug) << "Plane " << planeNum << " is a Shower Plane" 
                                //      << " with " << numDigits << " digits" << endl;
                        }
                        else if(numDigits<=2 && planeNum>-1){ 
                                planeArray->AddLast(new DmxMuonPlane(acd, planeCdhit, planeNum));
                                //MSG("DmxUtil", Msg::kDebug) << "Plane " << planeNum << " is a Muon Plane" 
                                //      << " with " << numDigits << " digits" << endl;
                        }

                        //reset the iterator over each plane's digits
                        planeCdhit.GetSet()->ClearSlice(false);
                        
                }//end if on a new plane

        }//end loop over all digits

        //cout << planeArray->GetLast() << endl;
        status->SetPlaneArray(planeArray);

        //zero the pointer to the planeArray because the DmxStatus object now owns the 
        //array
        planeArray = 0;

        //MSG("DmxUtil", Msg::kDebug) << "done filling plane array" << endl;
        return;
}

//-----------------------------------------------------------------------
DmxUtilities::~DmxUtilities()
{
  //MSG("Dmx", Msg::kDebug) << "deleting DmxUtilities object " << endl;
}

//-----------------------------------------------------------------------
Float_t DmxUtilities::CheckForXTalk(Int_t pixel, Float_t *planeSet)
{
  
  //go through and look to see if the digit is from cross talk
  Float_t crossTalk = -1.;

  if(pixel == 0){
    crossTalk = planeSet[1] + planeSet[4];
  }
  if(pixel == 1) {
    crossTalk = planeSet[0] + planeSet[2] + planeSet[5];
  }
  if(pixel == 2){
    crossTalk = planeSet[1] + planeSet[3] + planeSet[6];
  }
  if(pixel == 3){
    crossTalk = planeSet[2] + planeSet[7];
  }
  if(pixel == 4){
    crossTalk = planeSet[0] + planeSet[5] + planeSet[8];
  }
  if(pixel == 5){
    crossTalk = planeSet[1] + planeSet[4] + planeSet[6] + planeSet[9];
  }
  if(pixel == 6){
    crossTalk = planeSet[2] + planeSet[5] + planeSet[7] + planeSet[10];
  }
  if(pixel == 7){
    crossTalk = planeSet[3] + planeSet[6] + planeSet[11];
  }
  if(pixel == 8){
    crossTalk = planeSet[4] + planeSet[9] + planeSet[12];
  }
  if(pixel == 9){
    crossTalk = planeSet[5] + planeSet[8] + planeSet[10] + planeSet[13];
  }
  if(pixel == 10){
    crossTalk = planeSet[6] + planeSet[9] + planeSet[11] + planeSet[14];
  }
  if(pixel == 11){
    crossTalk = planeSet[7] + planeSet[10] + planeSet[15];
    //MSG("DmxUtil", Msg::kDebug)<< fUtilities.GetEventNumber() << "\t" 
    //               << pixel << "\t" << "\t" << planeSet[11] << "\t" 
    //               <<planeSet[7]/planeSet[11] << "\t" 
    //               << planeSet[10]/planeSet[11] << "\t" 
    //               << planeSet[15]/planeSet[11] << endl; 
  }
  if(pixel == 12){
    crossTalk = planeSet[8] + planeSet[13];
  }
  if(pixel == 13){
    crossTalk = planeSet[9] + planeSet[12] + planeSet[14];
  }
  if(pixel == 14){
    crossTalk = planeSet[10] + planeSet[13] + planeSet[15];
    //MSG("DmxUtil", Msg::kDebug)<< fUtilities.GetEventNumber() << "\t" 
    //               << pixel << "\t" << "\t" << planeSet[14] << "\t" 
    //               <<planeSet[13]/planeSet[14] << "\t" 
    //               << planeSet[10]/planeSet[14] << "\t" 
    //               << planeSet[15]/planeSet[14] << endl; 
  }
  if(pixel == 15){
    crossTalk = planeSet[11] + planeSet[14];
  }
  
  return crossTalk;
}

//-----------------------------------------------------------------------
Int_t DmxUtilities::GetDigitPixel(Int_t channel)
{
        
  //find the pixel that the digit came from.  make stupid arrays
  //to hold the pixel to va channel info, as for some reason it 
  //isnt available from the framework
  Int_t vaChip0[] = {-1, -1, 14, 0, 15, 1, 10, 4, 11, 
                     5, 6, 8, 7, 9, 2, 13, 3, 12};
  //Int_t vaChip1[] = {-1, -1, 22, 8, 23, 9, 18, 12, 19, 
  //         13, 14, 16, 15, 17, 10, 21, 11, 20};
  //Int_t vaChip2[] = {-1, -1, 6, 16, 7, 17, 2, 20, 3, 21, 
  //         22, 0, 23, 1, 18, 5, 19, 4};
    
  Int_t pixel = vaChip0[channel];
  //if(chip == 0){ pixel = vaChip0[channel];}
  //if(chip == 1){ pixel = vaChip1[channel];}
  //if(chip == 2){ pixel = vaChip2[channel];}
 
  //MSG("Dmx2", Msg::kDebug)<< "pixel = " << pixel << endl;
  return pixel;
}

//---------------------------------------------------------------------------
void DmxUtilities::FillHitPixels(CandDeMuxDigitHandleItr crdhi, Int_t chip,
                                 Float_t *eastSet, Float_t *westSet)
{
  Int_t pixel = -1;
  //loop over the digits and see if they have been set
  while( crdhi.IsValid() ){
            
    //get the digit and it's strip, plane, and side
    CandDeMuxDigitHandle *digit = crdhi.Ptr();
    Int_t digitChip = digit->GetChannelId().GetVaChip();
    Int_t channel = digit->GetChannelId().GetVaChannel();
    pixel = GetDigitPixel(channel);
    Float_t charge = digit->GetCharge();

    if( digitChip == chip){
      if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
        eastSet[pixel] += charge;
      }
      if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
        westSet[pixel] += charge;
      }
    }
    //MSG("Dmx2", Msg::kDebug)<< pixel << "\t" << channel<< "\t" << chip << "\t" 
    //             << planeSet[pixel] << endl;
   
    crdhi.Next();
  }
  
  crdhi.ResetFirst();
  return;
}

//method for finding the vertex of an event
//----------------------------------------------------------------------------
Int_t DmxUtilities::FindVertexPlane(DmxPlaneItr &planeItr)
{
  MSG("DmxUtil", Msg::kDebug) << "In FindVertexPlane" << endl;

  planeItr.ResetFirst();

  Int_t vertexPlaneNumber = -1;
  Int_t firstPlaneNumber = 0;
  if( planeItr.IsValid() )
    firstPlaneNumber = planeItr.Ptr()->GetPlaneNumber();

  planeItr.ResetFirst();

  //find the vertex.  make it the first plane in a set of 5 planes where at least 3
  //of the planes are valid.  use the number 3 since it is the least number of planes
  //you can have and still drop one to refine a fit.  use 5 because the trigger is 4/5
  Int_t beginPlaneNumber = 0;
  Int_t endPlaneNumber = 0;
  Int_t validCtr = 0;
  Int_t cnt = 0;

  while( planeItr.IsValid() ){ 

    beginPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
    endPlaneNumber = 0;
    validCtr = 0;
    
    cnt = 0;
    while( planeItr.IsValid() && cnt < 5 && validCtr < 3){
      if( planeItr.Ptr()->IsValid() ){ ++validCtr;}
      ++cnt;
      endPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
      MSG("DmxUtil", Msg::kDebug)<<"Start Plane = " << beginPlaneNumber 
                                << "\tThis Plane = " 
                                << endPlaneNumber << "\tValid = " 
                                << (Int_t)planeItr.Ptr()->IsValid() 
                                << "\tValidCtr = " << validCtr << endl;
      planeItr.Next();
    }
        
    if( validCtr >= 3 && (endPlaneNumber-beginPlaneNumber) <= 5 ){  
      vertexPlaneNumber = beginPlaneNumber;
      MSG("DmxUtil", Msg::kDebug)<<"\tVertex Plane = " << vertexPlaneNumber << endl;
      break;
    }
    else{
      for(Int_t a = 0; a < cnt-1; a++){
        planeItr.Prev();
      }
    }
  } // end loop to find vertex
  
  //see if you can find any contiguos planes before the vertex, if so, make them the vertex
  //slice from 0 to the vertex and then iterate backwards.  if there is a gap in signal
  //of more than 2 planes, stop
  
  planeItr.Reset();
  planeItr.GetSet()->ClearSlice(false);

  Int_t nextPlane = 0;

  if( vertexPlaneNumber > -1 && firstPlaneNumber < vertexPlaneNumber){
    planeItr.GetSet()->Slice(firstPlaneNumber, vertexPlaneNumber);
    planeItr.GetSet()->Update();

    planeItr.ResetLast();

    while( planeItr.IsValid() ){
    
      Int_t currentPlane = planeItr.Ptr()->GetPlaneNumber();
      MSG("DmxUtil", Msg::kDebug) << "\tcurrentPlaneNumber = " << currentPlane;
      planeItr.Prev();
      nextPlane = 0;
      if(planeItr.IsValid()){
        nextPlane = planeItr.Ptr()->GetPlaneNumber();
        MSG("DmxUtil", Msg::kDebug) << "\tnextPlaneNumber = " << nextPlane << endl;
      
        if( vertexPlaneNumber - nextPlane < 3)
          vertexPlaneNumber = nextPlane;
        else break;
      }
      
    }//end loop over planes

    planeItr.GetSet()->ClearSlice(false);
    planeItr.ResetFirst();
  }
  
  planeItr.GetSet()->ClearSlice(false);
  planeItr.ResetFirst();
  
  MSG("DmxUtil", Msg::kDebug) << "\tvertexPlaneNumber = " << vertexPlaneNumber << endl;
  return vertexPlaneNumber;
}

//--------------------------------------------------------------------------------
Int_t DmxUtilities::FindEndPlane(DmxPlaneItr &planeItr)
{
  Int_t endPlaneNumber = 500;

  while( planeItr.IsValid() && endPlaneNumber == 500){
    Int_t planeNumber = planeItr.Ptr()->GetPlaneNumber();
    
    planeItr.Next();
    if( planeItr.IsValid() ){
      //if the difference between the current plane number and the previous is >= 10
      //call the previous plane the last one in the event
      if( TMath::Abs(planeItr.Ptr()->GetPlaneNumber() 
                     - planeNumber) >= 5 ){ endPlaneNumber = planeNumber; }
    }
    else{
      //didn't find a cutoff plane so call the last plane it
      endPlaneNumber = planeNumber;
    }
  }

  return endPlaneNumber;
}

//-----------------------------------------------------------------------------
Int_t DmxUtilities::FindMuonStartPlane(DmxPlaneItr &planeItr, Float_t peGainConversion)
{
        //now locate where the muon starts, if anywhere - define the start
        //plane to be the first muon plane in a set of 4/5 muon planes with 
        //3 of the 4 having at least 3 pe's each
  
        //this method assumes that the planeItr was sliced to be from the vertex to the 
        //end planes of the event

        Int_t fMuonStartPlaneNumber = -1;

        Int_t ctr = 0;
        Int_t muonCtr = 0;
        Int_t startPlane = 0;
        Int_t endPlane = 0;
        Int_t peCtr = 0;
        
        DmxPlane *currentPlane = 0;
        
        while( (currentPlane = planeItr()) ){
                
                MSG("DmxUtilities", Msg::kDebug) << "plane " << currentPlane->GetPlaneNumber()
                                                                                << " type " << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
                                                                                << " charge " << currentPlane->GetPlaneCharge()
                                                                                << " strips " << currentPlane->GetNumberOfStrips()
                                                                                << " is valid " << (Int_t)currentPlane->IsValid() << endl;
        }

        planeItr.ResetFirst();
        
        while(planeItr.IsValid() && fMuonStartPlaneNumber == -1){
                ctr = 0;

                muonCtr = 0;
                startPlane = planeItr.Ptr()->GetPlaneNumber();
                endPlane = startPlane+10;
                peCtr = 0;
    
                MSG("DmxUtilities", Msg::kDebug) << "plane number " << startPlane 
                                                                                << " charge " << planeItr.Ptr()->GetPlaneCharge()
                                                                                << " type " 
                                                                            << DmxPlaneTypes::GetPlaneTypeString(planeItr.Ptr()->GetPlaneType()) 
                                                                                << endl;
                planeItr.Next();
                        
                //check 4 planes at a time
                while(ctr < 4 && planeItr.IsValid() ){
                        currentPlane = planeItr.Ptr();

                        MSG("DmxUtilities", Msg::kDebug) << "plane number " << currentPlane->GetPlaneNumber()
                                                                                        << " charge " << currentPlane->GetPlaneCharge()
                                                                                        << " type " 
                                                                                        << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType()) 
                                                                                        << " num strips " << currentPlane->GetNumberOfStrips()
                                                                                        << endl;
                        endPlane = currentPlane->GetPlaneNumber();

                        //check to see that its a muon plane or only has 2 strips hit
                        if( currentPlane->GetNumberOfStrips() < 3 )     ++muonCtr;
                        
                        if(currentPlane->GetPlaneCharge() > (peGainConversion * 3.)) ++peCtr;
                        
                        MSG("DmxUtilities", Msg::kDebug) << " muon planes = " << muonCtr 
                                                                                        << " peCtr = " << peCtr << endl;
                        planeItr.Next();
                        ++ctr;
                }
    
                if(muonCtr == 4 && peCtr >= 3 && endPlane-startPlane<=5) 
                        fMuonStartPlaneNumber = startPlane;
    
                //if there are still planes left to examine, the iterator will be valid
                //otherwise, you ran out of planes, so break the cycle
                if(planeItr.IsValid() ){
                        planeItr.Prev();
                        planeItr.Prev();
                        planeItr.Prev();
                        planeItr.Prev();
                }
        }//end loop to find the muon track starting plane
  
        planeItr.GetSet()->ClearSlice(false);
        planeItr.Reset();
  
        return fMuonStartPlaneNumber;
}

//-----------------------------------------------------------------------------
Int_t DmxUtilities::FindBeamMuonStartPlane(DmxPlaneItr &planeItr, Float_t peGainConversion)
{
        //now locate where the muon starts, if anywhere - define the start
        //plane to be the first muon plane in a set of 4/5 muon planes with 
        //3 of the 4 having at least 3 pe's each
        
        //this method assumes that the planeItr was sliced to be from the vertex to the 
        //end planes of the event
        
        Int_t fMuonStartPlaneNumber = -1;
        
        Int_t ctr = 0;
        Int_t muonCtr = 0;
        Int_t startPlane = 0;
        Int_t endPlane = 0;
        Int_t peCtr = 0;
        Int_t lastEventPlane = 0;
        
        DmxPlane *currentPlane = 0;
        
        while( (currentPlane = planeItr()) ){
                
                MSG("DmxUtilities", Msg::kDebug) << "plane " << currentPlane->GetPlaneNumber()
                << " type " << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
                << " charge " << currentPlane->GetPlaneCharge()
                << " strips " << currentPlane->GetNumberOfStrips()
                << " is valid " << (Int_t)currentPlane->IsValid() << endl;
                
                lastEventPlane = currentPlane->GetPlaneNumber();
        }
        
        planeItr.ResetFirst();
        
        while(planeItr.IsValid() && fMuonStartPlaneNumber == -1){
                ctr = 0;
                
                muonCtr = 0;
                startPlane = planeItr.Ptr()->GetPlaneNumber();
                endPlane = startPlane+10;
                peCtr = 0;
                
                MSG("DmxUtilities", Msg::kDebug) << "plane number " << startPlane 
                        << " charge " << planeItr.Ptr()->GetPlaneCharge()
                        << " type " 
                        << DmxPlaneTypes::GetPlaneTypeString(planeItr.Ptr()->GetPlaneType()) 
                        << endl;
                planeItr.Next();
                
                //check 4 planes at a time
                while(ctr < 4 && planeItr.IsValid() ){
                        currentPlane = planeItr.Ptr();
                        
                        MSG("DmxUtilities", Msg::kDebug) << "plane number " << currentPlane->GetPlaneNumber()
                                << " charge " << currentPlane->GetPlaneCharge()
                                << " type " 
                                << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType()) 
                                << " num strips " << currentPlane->GetNumberOfStrips()
                                << endl;
                        endPlane = currentPlane->GetPlaneNumber();
                        
                        //check to see that its a muon plane or only has 2 strips hit
                        if( currentPlane->GetNumberOfStrips() < 3 )     ++muonCtr;
                        
                        if(currentPlane->GetPlaneCharge() > (peGainConversion * 3.)) ++peCtr;
                        
                        MSG("DmxUtilities", Msg::kDebug) << " muon planes = " << muonCtr 
                                << " peCtr = " << peCtr << endl;
                        planeItr.Next();
                        ++ctr;
                }
                
                if(muonCtr == 4 && peCtr >= 3 && endPlane-startPlane<=5) 
                        fMuonStartPlaneNumber = startPlane;
                
                //if there are still planes left to examine, the iterator will be valid
                //otherwise, you ran out of planes, so break the cycle
                if(planeItr.IsValid() ){
                        planeItr.Prev();
                        planeItr.Prev();
                        planeItr.Prev();
                        planeItr.Prev();
                }
        }//end loop to find the muon track starting plane
        
        //found the most likely start plane for the muon.  but lets make sure 
        //we arent being fooled into thinking a shower event is really a muon event.  
        //loop over the planes and see the number of planes that are shower like and if 
        //there are lots more shower planes than muon planes
        planeItr.GetSet()->ClearSlice(false);
        
        planeItr.GetSet()->Slice(fMuonStartPlaneNumber, lastEventPlane);
        planeItr.GetSet()->Update();
        planeItr.ResetFirst();
        
        Int_t numShowerPlanes = 0;
        Int_t numMuonPlanes = 0;
        while( (currentPlane = planeItr()) ){
                
                if(currentPlane->GetPlaneType() == DmxPlaneTypes::kShower
                   && currentPlane->GetNumberOfStrips()>2) ++numShowerPlanes;
                else ++numMuonPlanes;
                
        }
        
        if(1.*numShowerPlanes>0.25*numMuonPlanes) fMuonStartPlaneNumber = -1;

        planeItr.GetSet()->ClearSlice(false);
        
        return fMuonStartPlaneNumber;
}

//---------------------------------------------------------------------------
//function to see if a fit makes since physically - ie its in the detector!
//also check to see if the fit has at least 50% of the signal mated
Bool_t DmxUtilities::IsValidFit(DmxPlaneItr &planeItr, Double_t a1, Double_t a2, 
                                                                Double_t a3, Double_t a4, Float_t offset)
{
        Int_t nonPhysicalCtr = 0;
        bool validFit = true;
        Int_t planesInEvent = planeItr.SizeSelect();
        //MSG("DmxUtil", Msg::kDebug)<<"\tin IsValidFit " << a << "\t" << b << endl;
        planeItr.Reset();

        DmxPlane *plane = 0;
        Double_t fitCoG = 0.;
        
        while( planeItr.IsValid() ){
   
                plane = planeItr.Ptr();
    
                fitCoG = (a1 + a2*(plane->GetZPosition()-offset) 
                                  + a3*TMath::Power(plane->GetZPosition()-offset,2)
                                  + a4*TMath::Power(plane->GetZPosition()-offset,3));
      //MSG("DmxUtil", Msg::kDebug)<<"\t" << plane->GetZPosition() << endl;
      
      //if the fit CoG for the plane is outside of the detector, it is non Physical
      if( TMath::Abs(fitCoG) >= 4.5) ++nonPhysicalCtr;
  
    planeItr.Next(); 
  }

  planeItr.Reset();
  
  if(planesInEvent<8 && (1.*nonPhysicalCtr)/(1.*planesInEvent)>0.33) validFit = false;
  else if(planesInEvent>=8 && nonPhysicalCtr>=4) validFit = false;
  
  return validFit;
}
//---------------------------------------------------------------------------
//function to see if a fit makes since physically - ie its in the detector!
//also check to see if the fit has at least 50% of the signal mated
Bool_t DmxUtilities::IsValidFit(DmxPlaneItr &planeItr, Double_t a1, Double_t a2, 
                                                                Double_t a3, Double_t a4)
{
        Int_t nonPhysicalCtr = 0;
        bool validFit = true;
        Int_t planesInEvent = planeItr.SizeSelect();
        //MSG("DmxUtil", Msg::kDebug)<<"\tin IsValidFit " << a << "\t" << b << endl;
        planeItr.Reset();

        DmxPlane *plane = 0;
        Double_t fitCoG = 0.;

        while( planeItr.IsValid() ){
   
                plane = planeItr.Ptr();
    
                fitCoG = (a1 + a2*(plane->GetZPosition()) 
                                  + a3*TMath::Power(plane->GetZPosition(),2)
                                  + a4*TMath::Power(plane->GetZPosition(),3));
                //MSG("DmxUtil", Msg::kDebug)<<"\t" << plane->GetZPosition() << endl;
      
                //if the fit CoG for the plane is outside of the detector, it is non Physical
                if( TMath::Abs(fitCoG) >= 4.5) ++nonPhysicalCtr;
  
                planeItr.Next(); 
        }

        planeItr.Reset();
  
        if(planesInEvent<8 && (1.*nonPhysicalCtr)/(1.*planesInEvent)>0.33) validFit = false;
        else if(planesInEvent>=8 && nonPhysicalCtr>=4) validFit = false;
  
        return validFit;
}

//---------------------------------------------------------------------------
//function to fix the cases where an end point of the fit doesnt jive
//with the best hyps in a plane
Int_t DmxUtilities::CheckFit(DmxPlaneItr &planeItr)
{
  DmxPlane *plane = 0;

  //compare the best hyps in valid shower planes to the fit hyp for the amount of signal
  //mated in the plane.  pick the one with the most mated signal, as long as the difference
  //is significant - like 0.2 or more.  also only use the best hyps if they have more than
  //85% of the signal mated in them
  Float_t firstBestMated = 0.;
  Float_t firstSMated = 0.;
  Float_t firstTMated = 0.;
  Float_t firstMated = 0.;
  Float_t lastBestMated = 0.;
  Float_t lastSMated = 0.;
  Float_t lastTMated = 0.;
  Float_t lastMated = 0.;
  
  Int_t ctr = 0;
  Int_t totPlanes = planeItr.SizeSelect();
  Int_t fixedPlanes = 0;
  planeItr.ResetFirst();

  //make an array to hold the CoG positions of the planes 
  Float_t firstPlaneCoG = -1.;
  Float_t firstPlaneBCoG = -1.;
  Float_t firstPlaneSCoG = -1.;
  Float_t firstPlaneTCoG = -1.;
  Float_t lastPlaneCoG = -1.;
  Float_t lastPlaneBCoG = -1.;
  Float_t lastPlaneSCoG = -1.;
  Float_t lastPlaneTCoG = -1.;
  Float_t firstPlaneSpace = 0.;
  Float_t lastPlaneSpace = 0.;

  while(planeItr.IsValid()){
    
    plane = planeItr.Ptr();
    if(ctr == 0){
      firstPlaneCoG  = plane->GetCoG();
      firstPlaneSpace = plane->GetPlaneNumber();
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
   
        firstPlaneBCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
        firstBestMated = dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
      
        firstPlaneSCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
        firstSMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
        
        firstPlaneTCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
        firstTMated = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
      
        firstMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
      }
      else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
        firstPlaneBCoG = firstPlaneCoG;
        firstPlaneSCoG = firstPlaneCoG;
        firstPlaneTCoG = firstPlaneCoG;
        firstMated = 1.;
        firstBestMated = 1.;
        firstSMated = 1.;
        firstTMated = 1.;
      }
    }
    else if(ctr == 1) firstPlaneSpace -= plane->GetPlaneNumber();
    else if(ctr == totPlanes-2) lastPlaneSpace = plane->GetPlaneNumber();
    else if(ctr == totPlanes-1){
      lastPlaneCoG  = plane->GetCoG();
      lastPlaneSpace -= plane->GetPlaneNumber();
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
   
        lastPlaneBCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
        lastBestMated = dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
      
        lastPlaneSCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
        lastSMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
        
        lastPlaneTCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
        lastTMated = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
      
        lastMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
      }
      else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
        lastPlaneBCoG = lastPlaneCoG;
        lastPlaneSCoG = lastPlaneCoG;
        lastPlaneTCoG = lastPlaneCoG;
        lastMated = 1.;
        lastBestMated = 1.;
        lastSMated = 1.;
        lastTMated = 1.;
      }
    }
    ++ctr;
    planeItr.Next();
  }
  planeItr.ResetFirst();

  //which is the better hypothesis among the 3 best?
  bool useFirstB = false;
  bool useFirstS = false;
  bool useFirstT = false;
  bool useLastB = false;
  bool useLastS = false;
  bool useLastT = false;

  if(firstBestMated>=firstSMated && firstBestMated>=firstTMated 
     && TMath::Abs(firstBestMated-firstMated)>0.2) useFirstB = true;
  else if(firstSMated>firstBestMated && firstSMated>=firstTMated 
     && TMath::Abs(firstSMated-firstMated)>0.2) useFirstS = true;
  else if(firstTMated>firstSMated && firstTMated>firstBestMated 
     && TMath::Abs(firstTMated-firstMated)>0.2) useFirstT = true;
 
  if(lastBestMated>=lastSMated && lastBestMated>=lastTMated 
     && TMath::Abs(lastBestMated-lastMated)>0.2) useLastB = true;
  else if(lastSMated>lastBestMated && lastSMated>=lastTMated 
     && TMath::Abs(lastSMated-lastMated)>0.2) useLastS = true;
  else if(lastTMated>lastSMated && lastTMated>lastBestMated 
     && TMath::Abs(lastTMated-lastMated)>0.2) useLastT = true;


//    cout << "first plane: " << firstPlaneCoG << " " 
//      << firstPlaneBCoG << " " << firstPlaneSCoG 
//         << " " << firstPlaneTCoG << " " << firstMated 
//      << " " << firstBestMated << " " 
//         << firstSMated << " " << firstTMated << " " 
//      << useFirstB << " " << useFirstS << " " 
//         << useFirstT << " " << firstPlaneSpace << endl;

//    cout << "last plane: " << lastPlaneCoG << " " << lastPlaneBCoG 
//      << " " << lastPlaneSCoG << " " << lastPlaneTCoG << " " << lastMated 
//      << " " << lastBestMated << " " << lastSMated << " " << lastTMated 
//      << " " << useLastB << " " << useLastS << " " << useLastT 
//      << " " << lastPlaneSpace << endl;

  //the first and last planes in the view should be reset to one of the
  //best hypothesis
  //
  //1.  the difference between the best hypothesis and the nearest neighboring
  //    valid plane hypothesis is less than 0.25m/plane
  //
  //    and
  //
  //2.  the best hypothesis has a lot more signal mated

  //plane iterator was previously set to the first plane in the set
    
  if(useFirstB 
     && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneBCoG)){
    plane = planeItr.Ptr();
    if(plane->IsStray()){
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
        dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("best");
      else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon) plane->SetStrips();
      plane->SetStray(false);
      ++fixedPlanes;
    }
  }
  else if(useFirstS 
          && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneSCoG)){
    plane = planeItr.Ptr();
    if(plane->IsStray()){
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
        dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("second");
      plane->SetStray(false);
      ++fixedPlanes;
    }
  }
  else if(useFirstT 
          && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneTCoG)){
    plane = planeItr.Ptr();
    if(plane->IsStray()){
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
        dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("third");
      plane->SetStray(false);
      ++fixedPlanes;
    }
  }

  //set the iterator to the last plane
  planeItr.ResetLast();

  if(useLastB && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneBCoG)){
    plane = planeItr.Ptr();
    if(plane->IsStray()){
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
        dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("best");
      else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon) plane->SetStrips();
      plane->SetStray(false);
      ++fixedPlanes;
    }
  }
  else if(useLastS 
          && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneSCoG)){
    plane = planeItr.Ptr();
    if(plane->IsStray()){
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
        dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("second");
      plane->SetStray(false);
      ++fixedPlanes;
    }
  }
  else if(useLastT 
          && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneTCoG)){
    plane = planeItr.Ptr();
    if(plane->IsStray()){
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
        dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("third");
      plane->SetStray(false);
      ++fixedPlanes;
    }
  }
  
  return fixedPlanes;
}

//---------------------------------------------------------------------------
//method to check the demuxing solution
Bool_t DmxUtilities::IsOverlappingMultiple(DmxPlaneItr &planeItr, Float_t vertexZ, 
                                              Float_t slopeCutOff, Float_t interceptCutOff,
                                              Float_t peakCutOff)
{

  planeItr.ResetFirst();
  DmxPlane *plane = 0;

  //make a histogram for each supermodule.  this is because the method is trickier in 
  //tpos-z space than strip-plane space, so just deal with the gap by pretending its not
  //there
  TH2F houghSpaceSM1 = TH2F("houghSpaceSM1", "Hough Space for Multiples", 
                            30, -30, 30, 500, -500, 500);
  TH2F houghSpaceSM2 = TH2F("houghSpaceSM2", "Hough Space for Multiples", 
                            30, -30, 30, 500, -500, 500);
  TH2F houghSpaceCutSM1 = TH2F("hS75SM1", "Hough Space for U with cut", 
                               30, -30, 30, 500, -500, 500);
  TH2F houghSpaceCutSM2 = TH2F("hS75SM2", "Hough Space for U with cut", 
                               30, -30, 30, 500, -500, 500);
      
  Float_t planeZ = 0.;
  Float_t cog = 0.;
  while( (plane = planeItr()) ){
    
    planeZ = 1.*plane->GetPlaneNumber();
    
    if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
     
      for(Int_t j = 0; j <3; j++){
        if( j==0 && dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis() ){
          cog = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
        }
        else if( j==1 && dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis() ){
          cog = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
        }
        else if( j==2 && dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis() ){
          cog = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
        }
        //loop over possible slopes
        for(Int_t m = -29; m<30; m+=2){
          
          Float_t c = cog - m*1.*(planeZ - vertexZ);
          
          if(planeZ<249) houghSpaceSM1.Fill(m,c);
          else if(planeZ>249) houghSpaceSM2.Fill(m,c);
        
        }//end loop over slopes
      }//end loop over hypotheses
      
    }//end if shower plane
    else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
      cog = dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG();
      
      //loop over possible slopes
      for(Int_t m = -29; m<30; m+=2){
        
        Float_t c = cog - m*1.*(planeZ - vertexZ);

        //fill the space with this entry 3 times so it is weighted the same as a shower plane
        if(planeZ<249){
          houghSpaceSM1.Fill(m,c);
          houghSpaceSM1.Fill(m,c);
          houghSpaceSM1.Fill(m,c);
        }
        else if(planeZ>249){
          houghSpaceSM2.Fill(m,c);
          houghSpaceSM2.Fill(m,c);
          houghSpaceSM2.Fill(m,c);
        }
      }//end loop over slopes
            
    }//end if muon plane
  }//end loop over planes

  planeItr.ResetFirst();

  //now fill a second histogram with the bins that are at least peakCutOff of the maximum
  //bin in the previous histogram
  Stat_t binContent = 0;
  Double_t maximumSM1 = houghSpaceSM1.GetMaximum();     
  Double_t maximumSM2 = houghSpaceSM2.GetMaximum();     
        
  for(Int_t x = 1; x < 30+1; x++){
    //cout << x << endl;
    for(Int_t y = 1; y < 500+1; y++){
      binContent = houghSpaceSM1.GetBinContent(x,y);
      if(binContent>= peakCutOff*maximumSM1) houghSpaceCutSM1.SetBinContent(x,y,binContent);
      binContent = houghSpaceSM2.GetBinContent(x,y);
      if(binContent>= peakCutOff*maximumSM2) houghSpaceCutSM2.SetBinContent(x,y,binContent);
          
    }
  }

  //if the rms in the slope is less than the value determined by tests
  //and the rms in the intercept is greater than the cutoff, mark the event
  //as a possible overlapping multiple
  bool multiple = false;
  if(houghSpaceCutSM1.GetRMS(1)<=slopeCutOff 
     && houghSpaceCutSM1.GetRMS(2)>=interceptCutOff)
    multiple = true;
  else if(houghSpaceCutSM2.GetRMS(1)<=slopeCutOff 
          && houghSpaceCutSM2.GetRMS(2)>=interceptCutOff)
    multiple = true;

  return multiple;
}

//-------------------------------------------------------------------
//method to help determine the goodness of the demuxing solution
Int_t DmxUtilities::FindPlanesOffFit(DmxPlaneItr &planeItr, Int_t strayCut)
{
  
  planeItr.Reset();

  Int_t strayPlanes = 0;

  DmxPlane *plane = 0;

  //MSG("DmxUtil", Msg::kDebug) << "checking for stray planes in event " << endl;
  while( planeItr.IsValid() ){
    
    plane = planeItr.Ptr();
    
    if( plane->IsValid() ){
      //&& plane->GetPlaneNumber()>= firstPlane && plane->GetPlaneNumber()<= lastPlane){
      Double_t fitCoG = -1.;
      
      
      if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
        Double_t diff1 = 0.;
        Double_t diff2 = 0.;
        Double_t diff3 = 0.;
        fitCoG = plane->GetCoG();

        if( dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis() ){
          diff1 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG());
          //MSG("DmxUtil", Msg::kDebug)<<"best " << dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG()<<endl;
        }
        if( dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis() ){
          diff2 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG());
          //MSG("DmxUtil", Msg::kDebug)<<"2nd " <<dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG()<<endl;
        }
        if( dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis() ){
          diff3 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG());
          //MSG("DmxUtil", Msg::kDebug)<<"3rd " << dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG()<<endl;
        }
        
        //fStrayCut is in strip numbers, the 0.042 factor is the width of a strip
        if(diff1>strayCut*.042 && diff2>strayCut*.042 && diff3>strayCut*.042){ 
          //MSG("DmxUtil", Msg::kDebug) << "plane = " << plane->GetPlaneNumber() << " is stray "
          //              << diff1 << " " << diff2 << " " << diff3 << endl;
      
          ++strayPlanes; 
          plane->SetStray(true);
          // MSG("DmxUtil", Msg::kDebug) << "\tstray plane = " << plane->GetPlaneNumber()  
          //              << "\tfit cog = " << fitCoG << endl;
        }
      }//end if shower plane
     //  else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
//      fitCoG = plane->GetCoG();
//      if( TMath::Abs(fitCoG-dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()) > strayCut*0.042){ 
//        //MSG("DmxUtil", Msg::kDebug)<<"muon " << plane->GetCoG()<<endl;
//        ++strayPlanes; 
//        plane->SetStray(true);
//        MSG("DmxUtil", Msg::kDebug) << "\tstray plane = " << plane->GetPlaneNumber() 
//                                << "\tfit cog = " << fitCoG << " initial cog = " 
//                                << dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()
//                                <<endl;
//      }
//       }//end if muon plane
    }//end if it is a valid plane
    planeItr.Next();
  }//end loop over planes

  planeItr.Reset();

  return strayPlanes;
}

//--------------------------------------------------------------------------
Float_t DmxUtilities::CheckFitWithTiming(DmxPlaneItr &planeItr)
{
  
  //the idea here is to get the offset from the center line from timing for a plane
  //in one view and compare that to the center of gravity from the demuxer of the neighboring
  //plane in the other view.

  //take the average of the centers of gravity for the planes on either side of the plane you 
  //are interested in 

  planeItr.ResetFirst();

  Int_t planeCtr = 0;
  Int_t numPlanes = planeItr.SizeSelect();
  Int_t forwardCtr = 0;
  Int_t backwardCtr = 0;
  bool opposite = false;
  Float_t offset = 0.;
  Float_t avOffsetDiff = 0.;
  Float_t avCoG = 0.;
  Float_t avCoGCnt = 0.;

  DmxPlane *plane = 0;
  DmxPlane *nextPlane = 0;
  DmxPlane *prevPlane = 0;

  while( planeItr.IsValid() ){

    plane = planeItr.Ptr();
    //MSG("DmxUtil", Msg::kDebug) << "plane of offset = " << plane->GetPlaneNumber() << endl; 

    offset = plane->GetTimingOffset();
    //MSG("DmxUtil", Msg::kDebug) << " offset = " << offset << " is valid = " << (Int_t)plane->IsValid() << endl;   
    forwardCtr = 0;
    backwardCtr = 0;
    avCoG = 0.;
    avCoGCnt = 0.;
    //this is the first plane, so you can only look at the one after it
    if(planeCtr == 0 && offset != -10.){

      //advance to the next plane in the opposite view
      while(planeItr.IsValid() && opposite == false){
        planeItr.Next();
        ++forwardCtr;
        nextPlane = planeItr.Ptr();
        if(nextPlane->GetPlaneView() != plane->GetPlaneView()){
          //MSG("DmxUtil", Msg::kDebug) << "next opposite plane to vertex = " << nextPlane->GetPlaneNumber() 
          //                 << " cog = " << nextPlane->GetCoG() << endl;
          opposite = true;
          avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(nextPlane->GetCoG()));
        }
      }//end loop to find next plane in opposite view

      while(forwardCtr > 0){
        planeItr.Prev();
        --forwardCtr;
      }//end loop to return planeItr to original plane
    
    }//end if first plane in set
    else if(planeCtr > 0 && planeCtr<numPlanes-1 && offset != -10.){

      opposite = false;
      //advance to the next plane in the opposite view
      while(planeItr.IsValid() && !opposite){
        planeItr.Next();
        ++forwardCtr;
        if(planeItr.IsValid()){
          nextPlane = planeItr.Ptr();
          if(nextPlane->GetPlaneView() != plane->GetPlaneView()){
            //MSG("DmxUtil", Msg::kDebug) << "next opposite plane = " << nextPlane->GetPlaneNumber() 
            //                 << " cog = " << nextPlane->GetCoG() << endl;
            opposite = true;
            avCoG += nextPlane->GetCoG();
            avCoGCnt += 1.;
          }
        }
      }//end loop to find next plane in opposite view

      while(forwardCtr > 0){
        planeItr.Prev();
        --forwardCtr;
      }//end loop to return planeItr to original plane
    
      opposite = false;
      
      //if(planeItr.IsValid() )MSG("DmxUtil", Msg::kDebug) << "start backward look from plane " 
      //                                                << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber()
      //                                                << endl;
      while(planeItr.IsValid() && !opposite){
        planeItr.Prev();
        ++backwardCtr;
        if(planeItr.IsValid()){
          prevPlane = planeItr.Ptr();
          if(prevPlane->GetPlaneView() != plane->GetPlaneView()){
            //MSG("DmxUtil", Msg::kDebug) << "prev opposite plane = " << prevPlane->GetPlaneNumber() 
            //                 << " cog = " << prevPlane->GetCoG() << endl;
            opposite = true;
            avCoG += prevPlane->GetCoG();
            avCoGCnt += 1.;
          }
        }
      }//end loop to find prev plane in opposite view

      while(backwardCtr > 0){
        planeItr.Next();
        --backwardCtr;
      }//end loop to return planeItr to original plane
    
      //get the average CoG for the planes on either side and subtract the timing offset
      if(avCoGCnt > 0.) avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(avCoG/avCoGCnt));

    }//end if plane between end planes
    else if(planeCtr == numPlanes-1 && offset != -10.){
      opposite = false;
      while(planeItr.IsValid() && !opposite){
        planeItr.Prev();
        ++backwardCtr;
        if(planeItr.IsValid()){
          prevPlane = planeItr.Ptr();
          if(prevPlane->GetPlaneView() != plane->GetPlaneView()){
            //MSG("DmxUtil", Msg::kDebug) << "prev opposite plane to end = " <<prevPlane->GetPlaneNumber() 
            //                 << " cog = " << prevPlane->GetCoG() << endl;
            opposite = true;
            avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(prevPlane->GetCoG()));
          }
        }
      }//end loop to find prev plane in opposite view

      while(backwardCtr > 0){
        planeItr.Next();
        --backwardCtr;
      }//end loop to return planeItr to original plane

    }//end if plane is end plane

    ++planeCtr;
    planeItr.Next();
  }//end loop over planes

  if(numPlanes>0) avOffsetDiff /= (1.0*numPlanes);

  planeItr.ResetFirst();

  return avOffsetDiff;
}

//----------------------------------------------------------------------
//method to find a least squares fit to a cubic function of the form
//y = a1 + a2*x + a3*x^2 + a4*x^3
void DmxUtilities::FindCubicFit(Double_t *x, Double_t *y, Double_t *weight,  
                                Int_t nPoints, Double_t &a1, Double_t &a2, 
                                Double_t &a3, Double_t &a4, Double_t &chiSq)
{

  //use method of determinants to do the least squares fit
  //this is from bevington chapter 7
  //
  //f_1(x) = 1
  //f_2(x) = x
  //f_3(x) = x^2
  //f_4(x) = x^3
  
  //need to find the sums
  //SIGMA y_i*f_1(x_i)/weight_i
  //SIGMA y_i*f_2(x_i)/weight_i
  //SIGMA y_i*f_3(x_i)/weight_i
  //SIGMA y_i*f_4(x_i)/weight_i
  //SIGMA f_1(x_i)*f_1(x_i)/weight_i
  //SIGMA f_1(x_i)*f_2(x_i)/weight_i
  //SIGMA f_1(x_i)*f_3(x_i)/weight_i
  //SIGMA f_1(x_i)*f_4(x_i)/weight_i
  //SIGMA f_2(x_i)*f_1(x_i)/weight_i
  //SIGMA f_2(x_i)*f_2(x_i)/weight_i
  //SIGMA f_2(x_i)*f_3(x_i)/weight_i
  //SIGMA f_2(x_i)*f_4(x_i)/weight_i
  //SIGMA f_3(x_i)*f_1(x_i)/weight_i
  //SIGMA f_3(x_i)*f_2(x_i)/weight_i
  //SIGMA f_3(x_i)*f_3(x_i)/weight_i
  //SIGMA f_3(x_i)*f_4(x_i)/weight_i
  //SIGMA f_4(x_i)*f_1(x_i)/weight_i
  //SIGMA f_4(x_i)*f_2(x_i)/weight_i
  //SIGMA f_4(x_i)*f_3(x_i)/weight_i
  //SIGMA f_4(x_i)*f_4(x_i)/weight_i
  //
  //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
  //10 instead of 16 variables

  Double_t yf1divWeight = 0.;
  Double_t yf2divWeight = 0.;
  Double_t yf3divWeight = 0.;
  Double_t yf4divWeight = 0.;
  Double_t f1f1divWeight = 0.;
  Double_t f1f2divWeight = 0.;
  Double_t f1f3divWeight = 0.;
  Double_t f1f4divWeight = 0.;
  Double_t f2f2divWeight = 0.;
  Double_t f2f3divWeight = 0.;
  Double_t f2f4divWeight = 0.;
  Double_t f3f3divWeight = 0.;
  Double_t f3f4divWeight = 0.;
  Double_t f4f4divWeight = 0.;
  
  for(Int_t i=0; i<nPoints; i++){
    yf1divWeight += y[i]/(weight[i]*weight[i]);
    yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
    yf3divWeight += (y[i]*x[i]*x[i])/(weight[i]*weight[i]);
    yf4divWeight += (y[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f1f1divWeight += 1./(weight[i]*weight[i]);
    f1f2divWeight += x[i]/(weight[i]*weight[i]);
    f1f3divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
    f1f4divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
    f2f3divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f2f4divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f3f3divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f3f4divWeight += (x[i]*x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f4f4divWeight += (x[i]*x[i]*x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
  }

  //make a bunch of matricies to hold these values
  TMatrix deltaMatrix(4,4);
  TMatrix a1Matrix(4,4);
  TMatrix a2Matrix(4,4);
  TMatrix a3Matrix(4,4);
  TMatrix a4Matrix(4,4);
  
  //fill the matricies
  deltaMatrix(0,0) = f1f1divWeight;
  deltaMatrix(0,1) = f1f2divWeight;
  deltaMatrix(0,2) = f1f3divWeight;
  deltaMatrix(0,3) = f1f4divWeight;
  deltaMatrix(1,0) = f1f2divWeight;
  deltaMatrix(1,1) = f2f2divWeight;
  deltaMatrix(1,2) = f2f3divWeight;
  deltaMatrix(1,3) = f2f4divWeight;
  deltaMatrix(2,0) = f1f3divWeight;
  deltaMatrix(2,1) = f2f3divWeight;
  deltaMatrix(2,2) = f3f3divWeight;
  deltaMatrix(2,3) = f3f4divWeight;
  deltaMatrix(3,0) = f1f4divWeight;
  deltaMatrix(3,1) = f2f4divWeight;
  deltaMatrix(3,2) = f3f4divWeight;
  deltaMatrix(3,3) = f4f4divWeight;
  
  a1Matrix(0,0) = yf1divWeight;
  a1Matrix(0,1) = f1f2divWeight;
  a1Matrix(0,2) = f1f3divWeight;
  a1Matrix(0,3) = f1f4divWeight;
  a1Matrix(1,0) = yf2divWeight;
  a1Matrix(1,1) = f2f2divWeight;
  a1Matrix(1,2) = f2f3divWeight;
  a1Matrix(1,3) = f2f4divWeight;
  a1Matrix(2,0) = yf3divWeight;
  a1Matrix(2,1) = f2f3divWeight;
  a1Matrix(2,2) = f3f3divWeight;
  a1Matrix(2,3) = f3f4divWeight;
  a1Matrix(3,0) = yf4divWeight;
  a1Matrix(3,1) = f2f4divWeight;
  a1Matrix(3,2) = f3f4divWeight;
  a1Matrix(3,3) = f4f4divWeight;
  
  a2Matrix(0,0) = f1f1divWeight;
  a2Matrix(0,1) = yf1divWeight;
  a2Matrix(0,2) = f1f3divWeight;
  a2Matrix(0,3) = f1f4divWeight;
  a2Matrix(1,0) = f1f2divWeight;
  a2Matrix(1,1) = yf2divWeight;
  a2Matrix(1,2) = f2f3divWeight;
  a2Matrix(1,3) = f2f4divWeight;
  a2Matrix(2,0) = f1f3divWeight;
  a2Matrix(2,1) = yf3divWeight;
  a2Matrix(2,2) = f3f3divWeight;
  a2Matrix(2,3) = f3f4divWeight;
  a2Matrix(3,0) = f1f4divWeight;
  a2Matrix(3,1) = yf4divWeight;
  a2Matrix(3,2) = f3f4divWeight;
  a2Matrix(3,3) = f4f4divWeight;
  
  a3Matrix(0,0) = f1f1divWeight;
  a3Matrix(0,1) = f1f2divWeight;
  a3Matrix(0,2) = yf1divWeight;
  a3Matrix(0,3) = f1f4divWeight;
  a3Matrix(1,0) = f1f2divWeight;
  a3Matrix(1,1) = f2f2divWeight;
  a3Matrix(1,2) = yf2divWeight;
  a3Matrix(1,3) = f2f4divWeight;
  a3Matrix(2,0) = f1f3divWeight;
  a3Matrix(2,1) = f2f3divWeight;
  a3Matrix(2,2) = yf3divWeight;
  a3Matrix(2,3) = f3f4divWeight;
  a3Matrix(3,0) = f1f4divWeight;
  a3Matrix(3,1) = f2f4divWeight;
  a3Matrix(3,2) = yf4divWeight;
  a3Matrix(3,3) = f4f4divWeight;
  
  a4Matrix(0,0) = f1f1divWeight;
  a4Matrix(0,1) = f1f2divWeight;
  a4Matrix(0,2) = f1f3divWeight;
  a4Matrix(0,3) = yf1divWeight;
  a4Matrix(1,0) = f1f2divWeight;
  a4Matrix(1,1) = f2f2divWeight;
  a4Matrix(1,2) = f2f3divWeight;
  a4Matrix(1,3) = yf2divWeight;
  a4Matrix(2,0) = f1f3divWeight;
  a4Matrix(2,1) = f2f3divWeight;
  a4Matrix(2,2) = f3f3divWeight;
  a4Matrix(2,3) = yf3divWeight;
  a4Matrix(3,0) = f1f4divWeight;
  a4Matrix(3,1) = f2f4divWeight;
  a4Matrix(3,2) = f3f4divWeight;
  a4Matrix(3,3) = yf4divWeight;
  
  a1 = -10000.;
  a2 = -10000.;
  a3 = -10000.;
  a4 = -10000.;
  Double_t delta = deltaMatrix.Determinant();
  
  if(delta != 0.){
    a1 = a1Matrix.Determinant()/delta;
    a2 = a2Matrix.Determinant()/delta;
    a3 = a3Matrix.Determinant()/delta;
    a4 = a4Matrix.Determinant()/delta;
  }

  //find the chi^2 value for the fit
  chiSq = 0.;
  Double_t chi = 0.;
  for(Int_t j = 0; j<nPoints; j++){
    
    chi = (y[j] - (a1 + a2*x[j] + a3*(x[j]*x[j]) + a4*(x[j]*x[j]*x[j])))/weight[j];
    chiSq += chi*chi;
  
  }  
  
  //get the chiSq per degree of freedom

  chiSq /= (1.*nPoints);

  return;
}

//----------------------------------------------------------------------
//method to find a least squares fit to a quadratic function of the form
//y = a1 + a2*x + a3*x^2 
void DmxUtilities::FindQuadraticFit(Double_t *x, Double_t *y, Double_t *weight,  
                                    Int_t nPoints, Double_t &a1, Double_t &a2, 
                                    Double_t &a3, Double_t &chiSq)
{

  //use method of determinants to do the least squares fit
  //this is from bevington chapter 7
  //
  //f_1(x) = 1
  //f_2(x) = x
  //f_3(x) = x^2
  
  //need to find the sums
  //SIGMA y_i*f_1(x_i)/weight_i
  //SIGMA y_i*f_2(x_i)/weight_i
  //SIGMA y_i*f_3(x_i)/weight_i
  //SIGMA f_1(x_i)*f_1(x_i)/weight_i
  //SIGMA f_1(x_i)*f_2(x_i)/weight_i
  //SIGMA f_1(x_i)*f_3(x_i)/weight_i
  //SIGMA f_2(x_i)*f_1(x_i)/weight_i
  //SIGMA f_2(x_i)*f_2(x_i)/weight_i
  //SIGMA f_2(x_i)*f_3(x_i)/weight_i
  //SIGMA f_3(x_i)*f_1(x_i)/weight_i
  //SIGMA f_3(x_i)*f_2(x_i)/weight_i
  //SIGMA f_3(x_i)*f_3(x_i)/weight_i

  //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
  //6 instead of 9 variables

  Double_t yf1divWeight = 0.;
  Double_t yf2divWeight = 0.;
  Double_t yf3divWeight = 0.;
  Double_t f1f1divWeight = 0.;
  Double_t f1f2divWeight = 0.;
  Double_t f1f3divWeight = 0.;
  Double_t f2f2divWeight = 0.;
  Double_t f2f3divWeight = 0.;
  Double_t f3f3divWeight = 0.;
  
  for(Int_t i=0; i<nPoints; i++){
    yf1divWeight += y[i]/(weight[i]*weight[i]);
    yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
    yf3divWeight += (y[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f1f1divWeight += 1./(weight[i]*weight[i]);
    f1f2divWeight += x[i]/(weight[i]*weight[i]);
    f1f3divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
    f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
    f2f3divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
    f3f3divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
  }

  //make a bunch of matricies to hold these values
  TMatrix deltaMatrix(3,3);
  TMatrix a1Matrix(3,3);
  TMatrix a2Matrix(3,3);
  TMatrix a3Matrix(3,3);
    
  //fill the matricies
  deltaMatrix(0,0) = f1f1divWeight;
  deltaMatrix(0,1) = f1f2divWeight;
  deltaMatrix(0,2) = f1f3divWeight;
  deltaMatrix(1,0) = f1f2divWeight;
  deltaMatrix(1,1) = f2f2divWeight;
  deltaMatrix(1,2) = f2f3divWeight;
  deltaMatrix(2,0) = f1f3divWeight;
  deltaMatrix(2,1) = f2f3divWeight;
  deltaMatrix(2,2) = f3f3divWeight;
  
  a1Matrix(0,0) = yf1divWeight;
  a1Matrix(0,1) = f1f2divWeight;
  a1Matrix(0,2) = f1f3divWeight;
  a1Matrix(1,0) = yf2divWeight;
  a1Matrix(1,1) = f2f2divWeight;
  a1Matrix(1,2) = f2f3divWeight;
  a1Matrix(2,0) = yf3divWeight;
  a1Matrix(2,1) = f2f3divWeight;
  a1Matrix(2,2) = f3f3divWeight;
  
  a2Matrix(0,0) = f1f1divWeight;
  a2Matrix(0,1) = yf1divWeight;
  a2Matrix(0,2) = f1f3divWeight;
  a2Matrix(1,0) = f1f2divWeight;
  a2Matrix(1,1) = yf2divWeight;
  a2Matrix(1,2) = f2f3divWeight;
  a2Matrix(2,0) = f1f3divWeight;
  a2Matrix(2,1) = yf3divWeight;
  a2Matrix(2,2) = f3f3divWeight;
    
  a3Matrix(0,0) = f1f1divWeight;
  a3Matrix(0,1) = f1f2divWeight;
  a3Matrix(0,2) = yf1divWeight;
  a3Matrix(1,0) = f1f2divWeight;
  a3Matrix(1,1) = f2f2divWeight;
  a3Matrix(1,2) = yf2divWeight;
  a3Matrix(2,0) = f1f3divWeight;
  a3Matrix(2,1) = f2f3divWeight;
  a3Matrix(2,2) = yf3divWeight;
    
  a1 = -10000.;
  a2 = -10000.;
  a3 = -10000.;
  Double_t delta = deltaMatrix.Determinant();
  
  if(delta != 0.){
    a1 = a1Matrix.Determinant()/delta;
    a2 = a2Matrix.Determinant()/delta;
    a3 = a3Matrix.Determinant()/delta;
  }

  //find the chi^2 value for the fit
  chiSq = 0.;
  Double_t chi = 0.;
  for(Int_t j = 0; j<nPoints; j++){
    
    chi = (y[j] - (a1 + a2*x[j] + a3*(x[j]*x[j])))/weight[j];
    chiSq += chi*chi;
  
  }  
  
  //get the chiSq per degree of freedom

  chiSq /= (1.*nPoints);

  return;
}

//-----------------------------------------------------------------------------------------
//method to find a least squares fit to a linear function of the form
//y = a1 + a2*x
void DmxUtilities::FindLinearFit(Double_t *x, Double_t *y, Double_t *weight, 
                                                                 Int_t nPoints, Double_t &a1, Double_t &a2, 
                                                                 Double_t &chiSq)
{

        //use method of determinants to do the least squares fit
        //this is from bevington chapter 6
        //
        //f_1(x) = 1
        //f_2(x) = x
  
        //need to find the sums
        //SIGMA y_i*f_1(x_i)/weight_i
        //SIGMA y_i*f_2(x_i)/weight_i
        //SIGMA f_1(x_i)*f_1(x_i)/weight_i
        //SIGMA f_1(x_i)*f_2(x_i)/weight_i
        //SIGMA f_2(x_i)*f_1(x_i)/weight_i
        //SIGMA f_2(x_i)*f_2(x_i)/weight_i
        //
        //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
        //3 instead of 4 variables

        Double_t yf1divWeight = 0.;
        Double_t yf2divWeight = 0.;
        Double_t f1f1divWeight = 0.;
        Double_t f1f2divWeight = 0.;
        Double_t f2f2divWeight = 0.;
  
        for(Int_t i=0; i<nPoints; ++i){
                MSG("DmxUtilities", Msg::kDebug) << " x = " << x[i]
                                                                                << " y = " << y[i] 
                                                                                << " weight = " << weight[i] << endl;
                if(weight[i]>0.){
                        yf1divWeight += y[i]/(weight[i]*weight[i]);
                        yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
                        f1f1divWeight += 1./(weight[i]*weight[i]);
                        f1f2divWeight += x[i]/(weight[i]*weight[i]);
                        f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);                     
                }
        }

        //make a bunch of matricies to hold these values
        TMatrix deltaMatrix(2,2);
        TMatrix a1Matrix(2,2);
        TMatrix a2Matrix(2,2);
  
        //fill the matricies
        deltaMatrix(0,0) = f1f1divWeight;
        deltaMatrix(0,1) = f1f2divWeight;
        deltaMatrix(1,0) = f1f2divWeight;
        deltaMatrix(1,1) = f2f2divWeight;
  
        a1Matrix(0,0) = yf1divWeight;
        a1Matrix(0,1) = f1f2divWeight;
        a1Matrix(1,0) = yf2divWeight;
        a1Matrix(1,1) = f2f2divWeight;
  
        a2Matrix(0,0) = f1f1divWeight;
        a2Matrix(0,1) = yf1divWeight;
        a2Matrix(1,0) = f1f2divWeight;
        a2Matrix(1,1) = yf2divWeight;
  
  
        a1 = -10000.;
        a2 = -10000.;
        Double_t delta = deltaMatrix.Determinant();
  
        if(delta != 0.){
                a1 = a1Matrix.Determinant()/delta;
                a2 = a2Matrix.Determinant()/delta;
        }

        //find the chi^2 value for the fit
        chiSq = 0.;
        Double_t chi = 0.;
        for(Int_t j = 0; j<nPoints; ++j){
                
                chi = (y[j] - (a1 + a2*x[j]))/weight[j];
                chiSq += chi*chi;
  
        }  
  
        //get the chiSq per degree of freedom

        chiSq /= (1.*nPoints);

        return;
}

---