LmvmEventMix.cxx

Go to the documentation of this file.

/* Copyright (C) 2023-2025 UGiessen, Giessen
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Cornelius Feier-Riesen [committer] */
 
#include "LmvmEventMix.h"
 
#include "CbmHistManager.h"
#include "CbmUtils.h"
#include "LmvmDef.h"
#include "LmvmKinePar.h"
#include "LmvmSimParam.h"
#include "LmvmUtils.h"
#include "TF1.h"
#include "TFile.h"
#include "TH1.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "TSystemDirectory.h"
 
#include <iomanip>
#include <iostream>
#include <string>
 
using namespace std;
using namespace Cbm;
 
/*********************************************************************************************************************
* This class is for providing invariant mass histograms from same and mixed events for like-sign and unlike-sign     *
* pairs to be used for the calculation of the combinatorial background in LmvmDrawAll. For this, in the analysis     *
* (LmvmTask) a root file is created that contains one vector<LmvmCand> (stored in file "lmvmcand.taskId.root") filled*
* with all LmvmCand candidates obtained in one job (max. 1000 events). From LmvmCand all necessary information can be*
* obtained for the required histograms. The parameter 'fNofFilesToMix' determines how many files are used for event  *
* mixing. That means, if fNofFilesToMix=10 and each file contains 1,000 events, the mixing depth is 10,000 events.   *
* In the same-event loop a 'vector<LmvmCand> candsSame' is filled with the vector elements extracted from the file   *
* with number 'taskId' and handed over to the function FillMinvHistos() with option "same". In the mixed-event loop  *
* a 'vector<LmvmCand> candsMix' is filled with the vector elements of exactly fNofFilesToMix files, starting at file *
* #taskId until file #taskId+fNofFilesToMix-1. If one of these files can not be opened or does not exist, the range  *
* is extended by 1 to make sure to mix exactly fNofFilesToMix files. After that, candsMix is handed over to          *
* FillMinvHistos() with the option "mix".                                                                            *
* In FillMinvHistos() the invariant mass for each possible pair combination is calculated and filled into an         *
* according histogram (like-sign or unlike-sign pair, same or mixed events).                                         *
**********************************************************************************************************************/
 
void LmvmEventMix::DoEventMix(const string& dataDir, const string& outputDir, const string& outFile,
                              const string& taskIdStr, int nofFilesToMix)
{
  CreateHistos();
 
  fDataDir       = dataDir;
  fOutputDir     = outputDir;
  fFileName      = outFile;
  fNofFilesToMix = nofFilesToMix;
  int taskId     = std::stoi(taskIdStr);
  cout << endl;
  LOG(info) << "LmvmEventMix::DoEventMix: Energy string = " << fEnergy
            << ". Mind this for calculation of weighting values!" << endl;
  fHMean.FillH1("hmx_nofJobs", 0.5);
 
  TStopwatch timer;
  timer.Start();
 
  // Choose signal (only one of five signals per run)
  int iSig = taskId % fHMean.fNofSignals;
  fSignal  = fHMean.fSignalNames[iSig];
  LOG(info) << "taskId = " << taskId << ", iSig = " << iSig << ", fSignal = " << fSignal;
 
  // Correct assignment of folder names to number 'iSig' is essential!
  string folder = (iSig == 0)   ? "inmed/"
                  : (iSig == 1) ? "qgp/"
                  : (iSig == 2) ? "omegaepem/"
                  : (iSig == 3) ? "phi/"
                  : (iSig == 4) ? "omegadalitz/"
                                : "";
  if (folder == "") LOG(error) << "LmvmEventMix::DoEventMix(): No folder name specified. Check!";
 
  // Get lmvmcand file with highest number; needed in mixed-event loop
  int iFileMax = 0;
  for (int iF = 1; iF <= 1000; iF++) {
    string iFString = Cbm::NumberToString(iF, 0);
    string filename = fDataDir + folder + "lmvmcand." + iFString + ".root";
    TFile* file(TFile::Open(filename.c_str()));
    if (file) iFileMax = iF;
  }
  LOG(info) << "iFileMax = " << iFileMax;
 
  // Filling candidate vectors (same and mixed events)
 
  // Loop for same events (use each file exactly once)
  vector<LmvmCand> candsSame;
  fNofEvents = 0;
  for (const string name : {"inmed", "omegaepem", "omegadalitz", "phi", "qgp"}) {
    candsSame.clear();
    string filename = fDataDir + name + "/lmvmcand." + taskIdStr + ".root";
    TFile* file(TFile::Open(filename.c_str()));
    if (!file) {
      LOG(info) << "file " << filename << " cannot be opened.";
      continue;
    }
    vector<LmvmCand>* cands(file->Get<vector<LmvmCand>>("fCandsTotal"));
    if (cands != nullptr) {
      LOG(info) << "File taken for same events: " << filename;
      candsSame.insert(candsSame.end(), cands->begin(), cands->end());
      string filenameAna = fDataDir + name + "/analysis." + taskIdStr + ".root";
      TFile* anafile(TFile::Open(filenameAna.c_str()));
      if (anafile) {
        TH1D* h(anafile->Get<TH1D>("hEventNumber"));
        fNofEvents += h->GetEntries();
      }
      else
        LOG(warning) << "LmvmEventMix (loop 'same'): anafile " << filenameAna
                     << " could not be opened or does not exist.";
    }
    else
      LOG(warning) << "LmvmEventMix::DoEventMix(): Loop 'same' with filename = " << filename
                   << ": candsSame = nullptr!";
    FillMinvHistos(candsSame, "same");
  }
 
  // Fill histograms for number of events and number of files
  // NOTE: why do the first lines not return the right number after 'hadd'??
  //fHMean.FillH1("hmx_EventNumber", 0.5, static_cast<double>(fNofEvents));
  //fHMean.FillH1("hmx_NofFilesToMix", 0.5, static_cast<double>(fNofFilesToMix));
  //fHMean.FillH1("hmx_EventNumber", 0.5, fNofEvents);
  //fHMean.FillH1("hmx_NofFilesToMix", 0.5, fNofFilesToMix);
  for (int i = 1; i <= fNofEvents; i++) {
    fHMean.FillH1("hmx_EventNumber", 0.5);
  }
  for (int i = 1; i <= fNofFilesToMix; i++) {
    fHMean.FillH1("hmx_NofFilesToMix", 0.5);
  }
 
  // Loop for mixed events (use exactly fNofFilesToMix files)
  vector<LmvmCand> candsMix;
  candsMix.clear();
  int iter        = -1;
  int nMixedFiles = 0;
  while (nMixedFiles < fNofFilesToMix) {
    iter++;  // to stop loop if e.g. no file is to open and nMixedFiles would not increase
    if (iter > 1000) {
      LOG(error) << "LmvmEventMix::DoEventMix(): iter = " << iter << ". Is too large. Check!";
      return;
    }
    int iFile = (taskId + iter == iFileMax) ? taskId + iter : (taskId + iter) % iFileMax;  // because no file #0 exists
    string iFileString = Cbm::NumberToString(iFile, 0);
    string filename    = fDataDir + folder + "lmvmcand." + iFileString + ".root";
    TFile* file(TFile::Open(filename.c_str()));
    if (!file) continue;
    vector<LmvmCand>* cands(file->Get<vector<LmvmCand>>("fCandsTotal"));
    if (cands != nullptr) {
      LOG(info) << "File taken for mixing: " << filename;
      candsMix.insert(candsMix.end(), cands->begin(), cands->end());
      nMixedFiles++;
    }
    else
      LOG(warning) << "LmvmEventMix::DoEventMix(): Loop 'mix': filename = " << filename << ", candsMix = nullptr!";
  }
  if (nMixedFiles < fNofFilesToMix) LOG(warn) << "nMixedFiles < fNofFilesToMix. Check cause!";
  FillMinvHistos(candsMix, "mix");
 
  // Finish
  double bW = fHMean.H1("hmx_MinvPM_sameEv_elid")->GetBinWidth(1);
  fHMean.fHM.ScaleByPattern("hmx_Minv.*", 1. / bW);
 
  int nCandsMix  = candsMix.size();
  int nCandsSame = candsSame.size();
  LOG(info) << "fNofEvents = " << fNofEvents;
  LOG(info) << "nCandsSame = " << nCandsSame << ", nCandsMix = " << nCandsMix;
 
  timer.Stop();
  LOG(info) << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s.";
  fHMean.FillH1("hmx_time", timer.RealTime());
 
  CheckMemory("End");
  SaveHist();
}
 
void LmvmEventMix::FillMinvHistos(const vector<LmvmCand>& cands, const string& mode)
{
  if (mode != "same" && mode != "mix")
    LOG(error) << "LmvmEventMix::FillMinvHistos: event mode must be 'same' or 'mix'!";
  int nCands = cands.size();
  for (int iC1 = 0; iC1 < nCands - 1; iC1++) {
    const auto& cand1 = cands[iC1];
    if (cand1.fCharge == 0) LOG(error) << "LmvmEventMix::FillMinvHistos: Charge of cand1 is 0. Check!";
    for (int iC2 = iC1 + 1; iC2 < nCands; iC2++) {
      const auto& cand2 = cands[iC2];
      if (cand2.fCharge == 0) LOG(error) << "LmvmEventMix::FillMinvHistos: Charge of cand2 is 0. Check!";
      if (cand1.fTaskId == "" || cand2.fTaskId == "")
        LOG(error) << "LmvmEventMix::FillMinvHistos: String 'fTaskId' of at least one candidate is empty. Check!";
      bool isSameEvent =
        (cand1.fEventNumber == cand2.fEventNumber && cand1.fTaskId == cand2.fTaskId && cand1.fName == cand2.fName);
      if ((mode == "same" && !isSameEvent) || (mode == "mix" && isSameEvent)) continue;
      string evmode = "";
      if (mode == "same" && isSameEvent) evmode = "_sameEv";
      if (mode == "mix" && !isSameEvent) evmode = "_mixedEv";
      if (evmode == "") {
        LOG(error) << "LmvmEventMix::FillMinvHistos: string 'evmode' is not defined. Check!";
        continue;
      }  // TODO: should not occur; delete this line
      string comb = (cand1.fCharge * cand2.fCharge < 0)        ? "PM"
                    : (cand1.fCharge < 0 && cand2.fCharge < 0) ? "MM"
                    : (cand1.fCharge > 0 && cand2.fCharge > 0) ? "PP"
                                                               : "";
      if (comb == "")
        LOG(error)
          << "LmvmEventMix::FillMinvHistos: 'comb' is empty. Check!";  // TODO: should not occur; delete this line and shorten line before (without last condition)
 
      double w         = (isSameEvent) ? LmvmUtils::GetWeightPair(cand1, cand2) : cand1.fWeight * cand2.fWeight;
      LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2);
      for (auto step : fHMean.fAnaSteps) {
        if (step < ELmvmAnaStep::ElId) continue;
        if (!(cand1.IsCutTill(step) && cand2.IsCutTill(step))) continue;
        fHMean.FillH1("hmx_Minv" + comb + evmode, step, pRec.fMinv, w);
        if (std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11)
          fHMean.FillH1("hmx_Minv" + comb + "_trueEl" + evmode, step, pRec.fMinv, w);
        if (!cand1.IsMcSignal() && !cand2.IsMcSignal())
          fHMean.FillH1("hmx_Minv" + comb + "_urqmdAll" + evmode, step, pRec.fMinv, w);
        if (!cand1.IsMcSignal() && !cand2.IsMcSignal() && std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11)
          fHMean.FillH1("hmx_Minv" + comb + "_urqmdEl" + evmode, step, pRec.fMinv, w);
      }
    }
  }
}
 
void LmvmEventMix::CreateHistos()
{
  for (const string comb : {"PM", "PP", "MM"}) {
    string hName = "hmx_Minv" + comb;
    for (const string cat : {"", "_trueEl", "_urqmdAll", "_urqmdEl"}) {
      string hNameCat = hName + cat;
      fHMean.CreateH1(hNameCat, {"sameEv", "mixedEv"}, fHMean.fAnaStepNames, "M_{ee} [GeV/c^{2}]",
                      "dN/dM_{ee}/N_{Ev} [GeV/c^{2}]^{-1}", 250, 0., 2.5, true);
    }
  }
  fHMean.CreateH1("hmx_EventNumber", "", "", 1, 0, 1.);    // TODO: needed?
  fHMean.CreateH1("hmx_nofJobs", "", "", 1, 0, 1.);        // for correct scaling in LmvmDrawAll
  fHMean.CreateH1("hmx_NofFilesToMix", "", "", 1, 0, 1.);  // To get mixing depth in LmvmDrawAll
  fHMean.CreateH1("hmx_time", "Time [s]", "a.u.", 100, 0., 200.);
}
 
void LmvmEventMix::CheckMemory(const string& text)
{
  static ProcInfo_t info;
  const float toMB = 1.f / 1024.f;
  gSystem->GetProcInfo(&info);
  cout << Form((text + ": Memory res|vir: %g|%g Mbytes\n").c_str(), info.fMemResident * toMB, info.fMemVirtual * toMB);
  cout << endl;
}
 
void LmvmEventMix::SaveHist()
{
  if (fOutputDir != "") {
    gSystem->mkdir(fOutputDir.c_str(), true);
    TFile* f = TFile::Open(string(fFileName).c_str(), "RECREATE");
    fHMean.WriteToFile();
    f->Close();
  }
  fHMean.fHM.SaveCanvasToImage(fOutputDir, "png");
}
 
ClassImp(LmvmEventMix);