ChopModule.cxx

Go to the documentation of this file.

// $Id: ChopModule.cxx,v 1.7 2007/11/11 08:26:13 rhatcher Exp $
//
// FILL_IN: [Document your code!!]
//
// n.tagg1@physics.ox.ac.uk
#include "ChopModule.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "Algorithm/AlgFactory.h"
#include "CandData/CandRecord.h"
#include "Candidate/CandContext.h"
#include "CandDigit/CandDigitListHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "CandChopListHandle.h"
#include "CandChopList.h"
#include "DigitVector.h"

#include "DataUtil/Truthifier.h"
#include "DataUtil/GetRunSnarlEvent.h"
#include "UgliGeometry/UgliGeomHandle.h"
//#include "DetSim/SimEventResult.h"
#include "Calibrator/Calibrator.h"


#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH2.h>
#include <TGraph.h>
#include <TStyle.h>

JOBMODULE(ChopModule, "ChopModule",
          "Chops digits into Chops (digitslices)");
CVSID("$Id: ChopModule.cxx,v 1.7 2007/11/11 08:26:13 rhatcher Exp $");



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

ChopModule::ChopModule() :
  fFile(0),
  fTree(0),
  fNuTree(0),
  fCanvas1(0),
  fCanvas2(0),
  fCanvas3(0),
  fleaf(0),
  fnu(0)
{
}

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

ChopModule::~ChopModule()
{
  // Save evaluation data (if any);
  if(fFile) {
    fTree->Write();
    fNuTree->Write();
    fFile->Write();
    fFile->Close();
  }

}

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

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

  // config.
  Registry& cfg = this->GetConfig();
  const char* algName       = cfg.GetCharString("ChopAlgorithm");
  const char* algConfigName = cfg.GetCharString("ChopAlgConfig");
  const char* listIn = cfg.GetCharString("ListIn");
  const char* listOut = cfg.GetCharString("ListOut");

  
  // Find PrimaryCandidateRecord fragment in MOM.
  CandRecord *candrec = dynamic_cast<CandRecord *>
    (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if (candrec == 0) {
    MSG("Chop", Msg::kWarning) << "No PrimaryCandidateRecord in MOM."
                               << endl;
    result.SetWarning().SetFailed();
    return result;
  }

  CandDigitListHandle *cdlh = dynamic_cast<CandDigitListHandle *>
    (candrec->FindCandHandle("CandDigitListHandle", listIn));

  // Require number of CandDigits to be non-zero.
  if (!cdlh || cdlh->GetNDaughters() < 1) {
    MSG("StripSR", Msg::kDebug)
      << "Null CandDigit list.  Bail out of event." << endl;
    result.SetFailed();
    return result;
  }
  
  AlgFactory &af = AlgFactory::GetInstance();
  AlgHandle adlh = af.GetAlgHandle(algName,algConfigName);
  CandContext cx(this, mom);
  cx.SetCandRecord(candrec);
  cx.SetDataIn(cdlh);

  CandChopListHandle chopListH = CandChopList::MakeCandidate(adlh, cx);
  chopListH.SetName(listOut);
  chopListH.SetTitle(Form("Chop list created by alg %s config %s",
                          algName, algConfigName) );
  candrec->SecureCandHandle(chopListH);
  
  
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}



const RawChannelId kQieRcid(Detector::kNear,ElecType::kQIE,0,0,false,false);

class Hi
{ public:
  int    index;
  float  time; 
  Hi(int i, float t) : index(i), time(t) {};
  bool operator<(const class Hi& other) const { return time < other.time; };
};

void TimeRankHistos(UInt_t n, std::vector<TH1*> inHists, std::vector<int>& outRank)
{
  outRank.resize(n);
  
  std::vector<Hi> sorter;
  for(UInt_t i=0;i<n;i++) {
    sorter.push_back(Hi(i,inHists[i]->GetMean()));
  }
  std::sort(sorter.begin(),sorter.end());

  for(UInt_t i=0;i<n;i++) {
    outRank[sorter[i].index] = i;
  }
}

//......................................................................
JobCResult ChopModule::Ana(const MomNavigator* mom)
{  
  JobCResult result;

  // Find PrimaryCandidateRecord fragment in MOM.
  CandRecord *candrec = dynamic_cast<CandRecord *>
    (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if (candrec == 0) {
    MSG("Chop", Msg::kWarning) << "No PrimaryCandidateRecord in  MOM."
                               << endl;
    result.SetWarning().SetFailed();
    return result;
  }  
  // Get CandDigits.
  CandDigitListHandle *cdlh_ptr = dynamic_cast<CandDigitListHandle *>
    (candrec->FindCandHandle("CandDigitListHandle", "canddigitlist"));
  if(cdlh_ptr==0) {
    MSG("Chop", Msg::kWarning) << "Couldn't find CandDigitListHandle." << endl;
    return JobCResult::kError;
  }

  // Get context. Bail out if this isn't MC.
  VldContext context = *(cdlh_ptr->GetVldContext());
  if(context.GetSimFlag() != SimFlag::kMC) {
    // MSG("Chop", Msg::kWarning) << "Running ChopModule::Ana is usless on non-MC data." << endl;
    //return JobCResult::kWarning;
  }

  // Get list of Chops.
  CandChopListHandle *choplist_ptr = dynamic_cast<CandChopListHandle *>
    (candrec->FindCandHandle("CandChopListHandle", "candchoplist"));
  if(choplist_ptr==0) {
    MSG("Chop", Msg::kWarning) << "Can't find ChopListHandle." << endl;
  }
  // Get list of SR slices.
  CandSliceListHandle *srslicelist_ptr = dynamic_cast<CandSliceListHandle *>
    (candrec->FindCandHandle("CandSliceListHandle", "CandSliceList"));
  //if(srslicelist_ptr==0) return JobCResult::kError;

  // Get simsnarl. Bail if we can't find it.
  bool haveTruth = false;
  SimSnarlRecord *ssr =
    dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));
  if (!ssr){
    MSG("Chop", Msg::kWarning) << "ChopModule::Ana couldn't find SimSnarl data, which it needs." << endl;
    haveTruth = false;
  } else {
    haveTruth = true;
  }

  // Get config
  double particleTimeout = GetConfig().GetDouble("particleTimeout");
  const char* filename   = GetConfig().GetCharString("filename");
  int         plots      = GetConfig().GetInt("plots");

  // Set up files, trees, etc
  TDirectory* oldDir = gDirectory;

  if(fFile==0) {
    fFile = new TFile(filename,"RECREATE");
    fFile->cd();

    fleaf = new ChopTreeLeaf;    
    fTree = new TTree("mytree","mytree");
    fTree->Branch("chopbr","ChopTreeLeaf",&fleaf);
    fTree->Write();
    
    fnu = new ChopNeutrinoLeaf; 
    fNuTree = new TTree("nutree","nutree");
    fNuTree->Branch("nubr","ChopNeutrinoLeaf",&fnu);
    fNuTree->Write();
    
    fFile->Write();
  }
  assert(fFile);
  fFile->cd();


  // Get context info: truth data, geometry, run number etc.

  // Find run number, etc:
  int run = -1, snarl = -1, event = -1;
  DataUtil::GetRunSnarlEvent(mom,run,snarl,event);
  // Get truth info:
  const Truthifier& truth = Truthifier::Instance(mom);

  // Get geometry:
  UgliGeomHandle ugli(context);
  // Find global time offset:
  //double timeOffset = 0;
  //const SimEventResult* ser = dynamic_cast<const SimEventResult*>
  //  (ssr->FindComponent("SimEventResult"));
  
  // Prep plots
  // Delete old histograms, if around.
  for(UInt_t ih = 0; ih<fCleanup.size(); ih++) delete fCleanup[ih];
  fCleanup.clear();
  
  // Make list of digits.
  DigitVector allDigits(cdlh_ptr);

  // Make list of chops.
  std::vector<DigitVector> chopList;
  if(choplist_ptr) {
    chopList.resize(choplist_ptr->GetNDaughters());
    
    CandDigitListHandleItr chopItr(choplist_ptr->GetDaughterIterator());
    UInt_t ichop=0;
    while( CandDigitListHandle* cdlh= chopItr() ) {
      chopList[ichop++].FillFrom(cdlh);
    }
    MSG("Chop",Msg::kDebug) << "Found " << chopList.size() << " Chops." << endl;
  }
  // Make list of SR slices.
  std::vector<DigitVector> refSliceList;

  if(srslicelist_ptr) {
    refSliceList.resize(srslicelist_ptr->GetNDaughters());
    CandSliceHandleItr cshItr(srslicelist_ptr->GetDaughterIterator());

    UInt_t ichop=0;
    while(CandSliceHandle* slice = cshItr()) {
      refSliceList[ichop++].FillFrom(slice);
    }
  } 
  MSG("Chop",Msg::kDebug) << "Found " << refSliceList.size() << " SR Slices." << endl;
  
  // Find earliest digit, last digit
    
  int tfirst = 90000000;
  int tlast  = -1;
  for(UInt_t idig = 0; idig<allDigits.size(); idig++) {
    CandDigitHandle cdh = allDigits[idig];
    int tdc = Calibrator::Instance().GetTDCFromTime(cdh.GetTime(CalTimeType::kNone), kQieRcid);    
    //int tdc2 = truth.GetRawDigit(cdh)->GetTDC();
    //if(tdc!=tdc2) MSG("Chop",Msg::kDebug) << "Round error: " << tdc << "/t" << tdc2 << endl;    
    if(tdc<tfirst) tfirst = tdc;
    if(tdc>tlast)  tlast  = tdc;
  }

  tfirst -= 5; // Pad with extra bins.
  tlast  += 5; // Ditto.

  // Neutrinos
  std::vector<Int_t> neutrinos(0);
  std::vector<TH1*>  nuHist(neutrinos.size());
  std::vector<TH1*>  nuHistInt(neutrinos.size());
  std::vector<float>  nuGraphTdc;
  std::vector<float>  nuGraphPlane;
  std::vector<float>  nuGraphE;
  std::vector<int>    nuGraphNu;
  std::vector<UInt_t> nuNGraph(neutrinos.size(),0); // How many points in each neutrino graph

  if(haveTruth) { // Have truth info

    // Make list of true neutrinos.  
    neutrinos = truth.GetListOfNeutrinoIndecies();

    // TZero times for neutrinos

    vector<double>  nuTZero(neutrinos.size(),1e19);  
  
    // Look at hits, find start times.
    for(Int_t s=0;s<truth.GetScintHitListSize(); s++) {
      const DigiScintHit* hit = truth.GetScintHit(s);
      if(hit){
        int nu_index = truth.GetNeutrinoOfTrack(hit->TrackId());
        for(UInt_t inu=0;inu<neutrinos.size();inu++) {
          if(neutrinos[inu]==nu_index){
            double t = hit->T1();
            if(t<nuTZero[inu]) nuTZero[inu]=t;
          }
        }
      }
    }


    ChopEvaluation allEval;
    allEval.Evaluate(allDigits,
                     truth, neutrinos,nuTZero, particleTimeout);
  
    // Create lists of neutrino digits.
    std::vector< std::vector<CandDigitHandle> >nuDigits(neutrinos.size());
    std::vector<ChopEvaluation> nuEval(neutrinos.size());
  
    for(UInt_t idig=0;idig<allDigits.size(); idig++) {
      for(UInt_t inu=0; inu < neutrinos.size(); inu++) {
        //float frac = truth.IsDigitFromNeutrino(neutrinos[inu],allDigits[idig]);
        float fract = truth.IsDigitFromNeutrino(neutrinos[inu],allDigits[idig],nuTZero[inu]+particleTimeout);
        if(fract>0) nuDigits[inu].push_back(allDigits[idig]);
      }
    }
    // evaluate each neutrino list.
    for(UInt_t inu=0; inu < neutrinos.size(); inu++) {
      nuEval[inu].Evaluate(nuDigits[inu],
                           truth, neutrinos,nuTZero, particleTimeout);
    }

    // Truth info for each neutrino:
    std::vector<ChopNeutrinoInfo> nuInfo(neutrinos.size());

    for(UInt_t inu=0; inu < neutrinos.size(); inu++) {
      const REROOT_NeuKin* neukin = truth.GetNeuKin(inu);
      const TParticle*     part   = truth.GetStdHepParticle(neutrinos[inu]);
      if(!neukin) continue;
      if(!part) continue;
      nuInfo[inu].inu = inu;
      nuInfo[inu].vtxt   = part->T();
      nuInfo[inu].vtxtdc = Calibrator::Instance().GetTDCFromTime(part->T(),kQieRcid);
      nuInfo[inu].vtxx   = part->Vx();
      nuInfo[inu].vtxy   = part->Vy();
      nuInfo[inu].vtxz   = part->Vz();
      nuInfo[inu].enu    = neukin-> P4Neu()[3];
      nuInfo[inu].emu    = neukin-> P4Mu1()[3] + neukin->P4Mu2()[3];
      nuInfo[inu].eem    = neukin-> P4El1()[3] + neukin->P4El2()[3];
      nuInfo[inu].eshw   = neukin-> P4Shw()[3];
      nuInfo[inu].evis   = nuInfo[inu].emu+nuInfo[inu].eem+nuInfo[inu].eshw;
      nuInfo[inu].iaction= neukin-> IAction();
      nuInfo[inu].X      = neukin->X();
      nuInfo[inu].Y      = neukin->Y();
      nuInfo[inu].tzero  = nuTZero[inu]; 
    }
  
    // Evaluate each chop.
    std::vector<ChopEvaluation> chopEval;

    if(GetConfig().GetInt("eval_sr") ==0) {
      // Evaluate the chopper.
      chopEval.resize(chopList.size());
      for(UInt_t ichop=0; ichop<chopEval.size(); ichop++) {
        chopEval[ichop].Evaluate(chopList[ichop],
                                 truth, neutrinos,nuTZero, particleTimeout);
      }
    } else {
      // Evaluate the SR Slicer.
      chopEval.resize(refSliceList.size());
      for(UInt_t ichop=0; ichop<chopEval.size(); ichop++) {
        chopEval[ichop].Evaluate(refSliceList[ichop],
                                 truth, neutrinos,nuTZero, particleTimeout);
      }
    }

    // Start evaluation of my slicing


    for(UInt_t ichop = 0; ichop < chopEval.size(); ichop++) {
      fleaf->Reset();
      fleaf->snarl    = snarl;
      fleaf->chop    = ichop;
      (*(ChopEvaluation*)fleaf) = (chopEval[ichop]);

    
      // Best neutrino info:
      fleaf->complete1 = 0;
      fleaf->viscal1 = 0;
      fleaf->viscalt1 = 0;
      fleaf->visfidt1 = 0;
      fleaf->strips1 = 0;
      if(chopEval[ichop].fBestNeutrino1>=0) {
        (*(ChopNeutrinoInfo*)fleaf) = nuInfo[chopEval[ichop].fBestNeutrino1];
        fleaf->complete1 = chopEval[ichop].fNuMatch_calt[chopEval[ichop].fBestNeutrino1] / 
          allEval.fNuMatch_calt[chopEval[ichop].fBestNeutrino1];
        fleaf->viscal1 =  allEval.fNuMatch_cal[chopEval[ichop].fBestNeutrino1];
        fleaf->viscalt1 = allEval.fNuMatch_calt[chopEval[ichop].fBestNeutrino1];
        fleaf->visfidt1 = allEval.fNuMatch_caltfid[chopEval[ichop].fBestNeutrino1];
        fleaf->strips1  = allEval.fNuMatch_strips[chopEval[ichop].fBestNeutrino1];
      }

      // Second best neutrino info:
      fleaf->complete2 = 0;
      fleaf->viscal2 = 0;
      fleaf->viscalt2 = 0;
      fleaf->visfidt2 = 0;
      fleaf->strips2 = 0;
      if(chopEval[ichop].fBestNeutrino2>=0) {
        (*(ChopNeutrinoInfo*)fleaf) = nuInfo[chopEval[ichop].fBestNeutrino2];
        fleaf->complete2 = chopEval[ichop].fNuMatch_calt[chopEval[ichop].fBestNeutrino2] / 
          allEval.fNuMatch_calt[chopEval[ichop].fBestNeutrino2];
        fleaf->viscal2 =  allEval.fNuMatch_cal[chopEval[ichop].fBestNeutrino2];
        fleaf->viscalt2 = allEval.fNuMatch_calt[chopEval[ichop].fBestNeutrino2];
        fleaf->visfidt2 = allEval.fNuMatch_caltfid[chopEval[ichop].fBestNeutrino2];
        fleaf->strips2  = allEval.fNuMatch_strips[chopEval[ichop].fBestNeutrino2];
      }
    
      fTree->Fill();
    }


  
    // Start evaluation of neutrinos

    for(UInt_t inu=0; inu < neutrinos.size(); inu++) {
     
      // Find best and second-best chop for this neutrino

      int neutrinoToChop1 = -1;
      int neutrinoToChop2 = -1;
      float neutrinoComplete1 = 0;
      float neutrinoComplete2 = 0;
    
      for(UInt_t ichop = 0; ichop<chopEval.size(); ichop++) {
        float complete = chopEval[ichop].fNuMatch_cal[inu]/allEval.fNuMatch_cal[inu];

        if(complete>neutrinoComplete1) {
          neutrinoComplete1 = complete;
          neutrinoToChop1 = ichop;
        }
      }

      for(UInt_t ichop = 0; ichop<chopEval.size(); ichop++) {
        float complete = chopEval[ichop].fNuMatch_cal[inu]/allEval.fNuMatch_cal[inu];

        if((complete>neutrinoComplete2) && ((int)ichop!= neutrinoToChop1)) {
          neutrinoComplete2 = complete;
          neutrinoToChop2 = ichop;
        }
      }
      
      if(allEval.fNuMatch_cal[inu] > 0 ) {
        MSG("Chop",Msg::kDebug) << "Neutrino " << inu 
                                << "\tVisCal = " << allEval.fNuMatch_cal[inu]
                                << " \tTDC = " << nuEval[inu].fBigTdc
                                << " \tTrueStartTime = " << nuTZero[inu]*1e9 << endl;
        if(neutrinoToChop1<0) {
          MSG("Chop",Msg::kDebug) << "\t\t missed completely." << endl;
        } 
        else if(chopEval[neutrinoToChop1].fBestNeutrino1!=(int)inu) {
          MSG("Chop",Msg::kDebug) << "\t\t" 
                                  << " got absorbed ( completeness " << neutrinoComplete1 << ")"             
                                  << " by chop " << neutrinoToChop1
                                  << " time: " << chopEval[neutrinoToChop1].fBigTdc 
                                  << " which came from a neutrino " << chopEval[neutrinoToChop1].fBestNeutrino1
                                  << " viscalt: " << allEval.fNuMatch_calt[chopEval[neutrinoToChop1].fBestNeutrino1]
                                  << endl;
        
        } else {
          MSG("Chop",Msg::kDebug) << "\t\t" 
                                  << "matches chop " << neutrinoToChop1 
                                  << " completeness = " << neutrinoComplete1 
                                  << " Strips: nu/chop/match = " 
                                  << allEval.fNuMatch_strips[inu] << "/" 
                                  << chopEval[neutrinoToChop1].fStrips << "/"  
                                  << chopEval[neutrinoToChop1].fNuMatch_strips[inu]  
                                  << endl;
        }
        //for(UInt_t ichop = 0; ichop<chopEval.size(); ichop++) {
        //      float complete = chopEval[ichop].fNuMatch_calt[inu]/allEval.fNuMatch_calt[inu];
        //      MSG("Chop",Msg::kDebug) << "\t\tchop " << ichop << "\t" << complete << endl;
        //}
      }
    

      fnu->Reset();
      *((ChopNeutrinoInfo*)fnu) = nuInfo[inu];
      *((ChopEvaluation*)fnu)   = nuEval[inu];

      fnu->snarl = snarl;
      fnu->viscal = allEval.fNuMatch_cal[inu];
      fnu->viscalt= allEval.fNuMatch_calt[inu];
      fnu->visfidt= allEval.fNuMatch_caltfid[inu];
      fnu->visfid = allEval.fNuMatch_fid[inu];
      fnu->sigcorr= nuEval[inu].fSigcor;
      fnu->strips = allEval.fNuMatch_strips[inu];
      fnu->calstrips = allEval.fNuMatch_calstrips[inu];

      fnu->chop1 = neutrinoToChop1;
      if(neutrinoToChop1>=0) {
        fnu->complete1      = neutrinoComplete1;
        fnu->completet1     = chopEval[neutrinoToChop1].fNuMatch_calt[inu]/allEval.fNuMatch_calt[inu];
        fnu->complete_tot1  = chopEval[neutrinoToChop1].fNuMatch[inu]/allEval.fNuMatch[inu];
        fnu->strips1        = chopEval[neutrinoToChop1].fStrips;
        fnu->calstrips1     = chopEval[neutrinoToChop1].fCalStrips;
        fnu->stripsmatch1   = chopEval[neutrinoToChop1].fNuMatch_strips[inu];
        fnu->calstripsmatch1= chopEval[neutrinoToChop1].fNuMatch_calstrips[inu];
        fnu->chop1_sigcorr= chopEval[neutrinoToChop1].fSigcor;
        fnu->chop1_viscal = chopEval[neutrinoToChop1].fSigcor_cal;
        fnu->chop1_visfid = chopEval[neutrinoToChop1].fSigcor_fid;
        fnu->chop1_nu     = chopEval[neutrinoToChop1].fBestNeutrino1;
        fnu->chop1_tdc    = chopEval[neutrinoToChop1].fBigTdc;
        fnu->chop1_planebeg= chopEval[neutrinoToChop1].fPlaneBeg;
        fnu->chop1_planeend= chopEval[neutrinoToChop1].fPlaneEnd;
      }
    
      fnu->chop2 = neutrinoToChop2;
      if(neutrinoToChop2>=0) {
        fnu->complete2      = neutrinoComplete2;
        fnu->completet2     = chopEval[neutrinoToChop2].fNuMatch_calt[inu]/allEval.fNuMatch_calt[inu];
        fnu->complete_tot2  = chopEval[neutrinoToChop2].fNuMatch[inu]/allEval.fNuMatch[inu];
        fnu->strips2      = chopEval[neutrinoToChop2].fStrips;
        fnu->calstrips2      = chopEval[neutrinoToChop2].fCalStrips;
        fnu->stripsmatch2 = chopEval[neutrinoToChop2].fNuMatch_strips[inu];
        fnu->calstripsmatch2 = chopEval[neutrinoToChop2].fNuMatch_calstrips[inu];
        fnu->chop2_sigcorr= chopEval[neutrinoToChop2].fSigcor;
        fnu->chop2_viscal = chopEval[neutrinoToChop2].fSigcor_cal;
        fnu->chop2_visfid = chopEval[neutrinoToChop2].fSigcor_fid;
        fnu->chop2_nu     = chopEval[neutrinoToChop2].fBestNeutrino1;
        fnu->chop2_tdc    = chopEval[neutrinoToChop2].fBigTdc;
        fnu->chop2_planebeg= chopEval[neutrinoToChop2].fPlaneBeg;
        fnu->chop2_planeend= chopEval[neutrinoToChop2].fPlaneEnd;
      }


      fNuTree->Fill();
    }

    if(plots) {
      nuHist.resize(neutrinos.size(),0);
      nuHistInt.resize(neutrinos.size(),0);
      nuNGraph.resize(neutrinos.size(),0);
      for(UInt_t inu = 0;inu < neutrinos.size(); inu++) {
        nuHist[inu] = new TH1F(Form("nu_%d",inu),
                               Form("Hits Nu:%d ",inu ),
                               600, 0, 600);    
        nuHistInt[inu] = new TH1F(Form("nuint_%d ",inu ),
                                  Form("Hits Nu:%d ",inu),
                                  600, 0, 600);        
        fCleanup.push_back(nuHist[inu]);
        fCleanup.push_back(nuHistInt[inu]);
      }
    
      // 
      // run through digits for plotting
      //
    
      for(UInt_t idig = 0; idig<allDigits.size(); idig++) {
        CandDigitHandle cdh = allDigits[idig];
        int tdc =  Calibrator::Instance().GetTDCFromTime(cdh.GetTime(CalTimeType::kNone), kQieRcid);    
        float sigcor =  cdh.GetCharge(CalDigitType::kSigCorr);
        int t = (tdc-tfirst);
        int plane = cdh.GetPlexSEIdAltL().GetPlane();
        PlexStripEndId seid = cdh.GetPlexSEIdAltL().GetBestSEId();
        UgliStripHandle ustrip = ugli.GetStripHandle(seid);
   
        if(plane<121) {

          for(UInt_t inu=0;inu<neutrinos.size();inu++) {
            int neu = neutrinos[inu];
            float frac  = truth.IsDigitFromNeutrino(neu,cdh);
            float fract = truth.IsDigitFromNeutrino(neu,cdh,nuTZero[inu]+particleTimeout);
            if(frac>0) {         
              nuHist[inu]->Fill(t,fract*sigcor);

              /*
            // find timing offset.
            // Reduced time: plane*0.0594m/c * 53.1e6
            const double c_buckets_per_plane = 0.0594 * 53.1e6/3e8;

            double tred = t - plane*c_buckets_per_plane;
            // Take off clear fibre length.
            tred -= ustrip.ClearFiber(seid.GetEnd()) * 1.6*53.1e6/3e8;
            
            // Take off green fibre length.
            const double rsqrt2 = 1.0/sqrt(2.0);
            const TVector3 kU(rsqrt2,rsqrt2,0.);
            const TVector3 kV(-rsqrt2,rsqrt2,0.);
            const TVector3 kBeamSpot(1.4885, 0.1397,0);
            Float_t stripX = 0.;
            if (seid.GetPlaneView()==PlaneView::kU) {
            float tpos = kBeamSpot.Dot(kV);
            float lambda = kV.Dot(ustrip.GlobalPos(0.));            
            stripX = -tpos + lambda;
            } else {        
            float tpos = kBeamSpot.Dot(kU);
            float lambda = kU.Dot(ustrip.GlobalPos(0.));
            stripX = tpos - lambda;
            }       
            double halfgreen = ustrip.GetHalfLength() + ustrip.WlsPigtail(seid.GetEnd());
            double wlsLen = halfgreen;
            wlsLen -= stripX;
            
            tred -= wlsLen * 1.6*53.1e6/3e8;
            
            nuHist[inu]->Fill(tred+0.0,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.1,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.2,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.3,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.4,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.5,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.6,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.7,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.8,fract*sigcor/10.);
            nuHist[inu]->Fill(tred+0.9,fract*sigcor/10.);
              */
              //nuHist[inu]->Fill(t,fract*sigcor);

              //nuHist[inu+neutrinos.size()]->Fill(t,fract*sigcor);
              nuGraphTdc.push_back(t);
              nuGraphPlane.push_back(plane);
              nuGraphNu.push_back(inu);
              nuGraphE.push_back(frac*sigcor);      
              nuNGraph[inu]++;
            }
          }
        }
      }
    }
  } // End haveTruth

  // Draw plots.


  if(plots) {
    //
    // Set up plots and data lists
    //
    
    fFile->cd();
    fFile->mkdir(Form("snarl%d",snarl));
    fFile->cd(Form("snarl%d",snarl));

    
    TH1* hitHist = 0;
    TH2* hit2Hist = 0;
    TH2* hitUHist = 0;
    TH2* hitVHist = 0;

    hitHist = new TH1F(Form("hits"),
                       Form("Digits Snarl %d",snarl),
                       600, 0, 600);
    fCleanup.push_back(hitHist);
    
    hit2Hist = new TH2F(Form("hits2"),
                        Form("Hits Vs Plane"),
                        600,0, 600,
                        61, 0, 121);
    fCleanup.push_back(hit2Hist);
    
    hitUHist = new TH2F(Form("hitsu"),
                        Form("Hits Vs U TPos"),
                        600,0, 600,
                        100, -2, 2);
    fCleanup.push_back(hitUHist);
    
    hitVHist = new TH2F(Form("hitsv"),
                        Form("Hits Vs V TPos"),
                        600,0, 600,
                        100, -2, 2);
    fCleanup.push_back(hitVHist);
    
    for(UInt_t idig = 0; idig<allDigits.size(); idig++) {
      CandDigitHandle cdh = allDigits[idig];
      int tdc =  Calibrator::Instance().GetTDCFromTime(cdh.GetTime(CalTimeType::kNone), kQieRcid);    
      float sigcor =  cdh.GetCharge(CalDigitType::kSigCorr);
      int t = (tdc-tfirst);
      int plane = cdh.GetPlexSEIdAltL().GetPlane();
      PlaneView::PlaneView_t view = cdh.GetPlexSEIdAltL().GetBestSEId().GetPlaneView();
      PlexStripEndId seid = cdh.GetPlexSEIdAltL().GetBestSEId();
      UgliStripHandle ustrip = ugli.GetStripHandle(seid);
      float tpos = ustrip.GetTPos();

      if(plane<121) {
        hitHist->Fill (t, sigcor);
        hit2Hist->Fill(t, plane, sigcor);

        if(view==PlaneView::kU) hitUHist->Fill(t, tpos, sigcor);
        else                    hitVHist->Fill(t, tpos, sigcor);
      }
    }
  

    // Chops
    std::vector<TH1*>   chopHist(chopList.size());
    std::vector<TH1*>   chopHistInt(chopList.size());
    std::vector<float>  chopGraphTdc;
    std::vector<float>  chopGraphPlane;
    std::vector<float>  chopGraphE;
    std::vector<int>    chopGraphChop;
    std::vector<UInt_t> chopNGraph(chopList.size(),0);

    for(UInt_t ichop = 0;ichop < chopList.size(); ichop++) {
      chopHist[ichop] = new TH1F(Form("chp_%d",ichop),
                                 Form("Hits Chop:%d",ichop),
                                 600, 0, 600);    
      chopHistInt[ichop] = new TH1F(Form("chpint_%d",ichop),
                                    Form("Hits Chop:%d ",ichop),
                                    600, 0, 600);        
      fCleanup.push_back(chopHist[ichop]);
      fCleanup.push_back(chopHistInt[ichop]);

      for(UInt_t idig = 0; idig < chopList[ichop].size(); idig++) {
        CandDigitHandle cdh = chopList[ichop][idig];
        int tdc =  Calibrator::Instance().GetTDCFromTime(cdh.GetTime(CalTimeType::kNone), kQieRcid);    
        float sig =  cdh.GetCharge(CalDigitType::kSigCorr);
        int t = (tdc-tfirst);
        int plane = cdh.GetPlexSEIdAltL().GetPlane();
        if(plane<121) {
          chopHist[ichop]->Fill(t,sig);
         
          chopGraphTdc.push_back(t);
          chopGraphPlane.push_back(plane);
          chopGraphE.push_back(sig);
          chopGraphChop.push_back(ichop);
          chopNGraph[ichop]++;
        }
      }
    }
    
    // Referece slices
    std::vector<TH1*>   refSliceHist(refSliceList.size());
    std::vector<TH1*>   refSliceHistInt(refSliceList.size());
    std::vector<float>  refSliceGraphTdc;
    std::vector<float>  refSliceGraphPlane;
    std::vector<float>  refSliceGraphE;
    std::vector<int>    refSliceGraphSlice;
    std::vector<UInt_t> refSliceNGraph(refSliceList.size(),0);

    for(UInt_t ichop = 0;ichop < refSliceList.size(); ichop++) {
      refSliceHist[ichop] = new TH1F(Form("refchp_%d",ichop),
                                     Form("Hits RefSlice:%d",ichop),
                                     600, 0, 600);    
      refSliceHistInt[ichop] = new TH1F(Form("refchpint_%d",ichop),
                                        Form("Hits RefSlice:%d ",ichop),
                                        600, 0, 600);        
      fCleanup.push_back(refSliceHist[ichop]);
      fCleanup.push_back(refSliceHistInt[ichop]);

      for(UInt_t idig = 0; idig < refSliceList[ichop].size(); idig++) {
        CandDigitHandle cdh = refSliceList[ichop][idig];
        int tdc =  Calibrator::Instance().GetTDCFromTime(cdh.GetTime(CalTimeType::kNone), kQieRcid);    
        float sig =  cdh.GetCharge(CalDigitType::kSigCorr);
        int t = (tdc-tfirst);
        int plane = cdh.GetPlexSEIdAltL().GetPlane();
        if(plane<121) {
          refSliceHist[ichop]->Fill(t,sig);
         
          refSliceGraphTdc.push_back(t);
          refSliceGraphPlane.push_back(plane);
          refSliceGraphE.push_back(sig);
          refSliceGraphSlice.push_back(ichop);
          refSliceNGraph[ichop]++;
        }
      }
    }


    std::vector<int> chpRank;
    TimeRankHistos(chopHist.size(),chopHist, chpRank);
    for(UInt_t i = 0; i< chopHist.size(); i++) {
      for(UInt_t j=i; j<chopHist.size(); j++) {
        chopHistInt[i]->Add(chopHist[j]);
      }
      int color = (chpRank[i] % 7) + 1; // Rotate through colors.
      chopHist[i] -> SetLineColor(color);
      chopHist[i] -> SetFillColor(color);
      chopHistInt[i] -> SetLineColor(color);
      chopHistInt[i] -> SetFillColor(color);
    }
    
    std::vector<int> refChpRank;
    TimeRankHistos(refSliceHist.size(),refSliceHist, refChpRank);
    for(UInt_t i = 0; i< refSliceHist.size(); i++) {
      for(UInt_t j=i; j<refSliceHist.size(); j++) {
        refSliceHistInt[i]->Add(refSliceHist[j]);
      }
      int color = (refChpRank[i] % 7) + 1; // Rotate through colors.
      refSliceHist[i] -> SetLineColor(color);
      refSliceHist[i] -> SetFillColor(color);
      refSliceHistInt[i] -> SetLineColor(color);
      refSliceHistInt[i] -> SetFillColor(color);
    }


    std::vector<int> nuRank;
    TimeRankHistos(neutrinos.size(), nuHist, nuRank);
    for(UInt_t inu = 0; inu< neutrinos.size(); inu++) {

      for(UInt_t i=inu; i<nuHist.size(); i++) {
        nuHistInt[inu]->Add(nuHist[i]);
      }
      int color = (nuRank[inu] % 7) +1;
      nuHist[inu] -> SetLineColor(color);
      nuHist[inu] -> SetFillColor(color);
      nuHistInt[inu] -> SetLineColor(color);
      nuHistInt[inu] -> SetFillColor(color);      
    }

    // Make graphs.
    std::vector<TGraph*> chopGraphs;
    for(UInt_t ichop=0;ichop<chopList.size();ichop++) {
      int n = chopNGraph[ichop];
      if(n>0) {
        TGraph* g = new TGraph(n);
        fCleanup.push_back(g);
        chopGraphs.push_back(g);
        g->SetMarkerStyle(20);
        g->SetMarkerColor((chpRank[ichop] % 7) +1);
        
        Int_t i=0;
        for(UInt_t p=0;p<chopGraphTdc.size();p++) {
          if(chopGraphChop[p] == (int)ichop) {
            g->SetPoint(i,chopGraphTdc[p],chopGraphPlane[p]);
            i++;
          }
        }
      }
    }

    std::vector<TGraph*> refSliceGraphs;
    for(UInt_t ichop=0;ichop<refSliceList.size();ichop++) {
      int n = refSliceNGraph[ichop];
      if(n>0) {
        TGraph* g = new TGraph(n);
        fCleanup.push_back(g);
        refSliceGraphs.push_back(g);
        g->SetMarkerStyle(20);
        g->SetMarkerColor((refChpRank[ichop] % 7) +1);  
        Int_t i=0;
        for(UInt_t p=0;p<refSliceGraphTdc.size();p++) {
          if(refSliceGraphSlice[p] == (int)ichop) {
            g->SetPoint(i,refSliceGraphTdc[p],refSliceGraphPlane[p]);
            i++;
          }
        }
      }
    }

    std::vector<TGraph*> nuGraphs;
    for(UInt_t inu=0;inu<neutrinos.size();inu++) {
      int n = nuNGraph[inu];
      if(n>0) {
        TGraph* g = new TGraph(n);
        fCleanup.push_back(g);
        nuGraphs.push_back(g);
        g->SetMarkerStyle(20);
        g->SetMarkerColor((nuRank[inu] % 7) +1);        
        Int_t i=0;
        for(UInt_t p=0;p<nuGraphTdc.size();p++) {
          if(nuGraphNu[p] == (int)inu) {
            g->SetPoint(i,nuGraphTdc[p],nuGraphPlane[p]);
            i++;
          }
        }
      }
    }
    


    gStyle->SetOptStat(0);

    if(!fCanvas1) fCanvas1 = new TCanvas("TimeCanvas","Time Histograms");
    fCanvas1->cd();
    fCanvas1->Clear();
    fCanvas1->SetLogy(1);
    fCanvas1->Divide(1,3,0.001,0.001,kWhite);
    
    fCanvas1->cd(1);
    fCanvas1->GetPad(1)->SetLogy(1);
    hitHist->Draw();
    for(UInt_t i = 0; i < nuHistInt.size(); i++) {
      if(nuHistInt[i]->GetEntries()>0) nuHistInt[i]->Draw("same");
    }
 
    fCanvas1->cd(2);
    fCanvas1->GetPad(2)->SetLogy(1);
    hitHist->Draw();
    for(UInt_t i = 0;i <chopHistInt.size(); i++) {
      if(chopHistInt[i]->GetEntries()>0) chopHistInt[i]->Draw("same");
    }
        
    fCanvas1->cd(3);
    fCanvas1->GetPad(3)->SetLogy(1);
    hitHist->Draw();
    for(UInt_t i = 0;i < refSliceHistInt.size(); i++) {
      if(refSliceHistInt[i]->GetEntries()>0) refSliceHistInt[i]->Draw("same");
    }    

    
    if(!fCanvas2) fCanvas2 = new TCanvas("WholeCanvas","Whole Event");
    fCanvas2->Clear();
    fCanvas2->Divide(1,3,0.001,0.001,kWhite);
    fCanvas2->cd(1);
    hit2Hist->Draw("colz");
    fCanvas2->cd(2);
    hitUHist->Draw("colz");
    fCanvas2->cd(3);
    hitVHist->Draw("colz");


    if(!fCanvas3) fCanvas3 = new TCanvas("PosTimeCanvas","Position vs Time");
    fCanvas3->cd();
    fCanvas3->Clear();
    fCanvas3->Divide(1,3,0.001,0.001,kWhite);
    TH2* h2frame = new TH2F("h2frame","Plane vs Time",600,0,600,130,0,130);
    fCleanup.push_back(h2frame);

    
    fCanvas3->cd(1);
    h2frame->Draw();
    for(UInt_t i=0;i<nuGraphs.size();i++) {
      nuGraphs[i]->Draw("p");
    }

    fCanvas3->cd(2);
    h2frame->Draw();
    for(UInt_t i=0;i<chopGraphs.size();i++) {
      chopGraphs[i]->Draw("p");
    }

    fCanvas3->cd(3);
    h2frame->Draw();
    for(UInt_t i=0;i<refSliceGraphs.size();i++) {
      refSliceGraphs[i]->Draw("p");
    }



  }

  {
    fFile->cd();
    fTree->AutoSave("SaveSelf Overwrite"); // SaveSelf makes it useable. Overwrite saves time, but is dangerous.
    fNuTree->AutoSave("SaveSelf Overwrite");
  }
  
  if(oldDir) oldDir->cd(); // restore old directory.
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

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

const Registry& ChopModule::DefaultConfig() const
{
  //======================================================================
  // Supply the default configuration for the module
  //======================================================================
  static Registry r; // Default configuration for module
  
  // Set name of config
  std::string name = this->JobCModule::GetName();
  name += ".config.default";
  r.SetName(name.c_str());
  
  // Set values in configuration
  r.UnLockValues();
  r.Set("ChopAlgorithm","AlgChopListFar");
  r.Set("ChopAlgConfig","default");
  r.Set("ListIn",         "canddigitlist");
  r.Set("ListOut",        "candchoplist");
  
  r.Set("plots",1);
  r.Set("eval_sr",0);
  r.Set("filename","mychopper.root");
  r.Set("particleTimeout",100e-9);

  r.LockValues();

  return r;
}

---