// -----------------------------------------------------------------------------
// ----- ana_digi_qa.C                                                     -----
// -----                                                                   -----
// ----- adapted by C. Simon on 2018-07-27                                 -----
// -----                                                                   -----
// -----------------------------------------------------------------------------

void ana_digi_qa(Int_t iNEvents, Int_t iCalMode, Int_t iCalSel, Int_t iCalRPC,
                 Int_t iRefSel, Int_t iDut, Int_t iMRef, Int_t iBRef, Int_t iSel2,
                 Double_t dDeadTime, Bool_t bWriteHits,
                 const TString& tUnpackDir, const TString& tInputFilePath, Int_t iNDigiFiles,
                 const TString& tMCLinkFileName, Int_t iNLinkFiles, const TString& tTofGeoVersion,
                 Bool_t bHeavyIonCollisions, Bool_t bIdealClustering, Bool_t bClusteringQA,
                 Bool_t bAverageSpills, Double_t dSpillLength, Double_t dSpillBreakLength,
                 Int_t iNSerialSpills, const TString& tSetupName, Bool_t bHitBuilding,
                 Bool_t bIdealCalibration, Bool_t bSimData, const TString& tCalibDir,
                 Int_t iCalibrationMode)
{
  // ---------------------------------------------------------------------------

  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 tInputFile            = tUnpackDir + "/data/unpack.out.root";
  TString tCalibDigiInputFile   = tCalibDir + "/data/digi.in.root";
  TString tOutputFile           = "./data/digi.out.root";
  TString tHistoFile            = "./hst/calib_cluster.hst.root";
  TString tInputClustCalibFile  = "./hst/calib_cluster.cor_in.root";
  TString tOutputClustCalibFile = "./hst/calib_cluster.cor_out.root";
  TString tQAHistoFile          = "./qa/hist.qa.root";

  // ---------------------------------------------------------------------------
  if(bSimData)
  {
    TString tSetupFile  = tSrcDir + "/geometry/setup/setup_" + tSetupName + ".C";
    TString tSetupFunct = TString("setup_") + tSetupName + "()";

    gROOT->LoadMacro(tSetupFile);
    gInterpreter->ProcessLine(tSetupFunct);
  }

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

  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)
  {
    TObjString* tTofDigiTbFile     = new TObjString(tParamDir + "/tof/tof_" + tTofGeoVersion + ".digitb.par");
    tParFileList->Add(tTofDigiTbFile);
  }

  // ---------------------------------------------------------------------------
  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_digi_qa: FAIRGeom not found in geometry file!" << std::endl;
    return;
  }

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

  FairFileSource* tFileSource = new FairFileSource(tInputFile);

  if(0 == iCalibrationMode || 2 == iCalibrationMode)
  {
    tFileSource->AddFriend(tCalibDigiInputFile);
  }

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

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

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

  Int_t iSet = 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);
      }
    }
  }

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

  CbmTofEventClusterizer* tTofEventClu = new CbmTofEventClusterizer("TOF Event Clusterizer", 1, bWriteHits);
  if(bSimData)
  {
    tTofEventClu->SetPs2Ns(kFALSE);
  }
  else
  {
    tTofEventClu->SetPs2Ns(kTRUE);
  }
  tTofEventClu->SwapChannelSides(kTRUE);
  tTofEventClu->SetFileIndex(0);
  tTofEventClu->SetCalMode(iCalMode);
  tTofEventClu->SetCalSel(iCalSel);
  tTofEventClu->SetCaldXdYMax(2.);
//  tTofEventClu->SetCaldXdYMax(1000.); // TODO
  tTofEventClu->SetCalCluMulMax(10.);
//  tTofEventClu->SetCalCluMulMax(3.);
//  tTofEventClu->SetCalCluMulMax(0.);
  tTofEventClu->SetCalRpc(iCalRPC);
  tTofEventClu->SetTRefId(iRefSel);
  tTofEventClu->SetTotMax(10.);
  tTofEventClu->SetTotMin(0.01);
  tTofEventClu->SetTotPreRange(5.);
  tTofEventClu->SetTotMean(2.);
  if(0 == iCalibrationMode)
  {
    tTofEventClu->SetMaxTimeDist(0.250);
    tTofEventClu->SetSubtractPropagationTime(kTRUE);
    tTofEventClu->SetSameTimeSpaceDist(kTRUE);
  }
  else if(2 == iCalibrationMode)
  {
    tTofEventClu->SetMaxTimeDist(0.5);
    tTofEventClu->SetSubtractPropagationTime(kTRUE);
//    tTofEventClu->SetSameTimeSpaceDist(kTRUE); FIXME
    tTofEventClu->SetPositionInterpolation(kTRUE);
  }
  else
  {
    tTofEventClu->SetMaxTimeDist(0.5);
//    tTofEventClu->SetSubtractPropagationTime(kTRUE);
    tTofEventClu->SetPositionInterpolation(kTRUE); // TODO new
  }
  //tTofEventClu->SetMaxTimeDist(0.);
  tTofEventClu->SetDelTofMax(6.);
  tTofEventClu->SetBeamRefMulMax(4);
//  tTofEventClu->SetBeamRefMulMax(1);
  tTofEventClu->SetBeamAddRefMul(-1);
  if(bSimData)
  {
    tTofEventClu->SetChannelDeadtime(dDeadTime);
  }
  else
  {
    tTofEventClu->SetChannelDeadtime(1000.*dDeadTime);
  }
  tTofEventClu->SetYFitMin(1000.);
  tTofEventClu->SetMemoryTime(0.);

  tTofEventClu->SetDutId(iDutType);
  tTofEventClu->SetDutSm(iDutModule);
  tTofEventClu->SetDutRpc(iDutCounter);
  tTofEventClu->SetSelId(iMRefType);
  tTofEventClu->SetSelSm(iMRefModule);
  tTofEventClu->SetSelRpc(iMRefCounter);
  tTofEventClu->SetBeamRefId(iBRefType);
  tTofEventClu->SetBeamRefSm(iBRefModule);
  tTofEventClu->SetBeamRefDet(iBRefCounter);
  tTofEventClu->SetSel2Id(iSel2Type); 
  tTofEventClu->SetSel2Sm(iSel2Module); 
  tTofEventClu->SetSel2Rpc(iSel2Counter);

  tTofEventClu->SetCalParFileName(tInputClustCalibFile);
  tTofEventClu->SetOutHstFileName(tOutputClustCalibFile);
//FIXME
/*
//                                  10987654321098765432109876543210
  tTofEventClu->SetDeadStrips( 2, 0b10000000000000000000000000000001); // 3-0-0
  tTofEventClu->SetDeadStrips( 3, 0b00000000000000000000001000000001); // 4-0-0
  tTofEventClu->SetDeadStrips( 5, 0b00000000000000000000000000000000); // 5-1-0
  tTofEventClu->SetDeadStrips( 6, 0b00000000000000000000000000000000); // 5-2-0
  tTofEventClu->SetDeadStrips( 7, 0b00000000000000001000000000000001); // 7-0-0
  tTofEventClu->SetDeadStrips( 8, 0b00000000100000000000000000000000); // 9-0-0
*/
  tTofEventClu->SetdTRefMax(2.); // lower time-axis limit for diamonds

  switch (iCalMode)
  {
    // initial calibration
    case 0:
      tTofEventClu->SetTotMax(10.);
      //tTofEventClu->SetTotMin(1.);
      tTofEventClu->SetTRefDifMax(2000.);
      tTofEventClu->PosYMaxScal(2000.);
      tTofEventClu->SetMaxTimeDist(0.);
      break;

    // save offsets, update walks, for diamonds
    case 1:
      tTofEventClu->SetTRefDifMax(50.);
      tTofEventClu->PosYMaxScal(0.1);
      break;
    case 11:
      tTofEventClu->SetTRefDifMax(4.);
      tTofEventClu->PosYMaxScal(1.0);
      break;   
    case 21:
      tTofEventClu->SetTRefDifMax(3.);
      tTofEventClu->PosYMaxScal(1.0);
      break;
    case 31:
      tTofEventClu->SetTRefDifMax(2.);
      tTofEventClu->PosYMaxScal(1.0);
      break;
    case 41:
      tTofEventClu->SetTRefDifMax(1.);
      tTofEventClu->PosYMaxScal(0.8);
      break;   
    case 51:
      tTofEventClu->SetTRefDifMax(0.7);
      tTofEventClu->PosYMaxScal(0.7);
      break;
    case 61:
      tTofEventClu->SetTRefDifMax(0.5);
      tTofEventClu->PosYMaxScal(0.7);
      break;   
    case 71:
      tTofEventClu->SetTRefDifMax(0.4);
      tTofEventClu->PosYMaxScal(0.6);
      break;

    // time difference calibration
    case 2:
      tTofEventClu->SetTRefDifMax(300.);
      tTofEventClu->PosYMaxScal(1000.);
      break;

    // time offsets
    case 3:
      tTofEventClu->SetTRefDifMax(200.);
      tTofEventClu->PosYMaxScal(1000.);
      tTofEventClu->SetMaxTimeDist(0.);
      break;
    case 12:
    case 13:
      tTofEventClu->SetTRefDifMax(100.);
      tTofEventClu->PosYMaxScal(100.);
      break;
    case 22:
    case 23:
      tTofEventClu->SetTRefDifMax(50.);
      tTofEventClu->PosYMaxScal(50.);
      break;
    case 32:
    case 33:
      tTofEventClu->SetTRefDifMax(25.);
      tTofEventClu->PosYMaxScal(20.);
      break;
    case 42:
    case 43:
      tTofEventClu->SetTRefDifMax(12.);
      tTofEventClu->PosYMaxScal(10.);
      break;
    case 52:
    case 53:
      tTofEventClu->SetTRefDifMax(5.);
      tTofEventClu->PosYMaxScal(4.);
      break;
    case 62:
    case 63:
      tTofEventClu->SetTRefDifMax(2.5);
      tTofEventClu->PosYMaxScal(2.);
      break;
    case 72:
    case 73:
      tTofEventClu->SetTRefDifMax(2.);
      tTofEventClu->PosYMaxScal(1.0);
      break;
    case 82:
    case 83:
      tTofEventClu->SetTRefDifMax(1.);
      tTofEventClu->PosYMaxScal(0.8);
      break;
    case 92:
    case 93:
      tTofEventClu->SetTRefDifMax(0.6);
      tTofEventClu->PosYMaxScal(0.8);
      break;

    // velocity dependence
    case 4:
      tTofEventClu->SetTRefDifMax(6.);
      tTofEventClu->PosYMaxScal(1.5);
      break;
    case 14:
      tTofEventClu->SetTRefDifMax(3.);
      tTofEventClu->PosYMaxScal(1.);
      break;
    case 24:
      tTofEventClu->SetTRefDifMax(1.5);
      tTofEventClu->PosYMaxScal(0.8);
      break;
    case 44:
      tTofEventClu->SetTRefDifMax(1.);
      tTofEventClu->PosYMaxScal(0.8);
      break;
    case 54:
      tTofEventClu->SetTRefDifMax(0.7);
      tTofEventClu->PosYMaxScal(0.7);
      break;
    case 64:
      tTofEventClu->SetTRefDifMax(0.5);
      tTofEventClu->PosYMaxScal(0.7);
      break;

    default:
      std::cout << "-E- ana_digi_qa: calibration mode not implemented!" << std::endl;
      return;
  }

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

  CbmTofClusterMC* tTofClu = new CbmTofClusterMC();
  tTofClu->SetTimeSortHits(kTRUE);
//  tTofClu->BuildHitsFromPointClusters(kTRUE);
  tTofClu->BuildHitsFromRefTracks(kTRUE);
  tTofClu->SetHitCompatibilityMode(kTRUE);
  tTofClu->SetAlternativeBranchNames(kFALSE);
  tTofClu->SetFileIndex(0);

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

  CbmTofFindPointType* tTofPointType = new CbmTofFindPointType();
  tTofPointType->SetAlternativeBranchNames(kFALSE);
  tTofPointType->SetFileIndex(0);
//  tTofPointType->SetConsiderPointClusters(kTRUE);
  tTofPointType->SetConsiderRefTracks(kTRUE);

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

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

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

  CbmTofRecoMCQa* tTofQa = new CbmTofRecoMCQa();
  tTofQa->SetAverageSpills(bAverageSpills);
  tTofQa->SetLocalRadius(0.2); // in [cm]
  tTofQa->SetSpillLength(dSpillLength);
  tTofQa->SetSpillBreakLength(dSpillBreakLength);
  tTofQa->SetDataTimeOffset(0.);
  tTofQa->SetNSerialSpills(iNSerialSpills);
  tTofQa->SetTimeAxisBinWidth(0.05); // in [s]
  tTofQa->SetHistogramFileName(tQAHistoFile.Data());
  tTofQa->SetFillLocalHistos(kFALSE); // FIXME
  tTofQa->SetHitCompatibilityMode(kTRUE);
  tTofQa->SetAlternativeBranchNames(kFALSE);
//  tTofQa->SetConsiderPointClusters(kTRUE);
  tTofQa->SetConsiderRefTracks(kTRUE);
  if(bHitBuilding && !bIdealClustering)
  {
    tTofQa->SetRealClustering(kTRUE);
  }

  // ---------------------------------------------------------------------------
  if(bSimData)
  {
    tRun->AddTask(tDataManager);

    if(bHitBuilding)
    {
      if(!bIdealClustering)
      {
        tRun->AddTask(tTofEventClu);

        if(bClusteringQA)
        {
          tRun->AddTask(tTofPointType);
        }
      }
      else
      {
        if(!bIdealCalibration)
        {
          tTofEventClu->SetAlternativeBranchNames(kTRUE);
          tTofEventClu->SetWriteHitsInOut(kFALSE);
          tRun->AddTask(tTofEventClu);
        }
        else
        {
          tTofClu->SetDuplicateCounterSide(5, 1, 0, 0);
          tTofClu->SetDuplicateCounterSide(5, 2, 0, 0);
        }

        tRun->AddTask(tTofClu);
      }
    }
    else
    {
      if(bClusteringQA)
      {
        tRun->AddTask(tTofPointType);
      }
    }

    if(bClusteringQA)
    {
      if(bHitBuilding)
      {
        tRun->AddTask(tTofMatch);
      }
      tRun->AddTask(tTofQa);
    }
  }
  else
  {
    tRun->AddTask(tTofEventClu);
  }

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

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

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

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

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

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

  if((bHitBuilding && !bIdealClustering) || !bSimData)
  {
    gROOT->LoadMacro((tMacroDir + "/save_hst.C").Data());
    gROOT->LoadMacro((tMacroDir + "/fit_ybox.h").Data());
    gROOT->LoadMacro((tMacroDir + "/pl_all_CluMul.C").Data());
    gROOT->LoadMacro((tMacroDir + "/pl_all_CluRate.C").Data());
    gROOT->LoadMacro((tMacroDir + "/pl_over_cluSel.C").Data());
    gROOT->LoadMacro((tMacroDir + "/pl_over_clu.C").Data());
    gROOT->LoadMacro((tMacroDir + "/pl_all_dTSel.C").Data());

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

    switch(iSet)
    {
      case 0:
      case 300400:
      case 300700:
      case 300900:
      case 400300:
      case 400700:
      case 400900:
      case 700300:
      case 700400:
      case 700900:
      case 900300:
      case 900400:
      case 900700:
        gInterpreter->ProcessLine("pl_over_clu(5, 2, 0)");
        gInterpreter->ProcessLine("pl_over_clu(5, 1, 0)");
        gInterpreter->ProcessLine("pl_over_clu(3, 0, 0)");
        gInterpreter->ProcessLine("pl_over_clu(9, 0, 0)");
        gInterpreter->ProcessLine("pl_over_clu(7, 0, 0)");
        gInterpreter->ProcessLine("pl_over_clu(4, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(0, 5, 2, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(0, 5, 1, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(0, 3, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(0, 9, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(0, 7, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(0, 4, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(1, 5, 2, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(1, 5, 1, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(1, 3, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(1, 9, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(1, 7, 0, 0)");
        gInterpreter->ProcessLine("pl_over_cluSel(1, 4, 0, 0)");
        break;

      case 1:
        break;

      default:
        ;
    }

    gInterpreter->ProcessLine("pl_all_dTSel()");
    gInterpreter->ProcessLine("pl_all_CluMul()");
    gInterpreter->ProcessLine("pl_all_CluRate()");
  }

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

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

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