diff --git a/PWGMM/Lumi/Tasks/CMakeLists.txt b/PWGMM/Lumi/Tasks/CMakeLists.txt index 87b7e7fe919..4003fa5c2d9 100644 --- a/PWGMM/Lumi/Tasks/CMakeLists.txt +++ b/PWGMM/Lumi/Tasks/CMakeLists.txt @@ -39,4 +39,9 @@ o2physics_add_dpl_workflow(lumistab SOURCES lumiStability.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(lumi-stability-light-ions + SOURCES lumiStabilityLightIons.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCCDB O2Physics::AnalysisCore COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/PWGMM/Lumi/Tasks/lumiStabilityLightIons.cxx b/PWGMM/Lumi/Tasks/lumiStabilityLightIons.cxx new file mode 100644 index 00000000000..6bff87598a7 --- /dev/null +++ b/PWGMM/Lumi/Tasks/lumiStabilityLightIons.cxx @@ -0,0 +1,310 @@ +// 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 lumiStabilityLightIons.cxx +/// \brief Analysis over BCs to study the luminosity stability along time +/// +/// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt, Stefanie Mrozinski (stefanie.mrozinski@cern.ch) - Goethe University Frankfurt + +#include "Common/CCDB/ctpRateFetcher.h" +#include "Common/Core/MetadataHelper.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" + +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsParameters/AggregatedRunInfo.h" +#include "DataFormatsParameters/GRPLHCIFData.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +o2::common::core::MetadataHelper metadataInfo; // Metadata helper + +using MyBCs = soa::Join; + +struct LumiStabilityLightIons { + Configurable cfgRequireGoodRCTQuality{"cfgRequireGoodRCTQuality", false, "Only store BCs with good quality of FT0 in RCT"}; + Configurable cfgDoFT0Vtx{"cfgDoFT0Vtx", true, "Create and fill histograms for the FT0 vertex trigger"}; + Configurable cfgDoFT0CE{"cfgDoFT0CE", true, "Create and fill histograms for the FT0 centrality trigger"}; + Configurable cfgDoFDD{"cfgDoFDD", true, "Create and fill histograms for the FDD trigger"}; + Configurable cfgDo1ZNC{"cfgDo1ZNC", true, "Create and fill histograms for the 1ZNC trigger"}; + + Configurable cfgDoBCA{"cfgDoBCA", false, "Create and fill histograms for the BCs of type A"}; + Configurable cfgDoBCB{"cfgDoBCB", true, "Create and fill histograms for the BCs of type B"}; + Configurable cfgDoBCC{"cfgDoBCC", false, "Create and fill histograms for the BCs of type C"}; + Configurable cfgDoBCE{"cfgDoBCE", false, "Create and fill histograms for the BCs of type E"}; + Configurable cfgDoBCL{"cfgDoBCL", false, "Create and fill histograms for leading BCs of type B"}; + + Configurable cfgEmptyBCsBeforeLeadingBC{"cfgEmptyBCsBeforeLeadingBC", 5, "Minimum number of empty BCs before a leading BC to identify it as such"}; + + std::bitset beamPatternA, beamPatternC; + std::bitset bcPatternA, bcPatternC, bcPatternB, bcPatternE; + + std::string strLPMProductionTag = ""; // MC production tag to be retrieved from AO2D metadata + + const int nBCsPerOrbit = 3564; + + aod::rctsel::RCTFlagsChecker isFT0GoodRCTChecker{aod::rctsel::kFT0Bad}; + parameters::GRPLHCIFData* mLHCIFdata = nullptr; + int mRunNumber = -1; + ctpRateFetcher mRateFetcher; + bool isLeadingBC = false; + + HistogramRegistry mHistManager{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + const int nTriggers = 5; + enum TriggerAliases { kAllBCs = 0, + kFT0Vtx = 1, + kFT0CE = 2, + kFDD = 3, + k1ZNC = 4 }; + const int nBCCategories = 5; + enum BCCategories { kBCA = 0, + kBCB = 1, + kBCC = 2, + kBCE = 3, + kBCL = 4 }; + + static constexpr std::string_view NBCsVsTimeHistNames[5][5] = + {{"AllBCs/BC_A/nBCsVsTime", "AllBCs/BC_B/nBCsVsTime", "AllBCs/BC_C/nBCsVsTime", "AllBCs/BC_E/nBCsVsTime", "AllBCs/BC_L/nBCsVsTime"}, + {"FT0VTx/BC_A/nBCsVsTime", "FT0VTx/BC_B/nBCsVsTime", "FT0VTx/BC_C/nBCsVsTime", "FT0VTx/BC_E/nBCsVsTime", "FT0VTx/BC_L/nBCsVsTime"}, + {"FT0CE/BC_A/nBCsVsTime", "FT0CE/BC_B/nBCsVsTime", "FT0CE/BC_C/nBCsVsTime", "FT0CE/BC_E/nBCsVsTime", "FT0CE/BC_L/nBCsVsTime"}, + {"FDD/BC_A/nBCsVsTime", "FDD/BC_B/nBCsVsTime", "FDD/BC_C/nBCsVsTime", "FDD/BC_E/nBCsVsTime", "FDD/BC_L/nBCsVsTime"}, + {"1ZNC/BC_A/nBCsVsTime", "1ZNC/BC_B/nBCsVsTime", "1ZNC/BC_C/nBCsVsTime", "1ZNC/BC_E/nBCsVsTime", "1ZNC/BC_L/nBCsVsTime"}}; + + static constexpr std::string_view NBCsVsBCIDHistNames[5][5] = + {{"AllBCs/BC_A/nBCsVsBCID", "AllBCs/BC_B/nBCsVsBCID", "AllBCs/BC_C/nBCsVsBCID", "AllBCs/BC_E/nBCsVsBCID", "AllBCs/BC_L/nBCsVsBCID"}, + {"FT0VTx/BC_A/nBCsVsBCID", "FT0VTx/BC_B/nBCsVsBCID", "FT0VTx/BC_C/nBCsVsBCID", "FT0VTx/BC_E/nBCsVsBCID", "FT0VTx/BC_L/nBCsVsBCID"}, + {"FT0CE/BC_A/nBCsVsBCID", "FT0CE/BC_B/nBCsVsBCID", "FT0CE/BC_C/nBCsVsBCID", "FT0CE/BC_E/nBCsVsBCID", "FT0CE/BC_L/nBCsVsBCID"}, + {"FDD/BC_A/nBCsVsBCID", "FDD/BC_B/nBCsVsBCID", "FDD/BC_C/nBCsVsBCID", "FDD/BC_E/nBCsVsBCID", "FDD/BC_L/nBCsVsBCID"}, + {"1ZNC/BC_A/nBCsVsBCID", "1ZNC/BC_B/nBCsVsBCID", "1ZNC/BC_C/nBCsVsBCID", "1ZNC/BC_E/nBCsVsBCID", "1ZNC/BC_L/nBCsVsBCID"}}; + + int64_t bcSOR; + int nBCsPerTF; + int64_t currentTFid = -1; + + void init(InitContext&) + { + mHistManager.add("hMu", "hMu", HistType::kTH1F, {{2000, 0., 0.2}}); + + strLPMProductionTag = metadataInfo.get("LPMProductionTag"); // to extract info from ccdb by the tag + + LOG(info) << "strLPMProductionTag: " << strLPMProductionTag; + + AxisSpec timeAxis{1200, 0., 1200., "#bf{t-t_{SOF} (min)}"}, bcIDAxis{3600, 0., 3600., "#bf{BC ID in orbit}"}; + + for (int iTrigger = 0; iTrigger < nTriggers; iTrigger++) { + if ((iTrigger == kAllBCs) || (iTrigger == kFT0Vtx && cfgDoFT0Vtx) || (iTrigger == kFT0CE && cfgDoFT0CE) || (iTrigger == kFDD && cfgDoFDD) || (iTrigger == k1ZNC && cfgDo1ZNC)) { + for (int iBCCategory = 0; iBCCategory < nBCCategories; iBCCategory++) { + if ((iBCCategory == kBCA && cfgDoBCA) || (iBCCategory == kBCB && cfgDoBCB) || (iBCCategory == kBCC && cfgDoBCC) || (iBCCategory == kBCE && cfgDoBCE) || (iBCCategory == kBCL && cfgDoBCL)) { + mHistManager.add(Form("%s", std::string(NBCsVsTimeHistNames[iTrigger][iBCCategory]).c_str()), "Time of triggered BCs since the start of fill;#bf{t-t_{SOF} (min)};#bf{#it{N}_{BC}}", HistType::kTH1F, {timeAxis}); + mHistManager.add(Form("%s", std::string(NBCsVsBCIDHistNames[iTrigger][iBCCategory]).c_str()), "BC ID of triggered BCs;#bf{BC ID in orbit};#bf{#it{N}_{BC}}", HistType::kTH1F, {bcIDAxis}); + } + } + } + } + + mHistManager.add("FT0Vtx_EvSel/nBCsVsTime", "Time of TVX triggered BCs since the start of fill;;#bf{#it{N}_{BC}}", HistType::kTH1F, {timeAxis}); + mHistManager.add("nBCsVsBCID", "Time of TVX triggered BCs since the start of fill;#bf{t-t_{SOF} (min)};#bf{#it{N}_{BC}}", HistType::kTH1F, {bcIDAxis}); + mHistManager.add("InteractionRateVsTime", "IR from CTP vs time since SOF;#bf{t-t_{SOF} (min)};#bf{#it{N}_{BC}}", HistType::kTH1F, {timeAxis}); + mHistManager.add("TFsPerMinute", "TFs seen in this minute (to account for failed jobs);#bf{t-t_{SOF} (min)};#bf{#it{N}_{TFs}}", HistType::kTH1F, {timeAxis}); + + if (cfgRequireGoodRCTQuality) + isFT0GoodRCTChecker.init({aod::rctsel::kFT0Bad}); + } + + void setLHCIFData(const auto& bc) + { + if (mRunNumber == bc.runNumber()) + return; + + auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance(); + uint64_t timeStamp = bc.timestamp(); + + std::map metadata; + mLHCIFdata = ccdbMgr.getSpecific("GLO/Config/GRPLHCIF", timeStamp, metadata); + if (mLHCIFdata == nullptr) + LOG(fatal) << "GRPLHCIFData not in database, timestamp:" << timeStamp; + mRunNumber = bc.runNumber(); + LOG(info) << "LHCIF data fetched for run " << mRunNumber << " and timestamp " << timeStamp; + + beamPatternA = mLHCIFdata->getBunchFilling().getBeamPattern(0); + beamPatternC = mLHCIFdata->getBunchFilling().getBeamPattern(1); + bcPatternA = beamPatternA & ~beamPatternC; + bcPatternC = ~beamPatternA & beamPatternC; + bcPatternB = beamPatternA & beamPatternC; + bcPatternE = ~beamPatternA & ~beamPatternC; + + auto runInfo = o2::parameters::AggregatedRunInfo::buildAggregatedRunInfo(o2::ccdb::BasicCCDBManager::instance(), mRunNumber, strLPMProductionTag); + bcSOR = runInfo.orbitSOR * nBCsPerOrbit; // first bc of the first orbit + LOG(info) << "BC SOR: " << bcSOR << " (orbit SOR: " << runInfo.orbitSOR << ") NBCs per orbit: " << nBCsPerOrbit; + nBCsPerTF = runInfo.orbitsPerTF * nBCsPerOrbit; // duration of TF in bcs + + return; + } + + double getTVXRate(const auto& bc) + { + auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance(); + double tvxRate = mRateFetcher.fetch(&ccdbMgr, bc.timestamp(), bc.runNumber(), "T0VTX"); + + return tvxRate; + } + + double calculateMu(const auto& bc) + { + + auto bfilling = mLHCIFdata->getBunchFilling(); + double nbc = bfilling.getFilledBCs().size(); + double nTriggersPerFilledBC = getTVXRate(bc) / nbc / o2::constants::lhc::LHCRevFreq; + double mu = -std::log(1 - nTriggersPerFilledBC); + + return mu; + } + + float getTimeSinceSOF(const auto& bc) + { + return (bc.timestamp() - mLHCIFdata->getFillNumberTime()) / 1e3 / 60; // Convert to minutes + } + + template + void fillHistograms(float timeSinceSOF, int64_t localBC) + { + mHistManager.fill(HIST(NBCsVsTimeHistNames[iTrigger][iBCCategory]), timeSinceSOF); + mHistManager.fill(HIST(NBCsVsBCIDHistNames[iTrigger][iBCCategory]), localBC); + } + + void process(MyBCs const& bcs, aod::FT0s const&) + { + int nEmptyBCs = 0; + for (const auto& bc : bcs) { + + if (bc.timestamp() == 0) + continue; + + setLHCIFData(bc); + + mHistManager.fill(HIST("hMu"), calculateMu(bc)); + + float timeSinceSOF = getTimeSinceSOF(bc); + + auto hRateHist = mHistManager.get(HIST("InteractionRateVsTime")); + hRateHist->SetBinContent(hRateHist->FindBin(timeSinceSOF), getTVXRate(bc)); + + if (bc.selection_bit(aod::evsel::kIsTriggerTVX)) + mHistManager.fill(HIST("FT0Vtx_EvSel/nBCsVsTime"), timeSinceSOF); + + int64_t globalBC = bc.globalBC(); + int localBC = globalBC % nBCsPerOrbit; + + int64_t thisTFid = (globalBC - bcSOR) / nBCsPerTF; + + if (thisTFid != currentTFid) { + currentTFid = thisTFid; + mHistManager.fill(HIST("TFsPerMinute"), timeSinceSOF); + } + + if (bcPatternB[localBC] && nEmptyBCs >= cfgEmptyBCsBeforeLeadingBC) { + isLeadingBC = true; + nEmptyBCs = 0; + } else { + isLeadingBC = false; + nEmptyBCs++; + } + + std::bitset<64> ctpInputMask(bc.inputMask()); + + for (int iTrigger = 0; iTrigger < nTriggers; iTrigger++) { + if ((iTrigger == kAllBCs) || (iTrigger == kFT0Vtx && cfgDoFT0Vtx) || (iTrigger == kFT0CE && cfgDoFT0CE) || (iTrigger == kFDD && cfgDoFDD) || (iTrigger == k1ZNC && cfgDo1ZNC)) { + for (int iBCCategory = 0; iBCCategory < nBCCategories; iBCCategory++) { + if ((iBCCategory == kBCA && cfgDoBCA) || (iBCCategory == kBCB && cfgDoBCB) || (iBCCategory == kBCC && cfgDoBCC) || (iBCCategory == kBCE && cfgDoBCE) || (iBCCategory == kBCL && cfgDoBCL)) { + if (iTrigger == kAllBCs) { + if (iBCCategory == kBCA && bcPatternA[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCB && bcPatternB[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCC && bcPatternC[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCE && bcPatternE[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCL && isLeadingBC) + fillHistograms(timeSinceSOF, localBC); + } + if (iTrigger == kFT0Vtx && ctpInputMask.test(2)) { + if (iBCCategory == kBCA && bcPatternA[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCB && bcPatternB[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCC && bcPatternC[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCE && bcPatternE[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCL && isLeadingBC) + fillHistograms(timeSinceSOF, localBC); + } + if (iTrigger == kFT0CE && ctpInputMask.test(4)) { + if (iBCCategory == kBCA && bcPatternA[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCB && bcPatternB[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCC && bcPatternC[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCE && bcPatternE[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCL && isLeadingBC) + fillHistograms(timeSinceSOF, localBC); + } + if (iTrigger == kFDD && ctpInputMask.test(15)) { + if (iBCCategory == kBCA && bcPatternA[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCB && bcPatternB[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCC && bcPatternC[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCE && bcPatternE[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCL && isLeadingBC) + fillHistograms(timeSinceSOF, localBC); + } + if (iTrigger == k1ZNC && ctpInputMask.test(25)) { + if (iBCCategory == kBCA && bcPatternA[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCB && bcPatternB[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCC && bcPatternC[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCE && bcPatternE[localBC]) + fillHistograms(timeSinceSOF, localBC); + if (iBCCategory == kBCL && isLeadingBC) + fillHistograms(timeSinceSOF, localBC); + } + } + } + } + } + mHistManager.fill(HIST("nBCsVsBCID"), localBC); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + metadataInfo.initMetadata(cfgc); + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}