From 5e2f410cabba2a530cecc2420301ab15903888c8 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 24 Mar 2026 19:28:59 +0000 Subject: [PATCH] Please consider the following formatting changes --- PWGUD/Core/UDHelpers.h | 2 +- PWGUD/TableProducer/SGCandProducer.cxx | 1584 ++++++++++++------------ 2 files changed, 793 insertions(+), 793 deletions(-) diff --git a/PWGUD/Core/UDHelpers.h b/PWGUD/Core/UDHelpers.h index ea137fbc17a..9dbd70c6897 100644 --- a/PWGUD/Core/UDHelpers.h +++ b/PWGUD/Core/UDHelpers.h @@ -677,7 +677,7 @@ inline bool getPhiEtaFromFitBit(FT0DetT& ft0Det, double& phi, double& eta) { - auto ref = decodeFitBit(bit); + auto ref = decodeFitBit(bit); if (ref.det != FitBitRef::Det::FT0) { return false; } diff --git a/PWGUD/TableProducer/SGCandProducer.cxx b/PWGUD/TableProducer/SGCandProducer.cxx index 80ec14a0b1f..da3c052c036 100644 --- a/PWGUD/TableProducer/SGCandProducer.cxx +++ b/PWGUD/TableProducer/SGCandProducer.cxx @@ -57,803 +57,803 @@ using namespace o2::aod::rctsel; #define getHist(type, name) std::get>(histPointers[name]) struct SGCandProducer { - Service ccdb; - // data inputs - using CCs = soa::Join; - using CC = CCs::iterator; - using BCs = soa::Join; - using BC = BCs::iterator; - using TCs = soa::Join; - using FWs = aod::FwdTracks; - - using MCCCs = soa::Join; - using MCCC = MCCCs::iterator; - - // get an SGCutparHolder - SGCutParHolder sameCuts = SGCutParHolder(); // SGCutparHolder - Configurable SGCuts{"SGCuts", {}, "SG event cuts"}; - FITCutParHolder fitCuts = FITCutParHolder(); - Configurable FITCuts{"FITCuts", {}, "FIT bitset cuts"}; - Configurable verboseInfo{"verboseInfo", false, "Print general info to terminal; default it false."}; - Configurable saveAllTracks{"saveAllTracks", true, "save only PV contributors or all tracks associated to a collision"}; - Configurable savenonPVCITSOnlyTracks{"savenonPVCITSOnlyTracks", false, "save non PV contributors with ITS only information"}; - Configurable rejectAtTFBoundary{"rejectAtTFBoundary", true, "reject collisions at a TF boundary"}; - Configurable noITSROFrameBorder{"noITSROFrameBorder", true, "reject ITS RO Frame Border"}; - Configurable noSameBunchPileUp{"noSameBunchPileUp", true, "reject SameBunchPileUp"}; - Configurable IsGoodVertex{"IsGoodVertex", false, "Select FT0 PV vertex matching"}; - Configurable ITSTPCVertex{"ITSTPCVertex", true, "reject ITS-only vertex"}; // if one wants to look at Single Gap pp events - Configurable> generatorIds{"generatorIds", std::vector{-1}, "MC generatorIds to process"}; - Configurable storeSG{"storeSG", true, "Store SG events in the output"}; - Configurable storeDG{"storeDG", true, "Store DG events in the output"}; - - Configurable saveFITbitsets{"saveFITbitsets", true, "Write FT0 and FV0 bitset tables to output"}; - - Configurable isGoodRCTCollision{"isGoodRCTCollision", true, "Check RCT flags for FT0,ITS,TPC and tracking"}; - Configurable isGoodRCTZdc{"isGoodRCTZdc", false, "Check RCT flags for ZDC if present in run"}; - - // Configurables to decide which tables are filled - Configurable fillTrackTables{"fillTrackTables", true, "Fill track tables"}; - Configurable fillFwdTrackTables{"fillFwdTrackTables", true, "Fill forward track tables"}; - - // SG selector - SGSelector sgSelector; - ctpRateFetcher mRateFetcher; - - // initialize RCT flag checker - RCTFlagsChecker myRCTChecker{"CBT"}; - - // data tables - Produces outputSGCollisions; - Produces outputCollisions; - Produces outputCollisionsSels; - Produces outputCollisionSelExtras; - Produces outputCollsLabels; - Produces outputZdcs; - Produces udZdcsReduced; - Produces outputTracks; - Produces outputTracksCov; - Produces outputTracksDCA; - Produces outputTracksPID; - Produces outputTracksPIDExtra; - Produces outputTracksExtra; - Produces outputTracksFlag; - Produces outputFwdTracks; - Produces outputFwdTracksExtra; - Produces outputTracksLabel; - Produces outputFITBits; - - // initialize histogram registry - HistogramRegistry registry{ - "registry", - {}}; - std::map histPointers; - - int runNumber = -1; - - // function to update UDFwdTracks, UDFwdTracksExtra - template - void updateUDFwdTrackTables(TFwdTrack const& fwdtrack, uint64_t const& bcnum) - { - outputFwdTracks(outputCollisions.lastIndex(), - fwdtrack.px(), fwdtrack.py(), fwdtrack.pz(), fwdtrack.sign(), - bcnum, fwdtrack.trackTime(), fwdtrack.trackTimeRes()); - outputFwdTracksExtra(fwdtrack.trackType(), - fwdtrack.nClusters(), - fwdtrack.pDca(), - fwdtrack.rAtAbsorberEnd(), - fwdtrack.chi2(), - fwdtrack.chi2MatchMCHMID(), - fwdtrack.chi2MatchMCHMFT(), - fwdtrack.mchBitMap(), - fwdtrack.midBitMap(), - fwdtrack.midBoards()); - } - - // function to update UDTracks, UDTracksCov, UDTracksDCA, UDTracksPID, UDTracksExtra, UDTracksFlag, - // and UDTrackCollisionIDs - template - void updateUDTrackTables(int64_t lastIndex, TTrack const& track, uint64_t const& bcnum) - { - outputTracks(lastIndex, - track.px(), track.py(), track.pz(), track.sign(), - bcnum, track.trackTime(), track.trackTimeRes()); - - // float sigmaY = track.sigmaY(); - // float sigmaZ = track.sigmaZ(); - float sigmaY = -1.; - float sigmaZ = -1.; - outputTracksCov(track.x(), track.y(), track.z(), sigmaY, sigmaZ); - - outputTracksDCA(track.dcaZ(), track.dcaXY()); - outputTracksPID(track.tpcNSigmaEl(), - track.tpcNSigmaMu(), - track.tpcNSigmaPi(), - track.tpcNSigmaKa(), - track.tpcNSigmaPr(), - track.beta(), - track.betaerror(), - track.tofNSigmaEl(), - track.tofNSigmaMu(), - track.tofNSigmaPi(), - track.tofNSigmaKa(), - track.tofNSigmaPr()); - outputTracksPIDExtra(track.tpcNSigmaDe(), - track.tpcNSigmaTr(), - track.tpcNSigmaHe(), - track.tpcNSigmaAl(), - track.tofNSigmaDe(), - track.tofNSigmaTr(), - track.tofNSigmaHe(), - track.tofNSigmaAl()); - outputTracksExtra(track.tpcInnerParam(), - track.itsClusterSizes(), - track.tpcNClsFindable(), - track.tpcNClsFindableMinusFound(), - track.tpcNClsFindableMinusCrossedRows(), - track.tpcNClsShared(), - track.trdPattern(), - track.itsChi2NCl(), - track.tpcChi2NCl(), - track.trdChi2(), - track.tofChi2(), - track.tpcSignal(), - track.tofSignal(), - track.trdSignal(), - track.length(), - track.tofExpMom(), - track.detectorMap()); - outputTracksFlag(track.has_collision(), - track.isPVContributor()); - outputTracksLabel(track.globalIndex()); - } - - // function to process reconstructed data - template - void processReco(std::string histdir, TCol const& collision, BCs const& bcs, - TCs const& tracks, FWs const& fwdtracks, - aod::FV0As const& fv0as, aod::FT0s const& ft0s, aod::FDDs const& fdds) - { - if (verboseInfo) - LOGF(debug, " collision %d", collision.globalIndex()); - getHist(TH1, histdir + "/Stat")->Fill(0., 1.); - // reject collisions at TF boundaries - if (rejectAtTFBoundary && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(1., 1.); - // reject collisions at ITS RO TF boundaries - if (noITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(2., 1.); - // reject Same Bunch PileUp - if (noSameBunchPileUp && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(3., 1.); - // check vertex matching to FT0 - if (IsGoodVertex && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(4., 1.); - // reject ITS Only vertices - if (ITSTPCVertex && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(5., 1.); - // nominal BC - if (!collision.has_foundBC()) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(6., 1.); - // RCT CBT for collision check - if (isGoodRCTCollision && !myRCTChecker(collision)) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(7., 1.); - // RCT CBT+ZDC for collision check - if (isGoodRCTZdc && !myRCTChecker(collision)) { - return; - } - getHist(TH1, histdir + "/Stat")->Fill(8., 1.); - - // - const int trs = collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard) ? 1 : 0; - const int trofs = collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard) ? 1 : 0; - const int hmpr = collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof) ? 1 : 0; - const int tfb = collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) ? 1 : 0; - const int itsROFb = collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder) ? 1 : 0; - const int sbp = collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup) ? 1 : 0; - const int zVtxFT0vPv = collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV) ? 1 : 0; - const int vtxITSTPC = collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC) ? 1 : 0; - auto bc = collision.template foundBC_as(); - double ir = 0.; - const uint64_t ts = bc.timestamp(); - const int runnumber = bc.runNumber(); - if (bc.has_zdc()) { - ir = mRateFetcher.fetch(ccdb.service, ts, runnumber, "ZNC hadronic") * 1.e-3; - } - auto newbc = bc; - - // obtain slice of compatible BCs - auto bcRange = udhelpers::compatibleBCs(collision, sameCuts.NDtcoll(), bcs, sameCuts.minNBCs()); - auto isSGEvent = sgSelector.IsSelected(sameCuts, collision, bcRange, bc); - // auto isSGEvent = sgSelector.IsSelected(sameCuts, collision, bcRange, tracks); - int issgevent = isSGEvent.value; - if (isSGEvent.bc && issgevent < 2) { - newbc = *(isSGEvent.bc); - } else { - if (verboseInfo) - LOGF(info, "No Newbc %i", bc.globalBC()); - } - getHist(TH1, histdir + "/Stat")->Fill(issgevent + 10, 1.); - if ((storeDG && issgevent == o2::aod::sgselector::DoubleGap) || (storeSG && (issgevent == o2::aod::sgselector::SingleGapA || issgevent == o2::aod::sgselector::SingleGapC))) { - if (verboseInfo) - LOGF(info, "Current BC: %i, %i, %i", bc.globalBC(), newbc.globalBC(), issgevent); - if (sameCuts.minRgtrwTOF()) { - if (udhelpers::rPVtrwTOF(tracks, collision.numContrib()) < sameCuts.minRgtrwTOF()) - return; - } - upchelpers::FITInfo fitInfo{}; - const uint8_t chFT0A = 0; - const uint8_t chFT0C = 0; - const uint8_t chFDDA = 0; - const uint8_t chFDDC = 0; - const uint8_t chFV0A = 0; - const int occ = collision.trackOccupancyInTimeRange(); - udhelpers::getFITinfo(fitInfo, newbc, bcs, ft0s, fv0as, fdds); - const int upc_flag = (collision.flags() & dataformats::Vertex>::Flags::UPCMode) ? 1 : 0; - // update SG candidates tables - outputCollisions(bc.globalBC(), bc.runNumber(), - collision.posX(), collision.posY(), collision.posZ(), upc_flag, - collision.numContrib(), udhelpers::netCharge(tracks), - 1.); // rtrwTOF); //omit the calculation to speed up the things while skimming - - outputSGCollisions(issgevent); - outputCollisionsSels(fitInfo.ampFT0A, fitInfo.ampFT0C, fitInfo.timeFT0A, fitInfo.timeFT0C, - fitInfo.triggerMaskFT0, - fitInfo.ampFDDA, fitInfo.ampFDDC, fitInfo.timeFDDA, fitInfo.timeFDDC, - fitInfo.triggerMaskFDD, - fitInfo.ampFV0A, fitInfo.timeFV0A, fitInfo.triggerMaskFV0A, - fitInfo.BBFT0Apf, fitInfo.BBFT0Cpf, fitInfo.BGFT0Apf, fitInfo.BGFT0Cpf, - fitInfo.BBFV0Apf, fitInfo.BGFV0Apf, - fitInfo.BBFDDApf, fitInfo.BBFDDCpf, fitInfo.BGFDDApf, fitInfo.BGFDDCpf); - outputCollisionSelExtras(chFT0A, chFT0C, chFDDA, chFDDC, chFV0A, occ, ir, trs, trofs, hmpr, tfb, itsROFb, sbp, zVtxFT0vPv, vtxITSTPC, collision.rct_raw()); - outputCollsLabels(collision.globalIndex()); - - uint64_t w1[4] = {0ull, 0ull, 0ull, 0ull}; - uint64_t w2[4] = {0ull, 0ull, 0ull, 0ull}; - - if (fitCuts.saveFITbitsets() && newbc.has_foundFT0() && newbc.has_fv0a()) { - udhelpers::buildFT0FV0Words(newbc.ft0(), newbc.fv0a(), w1, w2, - fitCuts.thr1_FT0A(), - fitCuts.thr1_FT0C(), - fitCuts.thr1_FV0A(), - fitCuts.thr2_FT0A(), - fitCuts.thr2_FT0C(), - fitCuts.thr2_FV0A()); - } - - outputFITBits(w1[0], w1[1], w1[2], w1[3], - w2[0], w2[1], w2[2], w2[3]); - - if (newbc.has_zdc()) { - auto zdc = newbc.zdc(); - udZdcsReduced(outputCollisions.lastIndex(), zdc.timeZNA(), zdc.timeZNC(), zdc.energyCommonZNA(), zdc.energyCommonZNC()); - } else { - udZdcsReduced(outputCollisions.lastIndex(), -999, -999, -999, -999); - } - // update SGTracks tables - if (fillTrackTables) { - for (const auto& track : tracks) { - if (track.pt() > sameCuts.minPt() && track.eta() > sameCuts.minEta() && track.eta() < sameCuts.maxEta()) { - if (track.isPVContributor()) { - updateUDTrackTables(outputCollisions.lastIndex(), track, bc.globalBC()); - } else if (saveAllTracks) { - if (track.itsClusterSizes() && track.itsChi2NCl() > 0 && ((track.tpcNClsFindable() == 0 && savenonPVCITSOnlyTracks) || track.tpcNClsFindable() > 50)) - updateUDTrackTables(outputCollisions.lastIndex(), track, bc.globalBC()); - // if (track.isPVContributor()) updateUDTrackTables(outputCollisions.lastIndex(), track, bc.globalBC()); - } - } - } - } - // update SGFwdTracks tables - if (fillFwdTrackTables) { - if (sameCuts.withFwdTracks()) { - for (const auto& fwdtrack : fwdtracks) { - if (!sgSelector.FwdTrkSelector(fwdtrack)) - updateUDFwdTrackTables(fwdtrack, bc.globalBC()); - } - } - } - } - } - - void init(InitContext& context) - { - ccdb->setURL("http://alice-ccdb.cern.ch"); - ccdb->setCaching(true); - ccdb->setFatalWhenNull(false); - sameCuts = (SGCutParHolder)SGCuts; - fitCuts = (FITCutParHolder)FITCuts; - - // add histograms for the different process functions - histPointers.clear(); - if (context.mOptions.get("processData")) { - histPointers.insert({"reco/Stat", registry.add("reco/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{14, -0.5, 13.5}}})}); - - const AxisSpec axisCountersTrg{10, 0.5, 10.5, ""}; - histPointers.insert({"reco/hCountersTrg", registry.add("reco/hCountersTrg", "Trigger counts before selections; Trigger; Counts", {HistType::kTH1F, {axisCountersTrg}})}); - histPointers.insert({"reco/hCountersTrgBcSel", registry.add("reco/hCountersTrgSel", "Trigger counts after BC selections; Trigger; Counts", {HistType::kTH1F, {axisCountersTrg}})}); - histPointers.insert({"reco/hLumi", registry.add("reco/hLumi", "Integrated luminosity before selections; Trigger; Luminosity, 1/#mub", {HistType::kTH1F, {axisCountersTrg}})}); - histPointers.insert({"reco/hLumiBcSel", registry.add("reco/hLumiBcSel", "Integrated luminosity before selections; Trigger; Luminosity, 1/#mub", {HistType::kTH1F, {axisCountersTrg}})}); - auto hCountersTrg = getHist(TH1, "reco/hCountersTrg"); - auto hCountersTrgBcSel = getHist(TH1, "reco/hCountersTrgBcSel"); - auto hLumi = getHist(TH1, "reco/hLumi"); - auto hLumiBcSel = getHist(TH1, "reco/hLumiBcSel"); - for (const auto& h : {hCountersTrg, hCountersTrgBcSel, hLumi, hLumiBcSel}) { - h->GetXaxis()->SetBinLabel(1, "TVX"); - h->GetXaxis()->SetBinLabel(2, "TCE"); - h->GetXaxis()->SetBinLabel(3, "ZEM"); - h->GetXaxis()->SetBinLabel(4, "ZNC"); - } - } - if (context.mOptions.get("processMcData")) { - histPointers.insert({"MCreco/Stat", registry.add("MCreco/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{14, -0.5, 13.5}}})}); - } - - if (isGoodRCTZdc) { - myRCTChecker.init("CBT", true); - } - } - - // process function for reconstructed data - void processData(CC const& collision, BCs const& bcs, TCs const& tracks, FWs const& fwdtracks, - aod::Zdcs const& /*zdcs*/, aod::FV0As const& fv0as, aod::FT0s const& ft0s, aod::FDDs const& fdds) - { - processReco(std::string("reco"), collision, bcs, tracks, fwdtracks, fv0as, ft0s, fdds); - } - PROCESS_SWITCH(SGCandProducer, processData, "Produce UD table with data", true); - - // process function for reconstructed MC data - void processMcData(MCCC const& collision, aod::McCollisions const& /*mccollisions*/, BCs const& bcs, - TCs const& tracks, FWs const& fwdtracks, aod::Zdcs const& /*zdcs*/, aod::FV0As const& fv0as, - aod::FT0s const& ft0s, aod::FDDs const& fdds) - { - // select specific processes with the GeneratorID - if (!collision.has_mcCollision()) - return; - auto mccol = collision.mcCollision(); - if (verboseInfo) - LOGF(info, "GeneratorId %d (%d)", mccol.getGeneratorId(), generatorIds->size()); - - if (std::find(generatorIds->begin(), generatorIds->end(), mccol.getGeneratorId()) != generatorIds->end()) { - if (verboseInfo) - LOGF(info, "Event with good generatorId"); - processReco(std::string("MCreco"), collision, bcs, tracks, fwdtracks, fv0as, ft0s, fdds); - } - } - PROCESS_SWITCH(SGCandProducer, processMcData, "Produce UD tables with MC data", false); + Service ccdb; + // data inputs + using CCs = soa::Join; + using CC = CCs::iterator; + using BCs = soa::Join; + using BC = BCs::iterator; + using TCs = soa::Join; + using FWs = aod::FwdTracks; + + using MCCCs = soa::Join; + using MCCC = MCCCs::iterator; + + // get an SGCutparHolder + SGCutParHolder sameCuts = SGCutParHolder(); // SGCutparHolder + Configurable SGCuts{"SGCuts", {}, "SG event cuts"}; + FITCutParHolder fitCuts = FITCutParHolder(); + Configurable FITCuts{"FITCuts", {}, "FIT bitset cuts"}; + Configurable verboseInfo{"verboseInfo", false, "Print general info to terminal; default it false."}; + Configurable saveAllTracks{"saveAllTracks", true, "save only PV contributors or all tracks associated to a collision"}; + Configurable savenonPVCITSOnlyTracks{"savenonPVCITSOnlyTracks", false, "save non PV contributors with ITS only information"}; + Configurable rejectAtTFBoundary{"rejectAtTFBoundary", true, "reject collisions at a TF boundary"}; + Configurable noITSROFrameBorder{"noITSROFrameBorder", true, "reject ITS RO Frame Border"}; + Configurable noSameBunchPileUp{"noSameBunchPileUp", true, "reject SameBunchPileUp"}; + Configurable IsGoodVertex{"IsGoodVertex", false, "Select FT0 PV vertex matching"}; + Configurable ITSTPCVertex{"ITSTPCVertex", true, "reject ITS-only vertex"}; // if one wants to look at Single Gap pp events + Configurable> generatorIds{"generatorIds", std::vector{-1}, "MC generatorIds to process"}; + Configurable storeSG{"storeSG", true, "Store SG events in the output"}; + Configurable storeDG{"storeDG", true, "Store DG events in the output"}; + + Configurable saveFITbitsets{"saveFITbitsets", true, "Write FT0 and FV0 bitset tables to output"}; + + Configurable isGoodRCTCollision{"isGoodRCTCollision", true, "Check RCT flags for FT0,ITS,TPC and tracking"}; + Configurable isGoodRCTZdc{"isGoodRCTZdc", false, "Check RCT flags for ZDC if present in run"}; + + // Configurables to decide which tables are filled + Configurable fillTrackTables{"fillTrackTables", true, "Fill track tables"}; + Configurable fillFwdTrackTables{"fillFwdTrackTables", true, "Fill forward track tables"}; + + // SG selector + SGSelector sgSelector; + ctpRateFetcher mRateFetcher; + + // initialize RCT flag checker + RCTFlagsChecker myRCTChecker{"CBT"}; + + // data tables + Produces outputSGCollisions; + Produces outputCollisions; + Produces outputCollisionsSels; + Produces outputCollisionSelExtras; + Produces outputCollsLabels; + Produces outputZdcs; + Produces udZdcsReduced; + Produces outputTracks; + Produces outputTracksCov; + Produces outputTracksDCA; + Produces outputTracksPID; + Produces outputTracksPIDExtra; + Produces outputTracksExtra; + Produces outputTracksFlag; + Produces outputFwdTracks; + Produces outputFwdTracksExtra; + Produces outputTracksLabel; + Produces outputFITBits; + + // initialize histogram registry + HistogramRegistry registry{ + "registry", + {}}; + std::map histPointers; + + int runNumber = -1; + + // function to update UDFwdTracks, UDFwdTracksExtra + template + void updateUDFwdTrackTables(TFwdTrack const& fwdtrack, uint64_t const& bcnum) + { + outputFwdTracks(outputCollisions.lastIndex(), + fwdtrack.px(), fwdtrack.py(), fwdtrack.pz(), fwdtrack.sign(), + bcnum, fwdtrack.trackTime(), fwdtrack.trackTimeRes()); + outputFwdTracksExtra(fwdtrack.trackType(), + fwdtrack.nClusters(), + fwdtrack.pDca(), + fwdtrack.rAtAbsorberEnd(), + fwdtrack.chi2(), + fwdtrack.chi2MatchMCHMID(), + fwdtrack.chi2MatchMCHMFT(), + fwdtrack.mchBitMap(), + fwdtrack.midBitMap(), + fwdtrack.midBoards()); + } + + // function to update UDTracks, UDTracksCov, UDTracksDCA, UDTracksPID, UDTracksExtra, UDTracksFlag, + // and UDTrackCollisionIDs + template + void updateUDTrackTables(int64_t lastIndex, TTrack const& track, uint64_t const& bcnum) + { + outputTracks(lastIndex, + track.px(), track.py(), track.pz(), track.sign(), + bcnum, track.trackTime(), track.trackTimeRes()); + + // float sigmaY = track.sigmaY(); + // float sigmaZ = track.sigmaZ(); + float sigmaY = -1.; + float sigmaZ = -1.; + outputTracksCov(track.x(), track.y(), track.z(), sigmaY, sigmaZ); + + outputTracksDCA(track.dcaZ(), track.dcaXY()); + outputTracksPID(track.tpcNSigmaEl(), + track.tpcNSigmaMu(), + track.tpcNSigmaPi(), + track.tpcNSigmaKa(), + track.tpcNSigmaPr(), + track.beta(), + track.betaerror(), + track.tofNSigmaEl(), + track.tofNSigmaMu(), + track.tofNSigmaPi(), + track.tofNSigmaKa(), + track.tofNSigmaPr()); + outputTracksPIDExtra(track.tpcNSigmaDe(), + track.tpcNSigmaTr(), + track.tpcNSigmaHe(), + track.tpcNSigmaAl(), + track.tofNSigmaDe(), + track.tofNSigmaTr(), + track.tofNSigmaHe(), + track.tofNSigmaAl()); + outputTracksExtra(track.tpcInnerParam(), + track.itsClusterSizes(), + track.tpcNClsFindable(), + track.tpcNClsFindableMinusFound(), + track.tpcNClsFindableMinusCrossedRows(), + track.tpcNClsShared(), + track.trdPattern(), + track.itsChi2NCl(), + track.tpcChi2NCl(), + track.trdChi2(), + track.tofChi2(), + track.tpcSignal(), + track.tofSignal(), + track.trdSignal(), + track.length(), + track.tofExpMom(), + track.detectorMap()); + outputTracksFlag(track.has_collision(), + track.isPVContributor()); + outputTracksLabel(track.globalIndex()); + } + + // function to process reconstructed data + template + void processReco(std::string histdir, TCol const& collision, BCs const& bcs, + TCs const& tracks, FWs const& fwdtracks, + aod::FV0As const& fv0as, aod::FT0s const& ft0s, aod::FDDs const& fdds) + { + if (verboseInfo) + LOGF(debug, " collision %d", collision.globalIndex()); + getHist(TH1, histdir + "/Stat")->Fill(0., 1.); + // reject collisions at TF boundaries + if (rejectAtTFBoundary && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(1., 1.); + // reject collisions at ITS RO TF boundaries + if (noITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(2., 1.); + // reject Same Bunch PileUp + if (noSameBunchPileUp && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(3., 1.); + // check vertex matching to FT0 + if (IsGoodVertex && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(4., 1.); + // reject ITS Only vertices + if (ITSTPCVertex && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(5., 1.); + // nominal BC + if (!collision.has_foundBC()) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(6., 1.); + // RCT CBT for collision check + if (isGoodRCTCollision && !myRCTChecker(collision)) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(7., 1.); + // RCT CBT+ZDC for collision check + if (isGoodRCTZdc && !myRCTChecker(collision)) { + return; + } + getHist(TH1, histdir + "/Stat")->Fill(8., 1.); + + // + const int trs = collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard) ? 1 : 0; + const int trofs = collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard) ? 1 : 0; + const int hmpr = collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof) ? 1 : 0; + const int tfb = collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) ? 1 : 0; + const int itsROFb = collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder) ? 1 : 0; + const int sbp = collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup) ? 1 : 0; + const int zVtxFT0vPv = collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV) ? 1 : 0; + const int vtxITSTPC = collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC) ? 1 : 0; + auto bc = collision.template foundBC_as(); + double ir = 0.; + const uint64_t ts = bc.timestamp(); + const int runnumber = bc.runNumber(); + if (bc.has_zdc()) { + ir = mRateFetcher.fetch(ccdb.service, ts, runnumber, "ZNC hadronic") * 1.e-3; + } + auto newbc = bc; + + // obtain slice of compatible BCs + auto bcRange = udhelpers::compatibleBCs(collision, sameCuts.NDtcoll(), bcs, sameCuts.minNBCs()); + auto isSGEvent = sgSelector.IsSelected(sameCuts, collision, bcRange, bc); + // auto isSGEvent = sgSelector.IsSelected(sameCuts, collision, bcRange, tracks); + int issgevent = isSGEvent.value; + if (isSGEvent.bc && issgevent < 2) { + newbc = *(isSGEvent.bc); + } else { + if (verboseInfo) + LOGF(info, "No Newbc %i", bc.globalBC()); + } + getHist(TH1, histdir + "/Stat")->Fill(issgevent + 10, 1.); + if ((storeDG && issgevent == o2::aod::sgselector::DoubleGap) || (storeSG && (issgevent == o2::aod::sgselector::SingleGapA || issgevent == o2::aod::sgselector::SingleGapC))) { + if (verboseInfo) + LOGF(info, "Current BC: %i, %i, %i", bc.globalBC(), newbc.globalBC(), issgevent); + if (sameCuts.minRgtrwTOF()) { + if (udhelpers::rPVtrwTOF(tracks, collision.numContrib()) < sameCuts.minRgtrwTOF()) + return; + } + upchelpers::FITInfo fitInfo{}; + const uint8_t chFT0A = 0; + const uint8_t chFT0C = 0; + const uint8_t chFDDA = 0; + const uint8_t chFDDC = 0; + const uint8_t chFV0A = 0; + const int occ = collision.trackOccupancyInTimeRange(); + udhelpers::getFITinfo(fitInfo, newbc, bcs, ft0s, fv0as, fdds); + const int upc_flag = (collision.flags() & dataformats::Vertex>::Flags::UPCMode) ? 1 : 0; + // update SG candidates tables + outputCollisions(bc.globalBC(), bc.runNumber(), + collision.posX(), collision.posY(), collision.posZ(), upc_flag, + collision.numContrib(), udhelpers::netCharge(tracks), + 1.); // rtrwTOF); //omit the calculation to speed up the things while skimming + + outputSGCollisions(issgevent); + outputCollisionsSels(fitInfo.ampFT0A, fitInfo.ampFT0C, fitInfo.timeFT0A, fitInfo.timeFT0C, + fitInfo.triggerMaskFT0, + fitInfo.ampFDDA, fitInfo.ampFDDC, fitInfo.timeFDDA, fitInfo.timeFDDC, + fitInfo.triggerMaskFDD, + fitInfo.ampFV0A, fitInfo.timeFV0A, fitInfo.triggerMaskFV0A, + fitInfo.BBFT0Apf, fitInfo.BBFT0Cpf, fitInfo.BGFT0Apf, fitInfo.BGFT0Cpf, + fitInfo.BBFV0Apf, fitInfo.BGFV0Apf, + fitInfo.BBFDDApf, fitInfo.BBFDDCpf, fitInfo.BGFDDApf, fitInfo.BGFDDCpf); + outputCollisionSelExtras(chFT0A, chFT0C, chFDDA, chFDDC, chFV0A, occ, ir, trs, trofs, hmpr, tfb, itsROFb, sbp, zVtxFT0vPv, vtxITSTPC, collision.rct_raw()); + outputCollsLabels(collision.globalIndex()); + + uint64_t w1[4] = {0ull, 0ull, 0ull, 0ull}; + uint64_t w2[4] = {0ull, 0ull, 0ull, 0ull}; + + if (fitCuts.saveFITbitsets() && newbc.has_foundFT0() && newbc.has_fv0a()) { + udhelpers::buildFT0FV0Words(newbc.ft0(), newbc.fv0a(), w1, w2, + fitCuts.thr1_FT0A(), + fitCuts.thr1_FT0C(), + fitCuts.thr1_FV0A(), + fitCuts.thr2_FT0A(), + fitCuts.thr2_FT0C(), + fitCuts.thr2_FV0A()); + } + + outputFITBits(w1[0], w1[1], w1[2], w1[3], + w2[0], w2[1], w2[2], w2[3]); + + if (newbc.has_zdc()) { + auto zdc = newbc.zdc(); + udZdcsReduced(outputCollisions.lastIndex(), zdc.timeZNA(), zdc.timeZNC(), zdc.energyCommonZNA(), zdc.energyCommonZNC()); + } else { + udZdcsReduced(outputCollisions.lastIndex(), -999, -999, -999, -999); + } + // update SGTracks tables + if (fillTrackTables) { + for (const auto& track : tracks) { + if (track.pt() > sameCuts.minPt() && track.eta() > sameCuts.minEta() && track.eta() < sameCuts.maxEta()) { + if (track.isPVContributor()) { + updateUDTrackTables(outputCollisions.lastIndex(), track, bc.globalBC()); + } else if (saveAllTracks) { + if (track.itsClusterSizes() && track.itsChi2NCl() > 0 && ((track.tpcNClsFindable() == 0 && savenonPVCITSOnlyTracks) || track.tpcNClsFindable() > 50)) + updateUDTrackTables(outputCollisions.lastIndex(), track, bc.globalBC()); + // if (track.isPVContributor()) updateUDTrackTables(outputCollisions.lastIndex(), track, bc.globalBC()); + } + } + } + } + // update SGFwdTracks tables + if (fillFwdTrackTables) { + if (sameCuts.withFwdTracks()) { + for (const auto& fwdtrack : fwdtracks) { + if (!sgSelector.FwdTrkSelector(fwdtrack)) + updateUDFwdTrackTables(fwdtrack, bc.globalBC()); + } + } + } + } + } + + void init(InitContext& context) + { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setFatalWhenNull(false); + sameCuts = (SGCutParHolder)SGCuts; + fitCuts = (FITCutParHolder)FITCuts; + + // add histograms for the different process functions + histPointers.clear(); + if (context.mOptions.get("processData")) { + histPointers.insert({"reco/Stat", registry.add("reco/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{14, -0.5, 13.5}}})}); + + const AxisSpec axisCountersTrg{10, 0.5, 10.5, ""}; + histPointers.insert({"reco/hCountersTrg", registry.add("reco/hCountersTrg", "Trigger counts before selections; Trigger; Counts", {HistType::kTH1F, {axisCountersTrg}})}); + histPointers.insert({"reco/hCountersTrgBcSel", registry.add("reco/hCountersTrgSel", "Trigger counts after BC selections; Trigger; Counts", {HistType::kTH1F, {axisCountersTrg}})}); + histPointers.insert({"reco/hLumi", registry.add("reco/hLumi", "Integrated luminosity before selections; Trigger; Luminosity, 1/#mub", {HistType::kTH1F, {axisCountersTrg}})}); + histPointers.insert({"reco/hLumiBcSel", registry.add("reco/hLumiBcSel", "Integrated luminosity before selections; Trigger; Luminosity, 1/#mub", {HistType::kTH1F, {axisCountersTrg}})}); + auto hCountersTrg = getHist(TH1, "reco/hCountersTrg"); + auto hCountersTrgBcSel = getHist(TH1, "reco/hCountersTrgBcSel"); + auto hLumi = getHist(TH1, "reco/hLumi"); + auto hLumiBcSel = getHist(TH1, "reco/hLumiBcSel"); + for (const auto& h : {hCountersTrg, hCountersTrgBcSel, hLumi, hLumiBcSel}) { + h->GetXaxis()->SetBinLabel(1, "TVX"); + h->GetXaxis()->SetBinLabel(2, "TCE"); + h->GetXaxis()->SetBinLabel(3, "ZEM"); + h->GetXaxis()->SetBinLabel(4, "ZNC"); + } + } + if (context.mOptions.get("processMcData")) { + histPointers.insert({"MCreco/Stat", registry.add("MCreco/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{14, -0.5, 13.5}}})}); + } + + if (isGoodRCTZdc) { + myRCTChecker.init("CBT", true); + } + } + + // process function for reconstructed data + void processData(CC const& collision, BCs const& bcs, TCs const& tracks, FWs const& fwdtracks, + aod::Zdcs const& /*zdcs*/, aod::FV0As const& fv0as, aod::FT0s const& ft0s, aod::FDDs const& fdds) + { + processReco(std::string("reco"), collision, bcs, tracks, fwdtracks, fv0as, ft0s, fdds); + } + PROCESS_SWITCH(SGCandProducer, processData, "Produce UD table with data", true); + + // process function for reconstructed MC data + void processMcData(MCCC const& collision, aod::McCollisions const& /*mccollisions*/, BCs const& bcs, + TCs const& tracks, FWs const& fwdtracks, aod::Zdcs const& /*zdcs*/, aod::FV0As const& fv0as, + aod::FT0s const& ft0s, aod::FDDs const& fdds) + { + // select specific processes with the GeneratorID + if (!collision.has_mcCollision()) + return; + auto mccol = collision.mcCollision(); + if (verboseInfo) + LOGF(info, "GeneratorId %d (%d)", mccol.getGeneratorId(), generatorIds->size()); + + if (std::find(generatorIds->begin(), generatorIds->end(), mccol.getGeneratorId()) != generatorIds->end()) { + if (verboseInfo) + LOGF(info, "Event with good generatorId"); + processReco(std::string("MCreco"), collision, bcs, tracks, fwdtracks, fv0as, ft0s, fdds); + } + } + PROCESS_SWITCH(SGCandProducer, processMcData, "Produce UD tables with MC data", false); }; struct McSGCandProducer { - // MC tables - Configurable verboseInfoMC{"verboseInfoMC", false, "Print general info to terminal; default it false."}; - Produces outputMcCollisions; - Produces outputMcParticles; - Produces outputMcCollsLabels; - Produces outputMcTrackLabels; - - // save all McTruth, even if the collisions is not reconstructed - Configurable> generatorIds{"generatorIds", std::vector{-1}, "MC generatorIds to process"}; - Configurable saveAllMcCollisions{"saveAllMcCollisions", true, "save all McCollisions"}; - - using CCs = soa::Join; - using BCs = soa::Join; - using TCs = soa::Join; - using UDCCs = soa::Join; - using UDTCs = soa::Join; - - // prepare slices - SliceCache cache; - PresliceUnsorted mcPartsPerMcCollision = aod::mcparticle::mcCollisionId; - Preslice udtracksPerUDCollision = aod::udtrack::udCollisionId; - - // initialize histogram registry - HistogramRegistry registry{ - "registry", - {}}; - - template - void updateUDMcCollisions(TMcCollision const& mccol, uint64_t globBC) - { - // save mccol - outputMcCollisions(globBC, - mccol.generatorsID(), - mccol.posX(), - mccol.posY(), - mccol.posZ(), - mccol.t(), - mccol.weight(), - mccol.impactParameter()); - } - - template - void updateUDMcParticle(TMcParticle const& McPart, int64_t McCollisionId, std::map& mcPartIsSaved) - { - // save McPart - // mother and daughter indices are set to -1 - // ATTENTION: this can be improved to also include mother and daughter indices - std::vector newmids; - int32_t newdids[2] = {-1, -1}; - - // update UDMcParticles - if (mcPartIsSaved.find(McPart.globalIndex()) == mcPartIsSaved.end()) { - outputMcParticles(McCollisionId, - McPart.pdgCode(), - McPart.statusCode(), - McPart.flags(), - newmids, - newdids, - McPart.weight(), - McPart.px(), - McPart.py(), - McPart.pz(), - McPart.e()); - mcPartIsSaved[McPart.globalIndex()] = outputMcParticles.lastIndex(); - } - } - - template - void updateUDMcParticles(TMcParticles const& McParts, int64_t McCollisionId, std::map& mcPartIsSaved) - { - // save McParts - // new mother and daughter ids - std::vector newmids; - int32_t newdids[2] = {-1, -1}; - int64_t newval = -1; - - // Determine the particle indices within the UDMcParticles table - // before filling the table - // This is needed to be able to assign the new daughter indices - std::map oldnew; - auto lastId = outputMcParticles.lastIndex(); - for (const auto& mcpart : McParts) { - auto oldId = mcpart.globalIndex(); - if (mcPartIsSaved.find(oldId) != mcPartIsSaved.end()) { - oldnew[oldId] = mcPartIsSaved[oldId]; - } else { - lastId++; - oldnew[oldId] = lastId; - } - } - - // all particles of the McCollision are saved - for (const auto& mcpart : McParts) { - if (mcPartIsSaved.find(mcpart.globalIndex()) == mcPartIsSaved.end()) { - // mothers - newmids.clear(); - auto oldmids = mcpart.mothersIds(); - for (const auto& oldmid : oldmids) { - auto m = McParts.rawIteratorAt(oldmid); - if (verboseInfoMC) - LOGF(debug, " m %d", m.globalIndex()); - if (mcPartIsSaved.find(oldmid) != mcPartIsSaved.end()) { - newval = mcPartIsSaved[oldmid]; - } else { - newval = -1; - } - newmids.push_back(newval); - } - // daughters - auto olddids = mcpart.daughtersIds(); - for (uint ii = 0; ii < olddids.size(); ii++) { - if (oldnew.find(olddids[ii]) != oldnew.end()) { - newval = oldnew[olddids[ii]]; - } else { - newval = -1; - } - newdids[ii] = newval; - } - if (verboseInfoMC) - LOGF(debug, " ms %i ds %i", oldmids.size(), olddids.size()); - - // update UDMcParticles - outputMcParticles(McCollisionId, - mcpart.pdgCode(), - mcpart.statusCode(), - mcpart.flags(), - newmids, - newdids, - mcpart.weight(), - mcpart.px(), - mcpart.py(), - mcpart.pz(), - mcpart.e()); - mcPartIsSaved[mcpart.globalIndex()] = outputMcParticles.lastIndex(); - } - } - } - - template - void updateUDMcTrackLabel(TTrack const& udtrack, std::map& mcPartIsSaved) - { - // udtrack (UDTCs) -> track (TCs) -> mcTrack (McParticles) -> udMcTrack (UDMcParticles) - auto trackId = udtrack.trackId(); - if (trackId >= 0) { - auto track = udtrack.template track_as(); - auto mcTrackId = track.mcParticleId(); - if (mcTrackId >= 0) { - if (mcPartIsSaved.find(mcTrackId) != mcPartIsSaved.end()) { - outputMcTrackLabels(mcPartIsSaved[mcTrackId], track.mcMask()); - } else { - outputMcTrackLabels(-1, track.mcMask()); - } - } else { - outputMcTrackLabels(-1, track.mcMask()); - } - } else { - outputMcTrackLabels(-1, -1); - } - } - - template - void updateUDMcTrackLabels(TTrack const& udtracks, std::map& mcPartIsSaved) - { - // loop over all tracks - for (const auto& udtrack : udtracks) { - // udtrack (UDTCs) -> track (TCs) -> mcTrack (McParticles) -> udMcTrack (UDMcParticles) - auto trackId = udtrack.trackId(); - if (trackId >= 0) { - auto track = udtrack.template track_as(); - auto mcTrackId = track.mcParticleId(); - if (mcTrackId >= 0) { - if (mcPartIsSaved.find(mcTrackId) != mcPartIsSaved.end()) { - outputMcTrackLabels(mcPartIsSaved[mcTrackId], track.mcMask()); - } else { - outputMcTrackLabels(-1, track.mcMask()); - } - } else { - outputMcTrackLabels(-1, track.mcMask()); - } - } else { - outputMcTrackLabels(-1, -1); - } - } - } - - // updating McTruth data and links to reconstructed data - void procWithSgCand(aod::McCollisions const& mccols, aod::McParticles const& mcparts, - UDCCs const& sgcands, UDTCs const& udtracks) - { - // use a hash table to keep track of the McCollisions which have been added to the UDMcCollision table - // {McCollisionId : udMcCollisionId} - // similar for the McParticles which have been added to the UDMcParticle table - // {McParticleId : udMcParticleId} - std::map mcColIsSaved; - std::map mcPartIsSaved; - - // loop over McCollisions and UDCCs simultaneously - auto mccol = mccols.iteratorAt(0); - auto mcOfInterest = std::find(generatorIds->begin(), generatorIds->end(), mccol.getGeneratorId()) != generatorIds->end(); - auto lastmccol = mccols.iteratorAt(mccols.size() - 1); - auto mccolAtEnd = false; - - auto sgcand = sgcands.iteratorAt(0); - auto lastsgcand = sgcands.iteratorAt(sgcands.size() - 1); - auto sgcandAtEnd = false; - - // advance dgcand and mccol until both are AtEnd - int64_t mccolId = mccol.globalIndex(); - int64_t mcsgId = -1; - bool goon = true; - while (goon) { - auto globBC = mccol.bc_as().globalBC(); - // check if dgcand has an associated McCollision - if (sgcand.has_collision()) { - auto sgcandCol = sgcand.collision_as(); - // colId = sgcandCol.globalIndex(); - if (sgcandCol.has_mcCollision()) { - mcsgId = sgcandCol.mcCollision().globalIndex(); - } else { - mcsgId = -1; - } - } else { - mcsgId = -1; - } - if (verboseInfoMC) - LOGF(info, "\nStart of loop mcsgId %d mccolId %d", mcsgId, mccolId); - - // two cases to consider - // 1. mcdgId <= mccolId: the event to process is a dgcand. In this case the Mc tables as well as the McLabel tables are updated - // 2. mccolId < mcdgId: the event to process is an MC event of interest without reconstructed dgcand. In this case only the Mc tables are updated - if ((!sgcandAtEnd && !mccolAtEnd && (mcsgId <= mccolId)) || mccolAtEnd) { - // this is case 1. - if (verboseInfoMC) - LOGF(info, "Doing case 1 with mcsgId %d", mcsgId); - - // update UDMcCollisions and UDMcColsLabels (for each UDCollision -> UDMcCollisions) - // update UDMcParticles and UDMcTrackLabels (for each UDTrack -> UDMcParticles) - // get dgcand tracks - auto sgTracks = udtracks.sliceByCached(aod::udtrack::udCollisionId, sgcand.globalIndex(), cache); - - // If the sgcand has an associated McCollision then the McCollision and all associated - // McParticles are saved - // but only consider generated events of interest - if (mcsgId >= 0 && mcOfInterest) { - if (mcColIsSaved.find(mcsgId) == mcColIsSaved.end()) { - if (verboseInfoMC) - LOGF(info, " Saving McCollision %d", mcsgId); - // update UDMcCollisions - auto sgcandMcCol = sgcand.collision_as().mcCollision(); - updateUDMcCollisions(sgcandMcCol, globBC); - mcColIsSaved[mcsgId] = outputMcCollisions.lastIndex(); - } - - // update UDMcColsLabels (for each UDCollision -> UDMcCollisions) - outputMcCollsLabels(mcColIsSaved[mcsgId]); - - // update UDMcParticles - auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mcsgId); - updateUDMcParticles(mcPartsSlice, mcColIsSaved[mcsgId], mcPartIsSaved); - - // update UDMcTrackLabels (for each UDTrack -> UDMcParticles) - updateUDMcTrackLabels(sgTracks, mcPartIsSaved); - - } else { - // If the sgcand has no associated McCollision then only the McParticles which are associated - // with the tracks of the sgcand are saved - if (verboseInfoMC) - LOGF(info, " Saving McCollision %d", -1); - - // update UDMcColsLabels (for each UDCollision -> UDMcCollisions) - outputMcCollsLabels(-1); - - // update UDMcParticles and UDMcTrackLabels (for each UDTrack -> UDMcParticles) - // loop over tracks of dgcand - for (const auto& sgtrack : sgTracks) { - if (sgtrack.has_track()) { - auto track = sgtrack.track_as(); - if (track.has_mcParticle()) { - auto mcPart = track.mcParticle(); - auto mcCol = mcPart.mcCollision(); - if (mcColIsSaved.find(mcCol.globalIndex()) == mcColIsSaved.end()) { - updateUDMcCollisions(mcCol, globBC); - mcColIsSaved[mcCol.globalIndex()] = outputMcCollisions.lastIndex(); - } - updateUDMcParticle(mcPart, mcColIsSaved[mcCol.globalIndex()], mcPartIsSaved); - updateUDMcTrackLabel(sgtrack, mcPartIsSaved); - } else { - outputMcTrackLabels(-1, track.mcMask()); - } - } else { - outputMcTrackLabels(-1, -1); - } - } - } - // advance sgcand - if (sgcand != lastsgcand) { - sgcand++; - } else { - sgcandAtEnd = true; - } - } else { - // this is case 2. - if (verboseInfoMC) - LOGF(info, "Doing case 2"); - - // update UDMcCollisions and UDMcParticles - // but only consider generated events of interest - if (mcOfInterest && mcColIsSaved.find(mccolId) == mcColIsSaved.end()) { - if (verboseInfoMC) - LOGF(info, " Saving McCollision %d", mccolId); - // update UDMcCollisions - updateUDMcCollisions(mccol, globBC); - mcColIsSaved[mccolId] = outputMcCollisions.lastIndex(); - - // update UDMcParticles - auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mccolId); - updateUDMcParticles(mcPartsSlice, mcColIsSaved[mccolId], mcPartIsSaved); - } - - // advance mccol - if (mccol != lastmccol) { - mccol++; - mcOfInterest = std::find(generatorIds->begin(), generatorIds->end(), mccol.getGeneratorId()) != generatorIds->end(); - mccolId = mccol.globalIndex(); - } else { - mccolAtEnd = true; - } - } - - goon = !sgcandAtEnd || !mccolAtEnd; - if (verboseInfoMC) - LOGF(info, "End of loop mcsgId %d mccolId %d", mcsgId, mccolId); - } - } - - // updating McTruth data only - void procWithoutSgCand(aod::McCollisions const& mccols, aod::McParticles const& mcparts) - { - // use a hash table to keep track of the McCollisions which have been added to the UDMcCollision table - // {McCollisionId : udMcCollisionId} - // similar for the McParticles which have been added to the UDMcParticle table - // {McParticleId : udMcParticleId} - std::map mcColIsSaved; - std::map mcPartIsSaved; - - // loop over McCollisions - for (auto const& mccol : mccols) { - int64_t mccolId = mccol.globalIndex(); - uint64_t globBC = mccol.bc_as().globalBC(); - - // update UDMcCollisions and UDMcParticles - if (mcColIsSaved.find(mccolId) == mcColIsSaved.end()) { - if (verboseInfoMC) - LOGF(info, " Saving McCollision %d", mccolId); - - // update UDMcCollisions - updateUDMcCollisions(mccol, globBC); - mcColIsSaved[mccolId] = outputMcCollisions.lastIndex(); - - // update UDMcParticles - auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mccolId); - updateUDMcParticles(mcPartsSlice, mcColIsSaved[mccolId], mcPartIsSaved); - } - } - } - - void init(InitContext& context) - { - // add histograms for the different process functions - if (context.mOptions.get("processMC")) { - registry.add("mcTruth/collisions", "Number of associated collisions", {HistType::kTH1F, {{11, -0.5, 10.5}}}); - registry.add("mcTruth/collType", "Collision type", {HistType::kTH1F, {{5, -0.5, 4.5}}}); - registry.add("mcTruth/IVMpt", "Invariant mass versus p_{T}", {HistType::kTH2F, {{150, 0.0, 3.0}, {150, 0.0, 3.0}}}); - } - } - - // process function for MC data - // save the MC truth of all events of interest and of the DG events - void processMC(aod::McCollisions const& mccols, aod::McParticles const& mcparts, - UDCCs const& sgcands, UDTCs const& udtracks, - CCs const& /*collisions*/, BCs const& /*bcs*/, TCs const& /*tracks*/) - { - if (verboseInfoMC) { - LOGF(info, "Number of McCollisions %d", mccols.size()); - LOGF(info, "Number of SG candidates %d", sgcands.size()); - LOGF(info, "Number of UD tracks %d", udtracks.size()); - } - if (mccols.size() > 0) { - if (sgcands.size() > 0) { - procWithSgCand(mccols, mcparts, sgcands, udtracks); - } else { - if (saveAllMcCollisions) { - procWithoutSgCand(mccols, mcparts); - } - } - } - } - PROCESS_SWITCH(McSGCandProducer, processMC, "Produce MC tables", false); - - void processDummy(aod::Collisions const& /*collisions*/) - { - // do nothing - if (verboseInfoMC) - LOGF(info, "Running dummy process function!"); - } - PROCESS_SWITCH(McSGCandProducer, processDummy, "Dummy function", true); + // MC tables + Configurable verboseInfoMC{"verboseInfoMC", false, "Print general info to terminal; default it false."}; + Produces outputMcCollisions; + Produces outputMcParticles; + Produces outputMcCollsLabels; + Produces outputMcTrackLabels; + + // save all McTruth, even if the collisions is not reconstructed + Configurable> generatorIds{"generatorIds", std::vector{-1}, "MC generatorIds to process"}; + Configurable saveAllMcCollisions{"saveAllMcCollisions", true, "save all McCollisions"}; + + using CCs = soa::Join; + using BCs = soa::Join; + using TCs = soa::Join; + using UDCCs = soa::Join; + using UDTCs = soa::Join; + + // prepare slices + SliceCache cache; + PresliceUnsorted mcPartsPerMcCollision = aod::mcparticle::mcCollisionId; + Preslice udtracksPerUDCollision = aod::udtrack::udCollisionId; + + // initialize histogram registry + HistogramRegistry registry{ + "registry", + {}}; + + template + void updateUDMcCollisions(TMcCollision const& mccol, uint64_t globBC) + { + // save mccol + outputMcCollisions(globBC, + mccol.generatorsID(), + mccol.posX(), + mccol.posY(), + mccol.posZ(), + mccol.t(), + mccol.weight(), + mccol.impactParameter()); + } + + template + void updateUDMcParticle(TMcParticle const& McPart, int64_t McCollisionId, std::map& mcPartIsSaved) + { + // save McPart + // mother and daughter indices are set to -1 + // ATTENTION: this can be improved to also include mother and daughter indices + std::vector newmids; + int32_t newdids[2] = {-1, -1}; + + // update UDMcParticles + if (mcPartIsSaved.find(McPart.globalIndex()) == mcPartIsSaved.end()) { + outputMcParticles(McCollisionId, + McPart.pdgCode(), + McPart.statusCode(), + McPart.flags(), + newmids, + newdids, + McPart.weight(), + McPart.px(), + McPart.py(), + McPart.pz(), + McPart.e()); + mcPartIsSaved[McPart.globalIndex()] = outputMcParticles.lastIndex(); + } + } + + template + void updateUDMcParticles(TMcParticles const& McParts, int64_t McCollisionId, std::map& mcPartIsSaved) + { + // save McParts + // new mother and daughter ids + std::vector newmids; + int32_t newdids[2] = {-1, -1}; + int64_t newval = -1; + + // Determine the particle indices within the UDMcParticles table + // before filling the table + // This is needed to be able to assign the new daughter indices + std::map oldnew; + auto lastId = outputMcParticles.lastIndex(); + for (const auto& mcpart : McParts) { + auto oldId = mcpart.globalIndex(); + if (mcPartIsSaved.find(oldId) != mcPartIsSaved.end()) { + oldnew[oldId] = mcPartIsSaved[oldId]; + } else { + lastId++; + oldnew[oldId] = lastId; + } + } + + // all particles of the McCollision are saved + for (const auto& mcpart : McParts) { + if (mcPartIsSaved.find(mcpart.globalIndex()) == mcPartIsSaved.end()) { + // mothers + newmids.clear(); + auto oldmids = mcpart.mothersIds(); + for (const auto& oldmid : oldmids) { + auto m = McParts.rawIteratorAt(oldmid); + if (verboseInfoMC) + LOGF(debug, " m %d", m.globalIndex()); + if (mcPartIsSaved.find(oldmid) != mcPartIsSaved.end()) { + newval = mcPartIsSaved[oldmid]; + } else { + newval = -1; + } + newmids.push_back(newval); + } + // daughters + auto olddids = mcpart.daughtersIds(); + for (uint ii = 0; ii < olddids.size(); ii++) { + if (oldnew.find(olddids[ii]) != oldnew.end()) { + newval = oldnew[olddids[ii]]; + } else { + newval = -1; + } + newdids[ii] = newval; + } + if (verboseInfoMC) + LOGF(debug, " ms %i ds %i", oldmids.size(), olddids.size()); + + // update UDMcParticles + outputMcParticles(McCollisionId, + mcpart.pdgCode(), + mcpart.statusCode(), + mcpart.flags(), + newmids, + newdids, + mcpart.weight(), + mcpart.px(), + mcpart.py(), + mcpart.pz(), + mcpart.e()); + mcPartIsSaved[mcpart.globalIndex()] = outputMcParticles.lastIndex(); + } + } + } + + template + void updateUDMcTrackLabel(TTrack const& udtrack, std::map& mcPartIsSaved) + { + // udtrack (UDTCs) -> track (TCs) -> mcTrack (McParticles) -> udMcTrack (UDMcParticles) + auto trackId = udtrack.trackId(); + if (trackId >= 0) { + auto track = udtrack.template track_as(); + auto mcTrackId = track.mcParticleId(); + if (mcTrackId >= 0) { + if (mcPartIsSaved.find(mcTrackId) != mcPartIsSaved.end()) { + outputMcTrackLabels(mcPartIsSaved[mcTrackId], track.mcMask()); + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, -1); + } + } + + template + void updateUDMcTrackLabels(TTrack const& udtracks, std::map& mcPartIsSaved) + { + // loop over all tracks + for (const auto& udtrack : udtracks) { + // udtrack (UDTCs) -> track (TCs) -> mcTrack (McParticles) -> udMcTrack (UDMcParticles) + auto trackId = udtrack.trackId(); + if (trackId >= 0) { + auto track = udtrack.template track_as(); + auto mcTrackId = track.mcParticleId(); + if (mcTrackId >= 0) { + if (mcPartIsSaved.find(mcTrackId) != mcPartIsSaved.end()) { + outputMcTrackLabels(mcPartIsSaved[mcTrackId], track.mcMask()); + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, -1); + } + } + } + + // updating McTruth data and links to reconstructed data + void procWithSgCand(aod::McCollisions const& mccols, aod::McParticles const& mcparts, + UDCCs const& sgcands, UDTCs const& udtracks) + { + // use a hash table to keep track of the McCollisions which have been added to the UDMcCollision table + // {McCollisionId : udMcCollisionId} + // similar for the McParticles which have been added to the UDMcParticle table + // {McParticleId : udMcParticleId} + std::map mcColIsSaved; + std::map mcPartIsSaved; + + // loop over McCollisions and UDCCs simultaneously + auto mccol = mccols.iteratorAt(0); + auto mcOfInterest = std::find(generatorIds->begin(), generatorIds->end(), mccol.getGeneratorId()) != generatorIds->end(); + auto lastmccol = mccols.iteratorAt(mccols.size() - 1); + auto mccolAtEnd = false; + + auto sgcand = sgcands.iteratorAt(0); + auto lastsgcand = sgcands.iteratorAt(sgcands.size() - 1); + auto sgcandAtEnd = false; + + // advance dgcand and mccol until both are AtEnd + int64_t mccolId = mccol.globalIndex(); + int64_t mcsgId = -1; + bool goon = true; + while (goon) { + auto globBC = mccol.bc_as().globalBC(); + // check if dgcand has an associated McCollision + if (sgcand.has_collision()) { + auto sgcandCol = sgcand.collision_as(); + // colId = sgcandCol.globalIndex(); + if (sgcandCol.has_mcCollision()) { + mcsgId = sgcandCol.mcCollision().globalIndex(); + } else { + mcsgId = -1; + } + } else { + mcsgId = -1; + } + if (verboseInfoMC) + LOGF(info, "\nStart of loop mcsgId %d mccolId %d", mcsgId, mccolId); + + // two cases to consider + // 1. mcdgId <= mccolId: the event to process is a dgcand. In this case the Mc tables as well as the McLabel tables are updated + // 2. mccolId < mcdgId: the event to process is an MC event of interest without reconstructed dgcand. In this case only the Mc tables are updated + if ((!sgcandAtEnd && !mccolAtEnd && (mcsgId <= mccolId)) || mccolAtEnd) { + // this is case 1. + if (verboseInfoMC) + LOGF(info, "Doing case 1 with mcsgId %d", mcsgId); + + // update UDMcCollisions and UDMcColsLabels (for each UDCollision -> UDMcCollisions) + // update UDMcParticles and UDMcTrackLabels (for each UDTrack -> UDMcParticles) + // get dgcand tracks + auto sgTracks = udtracks.sliceByCached(aod::udtrack::udCollisionId, sgcand.globalIndex(), cache); + + // If the sgcand has an associated McCollision then the McCollision and all associated + // McParticles are saved + // but only consider generated events of interest + if (mcsgId >= 0 && mcOfInterest) { + if (mcColIsSaved.find(mcsgId) == mcColIsSaved.end()) { + if (verboseInfoMC) + LOGF(info, " Saving McCollision %d", mcsgId); + // update UDMcCollisions + auto sgcandMcCol = sgcand.collision_as().mcCollision(); + updateUDMcCollisions(sgcandMcCol, globBC); + mcColIsSaved[mcsgId] = outputMcCollisions.lastIndex(); + } + + // update UDMcColsLabels (for each UDCollision -> UDMcCollisions) + outputMcCollsLabels(mcColIsSaved[mcsgId]); + + // update UDMcParticles + auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mcsgId); + updateUDMcParticles(mcPartsSlice, mcColIsSaved[mcsgId], mcPartIsSaved); + + // update UDMcTrackLabels (for each UDTrack -> UDMcParticles) + updateUDMcTrackLabels(sgTracks, mcPartIsSaved); + + } else { + // If the sgcand has no associated McCollision then only the McParticles which are associated + // with the tracks of the sgcand are saved + if (verboseInfoMC) + LOGF(info, " Saving McCollision %d", -1); + + // update UDMcColsLabels (for each UDCollision -> UDMcCollisions) + outputMcCollsLabels(-1); + + // update UDMcParticles and UDMcTrackLabels (for each UDTrack -> UDMcParticles) + // loop over tracks of dgcand + for (const auto& sgtrack : sgTracks) { + if (sgtrack.has_track()) { + auto track = sgtrack.track_as(); + if (track.has_mcParticle()) { + auto mcPart = track.mcParticle(); + auto mcCol = mcPart.mcCollision(); + if (mcColIsSaved.find(mcCol.globalIndex()) == mcColIsSaved.end()) { + updateUDMcCollisions(mcCol, globBC); + mcColIsSaved[mcCol.globalIndex()] = outputMcCollisions.lastIndex(); + } + updateUDMcParticle(mcPart, mcColIsSaved[mcCol.globalIndex()], mcPartIsSaved); + updateUDMcTrackLabel(sgtrack, mcPartIsSaved); + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, -1); + } + } + } + // advance sgcand + if (sgcand != lastsgcand) { + sgcand++; + } else { + sgcandAtEnd = true; + } + } else { + // this is case 2. + if (verboseInfoMC) + LOGF(info, "Doing case 2"); + + // update UDMcCollisions and UDMcParticles + // but only consider generated events of interest + if (mcOfInterest && mcColIsSaved.find(mccolId) == mcColIsSaved.end()) { + if (verboseInfoMC) + LOGF(info, " Saving McCollision %d", mccolId); + // update UDMcCollisions + updateUDMcCollisions(mccol, globBC); + mcColIsSaved[mccolId] = outputMcCollisions.lastIndex(); + + // update UDMcParticles + auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mccolId); + updateUDMcParticles(mcPartsSlice, mcColIsSaved[mccolId], mcPartIsSaved); + } + + // advance mccol + if (mccol != lastmccol) { + mccol++; + mcOfInterest = std::find(generatorIds->begin(), generatorIds->end(), mccol.getGeneratorId()) != generatorIds->end(); + mccolId = mccol.globalIndex(); + } else { + mccolAtEnd = true; + } + } + + goon = !sgcandAtEnd || !mccolAtEnd; + if (verboseInfoMC) + LOGF(info, "End of loop mcsgId %d mccolId %d", mcsgId, mccolId); + } + } + + // updating McTruth data only + void procWithoutSgCand(aod::McCollisions const& mccols, aod::McParticles const& mcparts) + { + // use a hash table to keep track of the McCollisions which have been added to the UDMcCollision table + // {McCollisionId : udMcCollisionId} + // similar for the McParticles which have been added to the UDMcParticle table + // {McParticleId : udMcParticleId} + std::map mcColIsSaved; + std::map mcPartIsSaved; + + // loop over McCollisions + for (auto const& mccol : mccols) { + int64_t mccolId = mccol.globalIndex(); + uint64_t globBC = mccol.bc_as().globalBC(); + + // update UDMcCollisions and UDMcParticles + if (mcColIsSaved.find(mccolId) == mcColIsSaved.end()) { + if (verboseInfoMC) + LOGF(info, " Saving McCollision %d", mccolId); + + // update UDMcCollisions + updateUDMcCollisions(mccol, globBC); + mcColIsSaved[mccolId] = outputMcCollisions.lastIndex(); + + // update UDMcParticles + auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mccolId); + updateUDMcParticles(mcPartsSlice, mcColIsSaved[mccolId], mcPartIsSaved); + } + } + } + + void init(InitContext& context) + { + // add histograms for the different process functions + if (context.mOptions.get("processMC")) { + registry.add("mcTruth/collisions", "Number of associated collisions", {HistType::kTH1F, {{11, -0.5, 10.5}}}); + registry.add("mcTruth/collType", "Collision type", {HistType::kTH1F, {{5, -0.5, 4.5}}}); + registry.add("mcTruth/IVMpt", "Invariant mass versus p_{T}", {HistType::kTH2F, {{150, 0.0, 3.0}, {150, 0.0, 3.0}}}); + } + } + + // process function for MC data + // save the MC truth of all events of interest and of the DG events + void processMC(aod::McCollisions const& mccols, aod::McParticles const& mcparts, + UDCCs const& sgcands, UDTCs const& udtracks, + CCs const& /*collisions*/, BCs const& /*bcs*/, TCs const& /*tracks*/) + { + if (verboseInfoMC) { + LOGF(info, "Number of McCollisions %d", mccols.size()); + LOGF(info, "Number of SG candidates %d", sgcands.size()); + LOGF(info, "Number of UD tracks %d", udtracks.size()); + } + if (mccols.size() > 0) { + if (sgcands.size() > 0) { + procWithSgCand(mccols, mcparts, sgcands, udtracks); + } else { + if (saveAllMcCollisions) { + procWithoutSgCand(mccols, mcparts); + } + } + } + } + PROCESS_SWITCH(McSGCandProducer, processMC, "Produce MC tables", false); + + void processDummy(aod::Collisions const& /*collisions*/) + { + // do nothing + if (verboseInfoMC) + LOGF(info, "Running dummy process function!"); + } + PROCESS_SWITCH(McSGCandProducer, processDummy, "Dummy function", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - WorkflowSpec workflow{ - adaptAnalysisTask(cfgc, TaskName{"sgcandproducer"}), - adaptAnalysisTask(cfgc, TaskName{"mcsgcandproducer"})}; - return workflow; + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, TaskName{"sgcandproducer"}), + adaptAnalysisTask(cfgc, TaskName{"mcsgcandproducer"})}; + return workflow; }