// $Id: Fun4Muons_Pisa.C,v 1.16 2004/04/23 12:49:14 hpereira Exp $
/*! 
	runs mutoo slow simulator on a PisaFile, create a Simulated DST
*/
void Fun4Muons_Pisa(
	int nEvents=0, 
	const char *filein = "pisa.root",
	const char *dstout = "dst_out.root",
	char *prdfout = "data_out.prdf",
	int run_number = 127535	)
{

	cout << "Fun4Muons_Pisa - filein: " << filein << endl;
	cout << "Fun4Muons_Pisa - dstout: " << dstout << endl;
	cout << "Fun4Muons_Pisa - prdfout: " << prdfout << endl;
	cout << "Fun4Muons_Pisa - run_number: " << run_number << endl;
	
	bool write_dst = true;			 
	bool write_prdf = false;
	bool add_disalignment = false;
	bool use_rpc = false;					
	bool use_fvtx = true;					
	bool do_eval = true;
	
	// load libraries
	gSystem->Load("libfun4all.so");
	gSystem->Load("libfun4allfuncs.so");
	gSystem->Load("libsimreco.so");
	gSystem->Load("libmutoo_subsysreco.so");
	gSystem->Load("libfvtx_subsysreco.so");

	// Read in MuTr alignment file to get proper "misalignment" of muon tracker stations
	// to match Run 4 real data:
	if( add_disalignment )
	TMutDatabaseCntrl::set_database_access("use_local_internal_align_file",true);

	///////////////////////////////////////////
	// recoConsts setup
	//////////////////////////////////////////
	recoConsts *rc = recoConsts::instance();

	// simVertexFlag = 0 (default) means that the BBC Z0 value will be used
	// simVertexFlag = 1 means that the same simZ0Vertex centroid value is used for all events
	// simVertexFlag = 2 means that the Z0 centroid is taken from the PISA event header for each event
	// The centroid values are modified by the Width values which are Gaussian sigma values
	rc->set_IntFlag("SIMVERTEXFLAG",2);	
	rc->set_FloatFlag("SIMZ0VERTEX",0);
	rc->set_FloatFlag("SIMT0VERTEX",0);	
	rc->set_FloatFlag("SIMZ0VERTEXWIDTH",2.0);
	rc->set_FloatFlag("SIMT0VERTEXWIDTH",0.0);
	
	rc->set_IntFlag("SIMULATIONFLAG",2);
	rc->set_IntFlag("EMBEDFLAG",0);
	rc->set_IntFlag("RUNNUMBER", run_number );
	
	// muid tube simulation efficiency
	rc->set_DoubleFlag("MUIOO_TUBE_EFF",1.0); 

	/* 
		disable check of mapFileScale since its decided by pisa, 
		not by the run number
	*/
	mMfmMT::setMapFileFlag( 4 );
	mMfmMT::setMapFileScale( 1.0 );
	MuonUtil::set_check_mapfile_scale( false );
	
	///////////////////////////////////////////
	// Make the Server
	///////////////////////////////////////////
	Fun4AllServer *se = Fun4AllServer::instance();

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

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

	// BBC/ZDC reconstruction
	// se->registerSubsystem( new BbcSimreco("BBC") );
	// se->registerSubsystem( new VtxReco("VTX") );
	
	// muon reconstruction
	MuonUnpackPisa *Muon_pisa = new MuonUnpackPisa();
	se->registerSubsystem(Muon_pisa);  //xxx was disabled
	Muon_pisa->Verbosity( 1 );

	// Rpc reconstruction
	if( use_rpc ) se->registerSubsystem(new RpcUnpackPisa()); 
	
	// Fvtx 
	if( use_fvtx ) {
		FvtxUnpackPisa *fvtx_unpack_pisa = new FvtxUnpackPisa();
		fvtx_unpack_pisa->Verbosity(1);
		fvtx_unpack_pisa->set_do_response( false );
		se->registerSubsystem( fvtx_unpack_pisa );
		if( do_eval ) se->registerSubsystem( new FvtxMCEval() );
	}
	
	if( prdfout && write_prdf ) 
	{
		
		// tell module to run the response and write it to the DST
		Muon_pisa->set_do_response( true );

		// tell module to write PRDF output if required
		se->registerSubsystem( new MuonPackPRDF() );
	
	}
	
	/* 
		Read in alignment correction file to get MuID misaligned to match Run 4
		real data. Update muid geometry consequently.
	*/
	if( add_disalignment )
	se->registerSubsystem( new MuonDisalign("MUONDISALIGN", "alignment_corrections.txt") );
	
	// map dumper
	MuonDumpMap *dump_map = new MuonDumpMap();
	dump_map->add_map( "TFvtxMCHitMap" );
	// se->registerSubsystem( dump_map );
	
	///////////////////////////////////////////
	// Input manager
	///////////////////////////////////////////
	Fun4AllInputManager *input_manager = new Fun4AllPisaInputManager("PisaIn");
	se->registerInputManager(input_manager);
	se->fileopen(input_manager->Name(),filein);
	
	///////////////////////////////////////////
	// Output manager
	///////////////////////////////////////////
	
	// DST output
	if( dstout && write_dst ) {
		
		Fun4AllDstOutputManager *dst_manager	= new Fun4AllDstOutputManager("DSTOUT",dstout);
		dst_manager->AddEventSelector(Muon_pisa->Name());	
		se->registerOutputManager(dst_manager);
		
		// Header and vertex nodes
		dst_manager->AddNode("RunHeader");	
		dst_manager->AddNode("EventHeader");
		dst_manager->AddNode("VtxOut");
		dst_manager->AddNode("BbcOut");
		dst_manager->AddNode("BbcRaw");
		dst_manager->AddNode("ZdcOut");
		dst_manager->AddNode("ZdcRaw");
		
		// RPC geometry
		if( use_rpc )
		{
			dst_manager->AddNode("RpcStationType1");	
			dst_manager->AddNode("RpcStationType2");	
		}
		
		// mutoo/muioo nodes
		dst_manager->AddNode("TMCPrimary");
		dst_manager->AddNode("TMuiMCHitO");
		dst_manager->AddNode("TMutMCHit");
		dst_manager->AddNode("TMutMCTrk");
		
		// RPC geometry
		if( use_fvtx )
		dst_manager->AddNode("TFvtxMCHit");	
		
		// RPC node
		if( use_rpc )
		dst_manager->AddNode("TRpcMCHit");
		 
		// From EVA node
		dst_manager->AddNode("header");
		dst_manager->AddNode("fkin");
		dst_manager->AddNode("primary");
		dst_manager->AddNode("pythia");
	}
	
	/* 
		prdf output
		up to now only the muon tracker/muid will be written
		to the prdf, since the other detectors (simreco modules)
		cannot write both on DST and PRDF at the same time 
		(this is controlled by the SIMULATION flag)
	*/
	if( prdfout && write_prdf ) {
		Fun4AllPrdfOutputManager *prdf_io	= new Fun4AllPrdfOutputManager("PRDFOUT",prdfout);
		prdf_io->AddEventSelector(Muon_pisa->Name());	
		prdf_io->InitPrdfNode( se->topNode() );
		se->registerOutputManager(prdf_io);
	}
	 
	///////////////////////////////////////////
	// Analyze the Data.
	//////////////////////////////////////////
	TFile f(filein);
	int jEvents = T->GetEntries();
	if(nEvents==0 || nEvents > jEvents ) {
		cout << "Fun4Muons_PISA - changing nEvents to " << jEvents << endl;
		nEvents = jEvents;
	}
	se->run(nEvents);
	se->End();
	
	cout << "Completed reconstruction." << endl;
	
}