/**
 * \file global_reco.C
 * \brief Event reconstruction macro.
 * \author Andrey Lebedev <andrey.lebedev@gsi.de>
 * \date 2010
 *
 * Macro runs event reconstruction in MVD, STS, TRD, MUCH and TOF.
 * Macro has 3 options "all", "hits" and "tracking".
 */

#include <iostream>
using std::cout;
using std::endl;

void global_reco(Int_t nEvents = 10000, // number of events
		TString opt = "all")
// if opt == "all" STS + hit producers + global tracking are executed
// if opt == "hits" STS + hit producers are executed
// if opt == "tracking" global tracking is executed
{
   TTree::SetMaxTreeSize(90000000000);
	TString script = TString(gSystem->Getenv("LIT_SCRIPT"));
	TString parDir = TString(gSystem->Getenv("VMCWORKDIR")) + TString("/parameters");

   // Input and output data
	TString dir = "events/much_anna_omega_8gev_10k/"; // Output directory
   TString mcFile = dir + "mc.0000.root"; // MC transport file
   TString parFile = dir + "param.0000.root"; // Parameters file
   TString globalRecoFile = dir + "global.reco.0000.root"; // Output file with reconstructed tracks and hits
   TString globalHitsFile = dir + "global.hits.0000.root"; // File with reconstructed STS tracks, STS, MUCH, TRD and TOF hits and digis
   TString globalTracksFile = dir + "global.tracks.0000.root";// Output file with global tracks

   // Digi files
   TList* parFileList = new TList();
   TObjString stsDigiFile = parDir + "/sts/sts_v12b_std.digi.par"; // STS digi file
   TObjString trdDigiFile = parDir + "/trd/trd_v13g.digi.par"; // TRD digi file
   TString muchDigiFile = parDir + "/much/much_v12c.digi.root"; // MUCH digi file
   TString stsMatBudgetFile = parDir + "/sts/sts_matbudget_v12b.root";
   TObjString tofDigiFile = parDir + "/tof/tof_v13b.digi.par";// TOF digi file

//   // Directory for output results
//   TString resultDir = "./test/";

   // Reconstruction parameters
   TString globalTrackingType = "branch"; // Global tracking type
   TString stsHitProducerType = "real"; // STS hit producer type: real, ideal
   TString trdHitProducerType = "smearing"; // TRD hit producer type: smearing, digi, clustering
   TString muchHitProducerType = "advanced"; // MUCH hit producer type: simple, advanced

//   // Normalization for efficiency
//   Int_t normStsPoints = 4;
//   Int_t normTrdPoints = 6;
//   Int_t normMuchPoints = 16;
//   Int_t normTofPoints = 1;
//   Int_t normTrdHits = 6;
//   Int_t normMuchHits = 16;
//   Int_t normTofHits = 1;

	if (script == "yes") {
		mcFile = TString(gSystem->Getenv("LIT_MC_FILE"));
		parFile = TString(gSystem->Getenv("LIT_PAR_FILE"));
		globalRecoFile = TString(gSystem->Getenv("LIT_GLOBAL_RECO_FILE"));
		globalHitsFile = TString(gSystem->Getenv("LIT_GLOBAL_HITS_FILE"));
		globalTracksFile = TString(gSystem->Getenv("LIT_GLOBAL_TRACKS_FILE"));

//		resultDir = TString(gSystem->Getenv("LIT_RESULT_DIR"));

		globalTrackingType = TString(gSystem->Getenv("LIT_GLOBAL_TRACKING_TYPE"));
		stsHitProducerType = TString(gSystem->Getenv("LIT_STS_HITPRODUCER_TYPE"));
		trdHitProducerType = TString(gSystem->Getenv("LIT_TRD_HITPRODUCER_TYPE"));
		muchHitProducerType = TString(gSystem->Getenv("LIT_MUCH_HITPRODUCER_TYPE"));

		stsDigiFile = TString(gSystem->Getenv("LIT_STS_DIGI"));
		trdDigiFile = TString(gSystem->Getenv("LIT_TRD_DIGI"));
		muchDigiFile = TString(gSystem->Getenv("LIT_MUCH_DIGI"));
		tofDigiFile = TString(gSystem->Getenv("LIT_TOF_DIGI"));
		stsMatBudgetFile = TString(gSystem->Getenv("LIT_STS_MAT_BUDGET_FILE"));

//		normStsPoints = TString(gSystem->Getenv("LIT_NORM_STS_POINTS")).Atoi();
//		normTrdPoints = TString(gSystem->Getenv("LIT_NORM_TRD_POINTS")).Atoi();
//		normMuchPoints = TString(gSystem->Getenv("LIT_NORM_MUCH_POINTS")).Atoi();
//		normTofPoints = TString(gSystem->Getenv("LIT_NORM_TOF_POINTS")).Atoi();
//		normTrdHits = TString(gSystem->Getenv("LIT_NORM_TRD_HITS")).Atoi();
//		normMuchHits = TString(gSystem->Getenv("LIT_NORM_MUCH_HITS")).Atoi();
//		normTofHits = TString(gSystem->Getenv("LIT_NORM_TOF_HITS")).Atoi();
	}

   parFileList->Add(&stsDigiFile);
   parFileList->Add(&trdDigiFile);
   parFileList->Add(&tofDigiFile);

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

	gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/loadlibs.C");
	loadlibs();

	FairRunAna *run = new FairRunAna();
	if (opt == "all") {
		run->SetInputFile(mcFile);
		run->SetOutputFile(globalRecoFile);
	} else if (opt == "hits") {
		run->SetInputFile(mcFile);
		run->SetOutputFile(globalHitsFile);
	} else if (opt == "tracking") {
		run->SetInputFile(mcFile);
		run->AddFriend(globalHitsFile);
		run->SetOutputFile(globalTracksFile);
	} else {
		std::cout << "-E- Incorrect opt parameter" << std::endl;
		exit(0);
	}

	if (opt == "all" || opt == "hits") {

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

		if (stsHitProducerType == "real") {
		// ----- STS REAL reconstruction -----------------------------------------------
         // -----   STS digitizer   -------------------------------------------------
         Double_t threshold = 4;
         Double_t noiseWidth = 0.01;
         Int_t nofBits = 12;
         Double_t ElectronsPerAdc =  10.;
         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->SetFrontNofElPerAdc(ElectronsPerAdc);
         stsDigitize->SetBackNofElPerAdc(ElectronsPerAdc);
         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);

		} else if (stsHitProducerType == "ideal") {

         // ----- STS IDEAL reconstruction   ---------------------------------------------
         FairTask* stsDigitize =	new CbmStsIdealDigitize("STSDigitize", iVerbose);
         run->AddTask(stsDigitize);

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

         FairTask* stsFindHits =	new CbmStsIdealFindHits("STSFindHits", iVerbose);
         run->AddTask(stsFindHits);

         FairTask* stsMatchHits = new CbmStsIdealMatchHits("STSMatchHits", iVerbose);
         run->AddTask(stsMatchHits);
		}

		FairTask* kalman = new CbmKF();
		run->AddTask(kalman);
		CbmL1* l1 = new CbmL1();
		//l1->SetExtrapolateToTheEndOfSTS(true);
		l1->SetMaterialBudgetFileName(stsMatBudgetFile);
		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);
		// ------------------------------------------------------------------------

		if (IsMuch(parFile)) {
			// -------- MUCH digitization ------------
			CbmMuchDigitizeGem* digitize = new CbmMuchDigitizeGem(muchDigiFile.Data());
			if (muchHitProducerType == "simple") {
				digitize->SetAlgorithm(0);
			} else if (muchHitProducerType == "advanced") {
				digitize->SetAlgorithm(1);
			}
			run->AddTask(digitize);
			CbmMuchFindHitsGem* findHits = new CbmMuchFindHitsGem(muchDigiFile.Data());
			run->AddTask(findHits);

		//	CbmMuchDigitizeStraws* strawDigitize = new CbmMuchDigitizeStraws(
		//	      "MuchDigitizeStraws", muchDigiFile.Data(), iVerbose);
		//	run->AddTask(strawDigitize);
		//	CbmMuchFindHitsStraws* strawFindHits = new CbmMuchFindHitsStraws(
		//			"MuchFindHitsStraws", muchDigiFile.Data(), iVerbose);
		//	strawFindHits->SetMerge(1);
		//	run->AddTask(strawFindHits);
			// -----------------------------------------------------------------
		}

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

			if (trdHitProducerType == "smearing") {
				CbmTrdHitProducerSmearing* trdHitProd = new CbmTrdHitProducerSmearing(radiator);
				trdHitProd->SetUseDigiPar(false);
				run->AddTask(trdHitProd);
			} else if (trdHitProducerType == "digi") {
				CbmTrdDigitizer* trdDigitizer = new CbmTrdDigitizer(radiator);
				run->AddTask(trdDigitizer);

				CbmTrdHitProducerDigi* trdHitProd = new CbmTrdHitProducerDigi();
				run->AddTask(trdHitProd);
			} else if (trdHitProducerType == "clustering") {
				// ----- TRD clustering -----
				CbmTrdDigitizerPRF* trdClustering = new CbmTrdDigitizerPRF("TRD Clusterizer", "TRD task", radiator, false, true);
				run->AddTask(trdClustering);

				CbmTrdClusterFinderFast* trdClusterfindingfast = new CbmTrdClusterFinderFast(true, true, false, 5.0e-7);
				run->AddTask(trdClusterfindingfast);

				CbmTrdHitProducerCluster* trdClusterHitProducer = new CbmTrdHitProducerCluster();
				run->AddTask(trdClusterHitProducer);
				// ----- End TRD Clustering -----
			}
			// ------------------------------------------------------------------------
		}

		if (IsTof(parFile)) {
			// ------ TOF hits --------------------------------------------------------
		   CbmTofHitProducerNew* tofHitProd = new CbmTofHitProducerNew("TOF HitProducerNew",iVerbose);
		   tofHitProd->SetInitFromAscii(kFALSE);
		   run->AddTask(tofHitProd);
			// ------------------------------------------------------------------------
		}
	}

	if (opt == "all" || opt == "tracking") {
		// ------ Global track reconstruction -------------------------------------
		CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks();
		//CbmLitFindGlobalTracksParallel* finder = new CbmLitFindGlobalTracksParallel();
		// Tracking method to be used
		// "branch" - branching tracking
		// "nn" - nearest neighbor tracking
		// "nn_parallel" - nearest neighbor parallel tracking
		finder->SetTrackingType(std::string(globalTrackingType));

		// 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();
			run->AddTask(trdMatchTracks);
		}

		if (IsMuch(parFile)) {
			CbmMuchMatchTracks* muchMatchTracks = new CbmMuchMatchTracks();
			run->AddTask(muchMatchTracks);
		}
	}

	if (opt == "all" || opt == "tracking") {
	  if (IsRich(parFile)) {
	     if (opt == "tracking") {
	        FairTask* kalman = new CbmKF();
	        run->AddTask(kalman);
	        FairTask* l1 = new CbmL1();
	        run->AddTask(l1);
	     }

	     CbmRichHitProducer* richHitProd  = new CbmRichHitProducer();
	     run->AddTask(richHitProd);

		  CbmRichReconstruction* richReco = new CbmRichReconstruction();
		  run->AddTask(richReco);

		  CbmRichMatchRings* matchRings = new CbmRichMatchRings();
		  run->AddTask(matchRings);
	   }
    }

    if (opt == "all" || opt == "tracking") {
      if (opt == "tracking") {
         FairTask* kalman = new CbmKF();
         run->AddTask(kalman);
         FairTask* l1 = new CbmL1();
         run->AddTask(l1);
      }
		// -----   Primary vertex finding   --------------------------------------
		CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
		CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
		run->AddTask(findVertex);
		// -----------------------------------------------------------------------
	}

//	if (opt == "all" || opt == "tracking") {
//		// ----- Reconstruction QA tasks -----------------------------------------
//		CbmLitTrackingQa* trackingQa = new CbmLitTrackingQa();
//		trackingQa->SetMinNofPointsSts(normStsPoints);
//		trackingQa->SetMinNofPointsTrd(normTrdPoints);
//		trackingQa->SetMinNofPointsMuch(normMuchPoints);
//		trackingQa->SetMinNofPointsTof(normTofPoints);
//		trackingQa->SetQuota(0.7);
//		trackingQa->SetMinNofHitsTrd(normTrdHits);
//		trackingQa->SetMinNofHitsMuch(normMuchHits);
//		trackingQa->SetMinNofHitsRich(7);
//		trackingQa->SetQuotaRich(0.6);
//		trackingQa->SetVerbose(normTofHits);
//		trackingQa->SetOutputDir(std::string(resultDir));
//		run->AddTask(trackingQa);
//
//		CbmLitFitQa* fitQa = new CbmLitFitQa();
//		fitQa->SetMvdMinNofHits(0);
//		fitQa->SetStsMinNofHits(normStsPoints);
//		fitQa->SetMuchMinNofHits(normMuchPoints);
//		fitQa->SetTrdMinNofHits(normTrdPoints);
//		fitQa->SetOutputDir(std::string(resultDir));
//		run->AddTask(fitQa);
//
//		CbmLitClusteringQa* clusteringQa = new CbmLitClusteringQa();
//		clusteringQa->SetOutputDir(std::string(resultDir));
//		run->AddTask(clusteringQa);
//		// -----------------------------------------------------------------------
//	}

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