void Fun4All_SIM(
		 int nEvents = 0,
		 int run = 427673, //test
		 std::string pisaFile = "PISAEvent.root",
		 std::string outputName  = "singlemuon_pdst.root"
		)
{

  /////////////////////////
  bool use_deadmaps = false;
  bool use_sim_vtx = true;
  bool add_noise = false;
  bool use_vtx_deadmaps = false;
  bool use_fvtx = true;
  bool use_svx = true;

  bool fill_mc_info = true;
  bool fill_MutrRefit_info = true;
 
  bool write_dst = false;
  bool write_pdst = false;

  bool use_eval = true;

	bool set_beam_seed = false;
	float beam_x = 0.1678; //Run15pp200 beam avg
	float beam_y = 0.0640; //Run15pp200 beam avg
  
  bool debug = false;
  int verbosity = 0;
  if(debug) verbosity = 5;
  /////////////////////////

  if(outputName.empty())
	{
		cout << "outputName is empty! Exiting..." << endl;
		return;
	}

  //std::string outputDST = outputName+Form("_%d_dst.root",run);
  //std::string outputNDST = outputName+Form("_%d_ndst.root",run);
  //std::string outputPDST = outputName+Form("_run%d_seg%05d.root",run,seg);
  std::string outputPDST = outputName;

  // load libraries
  gSystem->Load("libfun4all.so");
  gSystem->Load("libPISARoot.so");
  gSystem->Load("libfun4allfuncs.so");
  gSystem->Load("libmutoo_subsysreco.so");
  gSystem->Load("libfvtx_subsysreco.so");
  gSystem->Load("libMWG_interface.so");
  gSystem->Load("libPHGeant" );
  gSystem->Load("libsimreco_base");
  gSystem->Load("libmuon_subsysreco");
  gSystem->Load("librpc_subsysreco");
  gSystem->Load("librpc_muotrackreco");
  gSystem->Load("libfun4allfuncs_muons");
  gSystem->Load("libsvx");
  gSystem->Load("liblvl2");
  gSystem->Load("libPythia6.so");
  gSystem->Load("libPHPythia.so");
  gSystem->Load("libPHPythiaEventGen.so");
  gSystem->Load("libpicodst_object.so");
  
  cout << ">>>> libraries loaded" << endl;
  
  ///////////////////////////////////////////
  // recoConsts setup
  //////////////////////////////////////////
  recoConsts *rc = recoConsts::instance();
  rc->set_IntFlag("SVXACTIVE", use_svx);
  rc->set_IntFlag("PRINT_MUTOO_PARAMETERS",1);
  rc->set_IntFlag("RUNNUMBER", run);
  rc->set_IntFlag("MUONFUN4SIM",1);
  rc->set_IntFlag("SIMULATIONFLAG",2);
  rc->set_IntFlag("EMBEDFLAG",0);
  rc->set_CharFlag("EMBED_MC_TOPNODE","TOP");  // need to use "TOP" to get MC inoformation                                                                                                                                                                                                   

  ///////////////////////////////////////////
  // FVTX+MUTR setup
  //////////////////////////////////////////
  TFvtxGlobalParCntrl::set_bool_par("is_pp",true);
  TFvtxGlobalParCntrl::set_bool_par("use_svx",use_svx);
  TFvtxGlobalParCntrl::set_bool_par("is_sim",true);
  TFvtxGlobalParCntrl::set_bool_par("deadmap_use_calibration_database", use_deadmaps);
  TFvtxGlobalParCntrl::set_bool_par("deadmap_use_production_map_for_sim", use_deadmaps);
  TFvtxGlobalParCntrl::set_bool_par("geom_use_calibration_database", false); 
  TFvtxGlobalParCntrl::set_string_par("geom_root_file_path","./");
  TFvtxGlobalParCntrl::set_string_par("geom_root_file_name", "fvtx_geometry.root");

	if ( set_beam_seed ) TFvtxGlobalParCntrl::set_float_par("beam_x_seed", beam_x);
	if ( set_beam_seed ) TFvtxGlobalParCntrl::set_float_par("beam_y_seed", beam_y);
	TFvtxGlobalParCntrl::set_bool_par("beam_use_average_xy", set_beam_seed);

	TMutDatabaseCntrl::set_database_access("use_local_dead_HV_file",true); //mut.disabledAnodes.dat
	TString hv_file = Form("/direct/phenix+hhj/shlim/work/07.genericQA/01.mutr_database/Run15pp200/mut.disabledAnodes_Run15pp200/mut.disabledAnodes.dat_run%d",run);
  TMutDatabaseCntrl::set_filename("use_local_dead_HV_file", hv_file.Data());
	TMutDatabaseCntrl::set_database_access("use_local_dead_wire_file",false); //mut.disabledWires.dat
	TMutDatabaseCntrl::set_database_access("use_local_attenuated_chnl_file",false); //AttenuatedChannels.txt
	TMutDatabaseCntrl::set_database_access("use_local_dead_channel_file",false); //DeadChannels.dat

  ///////////////////////////////////////////
  // MUID setup
  //////////////////////////////////////////
	std::string muid_eff_south("muid_tube_eff_south_Run15pp200.txt");
	std::string muid_eff_north("muid_tube_eff_north_Run15pp200.txt");

	{
		std::cout << "using local two pack efficiciency files:" << muid_eff_south << ", "  << muid_eff_north << std::endl;
		TMuiHVMask::set_mode(TMuiHVMask::FROM_FILE);
		TMuiHVMask::set_filename_south(muid_eff_south.c_str());
		TMuiHVMask::set_filename_north(muid_eff_north.c_str());
	}


  TMutExtVtx::get().set_verbosity( MUTOO::SOME );
  if(!use_sim_vtx)
	{
		TMutExtVtx::get().set_priorities( "FVTX", 0 );
		TMutExtVtx::get().set_priorities( "SVX_PRECISE", 1 );
		TMutExtVtx::get().set_priorities( "BBC", 2 );
	}
	else{
		TMutExtVtx::get().set_priorities( "SIM", 0 );
	}
  //mMfmMT::setMapFileFlag( mMfmMT::MAP_3D_PLUS_PLUS );
  //mMfmMT::setMapFileScale( 1.0 );
  //MuonUtil::set_check_mapfile_scale( false );

  ///////////////////////////////////////////
  // Make the Server
  //////////////////////////////////////////
  Fun4AllServer *se = Fun4AllServer::instance();
  se->Verbosity(verbosity);


  ///////////////////////////////////////////
  // Subsystems
  //////////////////////////////////////////

  // run header and trigger setting
  se->registerSubsystem( new HeadSimreco() );
  //se->registerSubsystem( new TrigSimreco() );

  // vertex simulation
  // puts the vertex from the pisa header node into vtxOut object
  VtxSimreco* vtx_reco = new VtxSimreco();
  vtx_reco->SmearZ( false );
  vtx_reco->SmearY( false );
  vtx_reco->SmearX( false );
  vtx_reco->UseXY( true );
  vtx_reco->OverwriteBBC( true );
  vtx_reco->XVertexSigma( 0.01 ); // VTX resolution of 20 microns - current value is arbitrary (error doesn't affect DCA calculation)
  vtx_reco->YVertexSigma( 0.01 );
  vtx_reco->ZVertexSigma( 0.02 );
  se->registerSubsystem( vtx_reco );

  // Counter
  se->registerSubsystem( new MuonCounter() );

  // global detectors subsystem
  //se->registerSubsystem( new HeadReco() );
  //se->registerSubsystem( new TrigReco( ));
  //se->registerSubsystem( new BbcReco() );
  //se->registerSubsystem( new ZdcReco() );
  se->registerSubsystem( new VtxReco() );


  if(use_svx)
	{
		// SVX reconstruction
		SvxParManager *svxpar = new SvxParManager();
		svxpar->Verbosity(verbosity);
		svxpar->set_UseRefDiffPixelMap(1);

		//svxpar->set_BeamCenter(0.1247,-0.1601); //Run12pp510
		//svxpar->set_BeamCenter(beam_x,beam_y); //Run15pp200
		svxpar->set_OffsetVtxToCnt(0.0, 0.0, 0.0);
		svxpar->set_OffsetEastToWest(0.0, 0.0, 0.0);
		svxpar->set_ReadGeoParFromFile(1);
		svxpar->set_GeometryFileName("svxPISA.par"); //Run15pp200 - temporary
		svxpar->set_UseStripThresholdDatbase(true);
		se->registerSubsystem(svxpar);

		SvxSimulator *svxsim = new SvxSimulator();
		svxsim->Verbosity(10);
		se->registerSubsystem(svxsim);

		SvxApplyHotDead *svxhotdead = new SvxApplyHotDead();
		svxhotdead->Verbosity(verbosity);
		if(use_vtx_deadmaps) se->registerSubsystem(svxhotdead);

		SvxReco *svxrec = new SvxReco();
		svxrec->Verbosity(verbosity);
		svxrec->set_ThisIsSimulation();
		svxrec->set_StripixelAdcSumThreshold(0);
		se->registerSubsystem(svxrec);

		SvxPriVertexSeedFinder *svxvtxseedfinder = new SvxPriVertexSeedFinder();
		svxvtxseedfinder->Verbosity(verbosity);
		se->registerSubsystem(svxvtxseedfinder);

		SvxStandAloneReco *svxstandalone = new SvxStandAloneReco();
		svxstandalone->Verbosity(verbosity);
		svxstandalone->setVertexRecoFlag(2);
		svxstandalone->setPPFlag(true);
		se->registerSubsystem(svxstandalone);

		SvxPrimVertexFinder *svxprimvtxfinder = new SvxPrimVertexFinder();
		svxprimvtxfinder->Verbosity(verbosity);
		se->registerSubsystem(svxprimvtxfinder);
	}

  // muon prdf unpacker
  MuonUnpackPisa *muon_unpack_pisa =  new MuonUnpackPisa();
  muon_unpack_pisa->Verbosity( 1 );
  muon_unpack_pisa->set_flag(MuonUnpackPisa::DO_RESPONSE,1);
  //if(add_noise)muon_unpack_pisa->set_flag(MuonUnpackPisa::ADD_NOISE,1);
  //muon_unpack_pisa->set_flag(MuonUnpackPRDF::SKIP_ZERO_SUPPRESSION,1);
  se->registerSubsystem( muon_unpack_pisa );

  // fvtx prdf unpacker
  if (use_fvtx)
	{
		FvtxUnpackPisa *fvtx_unpack = new FvtxUnpackPisa();
		//fvtx_unpack->set_do_response( true );
		//fvtx_unpack->Verbosity(verbosity);
		se->registerSubsystem( fvtx_unpack );

		FvtxReco* fvtxreco = new FvtxReco();
		//fvtxreco->set_do_mutkalfiteval(true);         // Adds VTX hits to tracking
		//fvtxreco->run_after_burner(false);         // Adds VTX hits to tracking
		//fvtxreco->set_use_svx_cluster(true);     // Uses smeared PISA hits if false
		//fvtxreco->set_do_mutr_matching(false);     // Match FVTX-MuTr tracks and refit or not
		//fvtxreco->set_finder(3);                  // Use Hough track finder
		se->registerSubsystem(fvtxreco);

		TMutNode<mFvtxFindTrackPar>::find_node(se->topNode(),"mFvtxFindTrackPar")->print();

		FvtxPrimVertex* fvtxprimvtx = new FvtxPrimVertex();
		fvtxprimvtx->set_source(FvtxPrimVertex::Tracks,FvtxPrimVertex::Segments);
		if(!use_sim_vtx)se->registerSubsystem(fvtxprimvtx);
	}

  // mutoo reconstruction
  se->registerSubsystem( new MuiooReco() );

  MuonDev *muon_dev = new MuonDev();
  muon_dev->set_flag(MuonDev::SKIP_VERTEX_FIT, 1);
  muon_dev->set_flag(MuonDev::USE_MUID, false);
  //Option to force the cluster fitting to use a particular method:
  //muon_dev->set_flags(MuonDev::FORCE_LOOKUP);
  se->registerSubsystem( muon_dev );

  if(use_fvtx)
	{
		// Perform FVTX-Mutr track matching and refit track:
		FvtxRecoWithMut *fvtx_reco_withmut = new FvtxRecoWithMut();
		se->registerSubsystem( fvtx_reco_withmut );
	}

  if(use_eval)
	{
		MuonEval* mueval = new MuonEval();
		mueval->set_flags(0); // no ntuple output needed
		se->registerSubsystem (mueval);
	}

  if (use_fvtx && use_eval)
	{
		FvtxEval* fvtxeval = new FvtxEval("FvtxEval","fvtx_eval_pisa.root");
		se->registerSubsystem(fvtxeval);
	}

	se->registerSubsystem( new GlobalReco() );
	se->registerSubsystem( new GlobalReco_muons() );
  
  
  ///////////////////////////////////////////
  // IOManagers...
  ///////////////////////////////////////////

  // dst
  if( write_dst ) 
	{
		Fun4AllDstOutputManager *dstManager  = new Fun4AllDstOutputManager("DSTOUT",outputDST.c_str());

		dstManager->AddNode("Sync");
		dstManager->AddNode("RunHeader");
		dstManager->AddNode("EventHeader");
		dstManager->AddNode("PreviousEvent");
		dstManager->AddNode("TrigRunLvl1");
		dstManager->AddNode("VtxOut");
		dstManager->AddNode("BbcOut");
		//dstManager->AddNode("BbcRaw");
		dstManager->AddNode("ZdcOut");
		//dstManager->AddNode("ZdcRaw");

		dstManager->AddNode("TMCPrimary");
		dstManager->AddNode("PHPythiaHeader");
		dstManager->AddNode("PHPythia");

		dstManager->AddNode("TrigLvl1");
		//dstManager->AddNode("L2Decision");
		//dstManager->AddNode("Lvl2OutArray");

		// Muioo nodes
		dstManager->AddNode("TMuiHitO");
		dstManager->AddNode("TMuiMCHitO");
		dstManager->AddNode("TMuiMCHitMapO");


		// Mutoo nodes
		dstManager->AddNode("TMutHit");
		dstManager->AddNode("TMutMCHit");
		dstManager->AddNode("TMutMCTrk");
		dstManager->AddNode("TMutMCTrkMap");
		dstManager->AddNode("TMutHitMap");

		dstManager->AddNode("TRpcMCHitMap");
		dstManager->AddNode("TRpcHitMap");

		dstManager->AddNode("PHMuoTracksOO");

		// FVTX nodes
		dstManager->AddNode("TFvtxHit");
		dstManager->AddNode("TFvtxHitMap");
		dstManager->AddNode("TFvtxMCHit");
		dstManager->AddNode("TFvtxPisaHit");
		dstManager->AddNode("TFvtxCompactTrk");

		// SVX nodes
		dstManager->AddNode("McSingle");
		dstManager->AddNode("SvxPisaHit");
		dstManager->AddNode("SvxGhitList");
		dstManager->AddNode("SvxRawhitList");
		dstManager->AddNode("SvxGhitRawhitList");
		//dstManager->AddNode("PHCentralTrack");
		//dstManager->AddNode("SvxCentralTrackList");

		// From EVA node
		dstManager->AddNode("header");
		dstManager->AddNode("fkin");
		dstManager->AddNode("primary");
		//dstManager->AddNode("pythia");

		// PHGlobal
		dstManager->AddNode("PHGlobal");
		dstManager->AddNode("PHGlobal_MUON");

		dstManager->AddNode("PHPythiaHeader");
		dstManager->AddNode("PHPythia");

		se->registerOutputManager(dstManager);
	}


  // picoDST
  if( write_pdst )
	{
		// MWG
		gSystem->Load("libMWGOO.so");
		PHInclusiveNanoCuts *MWGcuts = new MWGInclusiveNanoCutsv2();
		se->registerSubsystem(new MWGFvtxReco(MWGcuts));

		mFillSingleMuonContainer* msngl = new mFillSingleMuonContainer();
		msngl->set_lastgap_cut(2);
		//msngl->set_bbcz_cut(50.0);
		//msngl->set_pz_cut(3.0);
		//msngl->set_DG0_cut(60);
		//msngl->set_DDG0_cut(40);
		//msngl->set_chi2_cut(15);
		msngl->set_is_sim(use_sim_vtx);
		se->registerSubsystem(msngl);

		if(fill_mc_info)
		{

			//********************* Mutr MC information Module ******************//
			//mFillMCSingleMuonContainer* msngl_mc = new mFillMCSingleMuonContainer();
			//se->registerSubsystem(msngl_mc);

			mFillMCSingleMuonFvtxContainer* msngl_fvtx_mc = new mFillMCSingleMuonFvtxContainer();
			se->registerSubsystem(msngl_fvtx_mc);

			//mFillMCDiMuonContainer* mdi_mc = new mFillMCDiMuonContainer();
			//se->registerSubsystem(mdi_mc);

		}

		Fun4AllOutputManager *outsmu = new Fun4AllDstOutputManager("Outsmu",outputPDST.c_str());
		outsmu->AddNode("SingleMuonContainer");
		if (fill_MutrRefit_info) outsmu->AddNode("MutrRefitSingleMuonContainer");
		//if (fill_mc_info) outsmu->AddNode("MCSingleMuonContainer");
		if(fill_mc_info && use_fvtx) outsmu->AddNode("MCSingleMuonFvtxContainer");
		outsmu->AddNode("Sync");
		outsmu->AddNode("TrigLvl1");
		outsmu->AddNode("VtxOut");
		//outsmu->AddNode("PHGlobal");
		//outsmu->AddNode("RunHeader");
		//outsmu->AddNode("TFvtxCompactTrk");
		//outsmu->AddNode("TFvtxCompactTrkMap");
		//outsmu->AddNode("FvtxPrimVtxContainer");
		//outsmu->AddNode("PHPythiaHeader");
		//outsmu->AddNode("PHPythia");
		//outsmu->AddNode("DiMuonContainer");
		//outsmu->AddNode("TriggerEmulatorContainer");
		//if(fill_mc_info) outsmu->AddNode("MCHepMCParticleContainer");
		//if(fill_rpc_info) outsmu->AddNode("RpcSingleMuonContainer");
		outsmu->AddEventSelector("mFillSingleMuonContainer");
		se->registerOutputManager(outsmu);

	}


  //VertexPrinter *printer = new VertexPrinter();
  //se->registerSubsystem(printer);
  
  Fun4AllPisaInputManager *inMan = new Fun4AllPisaInputManager("PisaIn");
  se->registerInputManager(inMan);
  TString tmp = pisaFile;
  std::string suffix = ".txt";
	if(tmp.EndsWith(suffix.c_str()))
	{
		int counter = 0;
		std::string theFile;
		ifstream inFile;
		inFile.open(pisaFile.c_str());
		while(inFile >> theFile)
		{
			inMan->AddFile(theFile.c_str());
		}
		inFile.close();
	}
	else
	{
		se->fileopen(inMan->Name(),pisaFile.c_str());
	}


	gSystem->ListLibraries();
  se->run(nEvents);

	/*
	Fun4AllHistoManager *hm= se->getHistoManager("mFvtxPostProductionAlignment");
	if ( hm ) hm->dumpHistos("histogram_local.root");
	*/

  se->End();


  //Fun4AllHistoManager *hm = se->getHistoManager("QA");
  //if( hm ) hm->setOutfileName( "QA.root" );
  //else cout << "Fun4Muons_Pisa - unable to find QA histograms" << endl;

  cout << ">>>> Completed reconstruction." << endl;
}