/**
 * \file global_sim.C
 *
 * \brief Simulation macro for tutorial on the CBM Software Workshop.
 *
 * \author Andrey Lebedev <andrey.lebedev@gsi.de>
 * \date 2011
 **/

void global_reco(Int_t nEvents = 100)
{
	TString script = TString(gSystem->Getenv("SCRIPT"));
	TString parDir = TString(gSystem->Getenv("VMCWORKDIR")) + TString("/parameters");

	TList *parFileList = new TList();

	// Output directory
	TString dir = "/data.local1/andrey/events/std_muon/";
	// MC transport file
	TString mcFile = dir + "mc.0000.root";
	// Parameters file
	TString parFile = dir + "param.0000.root";
	// Output file with reconstructed tracks and hits
	TString globalRecoFile = dir + "global.reco.0000.root";
	// Digi scheme for STS
	//TObjString stsDigiFile = parDir + "/sts/sts_v11a.digi.par";
	//parFileList->Add(&stsDigiFile);
	// Directory for output results
	TString resultsDir = "./results/";

	if (script == "yes") {
		dir = TString(gSystem->Getenv("OUTPUT_DIR"));
		mcFile = dir + "mc.0000.root";
		parFile = dir + "param.0000.root";
		globalRecoFile = dir + "global.reco.0000.root";
		resultsDir = TString(gSystem->Getenv("RESULTS_DIR"));
		TObjString stsDigiFile = parDir + "/" + TString(gSystem->Getenv("STS_DIGI"));
		parFileList->Add(&stsDigiFile);
	}

	Int_t iVerbose = 1;
	TStopwatch timer;
	timer.Start();

	gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
	basiclibs();
	gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/cbmrootlibs.C");
	cbmrootlibs();
	gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/determine_setup.C");

	FairRunAna* run = new FairRunAna();
	run->SetInputFile(mcFile);
	run->SetOutputFile(globalRecoFile);


	if (IsMvd(parFile)) {
		// ----- MVD reconstruction    --------------------------------------------
		CbmMvdDigitizeL* mvdDigi = new CbmMvdDigitizeL("MVD Digitiser", 0, iVerbose);
		run->AddTask(mvdDigi);

		CbmMvdFindHits* mvdHitFinder = new CbmMvdFindHits("MVD Hit Finder", 0, iVerbose);
		run->AddTask(mvdHitFinder);
		// -------------------------------------------------------------------------
	}

	// ----- STS REAL reconstruction -----------------------------------------------
	Double_t threshold  =  4;
	Double_t noiseWidth =  0.01;
	Int_t    nofBits    = 20;
	Double_t minStep    =  0.01;
	Double_t StripDeadTime = 0.1;
	CbmStsDigitize* stsDigitize = new CbmStsDigitize("STS Digitiser", iVerbose);
	stsDigitize->SetRealisticResponse();
	stsDigitize->SetFrontThreshold (threshold);
	stsDigitize->SetBackThreshold  (threshold);
	stsDigitize->SetFrontNoiseWidth(noiseWidth);
	stsDigitize->SetBackNoiseWidth (noiseWidth);
	stsDigitize->SetFrontNofBits   (nofBits);
	stsDigitize->SetBackNofBits    (nofBits);
	stsDigitize->SetFrontMinStep   (minStep);
	stsDigitize->SetBackMinStep    (minStep);
	stsDigitize->SetStripDeadTime  (StripDeadTime);
	run->AddTask(stsDigitize);

	FairTask* stsClusterFinder = new CbmStsClusterFinder("STS Cluster Finder",iVerbose);
	run->AddTask(stsClusterFinder);

	FairTask* stsFindHits = new CbmStsFindHits("STS Hit Finder", iVerbose);
	run->AddTask(stsFindHits);

	FairTask* stsMatchHits = new CbmStsMatchHits("STS Hit Matcher", iVerbose);
	run->AddTask(stsMatchHits);

	FairTask* kalman = new CbmKF();
	run->AddTask(kalman);
	FairTask* l1 = new CbmL1();
	run->AddTask(l1);
	CbmStsTrackFinder* trackFinder = new CbmL1StsTrackFinder();
	FairTask* findTracks = new CbmStsFindTracks(iVerbose, trackFinder);
	run->AddTask(findTracks);

	FairTask* stsMatchTracks = new CbmStsMatchTracks("STSMatchTracks", iVerbose);
	run->AddTask(stsMatchTracks);

	CbmStsTrackFitter* trackFitter = new CbmStsKFTrackFitter();
	FairTask* fitTracks = new CbmStsFitTracks("STS Track Fitter", trackFitter, iVerbose);
	run->AddTask(fitTracks);

	//	FairTask* stsFHQa = new CbmStsFindHitsQa("STSFindHitsQA",iVerbose);
	//	run->AddTask(stsFHQa);

	//  FairTask* stsRecoQa = new CbmStsReconstructionQa(kFALSE, 4, 0.7, 1);
	//	run->AddTask(stsRecoQa);
	// ------------------------------------------------------------------------

	if (IsTrd(parFile)) {
		// ----- TRD reconstruction-----------------------------------------
		// Update of the values for the radiator F.U. 17.08.07
		Int_t trdNFoils = 130; // number of polyetylene foils
		Float_t trdDFoils = 0.0013; // thickness of 1 foil [cm]
		Float_t trdDGap = 0.02; // thickness of gap between foils [cm]
		Bool_t simpleTR = kTRUE; // use fast and simple version for TR production

		CbmTrdRadiator *radiator = new CbmTrdRadiator(simpleTR, trdNFoils, trdDFoils, trdDGap);

		// ----- TRD hit smearing -----
		Double_t trdSigmaX[] = { 300, 400, 500 }; // Resolution in x [mum]
		// Resolutions in y - station and angle dependent [mum]
		Double_t trdSigmaY1[] = { 2700, 3700, 15000, 27600, 33000, 33000, 33000 };
		Double_t trdSigmaY2[] = { 6300, 8300, 33000, 33000, 33000, 33000, 33000 };
		Double_t trdSigmaY3[] = { 10300, 15000, 33000, 33000, 33000, 33000, 33000 };

		CbmTrdHitProducerSmearing* trdHitProd =
			  new CbmTrdHitProducerSmearing("TRD Hitproducer", "TRD task", radiator);
		trdHitProd->SetSigmaX(trdSigmaX);
		trdHitProd->SetSigmaY(trdSigmaY1, trdSigmaY2, trdSigmaY3);
		run->AddTask(trdHitProd);
		// ----- End TRD hit smearing -----
	}

	if (IsTof(parFile)) {
		// ------ TOF hits --------------------------------------------------------
		CbmTofHitProducer* tofHitProd = new CbmTofHitProducer(
				"TOF HitProducer", 1);
		run->AddTask(tofHitProd);
		// ------------------------------------------------------------------------
	}

	CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks();
	// Tracking method to be used
	// "branch" - branching tracking
	// "nn" - nearest neighbor tracking
	// "weight" - weighting tracking
	finder->SetTrackingType("branch");
	// Hit-to-track merger method to be used
	// "nearest_hit" - assigns nearest hit to the track
	finder->SetMergerType("nearest_hit");
	run->AddTask(finder);

	if (IsTrd(parFile)) {
		CbmTrdMatchTracks* trdMatchTracks = new CbmTrdMatchTracks(1);
		run->AddTask(trdMatchTracks);
	}

	if (IsRich(parFile)) {
		// ---------------------RICH Hit Producer ----------------------------------
		Double_t richPmtRad = 0.4; // PMT radius [cm]
		Double_t richPmtDist = 0.; // Distance between PMTs [cm]
		Int_t richDetType = 4; // Detector type Hamamatsu H8500-03 (no WLS)
		Int_t richNoise = 220; // Number of noise points per event
		Double_t richCollEff = 1.0; // Collection Efficiency of PMT electron optics
		Double_t richSMirror = 0.06; // Sigma for additional point smearing due to light scattering in mirror

		CbmRichHitProducer* richHitProd = new CbmRichHitProducer(richPmtRad,
			richPmtDist, richDetType, richNoise, iVerbose, richCollEff, richSMirror);
		run->AddTask(richHitProd);
		//--------------------------------------------------------------------------

		//----------------------RICH Track Extrapolation ---------------------------
		Int_t richNSts = 4; // minimum number of STS hits for extrapolation
		Double_t richZPos = 300.; // z position for extrapolation [cm]
		CbmRichTrackExtrapolation* richExtra = new CbmRichTrackExtrapolationKF(
			richNSts, iVerbose);
		CbmRichExtrapolateTracks* richExtrapolate = new CbmRichExtrapolateTracks();
		richExtrapolate->UseExtrapolation(richExtra, richZPos);
		run->AddTask(richExtrapolate);
		//--------------------------------------------------------------------------

		//--------------------- Rich Track Projection to photodetector -------------
		Int_t richZFlag = 1; // Projetion from IM plane (default)
		CbmRichProjectionProducer* richProj = new CbmRichProjectionProducer(iVerbose, richZFlag);
		run->AddTask(richProj);
		//--------------------------------------------------------------------------

		//--------------------- RICH Ring Finding ----------------------------------
		TString richGeoType = "compact";//choose between compact or large
		CbmRichRingFinderHough* richFinder = new CbmRichRingFinderHough(iVerbose,	richGeoType);
		CbmRichFindRings* richFindRings = new CbmRichFindRings();
		richFindRings->UseFinder(richFinder);
		run->AddTask(richFindRings);
		//--------------------------------------------------------------------------

		//-------------------- RICH Ring Fitting -----------------------------------
		CbmRichRingFitter* richFitter = new CbmRichRingFitterEllipseTau(iVerbose,	1, richGeoType);
		CbmRichFitRings* fitRings = new CbmRichFitRings("", "", richFitter);
		run->AddTask(fitRings);
		//--------------------------------------------------------------------------

		// ------------------- RICH Ring matching  ---------------------------------
		CbmRichMatchRings* matchRings = new CbmRichMatchRings(iVerbose);
		run->AddTask(matchRings);
		// -------------------------------------------------------------------------

		//--------------------- RICH ring-track assignment ------------------------
		Double_t richDistance = 10.; // Max. dist. ring centre to track [cm]
		Int_t richNPoints = 5; // Minmum number of hits on ring
		CbmRichRingTrackAssign* richAssign = new CbmRichRingTrackAssignClosestD(
			richDistance, richNPoints, iVerbose);
		CbmRichAssignTrack* assignTrack = new CbmRichAssignTrack();
		assignTrack->UseAssign(richAssign);
		run->AddTask(assignTrack);
		// ------------------------------------------------------------------------
	}

	// -----   Primary vertex finding   ---------------------------------------
	CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
	CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
	run->AddTask(findVertex);
	// ------------------------------------------------------------------------

	// -----   Track finding QA check   ------------------------------------
	CbmLitReconstructionQa* reconstructionQa = new CbmLitReconstructionQa();
	reconstructionQa->SetMinNofPointsSts(4);
	reconstructionQa->SetMinNofPointsTrd(8);
	reconstructionQa->SetMinNofPointsMuch(0);
	reconstructionQa->SetMinNofPointsTof(1);
	reconstructionQa->SetQuota(0.7);
	reconstructionQa->SetMinNofHitsTrd(8);
	reconstructionQa->SetMinNofHitsMuch(0);
	reconstructionQa->SetMinNofHitsRich(7);
	reconstructionQa->SetQuotaRich(0.6);
	reconstructionQa->SetVerbose(1);
	reconstructionQa->SetMomentumRange(0, 12);
	reconstructionQa->SetNofBinsMom(12);
	reconstructionQa->SetOutputDir(std::string(resultsDir));
	run->AddTask(reconstructionQa);
	// ------------------------------------------------------------------------


	// -----  Parameter database   --------------------------------------------
	FairRuntimeDb* rtdb = run->GetRuntimeDb();
	FairParRootFileIo* parIo1 = new FairParRootFileIo();
	FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
	parIo1->open(parFile.Data());
	parIo2->open(parFileList, "in");
	rtdb->setFirstInput(parIo1);
	rtdb->setSecondInput(parIo2);
	rtdb->setOutput(parIo1);
	rtdb->saveOutput();
	// ------------------------------------------------------------------------

	// -----   Initialize and run   --------------------------------------------
	run->Init();
	run->Run(0, nEvents);
	// ------------------------------------------------------------------------

	// -----   Finish   -------------------------------------------------------
	timer.Stop();
	cout << endl << endl;
 	cout << "Macro finished successfully." << endl;
        cout << "Test passed"<< endl;
        cout << " All ok " << endl;
	cout << "Output file is " << globalRecoFile << endl;
	cout << "Parameter file is " << parFile << endl;
	cout << "Real time " << timer.RealTime() << " s, CPU time "	<< timer.CpuTime() << " s" << endl;
	cout << endl;
	// ------------------------------------------------------------------------
}