MergeEvent.cxx

Go to the documentation of this file.

#include <cassert>

#include "MuonRemoval/MergeEvent.h"
#include "MuonRemoval/RawDigitInfo.h"
#include "MuonRemoval/SelectEvent.h"
#include "MuonRemoval/CandRmMuListHandle.h"
#include "MuonRemoval/CandRmMuHandle.h"
#include "MuonRemoval/AlgRmMu.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "Calibrator/Calibrator.h"

#include <CandData/CandHeader.h>
#include "CandData/CandRecord.h"
#include "CandDigit/CandDeMuxDigitListHandle.h"
#include <RecoBase/CandTrackListHandle.h>
#include <RecoBase/CandFitTrackHandle.h>
#include <RecoBase/CandTrackHandle.h>
#include <RecoBase/CandShowerHandle.h>
#include <RecoBase/CandShowerListHandle.h>
#include <RecoBase/CandStripHandle.h>
#include <RecoBase/CandStripListHandle.h>
#include <RecoBase/CandEventListHandle.h>
#include <RecoBase/CandEventHandle.h>

#include "RawData/RawRecord.h"
#include "RawData/RawDigit.h"
#include "RawData/RawHeader.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawChannelId.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "RawData/RawDaqHeaderBlock.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawVarcErrorInTfBlock.h"

#include "Algorithm/AlgConfig.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"

#include "Record/SimSnarlRecord.h"
#include "REROOT_Classes/REROOT_NeuKin.h"

#include "Candidate/CandContext.h"

#include "TFile.h"
#include "TTree.h"
#include "TObjArray.h"
#include "TParticle.h"

#include "Conventions/ReleaseType.h"
#include "DataUtil/EnergyCorrections.h"

#include <cmath>

JOBMODULE(MergeEvent, "MergeEvent","");
CVSID("$Id: MergeEvent.cxx,v 1.14 2009/02/15 23:04:54 boehm Exp $");
//......................................................................
MergeEvent::MergeEvent() :
  fInputFileName(""), fInputFile(NULL), fInputTree(NULL),
  fInputEventNumber(0), fInRawDigits(NULL),
  fOutputFileName(""), fOutputFile(NULL), fOutputTree(NULL),
  fOutputEventNumber(0), fOutRawDigits(NULL),
  fMergeAlg(""), fMergeConfig(""), fMergeListOut(""), fUseTrkForTiming(1)
{
  fInputP4[0] = 0; fInputP4[1] = 0; fInputP4[2] = 0; fInputP4[3] = 0;  
  fOutputP4[0] = 0; fOutputP4[1] = 0; fOutputP4[2] = 0; fOutputP4[3] = 0;  
  for (int i = 0; i<6; i++){
    fInputMRM[i] = 0;
    fOutputMRM[i] = 0;
  }
}

void MergeEvent::EndJob()
{
  if(fOutputEventNumber>0) {
    fOutputFile->cd();
    fOutputTree->Write();
    fOutputFile->Close();
    MSG("EvtMrg",Msg::kInfo) << " Wrote digits from "<< fOutputEventNumber 
                             << " electron events to summary file: "
                             << fOutputFileName.c_str() << endl;
  }
  if(fInputEventNumber>0) {
    fInputFile->Close();
    MSG("EvtMrg",Msg::kInfo) << " Merged digits from "<< fInputEventNumber 
                             << " electron events from summary file: "
                             << fInputFileName.c_str() << endl;
  }
}

MergeEvent::~MergeEvent() 
{
  if(fInRawDigits) delete fInRawDigits;
  if(fOutRawDigits) delete fOutRawDigits;
}

//......................................................................
void MergeEvent::OutfileSetUp() 
{
  fOutputFile = new TFile(fOutputFileName.c_str(),"RECREATE");  
  fOutputTree = new TTree("ElectronTree","Electron Tree");
  fOutputTree->Branch("true_p4",&fOutputP4,"true_p4[4]/F",32000);
  fOutputTree->Branch("mrminfo",&fOutputMRM,"mrminfo[6]/F",32000);
  fOutRawDigits = new TClonesArray("RawDigitInfo",1000);
  fOutputTree->Branch("rawdigits","TClonesArray",&fOutRawDigits,8000,1);  
  MSG("EvtMrg",Msg::kInfo) << " Opened summary file "<< fOutputFileName
                           << " for writing electron events" << endl;
}

//......................................................................
void MergeEvent::InfileSetUp()
{
  TDirectory* current_directory = gDirectory;
  fInputFile = TFile::Open(fInputFileName.c_str(),"READ");
  if(!fInputFile || !(fInputFile->IsOpen()) ){
    MSG("EvtMrg",Msg::kFatal) << " Unable to open: "<< fInputFileName << endl;
    return;
  }
  fInputTree = dynamic_cast<TTree*>(fInputFile->Get("ElectronTree"));
  if(!fInputTree){
    MSG("EvtMrg",Msg::kFatal) << " Unable to get tree : ElectronTree from  "
                              << fInputFileName << endl;
    return;
  }
  fInRawDigits = new TClonesArray("RawDigitInfo",1000);
  fInputTree->SetBranchAddress("true_p4",fInputP4);
  if (fInputTree->GetBranchStatus("mrminfo")) fInputTree->SetBranchAddress("mrminfo",fInputMRM);
  fInputTree->SetBranchAddress("rawdigits",&fInRawDigits);
  current_directory->cd();
  MSG("EvtMrg",Msg::kDebug) << " Opened input file: " << fInputFileName <<endl;
}

//......................................................................
JobCResult MergeEvent::Ana(const MomNavigator* mom)
{
  //set up output tree if first time:
  if(!fOutputTree) this->OutfileSetUp();
  
  JobCResult result(JobCResult::kPassed);
  
  // make sure raw record is available
  RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
  if (rr == 0) {
    MSG("EventSR", Msg::kWarning) << "No RawRecord in MOM." << endl;
    return result;
  }

  const RawDaqSnarlHeader* snarlHdr =
    dynamic_cast<const RawDaqSnarlHeader*>(rr->GetRawHeader());
  if (snarlHdr) {
    if (snarlHdr->GetSnarl()%1000==0) cout<<"Run "<<snarlHdr->GetRun()<<" Snarl "<<snarlHdr->GetSnarl()<<endl;
  }
  else cout<<"No snarl info"<<endl;  

  //fill raw digit array
  fOutRawDigits->Clear();
  TClonesArray &rawdigitl = *fOutRawDigits;
  const RawDigitDataBlock *rddb=dynamic_cast<const RawDigitDataBlock *>
    (rr->FindRawBlock("RawDigitDataBlock"));
  Int_t rawDigitCounter = 0;
  if (rddb) {
    TIter it=rddb->GetDatumIter();
    while (TObject *obj=it.Next()) {
      RawDigit *rd=dynamic_cast<RawDigit *>(obj);
      new(rawdigitl[rawDigitCounter]) RawDigitInfo(rd->GetChannel(),
                                                   rd->GetADC(),
                                                   rd->GetTDC(),
                                                   rd->GetCrateT0());
      rawDigitCounter++;
    }
    rawdigitl.Compress();
  }

  //get truth:
  SimSnarlRecord* simrec = dynamic_cast<SimSnarlRecord*>
    (mom->GetFragment("SimSnarlRecord"));
  if ( !simrec ) {
    MSG("NtpMC",Msg::kWarning) << "No SimSnarlRecord in Mom" << endl;
    result.SetWarning().SetFailed();
    return result;
  }
  const TClonesArray* neukinarray = dynamic_cast<const TClonesArray*>
    (simrec->FindComponent("TClonesArray","NeuKinList"));
  REROOT_NeuKin* rneukin = dynamic_cast<REROOT_NeuKin*>(neukinarray->At(0));
  if (rneukin){
    for(int i=0;i<4;i++) fOutputP4[i] = rneukin->P4Shw()[i];
  }
  else {
    const TClonesArray* parts = 
      (const TClonesArray*)simrec->FindComponent("TClonesArray","StdHep");
    
    TIter next(parts);
    TObject* obj = 0;
    while ( ( obj = next() ) ) {
      TParticle* part = (TParticle*)obj;
      if (part->GetStatusCode()==1){
        fOutputMRM[0] = part->Px();
        fOutputMRM[1] = part->Py();
        fOutputMRM[2] = part->Pz();
        fOutputMRM[3] = part->Energy();
      }
      if (part->GetStatusCode()==800){
        fOutputP4[0] = part->Px();
        fOutputP4[1] = part->Py();
        fOutputP4[2] = part->Pz();
        fOutputP4[3] = part->Energy();
        fOutputMRM[4] = part->Vx();
        fOutputMRM[5] = part->Vy();
        fOutputTree->Fill();
        fOutputEventNumber++;
        return result;
      }
    }
  }
  //fill tree
  fOutputTree->Fill();
  fOutputEventNumber++;
  return result;
}

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

JobCResult MergeEvent::Reco(MomNavigator* mom)
{
  JobCResult result(JobCResult::kPassed);
  MSG("EvtMrg", Msg::kDebug) << " Starting RemoveMuon::Reco() " <<endl;
  MSG("EvtMrg", Msg::kDebug) << "  Alg        : " << fMergeAlg<< endl;
  MSG("EvtMrg", Msg::kDebug) << "  AlgConfig  : " << fMergeConfig<< endl;
  MSG("EvtMrg", Msg::kDebug) << "  NameListOut: " << fMergeListOut<< endl;

  //set up input tree if first time:
  if(!fInputTree) this->InfileSetUp();
  
  // Find PrimaryCandidateRecord fragment in MOM.
  CandRecord *record = dynamic_cast<CandRecord *>(mom->GetFragment("CandRecord"));
  if (record == 0) {
    MSG("EvtMrg", Msg::kError) << "No PrimaryCandidateRecord in MOM."<< endl;
    result.SetWarning().SetFailed();
    return result;
  }
  const CandHeader* header  = dynamic_cast<const CandHeader*>(record->GetHeader());  
  
  //Find the raw data in MOM
  RawRecord *rawrecord = 
    dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord",0,"DaqSnarl"));
  assert(rawrecord);
  RawDigitDataBlock *input_datablock =  NULL;
  if(rawrecord){
    TIter rdbit = rawrecord->GetRawBlockIter();
    TObject *tob;
    while ((tob = rdbit())) {
      if (tob->InheritsFrom("RawDigitDataBlock")) {
        input_datablock = (RawDigitDataBlock *) tob;
      }
    }
  }

  CandEventListHandle * eventlist = 
    dynamic_cast<CandEventListHandle*>(record->FindCandHandle("CandEventListHandle"));
  if(eventlist==NULL) { 
    MSG("EvtMrg",Msg::kWarning) << " Bailing out of Event eventlist = " << eventlist 
                                << endl; 
    result.SetWarning().SetFailed();
    return result;
  }

  //get muon removed digit list:
  CandDigitListHandle* digitlist = dynamic_cast<CandDigitListHandle*>
    (record->FindCandHandle("CandDigitListHandle","stripdigitlist"));

  //make a new TClonesArray to hold track time corrected 
  //RawDigitInfo objects for whole snarl
  TClonesArray *fRawDigits = new TClonesArray("RawDigitInfo",5000);

  //get rmmulist to check that input electrons match with removed muons
  CandRmMuListHandle *rmmulist = 
    dynamic_cast<CandRmMuListHandle*>(record->FindCandHandle("CandRmMuListHandle"));  

  //map to link MC electron RawDigits to RmMu objects
  std::map<Int_t,Int_t> linkRmMuRawDig;

  //Set up calibrator
  Calibrator& calibrator = Calibrator::Instance();
  calibrator.ReInitialise(*(record->GetVldContext()));

  Detector::Detector_t det = record->GetVldContext()->GetDetector();

  //loop over removed muons:
  Int_t rmmuIndex = 0;
  TIter rmmuItr(rmmulist->GetDaughterIterator());
  while(CandRmMuHandle* rmmu = dynamic_cast<CandRmMuHandle*>(rmmuItr())){
    //read in first electron:
    if(fInputTree->GetEntry(fInputEventNumber++)<=0){
      MSG("EvtMrg", Msg::kFatal) << "Unable to read input event number: " 
                                 << fInputEventNumber-1 <<endl;
      result.SetFatal();
      return result;
    }
    MSG("EvtMrg", Msg::kDebug)  << " Cand:" << header->GetRun() << "/" 
                                << header->GetSnarl() << "/" 
                                << rmmu->GetOrigEvtIndex() << endl;
    MSG("EvtMrg",Msg::kDebug) << " Read New event P4: " 
                              << fInputP4[0] <<", " 
                              << fInputP4[1] <<", "
                              << fInputP4[2] <<", "
                              << fInputP4[3] <<endl;

    //Check that inputs and rmmu values agree
    const Float_t emass  = 0.000511; //GeV
    const Float_t mumass = 0.105658; //GeV
    
    Float_t etot = 0;
    Float_t pelec = 0;

    double rangemom = rmmu->GetMomRange();
    double curvemom = rmmu->GetMomCurv();

    const VldContext vc = *(record->GetVldContext());
    ReleaseType::Release_t release = ReleaseType::kCedarPhyDaikon;

    EnergyCorrections::WhichCorrection_t corrver = EnergyCorrections::kDefault;
    if (rangemom>0) rangemom=EnergyCorrections::FullyCorrectMomentumFromRange(rangemom,vc,release,corrver);
    curvemom = EnergyCorrections::FullyCorrectSignedMomentumFromCurvature(curvemom,vc,release,corrver);


    if(rmmu->GetIsCont()==1){
      etot =  sqrt(rangemom*rangemom + pow(mumass,2));
      pelec  = sqrt(etot*etot - emass*emass);
    }
    else { //use range:
      etot = sqrt(pow(curvemom,2)+mumass*mumass);
      pelec  = sqrt(etot*etot - emass*emass);
    }

    if(etot > 100) {etot = emass; pelec = 0.0;}
    if(fInputP4[3] > 100) {
      fInputP4[0] = 0;
      fInputP4[1] = 0;
      fInputP4[2] = 0;
      fInputP4[3] = emass;
    } 
   float cut = 1e-3;

    if( TMath::Abs(rmmu->GetVtxDCosX()*pelec - fInputP4[0])> cut ||
        TMath::Abs(rmmu->GetVtxDCosY()*pelec - fInputP4[1])> cut ||
        TMath::Abs(rmmu->GetVtxDCosZ()*pelec - fInputP4[2])> cut ||
        TMath::Abs(etot - fInputP4[3])> cut ) {
      MSG("EvtMrg", Msg::kInfo) << " Input momentum: " << fInputP4[0] << " "
                                << fInputP4[1] << " " << fInputP4[2] << " " 
                                << fInputP4[3] << endl
                                << " RmMu momentum: " << rmmu->GetVtxDCosX()*pelec 
                                << " " << rmmu->GetVtxDCosY()*pelec << " " 
                                << rmmu->GetVtxDCosZ()*pelec << " " << etot << endl;
      MSG("EvtMrg", Msg::kError)
        << " Input electron file out of synch with muon removal..." 
        << " look for a match nearby." << endl;
      
      Bool_t gotMatch = false;
      for(int ii=1;ii<50;ii++){
        for(int jj=-1;jj<=1;jj+=2){       
          if(fInputTree->GetEntry(fInputEventNumber+(ii*jj))) {
            if( TMath::Abs(rmmu->GetVtxDCosX()*pelec - fInputP4[0])< cut &&
                TMath::Abs(rmmu->GetVtxDCosY()*pelec - fInputP4[1])< cut &&
                TMath::Abs(rmmu->GetVtxDCosZ()*pelec - fInputP4[2])< cut &&
                TMath::Abs(etot - fInputP4[3])<1e-4 ) {
              //got a match!
              fInputEventNumber+=(ii*jj)+1;
              gotMatch = true;
              MSG("EvtMrg", Msg::kInfo) 
                << "Found an input momentum match; "
                "ii = " << ii << " jj = " << jj << " "
                << rmmu->GetVtxDCosX()*pelec - fInputP4[0] << " " 
                << rmmu->GetVtxDCosY()*pelec - fInputP4[1] << " " 
                << rmmu->GetVtxDCosZ()*pelec - fInputP4[2] << " " 
                << etot - fInputP4[3] << endl;
              break;
            }
          }       
        }
        if(gotMatch) break;
      }

      if(!gotMatch) {
        MSG("EvtMrg", Msg::kError) << "Can't find a match... quitting" << endl;
        result.SetWarning().SetFailed();
        return result;
      }
    }

    //We have a rmmu entry and a set of electron hits that match
    //check if the MC electron/muon digit list is empty
    //This might not be a problem for electrons. But sometimes the muon event has no hits since the MC muon momentum is independent of the original muon momentum and can be very small. TY 11/16/07
    if (!fInRawDigits->GetEntries()){
      MSG("EvtMrg", Msg::kWarning) << "The lepton has no hits, don't add it to the muon removal event." << endl;
      continue;
    }
    
    //save MRM information: pmu, Q2 and Eshw
    rmmu->SetMRMX(fInputMRM[0]);
    rmmu->SetMRMY(fInputMRM[1]);
    rmmu->SetMRMZ(fInputMRM[2]);
    rmmu->SetMRMQ2(fInputMRM[4]);
    rmmu->SetMRMEshw(fInputMRM[5]);

    //now find original event and track to get the timing offsets
    TIter evtItr(eventlist->GetDaughterIterator());
    const CandEventHandle *event = NULL;
    for(int i=0;i<rmmu->GetOrigEvtIndex()+1;i++) 
      event = dynamic_cast<const CandEventHandle*>(evtItr());      
    if(!SelectEvent(event)) {
      MSG("EvtMrg", Msg::kError) << " This event was selected for muon removal " 
                                 << " Now it is not passing... quitting"
                                 << endl;
      result.SetWarning().SetFailed();
      return result;
    }
    const CandTrackHandle* track = GetRemovableTrack(event);
    if(!track) {
      MSG("EvtMrg", Msg::kError) << " Cannot find removable track from event" 
                                 << " ... quitting"
                                 << endl;
      result.SetWarning().SetFailed();
      return result;
    }
    //check that track matches rmmu:
    if(rmmu->GetVtxPlane()!=track->GetVtxPlane() ||
       rmmu->GetNPlane()!= TMath::Abs(track->GetEndPlane()-
                                      track->GetVtxPlane())+1) {
      MSG("EvtMrg", Msg::kError) << " Removable track properties do not match"
                                 << " rmmu properties... quitting"
                                 << endl;
      result.SetWarning().SetFailed();
      return result;
    }

    //now calculate mean offsets between track/event and MC electron:
    int tracktdc = 0;
    int tracktdc_ndigit = 0 ;
    int track_offset = -1;
    if(fUseTrkForTiming){
      MSG("EvtMrg", Msg::kDebug) << " Using Track for timing " << endl;
      TIter trkstpIter(track->GetDaughterIterator());
      while(const CandStripHandle* strip = 
            dynamic_cast<const CandStripHandle*>(trkstpIter())){
        TIter trkdigIter(strip->GetDaughterIterator());
        if(strip->GetPlane()<track->GetVtxPlane()+10){
          while(const CandDigitHandle* digit = 
                dynamic_cast<const CandDigitHandle*>(trkdigIter())){
            const RawDigit* raw_digit = input_datablock->At(digit->GetRawDigitIndex());
            if(raw_digit){
              if(track_offset==-1) track_offset = raw_digit->GetTDC();    
              tracktdc += raw_digit->GetTDC() - track_offset;
              tracktdc_ndigit++;
              MSG("EvtMrg", Msg::kDebug) << " TDC:  " << raw_digit->GetTDC() 
                                         << " ADC:  " << raw_digit->GetADC() << endl;
            }
            else MSG("EvtMrg", Msg::kError) << "Digit with no raw digit"<<endl;   
          }
        }
      }
    }
    else {
      //Use all digts for finding the mean time of the event. 
      //  NB: You want to be careful that you dont have outliers screwing it up so 
      //      use a 5.0PE cut rather than the usual 3PE cut when doing timing 
      //      along tracks. 
      TIter digIter(digitlist->GetDaughterIterator());
      while(const CandDigitHandle* digit = 
            dynamic_cast<const CandDigitHandle*>(digIter())){
        if(digit->GetCharge(CalDigitType::kPE)>5.0){
          const RawDigit* raw_digit = input_datablock->At(digit->GetRawDigitIndex());   
          if(raw_digit){
            if(track_offset==-1) track_offset = raw_digit->GetTDC();      
            tracktdc+=raw_digit->GetTDC() - track_offset;
            tracktdc_ndigit++;
            MSG("EvtMrg", Msg::kDebug) << " TDC:  " << raw_digit->GetTDC() 
                                       << " ADC:  " << raw_digit->GetADC() << endl;
          }
          else {
            MSG("EvtMrg", Msg::kError) << "Digit with no raw digit"<<endl;
          }
        }
      }    
    }

    int meantracktdc = (track_offset!=-1)?((int)(round(((float)tracktdc)/ 
                                                       tracktdc_ndigit))+track_offset):0;
    MSG("EvtMrg", Msg::kDebug) << " Mean Track TDC:  " << meantracktdc << endl;

    
    //Calc mean and mode timing offset for the "MC" events.
    int rawtdc = 0;
    int rawtdc_ndigit = 0;
    TIter rawdigitIter(fInRawDigits);
    RawDigitInfo *tmp_digit = 0;
    std::map<Int_t,Int_t> freqTDC;
    while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawdigitIter()))){
      rawtdc_ndigit++;
      rawtdc += (int)(tmp_digit->tdc);
      if(freqTDC.find(int(tmp_digit->tdc))!=freqTDC.end()) 
        freqTDC[int(tmp_digit->tdc)] += 1;
      else freqTDC[int(tmp_digit->tdc)] = 1;
    }
    int moderawtdc = 0;
    int maxrawtdc = 0;
    std::map<Int_t,Int_t>::iterator begTDC = freqTDC.begin();
    std::map<Int_t,Int_t>::iterator endTDC = freqTDC.end();
    while(begTDC!=endTDC){
      if(begTDC->second>maxrawtdc) {
        moderawtdc = begTDC->first;
        maxrawtdc = begTDC->second;
      }
      begTDC++;
    }
    int meanrawtdc = 0;
    if(rawtdc_ndigit>0) meanrawtdc = (int)(round(((float)rawtdc)/rawtdc_ndigit));
    MSG("EvtMrg", Msg::kDebug) << "MODE, MAX, MEAN: " << moderawtdc << " " 
                               << maxrawtdc << " " << meanrawtdc << endl;

    //remove large outliers which skew the mean:
    if(maxrawtdc>4) meanrawtdc = moderawtdc; //trust the mode
    rawtdc = rawtdc_ndigit = 0;
    rawdigitIter.Reset();
    while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawdigitIter()))){   
      if(TMath::Abs(tmp_digit->tdc-meanrawtdc)<=5) { //ND: 5 buckets ~ 100ns
        rawtdc_ndigit++;
        rawtdc += (int)(tmp_digit->tdc);
        MSG("EvtMrg", Msg::kDebug) << " TDC:  " << tmp_digit->tdc 
                                   << " ADC:  " << tmp_digit->adc << endl;
      }
    }
    if(rawtdc_ndigit>0) meanrawtdc = (int)(round(((float)rawtdc)/rawtdc_ndigit));
    MSG("EvtMrg", Msg::kDebug) << " Mean Electron TDC:  " << meanrawtdc << endl;

    MSG("EvtMrg", Msg::kDebug)<< "Delta T muon and electron events: " << meantracktdc
                              <<" - "<< meanrawtdc <<" = " 
                              << meantracktdc - meanrawtdc <<endl;     
    
    //Now add RawDigitInfo objects to new TClonesArray for the snarl
    //and correct for the average delta time between the muon and electron
    rawdigitIter.Reset();
    TClonesArray &rawdigitl = *fRawDigits;
    Int_t rawDigitCounter = fRawDigits->GetEntries();
    tmp_digit = 0;
    while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawdigitIter()))) {
      Double_t TDC = tmp_digit->tdc - meanrawtdc + meantracktdc;
      if(TDC<0) TDC = 0;
      linkRmMuRawDig[rawDigitCounter] = rmmuIndex;
      new(rawdigitl[rawDigitCounter++]) RawDigitInfo(tmp_digit->rcid,tmp_digit->adc,
                                                     TDC,tmp_digit->t0);
    }
    rmmuIndex += 1;
  }
  
  // Add back in the electrons
  MSG("EvtMrg", Msg::kDebug) << " Running Event Merging:  " << fMergeAlg 
                             << ":"<< fMergeConfig<< " to produce " 
                             << fMergeListOut <<endl;   
  
  CandContext merge_cx(NULL, mom);
  TObjArray merge_input;
  merge_input.Add(fRawDigits);
  merge_input.Add(digitlist);
  merge_input.Add(rawrecord);
  merge_cx.SetDataIn(&merge_input);
  merge_cx.SetCandRecord(record);
  AlgHandle merge_algorithm = 
    AlgFactory::GetInstance().GetAlgHandle(fMergeAlg.c_str(), fMergeConfig.c_str());
  CandDigitListHandle merge_digitlist = 
    CandDigitList::MakeCandidate(merge_algorithm, merge_cx);
  merge_digitlist.SetName(fMergeListOut.c_str());
  merge_digitlist.SetTitle(fMergeListOut.c_str());
  record->SecureCandHandle(merge_digitlist);

  //now need to replace daughters in rmmu objects with those in new digitlist
  CandDigitListHandle* mergedlist = dynamic_cast<CandDigitListHandle*>
    (record->FindCandHandle("CandDigitListHandle",fMergeListOut.c_str()));
  TIter mergedlistIter(mergedlist->GetDaughterIterator());

  rmmuItr.Reset();
  while(CandRmMuHandle* rmmu = dynamic_cast<CandRmMuHandle*>(rmmuItr())){
    std::vector<CandDigitHandle*> mdigs(rmmu->GetNDaughters());
    std::vector<CandDigitHandle*> rmdigs(rmmu->GetNDaughters());
    std::vector<Int_t> rfk(rmmu->GetNDaughters());
    Int_t rmmuDigCounter = 0;
    TIter rmmuDigIter(rmmu->GetDaughterIterator());
    while(CandDigitHandle *rmmuDig = dynamic_cast<CandDigitHandle*>(rmmuDigIter())){      
      rfk[rmmuDigCounter] = rmmu->ReasonForKeeping(rmmuDig);
      mergedlistIter.Reset();
      while(CandDigitHandle *mergeDig = dynamic_cast<CandDigitHandle*>(mergedlistIter())){
        if( ( mergeDig->GetRawDigitIndex()==rmmuDig->GetRawDigitIndex() ) ){
          mdigs[rmmuDigCounter] = mergeDig;
          rmdigs[rmmuDigCounter] = rmmuDig;
          break;
        }
      }
      rmmuDigCounter += 1;
    }
    for(int i=0;i<rmmu->GetNDaughters();i++) {
      if(rmdigs[i]){
        rmmu->RemoveDaughter(rmdigs[i]);
        rmmu->AddDaughterLink(*mdigs[i]);
        rmmu->SetReasonForKeeping(mdigs[i],rfk[i]);
      }
    }
  }

  //now add electron digits to RmMu objects, and set ReasonForKeeping IsMCElec flag
  //loop over new raw digits:
  TIter rawDigIter(fRawDigits);
  Int_t rawDigitCounter = 0;
  RawDigitInfo *tmp_digit = 0;
  
  while((tmp_digit = dynamic_cast<RawDigitInfo*>(rawDigIter()))) {
    Bool_t gotDigit = false;
    //get rmmu index
    Int_t rmmuInd = linkRmMuRawDig[rawDigitCounter];
    MSG("EvtMrg", Msg::kDebug) << "RawDigit: " << rawDigitCounter 
                               << " rmmu index: " << rmmuInd << endl;
    //now loop over merged digit list
    mergedlistIter.Reset();
    while(CandDigitHandle *mergeDig = dynamic_cast<CandDigitHandle*>(mergedlistIter())){
      //look for equivalent channels:
      if(tmp_digit->rcid == mergeDig->GetChannelId()) {
        MSG("EvtMrg", Msg::kVerbose) << "Got channel match for electron rawdigit " 
                                     << rawDigitCounter << endl;
        //get raw digit and check that TDC is the same
        const RawDigit* raw_digit = input_datablock->At(mergeDig->GetRawDigitIndex());
        //also check against digit info since electron-only hits have no associated RawDigit
        Double_t rdTime = calibrator.GetTimeFromTDC(int(tmp_digit->tdc), tmp_digit->rcid);

        Double_t rdAdc = tmp_digit->adc - 50.;   // 50 is a qie offset
        if(det == Detector::kFar)  rdAdc = tmp_digit->adc;

        if( ( raw_digit!=NULL && tmp_digit->tdc==raw_digit->GetTDC() ) || 
            ( TMath::Abs(rdTime - mergeDig->GetTime())<1e-9 && 
              TMath::Abs(rdAdc  - mergeDig->GetCharge())<1e-6 ) ) {
          //found match; 
          MSG("EvtMrg", Msg::kDebug) << "Got TDC match for electron rawdigit " 
                                     << rawDigitCounter << endl;
          MSG("EvtMrg", Msg::kDebug) << "Time: " << rdTime*1e6 << " " 
                                     << mergeDig->GetTime()*1e6
                                     << " ADC: " << rdAdc << " " << mergeDig->GetCharge()
                                     << endl;
          rmmuItr.Reset();
          Int_t rmmuIndex = 0;
          while(CandRmMuHandle* rmmu = dynamic_cast<CandRmMuHandle*>(rmmuItr())){
            if(rmmuIndex==rmmuInd) {
              MSG("EvtMrg", Msg::kDebug) << "Got RmMu match for electron rawdigit " 
                                         << rawDigitCounter << endl;
              //got rmmu object, now look for match in rmmu digit list
              TIter rmmuDigIter(rmmu->GetDaughterIterator());
              while(CandDigitHandle *rmmuDig = dynamic_cast<CandDigitHandle*>(rmmuDigIter())){
                Int_t rfk = rmmu->ReasonForKeeping(rmmuDig);
                if(rmmuDig->IsEquivalent(mergeDig)) {
                  MSG("EvtMrg", Msg::kDebug) << "Got RmMu Digit match for electron rawdigit " 
                                             << rawDigitCounter << endl;
                  //this digit is already in list, set IsMCElec flag
                  rfk |= RmMuMask::kRMMU_ISMCELEC_MASK;
                  rmmu->SetReasonForKeeping(rmmuDig,rfk);
                  gotDigit = true;
                  break;
                }
              }
              if(!gotDigit) {
                //then add the daughter link and set the flag
                MSG("EvtMrg", Msg::kDebug) << "Adding electron digit for electron rawdigit " 
                                           << rawDigitCounter << endl;
                Int_t rfk = 0;
                rfk |= RmMuMask::kRMMU_ISMCELEC_MASK;
                rmmu->AddDaughterLink(*mergeDig);
                rmmu->SetReasonForKeeping(mergeDig,rfk);
              }
            }
            if(gotDigit) break;
            rmmuIndex++;
          } // end of loop over rmmulist
        } // found a raw digit match
      } // found a channel ID match
      if(gotDigit) break;
    } // end of loop over new merged digitlist
    rawDigitCounter++;
  } // end of loop over MC electron raw digits

  fRawDigits->Clear();
  delete fRawDigits;
  
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

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

const Registry& MergeEvent::DefaultConfig() const
{
  // Default configuration for module
  static Registry r;
  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());
  // Set values in configuration
  r.UnLockValues();
  // Config for merging events
  r.Set("InputFileName","input_events.root");
  r.Set("OutputFileName","output_events.root");
  r.Set("MergeAlg","AlgMergeEvent");
  r.Set("MergeConfig","default");
  r.Set("MergeListOut","mergedigitlist");  
  r.Set("UseTrkForTiming",1);
  r.LockValues();
  return r;
}

//......................................................................
void MergeEvent::Config(const Registry& r)
{
  fInputFileName  = r.GetCharString("InputFileName");
  fOutputFileName = r.GetCharString("OutputFileName");
  fMergeAlg       = r.GetCharString("MergeAlg");
  fMergeConfig    = r.GetCharString("MergeConfig");
  fMergeListOut   = r.GetCharString("MergeListOut");
  fUseTrkForTiming   = r.GetInt("UseTrkForTiming");
  MSG("EvtMrg",Msg::kInfo) << " **** " << fInputFileName << endl;
  MSG("EvtMrg",Msg::kInfo) << " **** " << fOutputFileName << endl;
  MSG("EvtMrg",Msg::kInfo) << " **** " << fMergeAlg << endl;
  MSG("EvtMrg",Msg::kInfo) << " **** " << fMergeConfig << endl;
  MSG("EvtMrg",Msg::kInfo) << " **** " << fMergeListOut << endl;
  if(fUseTrkForTiming) MSG("EvtMrg",Msg::kInfo) << " **** Using Trk For Timing" << endl;
}

---