#include "dreco.h"
#include "PHNodeIterator.h"
#include "RunToTimeObjy.hh"
#include "dMuiRoadsWrapper.h"
#include "dMuiRoadRawRelWrapper.h"
#include "dMutCalibCathodesWrapper.h"
#include "PHIODataNode.h"
#include "PHDataNode.h"
#include "PhHistogramFactory.hh"
#include "PHVector.h"
#include "Lvl2Event.h"
#include "Lvl2DecisionOutv1.h"
#include "EventHeaderv1.h"
#include "dMuoTracksWrapper.h"
#include "dMuoTrackRoadRelWrapper.h"
#include "dMuiRawWrapper.h"

// Display
#include "PhenixRun.hh"
#include "PhEventDisplay.hh"
#include "PhMuiDisplay.hh"
#include "PhCanvas.hh"
#include "mMuiInitModule.h"
// Level2 includes
#include "L2MuiTracks.h"  //This Must come after event display includes
#include "L2BbcZdcCentrality.h"
#include "dMuiClustersWrapper.h"
#include "dMuiRoadClusterRelWrapper.h"
// MicroDst
#include "TFile.h"
#include "TNtuple.h"

// Refit
#include "dMuiClusterControlWrapper.h"
#include "dMuiRoadControlWrapper.h"
#include "mMuiClusterFinderModule.h"
#include "mMuiRoadFinder.hh"
#include "mMuiRoadFinder2.hh"
#include "mMuiFillRoadRawRelModule.h"
#include "mMuiFillRawFkinRelModule.h"
#include "mMuiFillClusterFkinRelModule.h"
#include "mMuiFillRoadFkinRelModule.h"
#include "mfm_setfscale.h"
#include "mfm_getfscale.h"

#include "mMutFindCathodeClustsT.hh"
#include "mMutFitCathodeClustsTModule.h"
#include "mMuoFindTracks.hh"
#include "mMutCathodePulseParams.hh"
#include "mum_user.hh"
#include "mMutRootInitModule.h"
#include "mMfmMTModule.h"
#include "mMutdbInit.hh"
#include "dMuoFindTracksParWrapper.h"
#include "dMutHistosWrapper.h"
#include "dMutGlobalParWrapper.h"
#include "dMfmParWrapper.h"
#include "dMutCathodeClustersWrapper.h"
#include "mumdigiparWrapper.h"


#include "mMuiEndModule.h"
#include "mMutEndModule.h"
#include "mMutCloseRootFileModule.h"

#include "dMutCalibCathodesOut.h"

EventHeader* eventheader = NULL;

// Event Display
bool first_event = true;
bool use_embedding = false;

// Refit

mMuiInitModule* mMuiInit;
mMuiClusterFinderModule* mMuiClusterFinder;
mMuiRoadFinder* mMuiRoadFind;
mMuiRoadFinder2* mMuiRoadFind2;
mMuiFillRoadRawRelModule* mMuiFillRoadRawRel;

mMutFindCathodeClustsT* mMutFindCathodeClustsTM;
mMutFitCathodeClustsTModule* mMutFitCathodeClustsTM;
mMuoFindTracks* mMuoFindTracksM;
mMfmMTModule* mMfmMT;
mMutdbInit* mMutdbInitM;
mMutRootInitModule* mMutRootInit;
mum_user* mumUser;
mMuiEndModule* mMuiEnd;
mMutEndModule* mMutEnd;
mMutCloseRootFileModule* mMutCloseRootFile;

PHCompositeNode* motherNode;


PHNodeIOManager *dstOut = 0;

//Embedding variables
PHNodeIOManager *dst2IOManager;
PHCompositeNode* dst2Node;
int mergeMui(PHCompositeNode* mainNode, PHCompositeNode* embedNode);
int mergeMut(PHCompositeNode* mainNode, PHCompositeNode* embedNode);

// MicroDst
int nevt;
TFile* microdst;
TNtuple* Dimuons;
TNtuple* SingleMus;

int RefitFlag = 1;
bool simulationFlag = false;

int dinit(){
  cout<<endl<<"This analysis code refits muon tracks from the DST."<<endl;
}

void
setRefitFlag(int Flag=1)
{
  RefitFlag = Flag;
  cout << "\n setRefitFlag:  RefitFlag = " << RefitFlag << endl;
}

void
setEmbeddingFlag(bool Flag=true)
{
  use_embedding = Flag;
  cout << "\n Embed a hits from a secondary DST:  use_embedding = " << use_embedding << endl;
  cout << " Embedding from file dstin2.root" << endl;
}

void
setSimulationFlag(bool Flag=true)
{
  simulationFlag = Flag;
  cout << "\n setSimulationFlag:  simulationFlag = " << simulationFlag << endl;
  cout << "\n Reconstruction will assume primary dst to be simulated!\n";
}


int
dstOutFileOpen(PHString dstOutFile)
{
  dstOut = new PHNodeIOManager(dstOutFile, PHWrite);
  dstOut->SetCompressionLevel(3);
  return 0;
}


int close_all()
{

  if (RefitFlag){
    cout << "Calling close_all routines" << endl;

    mMuiEnd->event(motherNode);
    mMutEnd->event(motherNode);
    mMutCloseRootFile->event(motherNode);
  }

  if(microdst>0)
  {
    cout << "Closing MicroDst File" << endl;
    microdst->Write();
    microdst->Close();
  }

  delete dstOut;
}

int setup_all()
{
//  int deb;
//  cout << "gdb pause" << endl;
//  cin >> deb;

  nevt = 0;
  return 0;
}

int process_event(PHCompositeNode* topNode)
{
  nevt++;

  if(first_event){
//    setup_display(topNode);

  if (RefitFlag){         //Refind and fit muon tracks
    setup_refit(topNode);
  }

// Print once the dst contents
    topNode->print();
    first_event = false;
  }
// user Analysis Modules

  if (RefitFlag){         //Refind and fit muon tracks
    ana_refit(motherNode);
  }

  if (RefitFlag && dstOut){   //Write out new dst information
    dstOut->write(topNode);

  }
  return 0;
}
//-------------------------------------------------------------------- REFIT CODE
int setup_refit(PHCompositeNode* topNode)
{
// Create a mother node to fake the reconstruction structure
  motherNode = new PHCompositeNode("TOP");
  PHNodeIterator* mainIter = new PHNodeIterator(motherNode);
  PHCompositeNode* muiNode = new PHCompositeNode("MUI");
  motherNode->addNode(muiNode);
  PHCompositeNode* mutNode = new PHCompositeNode("MUT");
  motherNode->addNode(mutNode);
//  PHCompositeNode* prdfNode = new PHCompositeNode("PRDF");
//  motherNode->addNode(prdfNode);
  PHCompositeNode* evaNode = new PHCompositeNode("EVA");
  motherNode->addNode(evaNode);
  PHCompositeNode* parNode = new PHCompositeNode("PAR");
  motherNode->addNode(parNode);
// Add the topNode aka the dst node to the mother node
  motherNode->addNode(topNode);
  //Another DST node for embedding
  //Should not be added to the motherNode
  dst2Node = new PHCompositeNode("DST");


// Specify a new dst file:

PHNodeIOManager *dstOut = 0;

// Read the parameter file
  PHNodeIOManager * parIn = new PHNodeIOManager("mupar.root", PHReadOnly);
  parIn->read(parNode);

  PHIODataNode<TObject>* dn;

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

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

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

  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);
    }

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

  PHIODataNode<PHTable>* dMutCathodePulseParamsNode = (PHIODataNode<PHTable>*)mainIter->findFirst("PHIODataNode", "dMutCathodePulseParams");
  mMutCathodePulseParams* dMutCathodePulseParams = (mMutCathodePulseParams*)dMutCathodePulseParamsNode->getData();

  PHIODataNode<PHTable>* dMfmParNode = (PHIODataNode<PHTable>*)mainIter->findFirst("PHIODataNode", "dMfmPar");
  dMfmParWrapper* dMfmPar = (dMfmParWrapper*)dMfmParNode->getData();

// mumdigipar should be in my mupar.root
  PHIODataNode<PHTable>* MumDigiParNode = (PHIODataNode<PHTable>*)mainIter->findFirst("PHIODataNode", "mumdigipar");
  mumdigiparWrapper* mumdigipar;
  if (MumDigiParNode)
    {
      mumdigipar = (mumdigiparWrapper*)MumDigiParNode->getData();
    } else {
      cout << __FILE__ << ":" << __LINE__
	   << ": MumDigiParNode not found!" << endl;
      mumdigipar = 0;
    }

  PHIODataNode<EventHeader>::iterator eventheader(topNode);

  if(eventheader.find("EventHeader")==0)
  { cout<<"EventHeader not found.\n";
  //return -1;
  }
  

  
// MuID input parameters:

  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);
    }


  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);


// Monte Carlo, real data dependent variables for MuTr
  if (simulationFlag)
    {
      dMutGlobalPar->set_MCData(0, 1);
      mumdigipar->set_OfflineCathodeDCMScheme(0, 0);
      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);
      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, 2);
// Latest parameters
  dMuoFindTracksPar->set_MaxMuidChisqdf(0, 200);
  dMuoFindTracksPar->set_MinMuidNhits(0,6);
  dMuoFindTracksPar->set_MaxMuidPlaHits(0, 10);
  dMuoFindTracksPar->set_MinMuidTheta(0,0);

  dMuoFindTracksPar->set_muidfit(0, 1);
  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);
  // dMuoFindTracksPar->set_revbpkick(0,0);
  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);

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

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

  mumdigipar->SetRowCount(1);
  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);

  mumdigipar->set_hit_efficiency(0, 1.0);
  mumdigipar->set_noise_fraction(0, 0.0);
  mumdigipar->set_GainMin(0, 1.0);

  for (int i = 0; i < 3; i++)
    {
      mumdigipar->set_wire_spacing(i, 0, 0.5);
      mumdigipar->set_wire_anode_spacing(i, 0, 1.0);
      mumdigipar->set_wire_coarse_cath_spacing(0, 0, 1.0);
    }

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

  mumdigipar->set_new_stereo_angle1(0, 0, 22.5);
  mumdigipar->set_new_stereo_angle2(0, 0, 22.5);
  mumdigipar->set_new_stereo_angle3(0, 0, 22.5);
  mumdigipar->set_new_stereo_angle1(1, 0, 28.125);
  mumdigipar->set_new_stereo_angle2(1, 0, 26.25);
  mumdigipar->set_new_stereo_angle3(1, 0, 11.25);
  mumdigipar->set_new_stereo_angle1(2, 0, 22.5);
  mumdigipar->set_new_stereo_angle2(2, 0, 22.5);
  mumdigipar->set_new_stereo_angle3(2, 0, 22.5);
  mumdigipar->set_new_stereo_angle1(3, 0, 11.25);
  mumdigipar->set_new_stereo_angle2(3, 0, 30.00);
  mumdigipar->set_new_stereo_angle3(3, 0, 11.25);
  mumdigipar->set_new_stereo_angle1(4, 0, 22.5);
  mumdigipar->set_new_stereo_angle2(4, 0, 22.5);
  mumdigipar->set_new_stereo_angle3(4, 0, 22.5);
  mumdigipar->set_new_stereo_angle1(5, 0, 33.75);
  mumdigipar->set_new_stereo_angle2(5, 0, 33.75);
  mumdigipar->set_new_stereo_angle3(5, 0, 11.25);

  for (int i = 6; i < 10; i++)
    {
      mumdigipar->set_new_stereo_angle1(i, 0, 0.0);
      mumdigipar->set_new_stereo_angle2(i, 0, 0.0);
      mumdigipar->set_new_stereo_angle3(i, 0, 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);

  for (int i = 0; i < 5; i++)
    {
      mumdigipar->set_zmumoff(i, 0, 0.0);
      mumdigipar->set_xmumoff(i, 0, 0.0);
      mumdigipar->set_ymumoff(i, 0, 0.0);
    }

  mumdigipar->set_bzoff(0, 0.0);
  mumdigipar->set_bxoff(0, 0.0);
  mumdigipar->set_byoff(0, 0.0);
  mumdigipar->set_gridrsize(0, 1.0);
  mumdigipar->set_gridzsize(0, 1.0);
  mumdigipar->set_berr(0, 0.0);
  mumdigipar->set_rnfmin(0, 0.0);
  mumdigipar->set_rnfmax(0, 0.0);
  mumdigipar->set_znfmin(0, 0.0);
  mumdigipar->set_znfmax(0, 0.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);
  mumdigipar->set_chamber_type(0, 3);

  for (int i = 0; i < 6; i++)
    {
      mumdigipar->set_mt_coords_per_plane(i, 0, 2);
    }
  for (int i = 6; i < 10; i++)
    {
      mumdigipar->set_mt_coords_per_plane(i, 0, 0);
    }

  mumdigipar->set_tdcchans(0, 1024);
  mumdigipar->set_minhits(0, 10);
  mumdigipar->set_lgrid(0, 0);

  for (int i = 0; i < 6; i++)
    {
      mumdigipar->set_signal_type1(i, 0, 1);
      mumdigipar->set_signal_type2(i, 0, 1);
      mumdigipar->set_signal_type3(i, 0, 1);
    }
  for (int i = 6; i < 10; i++)
    {
      mumdigipar->set_signal_type1(i, 0, 0);
      mumdigipar->set_signal_type2(i, 0, 0);
      mumdigipar->set_signal_type3(i, 0, 0);
    }

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

  mumdigipar->set_drawevent(0, 0);
  mumdigipar->set_drawmuid(0, 0);
  mumdigipar->set_drawghit(0, 0);
  mumdigipar->set_drawchamber(0, 0);
  mumdigipar->set_drawarm(0, 0);

  for (int i = 0; i < 27; i++)
    {
      mumdigipar->set_drawplane(i, 0, 0);
    }
  mumdigipar->set_drawplane(0, 0, 1);
  mumdigipar->set_drawplane(1, 0, 1);
  for (int i = 0; i < 8; i++)
    {
      mumdigipar->set_drawoctant(i, 0, 1);
    }

  mumdigipar->set_PerfectCathFit(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);

  mumdigipar->set_cathode_strip_type(0, 0, 's');
  mumdigipar->set_cathode_strip_type(1, 0, 't');
  mumdigipar->set_cathode_strip_type(2, 0, 'r');
  mumdigipar->set_cathode_strip_type(3, 0, 'a');
  mumdigipar->set_cathode_strip_type(4, 0, 'i');
  mumdigipar->set_cathode_strip_type(5, 0, 'g');
  mumdigipar->set_cathode_strip_type(6, 0, 'h');
  mumdigipar->set_cathode_strip_type(7, 0, 't');

  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, 2);

  mMuiInit = new mMuiInitModule;
  mMuiClusterFinder = new mMuiClusterFinderModule;
  mMuiRoadFind = new mMuiRoadFinder;
  mMuiFillRoadRawRel = new mMuiFillRoadRawRelModule;
  mMuiRoadFind2 = new mMuiRoadFinder2;

  mumUser = new mum_user;
  mMfmMT = new mMfmMTModule;
  mMutdbInitM = new mMutdbInit;
  mMutRootInit = new mMutRootInitModule;
  mMutFindCathodeClustsTM = new mMutFindCathodeClustsT;
  mMutFitCathodeClustsTM = new mMutFitCathodeClustsTModule;
  mMuoFindTracksM = new mMuoFindTracks;
  mMuiEnd = new mMuiEndModule;
  mMutEnd = new mMutEndModule;
  mMutCloseRootFile = new mMutCloseRootFileModule;
//Magnetic field should come from rundb
  setMapFileScale(1.0);
  getMapFileScale();
  
  mMfmMT->event(motherNode);

// The time stamps should be taken from the dst !
// ... unless it is a simulated dst!
  PHTimeStamp startTime;

  if(simulationFlag){
    startTime.set(2001, 5, 4, 0, 0, 0);
  }else
  {
    time_t EventTime = (*eventheader).get_TimeStamp();
    startTime = EventTime;
  }
  cout << "startTime = " << startTime << endl;

  mMuiInit->SetSearchTimeStamp(startTime);
  mMuiInit->event(motherNode);


  mMutdbInitM->dbGetAll(motherNode, startTime);
  mMutRootInit->event(motherNode);
  mumUser->event(motherNode);
  if(use_embedding)
  {
    // Open the geant evaluation file in an elegant way
    dst2IOManager = new PHNodeIOManager("dstin2.root", PHReadOnly);
    if (dst2IOManager==0) {
      cout << "File 2nd Dst for embedding needed for embedding not exist" << endl;
      cout << "Embedding will not be performed!!!\n";
      use_embedding = false;
    }
  }

}

int ana_refit(PHCompositeNode *topNode)
{
  int TrackonDst = 0;
  int TrackRefit = 0;

  PHIODataNode<dMuoTracksWrapper>::iterator dMuoTracks(topNode);

  if(dMuoTracks.find("dMuoTracks")==0)
    { cout<<"dMuoTracks not found.\n";
    return -1;
    }

  PHIODataNode<dMuoTrackRoadRelWrapper>::iterator dMuoTrackRoadRel(topNode);

  if(dMuoTrackRoadRel.find("dMuoTrackRoadRel")==0)
    { cout<<"dMuoTrackRoadRel not found.\n";
    return -1;
    }

  PHIODataNode<dMuiClustersWrapper>::iterator dMuiClustersW(topNode);

  if(dMuiClustersW.find("dMuiClusters")==0)
    { cout<<"dMuiClusters not found.\n";
    return -1;
    }

  PHIODataNode<dMuiRoadsWrapper>::iterator dMuiRoads(topNode);

  if(dMuiRoads.find("dMuiRoads")==0)
    { cout<<"dMuiRoads not found.\n";
    return -1;
    }

  PHIODataNode<dMuiRoadRawRelWrapper>::iterator dMuiRoadRawRel(topNode);

  if(dMuiRoadRawRel.find("dMuiRoadRawRel")==0)
    { cout<<"dMuiRoadRawRel not found.\n";
    return -1;
    }

  if(use_embedding)
  {

    //Read in the file
    dst2IOManager->read(dst2Node);

    PHIODataNode<dMutCathodeClustersWrapper>::iterator dMutCathodeClustersW(topNode);

    if(dMutCathodeClustersW.find("dMutCathodeClusters")==0)
    { cout<<"dMutCathodeClusters not found.\n";
      return -1;
    }

    (*dMutCathodeClustersW).SetRowCount(0);

    mergeMui(topNode, dst2Node);
    mergeMut(topNode, dst2Node);
    
    mMutFindCathodeClustsTM->event(topNode);
    mMutFitCathodeClustsTM->event(topNode);
  }
  
  TrackonDst = (*dMuoTracks).RowCount();
  (*dMuoTracks).SetRowCount(0); // also done in dMuoTrackRoadRel
   // What about trackroadrel ?
  (*dMuoTrackRoadRel).SetRowCount(0);

  (*dMuiClustersW).SetRowCount(0); // 
  (*dMuiRoads).SetRowCount(0); // 
  (*dMuiRoadRawRel).SetRowCount(0); // 
  
// Refind muon identifier clusters and roads:
  mMuiClusterFinder->event(topNode);
  mMuiRoadFind->event(topNode);
  mMuiFillRoadRawRel->event(topNode);

// Try to find tracks from roads and clusters in dsts
  mMuoFindTracksM->event(topNode);

  mMuiRoadFind2->event(topNode);

  TrackRefit = (*dMuoTracks).RowCount();

  if ( (TrackonDst+TrackRefit) > 0 )
    cout << "Tracks ondst / refit : " << TrackonDst << " / " << TrackRefit << endl;
}
//------------------------------------------------------------------ DISPLAY CODE
int setup_display(PHCompositeNode* topNode)
{
  mMuiInitModule mMuiInit;
  PHTimeStamp t_search(2001,8,8,0,0,0);
  mMuiInit.SetSearchTimeStamp(t_search);
  mMuiInit.event(topNode);

  PhenixRun* phenix = new PhenixRun("Phenix-Run","Phenix-Run", topNode);
  gPhenix->GetDisplay()->AddSubsystem(new PhMuiDisplay("Mui-Display","Mui-Display"));
  gPhenix->StartRun();
  PhEventDisplay* dmuon = new PhEventDisplay("Main",gClient->GetRoot(),1000,1000,gPhenix->GetDisplay()->GetList());
  PhCanvas* can = (PhCanvas*)(dmuon->GetCanvas());
  can->DefaultRange(-800,-800,-800,800,800,800);//in cm
  dmuon->Popup();
  can->ResetRange();

  return 0;
}

int ana_display(PHCompositeNode* topNode)
{
  gPhenix->event(topNode);
  gPhenix->Draw();
  return 0;
}

int mergeMut(PHCompositeNode* mainNode, PHCompositeNode* embedNode)
{
  bool oldst = false;
  PHIODataNode<dMutCalibCathodesOut>::iterator rCalibCathodes1(mainNode);
  if(rCalibCathodes1.find("dMutCalibCathodesOut")==0)
  { cout<<"dMutCalibCathodesOut not found.\n";
    oldst = true;
  }
  
  PHIODataNode<dMutCalibCathodesOut>::iterator rCalibCathodes2(embedNode);
  if(rCalibCathodes2.find("dMutCalibCathodesOut")==0)
  { cout<<"Embedding dMutCalibCathodesOut not found.\n";
    return -1;
  }

  int beforeCount = 0;
  int afterCount  = 0;
  if(!oldst)
  {
    dMutCalibCathodesOut* fCalibCathodes1 = &(*rCalibCathodes1);
    dMutCalibCathodesOut* fCalibCathodes2 = &(*rCalibCathodes2);
  
    beforeCount = fCalibCathodes1->RowCount();
    
    for(int calibid2 = 0; calibid2 < fCalibCathodes2->RowCount(); calibid2++)
    {
      bool calibfound = false;
      for(int calibid1 = 0; calibid1 < fCalibCathodes1->RowCount(); calibid1++)
      {
        if(calibfound) continue;
        
        //Are these the same cathodes
        if(
          fCalibCathodes1->get_arm(calibid1) == fCalibCathodes2->get_arm(calibid2) &&
          fCalibCathodes1->get_station(calibid1) == fCalibCathodes2->get_station(calibid2) &&
          fCalibCathodes1->get_octant(calibid1) == fCalibCathodes2->get_octant(calibid2) &&
          fCalibCathodes1->get_halfoctant(calibid1) == fCalibCathodes2->get_halfoctant(calibid2) &&
          fCalibCathodes1->get_gap(calibid1) == fCalibCathodes2->get_gap(calibid2) &&
          fCalibCathodes1->get_icath(calibid1) == fCalibCathodes2->get_icath(calibid2) &&
          fCalibCathodes1->get_plane(calibid1) == fCalibCathodes2->get_plane(calibid2) &&
          fCalibCathodes1->get_strip(calibid1) ==fCalibCathodes2->get_strip(calibid2)
          )
        {
          //Our signal file contains a hit that overlaps with the real data
          // So lets add the total charge of the hits and for now do nothing with adc values
          calibfound = true;
          fCalibCathodes1->set_qstrip(calibid1,
                                      fCalibCathodes1->get_qstrip(calibid1)
                                      + fCalibCathodes2->get_qstrip(calibid2)
                                      );
          
        }
        
        
      }
  
      //If the cathode was not hit in first file add it from second file
      if(!calibfound)
      {
        fCalibCathodes1->Add();
        int newID = fCalibCathodes1->RowCount() - 1;
        fCalibCathodes1->set_qstrip(newID, fCalibCathodes2->get_qstrip(calibid2));
        fCalibCathodes1->set_qrms(newID, fCalibCathodes2->get_qrms(calibid2));
        fCalibCathodes1->set_tzero(newID, fCalibCathodes2->get_tzero(calibid2));
        for(int adcid = 0 ; adcid < 4; adcid++)
        {
    fCalibCathodes1->set_charge(adcid, newID,
              fCalibCathodes2->get_charge(adcid, calibid2));
        }
        for(int mcid = 0 ; mcid < 10; mcid++)
        {
    fCalibCathodes1->set_MutMCIndex(mcid, newID,
            fCalibCathodes2->get_MutMCIndex(mcid, calibid2));
        }
        fCalibCathodes1->set_arm(newID, fCalibCathodes2->get_arm(calibid2));
        fCalibCathodes1->set_station(newID, fCalibCathodes2->get_station(calibid2));
        fCalibCathodes1->set_octant(newID, fCalibCathodes2->get_octant(calibid2));
        fCalibCathodes1->set_halfoctant(newID, fCalibCathodes2->get_halfoctant(calibid2));
        fCalibCathodes1->set_gap(newID, fCalibCathodes2->get_gap(calibid2));
        fCalibCathodes1->set_icath(newID, fCalibCathodes2->get_icath(calibid2));
        fCalibCathodes1->set_plane(newID, fCalibCathodes2->get_plane(calibid2));
        fCalibCathodes1->set_strip(newID, fCalibCathodes2->get_strip(calibid2));        
        
      }
    }

    afterCount = fCalibCathodes1->RowCount();

  }else
  {
    //Before pro.19 old wrapped staf table was used for Calibrated Cathodes
    PHIODataNode<dMutCalibCathodesWrapper>::iterator rCalibCathodes1(mainNode);
    if(rCalibCathodes1.find("dMutCalibCathodes")==0)
    { cout<<"dMutCalibCathodes not found.\n";
      return -1;
    }
    
    
    dMutCalibCathodesWrapper* fCalibCathodes1 = &(*rCalibCathodes1);
    dMutCalibCathodesOut* fCalibCathodes2 = &(*rCalibCathodes2);

    int beforeCount = fCalibCathodes1->RowCount();

    for(int calibid2 = 0; calibid2 < fCalibCathodes2->RowCount(); calibid2++)
    {
      bool calibfound = false;
      for(int calibid1 = 0; calibid1 < fCalibCathodes1->RowCount(); calibid1++)
      {
        if(calibfound) continue;

        //Are these the same cathodes
        if(
           fCalibCathodes1->get_arm(calibid1) == fCalibCathodes2->get_arm(calibid2) &&
           fCalibCathodes1->get_station(calibid1) == fCalibCathodes2->get_station(calibid2) &&
           fCalibCathodes1->get_octant(calibid1) == fCalibCathodes2->get_octant(calibid2) &&
           fCalibCathodes1->get_halfoctant(calibid1) == fCalibCathodes2->get_halfoctant(calibid2) &&
           fCalibCathodes1->get_gap(calibid1) == fCalibCathodes2->get_gap(calibid2) &&
           fCalibCathodes1->get_icath(calibid1) == fCalibCathodes2->get_icath(calibid2) &&
           fCalibCathodes1->get_plane(calibid1) == fCalibCathodes2->get_plane(calibid2) &&
           fCalibCathodes1->get_strip(calibid1) ==fCalibCathodes2->get_strip(calibid2)
           )
        {
          //Our signal file contains a hit that overlaps with the real data
          // So lets add the total charge of the hits and for now do nothing with adc values
          calibfound = true;
          fCalibCathodes1->set_qstrip(calibid1,
                                      fCalibCathodes1->get_qstrip(calibid1)
                                      + fCalibCathodes2->get_qstrip(calibid2)
                                      );

        }


      }

      //If the cathode was not hit in first file add it from second file
      if(!calibfound)
      {
        fCalibCathodes1->SetRowCount(fCalibCathodes1->RowCount() + 1);
        int newID = fCalibCathodes1->RowCount() - 1;
        fCalibCathodes1->set_qstrip(newID, fCalibCathodes2->get_qstrip(calibid2));
        fCalibCathodes1->set_qrms(newID, fCalibCathodes2->get_qrms(calibid2));
        fCalibCathodes1->set_tzero(newID, fCalibCathodes2->get_tzero(calibid2));
        for(int adcid = 0 ; adcid < 4; adcid++)
        {
          fCalibCathodes1->set_charge(adcid, newID,
                                      fCalibCathodes2->get_charge(adcid, calibid2));
        }
        for(int mcid = 0 ; mcid < 10; mcid++)
        {
          fCalibCathodes1->set_MutMCIndex(mcid, newID,
                                          fCalibCathodes2->get_MutMCIndex(mcid, calibid2));
        }
        fCalibCathodes1->set_arm(newID, fCalibCathodes2->get_arm(calibid2));
        fCalibCathodes1->set_station(newID, fCalibCathodes2->get_station(calibid2));
        fCalibCathodes1->set_octant(newID, fCalibCathodes2->get_octant(calibid2));
        fCalibCathodes1->set_halfoctant(newID, fCalibCathodes2->get_halfoctant(calibid2));
        fCalibCathodes1->set_gap(newID, fCalibCathodes2->get_gap(calibid2));
        fCalibCathodes1->set_icath(newID, fCalibCathodes2->get_icath(calibid2));
        fCalibCathodes1->set_plane(newID, fCalibCathodes2->get_plane(calibid2));
        fCalibCathodes1->set_strip(newID, fCalibCathodes2->get_strip(calibid2));

      }
    }
    afterCount = fCalibCathodes1->RowCount();
  }    
  cout<<"MutMerge: Before/After "<<beforeCount <<"/"<< afterCount<<endl;

  return 0;
}

int mergeMui(PHCompositeNode* mainNode, PHCompositeNode* embedNode)
{
  PHIODataNode<dMuiRawWrapper>::iterator rMuiRaw1(mainNode);
  if(rMuiRaw1.find("dMuiRaw")==0)
  { cout<<"Main dMuiRaw not found.\n";
    return -1;
  }

  PHIODataNode<dMuiRawWrapper>::iterator rMuiRaw2(embedNode);
  if(rMuiRaw2.find("dMuiRaw")==0)
  { cout<<"Embedding dMuiRaw not found.\n";
    return -1;
  }
  
  dMuiRawWrapper& fMuiRaw1 = (*rMuiRaw1);
  dMuiRawWrapper& fMuiRaw2 = (*rMuiRaw2);

  int beforeCount = fMuiRaw1.RowCount();
  
  for(int rawid2 = 0; rawid2 < fMuiRaw2.RowCount(); rawid2++)
  {
    bool rawfound = false;
    for(int rawid1 = 0; rawid1 < fMuiRaw1.RowCount(); rawid1++)
    {
      if(
        fMuiRaw1.get_arm(rawid1) == fMuiRaw2.get_arm(rawid2) &&
        fMuiRaw1.get_plane(rawid1) == fMuiRaw2.get_plane(rawid2) &&
        fMuiRaw1.get_panel(rawid1) == fMuiRaw2.get_panel(rawid2) &&
        fMuiRaw1.get_orientation(rawid1) == fMuiRaw2.get_orientation(rawid2) &&
        fMuiRaw1.get_twopack(rawid1) == fMuiRaw2.get_twopack(rawid2)
        )
      {
        rawfound = true;
      }
    }
    
    if(!rawfound)
    {
      int newid = fMuiRaw1.RowCount();
      fMuiRaw1.SetRowCount(newid + 1);
      fMuiRaw1.set_dMuiRawIndex(newid, newid + 10000);
      fMuiRaw1.set_arm(newid, fMuiRaw2.get_arm(rawid2));
      fMuiRaw1.set_plane(newid, fMuiRaw2.get_plane(rawid2));
      fMuiRaw1.set_panel(newid, fMuiRaw2.get_panel(rawid2));
      fMuiRaw1.set_orientation(newid, fMuiRaw2.get_orientation(rawid2));
      fMuiRaw1.set_twopack(newid, fMuiRaw2.get_twopack(rawid2));
    }
  }

    cout<<"MuiMerge: Before/After "<<beforeCount <<"/"<< fMuiRaw1.RowCount()<<endl;

  
  return 0;
}