// -----------------------------------------------------------------------------
// ----- ana_hits.C                                                        -----
// -----                                                                   -----
// ----- adapted by C. Simon on 2018-06-04                                 -----
// -----                                                                   -----
// -----------------------------------------------------------------------------

void ana_hits(Int_t iNEvents, Int_t iGenCor, const TString& tUnpackDir, const TString& tCluDir,
              const TString& tInputFilePath, Int_t iNDigiFiles, const TString& tMCLinkFileName, Int_t iNLinkFiles,
              const TString& tTofGeoVersion, const TString& tSetupName, Bool_t bIdealCalibration,
              Bool_t bHeavyIonCollisions, Bool_t bIdealClustering, Bool_t bSeparateHitFile, Bool_t bAnalysisMCQA,
              Int_t iDut, Int_t iMRef, Int_t iBRef, Int_t iSel2, Double_t dScaleFactor, Bool_t bSimData,
              Int_t iCalibrationMode, Bool_t bFinalRun, Double_t dTimeInSpill, Int_t iNTISBins,
              Double_t dExtChi2Lim = 1000000., Double_t dExtChi2Lim2 = 1000000., Int_t iExtSigMode = 0,
              Double_t dExtSigX = 1., Double_t dExtSigY = 1., Double_t dExtSigT = 0.1,
              Bool_t bDifferentialTISAnalysis = kFALSE)
{
  // ---------------------------------------------------------------------------

  TString tLogLevel = "ERROR";
  TString tLogVerbosity = "LOW";
  gErrorIgnoreLevel = kWarning;
  gDebug = 0;

  TTree::SetMaxTreeSize(1000000000000LL); // [B] here: 1 TB

  // ---------------------------------------------------------------------------

  TStopwatch tTimer;
  tTimer.Start();

  // ---------------------------------------------------------------------------

  FairLogger::GetLogger()->SetLogScreenLevel(tLogLevel.Data());
  FairLogger::GetLogger()->SetLogVerbosityLevel(tLogVerbosity.Data());

  // ---------------------------------------------------------------------------

  TString tSrcDir   = gSystem->Getenv("VMCWORKDIR");
  TString tMacroDir = tSrcDir + "/macro/beamtime";
  TString tParamDir = tSrcDir + "/parameters";
  TString tGeoDir   = tSrcDir + "/geometry";

  TString tUnpackInputFile      = tUnpackDir + "/data/unpack.out.root";
  TString tCluInputFile         = tCluDir + "/data/digi.out.root";
  TString tOutputFile           = "./hits.out.root";
  TString tHistoFile            = "./ana_cluster.hst.root";
  TString tInputAnaCalibFile    = "./calib_ana.cor_in.root";
  TString tOutputAnaCalibFile   = "./calib_ana.cor_out.root";

  // ---------------------------------------------------------------------------

  TList* tParFileList = new TList();

  TObjString* tTofDigiParFile    = new TObjString(tParamDir + "/tof/tof_"+ tTofGeoVersion + ".digi.par");
  TObjString* tTofDigiBdfParFile = new TObjString(tParamDir + "/tof/tof_"+ tTofGeoVersion + ".digibdf.par");

  tParFileList->Add(tTofDigiParFile);
  tParFileList->Add(tTofDigiBdfParFile);

  // ---------------------------------------------------------------------------
  if(bSimData)
  {
    if(bHeavyIonCollisions)
    {
      TFile::Open(tInputFilePath + "/001/data/mc/00001/data/" + tSetupName + "_geofile_full.root");
    }
    else
    {
      TFile::Open(tInputFilePath + "/001/data/" + tSetupName + "_geofile_full.root");
    }
  }
  else
  {
    TFile::Open(tGeoDir + "/tof/geofile_tof_" + tTofGeoVersion + ".root");
  }

  gFile->Get("FAIRGeom");
  if(!gGeoManager)
  {
    std::cout << "-E- ana_hits: FAIRGeom not found in geometry file!" << std::endl;
    return;
  }

  // ---------------------------------------------------------------------------

  FairFileSource* tFileSource = new FairFileSource(tUnpackInputFile);
  if(bSeparateHitFile)
  {
    tFileSource->AddFriend(tCluInputFile);
  }

  // ---------------------------------------------------------------------------

  FairRunAna* tRun= new FairRunAna();
  tRun->SetContainerStatic(kTRUE);
  tRun->SetSource(tFileSource);
  tRun->SetOutputFile(tOutputFile);

/*
  FairMonitor::GetMonitor()->EnableMonitor(kTRUE);
*/
  // ---------------------------------------------------------------------------

  Int_t iSel = 1000*iDut + iMRef;

  Int_t iDutCounter = iDut%10;
  Int_t iDutModule  = ((iDut - iDutCounter)/10)%10;
  Int_t iDutType    = (iDut - iDutCounter - 10*iDutModule)/100;

  Int_t iMRefCounter = iMRef%10;
  Int_t iMRefModule  = ((iMRef - iMRefCounter)/10)%10;
  Int_t iMRefType    = (iMRef - iMRefCounter - 10*iMRefModule)/100;

  Int_t iBRefCounter = iBRef%10;
  Int_t iBRefModule  = ((iBRef - iBRefCounter)/10)%10;
  Int_t iBRefType    = (iBRef - iBRefCounter - 10*iBRefModule)/100;

  Int_t iSel2Counter = iSel2%10;
  Int_t iSel2Module  = ((iSel2 - iSel2Counter)/10)%10;
  Int_t iSel2Type    = (iSel2 - iSel2Counter - 10*iSel2Module)/100;

  // ---------------------------------------------------------------------------

  CbmMCDataManager* tDataManager = new CbmMCDataManager("MCManager");

  if(bSimData)
  {
    if(bHeavyIonCollisions)
    {
      for(Int_t iDigiFile = 1; iDigiFile <= iNDigiFiles; iDigiFile++)
      {
        for(Int_t iLinkFile = 1; iLinkFile <= iNLinkFiles; iLinkFile++)
        {
          tDataManager->AddFileToChain(tInputFilePath + TString::Format("/%03d/data", iDigiFile)
                                                      + TString::Format("/mc/%05d/data/", iLinkFile)
                                                      + tMCLinkFileName,
                                       iDigiFile - 1);
        }
      }
    }
    else
    {
      for(Int_t iDigiFile = 1; iDigiFile <= iNDigiFiles; iDigiFile++)
      {
        tDataManager->AddFileToChain(tInputFilePath + TString::Format("/%03d/data/", iDigiFile)
                                                    + tMCLinkFileName,
                                     iDigiFile - 1);
      }
    }
  }

  // ---------------------------------------------------------------------------

  CbmTofMatchReco* tTofMatch = new CbmTofMatchReco();
  tTofMatch->SetAlternativeBranchNames(kFALSE);
  tTofMatch->SetFileIndex(0);
  tTofMatch->SetMatchMCTracks(kTRUE);

  // ---------------------------------------------------------------------------

  CbmTofAnaTestbeam* tTofAnaTestbeam = new CbmTofAnaTestbeam("TOF TestBeam Analysis", 1);
  tTofAnaTestbeam->SetDut(iDutType);
  tTofAnaTestbeam->SetDutSm(iDutModule);
  tTofAnaTestbeam->SetDutRpc(iDutCounter);
  tTofAnaTestbeam->SetMrpcRef(iMRefType);
  tTofAnaTestbeam->SetMrpcRefSm(iMRefModule);
  tTofAnaTestbeam->SetMrpcRefRpc(iMRefCounter);
  if(0 != iSel2Type)
  {
    tTofAnaTestbeam->SetMrpcSel2(iSel2Type);
    tTofAnaTestbeam->SetMrpcSel2Sm(iSel2Module);
    tTofAnaTestbeam->SetMrpcSel2Rpc(iSel2Counter);
  }
  tTofAnaTestbeam->SetBeamRefSmType(iBRefType);
  tTofAnaTestbeam->SetBeamRefSmId(iBRefModule);
  tTofAnaTestbeam->SetBeamRefRpc(iBRefCounter);

  tTofAnaTestbeam->SetCorMode(iGenCor);
//  if(!bIdealCalibration)
  {
    tTofAnaTestbeam->SetCalParFileName(tInputAnaCalibFile);
  }
  tTofAnaTestbeam->SetCalOutFileName(tOutputAnaCalibFile);
  tTofAnaTestbeam->SetMonteCarloComparison(bAnalysisMCQA);
  tTofAnaTestbeam->SetMinMCRefTrackPoints(4);
  tTofAnaTestbeam->SetDelayMCPoints(!bIdealClustering);
  if(bFinalRun)
  {
    tTofAnaTestbeam->SetUseSigCalib(kTRUE);
  }

  tTofAnaTestbeam->SetSpillDuration(dTimeInSpill);
  tTofAnaTestbeam->SetNTISBins(iNTISBins);
  tTofAnaTestbeam->SetFitTISMeanDependency(kTRUE); // TODO: works?
  tTofAnaTestbeam->SetCorrectTISDependency(bDifferentialTISAnalysis);
  tTofAnaTestbeam->SetCorrectChi2LimitedResiduals(kFALSE);
//  tTofAnaTestbeam->SetSimData(bSimData); // currently no effect
  tTofAnaTestbeam->SetIgnoreStartTimeMUL(kTRUE);
  tTofAnaTestbeam->SetPositionInterpolation(kTRUE);

  tTofAnaTestbeam->SetExtSigMode(iExtSigMode);
  tTofAnaTestbeam->SetChi2Lim(dExtChi2Lim/3.);
  tTofAnaTestbeam->SetChi2Lim2(dExtChi2Lim2/3.);

  tTofAnaTestbeam->SetMulDMax(10);
  tTofAnaTestbeam->SetMul0Max(50);
  tTofAnaTestbeam->SetMul4Max(50);
  tTofAnaTestbeam->SetMulS2Max(dScaleFactor); // TODO: workaround to pass an external MUL(Sel2) limit to this macro
/*
// TODO: allows for additional suppression of multi-track events
         if MUL cuts on Sel2 and MRef do not suffice
  tTofAnaTestbeam->SetMrpcSel3(9);
  tTofAnaTestbeam->SetMrpcSel3Sm(0);
  tTofAnaTestbeam->SetMrpcSel3Rpc(0);
  tTofAnaTestbeam->SetMulS3Max(1);
//
*/
  if(10 == iGenCor)
  {
    tTofAnaTestbeam->SetDTD4MAX(100.);
  }
  else
  {
    tTofAnaTestbeam->SetDTD4MAX(6.);
  }

  tTofAnaTestbeam->SetDTDia(0.);
  tTofAnaTestbeam->SetHitDistMin(30.);
  tTofAnaTestbeam->SetHitDistAv(0.);
  tTofAnaTestbeam->SetTShift(0.);
  tTofAnaTestbeam->SetTOffD4(0.);

  tTofAnaTestbeam->SetCh4Sel(0.);
  tTofAnaTestbeam->SetDCh4Sel(100.);
  tTofAnaTestbeam->SetPosY4SelOff(0.);
  tTofAnaTestbeam->SetPosY4Sel(10.);
//  tTofAnaTestbeam->SetPosY4Sel(0.5*dScaleFactor);

  tTofAnaTestbeam->SetChS2Sel(0.);
  tTofAnaTestbeam->SetDChS2Sel(100.);
  tTofAnaTestbeam->SetPosYS2SelOff(0.);
  tTofAnaTestbeam->SetPosYS2Sel(10.);

  tTofAnaTestbeam->SetDXMean(0.);
  tTofAnaTestbeam->SetDYMean(0.);
  tTofAnaTestbeam->SetDTMean(0.);

  tTofAnaTestbeam->SetDXWidth(1.);
  tTofAnaTestbeam->SetDYWidth(1.);
  tTofAnaTestbeam->SetDTWidth(0.1);

  tTofAnaTestbeam->SetSel2XOff(0.);
  tTofAnaTestbeam->SetSel2YOff(0.);
  tTofAnaTestbeam->SetSel2TOff(0.);

  tTofAnaTestbeam->SetDXWidthSel24(1.);
  tTofAnaTestbeam->SetDYWidthSel24(1.);
  tTofAnaTestbeam->SetDTWidthSel24(0.1);

  // ---------------------------------------------------------------------------

  if(iExtSigMode)
  {
    // DUT - MRef
    if(1 == iExtSigMode)
    {
      tTofAnaTestbeam->SetDXWidth(dExtSigX);
      tTofAnaTestbeam->SetDYWidth(dExtSigY);
      tTofAnaTestbeam->SetDTWidth(dExtSigT);
    }
    // Sel2 - MRef
    if(2 == iExtSigMode)
    {
      tTofAnaTestbeam->SetDXWidthSel24(dExtSigX);
      tTofAnaTestbeam->SetDYWidthSel24(dExtSigY);
      tTofAnaTestbeam->SetDTWidthSel24(dExtSigT);
    }
  }

  // ---------------------------------------------------------------------------

  Double_t dFitRangeDX04(1.);
  Double_t dFitRangeDY04(1.);
  Double_t dFitRangeDT04(0.15);

  if(bSimData)
  {
    dFitRangeDX04 = 1.2;

    if(700 == iDut)
    {
      dFitRangeDY04 = 1.2;
    }
  }
  else
  {
    if(700 == iDut)
    {
      dFitRangeDX04 = 2.;
      dFitRangeDY04 = 1.3;
    }
  }

  tTofAnaTestbeam->SetFitRangeDXSel24(1.);
  tTofAnaTestbeam->SetFitRangeDYSel24(1.);
  tTofAnaTestbeam->SetFitRangeDTSel24(0.15);
  tTofAnaTestbeam->SetFitRangeDX04(dFitRangeDX04);
  tTofAnaTestbeam->SetFitRangeDY04(dFitRangeDY04);
  tTofAnaTestbeam->SetFitRangeDT04(dFitRangeDT04);
  tTofAnaTestbeam->SetCorRangeLimitDX04(-5., 5.);
  tTofAnaTestbeam->SetCorRangeLimitDY04(-5., 5.);
  tTofAnaTestbeam->SetCorRangeLimitDT04(-0.2, 0.2);

  // ---------------------------------------------------------------------------

  if(bSimData)
  {
    tRun->AddTask(tDataManager);
    if(bAnalysisMCQA)
    {
      tRun->AddTask(tTofMatch);
    }
  }

  tRun->AddTask(tTofAnaTestbeam);

  // ---------------------------------------------------------------------------

  FairRuntimeDb* tRuntimeDb = tRun->GetRuntimeDb();
  FairParAsciiFileIo* tParIo1 = new FairParAsciiFileIo();
  tParIo1->open(tParFileList, "in");
  tRuntimeDb->setFirstInput(tParIo1);

  // ---------------------------------------------------------------------------

  std::cout << "-I- ana_hits: initializing run..." << std::endl;
  tRun->Init();

  std::cout << "-I- ana_hits: starting run..." << std::endl;
  tRun->Run(0, iNEvents);

  // ---------------------------------------------------------------------------

  gROOT->LoadMacro((tMacroDir + "/save_hst.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_MatD4sel.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_Mat04D4best.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_TIS.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_eff_XY.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_MCEvent.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_MCHitSel.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_MCCounter.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_MCHitRes.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_SelTypes.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_all_SelTypes.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_TIS.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_CS.C").Data());
  gROOT->LoadMacro((tMacroDir + "/pl_over_CluArea.C").Data());

  gInterpreter->ProcessLine("save_hst(\"" + tHistoFile + "\", kTRUE)");

  gInterpreter->ProcessLine("pl_over_MatD4sel()");
  gInterpreter->ProcessLine("pl_over_Mat04D4best(1)");
  gInterpreter->ProcessLine("pl_TIS()");
  gInterpreter->ProcessLine("pl_eff_XY()");

  if(bAnalysisMCQA)
  {
    Double_t dTyMean = 66.8/367.4; // FIXME: hard-coded

    gInterpreter->ProcessLine(Form("pl_over_MCEvent(%f)", dTyMean));
    gInterpreter->ProcessLine("pl_over_MCHitSel()");
    gInterpreter->ProcessLine(Form("pl_over_MCHitRes(%f, %f, %f, %f)",
                                   dExtChi2Lim, dFitRangeDX04, dFitRangeDY04, dFitRangeDT04));
    gInterpreter->ProcessLine("pl_over_MCCounter(5, 2, 0)");
    gInterpreter->ProcessLine("pl_over_MCCounter(5, 1, 0)");
    gInterpreter->ProcessLine("pl_over_MCCounter(3, 0, 0)");
    gInterpreter->ProcessLine("pl_over_MCCounter(9, 0, 0)");
    gInterpreter->ProcessLine("pl_over_MCCounter(7, 0, 0)");
    gInterpreter->ProcessLine("pl_over_MCCounter(4, 0, 0)");
  }

  gInterpreter->ProcessLine("pl_over_SelTypes()");
  gInterpreter->ProcessLine("pl_all_SelTypes()");
  gInterpreter->ProcessLine(Form("pl_over_TIS(%f, %f, %f, %d, %d, %f, %f, %f)",
                                 dExtChi2Lim, dExtChi2Lim2,
                                 dTimeInSpill, iNTISBins, bDifferentialTISAnalysis,
                                 dFitRangeDX04, dFitRangeDY04, dFitRangeDT04));
  gInterpreter->ProcessLine(Form("pl_over_CS(%f, %f, %f, %f)",
                                 dExtChi2Lim, dFitRangeDX04, dFitRangeDY04, dFitRangeDT04));
  gInterpreter->ProcessLine("pl_over_CluArea()");

  // ---------------------------------------------------------------------------

  tTimer.Stop();
  Double_t dRealTime = tTimer.RealTime();
  Double_t dCPUTime = tTimer.CpuTime();

  std::cout << std::endl << std::endl;
  std::cout << "Macro finished successfully." << std::endl;
  std::cout << "Real time " << dRealTime << " s, CPU time " << dCPUTime << " s" << std::endl;
  std::cout << std::endl;

  // ---------------------------------------------------------------------------

  // Extract the maximal used memory an add is as Dart measurement
  // This line is filtered by CTest and the value send to CDash
  FairSystemInfo tSysInfo;
  Float_t fMaxMemory = tSysInfo.GetMaxMemory();
  std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
  std::cout << fMaxMemory;
  std::cout << "</DartMeasurement>" << endl;

  Float_t fCPUUsage = dCPUTime/dRealTime;
  std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
  std::cout << fCPUUsage;
  std::cout << "</DartMeasurement>" << endl;
/*
  FairMonitor::GetMonitor()->Print();
*/
  // ---------------------------------------------------------------------------

  FairRootManager::Instance()->CloseSink();

  // ---------------------------------------------------------------------------

  RemoveGeoManager();

  // ---------------------------------------------------------------------------
}