diff --git a/ALICE3/DataModel/OTFMID.h b/ALICE3/DataModel/OTFMID.h new file mode 100644 index 00000000000..363547294e7 --- /dev/null +++ b/ALICE3/DataModel/OTFMID.h @@ -0,0 +1,68 @@ +// 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 OTFMID.h +/// \author David Dobrigkeit Chinellato +/// \author Antonio Ortiz, UNAM +/// \since 11/05/2023 +/// \brief Set of tables for the ALICE3 OTFMID information +/// TODO: remove the columns that are not usefulf for me + +#ifndef ALICE3_DATAMODEL_MID_H_ +#define ALICE3_DATAMODEL_MID_H_ + +#include + +#include +#include + +namespace o2::aod +{ +namespace upgrade_mid +{ +DECLARE_SOA_COLUMN(MIDTrackTime, midTrackTime, float); //! Track time generated at the MID +DECLARE_SOA_COLUMN(MIDTrackLength, midTrackLength, float); //! track length for calculation of MID (generated) + +DECLARE_SOA_COLUMN(MIDEventTime, midEventTime, float); //! Event time reconstructed with the MID +DECLARE_SOA_COLUMN(MIDEventTimeErr, midEventTimeErr, float); //! Uncertainty on the event time reconstructed with the MID +DECLARE_SOA_COLUMN(MIDTrackPosReco, midTrackPosReco, std::vector); //! Measured position in MID +DECLARE_SOA_COLUMN(MIDTrackPos, midTrackPos, std::vector); //! Expected position in MID +DECLARE_SOA_COLUMN(MIDCandidate, midCandidate, bool); //! Muon candidate +DECLARE_SOA_COLUMN(MIDTrackTimeReco, midTrackTimeReco, float); //! Track time measured at MID +DECLARE_SOA_COLUMN(MIDTrackLengthReco, midTrackLengthReco, float); //! track length for calculation of MID (reconstructed) + +} // namespace upgrade_mid + +DECLARE_SOA_TABLE(UpgradeMidMCs, "AOD", "UPGRADEMIDMC", + upgrade_mid::MIDTrackTime, + upgrade_mid::MIDTrackLength); + +DECLARE_SOA_TABLE(UpgradeMids, "AOD", "UPGRADEMID", + upgrade_mid::MIDEventTime, + upgrade_mid::MIDEventTimeErr, + upgrade_mid::MIDTrackPosReco, + upgrade_mid::MIDTrackTimeReco, + upgrade_mid::MIDTrackLengthReco, + upgrade_mid::MIDCandidate); + +DECLARE_SOA_TABLE(UpgradeMidExpecteds, "AOD", "UPGRADEMIDEXPT", + upgrade_mid::MIDTrackPos, + upgrade_mid::MIDTrackTime, + upgrade_mid::MIDTrackLength); + +using UpgradeMidMC = UpgradeMidMCs::iterator; +using UpgradeMid = UpgradeMids::iterator; +using UpgradeMidExpected = UpgradeMidExpecteds::iterator; + +} // namespace o2::aod + +#endif // ALICE3_DATAMODEL_OTFMID_H_ diff --git a/ALICE3/DataModel/ReducedTablesAlice3.h b/ALICE3/DataModel/ReducedTablesAlice3.h index 363ba50842f..7ff870f9308 100644 --- a/ALICE3/DataModel/ReducedTablesAlice3.h +++ b/ALICE3/DataModel/ReducedTablesAlice3.h @@ -20,6 +20,7 @@ #include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "ALICE3/DataModel/OTFMID.h" #include "ALICE3/DataModel/OTFPIDTrk.h" #include "ALICE3/DataModel/OTFRICH.h" #include "ALICE3/DataModel/OTFTOF.h" @@ -275,6 +276,14 @@ DECLARE_SOA_TABLE(ReducedA3PIDTOF, "AOD", "REA3PIDTOF", upgrade_tof::NSigmaHelium3OuterTOF, upgrade_tof::NSigmaAlphaOuterTOF>); +DECLARE_SOA_TABLE(ReducedA3PIDMID, "AOD", "REA3PIDMID", + upgrade_mid::MIDEventTime, + upgrade_mid::MIDEventTimeErr, + upgrade_mid::MIDTrackPosReco, + upgrade_mid::MIDTrackTimeReco, + upgrade_mid::MIDTrackLengthReco, + upgrade_mid::MIDCandidate); + DECLARE_SOA_TABLE(ReducedA3PIDRich, "AOD", "REA3PIDRICH", upgrade_rich::NSigmaElectronRich, upgrade_rich::NSigmaMuonRich, diff --git a/ALICE3/TableProducer/OTF/CMakeLists.txt b/ALICE3/TableProducer/OTF/CMakeLists.txt index 6fe3f41a02f..d92809cae76 100644 --- a/ALICE3/TableProducer/OTF/CMakeLists.txt +++ b/ALICE3/TableProducer/OTF/CMakeLists.txt @@ -24,6 +24,11 @@ o2physics_add_dpl_workflow(onthefly-tofpid PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core O2Physics::FastTracker COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(onthefly-midpid + SOURCES onTheFlyMidPid.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core O2Physics::FastTracker + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(onthefly-richpid SOURCES onTheFlyRichPid.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core O2Physics::FastTracker diff --git a/ALICE3/TableProducer/OTF/onTheFlyMidPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyMidPid.cxx new file mode 100644 index 00000000000..bda0b32d5f2 --- /dev/null +++ b/ALICE3/TableProducer/OTF/onTheFlyMidPid.cxx @@ -0,0 +1,1124 @@ +// 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 onTheFlyMidPid.cxx +/// +/// \brief This task goes straight from a combination of track table and mcParticles +/// and a custom MID configuration to a table of muon id for the particles +/// being analysed. +/// +/// \author David Dobrigkeit Chinellato, UNICAMP +/// \author Antonio Ortiz, UNAM +/// \since May 22, 2024 +/// + +#include "GeometryContainer.h" + +#include "ALICE3/Core/FlatTrackSmearer.h" +#include "ALICE3/Core/TrackUtilities.h" +#include "ALICE3/DataModel/OTFCollision.h" +#include "ALICE3/DataModel/OTFMID.h" +#include "Common/Core/trackUtilities.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace o2; +using namespace o2::framework; + +struct OnTheFlyMidPid { + Produces upgradeMidMC; + Produces upgradeMid; + Produces upgradeMidExpected; + + // necessary for particle charges + Service pdg; + // Necessary for LUTs + Service ccdb; + + // these are the settings governing the MID layer to be used + // in the evolution of this effort + struct : ConfigurableGroup { + Configurable considerEventTime{"considerEventTime", true, "flag to consider event time"}; + // Parametrization of muon/pion absortion in iron + Configurable magRadiusC{"magRadiusC", 160., "Magnet coil radius (cm)"}; + Configurable absRadius1{"absRadius1", 220., "Absorber inner radius (cm)"}; + Configurable absRadius2{"absRadius2", 290., "Absorber outer radius (cm)"}; + Configurable absDensity{"absDensity", 7.874, "iron density (g/cm3)"}; + Configurable absLambdaPi{"absLambdaPi", 20.41, "illustrative effective interaction length (cm)"}; + Configurable absDeDx{"absDeDxMIP", 1.7, "dE/(rho dx) MIP MeV/g/cm2"}; + Configurable absFakeAt3{"absFakeAt3", 0.032, "fake-muon efficiency at 3 GeV/c from GEANT4+BDT abs(eta)<1.24, 70 cm"}; + Configurable absFakeAt10{"absFakeAt10", 0.08, "fake-muon efficiency at 10 GeV/c from GEANT4+BDT abs(eta)<1.24, 70 cm"}; + Configurable absSigmaEta{"absSigmaEta", 10.65, "user-configurable phenomenological eta width"}; + Configurable absPionFractionAt3{"absPionFractionAt3", 0.5, "parameter from GEANT4+BDT"}; // assumed, NOT measured + Configurable absEtaMax{"absEtaMax", 1.24, "max eta covered by detector"}; + Configurable midRadius{"midRadius", 302.5, "MID radius (cm)"}; + Configurable midLength{"midLength", 1000.0, "MID length (cm)"}; + Configurable> midCellDimension{"midCellDimension", {5.1, 5.1}, "MID cell dimension in Z and RPhi (cm)"}; + Configurable midFractionOfInactiveArea{"midFractionOfInactiveArea", (0.1 / 5.1), "MID fraction of inactive area"}; + Configurable midTimeReso{"midTimeReso", 1500, "MID time error (ps)"}; + Configurable nStepsLIntegrator{"nStepsLIntegrator", 200, "number of steps in length integrator"}; + Configurable multiplicityEtaRange{"multiplicityEtaRange", 0.800000012, "eta range to compute the multiplicity"}; + Configurable flagIncludeTrackTimeRes{"flagIncludeTrackTimeRes", true, "flag to include or exclude track time resolution"}; + Configurable flagMIDLoadDelphesLUTs{"flagMIDLoadDelphesLUTs", false, "flag to load Delphes LUTs for tracking correction (use recoTrack parameters if false)"}; + } simConfig; + + struct : ConfigurableGroup { + Configurable doQAplots{"doQAplots", true, "do basic velocity plot qa"}; + Configurable doSeparationVsPt{"doSeparationVsPt", true, "Produce plots vs pt or p"}; + Configurable nBinsBeta{"nBinsBeta", 2200, "number of bins in beta"}; + Configurable nBinsP{"nBinsP", 80, "number of bins in momentum/pT depending on doSeparationVsPt"}; + Configurable nBinsTrackLength{"nBinsTrackLength", 500, "number of bins in track length"}; + Configurable nBinsNsigmaWrongSpecies{"nBinsNsigmaWrongSpecies", 200, "number of bins in Nsigma plot (wrong species)"}; + Configurable minNsigmaRange{"minNsigmaRange", -10, "lower limit for the nsigma axis"}; + Configurable maxNsigmaRange{"maxNsigmaRange", +10, "upper limit for the nsigma axis"}; + Configurable nBinsDeltaCorrectSpecies{"nBinsDeltaCorrectSpecies", 200, "number of bins in Delta plot (correct speies)"}; + Configurable nBinsDeltaWrongSpecies{"nBinsDeltaWrongSpecies", 200, "number of bins in Delta plot (wrong species)"}; + Configurable nBinsTrackDeltaLength{"nBinsTrackDeltaLength", 100, "number of bins in delta track length"}; + Configurable nBinsNsigmaCorrectSpecies{"nBinsNsigmaCorrectSpecies", 200, "number of bins in Nsigma plot (correct speies)"}; + Configurable minDeltaRange{"minDeltaRange", -100, "lower limit for the nsigma axis"}; + Configurable maxDeltaRange{"maxDeltaRange", +100, "upper limit for the nsigma axis"}; + Configurable nBinsTimeRes{"nBinsTimeRes", 400, "number of bins plots time resolution"}; + Configurable nBinsRelativeEtaPt{"nBinsRelativeEtaPt", 400, "number of bins plots pt and eta relative errors"}; + Configurable nBinsEta{"nBinsEta", 400, "number of bins plot relative eta error"}; + Configurable nBinsMult{"nBinsMult", 200, "number of bins in multiplicity"}; + Configurable maxMultRange{"maxMultRange", 1000.f, "upper limit in multiplicity plots"}; + } plotsConfig; + + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; + + // Track smearer array, one per geometry + // std::vector> mSmearer; + std::vector> mSmearer; + // needed: random number generator for smearing + TRandom3 pRandomNumberGenerator; + + // for handling basic QA histograms if requested + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + OutputObj listEfficiency{"efficiency"}; + static constexpr int kParticles = 9; + + // Configuration defined at init time + o2::fastsim::GeometryContainer mGeoContainer; + float mMagneticField = 0.0f; + void init(o2::framework::InitContext& initContext) + { + mGeoContainer.setCcdbManager(ccdb.operator->()); + mGeoContainer.init(initContext); + + const int nGeometries = mGeoContainer.getNumberOfConfigurations(); + mMagneticField = mGeoContainer.getFloatValue(0, "global", "magneticfield"); + + pRandomNumberGenerator.SetSeed(0); // fully randomize + if (simConfig.flagMIDLoadDelphesLUTs) { + for (int icfg = 0; icfg < nGeometries; ++icfg) { + mSmearer.emplace_back(std::make_unique()); + mSmearer[icfg]->setCcdbManager(ccdb.operator->()); + std::map globalConfiguration = mGeoContainer.getConfiguration(icfg, "global"); + for (const auto& entry : globalConfiguration) { + int pdg = 0; + if (entry.first.find("lut") != 0) { + continue; + } + if (entry.first.find("lutEl") != std::string::npos) { + pdg = kElectron; + } else if (entry.first.find("lutMu") != std::string::npos) { + pdg = kMuonMinus; + } else if (entry.first.find("lutPi") != std::string::npos) { + pdg = kPiPlus; + } else if (entry.first.find("lutKa") != std::string::npos) { + pdg = kKPlus; + } else if (entry.first.find("lutPr") != std::string::npos) { + pdg = kProton; + } else if (entry.first.find("lutDe") != std::string::npos) { + pdg = o2::constants::physics::kDeuteron; + } else if (entry.first.find("lutTr") != std::string::npos) { + pdg = o2::constants::physics::kTriton; + } else if (entry.first.find("lutHe3") != std::string::npos) { + pdg = o2::constants::physics::kHelium3; + } else if (entry.first.find("lutAl") != std::string::npos) { + pdg = o2::constants::physics::kAlpha; + } + + std::string filename = entry.second; + if (pdg == 0) { + LOG(fatal) << "Unknown LUT entry " << entry.first << " for global configuration"; + } + LOG(info) << "Loading LUT for pdg " << pdg << " for config " << icfg << " from provided file '" << filename << "'"; + if (filename.empty()) { + LOG(warning) << "No LUT file passed for pdg " << pdg << ", skipping."; + } + // strip from leading/trailing spaces + filename.erase(0, filename.find_first_not_of(" ")); + filename.erase(filename.find_last_not_of(" ") + 1); + if (filename.empty()) { + LOG(warning) << "No LUT file passed for pdg " << pdg << ", skipping."; + } + bool success = mSmearer[icfg]->loadTable(pdg, filename.c_str()); + if (!success) { + LOG(fatal) << "Having issue with loading the LUT " << pdg << " " << filename; + } + } + } + } + + if (plotsConfig.doQAplots) { + const AxisSpec axisdNdeta{plotsConfig.nBinsMult, 0.0f, plotsConfig.maxMultRange, Form("dN/d#eta in |#eta| < %f", simConfig.multiplicityEtaRange.value)}; + + histos.add("h1dNdeta", "h2dNdeta", kTH1F, {axisdNdeta}); + histos.add("h2dEventTime", "h2dEventTime", kTH2F, {{2000, -10000, 10000, "computed"}, {2000, -10000, 10000, "generated"}}); + histos.add("h1dEventTimegen", "h1dEventTimegen", kTH1F, {{200, -1000, 1000, "generated"}}); + histos.add("h1dEventTimerec", "h1dEventTimerec", kTH1F, {{2000, -10000, 10000, "computed"}}); + histos.add("h1dEventTimedelta", "h1dEventTimedelta", kTH1F, {{200, -1000, 1000, "generated - computed"}}); + histos.add("h2dEventTimeres", "h2dEventTimeres", kTH2F, {axisdNdeta, {300, 0, 300, "resolution"}}); + listEfficiency.setObject(new THashList); + listEfficiency->Add(new TEfficiency("effEventTime", "effEventTime", plotsConfig.nBinsMult, 0.0f, plotsConfig.maxMultRange)); + + const AxisSpec axisMomentum{static_cast(plotsConfig.nBinsP), 0.0f, +10.0f, "#it{p} (GeV/#it{c})"}; + const AxisSpec axisRigidity{static_cast(plotsConfig.nBinsP), 0.0f, +10.0f, "#it{p} / |#it{z}| (GeV/#it{c})"}; + const AxisSpec axisVelocity{static_cast(plotsConfig.nBinsBeta), 0.0f, +1.1f, "Measured #beta"}; + const AxisSpec axisMomentumSmall{static_cast(plotsConfig.nBinsP), 0.0f, +1.0f, "#it{p} (GeV/#it{c})"}; + const AxisSpec axisTrackLength{static_cast(plotsConfig.nBinsTrackLength), 0.0f, 1000.0f, "Track length (cm)"}; + const AxisSpec axisTrackDeltaLength{static_cast(plotsConfig.nBinsTrackDeltaLength), 0.0f, 30.0f, "Delta Track length (cm)"}; + + histos.add("h2dVelocityVsMomentumInner", "h2dVelocityVsMomentumInner", kTH2F, {axisMomentum, axisVelocity}); + histos.add("h2dVelocityVsRigidityInner", "h2dVelocityVsRigidityInner", kTH2F, {axisRigidity, axisVelocity}); + histos.add("h2dTrackLengthInnerVsPt", "h2dTrackLengthInnerVsPt", kTH2F, {axisMomentum, axisTrackLength}); + histos.add("h2dTrackLengthInnerRecoVsPt", "h2dTrackLengthInnerRecoVsPt", kTH2F, {axisMomentum, axisTrackLength}); + + histos.add("h2dTrackLengthB2vsB0", "h2dDeltaTrackLengthInnerVsPt; for B=0; for B=2 T", kTH2F, {axisTrackLength, axisTrackLength}); + histos.add("h2HitMap", "h2HitMap", kTH2F, {{200, -simConfig.midLength / 2, simConfig.midLength / 2}, {1000, 0, simConfig.midRadius * 2 * M_PI}}); + histos.add("h2HitMapRec", "h2HitMapRec", kTH2F, {{200, -simConfig.midLength / 2, simConfig.midLength / 2}, {1000, 0, simConfig.midRadius * 2 * M_PI}}); + histos.add("h2HitDeltaXY", "h2HitDeltaXY", kTH2F, {{200, -100, 100, "#Delta#it{x} (cm)"}, {200, -100, 100, "#Delta#it{y} (cm)"}}); + + const AxisSpec axisPt{static_cast(plotsConfig.nBinsP), 0.0f, +10.0f, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec axisPtRes{1000, -1.0f, +1.0f, "(#it{p}_{T}^{gen}-#it{p}_{T}^{rec})/#it{p}_{T}^{gen}"}; + const AxisSpec axisEta{static_cast(plotsConfig.nBinsEta), -2.0f, +2.0f, "#eta"}; + const AxisSpec axisRelativePt{static_cast(plotsConfig.nBinsRelativeEtaPt), 0.0f, +10.0f, "#it{#sigma_{p_{T}}} / #it{p_{T}} (%)"}; + const AxisSpec axisRelativeEta{static_cast(plotsConfig.nBinsRelativeEtaPt), 0.0f, +10.0f, "#it{#sigma_{#eta}} / #it{#eta} (%)"}; + histos.add("h2dEtaVsPt", "h2dTrackLengthInnerRecoVsPt", kTH2F, {axisMomentum, axisEta}); + histos.add("h2dEtaVsP", "h2dTrackLengthInnerRecoVsP", kTH2F, {axisMomentum, axisEta}); + + const AxisSpec axisDelta{plotsConfig.nBinsDeltaCorrectSpecies, plotsConfig.minDeltaRange, plotsConfig.maxDeltaRange, "#Delta#{t} (ps)"}; + const AxisSpec axisNsigmaCorrect{plotsConfig.nBinsNsigmaCorrectSpecies, plotsConfig.minNsigmaRange, plotsConfig.maxNsigmaRange, "N#sigma - True "}; + const AxisSpec axisDeltaWrong{plotsConfig.nBinsDeltaWrongSpecies, plotsConfig.minDeltaRange, plotsConfig.maxDeltaRange, "#Delta - True "}; + const AxisSpec axisX{plotsConfig.doSeparationVsPt.value ? axisPt : axisMomentum}; + // const AxisSpec axisNSigma{axisNsigmaCorrect}; + + const AxisSpec axisTrackTimeRes{plotsConfig.nBinsTimeRes, 0.0f, +3000.0f, "Track time resolution (ps)"}; + const AxisSpec axisTotalTimeRes{plotsConfig.nBinsTimeRes, 0.0f, +3000.0f, "Total time resolution (ps)"}; + + histos.add("h1PtRes", "", kTH2F, {axisMomentum, axisPtRes}); + histos.add("h2EkinVsEloss", "", kTH2F, {{200, 0.0f, 2000.f, "#it{E}_{loss} (MeV)"}, {200, 0.0f, 20000.f, "KE (MeV)"}}); + } + } + + struct MIDEfficiency { + private: + const float layerRadius; + const float layerLength; + const float cellDimensionZ; + const float cellDimensionRPhi; + const float fractionInactive; + const float magField; + + TAxis* axisZ = nullptr; + TAxis* axisRPhi = nullptr; + TAxis* axisInCellZ = nullptr; + TAxis* axisInCellRPhi = nullptr; + + TH2F* hHitMapInCell = nullptr; + TH2F* hHitMapInCellBefore = nullptr; + TH2F* hHitMap = nullptr; + + public: + ~MIDEfficiency() + { + if (0) { + hHitMap->SaveAs(Form("/tmp/%s.png", hHitMap->GetName())); + hHitMapInCell->SaveAs(Form("/tmp/%s.png", hHitMapInCell->GetName())); + hHitMapInCellBefore->SaveAs(Form("/tmp/%s.png", hHitMapInCellBefore->GetName())); + } + + delete axisZ; + delete axisRPhi; + delete axisInCellZ; + delete axisInCellRPhi; + delete hHitMap; + delete hHitMapInCell; + delete hHitMapInCellBefore; + } + + MIDEfficiency(const MIDEfficiency&) = delete; + MIDEfficiency& operator=(const MIDEfficiency&) = delete; + + MIDEfficiency(float r, float l, std::array pDimensions, float fIA, float m) + : layerRadius(r), layerLength(l), cellDimensionZ(pDimensions[0]), cellDimensionRPhi(pDimensions[1]), fractionInactive(fIA), magField(m) + { + // Assuming square cells for simplicity + const float circumference = 2.0f * M_PI * layerRadius; + axisZ = new TAxis(static_cast(layerLength / cellDimensionZ), -layerLength / 2, layerLength); + axisRPhi = new TAxis(static_cast(circumference / cellDimensionRPhi), 0.f, circumference); + + const float inactiveBorderRPhi = cellDimensionRPhi * std::sqrt(fractionInactive) / 2; + const float inactiveBorderZ = cellDimensionZ * std::sqrt(fractionInactive) / 2; + const double arrayRPhi[4] = {-cellDimensionRPhi / 2, -cellDimensionRPhi / 2 + inactiveBorderRPhi, cellDimensionRPhi / 2 - inactiveBorderRPhi, cellDimensionRPhi / 2}; + for (int i = 0; i < 4; i++) { + LOG(info) << "arrayRPhi[" << i << "] = " << arrayRPhi[i]; + } + axisInCellRPhi = new TAxis(3, arrayRPhi); + const double arrayZ[4] = {-cellDimensionZ / 2, -cellDimensionZ / 2 + inactiveBorderZ, cellDimensionZ / 2 - inactiveBorderZ, cellDimensionZ / 2}; + for (int i = 0; i < 4; i++) { + LOG(info) << "arrayZ[" << i << "] = " << arrayZ[i]; + } + axisInCellZ = new TAxis(3, arrayZ); + + hHitMap = new TH2F(Form("hHitMap_R%.0f", layerRadius), "HitMap;z (cm); r#phi (cm)", 1000, -1000, 1000, 1000, -1000, 1000); + hHitMapInCell = new TH2F(Form("hHitMapInCell_R%.0f", layerRadius), "HitMapInCell;z (cm); r#phi (cm)", 1000, -10, 10, 1000, -10, 10); + hHitMapInCellBefore = new TH2F(Form("hHitMapInCellBefore_R%.0f", layerRadius), "HitMapInCell;z (cm); r#phi (cm)", 1000, -10, 10, 1000, -10, 10); + } + + bool isInMIDActiveArea(std::array hitPosition) + { + if (fractionInactive <= 0.0f) { + return true; + } + if (fractionInactive >= 1.0f) { + return false; + } + + // Convert 3D position to cylindrical coordinates for area calculation + const float phi = std::atan2(hitPosition[1], hitPosition[0]); + const float rphi = phi * layerRadius; + const float z = hitPosition[2]; + const float r = std::sqrt(hitPosition[0] * hitPosition[0] + hitPosition[1] * hitPosition[1]); + + // Check if hit is within layer geometric acceptance + if (std::abs(layerRadius - r) > 10.f) { + LOG(debug) << "Hit out of MID layer acceptance: r=" << r << " cm with respect to the layer radius " << layerRadius; + return false; + } + if (std::abs(z) > layerLength / 2.0f) { + LOG(debug) << "Hit out of MID layer acceptance: z=" << z << " cm with respect to the layer length " << layerLength; + return false; + } + + const int cellIndexPhi = axisRPhi->FindBin(rphi); + const int cellIndexZ = axisZ->FindBin(z); + + if (cellIndexPhi <= 0 || cellIndexPhi > axisRPhi->GetNbins() || cellIndexZ <= 0 || cellIndexZ > axisZ->GetNbins()) { + return false; + } + // Calculate local position within the pixel + const float localRPhi = (rphi - axisRPhi->GetBinCenter(cellIndexPhi)); + const float localZ = (z - axisZ->GetBinCenter(cellIndexZ)); + + // The difference between the hit and the pixel position cannot be greater than the size of the pixel + if (std::abs(localRPhi - axisRPhi->GetBinCenter(cellIndexPhi)) > axisRPhi->GetBinWidth(cellIndexPhi)) { + // LOG(warning) << "Local hit difference in phi is bigger than the pixel size"; + } + if (std::abs(localZ - axisZ->GetBinCenter(cellIndexZ)) > axisZ->GetBinWidth(cellIndexZ)) { + // LOG(warning) << "Local hit difference in z is bigger than the pixel size"; + } + hHitMapInCellBefore->Fill(localZ, localRPhi); + switch (axisInCellRPhi->FindBin(localRPhi)) { + case 0: + case 1: + case 3: + case 4: + return false; + } + switch (axisInCellZ->FindBin(localZ)) { + case 0: + case 1: + case 3: + case 4: + return false; + } + hHitMapInCell->Fill(localZ, localRPhi); + hHitMap->Fill(z, rphi); + return true; + } + + /// Check if a track hits the active area (convenience function) + /// \param track the track parameters for automatic hit position calculation + /// \return true if the hit is in the active area + bool isTrackInActiveArea(const o2::track::TrackParCov& track) + { + if (fractionInactive <= 0.f) { + return true; + } + float x, y, z; + if (!track.getXatLabR(layerRadius, x, magField)) { + LOG(debug) << "Could not propagate track to TOF layer at radius " << layerRadius << " cm"; + return false; + } + bool b; + ROOT::Math::PositionVector3D hit = track.getXYZGloAt(x, magField, b); + if (!b) { + LOG(debug) << "Could not get hit position at radius " << layerRadius << " cm"; + return false; + } + hit.GetCoordinates(x, y, z); + return isInMIDActiveArea(std::array{x, y, z}); + } + }; + + bool isInMIDActiveArea(const o2::track::TrackParCov& track) + { + static MIDEfficiency midLayer(simConfig.midRadius, simConfig.midLength, simConfig.midCellDimension, simConfig.midFractionOfInactiveArea, mMagneticField); + return midLayer.isTrackInActiveArea(track); + } + + struct TracksWithTime { + TracksWithTime(int pdgCode, + std::pair midTime, + std::pair trackLengthMID, + std::pair trackLengthInIron, + std::pair, std::array> posInMid, + std::pair momentum, + std::pair pseudorapidity, + float noSmearingPt) + : mPdgCode(pdgCode), + mMIDTime(midTime), + mTrackLengthMID(trackLengthMID), + mTrackLengthInIron(trackLengthInIron), + mPosInMid(posInMid), + mMomentum(momentum), + mPseudorapidity{pseudorapidity}, + mNoSmearingPt{noSmearingPt} {} + + int mPdgCode; + std::pair mMIDTime; // Measured time and expected resolution at the MID [ps] + std::pair mTrackLengthMID; // Measured and generated length at the MID [cm] + std::pair mTrackLengthInIron; // Measured and generated length within absorber [cm] + std::pair, std::array> mPosInMid; + std::pair mMomentum; // Measured momentum and uncertainty on the momentum [GeV/c] + std::pair mPseudorapidity; // Measured pseudorapidity and uncertainty on the pseudorapidity + float mNoSmearingPt; // No smearing pt + }; + + std::vector tracksWithTime; + bool eventTime(std::vector& tracks, + std::array& tzero) + { + + float sum = 0.; + float sumw = 0.; + + // Todo: check the different mass hypothesis iteratively + for (const auto& track : tracks) { + auto pdgInfo = pdg->GetParticle(track.mPdgCode); + if (pdgInfo == nullptr) { + continue; + } + const float mass = pdgInfo->Mass(); + const float mass2 = mass * mass; + const float tof = track.mMIDTime.first; // [ps] + if (tof <= 0) { + continue; + } + const float etof = track.mMIDTime.second; // [ps] + const float length = track.mTrackLengthMID.first; // [cm] + float p = track.mMomentum.first; // [GeV/c] + p *= std::abs(pdgInfo->Charge()) / 3.; // Total momentum + const float ep = track.mMomentum.second; // [GeV/c] + const float p2 = p * p; + const float lengthOverC = length * o2::constants::physics::invLightSpeedCm2PS; + const float texp = lengthOverC / p * std::sqrt(mass2 + p2); + // LOG(info) << "TOF: " << tof << " " << etof << " vs " << texp; + const float etexp = lengthOverC * mass2 / p2 / std::sqrt(mass2 + p2) * ep; + const float sigma = std::sqrt(etexp * etexp + etof * etof); + const float deltat = tof - texp; + + const float w = 1. / (sigma * sigma); + + sum += w * deltat; + sumw += w; + } + + static constexpr float kMaxEventTimeResolution = 200.f; + if (sumw <= 0. || tracks.size() <= 1 || std::sqrt(1. / sumw) > kMaxEventTimeResolution) { + tzero[0] = 0.; // [ps] + tzero[1] = kMaxEventTimeResolution; // [ps] + return false; + } + + tzero[0] = sum / sumw; + tzero[1] = std::sqrt(1. / sumw); + return true; + } + + /// returns track time resolution + /// \param pt the transverse momentum of the tarck + /// \param eta the pseudorapidity of the tarck + /// \param trackPtResolution the absolute resolution on pt + /// \param trackEtaResolution the absolute resolution on eta + /// \param mass the mass of the particle + /// \param detRadius the radius of the cylindrical layer + /// \param magneticField the magnetic field (along Z) + double calculateTrackTimeResolutionAdvanced(const float pt, + const float eta, + const float trackPtResolution, + const float trackEtaResolution, + const float mass, + const float detRadius, + const float magneticField) + { + // Compute tracking contribution to timing using the error propagation formula + // Uses light speed in m/ps, magnetic field in T (*0.1 for conversion kGauss -> T) + double a0 = mass * mass; + double a1 = 0.299792458 * (0.1 * magneticField) * (0.01 * o2::constants::physics::LightSpeedCm2NS / 1e+3); + double a2 = (detRadius * 0.01) * (detRadius * 0.01) * (0.299792458) * (0.299792458) * (0.1 * magneticField) * (0.1 * magneticField) / 2.0; + double dtofOndPt = (std::pow(pt, 4) * std::pow(std::cosh(eta), 2) * std::acos(1.0 - a2 / std::pow(pt, 2)) - 2.0 * a2 * std::pow(pt, 2) * (a0 + std::pow(pt * std::cosh(eta), 2)) / std::sqrt(a2 * (2.0 * std::pow(pt, 2) - a2))) / (a1 * std::pow(pt, 3) * std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2))); + double dtofOndEta = std::pow(pt, 2) * std::sinh(eta) * std::cosh(eta) * std::acos(1.0 - a2 / std::pow(pt, 2)) / (a1 * std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2))); + double trackTimeResolution = std::hypot(std::fabs(dtofOndPt) * trackPtResolution, std::fabs(dtofOndEta) * trackEtaResolution); + return trackTimeResolution; + } + + /* + // ============================================================ + // Track length in a B = 0 region + // + // Calculates the distance from the CURRENT track position + // to the cylinder of radius "radius". + // ============================================================ + */ + float computeTrackLengthB0(const o2::track::TrackParCov& track, + const float radius) + { + // Current global position + std::array pos{}; + track.getXYZGlo(pos); + + // Current global momentum + std::array mom{}; + if (!track.getPxPyPzGlo(mom)) { + return -100.f; + } + + const float px = mom[0]; + const float py = mom[1]; + // const float pz = mom[2]; + + const float pt = std::hypot(px, py); + + if (pt <= 0.f) { + return -100.f; + } + + // Unit vector in transverse plane + const float ux = px / pt; + const float uy = py / pt; + + const float b = pos[0] * ux + pos[1] * uy; + const float c = pos[0] * pos[0] + + pos[1] * pos[1] - + radius * radius; + + const float discriminant = b * b - c; + + if (discriminant < 0.f) { + return -100.f; + } + + const float sqrtDisc = std::sqrt(discriminant); + + const float t1 = -b + sqrtDisc; + const float t2 = -b - sqrtDisc; + + // We want the outward/forward intersection + float t = -1.f; + + if (t1 > 0.f && t2 > 0.f) { + t = std::min(t1, t2); + } else if (t1 > 0.f) { + t = t1; + } else if (t2 > 0.f) { + t = t2; + } else { + return -100.f; + } + + return t * std::sqrt(1.f + track.getTgl() * track.getTgl()); + } + + /* + // ============================================================ + // MID piecewise track length + // + // B = +20.0 kG R < R0 + // B = 0.0 kG R0 < R < R1 + // B = -14.3 kG R1 < R < R2 + // + // Returns total track length from current track position to R3. + // ============================================================ + */ + float computeTrackLengthMID(const o2::track::TrackParCov& inputTrack, + float R0, + float R1, + float R2, + float R3, + std::array& posMID, + std::array& deltaMID, + float& deltaposInIron) + { + float B0_m = +20.0f; // kG + float B1_m = 0.0f; // kG + float B2_m = 0.0f; + B2_m = -B0_m * pow(R0, 2) / (pow(R2, 2) - pow(R1, 2)); // kG + float B3_m = 0.0f; // kG + + // Work with a copy + auto track = inputTrack; + + float totalLength = 0.f; + float xAtR = 0.f; + + // ========================================================== + // REGION 1 + // current position -> R0 + // B = +20 kG + // ========================================================== + + const float L0 = + o2::upgrade::computeTrackLength(track, R0, B0_m); + + if (L0 < 0.f) { + return -100.f; + } + + totalLength += L0; + + // Find local x corresponding to global radius R0 + if (!track.getXatLabR(R0, + xAtR, + B0_m, + o2::track::DirOutward)) { + return -100.f; + } + + // Propagate track to R0 + if (!track.propagateTo(xAtR, B0_m)) { + return -100.f; + } + + // ========================================================== + // REGION 2 + // R0 -> R1 + // B = 0 + // ========================================================== + + const float L01 = computeTrackLengthB0(track, R1); + + if (L01 < 0.f) { + return -100.f; + } + + totalLength += L01; + + // Find local x at R1 assuming B = 0 + if (!track.getXatLabR(R1, + xAtR, + B1_m, + o2::track::DirOutward)) { + return -100.f; + } + + // Straight-line propagation to R1 + if (!track.propagateTo(xAtR, B1_m)) { + return -100.f; + } + + // ========================================================== + // REGION 3 (absorber) + // R1 -> R2 + // B = -14.3 kG + // ========================================================== + + std::array pos1{}; + track.getXYZGlo(pos1); + + float L12 = + o2::upgrade::computeTrackLength(track, R2, B2_m); + + if (L12 < 0.f) { + return -100.f; + } + totalLength += L12; + deltaposInIron = L12; + // Find local x corresponding to global radius R2 + if (!track.getXatLabR(R2, + xAtR, + B2_m, + o2::track::DirOutward)) { + return -100.f; + } + + // Propagate track to R2 + if (!track.propagateTo(xAtR, B2_m)) { + return -100.f; + } + + // ========================================================== + // REGION 2 + // R2 -> R3 + // B = 0 + // ========================================================== + + const float L23 = computeTrackLengthB0(track, R3); + + if (L23 < 0.f) { + return -100.f; + } + + totalLength += L23; + + // Find local x corresponding to global radius R3 + if (!track.getXatLabR(R3, + xAtR, + B3_m, + o2::track::DirOutward)) { + return -100.f; + } + + // Propagate track to R3 + if (!track.propagateTo(xAtR, B3_m)) { + return -100.f; + } + + // ========================================================== + // GLOBAL POSITION AT R3 + // ========================================================== + + std::array pos3{}; + track.getXYZGlo(pos3); + + posMID[0] = pos3[0]; + posMID[1] = pos3[1]; + posMID[2] = pos3[2]; + + // Displacement projections + deltaMID[0] = posMID[0] - pos1[0]; + deltaMID[1] = posMID[1] - pos1[1]; + deltaMID[2] = posMID[2] - pos1[2]; + + return totalLength; + } + // Multiple Coulomb scattering angle in iron + // + // momentum : particle momentum [GeV/c] + // thickness : iron thickness [cm] + // + // returns theta0 [rad] + + float moliereThetaIron(float momentum, float mass, float thickness) + { + float X0 = 1.757f; // radiation length of iron [cm] + // constexpr float mMuon = 0.105658f; // muon mass [GeV/c^2] + + if (momentum <= 0.f || thickness <= 0.f) { + return 0.f; + } + + // beta = p/E + const float energy = std::sqrt(momentum * momentum + + mass * mass); + + const float beta = momentum / energy; + + // Thickness in radiation lengths + const float xOverX0 = thickness / X0; + + // Highland approximation + const float theta0 = (0.0136f / (beta * momentum)) * sqrt(xOverX0) * (1.f + 0.038f * log(xOverX0)); + + return theta0; + } + float fPi(double p) + { + return p > 1. ? 1. - std::exp(-(p - 1.) / 1.5) : 0.; + } + float fMu(double p) + { + return p > 0.7 ? (1. - std::exp(-(p - 0.7) / 0.8)) * std::exp(-p / 3.) : 0.; + } + float refPi(double p) + { + if (p <= 1.) + return 0.; + float f3 = fPi(3.); + float mu10 = simConfig.absFakeAt3 * (1. - simConfig.absPionFractionAt3) * fMu(10.) / fMu(3.); + float slope = ((simConfig.absFakeAt10 - mu10) / (simConfig.absFakeAt3 * simConfig.absPionFractionAt3) / (fPi(10.) / f3) - 1.) / 7.; + return std::max(0., simConfig.absFakeAt3 * simConfig.absPionFractionAt3) * fPi(p) / f3 * (1. + slope * (p - 3.)); + } + float refMu(double p) + { + return simConfig.absFakeAt3 * (1. - simConfig.absPionFractionAt3) * fMu(p) / fMu(3.); + } + float decayAir(double p, double airCm) + { + if (p <= 0.) + return 0.; + return -std::expm1(-airCm * 0.13957039 / (p * 780.45)); // assume pion mass, ctau=780.45 cm + } + float survivalAir(double p, double airCm) + { + return 1. - decayAir(p, airCm); + } + // Range proxy: kinetic energy minus mean ionization loss in iron. + // Smoothed only to represent missing fluctuations; NOT a Geant4 range calculation. + float rangeTurnOn(double p, double ironCm, double mass) + { + if (p <= 0.) + return 0.; + float kinetic = std::sqrt(p * p + mass * mass) - mass; + float deficit = kinetic - (simConfig.absDensity * simConfig.absDeDx / 1000.) * ironCm; + float width = 0.10; // GeV, phenomenological + return 1. / (1. + std::exp(-std::max(-60., std::min(60., deficit / width)))); + } + float etaShape(double eta) + { + return std::exp(-0.5 * eta * eta / (simConfig.absSigmaEta * simConfig.absSigmaEta)); + } + float pion(float p, float eta, float airCm, float ironCm) + { + if (p <= 0.) + return 0.; + const float c = std::cosh(eta); + const float survivalRatio = survivalAir(p, airCm * c) / survivalAir(p, simConfig.absRadius1 * c); + const float hadronRatio = std::exp(-(ironCm - (simConfig.absRadius2 - simConfig.absRadius1)) * c / simConfig.absLambdaPi); + // Ratio to reference range at the SAME eta; reference shape is phenomenological. + const float rangeRatio = rangeTurnOn(p, ironCm * c, 0.13957039) / rangeTurnOn(p, (simConfig.absRadius2 - simConfig.absRadius1) * c, 0.13957039); // assume pion mass + return std::max(0., refPi(p) * etaShape(eta) * survivalRatio * hadronRatio * rangeRatio); + } + float muon(float p, float eta, float airCm, float ironCm) + { + if (p <= 0.) + return 0.; + const float c = std::cosh(eta); + const float d0 = decayAir(p, simConfig.absRadius1 * c); + const float decayRatio = d0 > 0. ? decayAir(p, airCm * c) / d0 : 0.; + // Effective secondary-muon momentum proxy: 0.8 * pion momentum. + // Secondary muons born in iron are NOT modeled; this is a phenomenological air-decay term. + const float pMu = 0.8 * p; + const float rangeRatio = rangeTurnOn(pMu, ironCm * c, 0.1056583755) / + rangeTurnOn(pMu, (simConfig.absRadius2 - simConfig.absRadius1) * c, 0.1056583755); + return std::max(0., refMu(p) * etaShape(eta) * decayRatio * rangeRatio); + } + float total(float p, float eta, float airCm, float ironCm) + { + return std::min(1., pion(p, eta, airCm, ironCm) + muon(p, eta, airCm, ironCm)); + } + float meanEta(float p, float airCm, float ironCm) + { + const int n = 1000; + float s = 0.; + for (int i = 0; i < n; ++i) { + float eta = -simConfig.absEtaMax + (i + .5) * 2. * simConfig.absEtaMax / n; + s += total(p, eta, airCm, ironCm); + } + return s / n; + } + + void process(soa::Join::iterator const& collision, + soa::Join const& tracks, + aod::McParticles const&, + aod::McCollisions const&) + { + o2::dataformats::VertexBase pvVtx({collision.posX(), collision.posY(), collision.posZ()}, + {collision.covXX(), collision.covXY(), collision.covYY(), + collision.covXZ(), collision.covYZ(), collision.covZZ()}); + + std::array mcPvCov = {0.}; + o2::dataformats::VertexBase mcPvVtx({0.0f, 0.0f, 0.0f}, mcPvCov); + const float eventCollisionTimePS = (simConfig.considerEventTime.value ? collision.collisionTime() * 1e3 : 0.f); // convert ns to ps + if (collision.has_mcCollision()) { + auto mcCollision = collision.mcCollision(); + mcPvVtx.setX(mcCollision.posX()); + mcPvVtx.setY(mcCollision.posY()); + mcPvVtx.setZ(mcCollision.posZ()); + } // else remains untreated for now + + // First we compute the number of charged particles in the event if LUTs are loaded + float dNdEta = 0.f; + for (const auto& track : tracks) { + if (!track.has_mcParticle()) + continue; + auto mcParticle = track.mcParticle(); + if (std::abs(mcParticle.eta()) > simConfig.multiplicityEtaRange) { + continue; + } + if (mcParticle.has_daughters()) { + continue; + } + const auto& pdgInfo = pdg->GetParticle(mcParticle.pdgCode()); + if (!pdgInfo) { + // LOG(warning) << "PDG code " << mcParticle.pdgCode() << " not found in the database"; + continue; + } + if (pdgInfo->Charge() == 0) { + continue; + } + dNdEta += 1.f; + } + if (plotsConfig.doQAplots) { + histos.fill(HIST("h1dNdeta"), dNdEta); + } + + tracksWithTime.clear(); // clear the vector of tracks with time to prepare the cache for the next event + tracksWithTime.reserve(tracks.size()); + // First loop to generate the arrival time of the particles to the MID layers + for (const auto& track : tracks) { + // first step: find precise arrival time (if any) + // --- convert track into perfect track + if (!track.has_mcParticle()) { // should always be OK but check please + upgradeMidMC(-999.f, -999.f); + continue; + } else { + LOG(debug) << "Track without mcParticle found!"; + } + const auto& mcParticle = track.mcParticle(); + float true_pt = mcParticle.pt(); + o2::track::TrackParCov o2track = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdg); + // get mass to calculate velocity + auto pdgInfo = pdg->GetParticle(mcParticle.pdgCode()); + if (pdgInfo == nullptr) { + LOG(error) << "PDG code " << mcParticle.pdgCode() << " not found in the database"; + upgradeMidMC(-999.f, -999.f); + continue; + } + + float xPv = -100.f; + static constexpr float kTrkXThreshold = -99.f; // Threshold to consider a good propagation of the track + if (o2track.propagateToDCA(mcPvVtx, mMagneticField)) { + xPv = o2track.getX(); + } + float trackLengthMID = -1; + float trackLengthMIDB2 = -1; + std::array posMid{-999., -999., -999.}; + std::array delaPosMid{-999., -999., -999.}; + float trackLengthInIron = -999.; + if (xPv > kTrkXThreshold) { + trackLengthMIDB2 = o2::upgrade::computeTrackLength(o2track, simConfig.magRadiusC, mMagneticField); + trackLengthMID = computeTrackLengthMID(o2track, simConfig.magRadiusC, simConfig.absRadius1, simConfig.absRadius2, simConfig.midRadius, posMid, delaPosMid, trackLengthInIron); + // thetamolInIron = moliereThetaIron(o2track.getP(), pdgInfo->Mass(), trackLengthInIron); + } + if (plotsConfig.doQAplots) { + histos.fill(HIST("h2dTrackLengthB2vsB0"), trackLengthMIDB2, trackLengthMID); + } + + // Check if the track hit a sensitive area of the MID + const bool activeMID = isInMIDActiveArea(o2track); + if (!activeMID) { + trackLengthMID = -999.f; + } else { + if (plotsConfig.doQAplots) { + float phi_tmp = std::atan2(posMid[1], posMid[0]); + if (phi_tmp < 0.f) { + phi_tmp += 2.f * M_PI; + } + histos.fill(HIST("h2HitMap"), posMid[2], simConfig.midRadius * phi_tmp); + } + } + + const float v = o2::upgrade::computeParticleVelocity(o2track.getP(), pdgInfo->Mass()); + const float expectedTimeMID = trackLengthMID > 0 ? trackLengthMID / v + eventCollisionTimePS : -999.f; // arrival time to the MID in ps + + upgradeMidMC(expectedTimeMID, trackLengthMID); + + // Smear with expected resolutions + const float measuredTimeMID = pRandomNumberGenerator.Gaus(expectedTimeMID, simConfig.midTimeReso); + + // Now we calculate the expected arrival time following certain mass hypotheses + // and the (imperfect!) reconstructed track parametrizations + float trackLengthRecoMID = -1; + auto recoTrack = getTrackParCov(track); + xPv = -100.f; + std::array posRecMid{-999., -999., -999.}; + std::array deltaPosRecMid{-999., -999., -999.}; + float trackLengthRecoInIron = -1.; + trackLengthRecoInIron = -1.; + trackLengthRecoMID = -1.; + if (recoTrack.propagateToDCA(pvVtx, mMagneticField)) { + xPv = recoTrack.getX(); + } + if (xPv > kTrkXThreshold) { + trackLengthRecoMID = computeTrackLengthMID(recoTrack, simConfig.magRadiusC, simConfig.absRadius1, simConfig.absRadius2, simConfig.midRadius, posRecMid, deltaPosRecMid, trackLengthRecoInIron); + // mimic scattering in iron + std::array mom{-999., -999., -999.}; + recoTrack.getPxPyPzGlo(mom); + for (int jj = 0; jj < 3; ++jj) { + float posRectmp = 0.0f; + float molieretmp = 0.0f; + molieretmp = moliereThetaIron(mom[jj], pdgInfo->Mass(), deltaPosRecMid[jj]); + posRectmp = pRandomNumberGenerator.Gaus(posRecMid[jj], molieretmp); + posRecMid[jj] = posRectmp; + } + if (plotsConfig.doQAplots) { + float phi_tmp = std::atan2(posRecMid[1], posRecMid[0]); + if (phi_tmp < 0.f) { + phi_tmp += 2.f * M_PI; + } + histos.fill(HIST("h2HitMapRec"), posRecMid[2], simConfig.midRadius * phi_tmp); + histos.fill(HIST("h2HitDeltaXY"), posRecMid[0] - posMid[0], posRecMid[1] - posMid[1]); + histos.fill(HIST("h1PtRes"), true_pt, (true_pt - recoTrack.getPt()) / true_pt); + } + } + // cache the track info needed for the event time calculation + // Reuse or emplace a new object in the vector + tracksWithTime.emplace_back(TracksWithTime{mcParticle.pdgCode(), + {measuredTimeMID, expectedTimeMID}, + {trackLengthRecoMID, trackLengthMID}, + {trackLengthRecoInIron, trackLengthInIron}, + {posRecMid, posMid}, + {recoTrack.getP(), recoTrack.getSigma1Pt2()}, + {recoTrack.getEta(), recoTrack.getSigmaTgl2()}, + o2track.getPt()}); + } + // Now we compute the event time for the tracks + + std::array tzero = {0.f, 0.f}; + bool etStatus = false; + if (simConfig.considerEventTime.value) { + etStatus = eventTime(tracksWithTime, tzero); + if (!etStatus) { + LOG(debug) << "Event time calculation failed with " << tracksWithTime.size() << " tracks with time and " << dNdEta << " charged particles"; + } + } + + if (plotsConfig.doQAplots) { + histos.fill(HIST("h2dEventTime"), tzero[0], eventCollisionTimePS); + histos.fill(HIST("h1dEventTimegen"), eventCollisionTimePS); + histos.fill(HIST("h1dEventTimerec"), tzero[0]); + histos.fill(HIST("h2dEventTimeres"), dNdEta, tzero[1]); + if (etStatus) { + histos.fill(HIST("h1dEventTimedelta"), eventCollisionTimePS - tzero[0]); + } + static_cast(listEfficiency->At(0))->Fill(etStatus, dNdEta); + } + + // Then we do a second loop to compute the measured quantities with the measured event time + int trackWithTimeIndex = 0; + for (const auto& track : tracks) { + if (!track.has_mcParticle()) { // should always be OK but check please + continue; + } + const auto& mcParticle = track.mcParticle(); + + const auto& trkWithTime = tracksWithTime[trackWithTimeIndex++]; + const float trackLengthRecoMID = trkWithTime.mTrackLengthMID.first; + const float trackLengthMID = trkWithTime.mTrackLengthMID.second; + const float measuredTimeMID = trkWithTime.mMIDTime.first; + const float expecTimeMID = trkWithTime.mMIDTime.second; + const float momentum = trkWithTime.mMomentum.first; + const float pseudorapidity = trkWithTime.mPseudorapidity.first; + const float noSmearingPt = trkWithTime.mNoSmearingPt; + std::array expectedPosMid = trkWithTime.mPosInMid.second; + std::array measuredPosMid = trkWithTime.mPosInMid.first; + float trackLengthInIronRec = trkWithTime.mTrackLengthInIron.first; + bool isMuCandidate = false; + if (abs(pseudorapidity) < simConfig.absEtaMax) { + if (trackLengthInIronRec > 0) { + float kineticEnergy = std::sqrt(pow(momentum, 2) + pow(0.1057, 2)) - 0.1057; + float energyloss = trackLengthInIronRec * simConfig.absDeDx * simConfig.absDensity / 1000.; + if (abs(trkWithTime.mPdgCode) == 13) { + if ((kineticEnergy - energyloss) > 0.1) { + if (plotsConfig.doQAplots) { + histos.fill(HIST("h2EkinVsEloss"), energyloss * 1000., kineticEnergy * 1000.); + } + isMuCandidate = true; + } + } else if (abs(trkWithTime.mPdgCode) == 11) { + isMuCandidate = false; + } else { // mostly hadrons, parametrization from simulations using pions, GEANT4+BDT + float eff_rndm = pRandomNumberGenerator.Uniform(0., 1.); + float eff_had = total(momentum, pseudorapidity, trackLengthRecoMID - trackLengthInIronRec, trackLengthInIronRec); + if (eff_rndm < eff_had && (kineticEnergy - energyloss) > 0.1) { + if (plotsConfig.doQAplots) { + histos.fill(HIST("h2EkinVsEloss"), energyloss * 1000., kineticEnergy * 1000.); + } + isMuCandidate = true; + } + } + } + } + auto truePdgInfo = pdg->GetParticle(mcParticle.pdgCode()); + float rigidity = momentum; // fallback to momentum if charge unknown + + // Use MC truth charge for rigidity calculation + if (truePdgInfo) { + const float trueCharge = std::abs(truePdgInfo->Charge()) / 3.0f; + if (trueCharge > 0) { + rigidity = momentum / trueCharge; + } + } + + if (plotsConfig.doQAplots) { + // unit conversion: length in cm, time in ps + const float midBeta = (trackLengthRecoMID / measuredTimeMID) / o2::constants::physics::LightSpeedCm2PS; + if (trackLengthRecoMID > 0) { + histos.fill(HIST("h2dVelocityVsMomentumInner"), momentum, midBeta); + histos.fill(HIST("h2dVelocityVsRigidityInner"), rigidity, midBeta); + histos.fill(HIST("h2dTrackLengthInnerVsPt"), noSmearingPt, trackLengthMID); + histos.fill(HIST("h2dTrackLengthInnerRecoVsPt"), noSmearingPt, trackLengthRecoMID); + histos.fill(HIST("h2dEtaVsPt"), noSmearingPt, pseudorapidity); + histos.fill(HIST("h2dEtaVsP"), momentum, pseudorapidity); + } + } + std::vector measuredPosMidVec(measuredPosMid.begin(), measuredPosMid.end()); + upgradeMid(tzero[0], tzero[1], measuredPosMidVec, measuredTimeMID, trackLengthRecoMID, isMuCandidate); + + std::vector expectedPosMidVec(expectedPosMid.begin(), expectedPosMid.end()); + upgradeMidExpected(expectedPosMidVec, expecTimeMID, trackLengthMID); + } + + if (trackWithTimeIndex != tracks.size()) { + LOG(fatal) << "Track with time index " << trackWithTimeIndex << " does not match the number of tracks " << tracks.size(); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; } diff --git a/ALICE3/Tasks/CMakeLists.txt b/ALICE3/Tasks/CMakeLists.txt index 6e6fbc37ea0..e86859ec374 100644 --- a/ALICE3/Tasks/CMakeLists.txt +++ b/ALICE3/Tasks/CMakeLists.txt @@ -44,6 +44,11 @@ o2physics_add_dpl_workflow(alice3-prefilterdilepton PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::FrameworkPhysicsSupport COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(alice3-prefilterdimuon + SOURCES alice3-prefilterdimuon.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::FrameworkPhysicsSupport + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(alice3-taskcorrelationddbar SOURCES alice3-taskcorrelationDDbar.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/ALICE3/Tasks/alice3-prefilterdimuon.cxx b/ALICE3/Tasks/alice3-prefilterdimuon.cxx new file mode 100644 index 00000000000..c797296b4a6 --- /dev/null +++ b/ALICE3/Tasks/alice3-prefilterdimuon.cxx @@ -0,0 +1,418 @@ +// 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 alice3-dimuon.cxx +/// \author s.scheid@cern.ch, daiki.sekihata@cern.ch +/// + +#include "ALICE3/DataModel/OTFMID.h" +#include "ALICE3/DataModel/prefilterDilepton.h" +#include "ALICE3/DataModel/tracksAlice3.h" +#include "Common/Core/RecoDecay.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include +#include + +#include +#include +#include +#include + +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::soa; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct Alice3DimuonEventCentSelection { + + Service inspdg; + + Produces eventCentSel; + Configurable minFITPart{"minFITPart", 0.f, "Minimum number of charged particles in the FIT acceptance"}; + Configurable maxFITPart{"maxFITPart", 0.f, "Maximum number of charged particles in the FIT acceptance"}; + + HistogramRegistry registry{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(InitContext&) + { + registry.add("Generated/Before/ParticlesFit", "Charged Particles in Fit acceptance per event", kTH1F, {{15000, 0, 15000}}); + registry.add("Generated/After/ParticlesFit", "Charged Particles in Fit acceptance per event", kTH1F, {{15000, 0, 15000}}); + registry.add("Generated/Before/ParticlesEta", "Charged Particles per event", kTH1F, {{1200, -6, 6}}); + registry.add("Generated/After/ParticlesEta", "Charged Particles per event", kTH1F, {{1200, -6, 6}}); + } + + using MyEvents = soa::Join; + Preslice perMCCollision = o2::aod::mcparticle::mcCollisionId; + + void processGen(MyEvents::iterator const& event, o2::aod::McCollisions const&, aod::McParticles const& mcParticles) + { + if (!event.has_mcCollision()) { + eventCentSel(0); + } else { + auto mccollision = event.mcCollision(); + auto mcParticles_per_coll = mcParticles.sliceBy(perMCCollision, mccollision.globalIndex()); + int nParticlesFIT = 0; + for (const auto& mcParticle : mcParticles_per_coll) { + if (mcParticle.isPhysicalPrimary()) { + auto pdgParticle = inspdg->GetParticle(mcParticle.pdgCode()); + if (pdgParticle) { + float charge = pdgParticle->Charge() / 3.f; // Charge in units of |e| + if (std::abs(charge) >= 1.) { + registry.fill(HIST("Generated/Before/ParticlesEta"), mcParticle.eta()); + if ((2.2 < mcParticle.eta() && mcParticle.eta() < 5.0) || (-3.4 < mcParticle.eta() && mcParticle.eta() < -2.3)) { + nParticlesFIT++; + } + } + } + } + } // end of mc particle loop + registry.fill(HIST("Generated/Before/ParticlesFit"), nParticlesFIT); + if (nParticlesFIT > minFITPart && nParticlesFIT < maxFITPart) { + registry.fill(HIST("Generated/After/ParticlesFit"), nParticlesFIT); + for (const auto& mcParticle : mcParticles_per_coll) { + if (mcParticle.isPhysicalPrimary()) { + auto pdgParticle = inspdg->GetParticle(mcParticle.pdgCode()); + if (pdgParticle) { + float charge = pdgParticle->Charge() / 3.f; // Charge in units of |e| + if (std::abs(charge) >= 1.) { + registry.fill(HIST("Generated/After/ParticlesEta"), mcParticle.eta()); + } + } + } + } // end of mc particle loop + eventCentSel(1); + } else { + eventCentSel(0); + } + } // if mc collision loop + } // end of process + + void processDummy(MyEvents::iterator const&) + { + eventCentSel(1); + } + + PROCESS_SWITCH(Alice3DimuonEventCentSelection, processGen, "Select centrality", false); + PROCESS_SWITCH(Alice3DimuonEventCentSelection, processDummy, "Dummy", true); +}; + +struct Alice3DimuonPrefilter { + + Produces pfb_derived; + + SliceCache cache_mc; + SliceCache cache_rec; + + Service inspdg; + + Configurable ptMin{"ptMin", 0.f, "Lower limit in pT"}; + Configurable ptMax{"ptMax", 5.f, "Upper limit in pT"}; + Configurable etaMin{"etaMin", -5.f, "Lower limit in eta"}; + Configurable etaMax{"etaMax", 5.f, "Upper limit in eta"}; + Configurable maxMass{"maxMass", 5.f, "Upper limit in mass"}; + Configurable maxOp{"maxOp", 5.f, "Upper limit in opening angle"}; + Configurable isMuMID{"isMuMID", true, "Muon candidate reaching MID"}; + Configurable nBinsP{"nBinsP", 80, "number of bins in momentum/pT depending on doSeparationVsPt"}; + Configurable nBinsEta{"nBinsEta", 400, "number of bins plot relative eta error"}; + std::unordered_map map_pfb; // map track.globalIndex -> prefilter bit + + HistogramRegistry registry{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + using MyTracksWithSmearings = soa::Join; + using MyTracksWithSmearing = MyTracksWithSmearings::iterator; + using MyTracksMCs = soa::Join; + using MyTracksMC = MyTracksMCs::iterator; + using DileptonCollisions = soa::Join; + using DileptonCollision = DileptonCollisions::iterator; + + Preslice perCollision = aod::track::collisionId; + Partition posTracks = o2::aod::track::signed1Pt > 0.f; + Partition negTracks = o2::aod::track::signed1Pt < 0.f; + Partition posTracksWithSmearing = o2::aod::track::signed1Pt > 0.f; + Partition negTracksWithSmearing = o2::aod::track::signed1Pt < 0.f; + + void init(InitContext&) + { + const AxisSpec axisM{500, 0, 5, "#it{m}_{ll} (GeV/#it{c}^{2})"}; + // const AxisSpec axisPt{1000, 0, 10, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec axisPt{nBinsP, 0.0f, +10.0f, "#it{p}_{T} (GeV/#it{c})"}; + // const AxisSpec axisSigmaEl{200, -10, 10, "n#sigma_{El}"}; + const AxisSpec axisTrackLengthMID{310, 0., 310., "Track length (cm)"}; + // const AxisSpec axisEta{1000, -5, 5, "#it{#eta}"}; + const AxisSpec axisEta{static_cast(nBinsEta), -2.0f, +2.0f, "#eta"}; + const AxisSpec axisPhi{360, 0, TMath::TwoPi(), "#it{#varphi} (rad.)"}; + const AxisSpec axisProdx{2000, -100, 100, "Prod. Vertex X (cm)"}; + const AxisSpec axisPrody{2000, -100, 100, "Prod. Vertex Y (cm)"}; + const AxisSpec axisProdz{2000, -100, 100, "Prod. Vertex Z (cm)"}; + + registry.add("Reconstructed/Pair/ULS/Mass_Pt", "Pair Mass vs. Pt", kTH2F, {axisM, axisPt}, true); + registry.add("ReconstructedFiltered/Pair/ULS/Mass_Pt", "Pair Mass vs. Pt", kTH2F, {axisM, axisPt}, true); + registry.add("Reconstructed/Track/Pt", "Track Pt", kTH1F, {axisPt}); + registry.add("Reconstructed/Track/Eta_P", "Track Eta vs Track P", kTH2F, {axisPt, axisEta}); + registry.add("Reconstructed/Track/Eta", "Track Eta", kTH1F, {axisEta}); + registry.add("Reconstructed/Track/Phi", "Track Phi", kTH1F, {axisPhi}); + registry.add("Reconstructed/Track/Eta_Pt", "Eta vs. Pt", kTH2F, {axisPt, axisEta}, true); + + if (doprocessPreFilter) { + // registry.add("Reconstructed/Track/SigmaOTofvspt", "Track #sigma oTOF", kTH2F, {axisPt, axisSigmaEl}); + // registry.add("Reconstructed/Track/SigmaITofvspt", "Track #sigma iTOF", kTH2F, {axisPt, axisSigmaEl}); + // registry.add("Reconstructed/Track/SigmaRichvspt", "Track #sigma RICH", kTH2F, {axisPt, axisSigmaEl}); + registry.add("Reconstructed/Track/midTrackLength", "Track length MID", kTH1F, {axisTrackLengthMID}); + } + } + + template + void PreFilter(TTracks const& posTracks_coll, TTracks const& negTracks_coll) + { + for (const auto& pos : posTracks_coll) { + if constexpr (isWithSmearing) { + if (!pos.isReconstructed() || !pos.selected()) { + continue; + } + if (pos.etaSmeared() < etaMin || etaMax < pos.etaSmeared()) { + continue; + } + if (pos.ptSmeared() < ptMin || ptMax < pos.ptSmeared()) { + continue; + } + registry.fill(HIST("Reconstructed/Track/Pt"), pos.ptSmeared()); + registry.fill(HIST("Reconstructed/Track/Eta"), pos.etaSmeared()); + // registry.fill(HIST("Reconstructed/Track/Eta_P"), pos.momentumSmeared(),pos.etaSmeared()); + registry.fill(HIST("Reconstructed/Track/Phi"), pos.phiSmeared()); + registry.fill(HIST("Reconstructed/Track/Eta_Pt"), pos.ptSmeared(), pos.etaSmeared()); + } else { + if (!pos.isReconstructed()) { + continue; + } + if (pos.eta() < etaMin || etaMax < pos.eta()) { + continue; + } + if (pos.pt() < ptMin || ptMax < pos.pt()) { + continue; + } + if (pos.midCandidate()) { + registry.fill(HIST("Reconstructed/Track/midTrackLength"), pos.midTrackLength()); + registry.fill(HIST("Reconstructed/Track/Pt"), pos.pt()); + registry.fill(HIST("Reconstructed/Track/Eta"), pos.eta()); + registry.fill(HIST("Reconstructed/Track/Eta_P"), pos.p(), pos.eta()); + registry.fill(HIST("Reconstructed/Track/Phi"), pos.phi()); + registry.fill(HIST("Reconstructed/Track/Eta_Pt"), pos.pt(), pos.eta()); + } + } + } + for (const auto& ele : negTracks_coll) { + if constexpr (isWithSmearing) { + if (!ele.isReconstructed() || !ele.selected()) { + continue; + } + if (ele.etaSmeared() < etaMin || etaMax < ele.etaSmeared()) { + continue; + } + if (ele.ptSmeared() < ptMin || ptMax < ele.ptSmeared()) { + continue; + } + registry.fill(HIST("Reconstructed/Track/Pt"), ele.ptSmeared()); + registry.fill(HIST("Reconstructed/Track/Eta"), ele.etaSmeared()); + // registry.fill(HIST("Reconstructed/Track/Eta_P"), ele.pSmeared(),ele.etaSmeared()); + registry.fill(HIST("Reconstructed/Track/Phi"), ele.phiSmeared()); + registry.fill(HIST("Reconstructed/Track/Eta_Pt"), ele.ptSmeared(), ele.etaSmeared()); + } else { + if (!ele.isReconstructed()) { + continue; + } + if (ele.eta() < etaMin || etaMax < ele.eta()) { + continue; + } + if (ele.pt() < ptMin || ptMax < ele.pt()) { + continue; + } + if (ele.midCandidate()) { + registry.fill(HIST("Reconstructed/Track/midTrackLength"), ele.midTrackLength()); + registry.fill(HIST("Reconstructed/Track/Pt"), ele.pt()); + registry.fill(HIST("Reconstructed/Track/Eta"), ele.eta()); + registry.fill(HIST("Reconstructed/Track/Eta_P"), ele.p(), ele.eta()); + registry.fill(HIST("Reconstructed/Track/Phi"), ele.phi()); + registry.fill(HIST("Reconstructed/Track/Eta_Pt"), ele.pt(), ele.eta()); + } + } + } + + for (const auto& [pos, ele] : combinations(CombinationsFullIndexPolicy(posTracks_coll, negTracks_coll))) { // ULS + if constexpr (isWithSmearing) { + if (!pos.isReconstructed() || !pos.selected()) { + continue; + } + if (pos.etaSmeared() < etaMin || etaMax < pos.etaSmeared()) { + continue; + } + if (pos.ptSmeared() < ptMin || ptMax < pos.ptSmeared()) { + continue; + } + if (!ele.isReconstructed() || !ele.selected()) { + continue; + } + if (ele.etaSmeared() < etaMin || etaMax < ele.etaSmeared()) { + continue; + } + if (ele.ptSmeared() < ptMin || ptMax < ele.ptSmeared()) { + continue; + } + } else { + if (!pos.isReconstructed()) { + continue; + } + if (pos.eta() < etaMin || etaMax < pos.eta()) { + continue; + } + if (pos.pt() < ptMin || ptMax < pos.pt()) { + continue; + } + if (!ele.isReconstructed()) { + continue; + } + if (ele.eta() < etaMin || etaMax < ele.eta()) { + continue; + } + if (ele.pt() < ptMin || ptMax < ele.pt()) { + continue; + } + } + if constexpr (isWithSmearing) { + + ROOT::Math::PtEtaPhiMVector v1(pos.ptSmeared(), pos.etaSmeared(), pos.phiSmeared(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(ele.ptSmeared(), ele.etaSmeared(), ele.phiSmeared(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float angle = RecoDecay::constrainAngle(ROOT::Math::VectorUtil::Angle(v1, v2), -M_PI); + // o2::math_utils::bringToPMPi(angle); + + registry.fill(HIST("Reconstructed/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt()); + + if (v12.M() < maxMass && std::abs(angle) < maxOp) { + map_pfb[pos.globalIndex()] = 1; + map_pfb[ele.globalIndex()] = 1; + registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt()); + } + + } else { + + if (pos.midCandidate()) { + if (ele.midCandidate()) { + + ROOT::Math::PtEtaPhiMVector v1(pos.pt(), pos.eta(), pos.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(ele.pt(), ele.eta(), ele.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float angle = RecoDecay::constrainAngle(ROOT::Math::VectorUtil::Angle(v1, v2), -M_PI); + // o2::math_utils::bringToPMPi(angle); + + registry.fill(HIST("Reconstructed/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt()); + + if (v12.M() < maxMass && std::abs(angle) < maxOp) { + map_pfb[pos.globalIndex()] = 1; + map_pfb[ele.globalIndex()] = 1; + registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt()); + } + } + } + } + } // combination + } + + void processPreFilter(DileptonCollisions const& collisions, MyTracksMCs const& tracks) + { + for (const auto& track : tracks) { + map_pfb[track.globalIndex()] = 0; + } // end of track loop + + for (const auto& collision : collisions) { + auto negTracks_coll = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache_rec); + auto posTracks_coll = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache_rec); + + if (collision.isEventCentSelected() == 0) { + for (const auto& pos : posTracks_coll) { + map_pfb[pos.globalIndex()] = 0; + } + for (const auto& neg : negTracks_coll) { + map_pfb[neg.globalIndex()] = 0; + } + continue; + } + PreFilter(posTracks_coll, negTracks_coll); + } + for (const auto& track : tracks) { + pfb_derived(map_pfb[track.globalIndex()]); + } // end of track loop + map_pfb.clear(); + } + + void processPreFilterWithSmearing(DileptonCollisions const& collisions, MyTracksWithSmearings const& tracks) + { + for (const auto& track : tracks) { + map_pfb[track.globalIndex()] = 0; + } // end of track loop + + for (const auto& collision : collisions) { + auto negTracks_coll = negTracksWithSmearing->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache_rec); + auto posTracks_coll = posTracksWithSmearing->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache_rec); + + if (collision.isEventCentSelected() == 0) { + for (const auto& pos : posTracks_coll) { + map_pfb[pos.globalIndex()] = 0; + } + for (const auto& neg : negTracks_coll) { + map_pfb[neg.globalIndex()] = 0; + } + continue; + } + PreFilter(posTracks_coll, negTracks_coll); + } + for (const auto& track : tracks) { + pfb_derived(map_pfb[track.globalIndex()]); + } // end of track loop + map_pfb.clear(); + } + + void processDummy(MyTracksMCs const& tracks) + { + for (int i = 0; i < tracks.size(); i++) { + pfb_derived(0); + } + } + + PROCESS_SWITCH(Alice3DimuonPrefilter, processPreFilter, "Run prefilter", false); + PROCESS_SWITCH(Alice3DimuonPrefilter, processPreFilterWithSmearing, "Run prefilter with smearing", false); + PROCESS_SWITCH(Alice3DimuonPrefilter, processDummy, "Dummy", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +}