DmxDeMuxCosmicsModule.cxx

Go to the documentation of this file.

//$Id: DmxDeMuxCosmicsModule.cxx,v 1.109 2005/03/26 21:41:48 gmieg Exp $
//
//DmxDeMuxCosmicsModule.cxx
//
//Module for demultiplexing cosmic ray events in the far detector
//
//Author:  B. Rebel 6/2001

#include "DeMux/DmxDeMuxCosmicsModule.h"
#include "DeMux/DmxChiSqrStat.h"
#include "DeMux/DmxUtilities.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCommand.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Algorithm/AlgConfig.h"
#include "CandDigit/CandDeMuxDigitListHandle.h"
#include "CandDigit/CandDigitHandle.h"
#include "CandDigit/CandDeMuxDigitHandle.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "CandData/CandRecord.h"
#include "CandData/CandHeader.h"
#include "Candidate/CandContext.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Conventions/PlaneView.h"
#include "CandTrackSR/CandTrackSR.h"
#include "Navigation/XxxItr.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "DeMux/DmxShowerPlane.h"
#include "DeMux/DmxMuonPlane.h"
#include "DeMux/DmxPlane.h"
#include "DeMux/DmxHypothesis.h"
#include "TH1.h"
#include "TH2.h"
#include "TClonesArray.h"
#include "RawData/RawRecord.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "Validity/VldContext.h"
#include "UgliGeometry/UgliGeomHandle.h"

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

#include <cassert>

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

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

//define the NavKey to sort the planes
static NavKey KeyOnPlaneNumber( const TObject *tobj){
  const DmxPlane *currentPlane = dynamic_cast<const DmxPlane *>(tobj);
  return currentPlane->GetPlaneNumber();
}

ClassImp(DmxDeMuxCosmicsModule)

CVSID("$Id: DmxDeMuxCosmicsModule.cxx,v 1.109 2005/03/26 21:41:48 gmieg Exp $");

// Declare this module to JobControll. Arguments are:
//  (1) The class name 
//  (2) The human-readable name 
//  (3) A short, human-readable description of what the module does
JOBMODULE(DmxDeMuxCosmicsModule, 
          "DeMuxCosmicsModule",
          "A module used for demultiplexing cosmics in the far detector");

//......................................................................
DmxDeMuxCosmicsModule::DmxDeMuxCosmicsModule() :
  fAverageTimingOffset(-5.),
  fCrate(-1),
  fDataType(0),
  fDeMuxedStrip(-1),
  fDeMuxedStripW(-1),
  fDigitCharge(-1.),
  fDigitChargeW(-1.),
  fDigitDeMuxedCorrectly(0),
  fDigitsInEvent(0),
  fDigitSide(0),
  fDirection(0),
  fDontUseCandDigitMasks(false),
  fEndPlane(-1),
  fEventCharge(0.),
  fEventDeMuxed(0),
  fEventLength(false),
  fEventNumber(0),
  fFilter(false),
  fHistoFileName("demux_cosmics.root"),
  fHoughSlope(5.),
  fHoughIntercept(10.),
  fHoughPeak(0.66),
  fHypothesisCoG(-1.),
  fHypothesisMatedSignalRatio(-1.),
  fHypothesisNumber(-1),
  fHypothesisRank(-1),
  fHypothesisStat(-1.),
  fHypothesisTieBreaker(-1.),
  fHypothesisValid(-1),
  fMakeTrees(true),
  fMated(0),
  fMatedSignalForValid(0.5),
  fMCDigitChargeE(0.),
  fMCDigitChargeW(0.),
  fMCVertexPlane(0),
  fMCEndPlane(0),
  fMultiple(0),
  fNonPhysicalFailure(0),
  fNonPhysicalPlanes(0),
  fNoVertexFailure(0),
  fNumberOfPlanesInEvent(0),
  fOneSidedUPlanes(0),
  fOneSidedVPlanes(0),
  fPlane(-1),
  fPlaneDeMuxedCoG(-1.),
  fPlaneStray(0),
  fPlanesInWindow(6),
  fPlaneType(-1),
  fSignalFractionLimit(0.1),
  fSingleModuleDeMux(false),
  fSingleModuleCharge(-1.),
  fSingleModuleCorrectCharge(-1.),
  fSingleModuleDigits(-1.),
  fSingleModuleCorrectDigits(-1.),
  fStrayDeltaStripLimit(6),
  fStrip0(0),
  fStrip1(0),
  fStrip2(0),
  fStrip3(0),
  fStrip4(0),
  fStrip5(0),
  fStrip6(0),
  fStrip7(0),
  fTime(-1.),
  fTruthStripE(-1),
  fTruthStripW(-1),
  fUStrayPlanes(0),
  fUValidPlanes(0),
  fVStrayPlanes(0),
  fVValidPlanes(0),
  fXTalkDigit(1),
  fUnsetDigitCheck(false),
  fValidPlanesFailure(0),
  fVertexPlane(-1),
  fWriteHistos(false),
  fXTalkSignalRatio(-1.),
  fUseStrayPlaneCheck(false),
  fStrayPlanesLimit(4),
  fDeMuxFile(0),
  fDeMuxTree(0),
  fDeMuxedDigitTree(0),
  fDmxPlaneTree(0),
  fSingleModuleTree(0),
  fTrackDigitsTree(0),
  fUnusedDigitsTree(0),
  fMCDeMuxTree(0),
  fMCDeMuxDigitTree(0),
  fTFolder(0)
{
 
  // Put DmxStatus object in /Loon/DeMux folder
  if (fTFolder==0) {
    TFolder *lf = dynamic_cast<TFolder *>
                         (gROOT->GetRootFolder()->FindObject("Loon"));
    if (lf==0) {
      lf = gROOT->GetRootFolder()->AddFolder("Loon", "Loon analysis");
      gROOT->GetListOfBrowsables()->Add(lf, "Loon");
    }
    fTFolder = lf->AddFolder("DeMux", "DeMux analysis");
    fTFolder->Add(&fStatus);
  }

   AlgFactory &af = AlgFactory::GetInstance();
   af.Register("AlgDeMuxCosmics", "default", "libDeMux.so", "AlgConfig");
   AlgHandle ah = af.GetAlgHandle("AlgDeMuxCosmics", "default");

   //get the AlgConfig class and set some variables for demuxing
   AlgConfig &acd = ah.GetAlgConfig();
   acd.UnLockValues();
   acd.Set("HoughInterceptRMSLimit",10.);
   acd.Set("HoughPeakLimit",0.60);
   acd.Set("HoughSlopeRMSLimit",10.);
   acd.Set("HypothesisSize",24);
   acd.Set("MuonTrackSlopeLimit",192.0);
   acd.Set("NumberOfHypotheses",169);
   acd.Set("NumberOfStrips", 192);
   acd.Set("PlanesInSet",fPlanesInWindow);
   acd.Set("RatioMatedSignalForValid",fMatedSignalForValid);
   acd.Set("UseCandDigitMask", 1);
   acd.Set("XTalkFractionLimit", 0.1);
   acd.Set("IgnoreSignalLimit", 1.5);
   if(fDataType == 1) acd.Set("AveragePEGainConversion", 60.);
   else acd.Set("AveragePEGainConversion", 60.);
   acd.Set("StrayDeltaStripLimit",fStrayDeltaStripLimit);
   acd.Set("StrayPlanesLimit",10);
   acd.Set("UseStrayPlaneCheck", 0);
   
   if(fUseStrayPlaneCheck){
     acd.Set("UseStrayPlaneCheck",(Int_t)fUseStrayPlaneCheck);
     //acd.SetStrayPlanesLimit(fStrayPlanesLimit);
   }
   
   if( fDontUseCandDigitMasks ){ acd.Set("UseCandDigitMask",0); }
   //MSG("DMX", Msg::kDebug) << fSignalFractionLimit << endl;
   if( fSignalFractionLimit > 0.){ acd.Set("XTalkFractionLimit",fSignalFractionLimit);}
   
   acd.LockValues();
   acd.LockKeys();
   
   MSG("JobC", Msg::kDebug) << "DmxDeMuxCosmicsModule::Constructor" << endl;
   //MSG("DMX", Msg::kDebug) << fSignalFractionLimit << endl;
   
}

//----------------------------------------------------------------------
DmxDeMuxCosmicsModule::~DmxDeMuxCosmicsModule()
{
  MSG("JobC", Msg::kDebug) << "DmxDeMuxCosmicsModule::Destructor" << endl;
  if(fDeMuxFile) delete fDeMuxFile;
}

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

void DmxDeMuxCosmicsModule::BeginJob()
{ 

  if(fMakeTrees && !fDeMuxFile){
    // save the current working directory
    TDirectory* savedir  = gDirectory;  

    //create the new TFile for demuxing
    // this changes the value of gDirectory
    fDeMuxFile = new TFile(fHistoFileName,"RECREATE");  
    
    // create TTree's, these will attach themselves to the current 
    //working directory (demux_cosmics.root)
   
    fDeMuxTree = new TTree("fDeMuxTree", "DeMux Info Tree");
    fDeMuxedDigitTree = new TTree("fDeMuxedDigitTree", "DeMuxed Digit Info");
    fDmxPlaneTree = new TTree("fDmxPlaneTree", "DeMux Plane Info");
    
    fDeMuxTree->Branch("Event", &fEventNumber, "fEventNumber/I"); 
    fDeMuxTree->Branch("DeMuxed", &fEventDeMuxed, "fEventDeMuxed/I");
    fDeMuxTree->Branch("Vertex", &fVertexPlane, "fVertexPlane/I");
    fDeMuxTree->Branch("VertexZ", &fVertexPlaneZ, "fVertexPlaneZ/F");
    fDeMuxTree->Branch("Direction", &fDirection, "fDirection/I");
    fDeMuxTree->Branch("End", &fEndPlane, "fEndPlane/I");
    fDeMuxTree->Branch("NumberOfPlanesHit", &fNumberOfPlanesInEvent, "fNumberOfPlanesInEvent/I");
    fDeMuxTree->Branch("EventCharge", &fEventCharge, "fEventCharge/F");
    fDeMuxTree->Branch("DigitsInEvent", &fDigitsInEvent, "fDigitsInEvent/I");
    fDeMuxTree->Branch("Multiple", &fMultiple, "fMultiple/I");
    fDeMuxTree->Branch("NonPhysicalPlanes", &fNonPhysicalPlanes, "fNonPhysicalPlanes/I");
    fDeMuxTree->Branch("ValidPlanesFailure", &fValidPlanesFailure, "fValidPlanesFailure/I");
    fDeMuxTree->Branch("NoVertexFailure", &fNoVertexFailure, "fNoVertexFailure/I");
    fDeMuxTree->Branch("NonPhysicalFailure", &fNonPhysicalFailure, "fNonPhysicalFailure/I");
    fDeMuxTree->Branch("OneSidedUPlanes", &fOneSidedUPlanes, "fOneSidedUPlanes/I");
    fDeMuxTree->Branch("OneSidedVPlanes", &fOneSidedVPlanes, "fOneSidedVPlanes/I");
    fDeMuxTree->Branch("SingleModuleCharge", &fSingleModuleCharge, "fSingleModuleCharge/F");
    fDeMuxTree->Branch("SingleModuleCorrectCharge", &fSingleModuleCorrectCharge, "fSingleModuleCorrectCharge/F");
    fDeMuxTree->Branch("SingleModuleDigits", &fSingleModuleDigits, "fSingleModuleDigits/F");
    fDeMuxTree->Branch("SingleModuleCorrectDigits", &fSingleModuleCorrectDigits, "fSingleModuleCorrectDigits/F");
    fDeMuxTree->Branch("UStrayPlanes", &fUStrayPlanes, "fUStrayPlanes/I");
    fDeMuxTree->Branch("VStrayPlanes", &fVStrayPlanes, "fVStrayPlanes/I");
    fDeMuxTree->Branch("UValidPlanes", &fUValidPlanes, "fUValidPlanes/I");
    fDeMuxTree->Branch("VValidPlanes", &fVValidPlanes, "fVValidPlanes/I");
    fDeMuxTree->Branch("UOverlappingMultiple", &fUOverlap, "fUOverlap/I");
    fDeMuxTree->Branch("VOverlappingMultiple", &fVOverlap, "fVOverlap/I");
    fDeMuxTree->Branch("AvTimingOffset", &fAverageTimingOffset, "fAverageTimingOffset/F");
    //fDeMuxTree->SetAutoSave(10000);
    
    fDeMuxedDigitTree->Branch("Event", &fEventNumber, "fEventNumber/I");
    fDeMuxedDigitTree->Branch("Plane", &fPlane, "fPlane/I");
    fDeMuxedDigitTree->Branch("Crate", &fCrate, "fCrate/I");
    fDeMuxedDigitTree->Branch("ADC", &fDigitCharge, "fDigitCharge/F");
    fDeMuxedDigitTree->Branch("Time", &fTime, "fTime/D");
    fDeMuxedDigitTree->Branch("Strip0", &fStrip0, "fStrip0/I");
    fDeMuxedDigitTree->Branch("Strip1", &fStrip1, "fStrip1/I");
    fDeMuxedDigitTree->Branch("Strip2", &fStrip2, "fStrip2/I");
    fDeMuxedDigitTree->Branch("Strip3", &fStrip3, "fStrip3/I");
    fDeMuxedDigitTree->Branch("Strip4", &fStrip4, "fStrip4/I");
    fDeMuxedDigitTree->Branch("Strip5", &fStrip5, "fStrip5/I");
    fDeMuxedDigitTree->Branch("Strip6", &fStrip6, "fStrip6/I");
    fDeMuxedDigitTree->Branch("Strip7", &fStrip7, "fStrip7/I");
    fDeMuxedDigitTree->Branch("DeMuxedStrip", &fDeMuxedStrip, "fDeMuxedStrip/I");
    fDeMuxedDigitTree->Branch("XTalkDigit", &fXTalkDigit, "fXTalkDigit/I");
    //fDeMuxedDigitTree->SetAutoSave(100000);
    
    fDmxPlaneTree->Branch("Event", &fEventNumber, "fEventNumber/I");
    fDmxPlaneTree->Branch("Plane", &fPlane, "fPlane/I");
    fDmxPlaneTree->Branch("PlaneZ", &fPlaneZ, "fPlaneZ/F");
    fDmxPlaneTree->Branch("PlaneType", &fPlaneType, "fPlaneType/I");
    fDmxPlaneTree->Branch("PlaneCoG", &fPlaneDeMuxedCoG, "fPlaneDeMuxedCoG/F");
    fDmxPlaneTree->Branch("Stray", &fPlaneStray, "fPlaneStray/I");
    fDmxPlaneTree->Branch("HypothesisCoG", &fHypothesisCoG, "fHypothesisCoG/F");
    fDmxPlaneTree->Branch("HypothesisCoGTPos", &fHypothesisCoGTPos, "fHypothesisCoGTPos/F");
    fDmxPlaneTree->Branch("Hypothesis", &fHypothesisNumber, "fHypothesisNumber/I");
    fDmxPlaneTree->Branch("MatedSignal", &fHypothesisMatedSignalRatio, "fHypothesisMatedSignalRatio/F");
    fDmxPlaneTree->Branch("Statistic", &fHypothesisStat, "fHypothesisStat/F");
    fDmxPlaneTree->Branch("TieBreaker", &fHypothesisTieBreaker, "fHypothesisTieBreaker/F");
    fDmxPlaneTree->Branch("Rank", &fHypothesisRank, "fHypothesisRank/I");
    fDmxPlaneTree->Branch("Valid", &fHypothesisValid, "fHypothesisValid/I");
    //fDmxPlaneTree->SetAutoSave(10000);

    fMCDeMuxDigitTree = new TTree("fMCDeMuxDigitTree", "MC DeMux Digit Tree");
    fMCDeMuxDigitTree->Branch("Event", &fEventNumber, "fEventNumber/I");
    fMCDeMuxDigitTree->Branch("TruthDigitChargeE", &fMCDigitChargeE, "fMCDigitChargeE/F");
    fMCDeMuxDigitTree->Branch("TruthDigitChargeW", &fMCDigitChargeW, "fMCDigitChargeW/F");
    fMCDeMuxDigitTree->Branch("DigitChargeE", &fDigitCharge, "fDigitCharge/F");
    fMCDeMuxDigitTree->Branch("DigitChargeW", &fDigitChargeW, "fDigitChargeW/F");
    fMCDeMuxDigitTree->Branch("Plane", &fPlane, "fPlane/I");
    fMCDeMuxDigitTree->Branch("Strip", &fDeMuxedStrip, "fDeMuxedStrip/I");
    fMCDeMuxDigitTree->Branch("MCStripSpread", &fMCStripSpread, "fMCStripSpread/I");
    //fMCDeMuxDigitTree->SetAutoSave(10000);
        
    fMCDeMuxTree = new TTree("fMCDeMuxTree", "MC DeMux Tree");
    fMCDeMuxTree->Branch("Event", &fEventNumber, "fEventNumber/I"); 
    fMCDeMuxTree->Branch("Vertex", &fMCVertexPlane, "fVertexPlane/I");
    fMCDeMuxTree->Branch("End", &fMCEndPlane, "fEndPlane/I");
    fMCDeMuxTree->Branch("ChargeCorrect", &fMCChargeCorrect, "fMCChargeCorrect/F");
    fMCDeMuxTree->Branch("ChargeCorrect1Hyp", &fMCChargeCorrect1Hyp, "fMCChargeCorrect1Hyp/F");
    fMCDeMuxTree->Branch("TotalCharge", &fMCTotalCharge, "fMCTotalCorrect/F");
    fMCDeMuxTree->Branch("TotalCharge1Hyp", &fMCTotalCharge1Hyp, "fMCTotalCharge1Hyp/F");
    fMCDeMuxTree->Branch("StripsCorrect", &fMCStripsCorrect, "fMCStripsCorrect/F");
    fMCDeMuxTree->Branch("StripSpread24Planes", &fMCStripSpread24Planes, "fMCStripSpread24Planes/I");
    fMCDeMuxTree->Branch("ShowerSpread", &fMCShowerSpread, "fMCShowerSpread/I");
    
    
    //fMCDeMuxTree->SetAutoSave(10000);
    
    // change back to current working directory before leaving constructor
    savedir->cd();   
  }

  //set up the TTrees for the various tests
  if( fSingleModuleDeMux && fMakeTrees){
    
    // save the current working directory
    TDirectory* savedir = gDirectory; 
    fDeMuxFile->cd();

    fSingleModuleTree = new TTree("fSingleModuleTree", "SingleModule Test Tree");
    fSingleModuleTree->Branch("Event", &fEventNumber, "fEventNumber/I"); 
    fSingleModuleTree->Branch("DigitCharge", &fDigitCharge, "fDigitCharge/F");
    fSingleModuleTree->Branch("XTalkSignalRatio", &fXTalkSignalRatio, "fXTalkSignalRatio/D");
    fSingleModuleTree->Branch("Strip", &fDeMuxedStrip, "fDeMuxedStrip/I");
    fSingleModuleTree->Branch("Mated", &fMated, "fMated/I");
    fSingleModuleTree->Branch("DigitDeMuxedCorrectly", &fDigitDeMuxedCorrectly, "fDigitDeMuxedCorrectly/I");
    
    // change back to current working directory
    savedir->cd(); 
  }
  if( fUnsetDigitCheck && fMakeTrees){

    // save the current working directory
    TDirectory* savedir = gDirectory; 
    fDeMuxFile->cd();

    fUnusedDigitsTree = new TTree("fUnusedDigitsTree", "UnsetDigitCheck Test Tree");
    fUnusedDigitsTree->Branch("Event", &fEventNumber, "fEventNumber/I"); 
    fUnusedDigitsTree->Branch("DigitCharge", &fDigitCharge, "fDigitCharge/F");
    fUnusedDigitsTree->Branch("XTalkSignalRatio", &fXTalkSignalRatio, "fXTalkSignalRatio/D");
    fUnusedDigitsTree->Branch("DigitSide", &fDigitSide, "fDigitSide/F");
    
    // change back to current working directory
    savedir->cd(); 
  }
  
  return;
}

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

JobCResult DmxDeMuxCosmicsModule::Reco(MomNavigator *mom)
{
  JobCResult result(JobCResult::kPassed);

  //MSG("JobCDmx", Msg::kDebug) << "DeMuxCosmicsModule::Reco starting" << endl;

  // Check that mom exists.
  assert(mom);

  // Find PrimaryCandidateRecord fragment in MOM.
  CandRecord *candrec = dynamic_cast<CandRecord *>
    (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if (candrec == 0) {
    result.SetError().SetFailed();
    return result;
  }

  //MSG("JobCDmx", Msg::kDebug) << "got cand record" << endl;
  //*-* Find Raw Data CandDeMuxDigitList fragment in PrimaryCandidateRecord.
  CandDeMuxDigitListHandle *crdlh = dynamic_cast<CandDeMuxDigitListHandle *>
    (candrec->FindCandHandle("CandDeMuxDigitListHandle", "canddigitlist"));
  if (crdlh == 0) {
    result.SetError().SetFailed();
    return result;
  }
    
  AlgFactory &af = AlgFactory::GetInstance();
  AlgHandle ah = af.GetAlgHandle("AlgDeMuxCosmics", "default");
  //AlgHandle ah = af.GetAlgHandle("AlgDeMuxGolden", "default");
    
  CandContext cx(this, mom);
  //cx.SetDataIn(status);
  
  //MSG("JobCDmx", Msg::kDebug) << "starting demux algorithm " << snarl << endl;
  ah.RunAlg(*crdlh, cx);
  //MSG("JobCDmx", Msg::kDebug) << "finished demux algorithm " << snarl << endl;
  
  
  //MSG("JobCDmx", Msg::kDebug) << "Event DeMuxed" << endl;
  
  return result;
}

//......................................................................
//this method fills various histograms used to test the demuxing algorithm
//efficiency

JobCResult DmxDeMuxCosmicsModule::Ana(const MomNavigator *mom)
{
  JobCResult result(JobCResult::kPassed);
  
  //MSG("JobCDmx", Msg::kDebug) << "DeMuxCosmicsModule::Ana" << endl;

  DmxStatus *status = dynamic_cast<DmxStatus *>(gROOT->GetRootFolder()->FindObject("Loon/DeMux/DmxStatus"));
  if (status == 0) {
    MSG("DMX", Msg::kWarning) << "//root/Loon/DeMux/DmxStatus not found." << endl;
    
    result.SetError().SetFailed();
    return result;
  }
  
  //if the event fails demuxing for any reason and filtering is turned on, 
  //dont pass the event on
  
  if(fFilter && fUseStrayPlaneCheck && (status->GetVOverlappingMultiple() || status->GetUOverlappingMultiple())){
    result.SetError().SetFailed();
    return result;
  }
  else if( !status->GetEventDeMuxed() && fFilter){
    result.SetError().SetFailed();
    return result;
  }
  
  //only do stuff in this method if you are making trees - otherwise there 
  //is no point to running the method.
  if( fMakeTrees ){
    
          DmxUtilities util;

          // Check that mom exists.
          assert(mom);
    
          // Find PrimaryCandidateRecord fragment in MOM.
          CandRecord *candrec = dynamic_cast<CandRecord *>
                  (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
          if(candrec == 0){
                  result.SetError().SetFailed();
                  return result;
          }
    
    //Find Raw Data CandDeMuxDigitList fragment in PrimaryCandidateRecord.
    CandDeMuxDigitListHandle *crdlh = dynamic_cast<CandDeMuxDigitListHandle *>
      (candrec->FindCandHandle("CandDeMuxDigitListHandle", "canddigitlist"));
    
    // Check that crdlh exists.
    assert(crdlh);
        
        //check that there are digits present
        //MSG("DeMuxCosmics", Msg::kWarning) << crdlh->GetNDaughters() << " digits present at start of ana" << endl;
        
    
    //get an iterator over the digits and one extra for the single module test
    CandDeMuxDigitHandleItr crdhi(crdlh->GetDaughterIterator());
    CandDeMuxDigitHandleItr smcrdhi(crdlh->GetDaughterIterator());

    //create the KeyFunc
    CandDeMuxDigitHandleKeyFunc *planeKF = crdhi.CreateKeyFunc();
    CandDeMuxDigitHandleKeyFunc *smPlaneKF = smcrdhi.CreateKeyFunc();
    
    //program the KeyFunc with the sort function
    planeKF->SetFun(KeyFromDigitPlane);
    smPlaneKF->SetFun(KeyFromDigitPlane);
    
    //get the NavSet from the iterator and pass the KeyFunc to it
    crdhi.GetSet()->AdoptSortKeyFunc(planeKF);
    smcrdhi.GetSet()->AdoptSortKeyFunc(smPlaneKF);
    
    //clear the KF pointer because we no longer own the KeyFunc
    planeKF = 0;
    smPlaneKF = 0;
    
    crdhi.ResetFirst();
    smcrdhi.ResetFirst();

    //get the TObjArray of DmxPlane objects
    const TObjArray *planeArray = status->GetPlaneArray();
    
    //create a TObjectItr
    TObjectItr planeItr(planeArray);
    
    //create a KeyFunc to sort planes by number
    TObjectKeyFunc *pnKF = planeItr.CreateKeyFunc();
    
    //program the KeyFunc with the sort function
    pnKF->SetFun(KeyOnPlaneNumber);
    
    //get the NavSet from the iterator and pass the KeyFunc to it
    planeItr.GetSet()->AdoptSortKeyFunc(pnKF);
    
    //clear the KF pointer because we no longer own the KeyFunc
    pnKF = 0;
    
    fNonPhysicalPlanes = 0;
    fHypothesisCoG = -1.;
    fPlaneStray = 0;
    fPlaneDeMuxedCoG = -1.;
    fHypothesisNumber = -1;
    fHypothesisStat = -1.;
    fHypothesisTieBreaker = -1.;
    fHypothesisMatedSignalRatio = -1.;
    fHypothesisValid = -1;
    fPlaneType = -1;
    fCrate = -1;
    fPlane = -1;
    fDeMuxedStrip = -1;
    fOneSidedUPlanes = 0;
    fOneSidedVPlanes = 0;
    fSingleModuleCharge = -1.;
    fSingleModuleDigits = -1.;
    fSingleModuleCorrectCharge = -1.;
    fSingleModuleCorrectDigits = -1.;
    fEventCharge = 0.;
    fDigitsInEvent = 0;
    fDigitCharge = -1.;
    fXTalkSignalRatio = -1.;
    fDigitDeMuxedCorrectly = 0;
    fDigitSide = 0;
    //get generic info for event
    fEventNumber = status->GetEventNumber();
    fEventDeMuxed = (Int_t)status->GetEventDeMuxed();
    fValidPlanesFailure = (Int_t)status->GetValidPlanesFailure();
    fNoVertexFailure = (Int_t)status->GetVertexPlaneFailure();
    fNonPhysicalFailure = (Int_t)status->GetNonPhysicalFailure();
    fUOverlap = (Int_t)status->GetUOverlappingMultiple();
    fVOverlap = (Int_t)status->GetVOverlappingMultiple();
    fAverageTimingOffset = status->GetAverageTimingOffset();
    fUSlope = status->GetUTrackSlope();
    fUIntercept = status->GetUTrackIntercept();
    fUSlopeRMS = status->GetUSlopeRMS();
    fUStrayPlanes = status->GetUStrayPlanes();
    fUInterceptRMS = -1.;
    fUChiSq = status->GetUTrackChiSq();
    fVSlope = status->GetVTrackSlope();
    fVIntercept = status->GetVTrackIntercept();
    fVSlopeRMS = status->GetVSlopeRMS();
    fVStrayPlanes = status->GetVStrayPlanes();
    fVInterceptRMS = -1.;
    fVChiSq = status->GetVTrackChiSq();
    fUValidPlanes = 0;
    fVValidPlanes = 0;
  
    //zero values that need to start over with each event
    fVertexPlane = status->GetVertexPlaneNumber();
    fVertexPlaneZ = status->GetVertexZPosition();
    fDirection = status->GetEventDirection();
    fEndPlane = status->GetEndPlaneNumber();
    fMultiple = (Int_t)status->GetMultipleMuon();
    
    if(fMultiple == 0){planeItr.GetSet()->Slice(fVertexPlane, fEndPlane);}
    else if(fMultiple == 1){planeItr.GetSet()->Slice(fVertexPlane, 500);}
    
    //-----------------------------get plane info--------------------------//
    
    //get the total charge and digits in the event between vertex and end planes
    while( planeItr.IsValid() ){
      DmxPlane *currentPlane = dynamic_cast<DmxPlane *>(planeItr.Ptr());
      fEventCharge += currentPlane->GetPlaneCharge();
      fPlane = currentPlane->GetPlaneNumber();
      fPlaneZ = currentPlane->GetZPosition();
      fPlaneDeMuxedCoG = currentPlane->GetStripCoG(); 
      fPlaneStray = (Int_t)currentPlane->IsStray();
      if(currentPlane->GetPlaneView() == PlaneView::kU){
        if(currentPlane->IsValid()){ ++fUValidPlanes; }
      }
      else if(currentPlane->GetPlaneView() == PlaneView::kV){
        if(currentPlane->IsValid()){ ++fVValidPlanes; }
      }
      //if(fPlaneFitCoG > 192. || fPlaneFitCoG < 0.){ ++fNonPhysicalPlanes; }
      if(TMath::Abs(fPlaneDeMuxedCoG) > 4.){ ++fNonPhysicalPlanes; }
      
      crdhi.GetSet()->Slice(fPlane);
      
      //fill the DmxPlane tree
      if( currentPlane->IsValid() && currentPlane->GetPlaneType() == DmxPlaneTypes::kShower ){
        fPlaneType = 0;
        
        //loop over the best hypotheses to get their stats
        if(dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()){
          fHypothesisNumber = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetLowerBound();
          fHypothesisRank = 1;
          fHypothesisCoG = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetCoG();
          fHypothesisCoGTPos = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestCoG();
          fHypothesisStat = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetStat();
          fHypothesisTieBreaker = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetTieBreakerStat();
          fHypothesisValid = (Int_t)dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->IsValid();
          fHypothesisMatedSignalRatio = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetBestHypothesis()->GetMatedSignalRatio();

          if(status->GetEventDeMuxed()){fDmxPlaneTree->Fill();}

        }
        if(dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis() ){
          fHypothesisNumber = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetLowerBound();
          fHypothesisRank = 2;
          fHypothesisCoG = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetCoG();
          fHypothesisCoGTPos = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestCoG();
          fHypothesisStat = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetStat();
          fHypothesisTieBreaker = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetTieBreakerStat();
          fHypothesisValid = (Int_t)dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->IsValid();
          fHypothesisMatedSignalRatio = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
          if(status->GetEventDeMuxed()){fDmxPlaneTree->Fill();}

        }
        if(dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis() ){
          fHypothesisNumber = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetLowerBound();
          fHypothesisRank = 3;
          fHypothesisCoG = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetCoG();
          fHypothesisCoGTPos = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestCoG();
          fHypothesisStat = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetStat();
          fHypothesisTieBreaker = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetTieBreakerStat();
          fHypothesisValid = (Int_t)dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->IsValid();
          fHypothesisMatedSignalRatio = dynamic_cast<DmxShowerPlane *>(currentPlane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
          if(status->GetEventDeMuxed()){fDmxPlaneTree->Fill();}
        }
        
        //fDmxPlaneTree->Fill();
        fHypothesisStat = -1.;
        fHypothesisMatedSignalRatio = -1.;
        fHypothesisTieBreaker = -1.;
        fHypothesisRank = -1;
        fHypothesisValid = -1;
        fHypothesisCoG = -1.;
      }
      else if( currentPlane->IsValid() && currentPlane->GetPlaneType() == DmxPlaneTypes::kMuon ){
        fPlaneType = 1;
        if((Int_t)currentPlane->GetStripCoG() > 11){
          fHypothesisNumber = (Int_t)currentPlane->GetStripCoG() - 11;
        }
        else{ fHypothesisNumber = 0; }
        
        fHypothesisTieBreaker = -1;
        fHypothesisRank = 1;
        fHypothesisCoGTPos = dynamic_cast<DmxMuonPlane *>(currentPlane)->GetInitialCoG();
        fHypothesisCoG = dynamic_cast<DmxMuonPlane *>(currentPlane)->GetInitialStripCoG();
        fHypothesisMatedSignalRatio = 1.;
        fHypothesisValid = (Int_t)currentPlane->IsValid();
        if(status->GetEventDeMuxed()){
          fDmxPlaneTree->Fill();
          if ( fEventNumber%1000 == 0 ) {
            // save current working directory
            TDirectory* savedir = gDirectory; 
            // move to file directory of tree, oddly SaveSelf 
            //doesn't work without it
            fDeMuxFile->cd();
            // save tree to file  
            fDmxPlaneTree->AutoSave();
            // save file directory containing this tree
            fDeMuxFile->SaveSelf();
            // change back to working to directory
            savedir->cd();  
          }
        }
      }
      
      //-------------------get digit info----------------------------------//
      
      fDigitsInEvent += crdhi.SizeSelect();
      
      //loop over the digits to fill the demuxed digit info
      while( crdhi.IsValid() ){
        CandDeMuxDigitHandle *digit = crdhi.Ptr();
        fStrip0 = 0;
        fStrip1 = 0;
        fStrip2 = 0;
        fStrip3 = 0;
        fStrip4 = 0;
        fStrip5 = 0;
        fStrip6 = 0;
        fStrip7 = 0;
        fXTalkDigit = 1;
        fTime = digit->GetTime();
        //get the strip, charge, crate for the digit
        if(digit->GetPlexSEIdAltL().GetBestWeight() == 1.){
          fDeMuxedStrip = digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
        }
        else{ fDeMuxedStrip = -1;}
        
        //set the alternatives list to the first one
        digit->GetPlexSEIdAltL().SetFirst();
        
        //loop over the strip end alt list
        Int_t cnt = 0;
        while( digit->GetPlexSEIdAltL().IsValid() ){
          if(cnt == 0){ fStrip0 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          else if(cnt == 1){ fStrip1 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          else if(cnt == 2){ fStrip2 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          else if(cnt == 3){ fStrip3 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          else if(cnt == 4){ fStrip4 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          else if(cnt == 5){ fStrip5 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          else if(cnt == 6){ fStrip6 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          else if(cnt == 7){ fStrip7 = digit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();}
          ++cnt;
          digit->GetPlexSEIdAltL().Next();
        }
        
        
        fDigitCharge = digit->GetCharge();
        
        //MSG("DMX", Msg::kDebug) << status->GetEventNumber() << "\t" << fPlane << "\t" << fDigitCharge 
        //                   << "\t" << neighborSignal << "\t" << neighborSignal*fSignalFractionLimit 
        //           << "\t" << fSignalFractionLimit << endl;
        
        if( digit->GetDeMuxDigitFlagWord() == 0) fXTalkDigit = 0;
        
        fCrate = digit->GetChannelId().GetCrate();
        if(status->GetEventDeMuxed()){
          fDeMuxedDigitTree->Fill();
          if ( fEventNumber%1000 == 0 ) {
            TDirectory* savedir = gDirectory; // save current working directory
            fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
            fDeMuxedDigitTree->AutoSave();  // save tree to file
            fDeMuxFile->SaveSelf();  // save file directory containing this tree
            savedir->cd();  // change back to working to directory
          }
        }
        crdhi.Next();
      }
      
      crdhi.GetSet()->Slice();
      crdhi.ResetFirst();
      
      planeItr.Next();
    }
    
    //if you arent using the cand digit masks to demux, then set them all to true
    //if(fDontUseCandDigitMasks){ crdhi.GetSet()->GetMasks().SetAllMasks(kTRUE); }
    
    planeItr.ResetFirst();
    planeItr.GetSet()->Slice();

    Int_t totalPlanes = planeItr.SizeSelect();
    fNumberOfPlanesInEvent = totalPlanes;
    
    //------------------start the efficiency tests----------------------------//
    
    if( fDataType == 0){
      if( fSingleModuleDeMux && status->GetEventDeMuxed() ){
        //get information about plane 10's shared tube 
        //to do some cross talk analysis later
        Float_t west10[16];
        Float_t east10[16];
        for(Int_t i = 0; i < 16; i++){
          west10[i] = 0.;
          east10[i] = 0.;
        }
        
        
        //slice the digit Itr to get plane 10
        crdhi.GetSet()->Slice(10);
        
        //fill the hit info
        util.FillHitPixels(crdhi, 2, east10, west10);
        
        //reset the itr
        crdhi.GetSet()->Slice();
        crdhi.ResetFirst();
        
        //slice digit Itr to get plane 12
        crdhi.GetSet()->Slice(12);
        crdhi.ResetFirst();
        smcrdhi.GetSet()->Slice(12);
        smcrdhi.ResetFirst();

        if(crdhi.SizeSelect()>0){
          fMated = 0;
          fSingleModuleCorrectCharge = 0.;
          fSingleModuleCorrectDigits = 0.;
          fSingleModuleCharge = 0.;
          fSingleModuleDigits = 0.;
        }
        
        //array to keep track of which pixels on which sides were hit for
        //the set digits//chip 2 can have xtalk from both planes 12 and 10, so put that into the signal ratio
        Float_t westSet1[16];
        Float_t eastSet1[16];
        Float_t westSet2[16];
        Float_t eastSet2[16];
        for(Int_t i = 0; i < 16; i++){
          westSet1[i] = 0.;
          eastSet1[i] = 0.;
          westSet2[i] = 0.;
          eastSet2[i] = 0.;
        }
        
        util.FillHitPixels(crdhi, 1, eastSet1, westSet1);
        util.FillHitPixels(crdhi, 2, eastSet2, westSet2);
        
        while( crdhi.IsValid() ){
          CandDeMuxDigitHandle *digit = crdhi.Ptr();
          
          //get the strip for the digit and check that it is within module 3 or module 6
          Int_t strip = digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
          fDigitCharge = digit->GetCharge()/55.;
          
          if( crdhi.Ptr()->GetPlexSEIdAltL().GetBestWeight() == 1. 
              && crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
            
              fSingleModuleCharge += fDigitCharge;
              fSingleModuleDigits += 1.;
              fDeMuxedStrip = strip;
              
              bool mated = false;
              //loop over the smcrdhi to see if there is an east side digit from this strip
              smcrdhi.ResetFirst();
              while( smcrdhi.IsValid() ){
                if(smcrdhi.Ptr()->GetPlexSEIdAltL().GetBestWeight() == 1. &&
                   smcrdhi.Ptr()->GetPlexSEIdAltL().GetBestSEId().GetStrip() == fDeMuxedStrip &&
                   smcrdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast)  mated = true;
                smcrdhi.Next();
              }
              
              if( mated ) fMated = 1;

              Int_t chip = digit->GetChannelId().GetVaChip();
              Int_t channel = digit->GetChannelId().GetVaChannel();
              Int_t pixel = util.GetDigitPixel(channel);
              
              //get the ratio of this digits charge to the sum of its 4 nearest neighbors
              Float_t nearestNeighborSignal = -1.;
              if(chip == 1){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet1);}
              if(chip == 2){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet2) 
                              + util.CheckForXTalk(pixel, west10);}
              
              if(nearestNeighborSignal > 0.){fXTalkSignalRatio = fDigitCharge / (nearestNeighborSignal/55.);}
              
              //MSG("DMX12", Msg::kDebug)<< fDigitCharge << "\t" << nearestNeighborSignal/55. 
              //                  << "\t" << fXTalkSignalRatio << endl;
              
              
              if( strip >= 56 && strip <= 75 ){
                fSingleModuleCorrectDigits += 1.;
                fSingleModuleCorrectCharge += fDigitCharge;
                fDigitDeMuxedCorrectly = 1;
              }
              else if( strip >= 116 && strip <= 135 ){
                fSingleModuleCorrectDigits += 1.;
                fSingleModuleCorrectCharge += fDigitCharge;
                fDigitDeMuxedCorrectly = 1;
              }
              else{
                fDigitDeMuxedCorrectly = 0;
                //MSG("DMX", Msg::kDebug)<<"XTalkSignalRatio = " << fXTalkSignalRatio << "\tNearestNeighborsSignal = "
                //              <<nearestNeighborSignal/55. << "\tDigitCharge = " << fDigitCharge << endl;
              } //end reconstructed in wrong module
              
              //MSG("DMX", Msg::kDebug)<<"SingleModuleCharge = " << fSingleModuleCharge 
              //                << "\tSingleModuleCorrectCharge = " <<fSingleModuleCorrectCharge 
              //                << "\tSingleModuleDigits = " << fSingleModuleDigits 
              //                << "\tSingleModuleCorrectDigits = " << fSingleModuleCorrectDigits 
              //                << endl;
              
              //fill the digit information in the tree
              fSingleModuleTree->Fill();
              if ( fEventNumber%1000 == 0 ) {
                TDirectory* savedir = gDirectory; // save current working directory
                fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
                fSingleModuleTree->AutoSave();  // save tree to file
                fDeMuxFile->SaveSelf();  // save file directory containing this tree
                savedir->cd();  // change back to working to directory
              }
              fXTalkSignalRatio = -1.;
              fMated = 0;
              fDigitDeMuxedCorrectly = 0;
            } //end if from west side and set
          crdhi.Next();
        } //end loop over digits
        
        crdhi.GetSet()->Slice();
        crdhi.ResetFirst();
        smcrdhi.GetSet()->Slice();
        smcrdhi.ResetFirst();
      }
      
      if( fUnsetDigitCheck && status->GetEventDeMuxed() ){
        
        //slice to get planes between the vertex and end
        planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
        planeItr.ResetFirst();
        
        //loop over the planes and print their reconstructions
        while( planeItr.IsValid() ){
          DmxPlane *plane = dynamic_cast<DmxPlane *>(planeItr.Ptr());
          
          Int_t planeNumber = plane->GetPlaneNumber();
          
          //slice the digit iterator to the current plane
          crdhi.GetSet()->Slice(planeNumber);
          
          //array to keep track of which pixels on which planes were hit for
          //the set digits
          Float_t westSet0[16];
          Float_t eastSet0[16];
          Float_t westSet1[16];
          Float_t eastSet1[16];
          Float_t westSet2[16];
          Float_t eastSet2[16];
          for(Int_t i = 0; i < 16; i++){
            westSet0[i] = 0.;
            eastSet0[i] = 0.;
            westSet1[i] = 0.;
            eastSet1[i] = 0.;
            westSet2[i] = 0.;
            eastSet2[i] = 0.;
          }
          
          util.FillHitPixels(crdhi, 0, eastSet0, westSet0);
          crdhi.ResetFirst();
          util.FillHitPixels(crdhi, 1, eastSet1, westSet1);
          crdhi.ResetFirst();
          util.FillHitPixels(crdhi, 2, eastSet2, westSet2);
          
          crdhi.ResetFirst();
          
          //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();
            fDigitCharge = digit->GetCharge()/55.;
            
            Int_t chip = digit->GetChannelId().GetVaChip();
            Int_t pixel = util.GetDigitPixel(digit->GetChannelId().GetVaChannel());
            
            //see if the best weight is < 1.0 - if so, the digit had no strip chosen
            if(digit->GetPlexSEIdAltL().GetBestWeight() < 1.0){
              
              Float_t nearestNeighborSignal = -1.;
              if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
                if(chip == 0){nearestNeighborSignal = util.CheckForXTalk(pixel, eastSet0);}
                if(chip == 1){nearestNeighborSignal = util.CheckForXTalk(pixel, eastSet1);}
                if(chip == 2){nearestNeighborSignal = util.CheckForXTalk(pixel, eastSet2);}
                fDigitSide = 0;
              }
              if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
                if(chip == 0){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet0);}
                if(chip == 1){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet1);}
                if(chip == 2){nearestNeighborSignal = util.CheckForXTalk(pixel, westSet2);}
                fDigitSide = 1;
              }
              
              if(nearestNeighborSignal > 0.){fXTalkSignalRatio = fDigitCharge / (nearestNeighborSignal/55.);}
              
              fUnusedDigitsTree->Fill();
              if ( fEventNumber%1000 == 0 ) {
                TDirectory* savedir = gDirectory; // save current working directory
                fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
                fUnusedDigitsTree->AutoSave();  // save tree to file
                fDeMuxFile->SaveSelf();  // save file directory containing this tree
                savedir->cd();  // change back to working to directory
              }

              fXTalkSignalRatio = -1.;
            }
            crdhi.Next();
          }
          
          //clear the digit slice
          crdhi.GetSet()->Slice();
          planeItr.Next();
        }
        
        //reset the digit itr
        crdhi.GetSet()->Slice();
        crdhi.ResetFirst();
        
        //reset the planeItr
        planeItr.GetSet()->Slice();
        planeItr.ResetFirst();
        
      } //end UnsetDigitCheck test
      
      if( fEventLength ){
        
        if( status->GetValidPlanesFailure() ){ 
          //MSG("DmxDigit", Msg::kDebug) << "Event = " << status->GetEventNumber() << endl;
          
          while( planeItr.IsValid() ){
            
            crdhi.GetSet()->Slice(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber());
            
            //declare variables to see if the digits from this plane come from each side
            Bool_t westDigits = false;
            Bool_t eastDigits = false;
            
            while( crdhi.IsValid() ){
              
              if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest ){
                //  MSG("DmxDigit", Msg::kDebug) << "\tPlane = " 
                //                                        << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber() 
                //                                        << "\tSide = West\t" << "chip = " 
                //                                        << crdhi.Ptr()->GetChannelId().GetVaChip()
                //                                        << "\tpixel = " 
                //                                        << status->GetDigitPixel(crdhi.Ptr()->GetChannelId().GetVaChannel()) 
                //                                        << "\tcharge = " << crdhi.Ptr()->GetCharge() << endl;
                westDigits = true;
              }
              if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast ){
                //  MSG("DmxDigit", Msg::kDebug) << "\tPlane = " 
                //                                        << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber() 
                //                                        << "\tSide = East\t" << "chip = " 
                //                                        << crdhi.Ptr()->GetChannelId().GetVaChip()
                //                                        << "\tpixel = " 
                //                                        << status->GetDigitPixel(crdhi.Ptr()->GetChannelId().GetVaChannel()) 
                //                                        << "\tcharge = " << crdhi.Ptr()->GetCharge() << endl;
                eastDigits = true;
              }
              crdhi.Next();
              
            } //end loop over digits
            
            crdhi.GetSet()->Slice();
            crdhi.ResetFirst();
            
            
            if( dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneView() == PlaneView::kU){
              if( !westDigits || !eastDigits){ fOneSidedUPlanes += 1; }
            }
            if( dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneView() == PlaneView::kV){
              if( !westDigits || !eastDigits){ fOneSidedVPlanes += 1; }
            }
            planeItr.Next();
            
          } //end loop over planes
          
          planeItr.ResetFirst();
          
        }//end valid plane failure
        else if( status->GetVertexPlaneFailure() ){ 
          if( fNumberOfPlanesInEvent > 3 ){
            //MSG("DmxVertex", Msg::kDebug) << "Event = " << status->GetEventNumber() << endl;
            
            bool westSide = false;
            bool eastSide = false;
            
            while( planeItr.IsValid() ){
              
              crdhi.GetSet()->Slice(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber());
              
              while( crdhi.IsValid() ){
                
                if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest ){
                  //  MSG("DmxVertex", Msg::kDebug) << "\tPlane = " 
                  //                                         << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber() 
                  //                                         << "\tIsValid = "
                  //                                         << (Int_t)dynamic_cast<DmxPlane *>(planeItr.Ptr())->IsValid() 
                  //                                         << "\tSide = West\t" << "chip = " 
                  //                                         << crdhi.Ptr()->GetChannelId().GetVaChip()
                  //                                         << "\tpixel = " 
                  //                                         << status->GetDigitPixel(crdhi.Ptr()->GetChannelId().GetVaChannel()) 
                  //                                         << "\tcharge = " << crdhi.Ptr()->GetCharge() << endl;
                  westSide = true;
                }
                if(crdhi.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast ){
                  //  MSG("DmxVertex", Msg::kDebug) << "\tPlane = " 
                  //                                         << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber() 
                  //                                         << "\tIsValid = "
                  //                                         << (Int_t)dynamic_cast<DmxPlane *>(planeItr.Ptr())->IsValid() 
                  //                                         << "\tSide = East\t" << "chip = " 
                  //                                         << crdhi.Ptr()->GetChannelId().GetVaChip()
                  //                                         << "\tpixel = " 
                  //                                         << status->GetDigitPixel(crdhi.Ptr()->GetChannelId().GetVaChannel()) 
                  //                                         << "\tcharge = " << crdhi.Ptr()->GetCharge() << endl;
                  eastSide = true;
                }
                crdhi.Next();
                
              } //end loop over digits
              
              crdhi.GetSet()->Slice();
              crdhi.ResetFirst();
              
              planeItr.Next();
              
            } //end loop over planes
            
            planeItr.ResetFirst();
          } //if more than 10 planes
          
        }// end vertex plane failure
      }// end if EventLength test
      
                //Find Raw Data CandDeMuxDigitList fragment in PrimaryCandidateRecord.
                CandDeMuxDigitListHandle *crdlh1 = dynamic_cast<CandDeMuxDigitListHandle *>
                (candrec->FindCandHandle("CandDeMuxDigitListHandle", "canddigitlist"));
                
                // Check that crdlh exists.
                assert(crdlh1);
                
                //check that there are digits present
                //MSG("DeMuxCosmics", Msg::kWarning) << crdlh1->GetNDaughters() << " digits present at end of real data" << endl;
                
                
                
    } //end if FarDetector data
    
    else if( fDataType == 1 ){

      //MSG("Dmx", Msg::kDebug) << "event " << fEventNumber << endl;
      //*-* Find Raw Data CandDigitList fragment in PrimaryCandidateRecord.
      CandDigitListHandle *cmdlh = dynamic_cast<CandDigitListHandle *>
        (candrec->FindCandHandle("CandDigitListHandle", "candmcdigitlist"));
      
      // Check that cmdlh exists.
      assert(cmdlh);
      
      // Get CandDigitHandleItr over MC Truth digits.
      CandDigitHandleItr cmdhi(cmdlh->GetDaughterIterator());
       
      //create the KeyFunc
      CandDigitHandleKeyFunc *mcplaneKF = cmdhi.CreateKeyFunc();
      
      //program the KeyFunc with the sort function
      mcplaneKF->SetFun(KeyFromDigitPlane);
            
      //get the NavSet from the iterator and pass the KeyFunc to it
      cmdhi.GetSet()->AdoptSortKeyFunc(mcplaneKF);
      
      //clear the KF pointer because we no longer own the KeyFunc
      mcplaneKF = 0;
       
      //reset the mc and raw digit iterators
      cmdhi.ResetFirst();
      crdhi.ResetFirst();
      
      //reset the plane iterator
      planeItr.ResetFirst();
      

      //make arrays to hold the strip hit information
      Float_t reconStripsW[192], reconStripsE[192], mcStripsW[192], mcStripsE[192];
      Float_t chargeCorrect = 0., stripsCorrect = 0.;
      Float_t totalCharge = 0., totalStrips = 0.;
      Float_t totalCharge1Hyp = 0.;//, chargeCorrect1Hyp = 0.;
      Float_t mcCoG = 0.;
      Float_t mcNum = 0.;
      Float_t mcDeNom = 0.;
      Int_t firstStrip = 0;
      Int_t lastStrip = 0;
      fMCStripSpread24Planes = 0;
      fMCShowerSpread = 0;
      Int_t lowShowerStrip = 192;
      Int_t highShowerStrip = -1;
      //loop over the planes in the event
      while( planeItr.IsValid() ){
        DmxPlane *currentPlane = dynamic_cast<DmxPlane *>(planeItr.Ptr());
        fPlane = currentPlane->GetPlaneNumber();
        mcNum = 0.;
        mcDeNom = 0.;
        mcCoG = 0.;
        firstStrip = 0;
        lastStrip = 0;

        //MSG("Dmx", Msg::kDebug) << "event " << fEventNumber << "plane " << fPlane << endl;
        //-------------------get digit info----------------------------------//
        //zero the arrays
        for(Int_t i = 0; i < 192; i++){
          reconStripsW[i] = 0;
          mcStripsW[i] = 0;
          reconStripsE[i] = 0;
          mcStripsE[i] = 0;
        }
      
        //slice for the digits in this plane
        crdhi.GetSet()->Slice(fPlane);
        cmdhi.GetSet()->Slice(fPlane);
        
        //loop over the digits to fill the demuxed digit info
        while( crdhi.IsValid() ){
          CandDeMuxDigitHandle *digit = crdhi.Ptr();
          //get the strip, charge, crate for the digit
          if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
            reconStripsW[digit->GetPlexSEIdAltL().GetBestSEId().GetStrip()] += digit->GetCharge();
          }
          else if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
            reconStripsE[digit->GetPlexSEIdAltL().GetBestSEId().GetStrip()] += digit->GetCharge();
          }
          
          crdhi.Next();
        }

        //now for the mc digits
        while( cmdhi.IsValid() ){
          CandDigitHandle *digit = cmdhi.Ptr();
          //get the strip, charge, crate for the digit
          if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
            mcStripsW[digit->GetPlexSEIdAltL().GetBestSEId().GetStrip()] += digit->GetCharge();
          }
          else if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
            mcStripsE[digit->GetPlexSEIdAltL().GetBestSEId().GetStrip()] += digit->GetCharge();
          }
          mcNum += digit->GetCharge() * digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
          mcDeNom += digit->GetCharge();
          
          cmdhi.Next();
        }

        if(mcDeNom > 0.) mcCoG = mcNum/mcDeNom;

        crdhi.GetSet()->Slice();
        crdhi.ResetFirst();
        cmdhi.GetSet()->Slice();
        cmdhi.ResetFirst();
        
        for(Int_t ss = 0; ss<192; ss++){
          if(mcStripsE[ss]>0. && firstStrip == 0) firstStrip = ss;
          else if(mcStripsW[ss]>0. && firstStrip == 0) firstStrip = ss;
          if(mcStripsE[ss]>0. || mcStripsW[ss]>0.) lastStrip = ss;

          //fill variables to check and see if the shower extends more than a meter in the detector
          if(currentPlane->GetPlaneType()==DmxPlaneTypes::kShower && (mcStripsE[ss]>0.||mcStripsW[ss]>0.)){
            if(ss>highShowerStrip) highShowerStrip = ss;
            if(ss<lowShowerStrip) lowShowerStrip = ss;
          }
        }

        fMCStripSpread = lastStrip - firstStrip;

        for(Int_t s = 0; s<192; s++){
          fDeMuxedStrip = s;
          fDigitCharge = -1.;
          fMCDigitChargeE = -1.;
          fDigitChargeW = -1.;
          fMCDigitChargeW = -1.;
          fDigitSide = -1;

          if(mcStripsE[s] != 0.){
            fMCDigitChargeE = mcStripsE[s];
            totalCharge += mcStripsE[s];
            totalStrips += 1.;
            if(TMath::Abs(mcCoG-s)<=23) totalCharge1Hyp += mcStripsE[s];
          }
          if(mcStripsW[s] != 0.){
            fMCDigitChargeW = mcStripsW[s];
            totalCharge += mcStripsW[s];
            totalStrips += 1.;
            if(TMath::Abs(mcCoG-s)<=23) totalCharge1Hyp += mcStripsW[s];
          }
          if(reconStripsE[s] != 0.) fDigitCharge = reconStripsE[s];
          if(reconStripsW[s] != 0.) fDigitChargeW = reconStripsW[s];
            

          if(reconStripsW[s] == mcStripsW[s] && reconStripsW[s] != 0.){
            chargeCorrect += mcStripsW[s];
            //if(TMath::Abs(mcCoG-s)<=23) chargeCorrect1Hyp += reconStripsW[s];
          }
          if(reconStripsE[s] == mcStripsE[s] && reconStripsE[s] != 0.){
            chargeCorrect += mcStripsE[s];
            //if(TMath::Abs(mcCoG-s)<=23) chargeCorrect1Hyp += reconStripsE[s];
          }
          if((reconStripsE[s] == mcStripsE[s]) && (reconStripsW[s] == mcStripsW[s])) stripsCorrect += 1.;
         
          
          fMCDeMuxDigitTree->Fill();
          if ( fEventNumber%1000 == 0 ) {
            TDirectory* savedir = gDirectory; // save current working directory
            fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
            fMCDeMuxDigitTree->AutoSave();  // save tree to file
            fDeMuxFile->SaveSelf();  // save file directory containing this tree
            savedir->cd();  // change back to working to directory
          }

        }//end loop to fill digit tree

        if(fMCStripSpread>23) ++fMCStripSpread24Planes;
             
        planeItr.Next();
      } //end loop over planes in event
      
      fMCChargeCorrect = -1.;
      fMCChargeCorrect1Hyp = -1.;
      fMCStripsCorrect = -1.;
     
      if(totalCharge != 0.){
        fMCChargeCorrect = chargeCorrect/totalCharge;
        fMCTotalCharge = totalCharge;
      }
      if(totalCharge1Hyp != 0.){
        fMCChargeCorrect1Hyp = chargeCorrect/totalCharge1Hyp;
        fMCTotalCharge1Hyp = totalCharge1Hyp;
      }
      if(totalStrips != 0.) fMCStripsCorrect = stripsCorrect/totalStrips;
      fMCShowerSpread = highShowerStrip-lowShowerStrip;

      //cout << "filling mc demux tree" << endl;
      fMCDeMuxTree->Fill();
      if ( fEventNumber%1000 == 0 ) {
        TDirectory* savedir = gDirectory; // save current working directory
        fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
        fMCDeMuxTree->AutoSave();  // save tree to file
        fDeMuxFile->SaveSelf();  // save file directory containing this tree
        savedir->cd();  // change back to working to directory
      }

    } //end if MonteCarlo data
    
    //fill the DeMux info Tree

    fDeMuxTree->Fill();
    if ( fEventNumber%1000 == 0 ) {
      TDirectory* savedir = gDirectory; // save current working directory
      fDeMuxFile->cd();  // move to file directory of tree, oddly SaveSelf doesn't work without it
      fDeMuxTree->AutoSave();  // save tree to file
      fDeMuxFile->SaveSelf();  // save file directory containing this tree
      savedir->cd();  // change back to working to directory
    }

  }
  //MSG("JobCDmx", Msg::kDebug) << "finished ana" << endl;
  
  // Find PrimaryCandidateRecord fragment in MOM.
  CandRecord *candrec = dynamic_cast<CandRecord *>
          (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if(candrec == 0){
          result.SetError().SetFailed();
          return result;
  }
  
  
  //Find Raw Data CandDeMuxDigitList fragment in PrimaryCandidateRecord.
  CandDeMuxDigitListHandle *crdlh = dynamic_cast<CandDeMuxDigitListHandle *>
      (candrec->FindCandHandle("CandDeMuxDigitListHandle", "canddigitlist"));
  
  // Check that crdlh exists.
  assert(crdlh);
  
  //check that there are digits present
  //MSG("DeMuxCosmics", Msg::kWarning) << crdlh->GetNDaughters() << " digits present at end of ana" << endl;
  
  return result;
}

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

void DmxDeMuxCosmicsModule::HandleCommand(JobCommand *command) 
{
  //
  //  Purpose:  Method to interpret module commands.
  //
  //  Arguments:
  //    command   in    Command to interpret.
  //
  //  Return:   n/a
  //
  //  Commands implemented:
  //    Set TestToMake  Set the test to make to obtain a figure of merit - fraction of
  //                    digits within n strips from a track, fraction of mated signal
  //                    DigitsAlongTrack, and MatedSignalFraction respectively
  //    Set DataType    Set the type of data to test - MonteCarlo or FarDetector
  //    Set MatedSignalLimit  Set the lower limit for fraction of signal coming from opposite
  //                          sides of the same strip in a plane.
  //    Set DeltaStripsFromTrack  Set the number of strips a digit can be off from a track and
  //                              still be considered along that track.
  //
  //  see http://beaker.astro.indiana.edu/brebel/demux_notes/how_to_demux.html for details
   MSG("JobC", Msg::kDebug) << "DeMuxCosmicsModule::HandleCommand" << endl;

   TString cmd = command->PopCmd();
   if(cmd == "Set"){   
     TString opt = command->PopOpt();
     if(opt == "DataType"){ fDataType = command->PopIntOpt();}
     else if(opt == "TestToMake"){ 
       TString testopt = command->PopOpt();
       if(testopt == "SingleModuleDeMux"){ fSingleModuleDeMux = true;}
       if(testopt == "EventLength"){ fEventLength = true;}
       if(testopt == "UnsetDigitCheck"){ fUnsetDigitCheck = true;}
     }
     else if(opt == "UseStrayPlanesCheck"){ fUseStrayPlaneCheck = true;}
     else if(opt == "DontUseCandDigitMasks"){ fDontUseCandDigitMasks = true;}
     else if(opt == "WriteHistos"){ fWriteHistos = true;}
     else if(opt == "DontMakeTrees"){ fMakeTrees = false;}
     else if(opt == "WindowSize"){ fPlanesInWindow = command->PopIntOpt();}
     else if(opt == "StrayDeltaStripLimit"){ fStrayDeltaStripLimit = command->PopIntOpt();}
     else if(opt == "StrayPlanesLimit"){ fStrayPlanesLimit = command->PopIntOpt();}
     else if(opt == "Filter"){ cout << fFilter << endl; fFilter = 1;}
     else if(opt == "HistoFileName"){ fHistoFileName = command->PopOpt();}
     else {
         MSG("JobC", Msg::kWarning)<< "DeMuxCosmicsModule: Unrecognized option " << opt << endl;
      }
   }
   else {
     MSG("JobC", Msg::kWarning)<< "DeMuxCosmicsModule: Unrecognized command " << cmd << endl;
   }
   return;
}

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

const Registry& DmxDeMuxCosmicsModule::DefaultConfig() const
{
  //  Commands implemented:
  //    Set TestToMake  Set the test to make to obtain a figure of merit - fraction of
  //                    digits within n strips from a track, fraction of mated signal
  //                    DigitsAlongTrack, and MatedSignalFraction respectively
  //    Set DataType    Set the type of data to test - MonteCarlo or FarDetector
  //    Set MatedSignalLimit  Set the lower limit for fraction of signal coming from opposite
  //                          sides of the same strip in a plane.
  //    Set DeltaStripsFromTrack  Set the number of strips a digit can be off from a track and
  //                              still be considered along that track.
  //
  //  see http://beaker.astro.indiana.edu/brebel/demux_notes/how_to_demux.html for details
   int itrue = 1;  // work around for lack of bool in registry
   int ifalse = 0; // work around for lack of bool in registry

   MSG("JobC", Msg::kDebug) << "DeMuxCosmicsModule::DefaultConfig" << endl;

   static Registry r;
 
   r.UnLockValues();

   r.Set("SingleModuleDeMux",ifalse);
   r.Set("EventLength",      ifalse);
   r.Set("UnsetDigitCheck",  ifalse);
   r.Set("DataType",         0);
   
   r.Set("WriteHistos",          itrue);
   r.Set("MakeTrees",            ifalse);
   r.Set("HistoFileName",        "demux_cosmics.root");
   r.Set("Filter",               ifalse);
   r.Set("UseCandDigitMasks",    itrue);
   r.Set("WindowSize",           6);
   r.Set("StrayDeltaStripLimit", 6);
   r.Set("StrayPlaneLimit",      4);
   r.Set("UseFigureOfMerit",     ifalse);
   r.Set("HoughSlopeRMS",        5.);
   r.Set("HoughInterceptRMS",    10.);
   r.Set("HoughPeak",            0.66);
   r.Set("XTalkFractionLimit",  0.1);
   r.Set("RatioMatedSignalForValid",  0.5);
   
   r.LockValues();
      
   return r;
}

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

void DmxDeMuxCosmicsModule::Config(const Registry& r)
{
  int         tmpb;  // a temp bool. See comment under DefaultConfig...
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string

  if (r.Get("SingleModuleDemux",   tmpb)) fSingleModuleDeMux     =  tmpb;
  if (r.Get("EventLength",         tmpb)) fEventLength           =  tmpb;
  if (r.Get("UnsetDigitCheck",     tmpb)) fUnsetDigitCheck       =  tmpb;
  if (r.Get("UseCandDigitMasks",   tmpb)) fDontUseCandDigitMasks = !tmpb;
  if (r.Get("WriteHistos",         tmpb)) fWriteHistos           =  tmpb;
  if (r.Get("MakeTrees",           tmpb)) fMakeTrees             =  tmpb;
  if (r.Get("HistoFileName",       tmps)) fHistoFileName         =  tmps;
  if (r.Get("DataType",            tmpi)) fDataType              =  tmpi;
  if (r.Get("WindowSize",          tmpi)) fPlanesInWindow        =  tmpi;
  if (r.Get("StrayDeltaStripLimit",tmpi)) fStrayDeltaStripLimit  =  tmpi;
  if (r.Get("StrayPlaneLimit",     tmpi)) fStrayPlanesLimit      =  tmpi;
  if (r.Get("UseFigureOfMerit",    tmpb)) fUseStrayPlaneCheck    =  tmpb;
  if (r.Get("Filter",              tmpb)) fFilter                =  tmpb;
  if (r.Get("HoughSlopeRMS",       tmpd)) fHoughSlope            =  tmpd;
  if (r.Get("HoughInterceptRMS",   tmpd)) fHoughIntercept        =  tmpd;
  if (r.Get("HoughPeak",           tmpd)) fHoughPeak             =  tmpd;
  if (r.Get("XTalkFractionLimit",  tmpd)) fSignalFractionLimit   =  tmpd;
  if (r.Get("RatioMatedSignalForValid",  tmpd)) fMatedSignalForValid   =  tmpd;
 
}

//----------------------------------------------------------------------------------
void DmxDeMuxCosmicsModule::Help() 
{
  MSG("JobC", Msg::kInfo) 
    << "DeMuxCosmicsModule::Help\n"
    <<"DmxDeMuxCosmicsModule is a module which demultiplexes events "
    <<"in the far detector." << endl
    << "see http://beaker.astro.indiana.edu/brebel/demux_notes/how_to_demux.html for details"
    << endl;
}

void DmxDeMuxCosmicsModule::EndJob() 
{
  //save the trees to a file
  if(fMakeTrees && fDeMuxFile){
    
//     if(fDataType==1){
//       fMCDeMuxTree->Write();
//       fMCDeMuxDigitTree->Write();
//     }
    fDeMuxFile->Write();
    fDeMuxFile->Close();
  }
  return;
}

void DmxDeMuxCosmicsModule::FindDigitsChiSq(Double_t &uChiSq, Double_t &vChiSq, CandDeMuxDigitHandleItr &digitItr, 
                                            Double_t uIntercept, Double_t uSlope, Double_t vIntercept, 
                                            Double_t vSlope)
{
  uChiSq = 0.;
  vChiSq = 0.;
  Double_t xi = 0.;
  Double_t weight = 0.;
  Double_t yi = 0.;
  Double_t yxi = 0.;
  Double_t quotient = 0.;
  Int_t uDigits = -1;
  Int_t vDigits = -1;
  digitItr.ResetFirst();
    
  while( digitItr.IsValid() ){

    CandDeMuxDigitHandle *digit = digitItr.Ptr();

    if( digitItr.GetSet()->GetMasks().GetMask(digit) ){
      xi = digit->GetPlexSEIdAltL().GetPlane();

      //make the weight the inverse of the charge - that means that digits with more
      //charge make a bigger contribution to the chi^2 and you are more interested 
      //in getting those right.
      weight = 1.; //10000. / digit->GetCharge();
      yi = 1. * digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
      
      if( digit->GetPlexSEIdAltL().GetPlaneView() == PlaneView::kU ){
        ++uDigits;
        yxi = uIntercept + (xi * uSlope);
        quotient = (yi - yxi) / weight;
      
        uChiSq += quotient * quotient;
      }
      else if( digit->GetPlexSEIdAltL().GetPlaneView() == PlaneView::kV ){
        ++vDigits;
        yxi = vIntercept + (xi * vSlope);
        quotient = (yi - yxi) / weight;
      
        vChiSq += quotient * quotient;
      }
    }
    digitItr.Next();
  }  

  digitItr.ResetFirst();
  
  if(vDigits > 0){ vChiSq /= 1.*vDigits; }
  if(uDigits > 0){ uChiSq /= 1.*uDigits; }
  return;
}

void DmxDeMuxCosmicsModule::FindPlanesChiSq(Double_t &uChiSq, Double_t &vChiSq, TObjectItr &planeItr, 
                                            Double_t uIntercept, Double_t uSlope, Double_t vIntercept, 
                                            Double_t vSlope)
{
  uChiSq = 0.;
  vChiSq = 0.;
  Double_t xi = 0.;
  Double_t weight = 0.;
  Double_t yi = 0.;
  Double_t yxi = 0.;
  Double_t quotient = 0.;
  Int_t uPlanes = -1;
  Int_t vPlanes = -1;
  planeItr.ResetFirst();
    
  while( planeItr.IsValid() ){

    DmxPlane *plane = dynamic_cast<DmxPlane *>(planeItr.Ptr());

    xi = plane->GetPlaneNumber();

    //weight each plane as 1.  this means that the chi^2 is the square of the 
    //difference between the fit and actual cog
    weight = 1.; 
    yi = plane->GetCoG();
      
    if( plane->GetPlaneView() == PlaneView::kU ){
        ++uPlanes;
        yxi = uIntercept + (xi * uSlope);
        quotient = (yi - yxi) / weight;
      
        uChiSq += quotient * quotient;
    }
    else if( plane->GetPlaneView() == PlaneView::kV ){
        ++vPlanes;
        yxi = vIntercept + (xi * vSlope);
        quotient = (yi - yxi) / weight;
      
        vChiSq += quotient * quotient;
    }
    
    planeItr.Next();
  }  

  planeItr.ResetFirst();
  
  if(vPlanes > 0){ vChiSq /= 1.*vPlanes; }
  if(uPlanes > 0){ uChiSq /= 1.*uPlanes; }
  return;
}

---