void run_reco(Int_t nEvents = 100)
{
   TTree::SetMaxTreeSize(90000000000);

	Int_t iVerbose = 0;

	TString script = TString(gSystem->Getenv("SCRIPT"));
	TString parDir = TString(gSystem->Getenv("VMCWORKDIR")) + TString("/parameters/");

	gRandom->SetSeed(10);

	TString mcFile = "/data/Sim_Outputs/Sim_Thesis/AuAu_10AGeV/setup_align_mc.root";
	TString parFile = "/data/Sim_Outputs/Sim_Thesis/AuAu_10AGeV/setup_align_param.root";
	TString recoFile ="/data/Sim_Outputs/Sim_Thesis/AuAu_10AGeV/setup_align_reco.root";
	TString geoFile = "/data/Sim_Outputs/Sim_Thesis/AuAu_10AGeV/setup_align__geofilefull.root";
	TString setupName = "setup_align";

	TString geoSetupFile = TString(gSystem->Getenv("VMCWORKDIR")) + "/macro/rich/alignment/misalignment_correction/Sim_Thesis/geosetup/setup_align.C";

	std::string resultDir = "cyl_rich_reco_qa/";

	if (script == "yes") {
		mcFile = TString(gSystem->Getenv("MC_FILE"));
		recoFile = TString(gSystem->Getenv("RECO_FILE"));
		parFile = TString(gSystem->Getenv("PAR_FILE"));
		resultDir = TString(gSystem->Getenv("LIT_RESULT_DIR"));
		geoSetupFile = TString(gSystem->Getenv("VMCWORKDIR")) + "/macro/rich/run/geosetup/" + TString(gSystem->Getenv("GEO_SETUP_FILE"));
	}

	remove(recoFile.Data());

    // -----   Load the geometry setup   --------------------------------------
    const char* setupName2 = setupName;
    TString setupFunct = "";
    setupFunct = setupFunct + setupName2 + "()";
    std::cout << "-I- setupFile: " << geoSetupFile << std::endl
	<< "-I- setupFunct: " << setupFunct << std::endl;
    gROOT->LoadMacro(geoSetupFile);
    gROOT->ProcessLine(setupFunct);
    //gInterpreter->ProcessLine(setupFunct);
    // ------------------------------------------------------------------------

    // ----    Debug option   -------------------------------------------------
    gDebug = 0;
    // ------------------------------------------------------------------------

    // -----   Timer   --------------------------------------------------------
    TStopwatch timer;
    timer.Start();
    // ------------------------------------------------------------------------

    // -----   Reconstruction run   -------------------------------------------
    FairRunAna *run= new FairRunAna();
    if (mcFile != "") run->SetInputFile(mcFile);
    if (recoFile != "") run->SetOutputFile(recoFile);
    run->SetGenerateRunInfo(kTRUE);
    run->SetGeomFile(geoFile);
    // ------------------------------------------------------------------------

    // ----- MC Data Manager   ------------------------------------------------
    CbmMCDataManager* mcManager=new CbmMCDataManager;
    mcManager->AddFile(mcFile);
    run->AddTask(mcManager);
    // ------------------------------------------------------------------------
    
    Bool_t isRich = IsRich(parFile);
    
	// =========================================================================
	// ===                      STS local reconstruction                     ===
	// =========================================================================
    Double_t dynRange       =   40960.;  // Dynamic range [e]
    Double_t threshold      =    4000.;  // Digitisation threshold [e]
    Int_t nAdc              =    4096;   // Number of ADC channels (12 bit)
    Double_t timeResolution =       5.;  // time resolution [ns]
    Double_t deadTime       = 9999999.;  // infinite dead time (integrate entire event)
    Double_t noise          =       0.;  // ENC [e]
    Int_t digiModel         =       1;   // User sensor type DSSD
    
    // The following settings correspond to a validated implementation.
    // Changing them is on your own risk.
    Int_t  eLossModel       = 1;         // Energy loss model: uniform
    Bool_t useLorentzShift  = kFALSE;    // Deactivate Lorentz shift
    Bool_t useDiffusion     = kFALSE;    // Deactivate diffusion
    Bool_t useCrossTalk     = kFALSE;    // Deactivate cross talk
    
    CbmStsDigitize* stsDigi = new CbmStsDigitize(digiModel);
    stsDigi->SetProcesses(eLossModel, useLorentzShift, useDiffusion, useCrossTalk);
    stsDigi->SetParameters(dynRange, threshold, nAdc, timeResolution, deadTime, noise);
    run->AddTask(stsDigi);


	FairTask* stsClusterFinder = new CbmStsFindClusters();
	run->AddTask(stsClusterFinder);

	FairTask* stsFindHits = new CbmStsFindHits();
	run->AddTask(stsFindHits);

	CbmKF* kalman = new CbmKF();
	run->AddTask(kalman);
	CbmL1* l1 = new CbmL1();
	l1->SetStsMaterialBudgetFileName(stsMatBudgetFileName.Data());
	if (mvdMatBudgetFileName != "") l1->SetMvdMaterialBudgetFileName(mvdMatBudgetFileName.Data());
	run->AddTask(l1);

	CbmStsTrackFinder* stsTrackFinder = new CbmL1StsTrackFinder();
	FairTask* stsFindTracks = new CbmStsFindTracks(1, stsTrackFinder);
	run->AddTask(stsFindTracks);

	// =========================================================================
	// ===                        Global tracking                            ===
	// =========================================================================

	CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks();
	finder->SetTrackingType(std::string("branch"));
	finder->SetMergerType("nearest_hit");
	run->AddTask(finder);

	CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
	CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
	run->AddTask(findVertex);


	if (isTrd) {
		CbmTrdSetTracksPidANN* trdSetTracksPidAnnTask = new CbmTrdSetTracksPidANN("CbmTrdSetTracksPidANN","CbmTrdSetTracksPidANN");
		trdSetTracksPidAnnTask->SetTRDGeometryType("h++");
		run->AddTask(trdSetTracksPidAnnTask);
	}//isTrd

    // =========================================================================
    // ===                        RICH reconstruction                        ===
    // =========================================================================
	if (isRich){
		CbmRichDigitizer* richDigitizer = new CbmRichDigitizer();
        //richDigitizer->SetNofNoiseHits(0);
		run->AddTask(richDigitizer);

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

		CbmRichReconstruction* richReco = new CbmRichReconstruction();
       // richReco->SetRunExtrapolation(true);
       // richReco->SetRunProjection(true);
       // richReco->SetRunTrackAssign(true);
        //richReco->SetFinderName("ideal");
        //richReco->SetProjectionName("tgeo");
       // richReco->SetFitterName("circle_cop");;
		run->AddTask(richReco);
	}//isRich


	CbmMatchRecoToMC* matchRecoToMc = new CbmMatchRecoToMC();
	run->AddTask(matchRecoToMc);
    
    // RICH reco QA
    CbmRichRecoQa* richRecoQa = new CbmRichRecoQa();
    richRecoQa->SetOutputDir(resultDir);
    run->AddTask(richRecoQa);

	// Reconstruction Qa
	CbmLitTrackingQa* trackingQa = new CbmLitTrackingQa();
	trackingQa->SetMinNofPointsSts(4);
	trackingQa->SetUseConsecutivePointsInSts(true);
	trackingQa->SetMinNofPointsTrd(minNofPointsTrd);
	trackingQa->SetMinNofPointsMuch(10);
	trackingQa->SetMinNofPointsTof(1);
	trackingQa->SetQuota(0.7);
	trackingQa->SetMinNofHitsTrd(minNofPointsTrd);
	trackingQa->SetMinNofHitsMuch(10);
	trackingQa->SetVerbose(0);
	trackingQa->SetMinNofHitsRich(7);
	trackingQa->SetQuotaRich(0.6);
	trackingQa->SetOutputDir(resultDir);
	trackingQa->SetPRange(20, 0., 10.);
	trackingQa->SetTrdAnnCut(trdAnnCut);
	std::vector<std::string> trackCat, richCat;
	trackCat.push_back("All");
	trackCat.push_back("Electron");
	richCat.push_back("Electron");
	richCat.push_back("ElectronReference");
	trackingQa->SetTrackCategories(trackCat);
	trackingQa->SetRingCategories(richCat);
	run->AddTask(trackingQa);

	CbmLitFitQa* fitQa = new CbmLitFitQa();
	fitQa->SetMvdMinNofHits(0);
	fitQa->SetStsMinNofHits(6);
	fitQa->SetMuchMinNofHits(10);
	fitQa->SetTrdMinNofHits(minNofPointsTrd);
	fitQa->SetOutputDir(resultDir);
	//run->AddTask(fitQa);

	CbmLitClusteringQa* clusteringQa = new CbmLitClusteringQa();
	clusteringQa->SetOutputDir(resultDir);
	//run->AddTask(clusteringQa);

	//CbmLitTofQa* tofQa = new CbmLitTofQa();
	//tofQa->SetOutputDir(std::string(resultDir));
	//run->AddTask(tofQa);

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

    run->Init();
    cout << "Starting run" << endl;
    run->Run(0,nEvents);

    // -----   Finish   -------------------------------------------------------
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished successfully." << endl;
    cout << "Output file is "    << recoFile << endl;
    cout << "Parameter file is " << parFile << endl;
    cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
    cout << endl;

    cout << " Test passed" << endl;
    cout << " All ok " << endl;
}