#include "muifuncs.h"
#include "recoConsts.h"
#include "RunToTimeObjy.hh"
#include "PHTimeStamp.h"
#include "dMutAnodesDCMWrapper.h"
#include "dMuiRoadRawRelWrapper.h"
#include "dMuiRawWrapper.h"
#include "precoSelect.h"
#include "getMapFileFlag.h"
#include "pisaReco.h"
#include "mfm_setfscale.h"
#include "mfm_getfscale.h"

#define DEBUG 

void 
setup_mui() // for MC tune
{

  PHNodeIterator* mainIter = new PHNodeIterator(topNode);
  PHCompositeNode* muiNode = new PHCompositeNode("MUI");
  topNode->addNode(muiNode);
  PHCompositeNode* mutNode = new PHCompositeNode("MUT");
  topNode->addNode(mutNode);
  PHBoolean IsFile;
  int msgnum = 0;
  size_t mr;

  // Read in par file used in response chain if we are doing reconstruction:

  if (simulationFlag != 3){
    PHNodeIOManager * parIn = new PHNodeIOManager();
    IsFile = parIn->setFile("mupar.root","",PHReadOnly);
      if (!IsFile) {
        cout << "EXIT : File mupar.root does not exist" << endl;
      }
    parIn->read(parNode);
  }

  PHIODataNode<TObject>* dn;

  if(simulationFlag == 2 || simulationFlag == 3) {
    //
    // Insertions here for the initializations needed for
    // the slow simulator part of the detector resposne
    //

    if(simulationFlag == 3) {
      //
      // Insertions here for the initializations needed for
      // the FEM and DCM modules which are needed for PRDF output filling
      //

      // Set up the tables needed for mui

      mr=60000;
      munhitsWrapper* munhits = new munhitsWrapper("munhits",mr);
      PHIODataNode<PHTable>* munhitsNode = new PHIODataNode<PHTable>(munhits,"munhits");
      geaNode->addNode(munhitsNode);

      mr=4000;
      dMuiRawGhitRelWrapper* dMuiRawGhitRel = new dMuiRawGhitRelWrapper("dMuiRawGhitRel",mr);
      PHIODataNode<PHTable>* dMuiRawGhitRelNode = new PHIODataNode<PHTable>(dMuiRawGhitRel,"dMuiRawGhitRel");
      evaNode->addNode(dMuiRawGhitRelNode); // correct the node in deftext

      mr=4000;
      dMuiDCMGhitRelWrapper* dMuiDCMGhitRel = new dMuiDCMGhitRelWrapper("dMuiDCMGhitRel",mr);
      PHIODataNode<PHTable>* dMuiDCMGhitRelNode = new PHIODataNode<PHTable>(dMuiDCMGhitRel,"dMuiDCMGhitRel");
      evaNode->addNode(dMuiDCMGhitRelNode); // correct the node in deftext

      mr=4;
      dMuiFEMWrapper* dMuiFEM = new dMuiFEMWrapper("dMuiFEM",mr);
      PHIODataNode<PHTable>* dMuiFEMNode = new PHIODataNode<PHTable>(dMuiFEM,"dMuiFEM");
      muiNode->addNode(dMuiFEMNode);

      mr=4;
      dMuiDCMWrapper* dMuiDCM = new dMuiDCMWrapper("dMuiDCM",mr);
      PHIODataNode<PHTable>* dMuiDCMNode = new PHIODataNode<PHTable>(dMuiDCM,"dMuiDCM");
      dcmNode->addNode(dMuiDCMNode);

// Set up the tables needed for mut

      mr=60000;
      mumhitsWrapper* mumhits = new mumhitsWrapper("mumhits",mr);
      PHIODataNode<PHTable>* mumhitsNode = new PHIODataNode<PHTable>(mumhits,"mumhits");
      geaNode->addNode(mumhitsNode);

      mr=1;
      mumgeoWrapper* mumgeo = new mumgeoWrapper("mumgeo",mr);
      PHIODataNode<PHTable>* mumgeoNode = new PHIODataNode<PHTable>(mumgeo,"mumgeo");
      parNode->addNode(mumgeoNode);

      mr=1;
      dMutDeadChannelsWrapper* dMutDeadChannels = new dMutDeadChannelsWrapper("dMutDeadChannels",mr);
      PHIODataNode<PHTable>* dMutDeadChannelsNode = new PHIODataNode<PHTable>(dMutDeadChannels,"dMutDeadChannels");
      parNode->addNode(dMutDeadChannelsNode);

     mr=5000;
     dMutCathodeChargeWrapper* dMutCathodeCharge = new dMutCathodeChargeWrapper("dMutCathodeCharge",mr);
     PHIODataNode<PHTable>* dMutCathodeChargeNode = new PHIODataNode<PHTable>(dMutCathodeCharge,"dMutCathodeCharge");
     mutNode->addNode(dMutCathodeChargeNode);

      mr=510;
     dMutCathodesFEMWrapper* dMutCathodesFEM = new dMutCathodesFEMWrapper("dMutCathodesFEM",mr);
     PHIODataNode<PHTable>* dMutCathodesFEMNode = new PHIODataNode<PHTable>(dMutCathodesFEM,"dMutCathodesFEM");
     mutNode->addNode(dMutCathodesFEMNode);

     mr=510;
     dMutCathodesFEMMCWrapper* dMutCathodesFEMMC = new dMutCathodesFEMMCWrapper("dMutCathodesFEMMC",mr);
     PHIODataNode<PHTable>* dMutCathodesFEMMCNode = new PHIODataNode<PHTable>(dMutCathodesFEMMC,"dMutCathodesFEMMC");
     evaNode->addNode(dMutCathodesFEMMCNode);

     mr=20;
     dMutAnodesFEMWrapper* dMutAnodesFEM  = new dMutAnodesFEMWrapper("dMutAnodesFEM ",mr);
     PHIODataNode<PHTable>* dMutAnodesFEMNode = new PHIODataNode<PHTable>(dMutAnodesFEM ,"dMutAnodesFEM ");
     mutNode->addNode(dMutAnodesFEMNode);

     mr=20;
     dMutAnodesFEMMCWrapper* dMutAnodesFEMMC = new dMutAnodesFEMMCWrapper("dMutAnodesFEMMC",mr);
     PHIODataNode<PHTable>* dMutAnodesFEMMCNode = new PHIODataNode<PHTable>(dMutAnodesFEMMC,"dMutAnodesFEMMC");
     evaNode->addNode(dMutAnodesFEMMCNode);
     // ends of the definition of output tables for module mMutXYtoFE

     // begining of table def for mMutZeroSup

     mr=510;
     dMutCathodesDCMWrapper* dMutCathodesDCM = new dMutCathodesDCMWrapper("dMutCathodesDCM",mr);
     PHIODataNode<PHTable>* dMutCathodesDCMNode = new PHIODataNode<PHTable>(dMutCathodesDCM,"dMutCathodesDCM");
     dcmNode->addNode(dMutCathodesDCMNode);

     mr=510;
     dMutCathodesDCMMCWrapper* dMutCathodesDCMMC = new dMutCathodesDCMMCWrapper("dMutCathodesDCMMC",mr);
     PHIODataNode<PHTable>* dMutCathodesDCMMCNode = new PHIODataNode<PHTable>(dMutCathodesDCMMC,"dMutCathodesDCMMC");
     evaNode->addNode(dMutCathodesDCMMCNode);

     mr=20;
     dMutAnodesDCMWrapper* dMutAnodesDCM = new dMutAnodesDCMWrapper("dMutAnodesDCM",mr);
     PHIODataNode<PHTable>* dMutAnodesDCMNode = new PHIODataNode<PHTable>(dMutAnodesDCM,"dMutAnodesDCM");
     dcmNode->addNode(dMutAnodesDCMNode);

     mr=20;
     dMutAnodesDCMMCWrapper* dMutAnodesDCMMC = new dMutAnodesDCMMCWrapper("dMutAnodesDCMMC",mr);
     PHIODataNode<PHTable>* dMutAnodesDCMMCNode = new PHIODataNode<PHTable>(dMutAnodesDCMMC,"dMutAnodesDCMMC");
     evaNode->addNode(dMutAnodesDCMMCNode);

     mr=1;
     dMutCSCGeomWrapper* dMutCSCGeom = new dMutCSCGeomWrapper("dMutCSCGeom",mr);
     PHIODataNode<PHTable>* dMutCSCGeomNode = new PHIODataNode<PHTable>(dMutCSCGeom,"dMutCSCGeom");
     geaNode->addNode(dMutCSCGeomNode);

    } // check on PISA-to-PRDF mode

    if(simulationFlag == 2) {
      // Insertions here for the initializations needed for
      // event reconstruction in PISA-to-DST mode when a PRDF is not written
      //

      cout << "ERROR:  simulatationFlag=2 not yet implemented for Muons" << endl;
      return;

    } // check on PISA-to-DST mode

  } // check on PISA-to-DST or PISA-to-PRDF mode

  // EVAL
  if (evaluationFlag) {

    // Open the geant evaluation file in an elegant way
    //    ioEval = new PHNodeIOManager("rawrel.root",PHReadOnly);
    ioEval = new PHNodeIOManager();
    IsFile = ioEval->setFile("rawrel.root","",PHReadOnly);
    if (!IsFile) {
      cout << "EXIT : File rawrel.root needed for Evaluation does not exist" << endl;
    }

    // Create GEA hits table
    mumhitsWrapper* mumhits = new mumhitsWrapper;
    dn = new PHIODataNode<TObject>(mumhits,"mumhits");
    geaNode->addNode(dn);
    munhitsWrapper* munhits = new munhitsWrapper;
    dn = new PHIODataNode<TObject>(munhits,"munhits");
    geaNode->addNode(dn);
    //EOEVAL
  }

  dMuiPseudoTriggerOut* PseudoTriggerOut = new dMuiPseudoTriggerOutv1();; 
  PHIODataNode<PHObject>* PseudoTriggerOutNode = new PHIODataNode<PHObject>(PseudoTriggerOut, "dMuiPseudoTriggerOut","PHObject");
  dstNode->addNode( PseudoTriggerOutNode );

  dMuiRawWrapper* dMuiRawW = new dMuiRawWrapper("dMuiRaw", 2000);

  PHIODataNode<PHTable>* dMuiRawNode =
    new PHIODataNode<PHTable>(dMuiRawW, "dMuiRaw");

  muiNode->addNode(dMuiRawNode);
  // Now let's try to cheat?
  dstNode->addNode(dMuiRawNode);

  dMuiRoadRawRelWrapper* dMuiRoadRawRelW = new dMuiRoadRawRelWrapper("dMuiRoadRawRel", 120000);

  PHIODataNode<PHTable>* dMuiRoadRawRelNode =
    new PHIODataNode<PHTable>(dMuiRoadRawRelW, "dMuiRoadRawRel");

  muiNode->addNode(dMuiRoadRawRelNode);
  // Now let's try to cheat?
  dstNode->addNode(dMuiRoadRawRelNode);

  mr=1;
  dMuiReadoutControlWrapper* dMuiReadoutControl = new dMuiReadoutControlWrapper("dMuiReadoutControl",mr);
  PHIODataNode<PHTable>* dMuiReadoutControlNode = new PHIODataNode<PHTable>(dMuiReadoutControl,"dMuiReadoutControl");
  parNode->addNode(dMuiReadoutControlNode);

  dMuiRoadControlWrapper* dMuiRoadControl = new
      dMuiRoadControlWrapper;
  dn = new PHIODataNode<TObject>(dMuiRoadControl, "dMuiRoadControl");
  parNode->addNode(dn);

  dMuiClusterControlWrapper* dMuiClusterControl = new
      dMuiClusterControlWrapper;
  dn = new
       PHIODataNode<TObject>(dMuiClusterControl, "dMuiClusterControl");
  parNode->addNode(dn);

  dMuiMuonIdProfileParWrapper* dMuiMuonIdProfilePar = new
      dMuiMuonIdProfileParWrapper;
  dn = new
       PHIODataNode<TObject>(dMuiMuonIdProfilePar, "dMuiMuonIdProfilePar");
  parNode->addNode(dn);

  dMuoFindTracksParWrapper* dMuoFindTracksPar = new
      dMuoFindTracksParWrapper;
  dn = new
       PHIODataNode<TObject>(dMuoFindTracksPar, "dMuoFindTracksPar");
  parNode->addNode(dn);

#ifdef DEBUG
  cout << "Message: " << msgnum++ << endl;
#endif

  PHIODataNode<PHTable>* DMUTGLOBALPARNode = (PHIODataNode<PHTable>*)mainIter->findFirst("PHIODataNode", "dMutGlobalPar");
  dMutGlobalParWrapper* dMutGlobalPar;
  if (DMUTGLOBALPARNode)
    {
      dMutGlobalPar = (dMutGlobalParWrapper*)DMUTGLOBALPARNode->getData();
    }
  else
    {
      dMutGlobalPar = new dMutGlobalParWrapper;
      dn = new PHIODataNode<TObject>(dMutGlobalPar, "dMutGlobalPar");
      parNode->addNode(dn);
    }

#ifdef DEBUG
  cout << "Message: " << msgnum++ << endl;
#endif

  dMutHistosWrapper* dMutHistos = new dMutHistosWrapper;
  dn = new PHIODataNode<TObject>(dMutHistos, "dMutHistos");
  parNode->addNode(dn);

#ifdef DEBUG
  cout << "Message: " << msgnum++ << endl;
#endif

  PHIODataNode<PHTable>* MumDigiParNode = (PHIODataNode<PHTable>*)mainIter->findFirst("PHIODataNode", "mumdigipar");
  mumdigiparWrapper* mumdigipar;
  if (MumDigiParNode)
    {
      mumdigipar = (mumdigiparWrapper*)MumDigiParNode->getData();
    }
  else
    {
      mumdigipar = new mumdigiparWrapper;
      dn = new PHIODataNode<TObject>(mumdigipar, "mumdigipar");
      parNode->addNode(dn);
    }

#ifdef DEBUG
  cout << "Message: " << msgnum++ << endl;
#endif

  PHIODataNode<PHTable>* dMutCathodePulseParamsNode = (PHIODataNode<PHTable>*)mainIter->findFirst("PHIODataNode", "dMutCathodePulseParams");
  mMutCathodePulseParams* dMutCathodePulseParams;
  if (dMutCathodePulseParamsNode){
      dMutCathodePulseParams = (mMutCathodePulseParams*)dMutCathodePulseParamsNode->getData();
  }
  else{
      dMutCathodePulseParams = new mMutCathodePulseParams;
      dn = new PHIODataNode<TObject>(dMutCathodePulseParams, "dMutCathodePulseParams");
      parNode->addNode(dn);
  }

  PHIODataNode<PHTable>* dMfmParNode = (PHIODataNode<PHTable>*)mainIter->findFirst("PHIODataNode", "dMfmPar");
  dMfmParWrapper* dMfmPar;
  if (dMfmParNode){
      dMfmPar = (dMfmParWrapper*)dMfmParNode->getData();
  }
  else{
      dMfmPar = new dMfmParWrapper;
      dn = new PHIODataNode<TObject>(dMfmPar, "dMfmPar");
      parNode->addNode(dn);
  }

  // Define tables which are used for evaluations, but also in the interface
  // of general reconstruction in the mut code:
  // These tables are defined by the evaluation routines, if they are called

  if (!evaluationFlag) {
    dMutAnodesDCMWrapper* dMutAnodesDCM = new dMutAnodesDCMWrapper;
    dn = new PHIODataNode<TObject>(dMutAnodesDCM, "dMutAnodesDCM");
    dcmNode->addNode(dn);

    dMutClusterGhitRelWrapper* dMutClusterGhitRel = new dMutClusterGhitRelWrapper;
    dn = new PHIODataNode<TObject>(dMutClusterGhitRel, "dMutClusterGhitRel");
    dstNode->addNode(dn);
  }

  dMuiRoadFkinRelWrapper* dMuiRoadFkinRel = new dMuiRoadFkinRelWrapper;
  dn = new PHIODataNode<TObject>(dMuiRoadFkinRel, "dMuiRoadFkinRel");
  dstNode->addNode(dn);

  dMuoTrackRoadRelWrapper* dMuoTrackRoadRel =
    new dMuoTrackRoadRelWrapper("dMuoTrackRoadRel", 1000);
  dn = new PHIODataNode<TObject>(dMuoTrackRoadRel, "dMuoTrackRoadRel");
  dstNode->addNode(dn);

#ifdef DEBUG
  cout << "Finished Instantiating Tables\n";
#endif

  dMuiReadoutControl->SetRowCount((size_t)1);
  Int_t i;
  for (i = 0; i < 5; i++) {
    dMuiReadoutControl->set_NeutronEffic(i,0,0.0);
    dMuiReadoutControl->set_GammaEffic(i,0,0.0);
    dMuiReadoutControl->set_ChargeEffic(i,0,97.0);
    dMuiReadoutControl->set_AveNoiseHits(i,0,0.0);
    dMuiReadoutControl->set_RmsNoiseHits(i,0,0.0);
  }


  dMuiRoadControl->SetRowCount((size_t)1);
  dMuiRoadControl->set_NumSeedLoops(0, 2);
  for (Int_t i = 0; i < 2; i++)
    {
      dMuiRoadControl->set_SearchLength(i, 0, 6);
    }
  dMuiRoadControl->set_SearchOrder(0, 0, 0, -1);
  dMuiRoadControl->set_SearchOrder(0, 1, 0, 1);
  dMuiRoadControl->set_SearchOrder(0, 2, 0, 0);
  dMuiRoadControl->set_SearchOrder(0, 3, 0, 2);
  dMuiRoadControl->set_SearchOrder(0, 4, 0, 3);
  dMuiRoadControl->set_SearchOrder(0, 5, 0, 4);
  dMuiRoadControl->set_SearchOrder(1, 0, 0, -1);
  dMuiRoadControl->set_SearchOrder(1, 1, 0, 2);
  dMuiRoadControl->set_SearchOrder(1, 2, 0, 1);
  dMuiRoadControl->set_SearchOrder(1, 3, 0, 0);
  dMuiRoadControl->set_SearchOrder(1, 4, 0, 3);
  dMuiRoadControl->set_SearchOrder(1, 5, 0, 4);
  dMuiRoadControl->set_ClusterCollectMode(0, 0);
  dMuiRoadControl->set_minLastGap1D(0, 2);
  dMuiRoadControl->set_minFiredGaps(0, 2);
  dMuiRoadControl->set_maxSkippedGaps(0, 2);
  dMuiRoadControl->set_minSharedHits1D(0, 5);
  dMuiRoadControl->set_maxXRef1D(0, 180.0);
  dMuiRoadControl->set_maxYRef1D(0, 180.0);
  dMuiRoadControl->set_minLastGap2D(0, 2);
  dMuiRoadControl->set_maxDelLastGap(0, 1);
  dMuiRoadControl->set_maxDelHitsPerGap(0, 1);
  dMuiRoadControl->set_maxDelTotalHits(0, 1);
  dMuiRoadControl->set_maxXRef2D(0, 180.0);
  dMuiRoadControl->set_maxYRef2D(0, 180.0);
  dMuiRoadControl->set_maxXChisq(0, 1000.0);
  dMuiRoadControl->set_maxYChisq(0, 1000.0);
  dMuiRoadControl->set_minSharedHits2D(0, 8);
  dMuiRoadControl->set_mvdvertex(0, 0);
  for (Int_t i = 0; i < 5; i++)
    {
      dMuiRoadControl->set_WeightPar1D(i, 0, 1.0);
      dMuiRoadControl->set_WeightPar2D(i, 0, 1.0);
      dMuiRoadControl->set_WindowParAlpha(i, 0, 0.0);
      dMuiRoadControl->set_WindowParBeta(i, 0, 0.0);
      dMuiRoadControl->set_WindowParGamma(i, 0, 15.0);
    }

  // yajun's new muid ghost/group parameters
  dMuiRoadControl->set_MutrWindow(0, 20.0);
  dMuiRoadControl->set_MuidWindow(0, 30.0);

  dMuiClusterControl->SetRowCount(1);
  dMuiClusterControl->set_xClusterSep(0, 1);
  dMuiClusterControl->set_yClusterSep(0, 1);
  
  for (Int_t i = 0;i < 5;i++)
    {
      dMuiClusterControl->set_xMaxClusterSize(i, 0, 1);
      dMuiClusterControl->set_yMaxClusterSize(i, 0, 1);
    }

  dMuiMuonIdProfilePar->SetRowCount(1);
  dMuiMuonIdProfilePar->set_pmax(0, 0, 0.0);
  dMuiMuonIdProfilePar->set_pmax(1, 0, 1.0);
  dMuiMuonIdProfilePar->set_pmax(2, 0, 1.3);
  dMuiMuonIdProfilePar->set_pmax(3, 0, 1.7);
  dMuiMuonIdProfilePar->set_pmax(4, 0, 200.0);
  dMuiMuonIdProfilePar->set_xVertexCut(0, 80.0);
  dMuiMuonIdProfilePar->set_yVertexCut(0, 80.0);

#ifdef DEBUG
  cout << "Finished setting MuID control parameters\n";
#endif

  // Monte Carlo, real data dependent variables for MuTr:

  if (simulationFlag)
    {
      dMutGlobalPar->set_MCData(0, 0); // 1 =use MC geometry; 0 = real one; default = 0 as of 04/10/2002

      mumdigipar->set_OfflineCathodeDCMScheme(0, 0);

      cout << "MUTR MC inputs: mutr_e_loss = " << mutr_e_loss 
	   << "\t mutr_e_base ="   << mutr_e_base << endl;

      //dedx = landau_distr*energy_loss+baseline (tur2/12 ~ 1/6)
      mumdigipar->set_LandauScale(0,0,mutr_e_loss); // default = 40.0 (energy_laoss) for all 3 stations
      mumdigipar->set_LandauScale(1,0,mutr_e_loss); // 0-2: south Staiotn-1,2,3
      mumdigipar->set_LandauScale(2,0,mutr_e_loss); // mutr_e_loss=7.0
      mumdigipar->set_LandauScale(3,0,mutr_e_loss);//3-5:north station-1-3
      mumdigipar->set_LandauScale(4,0,mutr_e_loss);
      mumdigipar->set_LandauScale(5,0,mutr_e_loss);
      
      mumdigipar->set_LandauOffset(0,0,mutr_e_base); // default baseline= 165 for all
      mumdigipar->set_LandauOffset(1,0,mutr_e_base); // mutr_e_base=28
      mumdigipar->set_LandauOffset(2,0,mutr_e_base);
      mumdigipar->set_LandauOffset(3,0,mutr_e_base);
      mumdigipar->set_LandauOffset(4,0,mutr_e_base);
      mumdigipar->set_LandauOffset(5,0,mutr_e_base);
      
      //      
      mumdigipar->set_resolution(0, 0, 0.02);
      mumdigipar->set_resolution(1, 0, 0.02);
      mumdigipar->set_resolution(2, 0, 0.02);
      mumdigipar->set_resolution(3, 0, 0.00);
      mumdigipar->set_resolution(4, 0, 0.00);

      dMutCathodePulseParams->setTrise(500.);
      dMutCathodePulseParams->setTfall(1500.);
      dMutCathodePulseParams->setFitRise(0);
      dMutCathodePulseParams->setFitFall(0);
      dMutCathodePulseParams->setFitType("exp");
    }
  else
    {
      dMutGlobalPar->set_MCData(0, 0);

      mumdigipar->set_OfflineCathodeDCMScheme(0, 0);

      mumdigipar->set_resolution(0, 0, 0.06);
      mumdigipar->set_resolution(1, 0, 0.06);
      mumdigipar->set_resolution(2, 0, 0.06);
      mumdigipar->set_resolution(3, 0, 0.06);
      mumdigipar->set_resolution(4, 0, 0.06);

      dMutCathodePulseParams->setTrise(4500.);
      dMutCathodePulseParams->setTfall(1000.);
      dMutCathodePulseParams->setFitRise(1);
      dMutCathodePulseParams->setFitFall(1);
      dMutCathodePulseParams->setFitType("avg");
    }

  dMuoFindTracksPar->SetRowCount(1);
  dMuoFindTracksPar->set_weed_chisq_imp(0, 0, 4.0);
  dMuoFindTracksPar->set_weed_chisq_imp(1, 0, 4.0);

  dMuoFindTracksPar->set_maskw_the(0, 0, 8.0);
  dMuoFindTracksPar->set_maskw_the(1, 0, 12.0);
  dMuoFindTracksPar->set_maskw_the(2, 0, 14.0);
  dMuoFindTracksPar->set_maskw_the(3, 0, 20.0);
  dMuoFindTracksPar->set_maskw_the(4, 0, 30.0);
  dMuoFindTracksPar->set_maskw_the(5, 0, 50.0);

  dMuoFindTracksPar->set_maskw_phi(0, 0, 10.0);
  dMuoFindTracksPar->set_maskw_phi(1, 0, 25.0);
  dMuoFindTracksPar->set_maskw_phi(2, 0, 40.0);
  dMuoFindTracksPar->set_maskw_phi(3, 0, 10.0);
  dMuoFindTracksPar->set_maskw_phi(4, 0, 25.0);
  dMuoFindTracksPar->set_maskw_phi(5, 0, 40.0);

  dMuoFindTracksPar->set_candtr_sigphi(0, 0, 0.0);
  dMuoFindTracksPar->set_candtr_sigphi(1, 0, 0.0);
  dMuoFindTracksPar->set_candtr_sigphi(2, 0, 2.0);
  dMuoFindTracksPar->set_candtr_sigphi(3, 0, 0.0);
  dMuoFindTracksPar->set_candtr_sigphi(4, 0, 0.0);
  dMuoFindTracksPar->set_candtr_sigphi(5, 0, 2.0);

  for (int i = 0; i < 9; i++)
    {
      dMuoFindTracksPar->set_projerr(i, 0, 0.6);
    }
  dMuoFindTracksPar->set_projerr(2, 0, 0.9);

  dMuoFindTracksPar->set_fqcut(0, 0.20);
  dMuoFindTracksPar->set_chifitcut(0, 0.0);
  dMuoFindTracksPar->set_qrescut(0, 0.15);
  dMuoFindTracksPar->set_awirecut(0, 9.00);
  dMuoFindTracksPar->set_muidmask(0, 1);
  //dMuoFindTracksPar->set_MinMuidDepth(0, 4);
  dMuoFindTracksPar->set_MinMuidDepth(0, 2);
  //dMuoFindTracksPar->set_MaxMuidChisqdf(0, 20.);
  dMuoFindTracksPar->set_MaxMuidChisqdf(0, 200.);
  dMuoFindTracksPar->set_MinMuidNhits(0, 5);
  dMuoFindTracksPar->set_MaxMuidPlaHits(0, 10);
  dMuoFindTracksPar->set_MinMuidTheta(0, 0.); //0 deg
  //dMuoFindTracksPar->set_MinMuidTheta(0, 0.209); //12 deg
  //dMuoFindTracksPar->set_MinMuidTheta(0, 0.349); //20 deg
  dMuoFindTracksPar->set_muidfit(0, 1);

  dMuoFindTracksPar->set_RealVertErr(0, 1.0);
  dMuoFindTracksPar->set_NoVertErr(0, 25.0);
  dMuoFindTracksPar->set_MCVertErr(0, 1.0);
  dMuoFindTracksPar->set_NoVertMCErr(0, 25.0);


  dMuoFindTracksPar->set_mvdvert(0, 0);
  dMuoFindTracksPar->set_bbcvert(0, 0);
  dMuoFindTracksPar->set_globalvert(0, 1); //use global vtx

  dMuoFindTracksPar->set_prfit(0, 1);
  dMuoFindTracksPar->set_min_stub_hits(0, 2);
  dMuoFindTracksPar->set_min_track23_hits(0, 6);
  dMuoFindTracksPar->set_min_track_hits(0, 10);
  dMuoFindTracksPar->set_nplskip(0, 2);

  dMuoFindTracksPar->set_ghost_nshare_min(0, 0, 0);
  dMuoFindTracksPar->set_ghost_nshare_min(1, 0, 5);
  dMuoFindTracksPar->set_ghost_nshare_min(2, 0, 0);

  dMuoFindTracksPar->set_max_good_hits(0, 4);

  // set sign of bend plane kicks (revbpkick)
  if(getMapFileFlag()==0) {
    dMuoFindTracksPar->set_revbpkick(0,1);
  } else {
    if(getMapFileScale() == 0) {
      dMuoFindTracksPar->set_revbpkick(0,0);
    } else if(getMapFileScale() == -1) {
      dMuoFindTracksPar->set_revbpkick(0,1);
    } else {
      dMuoFindTracksPar->set_revbpkick(0,-1);
    }
  }

  dMuoFindTracksPar->set_max_stub_chisq(0, 0, 100.0);
  dMuoFindTracksPar->set_max_stub_chisq(1, 0, 100.0);
  dMuoFindTracksPar->set_max_stub_chisq(2, 0, 50.0);
  dMuoFindTracksPar->set_max_stub_chisq(3, 0, 100.0);
  dMuoFindTracksPar->set_max_stub_chisq(4, 0, 100.0);
  dMuoFindTracksPar->set_max_stub_chisq(5, 0, 50.0);

  dMuoFindTracksPar->set_stub_fit_dwh(0, 0, 0.0);
  dMuoFindTracksPar->set_stub_fit_dwh(1, 0, 0.0);
  dMuoFindTracksPar->set_stub_fit_dwh(2, 0, 0.0);
  dMuoFindTracksPar->set_stub_fit_dwh(3, 0, 0.0);
  dMuoFindTracksPar->set_stub_fit_dwh(4, 0, 0.0);
  dMuoFindTracksPar->set_stub_fit_dwh(5, 0, 0.0);

  dMuoFindTracksPar->set_max_track23_chisq(0, 5000.0);
  dMuoFindTracksPar->set_max_track_chisq(0, 150.0);
  dMuoFindTracksPar->set_mscat_z(0, 85.0);
  dMuoFindTracksPar->set_max_geane_chisqr(0, 10000.0);

  dMuoFindTracksPar->set_ghost_qual_fac(0, 0, 2.0);
  dMuoFindTracksPar->set_ghost_qual_fac(1, 0, 6.0);
  dMuoFindTracksPar->set_ghost_qual_fac(2, 0, 7.5);

  dMuoFindTracksPar->set_ghost_qual_subset_fac(0, 0, 1.2);
  dMuoFindTracksPar->set_ghost_qual_subset_fac(1, 0, 1.2);
  dMuoFindTracksPar->set_ghost_qual_subset_fac(2, 0, 0.2);

  for (int i = 0; i < 3; i++)
    {
      dMuoFindTracksPar->set_ghost_chi_min(i, 0, 0.0);
    }

  for (int i = 0; i < 2; i++)
    {
      dMuoFindTracksPar->set_weed_chisq_bad(i, 0, 10000.0);
      dMuoFindTracksPar->set_weed_chisq_max(i, 0, 0.5);
    }

  dMuoFindTracksPar->set_max_xproj(0, 0, 300.0);
  dMuoFindTracksPar->set_max_yproj(0, 0, 300.0);
  dMuoFindTracksPar->set_max_xproj(1, 0, 250.0);
  dMuoFindTracksPar->set_max_yproj(1, 0, 250.0);
  dMuoFindTracksPar->set_max_xproj(2, 0, 150.0);
  dMuoFindTracksPar->set_max_yproj(2, 0, 150.0);
  dMuoFindTracksPar->set_max_xproj(3, 0, 300.0);
  dMuoFindTracksPar->set_max_yproj(3, 0, 300.0);
  dMuoFindTracksPar->set_max_xproj(4, 0, 250.0);
  dMuoFindTracksPar->set_max_yproj(4, 0, 250.0);
  dMuoFindTracksPar->set_max_xproj(5, 0, 150.0);
  dMuoFindTracksPar->set_max_yproj(5, 0, 150.0);

  dMuoFindTracksPar->set_max_xproj_tr(0, 0, 50.0);
  dMuoFindTracksPar->set_max_yproj_tr(0, 0, 50.0);
  dMuoFindTracksPar->set_max_xproj_tr(1, 0, 200.0);
  dMuoFindTracksPar->set_max_yproj_tr(1, 0, 200.0);
  dMuoFindTracksPar->set_max_xproj_tr(2, 0, 100.0);
  dMuoFindTracksPar->set_max_yproj_tr(2, 0, 100.0);
  dMuoFindTracksPar->set_max_xproj_tr(3, 0, 50.0);
  dMuoFindTracksPar->set_max_yproj_tr(3, 0, 50.0);
  dMuoFindTracksPar->set_max_xproj_tr(4, 0, 200.0);
  dMuoFindTracksPar->set_max_yproj_tr(4, 0, 200.0);
  dMuoFindTracksPar->set_max_xproj_tr(5, 0, 100.0);
  dMuoFindTracksPar->set_max_yproj_tr(5, 0, 100.0);

  dMuoFindTracksPar->set_awirecutchi(0, 0.1);
  dMuoFindTracksPar->set_StubFitVert(0, 0);
  //dMuoFindTracksPar->set_StubFitVert(0, 1);

  dMutGlobalPar->SetRowCount(1);
  dMutGlobalPar->set_MCAnalysis(0, 0);

  dMutHistos->SetRowCount(1);
  dMutHistos->set_Dimuons(0, 1); //all
  dMutHistos->set_Thrown(0, 0);
  dMutHistos->set_Recon(0, 1);
  dMutHistos->set_BPFit(0, 1); // Mike
  dMutHistos->set_Stub(0, 1);  //Yuki
  dMutHistos->set_NCandHits(0, 0); //
  dMutHistos->set_CandTrs(0, 0);  //
  dMutHistos->set_CheckHit(0, 0); //
  dMutHistos->set_CheckMask(0, 1); // Mike
  dMutHistos->set_RawData(0, 0); //
  dMutHistos->set_StripData(0, 1); // Ming
  dMutHistos->set_ClustData(0, 1); //Yuki
  dMutHistos->set_ClustData2(0, 1);

  // Change some flags for evaluation purposes
  if (evaluationFlag) {
    dMuoFindTracksPar->set_globalvert(0, 0);
    dMuoFindTracksPar->set_mcvert(0, 1);
    dMutGlobalPar->set_MCAnalysis(0, 1);
    dMutHistos->set_Thrown(0, 1);
    dMutHistos->set_Recon(0, 1);
    dMutHistos->set_ClustData2(0, 1);
  }

  mumdigipar->SetRowCount(1);

  mumdigipar->set_hit_efficiency(0, 1.0);

  for (int i = 0; i < 3; i++) {
    mumdigipar->set_st_angacc(i, 0, 10.0);
    mumdigipar->set_st_angacc(i + 3, 0, 12.0);
  }

  mumdigipar->set_baselor(0, 0.85);

  for (int i = 0; i < 9; i++) {
    mumdigipar->set_keepgap(i, 0, 1);
  }
  mumdigipar->set_keepgap(8, 0, 0);

  for (int i = 0; i < 10; i++) {
    mumdigipar->set_threshold(i, 0, 6.0);
    mumdigipar->set_noise_level(i, 0, 2.0);
    mumdigipar->set_noise_scale(i, 0, 1.0);
  }
  for (int i = 5; i < 10; i++) {
    mumdigipar->set_noise_scale(i, 0, 0.0);
  }
  mumdigipar->set_threshold(9, 0, 0.0);
  mumdigipar->set_noise_level(9, 0, 0.0);

  mumdigipar->set_GainMin(0, 1.0);

  for (int i = 0; i < 3; i++) {
    mumdigipar->set_wire_spacing(i, 0, 0.5);
  }
  for (int i = 3; i < 5; i++) {
    mumdigipar->set_wire_spacing(i, 0, 0.0);
  }

  mumdigipar->set_CathodeDCMScheme(0, 1111);

  mumdigipar->set_ThreshFactor(0, 3.0);
  mumdigipar->set_ThreshFactor2(0, -100.0);
  mumdigipar->set_MinClusCharge(0, 2.0);

  dMutCathodePulseParams->setNsamples(4);
  float tsamples[4];
  tsamples[0] = 100.;
  tsamples[1] = 600.;
  tsamples[2] = 700.;
  tsamples[3] = 800.;
  dMutCathodePulseParams->setTsamples(tsamples);

  dMfmPar->SetRowCount(1);
  dMfmPar->set_verbosity(0, 0);
  dMfmPar->set_mapFile(0, getMapFileFlag() );
  mMuiDCMUnpack = new mMuiDCMUnpackModule;
  mMuiClusterFinder = new mMuiClusterFinderModule;
  mMuiFillRoadRawRel = new mMuiFillRoadRawRelModule;

  if (evaluationFlag) {
    mMuiFillRawFkinRel = new mMuiFillRawFkinRelModule;
    mMuiFillClusterFkinRel = new mMuiFillClusterFkinRelModule;
    mMuiFillRoadFkinRel = new mMuiFillRoadFkinRelModule; 
    mMuiEvaluator = new mMuiEval;
    mMutFillRawCathGhitRel= new mMutFillRawCathGhitRelModule;
    mMutFillClusterGhitRel = new mMutFillClusterGhitRelModule;
    mMutCathFitEvaluator = new mMutCathFitEval;
  }

  mMuiChannelGhitRel = new mMuiChannelGhitRelModule;
  mMuiFEMToDCM = new mMuiFEMToDCMModule;
  mMuiRawToFEM = new mMuiRawToFEMModule;
  mMuiReadout = new mMuiReadoutModule;
  mMuiReadTubeEff = new mMuiReadTubeEffModule;
  mMuiRoadFind = new mMuiRoadFinder;
  mMuiRoadFind2 = new mMuiRoadFinder2;
  mMuiHistoBase = new mMuiHistoBaseModule;
  mMuiInit = new mMuiInitModule;
  mMuiEnd = new mMuiEndModule;
  mMuiMuonId = new mMuiMuonIdModule;
  mMuiPseudoTrig = new mMuiPseudoTrigger;

  mMfmMT = new mMfmMTModule;

  mMutdbInitM = new mMutdbInit;
  mMutRootInit = new mMutRootInitModule;
  mumUser = new mum_user;
  mMutXYtoFEEM = new mMutXYtoFEE;
  mMutZeroSup = new mMutZeroSupModule;
  mMutZeroSuptoRawDataM = new mMutZeroSuptoRawData;
  mMutApplyCalibM = new mMutApplyCalib;
  mMutFindCathodeClustsTM = new mMutFindCathodeClustsT;
  mMutFitCathodeClustsT = new mMutFitCathodeClustsTModule;
  mMuoFindTracksM = new mMuoFindTracks;
  mMuoFillDataMineVars = new mMuoFillDataMineVarsModule;
  mMutEnd = new mMutEndModule;
  mMutCloseRootFile = new mMutCloseRootFileModule;

#ifdef DEBUG
  cout << "Finished instantiating analysis modules\n";
#endif

  //Get gloabal TimeStamp from preco
#ifdef DEBUG
  cout << "PRECO Global TimeStamp = " << TimeStamp << endl;
#endif

  PHTimeStamp startTime = TimeStamp;

  if (simulationFlag)
    {
      PHTimeStamp t_search(2000, 4, 2, 0, 0, 0);  // April 2, 2000, midnight.
      mMuiInit->SetSearchTimeStamp(t_search);

      if( simulationFlag == 3 ) {
#ifdef DEBUG
	cout << " calling mMuiReadTubeEff " << endl;
#endif
	mMuiReadTubeEff->event(topNode);
      }

      mMuiPseudoTrig->MapSelect(2);
    }
  else
    {
#ifdef DEBUG
      cout << "--->MUID: Gloabal startTime for run " << runNumber 
	   << " is:  " << startTime << endl;
#endif

      mMuiInit->SetSearchTimeStamp(startTime);
      mMuiPseudoTrig->MapSelect(0);
    }
  mMuiInit->event(topNode);

#ifdef DEBUG
  cout << "murec: Calling mMfmMT\n";
#endif

  mMfmMT->event(topNode);

#ifdef DEBUG
  cout << "murec: Calling mMutdbInit\n";
#endif
  if (simulationFlag)
    { 
      cout << "Using real calib data for MC: --->MUTR: Global startTime for run " << runNumber 
	   << " is:  " << startTime << endl;
 
      if (runNumber>0) { // to use real calib data information for better mutr response 
	mMutdbInitM->dbGetAll(topNode, startTime);
      }
      else {
	PHTimeStamp tsearch(2001, 5, 4, 0, 0, 0);
	mMutdbInitM->dbGetAll(topNode, tsearch);
      }
    }
  else
    {
#ifdef DEBUG
      cout << "--->MUTR: Global startTime for run " << runNumber 
	   << " is:  " << startTime << endl;
#endif
      mMutdbInitM->dbGetAll(topNode, startTime);
    }

#ifdef DEBUG
  cout << "murec: Calling mMutRootInit\n";
#endif
  mMutRootInit->event(topNode);

#ifdef DEBUG
  cout << "murec: Calling mumUser\n";
#endif
  mumUser->event(topNode);

  return;
}

void 
ana_mui()
{

  if(simulationFlag == 2 || simulationFlag == 3) {
    //
    // insertions here for the event-by-event slow simulator code
    // which initially processes the PISA hits information
    //

  } // check on PISA-to-DST or PISA-to-PRDF mode

  if(simulationFlag == 3) {
    //
    // insertions here for the event-by-event FEM and DCM modules
    //
      mMuiReadout->event(topNode);  
      mMuiRawToFEM->event(topNode);
      mMuiFEMToDCM->event(topNode);
      mMuiChannelGhitRel->event(topNode);
      MuiPutDCM(topNode);        // put DCM table in PRDF node and On-Line buffers

      mMutXYtoFEEM->event(topNode);  //
      mMutZeroSup->event(topNode);
      MutPutDCM(topNode);

      return;  // simulationFlag = 3 means no event reconstruction

  } // check on PISA-to-PRDF mode

  if(simulationFlag == 2) {
    //
    // insertions here for any special PISA-to-DST event-by-event modules
    // needed for reconstruction not from a simulated PRDF
    //
    cout << "ERROR:  simulationFlag=2 not yet implemented for Muons" << endl;

  } // check on PISA-to-DST mode

  if (evaluationFlag) {

    ioEval->read(evaNode);
    
    pisarun->GetOneEvent(pisaEventArray, &kPisaEvent, TTArray, nFiles);

    KinGetGEA(topNode);  // fill the fkin table
    MutGetGEA(topNode);  // fill the mumhits table
    MuiGetGEA(topNode);  // fill the munhits table
  }

#ifdef DEBUG
  cout << "Calling mMuiPseudoTrigger\n";
#endif
  mMuiPseudoTrig->event(topNode);

#ifdef DEBUG
  cout << "Calling MuiGetDCM\n";
#endif
  MuiGetDCM(topNode);

#ifdef DEBUG
  cout << "Calling mMuiDCMUnpack\n";
#endif
  mMuiDCMUnpack->event(topNode);

#ifdef DEBUG
  cout << "Calling mMuiClusterFinder\n";
#endif
  mMuiClusterFinder->event(topNode);

#ifdef DEBUG
  cout << "Calling mMuiRoadFinder\n";
#endif
  mMuiRoadFind->event(topNode);

#ifdef DEBUG
  cout << "Calling mMuiFillRoadRawRel\n";
#endif
  mMuiFillRoadRawRel->event(topNode);

  if (evaluationFlag) {
#ifdef DEBUG
    cout << "Calling mMuiFillRawFkinRel\n";
#endif 
    mMuiFillRawFkinRel->event(topNode);
    
#ifdef DEBUG
    cout << "Calling mMuiFillClusterFkinRel\n";
#endif
    mMuiFillClusterFkinRel->event(topNode);

#ifdef DEBUG
    cout << "Calling mMuiFillRoadFkinRel\n";
#endif
    mMuiFillRoadFkinRel->event(topNode);
  }

#ifdef DEBUG
  cout << "Calling mMuiHistoBase\n";
#endif
  mMuiHistoBase->event(topNode);

  if (precoSelect::isMUTR())
    {
      
#ifdef DEBUG
      cout << "murec: Calling mMutZeroSuptoRawData\n";
#endif
      mMutZeroSuptoRawDataM->event(topNode);
      
      if (evaluationFlag) {
#ifdef DEBUG
	cout << "murec: Calling mMutFillRawCathGhitRel\n";
#endif
	mMutFillRawCathGhitRel->event(topNode);
      }

#ifdef DEBUG
      cout << "murec: Calling mMutApplyCalib\n";
#endif
      mMutApplyCalibM->event(topNode);

#ifdef DEBUG
      cout << "murec: Calling mMutFindCathodeClustsT\n";
#endif
      mMutFindCathodeClustsTM->event(topNode);

      if (evaluationFlag) {
#ifdef DEBUG
	cout << "murec: Calling mMutFillClusterGhitRel\n";
#endif
	mMutFillClusterGhitRel->event(topNode);
      }

#ifdef DEBUG
      cout << "murec: Calling mMutFitCathodeClustsT\n";
#endif
      mMutFitCathodeClustsT->event(topNode);

	if (evaluationFlag) {
#ifdef DEBUG
	  cout << "murec: Calling mMutCathFitEval\n";
#endif
	  mMutCathFitEvaluator->event(topNode);
	}

#ifdef DEBUG
      cout << "murec: Calling mMuoFindTracks\n";
#endif
      mMuoFindTracksM->event(topNode);

#ifdef DEBUG
      cout << "murec: Calling mMuiRoadFinder2\n";
#endif
      mMuiRoadFind2->event(topNode);

#ifdef DEBUG
      cout << "murec: Calling mMuiMuonId\n";
#endif
      mMuiMuonId->event(topNode);

#ifdef DEBUG
      cout << "murec: Calling mMuoFillDataMineVars\n";
#endif
      mMuoFillDataMineVars->event(topNode);
    }

  if (evaluationFlag) {
  // Clean up PISA memory use at end of full event
    pisarun->HitsClear(); // releases memory assigned to XxxPISAHit globals
  }

  return;
}

void 
end_mui()
{
#ifdef DEBUG
  cout << "murec: Calling mMuiEnd\n";
#endif
  mMuiEnd->event(topNode);

#ifdef DEBUG
  cout << "murec: Calling mMutEnd\n";
#endif
  mMutEnd->event(topNode);

#ifdef DEBUG
  cout << "murec: Calling mMutCloseRootFile\n";
#endif
  mMutCloseRootFile->event(topNode);

}