diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx deleted file mode 100644 index 3110006d11b..00000000000 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ /dev/null @@ -1,551 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// -/// @file DelphesO2LutWriter.cxx -/// @brief Porting to O2Physics of DelphesO2 code. -/// Minimal changes have been made to the original code for adaptation purposes, formatting and commented parts have been considered. -/// Relevant sources: -/// DelphesO2/src/lutWrite.cc https://github.com/AliceO2Group/DelphesO2/blob/master/src/lutWrite.cc -/// @author: Roberto Preghenella -/// @email: preghenella@bo.infn.it -/// - -#include "ALICE3/Core/DelphesO2LutWriter.h" - -#include "ALICE3/Core/DelphesO2TrackSmearer.h" -#include "ALICE3/Core/FastTracker.h" -#include "ALICE3/Core/TrackUtilities.h" - -#include -#include - -#include -#include -#include -#include -#include // IWYU pragma: keep (do not replace with TMatrixDfwd.h) -#include -#include -#include -#include -#include -#include - -#include - -#include -#include -#include -#include -#include - -// #define USE_FWD_PARAM -#ifdef USE_FWD_PARAM -#include "fwdRes.C" -#endif - -namespace o2::fastsim -{ - -void DelphesO2LutWriter::print() const -{ - LOG(info) << " --- Printing configuration of LUT writer --- "; - LOG(info) << " -> etaMaxBarrel = " << etaMaxBarrel; - LOG(info) << " -> usePara = " << usePara; - LOG(info) << " -> useDipole = " << useDipole; - LOG(info) << " -> useFlatDipole = " << useFlatDipole; - LOG(info) << " -> mAtLeastHits = " << mAtLeastHits; - LOG(info) << " -> mAtLeastCorr = " << mAtLeastCorr; - LOG(info) << " -> mAtLeastFake = " << mAtLeastFake; - LOG(info) << " -> Nch Binning: = " << mNchBinning.toString(); - LOG(info) << " -> Radius Binning: = " << mRadiusBinning.toString(); - LOG(info) << " -> Eta Binning: = " << mEtaBinning.toString(); - LOG(info) << " -> Pt Binning: = " << mPtBinning.toString(); - LOG(info) << " --- End of configuration --- "; -} - -std::string DelphesO2LutWriter::LutBinning::toString() const -{ - std::string str = ""; - str.append(log ? "log" : "lin"); - str.append(" nbins: "); - str.append(std::to_string(nbins)); - str.append(" min: "); - str.append(std::to_string(min)); - str.append(" max: "); - str.append(std::to_string(max)); - return str; -} - -bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry, - float pt, - float eta, - const float mass, - size_t itof, - size_t otof, - int q, - const float nch) -{ - lutEntry.valid = false; - - static TLorentzVector tlv; - tlv.SetPtEtaPhiM(pt, eta, 0., mass); - o2::track::TrackParCov trkIn; - o2::upgrade::convertTLorentzVectorToO2Track(q, tlv, {0., 0., 0.}, trkIn); - // tlv.Print(); - // return fmt::format("X:{:+.4e} Alp:{:+.3e} Par: {:+.4e} {:+.4e} {:+.4e} {:+.4e} {:+.4e} |Q|:{:d} {:s}\n", - // getX(), getAlpha(), getY(), getZ(), getSnp(), getTgl(), getQ2Pt(), getAbsCharge(), getPID().getName()); - // trkIn.print(); - o2::track::TrackParCov trkOut; - const int status = fat.FastTrack(trkIn, trkOut, nch); - if (status <= mAtLeastHits) { - LOGF(info, " --- fatSolve: FastTrack failed ---"); - // tlv.Print(); - return false; - } - LOGF(info, " --- fatSolve: FastTrack succeeded %d ---", status); - // trkOut.print(); - lutEntry.valid = true; - lutEntry.itof = fat.GetGoodHitProb(itof); - lutEntry.otof = fat.GetGoodHitProb(otof); - static constexpr int nCov = 15; - for (int i = 0; i < nCov; ++i) - lutEntry.covm[i] = trkOut.getCov()[i]; - - // define the efficiency - auto totfake = 0.; - lutEntry.eff = 1.; - for (size_t i = 1; i < fat.GetNLayers(); ++i) { - if (fat.IsLayerInert(i)) - continue; // skip inert layers - auto igoodhit = fat.GetGoodHitProb(i); - if (igoodhit <= 0. || i == itof || i == otof) - continue; - lutEntry.eff *= igoodhit; - auto pairfake = 0.; - for (size_t j = i + 1; j < fat.GetNLayers(); ++j) { - auto jgoodhit = fat.GetGoodHitProb(j); - if (jgoodhit <= 0. || j == itof || j == otof) - continue; - pairfake = (1. - igoodhit) * (1. - jgoodhit); - break; - } - totfake += pairfake; - } - lutEntry.eff2 = (1. - totfake); - - return true; -} - -#ifdef USE_FWD_PARAM -bool DelphesO2LutWriter::fwdSolve(float* covm, float pt, float eta, float mass) -{ - if (fwdRes(covm, pt, eta, mass) < 0) - return false; - return true; -} -#else -bool DelphesO2LutWriter::fwdSolve(float*, float, float, float) -{ - return false; -} -#endif - -bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, float mass, float Bfield) -{ - lutEntry.valid = false; - - // parametrised forward response; interpolates between FAT at eta = 1.75 and a fixed parametrisation at eta = 4; only diagonal elements - static constexpr float etaLimit = 4.0f; - if (std::fabs(eta) < etaMaxBarrel || std::fabs(eta) > etaLimit) - return false; - - if (!fatSolve(lutEntry, pt, etaMaxBarrel, mass)) - return false; - static constexpr int nCov = 15; - float covmbarrel[nCov] = {0}; - for (int i = 0; i < nCov; ++i) { - covmbarrel[i] = lutEntry.covm[i]; - } - - // parametrisation at eta = 4 - const double beta = 1. / std::sqrt(1 + mass * mass / pt / pt / std::cosh(eta) / std::cosh(eta)); - const float dcaPos = 2.5e-4 / std::sqrt(3); // 2.5 micron/sqrt(3) - const float r0 = 0.5; // layer 0 radius [cm] - const float r1 = 1.3; - const float r2 = 2.5; - const float x0layer = 0.001; // material budget (rad length) per layer - const double sigmaAlpha = 0.0136 / beta / pt * std::sqrt(x0layer * std::cosh(eta)) * (1 + 0.038 * std::log(x0layer * std::cosh(eta))); - const double dcaxyMs = sigmaAlpha * r0 * std::sqrt(1 + r1 * r1 / (r2 - r0) / (r2 - r0)); - const double dcaxy2 = dcaPos * dcaPos + dcaxyMs * dcaxyMs; - - const double dcazMs = sigmaAlpha * r0 * std::cosh(eta); - const double dcaz2 = dcaPos * dcaPos + dcazMs * dcazMs; - - const float Leta = 2.8 / std::sinh(eta) - 0.01 * r0; // m - const double relmomresPos = 10e-6 * pt / 0.3 / Bfield / Leta / Leta * std::sqrt(720. / 15.); - - const float relmomresBarrel = std::sqrt(covmbarrel[14]) * pt; - const float rOuter = 1; // m - const float relmomresPosBarrel = 10e-6 * pt / 0.3 / Bfield / rOuter / rOuter / std::sqrt(720. / 15.); - const float relmomresMSBarrel = std::sqrt(relmomresBarrel * relmomresBarrel - relmomresPosBarrel * relmomresPosBarrel); - - // interpolate MS contrib (rel resolution 0.4 at eta = 4) - const float relmomresMSEta4 = 0.4 / beta * 0.5 / Bfield; - const float relmomresMS = relmomresMSEta4 * std::pow(relmomresMSEta4 / relmomresMSBarrel, (std::fabs(eta) - 4.) / (4. - etaMaxBarrel)); - const float momresTot = pt * std::sqrt(relmomresPos * relmomresPos + relmomresMS * relmomresMS); // total absolute mom reso - - // Fill cov matrix diag - for (int i = 0; i < 15; ++i) - lutEntry.covm[i] = 0; - - lutEntry.covm[0] = covmbarrel[0]; - if (dcaxy2 > lutEntry.covm[0]) - lutEntry.covm[0] = dcaxy2; - lutEntry.covm[2] = covmbarrel[2]; - if (dcaz2 > lutEntry.covm[2]) - lutEntry.covm[2] = dcaz2; - lutEntry.covm[5] = covmbarrel[5]; // sigma^2 sin(phi) - lutEntry.covm[9] = covmbarrel[9]; // sigma^2 tanl - lutEntry.covm[14] = momresTot * momresTot / pt / pt / pt / pt; // sigma^2 1/pt - // Check that all numbers are numbers - for (int i = 0; i < 15; ++i) { - if (std::isnan(lutEntry.covm[i])) { - LOGF(info, " --- lutEntry.covm[%d] is NaN", i); - return false; - } - } - return true; -} - -void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, size_t itof, size_t otof) -{ - - if (useFlatDipole && useDipole) { - LOGF(info, "Both dipole and dipole flat flags are on, please use only one of them"); - return; - } - - // output file - std::ofstream lutFile(filename, std::ofstream::binary); - if (!lutFile.is_open()) { - LOGF(info, "Did not manage to open output file!!"); - return; - } - - // write header - lutHeader_t lutHeader; - // pid - lutHeader.pdg = pdg; - const TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdg); - if (!particle) { - LOG(fatal) << "Cannot find particle with PDG code " << pdg; - return; - } - lutHeader.mass = particle->Mass(); - const int q = std::abs(particle->Charge()) / 3; - if (q <= 0) { - LOGF(info, "Negative or null charge (%f) for pdg code %i. Fix the charge!", particle->Charge(), pdg); - return; - } - lutHeader.field = field; - auto setMap = [](map_t& map, LutBinning b) { - map.log = b.log; - map.nbins = b.nbins; - map.min = b.min; - map.max = b.max; - }; - // nch - setMap(lutHeader.nchmap, mNchBinning); - // radius - setMap(lutHeader.radmap, mRadiusBinning); - // eta - setMap(lutHeader.etamap, mEtaBinning); - // pt - setMap(lutHeader.ptmap, mPtBinning); - - lutFile.write(reinterpret_cast(&lutHeader), sizeof(lutHeader)); - - // entries - const int nnch = lutHeader.nchmap.nbins; - const int nrad = lutHeader.radmap.nbins; - const int neta = lutHeader.etamap.nbins; - const int npt = lutHeader.ptmap.nbins; - lutEntry_t lutEntry; - - // write entries - int nCalls = 0; - int successfullCalls = 0; - int failedCalls = 0; - for (int inch = 0; inch < nnch; ++inch) { - LOGF(info, " --- writing nch = %d/%d", inch, nnch); - auto nch = lutHeader.nchmap.eval(inch); - lutEntry.nch = nch; - fat.SetdNdEtaCent(nch); - for (int irad = 0; irad < nrad; ++irad) { - LOGF(info, " --- writing irad = %d/%d", irad, nrad); - for (int ieta = 0; ieta < neta; ++ieta) { - LOGF(info, " --- writing ieta = %d/%d", ieta, neta); - auto eta = lutHeader.etamap.eval(ieta); - lutEntry.eta = lutHeader.etamap.eval(ieta); - for (int ipt = 0; ipt < npt; ++ipt) { - nCalls++; - LOGF(info, " --- writing ipt = %d/%d", ipt, npt); - lutEntry.pt = lutHeader.ptmap.eval(ipt); - lutEntry.valid = true; - if (std::fabs(eta) <= etaMaxBarrel) { // full lever arm ends at etaMaxBarrel - LOGF(info, "Solving in the barrel"); - // LOGF(info, " --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); - successfullCalls++; - if (!fatSolve(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, itof, otof, q)) { - // LOGF(info, " --- fatSolve: error"); - lutEntry.valid = false; - lutEntry.eff = 0.; - lutEntry.eff2 = 0.; - for (int i = 0; i < 15; ++i) { - lutEntry.covm[i] = 0.; - } - successfullCalls--; - failedCalls++; - } - } else { - LOGF(info, "Solving outside the barrel"); - // LOGF(info, " --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); - lutEntry.eff = 1.; - lutEntry.eff2 = 1.; - bool retval = true; - successfullCalls++; - if (useFlatDipole) { // Using the parametrization at the border of the barrel - retval = fatSolve(lutEntry, lutEntry.pt, etaMaxBarrel, lutHeader.mass, itof, otof, q); - } else if (usePara) { - retval = fwdPara(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, field); - } else { - retval = fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass); - } - if (useDipole) { // Using the parametrization at the border of the barrel only for efficiency and momentum resolution - lutEntry_t lutEntryBarrel; - retval = fatSolve(lutEntryBarrel, lutEntry.pt, etaMaxBarrel, lutHeader.mass, itof, otof, q); - lutEntry.valid = lutEntryBarrel.valid; - lutEntry.covm[14] = lutEntryBarrel.covm[14]; - lutEntry.eff = lutEntryBarrel.eff; - lutEntry.eff2 = lutEntryBarrel.eff2; - } - if (!retval) { - LOGF(info, " --- fwdSolve: error"); - lutEntry.valid = false; - for (int i = 0; i < 15; ++i) { - lutEntry.covm[i] = 0.; - } - successfullCalls--; - failedCalls++; - } - } - LOGF(info, "Diagonalizing"); - diagonalise(lutEntry); - LOGF(info, "Writing"); - lutFile.write(reinterpret_cast(&lutEntry), sizeof(lutEntry_t)); - } - } - } - } - LOGF(info, " --- finished writing LUT file %s", filename); - LOGF(info, " --- successfull calls: %d/%d, failed calls: %d/%d", successfullCalls, nCalls, failedCalls, nCalls); - lutFile.close(); -} - -void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry) -{ - static constexpr int kEig = 5; - TMatrixDSym m(kEig); - for (int i = 0, k = 0; i < kEig; ++i) { - for (int j = 0; j < i + 1; ++j, ++k) { - m(i, j) = lutEntry.covm[k]; - m(j, i) = lutEntry.covm[k]; - } - } - - // m.Print(); - TMatrixDSymEigen eigen(m); - // eigenvalues vector - const TVectorD& eigenVal = eigen.GetEigenValues(); - for (int i = 0; i < kEig; ++i) - lutEntry.eigval[i] = eigenVal[i]; - // eigenvectors matrix - TMatrixD eigenVec = eigen.GetEigenVectors(); - for (int i = 0; i < kEig; ++i) - for (int j = 0; j < kEig; ++j) - lutEntry.eigvec[i][j] = eigenVec[i][j]; - // inverse eigenvectors matrix - eigenVec.Invert(); - for (int i = 0; i < kEig; ++i) - for (int j = 0; j < kEig; ++j) - lutEntry.eiginv[i][j] = eigenVec[i][j]; -} - -TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int vs, float nch, float radius, float eta, float pt) -{ - LOGF(info, " --- reading LUT file %s", filename); - // vs - static const int kNch = 0; - static const int kEta = 1; - static const int kPt = 2; - - // what - static const int kEfficiency = 0; - static const int kEfficiency2 = 1; - static const int kEfficiencyInnerTOF = 2; - static const int kEfficiencyOuterTOF = 3; - static const int kPtResolution = 4; - static const int kRPhiResolution = 5; - static const int kZResolution = 6; - - o2::delphes::DelphesO2TrackSmearer smearer; - smearer.loadTable(pdg, filename); - auto lutHeader = smearer.getLUTHeader(pdg); - lutHeader->print(); - map_t lutMap; - switch (vs) { - case kNch: - lutMap = lutHeader->nchmap; - break; - case kEta: - lutMap = lutHeader->etamap; - break; - case kPt: - lutMap = lutHeader->ptmap; - break; - } - auto nbins = lutMap.nbins; - auto g = new TGraph(); - g->SetName(Form("lut_%s_%d_vs_%d_what_%d", filename, pdg, vs, what)); - g->SetTitle(Form("LUT for %s, pdg %d, vs %d, what %d", filename, pdg, vs, what)); - switch (vs) { - case kNch: - LOGF(info, " --- vs = kNch"); - g->GetXaxis()->SetTitle("Nch"); - break; - case kEta: - LOGF(info, " --- vs = kEta"); - g->GetXaxis()->SetTitle("#eta"); - break; - case kPt: - LOGF(info, " --- vs = kPt"); - g->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - break; - default: - LOGF(info, " --- error: unknown vs %d", vs); - return nullptr; - } - switch (what) { - case kEfficiency: - LOGF(info, " --- what = kEfficiency"); - g->GetYaxis()->SetTitle("Efficiency (%)"); - break; - case kEfficiency2: - LOGF(info, " --- what = kEfficiency2"); - g->GetYaxis()->SetTitle("Efficiency2 (%)"); - break; - case kEfficiencyInnerTOF: - LOGF(info, " --- what = kEfficiencyInnerTOF"); - g->GetYaxis()->SetTitle("Inner TOF Efficiency (%)"); - break; - case kEfficiencyOuterTOF: - LOGF(info, " --- what = kEfficiencyOuterTOF"); - g->GetYaxis()->SetTitle("Outer TOF Efficiency (%)"); - break; - case kPtResolution: - LOGF(info, " --- what = kPtResolution"); - g->GetYaxis()->SetTitle("p_{T} Resolution (%)"); - break; - case kRPhiResolution: - LOGF(info, " --- what = kRPhiResolution"); - g->GetYaxis()->SetTitle("R#phi Resolution (#mum)"); - break; - case kZResolution: - LOGF(info, " --- what = kZResolution"); - g->GetYaxis()->SetTitle("Z Resolution (#mum)"); - break; - default: - LOGF(info, " --- error: unknown what %d", what); - return nullptr; - } - - bool canBeInvalid = true; - for (int i = 0; i < nbins; ++i) { - switch (vs) { - case kNch: - nch = lutMap.eval(i); - break; - case kEta: - eta = lutMap.eval(i); - break; - case kPt: - pt = lutMap.eval(i); - break; - } - float eff = 0.; - auto lutEntry = smearer.getLUTEntry(pdg, nch, radius, eta, pt, eff); - if (!lutEntry->valid || lutEntry->eff == 0.) { - if (!canBeInvalid) { - LOGF(info, " --- warning: it cannot be invalid"); - } - continue; - } - canBeInvalid = false; - - double cen = 0.; - switch (vs) { - case kNch: - cen = lutEntry->nch; - break; - case kEta: - cen = lutEntry->eta; - break; - case kPt: - cen = lutEntry->pt; - break; - } - double val = 0.; - switch (what) { - case kEfficiency: - val = lutEntry->eff * 100.; // efficiency (%) - break; - case kEfficiency2: - val = lutEntry->eff2 * 100.; // efficiency (%) - break; - case kEfficiencyInnerTOF: - val = lutEntry->itof * 100.; // efficiency (%) - break; - case kEfficiencyOuterTOF: - val = lutEntry->otof * 100.; // efficiency (%) - break; - case kPtResolution: - val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; // pt resolution (%) - break; - case kRPhiResolution: - val = std::sqrt(lutEntry->covm[0]) * 1.e4; // rphi resolution (um) - break; - case kZResolution: - val = std::sqrt(lutEntry->covm[1]) * 1.e4; // z resolution (um) - break; - default: - LOGF(info, " --- error: unknown what %d", what); - break; - } - g->AddPoint(cen, val); - } - - return g; -} -} // namespace o2::fastsim - -ClassImp(o2::fastsim::DelphesO2LutWriter); diff --git a/ALICE3/Core/DelphesO2LutWriter.h b/ALICE3/Core/DelphesO2LutWriter.h deleted file mode 100644 index e234efd3bdb..00000000000 --- a/ALICE3/Core/DelphesO2LutWriter.h +++ /dev/null @@ -1,97 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. -/// \file DelphesO2LutWriter.h -/// \brief Porting to O2Physics of DelphesO2 code. -/// Minimal changes have been made to the original code for adaptation purposes, formatting and commented parts have been considered. -/// Relevant sources: -/// DelphesO2/src/lutWrite.cc https://github.com/AliceO2Group/DelphesO2/blob/master/src/lutWrite.cc -/// \author: Roberto Preghenella -/// \email: preghenella@bo.infn.it -/// - -#ifndef ALICE3_CORE_DELPHESO2LUTWRITER_H_ -#define ALICE3_CORE_DELPHESO2LUTWRITER_H_ - -#include "ALICE3/Core/DelphesO2TrackSmearer.h" -#include "ALICE3/Core/FastTracker.h" - -#include - -#include - -#include - -#include -#include - -namespace o2::fastsim -{ -class DelphesO2LutWriter -{ - public: - DelphesO2LutWriter() = default; - virtual ~DelphesO2LutWriter() = default; - - // Setters - void setBinningNch(bool log, int nbins, float min, float max) { mNchBinning = {log, nbins, min, max}; } - void setBinningRadius(bool log, int nbins, float min, float max) { mRadiusBinning = {log, nbins, min, max}; } - void setBinningEta(bool log, int nbins, float min, float max) { mEtaBinning = {log, nbins, min, max}; } - void setBinningPt(bool log, int nbins, float min, float max) { mPtBinning = {log, nbins, min, max}; } - void setEtaMaxBarrel(float eta) { etaMaxBarrel = eta; } - void setAtLeastHits(int n) { mAtLeastHits = n; } - void setAtLeastCorr(int n) { mAtLeastCorr = n; } - void setAtLeastFake(int n) { mAtLeastFake = n; } - bool fatSolve(lutEntry_t& lutEntry, - float pt = 0.1, - float eta = 0.0, - const float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion], - size_t itof = 0, - size_t otof = 0, - int q = 1, - const float nch = 1); - - void print() const; - bool fwdSolve(float* covm, float pt = 0.1, float eta = 0.0, float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion]); - bool fwdPara(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion], float Bfield = 0.5); - void lutWrite(const char* filename = "lutCovm.dat", int pdg = 211, float field = 0.2, size_t itof = 0, size_t otof = 0); - TGraph* lutRead(const char* filename, int pdg, int what, int vs, float nch = 0., float radius = 0., float eta = 0., float pt = 0.); - - o2::fastsim::FastTracker fat; - - private: - void diagonalise(lutEntry_t& lutEntry); - float etaMaxBarrel = 1.75f; - bool usePara = true; // use fwd parameterisation - bool useDipole = false; // use dipole i.e. flat parametrization for efficiency and momentum resolution - bool useFlatDipole = false; // use dipole i.e. flat parametrization outside of the barrel - - int mAtLeastHits = 4; - int mAtLeastCorr = 4; - int mAtLeastFake = 0; - - // Binning of the LUT to make - struct LutBinning { - bool log; - int nbins; - float min; - float max; - std::string toString() const; - }; - LutBinning mNchBinning = {true, 20, 0.5f, 3.5f}; - LutBinning mRadiusBinning = {false, 1, 0.0f, 100.0f}; - LutBinning mEtaBinning = {false, 80, -4.0f, 4.0f}; - LutBinning mPtBinning = {true, 200, -2.0f, 2.0f}; - - ClassDef(DelphesO2LutWriter, 1); -}; -} // namespace o2::fastsim - -#endif // ALICE3_CORE_DELPHESO2LUTWRITER_H_ diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx deleted file mode 100644 index 7d7f6fc6077..00000000000 --- a/ALICE3/Core/DelphesO2TrackSmearer.cxx +++ /dev/null @@ -1,359 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// -/// \file DelphesO2TrackSmearer.cxx -/// \author Roberto Preghenella -/// \brief Porting to O2Physics of DelphesO2 code. -/// Minimal changes have been made to the original code for adaptation purposes, formatting and commented parts have been considered. -/// Relevant sources: -/// DelphesO2/src/lutCovm.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/lutCovm.hh -/// DelphesO2/src/TrackSmearer.cc https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.cc -/// DelphesO2/src/TrackSmearer.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.hh -/// @email: preghenella@bo.infn.it -/// - -////////////////////////////// -// DelphesO2/src/lutCovm.cc // -////////////////////////////// - -/// @author: Roberto Preghenella -/// @email: preghenella@bo.infn.it - -#include "ALICE3/Core/DelphesO2TrackSmearer.h" - -#include "ALICE3/Core/GeometryContainer.h" - -#include -#include - -#include - -#include -#include -#include -#include - -namespace o2 -{ -namespace delphes -{ - -/*****************************************************************/ - -bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) -{ - if (!filename || filename[0] == '\0') { - LOG(info) << " --- No LUT file provided for PDG " << pdg << ". Skipping load."; - return false; - } - const auto ipdg = getIndexPDG(pdg); - LOGF(info, "Will load %s lut file ..: '%s'", getParticleName(pdg), filename); - if (mLUTHeader[ipdg] && !forceReload) { - LOG(info) << " --- LUT table for PDG " << pdg << " has been already loaded with index " << ipdg << std::endl; - return false; - } - - const std::string localFilename = o2::fastsim::GeometryEntry::accessFile(filename, "./.ALICE3/LUTs/", mCcdbManager, 10); - mLUTHeader[ipdg] = new lutHeader_t; - - std::ifstream lutFile(localFilename, std::ifstream::binary); - if (!lutFile.is_open()) { - LOG(info) << " --- cannot open covariance matrix file for PDG " << pdg << ": " << localFilename << std::endl; - delete mLUTHeader[ipdg]; - mLUTHeader[ipdg] = nullptr; - return false; - } - lutFile.read(reinterpret_cast(mLUTHeader[ipdg]), sizeof(lutHeader_t)); - if (lutFile.gcount() != sizeof(lutHeader_t)) { - LOG(info) << " --- troubles reading covariance matrix header for PDG " << pdg << ": " << filename << std::endl; - LOG(info) << " --- expected/detected " << sizeof(lutHeader_t) << "/" << lutFile.gcount() << std::endl; - delete mLUTHeader[ipdg]; - mLUTHeader[ipdg] = nullptr; - return false; - } - if (mLUTHeader[ipdg]->version != LUTCOVM_VERSION) { - LOG(info) << " --- LUT header version mismatch: expected/detected = " << LUTCOVM_VERSION << "/" << mLUTHeader[ipdg]->version << std::endl; - delete mLUTHeader[ipdg]; - mLUTHeader[ipdg] = nullptr; - return false; - } - bool specialPdgCase = false; - switch (pdg) { // Handle special cases - case o2::constants::physics::kAlpha: // Special case: Allow Alpha particles to use He3 LUT - specialPdgCase = (mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3); - if (specialPdgCase) - LOG(info) - << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; - break; - default: - break; - } - if (mLUTHeader[ipdg]->pdg != pdg && !specialPdgCase) { - LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl; - delete mLUTHeader[ipdg]; - mLUTHeader[ipdg] = nullptr; - return false; - } - const int nnch = mLUTHeader[ipdg]->nchmap.nbins; - const int nrad = mLUTHeader[ipdg]->radmap.nbins; - const int neta = mLUTHeader[ipdg]->etamap.nbins; - const int npt = mLUTHeader[ipdg]->ptmap.nbins; - mLUTEntry[ipdg] = new lutEntry_t****[nnch]; - for (int inch = 0; inch < nnch; ++inch) { - mLUTEntry[ipdg][inch] = new lutEntry_t***[nrad]; - for (int irad = 0; irad < nrad; ++irad) { - mLUTEntry[ipdg][inch][irad] = new lutEntry_t**[neta]; - for (int ieta = 0; ieta < neta; ++ieta) { - mLUTEntry[ipdg][inch][irad][ieta] = new lutEntry_t*[npt]; - for (int ipt = 0; ipt < npt; ++ipt) { - mLUTEntry[ipdg][inch][irad][ieta][ipt] = new lutEntry_t; - lutFile.read(reinterpret_cast(mLUTEntry[ipdg][inch][irad][ieta][ipt]), sizeof(lutEntry_t)); - if (lutFile.gcount() != sizeof(lutEntry_t)) { - LOG(info) << " --- troubles reading covariance matrix entry for PDG " << pdg << ": " << localFilename << std::endl; - LOG(info) << " --- expected/detected " << sizeof(lutHeader_t) << "/" << lutFile.gcount() << std::endl; - return false; - } - } - } - } - } - LOG(info) << " --- read covariance matrix table for PDG " << pdg << ": " << filename << std::endl; - mLUTHeader[ipdg]->print(); - - lutFile.close(); - return true; -} - -/*****************************************************************/ - -lutEntry_t* TrackSmearer::getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff) -{ - const int ipdg = getIndexPDG(pdg); - if (!mLUTHeader[ipdg]) { - return nullptr; - } - - auto inch = mLUTHeader[ipdg]->nchmap.find(nch); - auto irad = mLUTHeader[ipdg]->radmap.find(radius); - auto ieta = mLUTHeader[ipdg]->etamap.find(eta); - auto ipt = mLUTHeader[ipdg]->ptmap.find(pt); - - // Interpolate if requested - auto fraction = mLUTHeader[ipdg]->nchmap.fracPositionWithinBin(nch); - if (mInterpolateEfficiency) { - static constexpr float kFractionThreshold = 0.5f; - if (fraction > kFractionThreshold) { - switch (mWhatEfficiency) { - case 1: - if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) { - interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff; - } else { - interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff; - } - break; - case 2: - if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) { - interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff2; - } else { - interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2; - } - break; - default: - LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency; - } - } else { - float comparisonValue = mLUTHeader[ipdg]->nchmap.log ? std::log10(nch) : nch; - switch (mWhatEfficiency) { - case 1: - if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) { - interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff; - } else { - interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff; - } - break; - case 2: - if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) { - interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff2; - } else { - interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2; - } - break; - default: - LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency; - } - } - } else { - switch (mWhatEfficiency) { - case 1: - interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff; - break; - case 2: - interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2; - break; - default: - LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency; - } - } - return mLUTEntry[ipdg][inch][irad][ieta][ipt]; -} //; - -/*****************************************************************/ - -bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float interpolatedEff) -{ - bool isReconstructed = true; - // generate efficiency - if (mUseEfficiency) { - auto eff = 0.; - switch (mWhatEfficiency) { - case 1: - eff = lutEntry->eff; - break; - case 2: - eff = lutEntry->eff2; - break; - } - if (mInterpolateEfficiency) - eff = interpolatedEff; - if (gRandom->Uniform() > eff) - isReconstructed = false; - } - - // return false already now in case not reco'ed - if (!isReconstructed && mSkipUnreconstructed) - return false; - - // transform params vector and smear - static constexpr int kParSize = 5; - double params[kParSize]; - for (int i = 0; i < kParSize; ++i) { - double val = 0.; - for (int j = 0; j < kParSize; ++j) - val += lutEntry->eigvec[j][i] * o2track.getParam(j); - params[i] = gRandom->Gaus(val, std::sqrt(lutEntry->eigval[i])); - } - // transform back params vector - for (int i = 0; i < kParSize; ++i) { - double val = 0.; - for (int j = 0; j < kParSize; ++j) - val += lutEntry->eiginv[j][i] * params[j]; - o2track.setParam(val, i); - } - // should make a sanity check that par[2] sin(phi) is in [-1, 1] - if (std::fabs(o2track.getParam(2)) > 1.) { - LOG(info) << " --- smearTrack failed sin(phi) sanity check: " << o2track.getParam(2) << std::endl; - } - // set covariance matrix - static constexpr int kCovMatSize = 15; - for (int i = 0; i < kCovMatSize; ++i) - o2track.setCov(lutEntry->covm[i], i); - return isReconstructed; -} - -/*****************************************************************/ - -bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch) -{ - auto pt = o2track.getPt(); - switch (pdg) { - case o2::constants::physics::kHelium3: - case -o2::constants::physics::kHelium3: - pt *= 2.f; - break; - } - auto eta = o2track.getEta(); - float interpolatedEff = 0.0f; - lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, interpolatedEff); - if (!lutEntry || !lutEntry->valid) - return false; - return smearTrack(o2track, lutEntry, interpolatedEff); -} - -/*****************************************************************/ -// relative uncertainty on pt -double TrackSmearer::getPtRes(const int pdg, const float nch, const float eta, const float pt) -{ - float dummy = 0.0f; - lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); - auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt; - return val; -} - -/*****************************************************************/ -// relative uncertainty on eta -double TrackSmearer::getEtaRes(const int pdg, const float nch, const float eta, const float pt) -{ - float dummy = 0.0f; - lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); - auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2 - auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty - etaRes /= lutEntry->eta; // relative uncertainty - return etaRes; -} -/*****************************************************************/ -// absolute uncertainty on pt -double TrackSmearer::getAbsPtRes(const int pdg, const float nch, const float eta, const float pt) -{ - float dummy = 0.0f; - lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); - auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * lutEntry->pt; - return val; -} - -/*****************************************************************/ -// absolute uncertainty on eta -double TrackSmearer::getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt) -{ - float dummy = 0.0f; - lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); - auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2 - auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty - return etaRes; -} -/*****************************************************************/ -// efficiency -double TrackSmearer::getEfficiency(const int pdg, const float nch, const float eta, const float pt) -{ - float efficiency = 0.0f; - getLUTEntry(pdg, nch, 0., eta, pt, efficiency); - return efficiency; -} -/*****************************************************************/ -// Only in DelphesO2 -// bool TrackSmearer::smearTrack(Track& track, bool atDCA) -// { - -// O2Track o2track; -// TrackUtils::convertTrackToO2Track(track, o2track, atDCA); -// int pdg = track.PID; -// float nch = mdNdEta; // use locally stored dNch/deta for the time being -// if (!smearTrack(o2track, pdg, nch)) -// return false; -// TrackUtils::convertO2TrackToTrack(o2track, track, atDCA); -// return true; - -// #if 0 -// lutEntry_t* lutEntry = getLUTEntry(track.PID, 0., 0., track.Eta, track.PT); -// if (!lutEntry) -// return; - -// O2Track o2track; -// TrackUtils::convertTrackToO2Track(track, o2track, atDCA); -// smearTrack(o2track, lutEntry); -// TrackUtils::convertO2TrackToTrack(o2track, track, atDCA); -// #endif -// } - -/*****************************************************************/ - -} // namespace delphes -} // namespace o2 diff --git a/ALICE3/Core/DelphesO2TrackSmearer.h b/ALICE3/Core/DelphesO2TrackSmearer.h deleted file mode 100644 index 6d941f336ca..00000000000 --- a/ALICE3/Core/DelphesO2TrackSmearer.h +++ /dev/null @@ -1,258 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// -/// @file DelphesO2TrackSmearer.h -/// @brief Porting to O2Physics of DelphesO2 code. -/// Minimal changes have been made to the original code for adaptation purposes, formatting and commented parts have been considered. -/// Relevant sources: -/// DelphesO2/src/lutCovm.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/lutCovm.hh -/// DelphesO2/src/TrackSmearer.cc https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.cc -/// DelphesO2/src/TrackSmearer.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.hh -/// @author: Roberto Preghenella -/// @email: preghenella@bo.infn.it -/// - -#ifndef ALICE3_CORE_DELPHESO2TRACKSMEARER_H_ -#define ALICE3_CORE_DELPHESO2TRACKSMEARER_H_ - -#include -#include - -#include -#include -#include - -/////////////////////////////// -/// DelphesO2/src/lutCovm.hh // -/////////////////////////////// - -#define LUTCOVM_VERSION 20210801 - -struct map_t { - int nbins = 1; - float min = 0.; - float max = 1.e6; - bool log = false; - float eval(int bin) - { - float width = (max - min) / nbins; - float val = min + (bin + 0.5) * width; - if (log) - return pow(10., val); - return val; - } - // function needed to interpolate some dimensions - float fracPositionWithinBin(float val) - { - float width = (max - min) / nbins; - int bin; - float returnVal = 0.5f; - if (log) { - bin = static_cast((log10(val) - min) / width); - returnVal = ((log10(val) - min) / width) - bin; - } else { - bin = static_cast((val - min) / width); - returnVal = val / width - bin; - } - return returnVal; - } - - int find(float val) - { - float width = (max - min) / nbins; - int bin; - if (log) - bin = static_cast((log10(val) - min) / width); // Changed due to MegaLinter error. - // bin = (int)((log10(val) - min) / width); // Original line. - else - bin = static_cast((val - min) / width); // Changed due to MegaLinter error. - // bin = (int)((val - min) / width); // Original line. - if (bin < 0) - return 0; - if (bin > nbins - 1) - return nbins - 1; - return bin; - } //; - void print() { printf("nbins = %d, min = %f, max = %f, log = %s \n", nbins, min, max, log ? "on" : "off"); } //; -}; - -struct lutHeader_t { - int version = LUTCOVM_VERSION; - int pdg = 0; - float mass = 0.; - float field = 0.; - map_t nchmap; - map_t radmap; - map_t etamap; - map_t ptmap; - bool check_version() - { - return (version == LUTCOVM_VERSION); - } //; - void print() - { - printf(" version: %d \n", version); - printf(" pdg: %d \n", pdg); - printf(" field: %f \n", field); - printf(" nchmap: "); - nchmap.print(); - printf(" radmap: "); - radmap.print(); - printf(" etamap: "); - etamap.print(); - printf(" ptmap: "); - ptmap.print(); - } //; -}; - -struct lutEntry_t { - float nch = 0.; - float eta = 0.; - float pt = 0.; - bool valid = false; - float eff = 0.; - float eff2 = 0.; - float itof = 0.; - float otof = 0.; - float covm[15] = {0.}; - float eigval[5] = {0.}; - float eigvec[5][5] = {{0.}}; - float eiginv[5][5] = {{0.}}; - void print() - { - printf(" --- lutEntry: pt = %f, eta = %f (%s)\n", pt, eta, valid ? "valid" : "not valid"); - printf(" efficiency: %f\n", eff); - printf(" covMatix: "); - int k = 0; - for (int i = 0; i < 5; ++i) { - for (int j = 0; j < i + 1; ++j) - printf("% e ", covm[k++]); - printf("\n "); - } - printf("\n"); - } -}; - -//////////////////////////////////// -/// DelphesO2/src/TrackSmearer.hh // -//////////////////////////////////// - -using O2Track = o2::track::TrackParCov; - -namespace o2 -{ -namespace delphes -{ - -class TrackSmearer -{ - - public: - TrackSmearer() = default; - ~TrackSmearer() = default; - - /** LUT methods **/ - bool loadTable(int pdg, const char* filename, bool forceReload = false); - bool hasTable(int pdg) { return (mLUTHeader[getIndexPDG(pdg)] != nullptr); } //; - void useEfficiency(bool val) { mUseEfficiency = val; } //; - void interpolateEfficiency(bool val) { mInterpolateEfficiency = val; } //; - void skipUnreconstructed(bool val) { mSkipUnreconstructed = val; } //; - void setWhatEfficiency(int val) { mWhatEfficiency = val; } //; - lutHeader_t* getLUTHeader(int pdg) { return mLUTHeader[getIndexPDG(pdg)]; } //; - lutEntry_t* getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff); - - bool smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float interpolatedEff); - bool smearTrack(O2Track& o2track, int pdg, float nch); - // bool smearTrack(Track& track, bool atDCA = true); // Only in DelphesO2 - double getPtRes(const int pdg, const float nch, const float eta, const float pt); - double getEtaRes(const int pdg, const float nch, const float eta, const float pt); - double getAbsPtRes(const int pdg, const float nch, const float eta, const float pt); - double getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt); - double getEfficiency(const int pdg, const float nch, const float eta, const float pt); - - int getIndexPDG(const int pdg) - { - switch (abs(pdg)) { - case 11: - return 0; // Electron - case 13: - return 1; // Muon - case 211: - return 2; // Pion - case 321: - return 3; // Kaon - case 2212: - return 4; // Proton - case 1000010020: - return 5; // Deuteron - case 1000010030: - return 6; // Triton - case 1000020030: - return 7; // Helium3 - case 1000020040: - return 8; // Alphas - default: - return 2; // Default: pion - } - } - - const char* getParticleName(int pdg) - { - switch (abs(pdg)) { - case 11: - return "electron"; - case 13: - return "muon"; - case 211: - return "pion"; - case 321: - return "kaon"; - case 2212: - return "proton"; - case 1000010020: - return "deuteron"; - case 1000010030: - return "triton"; - case 1000020030: - return "helium3"; - case 1000020040: - return "alpha"; - default: - return "pion"; // Default: pion - } - } - void setdNdEta(float val) { mdNdEta = val; } //; - void setCcdbManager(o2::ccdb::BasicCCDBManager* mgr) { mCcdbManager = mgr; } //; - - protected: - static constexpr unsigned int nLUTs = 9; // Number of LUT available - lutHeader_t* mLUTHeader[nLUTs] = {nullptr}; - lutEntry_t***** mLUTEntry[nLUTs] = {nullptr}; - bool mUseEfficiency = true; - bool mInterpolateEfficiency = false; - bool mSkipUnreconstructed = true; // don't smear tracks that are not reco'ed - int mWhatEfficiency = 1; - float mdNdEta = 1600.; - - private: - o2::ccdb::BasicCCDBManager* mCcdbManager = nullptr; -}; - -} // namespace delphes -} // namespace o2 - -namespace o2::delphes -{ -using DelphesO2TrackSmearer = TrackSmearer; -} - -#endif // ALICE3_CORE_DELPHESO2TRACKSMEARER_H_