FitTrackMSListModule.cxx

Go to the documentation of this file.

// $Id: FitTrackMSListModule.cxx,v 1.11 2010/01/06 18:45:41 rhatcher Exp $
//
// A JobControl Module that calculates a fit and creates a 
// CandFitTrackMSList from a CandTrackList
//
// Tom Bringley
// ttb2@duke.edu
// 6/13/2001

#include <cassert>
#include <cmath>

#include "CandFitTrackMS/FitTrackMSListModule.h"

#include "Algorithm/AlgConfig.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "CandData/CandHeader.h"
#include "CandData/CandRecord.h"
#include "CandFitTrackMS/CandFitTrackMSListHandle.h"
#include "CandFitTrackMS/CandFitTrackMSList.h"
#include "CandFitTrackMS/CandFitTrackMSHandle.h"
#include "Candidate/CandContext.h"
#include "Conventions/Munits.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCommand.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "RawData/RawHeader.h"
#include "RawData/RawRecord.h"
#include "RecoBase/CandClusterListHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "Validity/VldContext.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "Calibrator/CalMapperFits.h"
#include "CandDigit/CandDigitListHandle.h"
#include "RawData/RawDigit.h"
#include "RawData/RawChannelId.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "RawData/RawDaqHeaderBlock.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawVarcErrorInTfBlock.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/LinearFit.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliStripHandle.h"

#include "MINF_Classes/MINFast.h"
#include "REROOT_Classes/REROOT_NeuKin.h"
#include "REROOT_Classes/REROOT_NeuVtx.h"

ClassImp(FitTrackMSListModule)

CVSID("$Id: FitTrackMSListModule.cxx,v 1.11 2010/01/06 18:45:41 rhatcher Exp $");

JOBMODULE(FitTrackMSListModule,"FitTrackMSListModule",
          "Builds CandFitTrackMSList from CandTrackList");


FitTrackMSListModule::FitTrackMSListModule() :
    fListIn("CandTrackList"),fListOut("CandFitTrackMSList"),fFile(0)
{
   MSG("FitTrackMS", Msg::kVerbose) << "FitTrackMSListModule::Constructor\n";

  // Get Singleton instance of AlgFactory.
  AlgFactory &af = AlgFactory::GetInstance();

  // Register (Algorithm, configset) by names ("AlgFitTrackMS", "default")
  af.Register("AlgFitTrackMS", "default", "libCandFitTrackMS.so", "AlgConfig");

  // Register (Algorithm, configset) by names ("AlgFitTrackMSList", "default")
  af.Register("AlgFitTrackMSList", "default", "libCandFitTrackMS.so",
              "AlgConfig");
}


FitTrackMSListModule::~FitTrackMSListModule() 
{
   MSG("FitTrackMS", Msg::kVerbose) << "FitTrackMSListModule::Destructor\n";
}

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

   MSG("FitTrackMS", Msg::kVerbose) << "FitTrackMSListModule::Reco\n";

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

// Find CandSliceList fragment in PrimaryCandidateRecord.
   MSG("FitTrackMS", Msg::kVerbose)
      << "CandSliceListHandle *cslh = (CandSliceListHandle *) "
      << "candrec->FindCandHandle(\"CandSliceListHandle\", "
      << "fListIn.Data());" << endl;
   /* create a pointer which is a slice list handle obtained from the
      candidate record.  dynamic cast deals with the function
      FindCANDIDATEHandle ; cast casts the type from generic handle to 
      CandidateSliceListHandle -- this makes something which is a
      pointer to a handle*/
   CandSliceListHandle *cslh = dynamic_cast<CandSliceListHandle *>
      (candrec->FindCandHandle("CandSliceListHandle"));
   CandTrackListHandle *ctlh = dynamic_cast<CandTrackListHandle *>
      (candrec->FindCandHandle("CandTrackListHandle"));

   TObjArray cxin;
   cxin.Add(cslh);
   cxin.Add(ctlh);


// Get Singleton instance of AlgFactory.
   MSG("FitTrackMS", Msg::kDebug)
                      << "Get Singleton instance of AlgFactory." << endl
               << "AlgFactory &af = AlgFactory::GetInstance();" << endl;
   /* &af is a reference to an algfactory object, GetInstance returns
      a reference to an AlgFactory*/
   AlgFactory &af = AlgFactory::GetInstance();

// Build a CandFitTrackMSList containing all CandFitTrackMSs in Frame.
   MSG("FitTrackMS", Msg::kDebug)
        << "Ask AlgFactory for Singleton AlgFitTrackMSList instance." << endl
   << "AlgHandle adlh = af.GetAlgHandle(\"AlgFitTrackMSList\", \"default\");"
                                                                << endl;
   AlgHandle adlh = af.GetAlgHandle("AlgFitTrackMSList", "default");

   MSG("FitTrackMS", Msg::kDebug) << "CandContext cx(this);" << endl;
   CandContext cx(this, mom);

   MSG("FitTrackMS", Msg::kDebug) << "cx.SetDataIn(cslh);" << endl;
   cx.SetDataIn(&cxin);
   cx.SetCandRecord(candrec);


   MSG("FitTrackMS", Msg::kDebug)
            << "ctllh = CandFitTrackMSList::MakeCandidate(adlh, cx);" << endl;
   CandFitTrackMSListHandle ctllh = CandFitTrackMSList::MakeCandidate(adlh, cx);
   ctllh.SetName(fListOut.Data());
   ctllh.SetTitle(TString("Created by FitTrackMSListModule from ").
                 Append(ctlh->GetName()));

// Give the CandHandle to the CandRecord
   MSG("FitTrackMS", Msg::kDebug) << "candrec->SecureCandHandle(ctllh);"
                                                                << endl;
   candrec->SecureCandHandle(ctllh);

   return result;
}

const Registry &FitTrackMSListModule::DefaultConfig() const
{
  // Returns registry item containing default configuration
  // for the module

  MSG("FitTrackMS",Msg::kDebug) <<
    "FitTrackMSListModule::DefaultConfig" << endl;

  static Registry r;

  std::string name = this->JobCModule::GetName();
  name += ".config.default";
  r.SetName(name.c_str());
  r.UnLockValues();

  // Set Nofit to nonzero for testing purposes to skip fitting step.
  r.Set("Nofit",0);
  r.Set("NoBField",0);
  r.Set("NoMS",0);
  r.Set("BothFit",1);
  r.Set("FullAna",0);
  r.Set("BFisFlipped",1);
  
  // physical constants and data
  r.Set("PosErr",1.18*Munits::cm);
  r.Set("XZero",1.76*Munits::cm);
  r.Set("Dedx",1.454*Munits::GeV / Munits::m);
  r.Set("SuperModGapSize",1.593*Munits::m);
  r.Set("SuperModSkippedPlane",250);

  // parameters
  r.Set("MaxHits",400);
  r.Set("MinHits",5);
  r.Set("MaxP",100*Munits::GeV);
  r.Set("MinP",0.01*Munits::GeV);
  r.Set("MaxIter",30);
  r.Set("InShower",99);
  r.Set("LTolerance",.01);
  r.Set("PTolerance",.01);

  r.LockValues();

  return r;
}


void FitTrackMSListModule::Config(const Registry &r)
{
  // Configures the module, given a registry.  Uses AlgConfig

  MSG("FitTrackMS", Msg::kDebug) << "FitTrackMSListModule::Config" <<
    endl;
  
  // Get Singleton instance of AlgFactory.
  AlgFactory &af = AlgFactory::GetInstance();
  AlgHandle adf = af.GetAlgHandle("AlgFitTrackMS","default");
  AlgConfig &acdf = adf.GetAlgConfig();

  acdf.UnLockValues();
  acdf.Merge(r);
  acdf.LockValues();
}


JobCResult FitTrackMSListModule::Ana(const MomNavigator *mom)
{

  JobCResult result(JobCResult::kPassed);


  MSG("FitTrackMS", Msg::kDebug) << "FitTrackMSListModule::Ana\n";

  if (!fFile) {
    fFile = new TFile("fittrackms.root","RECREATE");
  
    fStripNt = new TTree("stripnt","StripSR Tree");
    fStripNt->SetAutoSave(100000);
    fStripNt->Branch("run",&nt_run,"run/I");
    fStripNt->Branch("snarl",&nt_snarl,"snarl/I");
    fStripNt->Branch("timesec",&nt_timesec,"timesec/I");
    fStripNt->Branch("timens",&nt_timens,"timens/D");
    fStripNt->Branch("track",&nt_track,"track/I");
    fStripNt->Branch("plane",&nt_plane,"plane/i2");
    fStripNt->Branch("planeview",&nt_planeview,"planeview/I1");
    fStripNt->Branch("strip",&nt_strip,"strip/i1");
    fStripNt->Branch("tpos",&nt_tpos,"tpos/F");
    fStripNt->Branch("inshower",&nt_inshower,"inshower/i1");
    fStripNt->Branch("adc",&nt_adc,"adc[2]/I2");
    fStripNt->Branch("adcplane",&nt_adcplane,"adcplane/F");
    fStripNt->Branch("time",&nt_time,"time[2]/D");
    fStripNt->Branch("corrtime",&nt_corrtime,"corrtime[2]/D");
    fStripNt->Branch("residtime",&nt_residtime,"residtime[2]/D");
    fStripNt->Branch("crate",&nt_crate,"crate[2]/I1");
    fStripNt->Branch("varc",&nt_varc,"varc[2]/I1");
    fStripNt->Branch("vmm",&nt_vmm,"vmm[2]/I1");
    fStripNt->Branch("vaadc",&nt_vaadc,"vaadc[2]/I1");
    fStripNt->Branch("vachip",&nt_vachip,"vachip[2]/I1");
    fStripNt->Branch("vachannel",&nt_vachannel,"vachannel[2]/I1");
    fStripNt->Branch("pixel",&nt_pixel,"pixel[2]/I1");
    fStripNt->Branch("uvz",&nt_uvz,"uvz[3]/F");
    fStripNt->Branch("vtxuvz",&nttrack_uvz,"vtxuvz[3]/F");
    fStripNt->Branch("vtxdircos",&nttrack_dircos,"vtxdircos[3]/F");
    fStripNt->Branch("houghdircos",&nttrack_houghdircos,"houghdircos[3]/F");
    fStripNt->Branch("fitdircos",&nttrack_fitdircos,"fitdircos[3]/F");
    fStripNt->Branch("fitintercept",&nttrack_fitintercept,"fitintercept[2]/F");
    fStripNt->Branch("enduvz",&nttrack_uvzend,"enduvz[3]/F");
    fStripNt->Branch("enddircos",&nttrack_dircosend,"enddircos[3]/F");
    fStripNt->Branch("dircos",&nt_dircos,"dircos[3]/F");
    fStripNt->Branch("clearfiber",&nt_clearfiber,"clearfiber[2]/F");
    fStripNt->Branch("wlspigtail",&nt_wlspigtail,"wlspigtail[2]/F");
    fStripNt->Branch("halflength",&nt_halflength,"halflength/F");
    fStripNt->Branch("timeslope",&nttrack_timeslope,"timeslope/F");
    fStripNt->Branch("pmtph",&nt_pmtph,"pmtph[2]/F");
    fStripNt->Branch("houghchi2",&nttrack_houghchi2,"houghchi2/F");
    fStripNt->Branch("attndist",&nt_attndist,"attndist[2]/F");
    fStripNt->Branch("c1",&nt_attnc1,"c1[2]/F");
    fStripNt->Branch("c2",&nt_attnc2,"c2[2]/F");
    fStripNt->Branch("l1",&nt_attnl1,"l1[2]/F");
    fStripNt->Branch("l2",&nt_attnl2,"l2[2]/F");
    fStripNt->Branch("attnnorm",&nt_attnnorm,"attnnorm[2]/F");
    fStripNt->Branch("attncorr",&nt_attncorr,"attncorr[2]/F");
    fStripNt->Branch("ds",&nt_ds,"ds/F");
  
    fPlaneNt = new TTree("planent","PlaneSR Tree");
    fPlaneNt->SetAutoSave(10000);
    fPlaneNt->Branch("run",&nt_run,"run/I");
    fPlaneNt->Branch("snarl",&nt_snarl,"snarl/I");
    fPlaneNt->Branch("timesec",&nt_timesec,"timesec/I");
    fPlaneNt->Branch("timens",&nt_timens,"timens/D");
    fPlaneNt->Branch("track",&nt_track,"track/I");
    fPlaneNt->Branch("nstrip",&ntplane_nstrip,"nstrip/I");
    fPlaneNt->Branch("adc",&ntplane_adc,"adc/I");
    fPlaneNt->Branch("plane",&ntplane_plane,"plane/i2");
    fPlaneNt->Branch("uvz",&ntplane_uvz,"uvz[3]/F");
    fPlaneNt->Branch("dircos",&ntplane_dircos,"dircos[3]/F");
  
    fTrackNt = new TTree("tracknt","FitTrackMS Tree");
    fTrackNt->SetAutoSave(10000);
    fTrackNt->Branch("run",&nt_run,"run/I");
    fTrackNt->Branch("snarl",&nt_snarl,"snarl/I");
    fTrackNt->Branch("timesec",&nt_timesec,"timesec/I");
    fTrackNt->Branch("timens",&nt_timens,"timens/D");
    fTrackNt->Branch("track",&nt_track,"track/I");
    fTrackNt->Branch("nstrip",&nttrack_nstrip,"nstrip/I");
    fTrackNt->Branch("ntrackstrip",&nttrack_ntrackstrip,"ntrackstrip/I");
    fTrackNt->Branch("ndigit",&nttrack_ndigit,"ndigit/I");
    fTrackNt->Branch("ntrackdigit",&nttrack_ntrackdigit,"ntrackdigit/I");
    fTrackNt->Branch("ntimefitdigit",&nttrack_ntimefitdigit,"ntimefitdigit/I");
    fTrackNt->Branch("adc",&nttrack_adc,"adc/I");
    fTrackNt->Branch("nplane",&nttrack_nplane,"nplane/i2");
    fTrackNt->Branch("nplanetrack",&nttrack_nplane_track,"nplanetrack/i2");
    fTrackNt->Branch("uvz",&nttrack_uvz,"uvz[3]/F");
    fTrackNt->Branch("dircos",&nttrack_dircos,"dircos[3]/F");
    fTrackNt->Branch("houghdircos",&nttrack_houghdircos,"houghdircos[3]/F");
    fTrackNt->Branch("fitdircos",&nttrack_fitdircos,"fitdircos[3]/F");
    fTrackNt->Branch("uvzend",&nttrack_uvzend,"uvzend[3]/F");
    fTrackNt->Branch("dircosend",&nttrack_dircosend,"dircosend[3]/F");
    fTrackNt->Branch("timeslope",&nttrack_timeslope,"timeslope/F");
    fTrackNt->Branch("timefitchi2",&nttrack_timefitchi2,"timefitchi2/F");
    fTrackNt->Branch("houghchi2",&nttrack_houghchi2,"houghchi2/F");
    fTrackNt->Branch("year",&nttrack_date_year,"year/i3");
    fTrackNt->Branch("month",&nttrack_date_month,"month/i2");
    fTrackNt->Branch("day",&nttrack_date_day,"day/i2");
    fTrackNt->Branch("hour",&nttrack_date_hour,"hour/i2");
    fTrackNt->Branch("minute",&nttrack_date_minute,"minute/i2");
    fTrackNt->Branch("second",&nttrack_date_second,"second/i2");
    fTrackNt->Branch("nanosecond",&nttrack_date_nanosecond,"nanosecond/i");
    fTrackNt->Branch("snarl_nstrip",&nttrack_snarl_nstrip,"snarl_nstrip/i");
    fTrackNt->Branch("snarl_begplane",&nttrack_snarl_begplane,"snarl_begplane/i2");
    fTrackNt->Branch("snarl_endplane",&nttrack_snarl_endplane,"snarl_endplane/i2");
    fTrackNt->Branch("snarl_nplane",&nttrack_snarl_nplane,"snarl_nplane/i2");
    fTrackNt->Branch("snarl_adc",&nttrack_snarl_adc,"snarl_adc/F");
    fTrackNt->Branch("utime",&nttrack_utime,"utime[2]/F");
    fTrackNt->Branch("vtime",&nttrack_vtime,"vtime[2]/F");
    fTrackNt->Branch("nplanefaildemux",&nttrack_nplanefaildemux,"nplanefaildemux/i2");
    fTrackNt->Branch("begdigitstrip",&nttrack_begdigitstrip,"begdigitstrip[2]/F");
    fTrackNt->Branch("dplanebeg",&nttrack_dplanebeg,"dplanebeg/I");
    fTrackNt->Branch("dplaneend",&nttrack_dplaneend,"dplaneend/I");

    fTrackNt->Branch("momentum",&nttrack_momentum,"momentum/F");
    fTrackNt->Branch("momentumL",&nttrack_momentumL,"momentumL/F");
    fTrackNt->Branch("momentumBF",&nttrack_momentumBF,"momentumBF/F");
    fTrackNt->Branch("momentumMS",&nttrack_momentumMS,"momentumMS/F");
    fTrackNt->Branch("momentumBoth",&nttrack_momentumBoth,"momentumBoth/F");
    fTrackNt->Branch("momentumAlt",&nttrack_momentumAlt,"momentumAlt/F");
    fTrackNt->Branch("chi2",&nttrack_chi2,"chi2/F");
    fTrackNt->Branch("chi2L",&nttrack_chi2L,"chi2L/F");
    fTrackNt->Branch("chi2BF",&nttrack_chi2BF,"chi2BF/F");
    fTrackNt->Branch("chi2MS",&nttrack_chi2MS,"chi2MS/F");
    fTrackNt->Branch("chi2Both",&nttrack_chi2Both,"chi2Both/F");
    fTrackNt->Branch("chi2Alt",&nttrack_chi2Alt,"chi2Alt/F");
 
    fTrackNt->Branch("ds",&nttrack_ds,"ds/F");
    fTrackNt->Branch("maxradius",&nttrack_maxradius,"maxradius/F");

    fTrackNt->Branch("charge",&nttrack_charge,"charge/F");
    fTrackNt->Branch("charged",&nttrack_charged,"charged/F");
    fTrackNt->Branch("flag",&nttrack_flag,"flag/I");
    fTrackNt->Branch("iter",&nttrack_iter,"iter/I");

    if (gMINFast) {
      fTrackNt->Branch("inu",&ntreroot_inu,"inu/I1");
      fTrackNt->Branch("iaction",&ntreroot_iaction,"inu/I1");
      fTrackNt->Branch("p4nu",&ntreroot_p4nu,"p4nu[4]/F");
      fTrackNt->Branch("p4sh",&ntreroot_p4sh,"p4sh[4]/F");
      fTrackNt->Branch("p4mu",&ntreroot_p4mu,"p4mu[4]/F");
      fTrackNt->Branch("p4el",&ntreroot_p4el,"p4el[4]/F");
      fTrackNt->Branch("xyz",&ntreroot_xyz,"xyz[3]/F");
    }
  
  }


  REROOT_Event *ev = 0;
  REROOT_NeuKin *rneukin = 0;
  REROOT_NeuVtx *rneuvtx = 0;
  if (gMINFast) {
    ev = gMINFast->GetREROOTEvent();
    rneukin = dynamic_cast<REROOT_NeuKin*>(ev->neukins()->First());
    rneuvtx = dynamic_cast<REROOT_NeuVtx*>(ev->neuvtxs()->First());
  }

  ntreroot_inu = 0;
  ntreroot_iaction = 0;
  for (int i=0; i<4; i++) {
    ntreroot_p4nu[i] = 0.;
    ntreroot_p4sh[i] = 0.;
    ntreroot_p4mu[i] = 0.;
    ntreroot_p4el[i] = 0.;
  }
  for (int i=0; i<3; i++) {
    ntreroot_xyz[i] = 0.;
  }

  if (rneuvtx && rneukin) {
    ntreroot_inu = rneukin->INu();
    ntreroot_iaction = rneukin->IAction();
    for (int i=0; i<4; i++) {
      ntreroot_p4nu[i] = rneukin->P4Neu()[i];
      ntreroot_p4sh[i] = rneukin->P4Shw()[i];
      ntreroot_p4mu[i] = rneukin->P4Mu1()[i];
      ntreroot_p4el[i] = rneukin->P4El1()[i];
    }
    ntreroot_xyz[0] = rneuvtx->X();
    ntreroot_xyz[1] = rneuvtx->Y();
    ntreroot_xyz[2] = rneuvtx->Z();
  }

  
//  Int_t pixel2vach[17] = {0,3,5,14,16,7,9,10,12,11,13,6,8,17,15,2,4};
  Int_t vach2pixel[18] = {0,0,15,1,16,2,11,5,12,6,7,9,8,10,3,14,4,13};

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

  CandFitTrackListHandle *tracklist = dynamic_cast<CandFitTrackListHandle*>
    (candrec->FindCandHandle("CandFitTrackListHandle"));
  if (!tracklist) {
    MSG("FitTrackMS", Msg::kWarning)
      << "No CandFitTrackListHandle in CandRecord." << endl;
    return result;
  }

  CandDigitListHandle *digitlist = dynamic_cast<CandDigitListHandle*>
    (candrec->FindCandHandle("CandDigitListHandle"));

  CandSliceListHandle *slicelist = dynamic_cast<CandSliceListHandle*>
    (candrec->FindCandHandle("CandSliceListHandle"));

  nttrack_snarl_nstrip = 0;
  nttrack_snarl_begplane = 0;
  nttrack_snarl_endplane = 0;
  nttrack_snarl_nplane = 0;
  nttrack_snarl_adc = 0.;

  Int_t planedigit[1000];
  Int_t planestrip[1000];
  Float_t planecharge[1000];
  for (int i=0; i<1000; i++) {
    planedigit[i] = 0;
    planestrip[i] = 0;
    planecharge[i] = 0.;
  }

  if (slicelist && slicelist->GetNDaughters()>0) {
    TIter sliceItr(slicelist->GetDaughterIterator());
    CandSliceHandle *slice = dynamic_cast<CandSliceHandle*>(sliceItr());
    nttrack_snarl_nstrip = slice->GetNStrip();
    nttrack_snarl_begplane = slice->GetBegPlane();
    nttrack_snarl_endplane = slice->GetEndPlane();
    nttrack_snarl_adc = slice->GetCharge();
    TIter stripItr(slice->GetDaughterIterator());
    while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
      Int_t iplane = strip->GetPlane();
      if (iplane>=0 && iplane<1000) {
        planestrip[iplane]++;
        planedigit[iplane] += strip->GetNDigit();
        planecharge[iplane] += strip->GetCharge();
      }
    }
  }

  for (int i=0; i<1000; i++) {
    if (planecharge[i]>200.) {
      nttrack_snarl_nplane++;
    }
  }

  nttrack_nplanefaildemux = 0;
  for (int i=0; i<1000; i++) {
    if (planestrip[i]>0 && (Float_t)(planedigit[i])/(Float_t)(planestrip[i])<1.2) {
      nttrack_nplanefaildemux++;
    }
  }


  // we require a DaqSnarl (ie. has RawDigitDataBlock) RawRecord
  RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord",0,"DaqSnarl"));
  if (rr == 0) {
    MSG("FitTrackMS", Msg::kWarning) << "No RawRecord in MOM." << endl;
    return result;
  }

  const RawDaqSnarlHeader* snarlHdr =
     dynamic_cast<const RawDaqSnarlHeader*>(rr->GetRawHeader());
  if (snarlHdr) {
     nt_run = snarlHdr->GetRun();
     nt_snarl = snarlHdr->GetSnarl();
  }


  const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *>
                        (rr->FindRawBlock("RawDigitDataBlock"));
  if (rddb == 0) {
    MSG("FitTrackMS", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
    return result;

  }

  TIter trackItr(tracklist->GetDaughterIterator());

  nt_timesec = -1;
  nt_timens = -1.;
  if (digitlist) {
    nt_timens = digitlist->GetAbsTime();
  }

  candrec->GetVldContext()->GetTimeStamp().GetTime(kTRUE, 0, &nttrack_date_hour, &nttrack_date_minute, &nttrack_date_second);
  candrec->GetVldContext()->GetTimeStamp().GetDate(kTRUE, 0, &nttrack_date_year, &nttrack_date_month, &nttrack_date_day);

  nttrack_date_nanosecond = candrec->GetVldContext()->GetTimeStamp().GetNanoSec();

  nttrack_nstrip = 0;
  nttrack_ndigit = 0;
  nttrack_nplane = 0;
  nttrack_nplane_track = 0;
  nttrack_adc = 0;
  nttrack_uvz[0] = 0.;
  nttrack_uvz[1] = 0.;
  nttrack_uvz[2] = 0.;
  nttrack_dircos[0] = 0.;
  nttrack_dircos[1] = 0.;
  nttrack_dircos[2] = 0.;
  nttrack_uvzend[0] = 0.;
  nttrack_uvzend[1] = 0.;
  nttrack_uvzend[2] = 0.;
  nttrack_dircosend[0] = 0.;
  nttrack_dircosend[1] = 0.;
  nttrack_dircosend[2] = 0.;
  nttrack_houghdircos[0] = 0.;
  nttrack_houghdircos[1] = 0.;
  nttrack_houghdircos[2] = 0.;
  nttrack_fitdircos[0] = 0.;
  nttrack_fitdircos[1] = 0.;
  nttrack_fitdircos[2] = 0.;
  nttrack_fitintercept[0] = 0.;
  nttrack_fitintercept[1] = 0.;
  nttrack_houghchi2 = 0.;
  nttrack_utime[0] = 0.;
  nttrack_utime[1] = 0.;
  nttrack_vtime[0] = 0.;
  nttrack_vtime[1] = 0.;
  nttrack_dplanebeg = 0;
  nttrack_dplaneend = 0;
  nttrack_momentum = 0.;
  nttrack_ds = 0.;
  nttrack_maxradius = 0.;
  nttrack_begdigitstrip[0] = 0.;
  nttrack_begdigitstrip[1] = 0.;
  nttrack_ntrackstrip = 0;
  nttrack_ntrackdigit = 0;
  nttrack_ntimefitdigit = 0;
  nttrack_dircosend[0] = 0.;
  nttrack_dircosend[1] = 0.;
  nttrack_dircosend[2] = 0.;
  nttrack_houghdircos[0] = 0.;
  nttrack_houghdircos[1] = 0.;
  nttrack_houghdircos[2] = 0.;
  nttrack_houghchi2 = 0.;
  nttrack_nplane_track = 0;
  nttrack_timefitchi2 = 0;
  nttrack_timeslope = 0.;

  nttrack_momentumL = 0.;
  nttrack_momentumMS = 0.;
  nttrack_momentumBF = 0.;
  nttrack_momentumBoth = 0.;
  nttrack_momentumAlt = 0.;
  nttrack_chi2 = 0.;
  nttrack_chi2L = 0.;
  nttrack_chi2MS = 0.;
  nttrack_chi2BF = 0.;
  nttrack_chi2Both = 0.;
  nttrack_chi2Alt = 0.;
  
  nttrack_flag = 0;
  nttrack_iter = 0;
  nttrack_charge = 0.;
  nttrack_charged = 0.;

  
  DbiResultPtr<CalMapperFits> fMapperFitsPtr;
  fMapperFitsPtr.NewQuery(*candrec->GetVldContext(),0);

  Int_t itrack=0;
  while (CandFitTrackHandle *track = dynamic_cast<CandFitTrackHandle*>(trackItr())) {
    CandFitTrackMSHandle *fittrackms = 0;
    if (track->InheritsFrom("CandFitTrackMSHandle")) {
      fittrackms = dynamic_cast<CandFitTrackMSHandle*>(track);
    }
    itrack++;
    nt_track = itrack;
    nttrack_nstrip = track->GetNStrip();
    nttrack_ndigit = track->GetNDigit();
    nttrack_nplane = track->GetNPlane();
    nttrack_nplane_track = 0;
    nttrack_adc = 0;
    nttrack_uvz[0] = track->GetVtxU();
    nttrack_uvz[1] = track->GetVtxV();
    nttrack_uvz[2] = track->GetVtxZ();
    nttrack_dircos[0] = track->GetDirCosU();
    nttrack_dircos[1] = track->GetDirCosV();
    nttrack_dircos[2] = track->GetDirCosZ();
    Int_t endplane = track->GetEndPlane();
    nttrack_uvzend[0] = track->GetU(endplane);
    nttrack_uvzend[1] = track->GetV(endplane);
    nttrack_uvzend[2] = track->GetZ(endplane);
    nttrack_dircosend[0] = 0.;
    nttrack_dircosend[1] = 0.;
    nttrack_dircosend[2] = 0.;
    nttrack_houghdircos[0] = 0.;
    nttrack_houghdircos[1] = 0.;
    nttrack_houghdircos[2] = 0.;
    nttrack_fitdircos[0] = 0.;
    nttrack_fitdircos[1] = 0.;
    nttrack_fitdircos[2] = 0.;
    nttrack_fitintercept[0] = 0.;
    nttrack_fitintercept[1] = 0.;
    nttrack_houghchi2 = 999999.;
    nttrack_utime[0] = 0.;
    nttrack_utime[1] = 0.;
    nttrack_vtime[0] = 0.;
    nttrack_vtime[1] = 0.;
    nttrack_dplanebeg = track->GetBegPlane(PlaneView::kU)-track->GetBegPlane(PlaneView::kV);
    nttrack_dplaneend = track->GetEndPlane(PlaneView::kU)-track->GetEndPlane(PlaneView::kV);
    nttrack_momentum = track->GetMomentum();
    nttrack_ds = 0;   //track->GetdS();
    nttrack_maxradius = 0.;
    Float_t totcharge[2][2] = {{0.,0.},{0.,0.}};
    nttrack_begdigitstrip[0] = (Float_t)(planedigit[track->GetBegPlane(PlaneView::kU)])/(Float_t)(planestrip[track->GetBegPlane(PlaneView::kU)]);
    nttrack_begdigitstrip[1] = (Float_t)(planedigit[track->GetBegPlane(PlaneView::kV)])/(Float_t)(planestrip[track->GetBegPlane(PlaneView::kV)]);

    if (fittrackms) {

      nttrack_momentumL = fittrackms->GetMomentumL();
      nttrack_momentumBF = fittrackms->GetMomentumBF();
      nttrack_momentumMS = fittrackms->GetMomentumMS();
      nttrack_momentumBoth = fittrackms->GetMomentumBoth();
      nttrack_momentumAlt = fittrackms->GetMomentumAlt();

      nttrack_chi2 = fittrackms->GetChi2();
      nttrack_chi2L = fittrackms->GetChi2L();
      nttrack_chi2BF = fittrackms->GetChi2BF();
      nttrack_chi2MS = fittrackms->GetChi2MS();
      nttrack_chi2Both = fittrackms->GetChi2Both();
      nttrack_chi2Alt = fittrackms->GetChi2Alt();

      nttrack_flag = fittrackms->GetFlag();
      nttrack_iter = fittrackms->GetIter();
      nttrack_charge = fittrackms->GetEMCharge();
      nttrack_charged = fittrackms->GetEMChargeD();
    }
    nttrack_timeslope = fabs(track->GetTimeSlope());

    Double_t timeslope = fabs(track->GetTimeSlope());
    Double_t timeoffset = track->GetTimeOffset();

    TIter stripItr(track->GetDaughterIterator());

    Double_t uzfit[1000],ufit[1000],uwfit[1000],uph[1000];
    Double_t vzfit[1000],vfit[1000],vwfit[1000],vph[1000];
  
    for (int i=0; i<1000; i++) {
      uzfit[i] = 0.;
      ufit[i] = 0.;
      uwfit[i] = 0.;
      uph[i] = 0.;
      vzfit[i] = 0.;
      vfit[i] = 0.;
      vwfit[i] = 0.;
      vph[i] = 0.;
    }

    stripItr.Reset();
    while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
      nt_plane = strip->GetPlane();
      nt_tpos = strip->GetTPos();
      if (strip->GetPlaneView()==PlaneView::kU) {
        if (strip->GetNDigit(StripEnd::kNegative)>0) {
          totcharge[0][0] += strip->GetCharge(StripEnd::kNegative);
          nttrack_utime[0] += strip->GetCharge(StripEnd::kNegative)*track->GetT(nt_plane,StripEnd::kNegative);
        }
        if (strip->GetNDigit(StripEnd::kPositive)>0) {
          totcharge[0][1] += strip->GetCharge(StripEnd::kPositive);
          nttrack_utime[1] += strip->GetCharge(StripEnd::kPositive)*track->GetT(nt_plane,StripEnd::kPositive);
        }
        uzfit[nt_plane] = strip->GetZPos();
        ufit[nt_plane] += nt_tpos*strip->GetCharge();
        uph[nt_plane] += strip->GetCharge();
        uwfit[nt_plane] = 1.;
      } else if (strip->GetPlaneView()==PlaneView::kV) {
        if (strip->GetNDigit(StripEnd::kNegative)>0) {
          totcharge[1][0] += strip->GetCharge(StripEnd::kNegative);
          nttrack_vtime[0] += strip->GetCharge(StripEnd::kNegative)*track->GetT(nt_plane,StripEnd::kNegative);
        }
        if (strip->GetNDigit(StripEnd::kPositive)>0) {
          totcharge[1][1] += strip->GetCharge(StripEnd::kPositive);
          nttrack_vtime[1] += strip->GetCharge(StripEnd::kPositive)*track->GetT(nt_plane,StripEnd::kPositive);
        }
        vzfit[nt_plane] = strip->GetZPos();
        vfit[nt_plane] += nt_tpos*strip->GetCharge();
        vph[nt_plane] += strip->GetCharge();
        vwfit[nt_plane] = 1.;
      }
    }

    for (int i=0; i<2; i++) {
      if (totcharge[0][i]>0.) {
        nttrack_utime[i] /= totcharge[0][i];
      }
      else {
        nttrack_utime[i] = -1.;
      }
      if (totcharge[1][i]>0.) {
        nttrack_vtime[i] /= totcharge[1][i];
      }
      else {
        nttrack_vtime[i] = -1.;
      }
    }

    for (int i=0; i<1000; i++) {
      if (uph[i]>0.) {
        ufit[i] /= uph[i];
      }
      if (vph[i]>0.) {
        vfit[i] /= vph[i];
      }
    }

    Double_t uparm[2],ueparm[2];
    Double_t vparm[2],veparm[2];

    LinearFit::Weighted(1000,uzfit,ufit,uwfit,uparm,ueparm);
    LinearFit::Weighted(1000,vzfit,vfit,vwfit,vparm,veparm);

    Double_t dudz = uparm[1];
    Double_t dvdz = vparm[1];
    nttrack_fitdircos[0] = dudz/sqrt(1.+dudz*dudz+dvdz*dvdz);
    nttrack_fitdircos[1] = dvdz/sqrt(1.+dudz*dudz+dvdz*dvdz);
    nttrack_fitdircos[2] = 1./sqrt(1.+dudz*dudz+dvdz*dvdz);

    nttrack_fitintercept[0] = uparm[0];
    nttrack_fitintercept[1] = vparm[0];

    if (track->GetTimeSlope()<0.) {
      nttrack_fitdircos[0] *= -1.;
      nttrack_fitdircos[1] *= -1.;
      nttrack_fitdircos[2] *= -1.;
    }

    stripItr.Reset();
    while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
      for (int i=0; i<2; i++) {
        nt_adc[i] = -1;
        nt_time[i] = -1.;
        nt_corrtime[i] = -1.;
        nt_residtime[i] = -1.;
        nt_crate[i] = -1;
        nt_varc[i] = -1;
        nt_vmm[i] = -1;
        nt_vaadc[i] = -1;
        nt_vachip[i] = -1;
        nt_vachannel[i] = -1;
        nt_pixel[i] = -1;
        nt_pmtph[i] = 0.;
      }
      nt_plane = strip->GetPlane();
      nt_planeview = (Int_t)strip->GetPlaneView();
      nt_strip = strip->GetStrip();
      nt_tpos = strip->GetTPos();
      nt_inshower = track->IsInShower(strip);
      nt_uvz[0] = track->GetU(nt_plane);
      nt_uvz[1] = track->GetV(nt_plane);
      nt_uvz[2] = track->GetZ(nt_plane);
      //      if (tracksr) {
      //        nt_dircos[0] = tracksr->GetDirCosU(nt_plane);
      //        nt_dircos[1] = tracksr->GetDirCosV(nt_plane);
      //        nt_dircos[2] = tracksr->GetDirCosZ(nt_plane);
      //      }
      nt_adcplane = planecharge[nt_plane];
      nt_ds = 0; //track->GetdS(nt_plane);
      TIter digitItr(strip->GetDaughterIterator());
      PlexStripEndId *plexstripendid=0;
      UgliGeomHandle *ugh = 0;
      while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())) {
        if (!plexstripendid) {
          plexstripendid = new PlexStripEndId(digit->GetPlexSEIdAltL().GetBestSEId());
          ugh = new UgliGeomHandle(*digit->GetVldContext());
        }
        const RawDigit *rd = rddb->At(digit->GetRawDigitIndex());
        if (rd) nttrack_adc += rd->GetADC();
        switch (digit->GetPlexSEIdAltL().GetBestSEId().GetEnd()) {
        case StripEnd::kNegative: case StripEnd::kWhole: case StripEnd::kUnknown:
          nt_time[0] = strip->GetBegTime(digit->GetPlexSEIdAltL().GetBestSEId().GetEnd());
          nt_corrtime[0] = track->GetT(strip->GetPlane(),StripEnd::kNegative);
          nt_residtime[0] = nt_corrtime[0]-(timeoffset+timeslope*nt_ds);
          if (rd) {
            RawChannelId rawid = rd->GetChannel();
            nt_timesec = rd->GetCrateT0().GetSec();
            nt_adc[0] = rd->GetADC();
            nt_crate[0] = rawid.GetCrate();
            nt_varc[0] = rawid.GetVarcId();
            nt_vmm[0] = rawid.GetVmm();
            nt_vaadc[0] = rawid.GetVaAdcSel();
            nt_vachip[0] = rawid.GetVaChip();
            nt_vachannel[0] = rawid.GetVaChannel();
            nt_pixel[0] = -1;
            if (nt_vachannel[0]>=2 && nt_vachannel[0]<=17) {
              nt_pixel[0] = vach2pixel[nt_vachannel[0]];
            }
          }
          break;
        case StripEnd::kPositive:
          nt_time[1] = strip->GetBegTime(digit->GetPlexSEIdAltL().GetBestSEId().GetEnd());
          nt_corrtime[1] = track->GetT(strip->GetPlane(),StripEnd::kPositive);
          nt_residtime[1] = nt_corrtime[1]-(timeoffset+timeslope*nt_ds);
          if (rd) {
            RawChannelId rawid = rd->GetChannel();
            nt_timesec = rd->GetCrateT0().GetSec();
            nt_adc[1] = rd->GetADC();
            nt_crate[1] = rawid.GetCrate();
            nt_varc[1] = rawid.GetVarcId();
            nt_vmm[1] = rawid.GetVmm();
            nt_vaadc[1] = rawid.GetVaAdcSel();
            nt_vachip[1] = rawid.GetVaChip();
            nt_vachannel[1] = rawid.GetVaChannel();
            nt_pixel[1] = -1;
            if (nt_vachannel[1]>=2 && nt_vachannel[1]<=17) {
              nt_pixel[1] = vach2pixel[nt_vachannel[1]];
            }
          }
          break;
        default:
          MSG("FitTrackMS",Msg::kError) << "undefined StripEnd " << digit->GetPlexSEIdAltL().GetBestSEId().GetEnd() << "\n";
          break;
        }
      }
      const CalMapperFits *mapcal[2];
      mapcal[0] = fMapperFitsPtr.GetRowByIndex(strip->GetStripEndId(StripEnd::kNegative).BuildPlnStripEndKey());
      mapcal[1] = fMapperFitsPtr.GetRowByIndex(strip->GetStripEndId(StripEnd::kPositive).BuildPlnStripEndKey());
      nt_attnc1[0] = 0.;
      nt_attnc2[0] = 0.;
      nt_attnl1[0] = 0.;
      nt_attnl2[0] = 0.;
      nt_attnnorm[0] = 0.;
      nt_attnc1[1] = 0.;
      nt_attnc2[1] = 0.;
      nt_attnl1[1] = 0.;
      nt_attnl2[1] = 0.;
      nt_attnnorm[1] = 0.;
      if (mapcal[0]) {
        nt_attnc1[0] = mapcal[0]->GetC1();
        nt_attnc2[0] = mapcal[0]->GetC2();
        nt_attnl1[0] = mapcal[0]->GetLambda1()/100.;
        nt_attnl2[0] = mapcal[0]->GetLambda2()/100.;
        nt_attnnorm[0] = mapcal[0]->GetNorm();
      }
      if (mapcal[1]) {
        nt_attnc1[1] = mapcal[1]->GetC1();
        nt_attnc2[1] = mapcal[1]->GetC2();
        nt_attnl1[1] = mapcal[1]->GetLambda1()/100.;
        nt_attnl2[1] = mapcal[1]->GetLambda2()/100.;
        nt_attnnorm[1] = mapcal[1]->GetNorm();
      }
      UgliStripHandle striphandle = ugh->GetStripHandle(*plexstripendid);
      nt_clearfiber[0] = striphandle.ClearFiber(StripEnd::kNegative);
      nt_clearfiber[1] = striphandle.ClearFiber(StripEnd::kPositive);
      nt_wlspigtail[0] = striphandle.WlsPigtail(StripEnd::kNegative);
      nt_wlspigtail[1] = striphandle.WlsPigtail(StripEnd::kPositive);
      nt_halflength = striphandle.GetHalfLength();
      nt_attndist[0] = 0.;
      nt_attndist[1] = 0.;
      if (strip->GetPlaneView()==PlaneView::kU) {
        nt_attndist[0] = nt_halflength+nt_wlspigtail[0]-nt_uvz[1];
        nt_attndist[1] = nt_halflength+nt_wlspigtail[1]+nt_uvz[1];
      }
      if (strip->GetPlaneView()==PlaneView::kV) {
        nt_attndist[0] = nt_halflength+nt_wlspigtail[0]+nt_uvz[0];
        nt_attndist[1] = nt_halflength+nt_wlspigtail[1]-nt_uvz[0];
      }
      for (int iside=0; iside<2; iside++) {
        nt_attncorr[iside] = nt_attnc1[iside]*exp(-nt_attndist[iside]/nt_attnl1[iside])+nt_attnc2[iside]*exp(-nt_attndist[iside]/nt_attnl2[iside]);
      }
      if (digitlist) {
        TIter alldigitItr(digitlist->GetDaughterIterator());
        while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(alldigitItr())) {
          const RawDigit *rd = rddb->At(digit->GetRawDigitIndex());
          if (rd) {
            for (int i=0; i<2; i++) {
              RawChannelId rawid = rd->GetChannel();
              if (rawid.GetCrate()==nt_crate[i] &&
                  rawid.GetVarcId()==nt_varc[i] &&
                  rawid.GetVmm()==nt_vmm[i] &&
                  rawid.GetVaAdcSel()==nt_vaadc[i] &&
                  rawid.GetVaChip()==nt_vachip[i]) {
                nt_pmtph[i] += digit->GetCharge();
              }
            }
          }
        }
      }
      delete plexstripendid;
      delete ugh;
      if (nt_uvz[0]*nt_uvz[0]+nt_uvz[1]*nt_uvz[1]>nttrack_maxradius*nttrack_maxradius) {
        nttrack_maxradius = sqrt(nt_uvz[0]*nt_uvz[0]+nt_uvz[1]*nt_uvz[1]);
      }
      fStripNt->Fill();
    }



    fTrackNt->Fill();

    for (Int_t iplane = track->GetBegPlane(); iplane<=track->GetEndPlane(); iplane++) {
      ntplane_nstrip = 0;
      ntplane_adc = 0;
      ntplane_plane = iplane;
      ntplane_uvz[0] = 0.;
      ntplane_uvz[1] = 0.;
      ntplane_uvz[2] = 0.;
      ntplane_dircos[0] = 0.;
      ntplane_dircos[1] = 0.;
      ntplane_dircos[2] = 0.;
      Float_t totalcharge=0.;
      stripItr.Reset();
      while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
        if (strip->GetPlane()==iplane) {
          ntplane_nstrip++;
          totalcharge += strip->GetCharge();
          ntplane_uvz[0] += strip->GetCharge()*track->GetU(iplane);
          ntplane_uvz[1] += strip->GetCharge()*track->GetV(iplane);
          ntplane_uvz[2] += strip->GetCharge()*track->GetZ(iplane);
          //          if (tracksr) {
          //            ntplane_dircos[0] += strip->GetCharge()*tracksr->GetDirCosU(iplane);
          //            ntplane_dircos[1] += strip->GetCharge()*tracksr->GetDirCosV(iplane);
          //            ntplane_dircos[2] += strip->GetCharge()*tracksr->GetDirCosZ(iplane);
          //          }
          TIter digitItr(strip->GetDaughterIterator());
          while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())) {
            const RawDigit *rd = rddb->At(digit->GetRawDigitIndex());
            if (rd) ntplane_adc += rd->GetADC();
          }
        }
      }
      if (ntplane_nstrip>0) {
        for (int i=0; i<3; i++) {
          ntplane_uvz[i] /= totalcharge;
          ntplane_dircos[i] /= totalcharge;
        }
        fPlaneNt->Fill();
      }
    }
 
  }

  if (!itrack) {
    nt_track = -1;
    fTrackNt->Fill();
  }

  fFile->Write("",TObject::kOverwrite);
  
  return result;
}

void FitTrackMSListModule::EndJob()
{
  if (fFile) {
    fFile->Write("",TObject::kOverwrite);
    fFile->Close();
  }
}

---