From e58155cda739071a4812a46f8fab8806ecc51fde Mon Sep 17 00:00:00 2001 From: Paul Date: Fri, 7 Apr 2023 16:45:39 +0200 Subject: [PATCH 1/4] . added tutorial examples . some updates DGPIDSelector . adaptions of DGCandAnalyzer, diffQA, diffMCQA --- PWGUD/Core/DGPIDSelector.cxx | 71 +-- PWGUD/Core/DGPIDSelector.h | 105 ++-- PWGUD/TableProducer/DGCandProducer.cxx | 16 +- PWGUD/Tasks/CMakeLists.txt | 17 +- PWGUD/Tasks/DGCandAnalyzer.cxx | 1 + PWGUD/Tasks/UDTutorial_01.cxx | 147 ++++++ PWGUD/Tasks/UDTutorial_02a.cxx | 183 +++++++ PWGUD/Tasks/UDTutorial_02b.cxx | 173 +++++++ PWGUD/Tasks/diffMCQA.cxx | 691 ++++++++++++------------- PWGUD/Tasks/diffQA.cxx | 454 +++++++++------- 10 files changed, 1193 insertions(+), 665 deletions(-) create mode 100644 PWGUD/Tasks/UDTutorial_01.cxx create mode 100644 PWGUD/Tasks/UDTutorial_02a.cxx create mode 100644 PWGUD/Tasks/UDTutorial_02b.cxx diff --git a/PWGUD/Core/DGPIDSelector.cxx b/PWGUD/Core/DGPIDSelector.cxx index a011d30461b..a992d984abb 100644 --- a/PWGUD/Core/DGPIDSelector.cxx +++ b/PWGUD/Core/DGPIDSelector.cxx @@ -122,16 +122,29 @@ DGAnaparHolder::~DGAnaparHolder() void DGAnaparHolder::Print() { LOGF(info, " DGAnaparHolder"); - LOGF(info, " nCombine: %i", mNCombine); - LOGF(info, " max dcaxy: %f", mMaxDCAxy); - LOGF(info, " max dcaz: %f", mMaxDCAz); + LOGF(info, " min number tracks: %d", mMinNTracks); + LOGF(info, " max number tracks: %d", mMaxNTracks); + LOGF(info, " min fraction of PV contr. with TOF: %f", mMinRgtrwTOF); + LOGF(info, " max dcaxy: %f", mMaxDCAxy); + LOGF(info, " max dcaz: %f", mMaxDCAz); + LOGF(info, " min dBC: %d", mdBCMin); + LOGF(info, " max dBC: %d", mdBCMax); + LOGF(info, " min track pT: %f", mMinpt); + LOGF(info, " max track pT: %f", mMaxpt); + LOGF(info, " min eta: %f", mMineta); + LOGF(info, " max eta: %f", mMaxeta); + LOGF(info, " min alpha: %f", mMinAlpha); + LOGF(info, " max alpha: %f", mMaxAlpha); + LOGF(info, " min system pT: %f", mMinptsys); + LOGF(info, " max system pT: %f", mMaxptsys); + LOGF(info, " nCombine: %d", mNCombine); LOGF(info, " net charges"); for (auto ch : mNetCharges) { LOGF(info, " %i", ch); } LOGF(info, " PIDs"); for (auto pid : mDGPIDs) { - LOGF(info, " %f", pid); + LOGF(info, " %d", pid); } PIDCuts().Print(); } @@ -145,8 +158,11 @@ DGPIDCuts DGAnaparHolder::PIDCuts() // ----------------------------------------------------------------------------- void DGAnaparHolder::makeUniquePermutations() { + // reset + muniquePerms.clear(); + // all permutations of mNCombine elements - std::vector> perms; + std::vector> perms; permutations(mNCombine, perms); // compute unique permutations @@ -160,7 +176,7 @@ void DGAnaparHolder::makeUniquePermutations() cnt = -1; for (auto ind : perm) { cnt++; - perminfo[cnt] = static_cast(mDGPIDs[ind]); + perminfo[cnt] = mDGPIDs[ind]; } hashstr = ""; for (auto tok : perminfo) { @@ -183,9 +199,8 @@ void DGAnaparHolder::makeUniquePermutations() std::vector DGAnaparHolder::uniquePermutations() { // create unique permutations if not done already - if (muniquePerms.size() == 0) { + if (muniquePerms.size() < mNCombine) { makeUniquePermutations(); - LOGF(info, "Number of unique permutations: %i", muniquePerms.size() / mNCombine); } return muniquePerms; @@ -193,9 +208,8 @@ std::vector DGAnaparHolder::uniquePermutations() // ----------------------------------------------------------------------------- // find all permutations of n0 elements -void DGAnaparHolder::permutations(std::vector& ref, int n0, int np, std::vector>& perms) +void DGAnaparHolder::permutations(std::vector& ref, int n0, int np, std::vector>& perms) { - // create local reference auto ref2u = ref; @@ -205,7 +219,7 @@ void DGAnaparHolder::permutations(std::vector& ref, int n0, int np, std::v // create a new permutation // copy first n0-np elements from ref // then rotate last np elements of ref - std::vector perm(n0, 0); + std::vector perm(n0, 0); for (auto ii = 0; ii < n0 - np; ii++) { perm[ii] = ref2u[ii]; } @@ -233,7 +247,7 @@ void DGAnaparHolder::permutations(std::vector& ref, int n0, int np, std::v //----------------------------------------------------------------------------- // find all permutations of n0 elements -int DGAnaparHolder::permutations(int n0, std::vector>& perms) +int DGAnaparHolder::permutations(int n0, std::vector>& perms) { // initialize with first trivial combination perms.clear(); @@ -241,7 +255,7 @@ int DGAnaparHolder::permutations(int n0, std::vector>& perms) return 0; } - std::vector ref(n0, 0); + std::vector ref(n0, 0); for (auto ii = 0; ii < n0; ii++) { ref[ii] = ii; } @@ -323,8 +337,8 @@ int DGPIDSelector::pid2ind(int pid) // ----------------------------------------------------------------------------- // find selections of np out of n0 -void DGPIDSelector::combinations(int n0, std::vector& pool, int np, std::vector& inds, int n, - std::vector>& combs) +void DGPIDSelector::combinations(int n0, std::vector& pool, int np, std::vector& inds, int n, + std::vector>& combs) { // loop over pool for (auto ii = 0; ii < n0 - n; ii++) { @@ -335,7 +349,7 @@ void DGPIDSelector::combinations(int n0, std::vector& pool, int np, std::v // else get next inds if (np == 1) { - std::vector comb(n + 1, 0); + std::vector comb(n + 1, 0); for (uint ii = 0; ii < inds.size(); ii++) { comb[ii] = inds[ii]; } @@ -344,7 +358,7 @@ void DGPIDSelector::combinations(int n0, std::vector& pool, int np, std::v } else { auto n0new = n0 - ii; - std::vector newpool(n0new, 0); + std::vector newpool(n0new, 0); for (auto kk = 0; kk < n0new; kk++) { newpool[kk] = pool[kk + ii + 1]; } @@ -358,7 +372,7 @@ void DGPIDSelector::combinations(int n0, std::vector& pool, int np, std::v // ----------------------------------------------------------------------------- // find all possible selections of np out of n0 -int DGPIDSelector::combinations(int n0, int np, std::vector>& combs) +int DGPIDSelector::combinations(int n0, int np, std::vector>& combs) { // initialisations combs.clear(); @@ -366,11 +380,11 @@ int DGPIDSelector::combinations(int n0, int np, std::vector>& return 0; } - std::vector pool(n0, 0); + std::vector pool(n0, 0); for (auto ii = 0; ii < n0; ii++) { pool[ii] = ii; } - std::vector inds(np, 0); + std::vector inds(np, 0); // iterate recursively combinations(n0, pool, np, inds, 0, combs); @@ -379,21 +393,21 @@ int DGPIDSelector::combinations(int n0, int np, std::vector>& } // ----------------------------------------------------------------------------- -std::vector> DGPIDSelector::combinations(int nPool) +std::vector> DGPIDSelector::combinations(int nPool) { // get the unique permutations auto uniquePerms = mAnaPars.uniquePermutations(); auto numUniquePerms = uniquePerms.size() / mAnaPars.nCombine(); // all selections of nCombine elements from nPool elements - std::vector> combs; + std::vector> combs; combinations(nPool, mAnaPars.nCombine(), combs); // permute the combinations - std::vector> copes; + std::vector> copes; for (auto comb : combs) { for (auto ii = 0u; ii < numUniquePerms; ii++) { - std::vector cope(mAnaPars.nCombine(), 0); + std::vector cope(mAnaPars.nCombine(), 0); for (auto jj = 0; jj < mAnaPars.nCombine(); jj++) { auto ind = ii * mAnaPars.nCombine() + jj; cope[uniquePerms[ind]] = comb[jj]; @@ -402,15 +416,6 @@ std::vector> DGPIDSelector::combinations(int nPool) } } - // print copes - LOGF(debug, "copes"); - for (auto cope : copes) { - LOGF(debug, " cope"); - for (auto ind : cope) { - LOGF(debug, " %i", ind); - } - } - return copes; } diff --git a/PWGUD/Core/DGPIDSelector.h b/PWGUD/Core/DGPIDSelector.h index 8d66ea9f8d2..12a1b627fa3 100644 --- a/PWGUD/Core/DGPIDSelector.h +++ b/PWGUD/Core/DGPIDSelector.h @@ -107,15 +107,17 @@ struct DGPIDCuts { struct DGAnaparHolder { public: // constructor - DGAnaparHolder(int nCombine = 2, float maxDCAxy = 100., float maxDCAz = 100, + DGAnaparHolder(int MinNTracks = 0, int MaxNTracks = 10000, float minrgtrwTOF = 0., + float maxDCAxy = 100., float maxDCAz = 100, int dBCMin = 0, int dBCMax = 0, - float minptsys = 0.0, float maxptsys = 100.0, float minpt = 0.0, float maxpt = 100.0, float mineta = -2.0, float maxeta = 2.0, float minalpha = 0.0, float maxalpha = 3.2, - std::vector netCharges = {-2, -1, 0, 1, 2}, - std::vector DGPIDs = {211, 211}, - std::vector DGPIDCutValues = {}) : mNCombine{nCombine}, mdBCMin{dBCMin}, mdBCMax{dBCMax}, mMaxDCAxy{maxDCAxy}, mMaxDCAz{maxDCAz}, mMinpt{minpt}, mMaxpt{maxpt}, mMineta{mineta}, mMaxeta{maxeta}, mMinptsys{minptsys}, mMaxptsys{maxptsys}, mMinAlpha{minalpha}, mMaxAlpha{maxalpha}, mNetCharges{netCharges}, mDGPIDs{DGPIDs}, mDGPIDCutValues{DGPIDCutValues} + float minptsys = 0.0, float maxptsys = 100.0, + int nCombine = 2, + std::vector netCharges = {0}, + std::vector DGPIDs = {211, 211}, + std::vector DGPIDCutValues = {}) : mMinNTracks{MinNTracks}, mMaxNTracks{MaxNTracks}, mMinRgtrwTOF{minrgtrwTOF}, mMaxDCAxy{maxDCAxy}, mMaxDCAz{maxDCAz}, mdBCMin{dBCMin}, mdBCMax{dBCMax}, mMinpt{minpt}, mMaxpt{maxpt}, mMineta{mineta}, mMaxeta{maxeta}, mMinAlpha{minalpha}, mMaxAlpha{maxalpha}, mMinptsys{minptsys}, mMaxptsys{maxptsys}, mNCombine{nCombine}, mNetCharges{netCharges}, mDGPIDs{DGPIDs}, mDGPIDCutValues{DGPIDCutValues} { if (mdBCMin < -16) { mdBCMin = -16; @@ -127,69 +129,60 @@ struct DGAnaparHolder { } else if (mdBCMax > 15) { mdBCMax = 15; } - makeUniquePermutations(); } ~DGAnaparHolder(); - // getter + // helper void Print(); - int nCombine() const { return mNCombine; } + + // getter + int minNTracks() const { return mMinNTracks; } + int maxNTracks() const { return mMaxNTracks; } + float minRgtrwTOF() const { return mMinRgtrwTOF; } + float maxDCAxy() const { return mMaxDCAxy; } + float maxDCAz() const { return mMaxDCAz; } int dBCMin() const { return mdBCMin; } int dBCMax() const { return mdBCMax; } - float maxDCAxy() { return mMaxDCAxy; } - float maxDCAz() { return mMaxDCAz; } - float minptsys() { return mMinptsys; } - float maxptsys() { return mMaxptsys; } - float minpt() { return mMinpt; } - float maxpt() { return mMaxpt; } - float mineta() { return mMineta; } - float maxeta() { return mMaxeta; } - float minAlpha() { return mMinAlpha; } - float maxAlpha() { return mMaxAlpha; } - std::vector netCharges() { return mNetCharges; } - std::vector PIDs() { return mDGPIDs; } + float minpt() const { return mMinpt; } + float maxpt() const { return mMaxpt; } + float mineta() const { return mMineta; } + float maxeta() const { return mMaxeta; } + float minAlpha() const { return mMinAlpha; } + float maxAlpha() const { return mMaxAlpha; } + float minptsys() const { return mMinptsys; } + float maxptsys() const { return mMaxptsys; } + int nCombine() const { return mNCombine; } + std::vector netCharges() const { return mNetCharges; } + std::vector PIDs() const { return mDGPIDs; } DGPIDCuts PIDCuts(); std::vector uniquePermutations(); private: // helper functions - void permutations(std::vector& ref, int n0, int np, std::vector>& perms); - int permutations(int n0, std::vector>& perms); + void permutations(std::vector& ref, int n0, int np, std::vector>& perms); + int permutations(int n0, std::vector>& perms); void makeUniquePermutations(); - // number of tracks to combine - int mNCombine; - - // BBFlags - int mdBCMin; - int mdBCMax; - - // dca of tracks + // arguments + int mMinNTracks; + int mMaxNTracks; + float mMinRgtrwTOF; float mMaxDCAxy; float mMaxDCAz; - - // pt-range of tracks + int mdBCMin; + int mdBCMax; float mMinpt; float mMaxpt; - - // eta range of tracks float mMineta; float mMaxeta; - - // pt-range of system - float mMinptsys; - float mMaxptsys; - - // alpha-range when 2 tracks float mMinAlpha; float mMaxAlpha; - - // net charge of all tracks + float mMinptsys; + float mMaxptsys; + int mNCombine; std::vector mNetCharges; - - // PID information - std::vector mDGPIDs; + std::vector mDGPIDs; std::vector mDGPIDCutValues; // unique permutations @@ -204,7 +197,7 @@ struct DGParticle { public: DGParticle(); template - DGParticle(TDatabasePDG* pdg, DGAnaparHolder anaPars, TTrack const& tracks, std::vector comb) + DGParticle(TDatabasePDG* pdg, DGAnaparHolder anaPars, TTrack const& tracks, std::vector comb) { // compute invariant mass TLorentzVector lvtmp; @@ -227,7 +220,7 @@ struct DGParticle { // getter void Print(); - std::vector trkinds() { return mtrkinds; } + std::vector trkinds() { return mtrkinds; } float M() { return mIVM.M(); } float Perp() { return mIVM.Perp(); } @@ -236,7 +229,7 @@ struct DGParticle { TLorentzVector mIVM; // indices of tracks included - std::vector mtrkinds; + std::vector mtrkinds; // ClassDefNV(DGParticle, 1); }; @@ -254,14 +247,13 @@ struct DGPIDSelector { // getters void Print(); template - bool isGoodCombination(std::vector comb, TTrack const& tracks) + bool isGoodCombination(std::vector comb, TTrack const& tracks) { // compute net charge of track combination int netCharge = 0.; for (auto const& ind : comb) { netCharge += (tracks.begin() + ind).sign(); } - LOGF(debug, "Net charge %i", netCharge); // is this in the list of accepted net charges? auto netCharges = mAnaPars.netCharges(); @@ -309,19 +301,16 @@ struct DGPIDSelector { for (auto pidcut : pidcuts) { // skip cut if it does not apply to this track - LOGF(debug, "nPart %i %i, Type %i Apply %i", pidcut.nPart(), cnt, pidcut.cutType(), pidcut.cutApply()); if (pidcut.nPart() != cnt || pidcut.cutApply() <= 0) { continue; } // check pt for pid cut - LOGF(debug, "pT %f %f %f", track.pt(), pidcut.cutPtMin(), pidcut.cutPtMax()); if (track.pt() < pidcut.cutPtMin() || track.pt() > pidcut.cutPtMax()) { continue; } // is detector information required - LOGF(debug, "TPC %i TOF %i", track.hasTPC(), track.hasTOF()); if (pidcut.cutApply() == 2) { if (pidcut.cutDetector() == 1 && !track.hasTPC()) { return false; @@ -332,7 +321,6 @@ struct DGPIDSelector { } // get detector value - LOGF(debug, "cutPID %i", pidcut.cutPID()); float detValue = 0.; if (pidcut.cutDetector() == 1) { if (!track.hasTPC()) { @@ -386,6 +374,7 @@ struct DGPIDSelector { for (auto comb : combs) { // is combination compatible with netCharge requirements? if (!isGoodCombination(comb, tracks)) { + LOGF(info, "Not a good combination"); continue; } // is tracks compatible with PID requirements? @@ -463,10 +452,10 @@ struct DGPIDSelector { TDatabasePDG* fPDG; // helper functions for computeIVMs - void combinations(int n0, std::vector& pool, int np, std::vector& inds, int n, - std::vector>& combs); - int combinations(int n0, int np, std::vector>& combs); - std::vector> combinations(int nPool); + void combinations(int n0, std::vector& pool, int np, std::vector& inds, int n, + std::vector>& combs); + int combinations(int n0, int np, std::vector>& combs); + std::vector> combinations(int nPool); // ClassDefNV(DGPIDSelector, 1); }; diff --git a/PWGUD/TableProducer/DGCandProducer.cxx b/PWGUD/TableProducer/DGCandProducer.cxx index 3213ceb69fc..fd7570880f2 100644 --- a/PWGUD/TableProducer/DGCandProducer.cxx +++ b/PWGUD/TableProducer/DGCandProducer.cxx @@ -82,11 +82,8 @@ struct DGCandProducer { using MCTC = MCTCs::iterator; // extract FIT information - upchelpers::FITInfo getFITinfo(uint64_t const& bcnum, BCs const& bcs, aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds) + void getFITinfo(upchelpers::FITInfo& info, uint64_t const& bcnum, BCs const& bcs, aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds) { - // FITinfo - upchelpers::FITInfo info{}; - // find bc with globalBC = bcnum Partition selbc = aod::bc::globalBC == bcnum; selbc.bindTable(bcs); @@ -177,8 +174,6 @@ struct DGCandProducer { if (bc2u.selection()[evsel::kIsBBFDC]) SETBIT(info.BBFDDCpf, bit); } - - return info; } // function to update UDTracks, UDTracksCov, UDTracksDCA, UDTracksPID, UDTracksExtra, UDTracksFlag, @@ -296,19 +291,22 @@ struct DGCandProducer { return; } auto bc = collision.foundBC_as(); + LOGF(debug, " BC id %d", bc.globalBC()); // obtain slice of compatible BCs auto bcRange = udhelpers::compatibleBCs(collision, diffCuts.NDtcoll(), bcs, diffCuts.minNBCs()); + LOGF(debug, " Size of bcRange %d", bcRange.size()); // apply DG selection auto isDGEvent = dgSelector.IsSelected(diffCuts, collision, bcRange, tracks, fwdtracks); // save DG candidates if (isDGEvent == 0) { - LOGF(debug, " Data: good collision!"); + LOGF(debug, " Data: good collision!"); // fill FITInfo - upchelpers::FITInfo fitInfo = getFITinfo(bc.globalBC(), bcs, ft0s, fv0as, fdds); + upchelpers::FITInfo fitInfo{}; + getFITinfo(fitInfo, bc.globalBC(), bcs, ft0s, fv0as, fdds); // update DG candidates tables auto rtrwTOF = udhelpers::rPVtrwTOF(tracks, collision.numContrib()); @@ -347,7 +345,7 @@ struct DGCandProducer { pt2 = tr.pt() * tr.sign(); signalTPC2 = tr.tpcSignal(); } - LOGF(debug, "track[%d] %d pT %f ITS %d TPC %d TRD %d TOF %d", + LOGF(debug, " track[%d] %d pT %f ITS %d TPC %d TRD %d TOF %d", cnt, tr.isGlobalTrack(), tr.pt(), tr.itsNCls(), tr.tpcNClsCrossedRows(), tr.hasTRD(), tr.hasTOF()); } } diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index b180cf40554..5c242b7db90 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -39,6 +39,21 @@ o2physics_add_dpl_workflow(dgcand-analyzer PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGCutparHolder O2Physics::UDGoodRunSelector O2Physics::DGPIDSelector O2Physics::UDFSParser COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(udtutorial-01 + SOURCES UDTutorial_01.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(udtutorial-02a + SOURCES UDTutorial_02a.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(udtutorial-02b + SOURCES UDTutorial_02b.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(upccand-analyzer SOURCES UPCCandidateAnalyzer.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase @@ -52,4 +67,4 @@ o2physics_add_dpl_workflow(upccand-producer-qa o2physics_add_dpl_workflow(upc-mft SOURCES upcMft.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsBase O2::DetectorsCommonDataFormats - COMPONENT_NAME Analysis) \ No newline at end of file + COMPONENT_NAME Analysis) diff --git a/PWGUD/Tasks/DGCandAnalyzer.cxx b/PWGUD/Tasks/DGCandAnalyzer.cxx index 231b54a25c9..a1e76cfb5d4 100644 --- a/PWGUD/Tasks/DGCandAnalyzer.cxx +++ b/PWGUD/Tasks/DGCandAnalyzer.cxx @@ -260,6 +260,7 @@ struct DGCandAnalyzer { // find track combinations which are compatible with PID cuts auto nIVMs = pidsel.computeIVMs(PVContributors); + LOGF(info, "Number of IVMs %d", nIVMs); // update candCase histogram if (nIVMs > 0) { diff --git a/PWGUD/Tasks/UDTutorial_01.cxx b/PWGUD/Tasks/UDTutorial_01.cxx new file mode 100644 index 00000000000..2ed5a158e14 --- /dev/null +++ b/PWGUD/Tasks/UDTutorial_01.cxx @@ -0,0 +1,147 @@ +// 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. +// +// \brief UD tutorial +// \author Paul Buehler, paul.buehler@oeaw.ac.at +// \since April 2023 + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" + +#include "TVector3.h" +#include "Common/DataModel/PIDResponse.h" +#include "PWGUD/DataModel/UDTables.h" +#include "PWGUD/Core/UDHelpers.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UDTutorial01 { + + // configurables + Configurable verbose{"Verbose", {}, "Additional print outs"}; + ConfigurableAxis ptAxis{"ptAxis", {250, 0.0, 2.5}, "p_T axis"}; + ConfigurableAxis etaAxis{"etaAxis", {300, -1.5, 1.5}, ""}; + ConfigurableAxis sigTPCAxis{"sigTPCAxis", {100, -100.0, 100.0}, ""}; + ConfigurableAxis sigTOFAxis{"sigTOFAxis", {100, -100.0, 100.0}, ""}; + + // initialize histogram registry + HistogramRegistry registry{ + "registry", + {}}; + + void init(InitContext&) + { + // Collision histograms + registry.add("collisions/BC", "Relative BC number; Relative BC; Collisions", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); + registry.add("collisions/multiplicityAll", "Multiplicity of all tracks; Tracks; Tracks", {HistType::kTH1F, {{201, -0.5, 200.5}}}); + registry.add("collisions/multiplicityPVC", "Multiplicity of PV contributors; PV contributors; Tracks", {HistType::kTH1F, {{201, -0.5, 200.5}}}); + + // track histograms + const AxisSpec axispt{ptAxis, "p_{T} axis"}; + const AxisSpec axiseta{etaAxis, "pseudo rapidity axis"}; + registry.add("tracks/QCAll", "Track QC of all tracks; Hit in detector; Tracks", {HistType::kTH1F, {{5, -0.5, 4.5}}}); + registry.add("tracks/QCPVC", "Track QC of PV contributors; Hit in detector; Tracks", {HistType::kTH1F, {{5, -0.5, 4.5}}}); + registry.add("tracks/ptAll", "track pt of all tracks; p_{T} [GeV/c]; Tracks", {HistType::kTH1F, {axispt}}); + registry.add("tracks/ptPVC", "track pt of PV contributors; p_{T} [GeV/c]; Tracks", {HistType::kTH1F, {axispt}}); + registry.add("tracks/etavsptAll", "track eta versus pt of all tracks; eta; p_{T} [GeV/c]; Tracks", {HistType::kTH2F, {axiseta, axispt}}); + registry.add("tracks/etavsptPVC", "track eta versus pt of PV contributors; eta; p_{T} [GeV/c]; Tracks", {HistType::kTH2F, {axiseta, axispt}}); + + const AxisSpec axisp{ptAxis, "momentum axis"}; + const AxisSpec axisTPCsig{sigTPCAxis, "TPC signal"}; + const AxisSpec axisTOFsig{sigTOFAxis, "TOF signal"}; + registry.add("tracks/TPCSignalvspAll", "TPC signal versus track momentum of all tracks; Track momentum [GeV/c]; TPC signal [arb. units]; Tracks", {HistType::kTH2F, {axisp, axisTPCsig}}); + registry.add("tracks/TPCSignalvspPVC", "TPC signal versus track momentum of PV contributors; Track momentum [GeV/c]; TPC signal [arb. units]; Tracks", {HistType::kTH2F, {axisp, axisTPCsig}}); + registry.add("tracks/TOFSignalvspAll", "TOF signal versus track momentum of all tracks; Track momentum [GeV/c]; TOF signal [arb. units]; Tracks", {HistType::kTH2F, {axisp, axisTOFsig}}); + registry.add("tracks/TOFSignalvspPVC", "TOF signal versus track momentum of PV contributors; Track momentum [GeV/c]; TOF signal [arb. units]; Tracks", {HistType::kTH2F, {axisp, axisTOFsig}}); + + // FIT histograms + registry.add("FIT/BBFV0A", "Beam-beam in V0A; BC relative to associated BC; Collisions", {HistType::kTH1F, {{32, -16.5, 15.5}}}); + registry.add("FIT/BBFT0A", "Beam-beam in T0A; BC relative to associated BC; Collisions", {HistType::kTH1F, {{32, -16.5, 15.5}}}); + registry.add("FIT/BBFT0C", "Beam-beam in T0C; BC relative to associated BC; Collisions", {HistType::kTH1F, {{32, -16.5, 15.5}}}); + registry.add("FIT/BBFDDA", "Beam-beam in FDDA; BC relative to associated BC; Collisions", {HistType::kTH1F, {{32, -16.5, 15.5}}}); + registry.add("FIT/BBFDDC", "Beam-beam in FDDA; BC relative to associated BC; Collisions", {HistType::kTH1F, {{32, -16.5, 15.5}}}); + } + + // define data types + using UDCollisionsFull = soa::Join; + using UDCollisionFull = UDCollisionsFull::iterator; + using UDTracksFull = soa::Join; + + void process(UDCollisionFull const& dgcand, UDTracksFull const& dgtracks) + { + if (verbose) { + LOGF(info, " DG candidate %d", dgcand.globalIndex()); + } + + // fill collision histograms + registry.get(HIST("collisions/multiplicityAll"))->Fill(dgtracks.size(), 1.); + // select PV contributors + Partition PVContributors = aod::udtrack::isPVContributor == true; + PVContributors.bindTable(dgtracks); + registry.get(HIST("collisions/multiplicityPVC"))->Fill(PVContributors.size(), 1.); + + // relative BC number + auto bcnum = dgcand.globalBC() % o2::constants::lhc::LHCMaxBunches; + registry.get(HIST("collisions/BC"))->Fill(bcnum, 1.); + + // fill track histograms + if (verbose) { + LOGF(info, " Number of tracks %d", dgtracks.size()); + LOGF(info, " Number of PV contributors %d", PVContributors.size()); + } + for (auto track : dgtracks) { + registry.get(HIST("tracks/QCAll"))->Fill(0., 1.); + registry.get(HIST("tracks/QCAll"))->Fill(1., track.hasITS() * 1.); + registry.get(HIST("tracks/QCAll"))->Fill(2., track.hasTPC() * 1.); + registry.get(HIST("tracks/QCAll"))->Fill(3., track.hasTRD() * 1.); + registry.get(HIST("tracks/QCAll"))->Fill(4., track.hasTOF() * 1.); + + auto vtrk = TVector3(track.px(), track.py(), track.pz()); + registry.get(HIST("tracks/ptAll"))->Fill(track.pt(), 1.); + registry.get(HIST("tracks/etavsptAll"))->Fill(vtrk.Eta(), track.pt(), 1.); + + auto signalTPC = track.tpcSignal() * track.sign(); + registry.get(HIST("tracks/TPCSignalvspAll"))->Fill(vtrk.Mag(), signalTPC, 1.); + auto signalTOF = track.tofSignal() / 1.E3; + registry.get(HIST("tracks/TOFSignalvspAll"))->Fill(vtrk.Mag(), signalTOF, 1.); + + if (track.isPVContributor()) { + registry.get(HIST("tracks/QCPVC"))->Fill(0., 1.); + registry.get(HIST("tracks/QCPVC"))->Fill(1., track.hasITS() * 1.); + registry.get(HIST("tracks/QCPVC"))->Fill(2., track.hasTPC() * 1.); + registry.get(HIST("tracks/QCPVC"))->Fill(3., track.hasTRD() * 1.); + registry.get(HIST("tracks/QCPVC"))->Fill(4., track.hasTOF() * 1.); + registry.get(HIST("tracks/ptPVC"))->Fill(track.pt(), 1.); + registry.get(HIST("tracks/etavsptPVC"))->Fill(vtrk.Eta(), track.pt(), 1.); + registry.get(HIST("tracks/TPCSignalvspPVC"))->Fill(vtrk.Mag(), signalTPC, 1.); + registry.get(HIST("tracks/TOFSignalvspPVC"))->Fill(vtrk.Mag(), signalTOF, 1.); + } + } + + // fill FIT histograms + for (auto bit = 0; bit < 33; bit++) { + registry.get(HIST("FIT/BBFV0A"))->Fill(bit - 16, TESTBIT(dgcand.bbFV0Apf(), bit)); + registry.get(HIST("FIT/BBFT0A"))->Fill(bit - 16, TESTBIT(dgcand.bbFT0Apf(), bit)); + registry.get(HIST("FIT/BBFT0C"))->Fill(bit - 16, TESTBIT(dgcand.bbFT0Cpf(), bit)); + registry.get(HIST("FIT/BBFDDA"))->Fill(bit - 16, TESTBIT(dgcand.bbFDDApf(), bit)); + registry.get(HIST("FIT/BBFDDC"))->Fill(bit - 16, TESTBIT(dgcand.bbFDDCpf(), bit)); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"udtutorial01"}), + }; +} diff --git a/PWGUD/Tasks/UDTutorial_02a.cxx b/PWGUD/Tasks/UDTutorial_02a.cxx new file mode 100644 index 00000000000..75aa83c87a7 --- /dev/null +++ b/PWGUD/Tasks/UDTutorial_02a.cxx @@ -0,0 +1,183 @@ +// 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. +// +// \brief UD tutorial +// \author Paul Buehler, paul.buehler@oeaw.ac.at +// \since April 2023 + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" + +#include "TDatabasePDG.h" +#include "TLorentzVector.h" +#include "Common/DataModel/PIDResponse.h" +#include "PWGUD/DataModel/UDTables.h" +#include "PWGUD/Core/UDHelpers.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UDTutorial02a { + + // configurables + Configurable verbose{"Verbose", {}, "Additional print outs"}; + ConfigurableAxis IVMAxis{"IVMAxis", {350, 0.0, 3.5}, "Invariant mass axis"}; + ConfigurableAxis ptAxis{"ptAxis", {250, 0.0, 2.5}, "p_T axis"}; + ConfigurableAxis nsTPCAxis{"nsTPCAxis", {100, -20.0, 20.0}, "nSigma TPC axis"}; + ConfigurableAxis nsTOFAxis{"nsTOFAxis", {100, -100.0, 100.0}, "nSigma TOF axis"}; + + // a pdg object + TDatabasePDG* pdg = nullptr; + + // initialize histogram registry + HistogramRegistry registry{ + "registry", + {}}; + + void init(InitContext&) + { + pdg = TDatabasePDG::Instance(); + + // dgcandidates histograms + const AxisSpec axisIVM{IVMAxis, "IVM axis"}; + const AxisSpec axispt{ptAxis, "pt axis"}; + registry.add("dgcandidates/IVMptsys", "IVM versus system pT; Invariant mass [GeV/c^{2}]; p_{T, system} {GeV/c]", {HistType::kTH2F, {axisIVM, axispt}}); + registry.add("dgcandidates/IVMpttrk", "IVM versus track pT; Invariant mass [GeV/c^{2}]; p_{T, track} {GeV/c]", {HistType::kTH2F, {axisIVM, axispt}}); + + const AxisSpec axisnsTPC{nsTOFAxis, "nSigma TPC axis"}; + registry.add("dgcandidates/nSigmaTPCEl", "TPC nSigma electrons versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{e}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCPi", "TPC nSigma pions versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{pi}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCMu", "TPC nSigma muons versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{mu}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCKa", "TPC nSigma kaon versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{K}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCPr", "TPC nSigma protons versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{p}", {HistType::kTH2F, {axispt, axisnsTPC}}); + + const AxisSpec axisnsTOF{nsTOFAxis, "nSigma TOF axis"}; + registry.add("dgcandidates/nSigmaTOFEl", "TOF nSigma electrons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{e}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFPi", "TOF nSigma pions versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{pi}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFMu", "TOF nSigma muons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{mu}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFKa", "TOF nSigma kaons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{K}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFPr", "TOF nSigma protons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{p}", {HistType::kTH2F, {axispt, axisnsTOF}}); + } + + using UDCollisionsFull = soa::Join; + using UDCollisionFull = UDCollisionsFull::iterator; + using UDTracksFull = soa::Join; + + void process(UDCollisionFull const& dgcand, UDTracksFull const& dgtracks) + { + + // skip events with too few/many tracks + Partition PVContributors = aod::udtrack::isPVContributor == true; + PVContributors.bindTable(dgtracks); + if (dgcand.numContrib() != 2) { + if (verbose) { + LOGF(info, " Candidate rejected: Number of PV contributors is %d", dgcand.numContrib()); + } + return; + } + + // skip events with net charge != 0 + if (dgcand.netCharge() != 0) { + if (verbose) { + LOGF(info, " Candidate rejected: Net charge is %d", dgcand.netCharge()); + } + return; + } + + // skip events with out-of-range rgtrwTOF (fraction-of-good-tracks-with-TOF-hit) + auto rtrwTOF = udhelpers::rPVtrwTOF(dgtracks, PVContributors.size()); + if (rtrwTOF < 0.5) { + if (verbose) { + LOGF(debug, " Candidate rejected: rtrwTOF is %f", rtrwTOF); + } + return; + } + + // check FIT information + auto bitMin = -1 + 16; + auto bitMax = 1 + 16; + for (auto bit = bitMin; bit <= bitMax; bit++) { + if (TESTBIT(dgcand.bbFT0Apf(), bit) || + TESTBIT(dgcand.bbFT0Cpf(), bit) || + TESTBIT(dgcand.bbFV0Apf(), bit) || + TESTBIT(dgcand.bbFDDApf(), bit) || + TESTBIT(dgcand.bbFDDCpf(), bit)) { + return; + } + } + + // check PID of tracks, use nSigmaTPC + // cut on track pT + for (auto trk : PVContributors) { + if (trk.tpcNSigmaPi() < -3. || trk.tpcNSigmaPi() > 3.) { + if (verbose) { + LOGF(info, " Candidate rejected: nSigmaTPC pion is %f", trk.tpcNSigmaPi()); + } + return; + } + if (trk.pt() < 0.1) { + if (verbose) { + LOGF(info, " Candidate rejected: Track pT is %f", trk.pt()); + } + return; + } + } + + // compute invariant mass + TParticlePDG* pion = pdg->GetParticle(211); + TLorentzVector lvtmp; + auto ivm = TLorentzVector(0., 0., 0., 0.); + for (auto trk : PVContributors) { + lvtmp.SetXYZM(trk.px(), trk.py(), trk.pz(), pion->Mass()); + ivm += lvtmp; + } + + // cut on system pT + if (ivm.Perp() < 0.1) { + if (verbose) { + LOGF(info, " Candidate rejected: System pT is %f", ivm.Perp()); + } + return; + } + + // update histograms + if (verbose) { + LOGF(info, " Candidate accepted!"); + } + registry.get(HIST("dgcandidates/IVMptsys"))->Fill(ivm.M(), ivm.Perp()); + for (auto trk : PVContributors) { + registry.get(HIST("dgcandidates/IVMpttrk"))->Fill(ivm.M(), trk.pt()); + + // fill nSigma histograms + auto mom = sqrt(pow(trk.px(), 2) + pow(trk.py(), 2) + pow(trk.pz(), 2)); + registry.get(HIST("dgcandidates/nSigmaTPCEl"))->Fill(mom, trk.tpcNSigmaEl()); + registry.get(HIST("dgcandidates/nSigmaTPCPi"))->Fill(mom, trk.tpcNSigmaPi()); + registry.get(HIST("dgcandidates/nSigmaTPCMu"))->Fill(mom, trk.tpcNSigmaMu()); + registry.get(HIST("dgcandidates/nSigmaTPCKa"))->Fill(mom, trk.tpcNSigmaKa()); + registry.get(HIST("dgcandidates/nSigmaTPCPr"))->Fill(mom, trk.tpcNSigmaPr()); + if (trk.hasTOF()) { + registry.get(HIST("dgcandidates/nSigmaTOFEl"))->Fill(mom, trk.tofNSigmaEl()); + registry.get(HIST("dgcandidates/nSigmaTOFPi"))->Fill(mom, trk.tofNSigmaPi()); + registry.get(HIST("dgcandidates/nSigmaTOFMu"))->Fill(mom, trk.tofNSigmaMu()); + registry.get(HIST("dgcandidates/nSigmaTOFKa"))->Fill(mom, trk.tofNSigmaKa()); + registry.get(HIST("dgcandidates/nSigmaTOFPr"))->Fill(mom, trk.tofNSigmaPr()); + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"udtutorial02a"}), + }; +} diff --git a/PWGUD/Tasks/UDTutorial_02b.cxx b/PWGUD/Tasks/UDTutorial_02b.cxx new file mode 100644 index 00000000000..57b4fc4b31b --- /dev/null +++ b/PWGUD/Tasks/UDTutorial_02b.cxx @@ -0,0 +1,173 @@ +// 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. +// +// \brief UD tutorial +// \author Paul Buehler, paul.buehler@oeaw.ac.at +// \since April 2023 + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" + +#include "Common/DataModel/PIDResponse.h" +#include "PWGUD/DataModel/UDTables.h" +#include "PWGUD/Core/UDHelpers.h" +#include "PWGUD/Core/DGPIDSelector.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UDTutorial02b { + + // configurables + Configurable verbose{"Verbose", {}, "Additional print outs"}; + ConfigurableAxis IVMAxis{"IVMAxis", {350, 0.0, 3.5}, "Invariant mass axis"}; + ConfigurableAxis ptAxis{"ptAxis", {250, 0.0, 2.5}, "p_T axis"}; + ConfigurableAxis nsTPCAxis{"nsTPCAxis", {100, -20.0, 20.0}, "nSigma TPC axis"}; + ConfigurableAxis nsTOFAxis{"nsTOFAxis", {100, -100.0, 100.0}, "nSigma TOF axis"}; + + // analysis cuts + DGAnaparHolder anaPars = DGAnaparHolder(); + Configurable DGPars{"anaPars", {}, "Analysis parameters"}; + + // PID selector + DGPIDSelector pidsel = DGPIDSelector(); + + // initialize histogram registry + HistogramRegistry registry{ + "registry", + {}}; + + void init(InitContext&) + { + anaPars = (DGAnaparHolder)DGPars; + pidsel.init(anaPars); + if (verbose) { + LOGF(info, ""); + pidsel.Print(); + } + + // dgcandidates histograms + const AxisSpec axisIVM{IVMAxis, "IVM axis"}; + const AxisSpec axispt{ptAxis, "pt axis"}; + registry.add("dgcandidates/IVMptsys", "IVM versus system pT; Invariant mass [GeV/c^{2}]; p_{T, system} {GeV/c]", {HistType::kTH2F, {axisIVM, axispt}}); + registry.add("dgcandidates/IVMpttrk", "IVM versus track pT; Invariant mass [GeV/c^{2}]; p_{T, track} {GeV/c]", {HistType::kTH2F, {axisIVM, axispt}}); + + const AxisSpec axisnsTPC{nsTOFAxis, "nSigma TPC axis"}; + registry.add("dgcandidates/nSigmaTPCEl", "TPC nSigma electrons versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{e}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCPi", "TPC nSigma pions versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{pi}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCMu", "TPC nSigma muons versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{mu}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCKa", "TPC nSigma kaon versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{K}", {HistType::kTH2F, {axispt, axisnsTPC}}); + registry.add("dgcandidates/nSigmaTPCPr", "TPC nSigma protons versus track pT; Track p_{T} [GeV/c]; TPC nSigma_{p}", {HistType::kTH2F, {axispt, axisnsTPC}}); + + const AxisSpec axisnsTOF{nsTOFAxis, "nSigma TOF axis"}; + registry.add("dgcandidates/nSigmaTOFEl", "TOF nSigma electrons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{e}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFPi", "TOF nSigma pions versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{pi}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFMu", "TOF nSigma muons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{mu}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFKa", "TOF nSigma kaons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{K}", {HistType::kTH2F, {axispt, axisnsTOF}}); + registry.add("dgcandidates/nSigmaTOFPr", "TOF nSigma protons versus track pT; Track p_{T} [GeV/c]; TOF nSigma_{p}", {HistType::kTH2F, {axispt, axisnsTOF}}); + } + + using UDCollisionsFull = soa::Join; + using UDCollisionFull = UDCollisionsFull::iterator; + using UDTracksFull = soa::Join; + + void process(UDCollisionFull const& dgcand, UDTracksFull const& dgtracks) + { + + // skip events with too few/many tracks + Partition PVContributors = aod::udtrack::isPVContributor == true; + PVContributors.bindTable(dgtracks); + if (dgcand.numContrib() < anaPars.minNTracks() || dgcand.numContrib() > anaPars.maxNTracks()) { + if (verbose) { + LOGF(info, " Candidate rejected: Number of PV contributors %d not in range [%d, %d].", dgcand.numContrib(), anaPars.minNTracks(), anaPars.maxNTracks()); + } + return; + } + + // skip events with out-of-range net charge + auto netChargeValues = anaPars.netCharges(); + if (std::find(netChargeValues.begin(), netChargeValues.end(), dgcand.netCharge()) == netChargeValues.end()) { + if (verbose) { + LOGF(info, " Candidate rejected: Net charge %d not in selected set.", dgcand.netCharge()); + } + return; + } + + // skip events with out-of-range rgtrwTOF (fraction-of-good-tracks-with-TOF-hit) + auto rtrwTOF = udhelpers::rPVtrwTOF(dgtracks, PVContributors.size()); + if (rtrwTOF < anaPars.minRgtrwTOF()) { + if (verbose) { + LOGF(debug, " Candidate rejected: rtrwTOF %f below threshold of %f.", rtrwTOF, anaPars.minRgtrwTOF()); + } + return; + } + + // check FIT information + auto bitMin = anaPars.dBCMin() + 16; + auto bitMax = anaPars.dBCMax() + 16; + for (auto bit = bitMin; bit <= bitMax; bit++) { + if (TESTBIT(dgcand.bbFT0Apf(), bit) || + TESTBIT(dgcand.bbFT0Cpf(), bit) || + TESTBIT(dgcand.bbFV0Apf(), bit) || + TESTBIT(dgcand.bbFDDApf(), bit) || + TESTBIT(dgcand.bbFDDCpf(), bit)) { + return; + } + } + + // find track combinations which are compatible with PID cuts + auto nIVMs = pidsel.computeIVMs(PVContributors); + if (verbose) { + LOGF(info, "< Found %d candidates", nIVMs); + } + + // update histograms + for (auto ivm : pidsel.IVMs()) { + // cut on pt-system + if (ivm.Perp() < anaPars.minptsys() || ivm.Perp() > anaPars.maxptsys()) { + continue; + } + + if (verbose) { + LOGF(info, " Candidate accepted!"); + } + registry.get(HIST("dgcandidates/IVMptsys"))->Fill(ivm.M(), ivm.Perp()); + for (auto ind : ivm.trkinds()) { + auto track = PVContributors.begin() + ind; + registry.get(HIST("dgcandidates/IVMpttrk"))->Fill(ivm.M(), track.pt()); + + // fill nSigma histograms + auto mom = sqrt(pow(track.px(), 2) + pow(track.py(), 2) + pow(track.pz(), 2)); + registry.get(HIST("dgcandidates/nSigmaTPCEl"))->Fill(mom, track.tpcNSigmaEl()); + registry.get(HIST("dgcandidates/nSigmaTPCPi"))->Fill(mom, track.tpcNSigmaPi()); + registry.get(HIST("dgcandidates/nSigmaTPCMu"))->Fill(mom, track.tpcNSigmaMu()); + registry.get(HIST("dgcandidates/nSigmaTPCKa"))->Fill(mom, track.tpcNSigmaKa()); + registry.get(HIST("dgcandidates/nSigmaTPCPr"))->Fill(mom, track.tpcNSigmaPr()); + if (track.hasTOF()) { + LOGF(debug, "tofNSigmaPi %f", track.tofNSigmaPi()); + registry.get(HIST("dgcandidates/nSigmaTOFEl"))->Fill(mom, track.tofNSigmaEl()); + registry.get(HIST("dgcandidates/nSigmaTOFPi"))->Fill(mom, track.tofNSigmaPi()); + registry.get(HIST("dgcandidates/nSigmaTOFMu"))->Fill(mom, track.tofNSigmaMu()); + registry.get(HIST("dgcandidates/nSigmaTOFKa"))->Fill(mom, track.tofNSigmaKa()); + registry.get(HIST("dgcandidates/nSigmaTOFPr"))->Fill(mom, track.tofNSigmaPr()); + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"udtutorial02b"}), + }; +} diff --git a/PWGUD/Tasks/diffMCQA.cxx b/PWGUD/Tasks/diffMCQA.cxx index c5fe025fea5..a54b0c29d81 100644 --- a/PWGUD/Tasks/diffMCQA.cxx +++ b/PWGUD/Tasks/diffMCQA.cxx @@ -10,36 +10,6 @@ // or submit itself to any jurisdiction. /// /// \brief A QA task for DG events -/// -/// options: -/// DiffCuts.mNDtcoll(4) -/// DiffCuts.mMinNBCs(7) -/// DiffCuts.mMinNTracks(0) -/// DiffCuts.mMaxNTracks(10000) -/// DiffCuts.mNetCharges({0}) -/// DiffCuts.mPidHypo(211) -/// DiffCuts.mMinPosz(-1000.) -/// DiffCuts.mMaxPosz(1000.) -/// DiffCuts.mMinPt(0.) -/// DiffCuts.mMaxPt(1000.) -/// DiffCuts.mMinEta(-1.) -/// DiffCuts.mMaxEta(1.) -/// DiffCuts.mMinIVM(0.) -/// DiffCuts.mMaxIVM(1000.) -/// DiffCuts.mMaxnSigmaTPC(1000.) -/// DiffCuts.mMaxnSigmaTOF(1000.) -/// DiffCutsX.mFITAmpLimits({0., 0., 0., 0., 0.}) -/// -/// usage: copts="--configuration json://DiffQAConfig.json -b" -/// -/// o2-analysis-timestamp $copts | -/// o2-analysis-track-propagation $copts | -/// o2-analysis-event-selection $copts | -/// o2-analysis-ft0-corrected-table $copts | -/// o2-analysis-trackextension $copts | -/// o2-analysis-trackselection $copts | -/// o2-analysis-ud-diff-mcqa $copts > diffQA.log -/// /// \author Paul Buehler, paul.buehler@oeaw.ac.at /// \since 01.10.2021 @@ -73,102 +43,17 @@ struct DiffMCQA { o2::dataformats::bcRanges abcrs = o2::dataformats::bcRanges("ambiguous_tracks"); o2::dataformats::bcRanges afbcrs = o2::dataformats::bcRanges("ambiguous_fwdtracks"); - // define histograms - // Stat: - // bin 0: all collisions - // bin 1: clean FIT - // bin 2: clean ZDC - // bin 3: clean Calo - // bin 4: no V0s - // bin 5: no Cascades - // bin 6: no FWD tracks - // bin 7: no global tracks which are no vtx tracks - // bin 8: no vtx tracks which are no global tracks - // bin 9: possible ambiguous tracks - // bin 10: possible ambiguous FwdTracks - // bin 11: fraction of PV tracks with TOF hit > minRgtrwTOF - // bin 12: number of tracks >= minimum number - // bin 13: number of tracks <= maximum number - // bin 14: minimum pt <= pt of vtx tracks <= maximum pt - // bin 15: minimum eta <= eta of vtx tracks <= maximum eta - // bin 16: net charge >= minimum net charge - // bin 17: net charge <= maximum net charge - // bin 18: IVM >= minimum IVM - // bin 19: IVM <= maximum IVM - // - // 3 diverent versions of histograms: - // Diff1: Pythia MBR - // Diff2: GRANIITTI - // : Rest (Pythia MB) + // 3 different versions of histograms: + // MBRDiff: Pythia MBR + // Diff : GRANIITTI + // nonDiff: Rest (Pythia MB) // + // initialize HistogramRegistry HistogramRegistry registry{ "registry", - { - // non diffractive events - {"Stat", "#Stat", {HistType::kTH1F, {{20, -0.5, 19.5}}}}, - {"cleanFIT", "#cleanFIT", {HistType::kTH2F, {{10, -0.5, 9.5}, {2, -0.5, 1.5}}}}, - {"Tracks", "#Tracks", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"vtxTracks", "#vtxTracks", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"globalTracks", "#globalTracks", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"rejectedTracks", "#rejectedTracks", {HistType::kTH1F, {{17, -0.5, 16.5}}}}, - {"tResvsrTOFTracks", "#tResvsrTOFTracks", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}}, - {"vtxPosxy", "#vtxPosxy", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}}, - {"vtxPosz", "#vtxPosz", {HistType::kTH1F, {{1000, -100., 100.}}}}, - {"etapt", "#etapt", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}}, - {"dEdxTPC", "#dEdxTPC", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}}, - {"dEdxTOF", "#dEdxTOF", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}}, - {"vtxPosxyDG", "#vtxPosxyDG", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}}, - {"vtxPoszDG", "#vtxPoszDG", {HistType::kTH1F, {{1000, -100., 100.}}}}, - {"etaptDG", "#etaptDG", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}}, - {"dEdxTPCDG", "#dEdxTPCDG", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}}, - {"dEdxTOFDG", "#dEdxTOFDG", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}}, - {"netChargeDG", "#netChargeDG", {HistType::kTH1F, {{21, -10.5, 10.5}}}}, - {"IVMptSysDG", "#IVMptSysDG", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}}, - {"IVMptTrkDG", "#IVMptTrkDG", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}}, - // PYTHIA8 diffractive events - {"StatDiff1", "#StatDiff1", {HistType::kTH1F, {{20, -0.5, 19.5}}}}, - {"cleanFITDiff1", "#cleanFITDiff1", {HistType::kTH2F, {{10, -0.5, 9.5}, {2, -0.5, 1.5}}}}, - {"TracksDiff1", "#TracksDiff1", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"vtxTracksDiff1", "#vtxTracksDiff1", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"globalTracksDiff1", "#globalTracksDiff1", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"rejectedTracksDiff1", "#rejectedTracksDiff1", {HistType::kTH1F, {{17, -0.5, 16.5}}}}, - {"tResvsrTOFTracksDiff1", "#tResvsrTOFTracksDiff1", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}}, - {"vtxPosxyDiff1", "#vtxPosxyDiff1", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}}, - {"vtxPoszDiff1", "#vtxPoszDiff1", {HistType::kTH1F, {{1000, -100., 100.}}}}, - {"etaptDiff1", "#etaptDiff1", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}}, - {"dEdxTPCDiff1", "#dEdxTPCDiff1", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}}, - {"dEdxTOFDiff1", "#dEdxTOFDiff1", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}}, - {"vtxPosxyDGDiff1", "#vtxPosxyDGDiff1", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}}, - {"vtxPoszDGDiff1", "#vtxPoszDGDiff1", {HistType::kTH1F, {{1000, -100., 100.}}}}, - {"etaptDGDiff1", "#etaptDGDiff1", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}}, - {"dEdxTPCDGDiff1", "#dEdxTPCDGDiff1", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}}, - {"dEdxTOFDGDiff1", "#dEdxTOFDGDiff1", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}}, - {"netChargeDGDiff1", "#netChargeDGDiff1", {HistType::kTH1F, {{21, -10.5, 10.5}}}}, - {"IVMptSysDGDiff1", "#IVMptSysDGDiff1", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}}, - {"IVMptTrkDGDiff1", "#IVMptTrkDGDiff1", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}}, - // GRANIITTI diffractive events - {"StatDiff2", "#StatDiff2", {HistType::kTH1F, {{20, -0.5, 19.5}}}}, - {"cleanFITDiff2", "#cleanFITDiff2", {HistType::kTH2F, {{10, -0.5, 9.5}, {2, -0.5, 1.5}}}}, - {"TracksDiff2", "#TracksDiff2", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"vtxTracksDiff2", "#vtxTracksDiff2", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"globalTracksDiff2", "#globalTracksDiff2", {HistType::kTH1F, {{50, 0.5, 50.5}}}}, - {"rejectedTracksDiff2", "#rejectedTracksDiff2", {HistType::kTH1F, {{17, -0.5, 16.5}}}}, - {"tResvsrTOFTracksDiff2", "#tResvsrTOFTracksDiff2", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}}, - {"vtxPosxyDiff2", "#vtxPosxyDiff2", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}}, - {"vtxPoszDiff2", "#vtxPoszDiff2", {HistType::kTH1F, {{1000, -100., 100.}}}}, - {"etaptDiff2", "#etaptDiff2", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}}, - {"dEdxTPCDiff2", "#dEdxTPCDiff2", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}}, - {"dEdxTOFDiff2", "#dEdxTOFDiff2", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}}, - {"vtxPosxyDGDiff2", "#vtxPosxyDGDiff2", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}}, - {"vtxPoszDGDiff2", "#vtxPoszDGDiff2", {HistType::kTH1F, {{1000, -100., 100.}}}}, - {"etaptDGDiff2", "#etaptDGDiff2", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}}, - {"dEdxTPCDGDiff2", "#dEdxTPCDGDiff2", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}}, - {"dEdxTOFDGDiff2", "#dEdxTOFDGDiff2", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}}, - {"netChargeDGDiff2", "#netChargeDGDiff2", {HistType::kTH1F, {{21, -10.5, 10.5}}}}, - {"IVMptSysDGDiff2", "#IVMptSysDGDiff2", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}}, - {"IVMptTrkDGDiff2", "#IVMptTrkDGDiff2", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}}, - }}; + {}}; + // define abbreviations using CCs = soa::Join; using CC = CCs::iterator; using BCs = soa::Join; @@ -177,11 +62,105 @@ struct DiffMCQA { using ATs = aod::AmbiguousTracks; using AFTs = aod::AmbiguousFwdTracks; - void init(InitContext&) + void init(InitContext& context) { maxdEdxTPC = 0.; maxdEdxTOF = 0.; diffCuts = (DGCutparHolder)DGCuts; + + // add histograms for the different process functions + if (context.mOptions.get("processMCTruth")) { + registry.add("MCTruth/Stat", "Simulated event type; event type; Entries", {HistType::kTH1F, {{4, 0.5, 4.5}}}); + } + + if (context.mOptions.get("processMain")) { + // non diffractive events + registry.add("nonDiff/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{20, -0.5, 19.5}}}); + registry.add("nonDiff/cleanFIT", "Statistics of collisions with empty FIT; Multiple of collision time resolution; FIT status; Collisions", {HistType::kTH2F, {{10, -0.5, 9.5}, {2, -0.5, 1.5}}}); + registry.add("nonDiff/Tracks", "Number of tracks; Number of tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("nonDiff/PVTracks", "Number of PV tracks; Number of PV tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("nonDiff/globalTracks", "Number of global tracks; Number of global tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("nonDiff/notPVTracks", "Not PV tracks; Track status bit; Not PV tracks", {HistType::kTH1F, {{17, -0.5, 16.5}}}); + registry.add("nonDiff/tResvsrTOFTracks", "#tResvsrTOFTracks", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}); + registry.add("nonDiff/PVposxy", "Vertex position in x and y direction; V_x; V_y; Collisions", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); + registry.add("nonDiff/PVposz", "Vertex position in z direction; V_z; Collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("nonDiff/etapt", "eta versus pT of all tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("nonDiff/dEdxTPC", "TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}); + registry.add("nonDiff/dEdxTOF", "TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); + registry.add("nonDiff/PVposxyDG", "DG: Vertex position in x and y direction; V_x; V_y; Collisions", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); + registry.add("nonDiff/PVposzDG", "DG: Vertex position in z direction; V_z; Collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("nonDiff/etaptDG", "DG: eta versus pT of all tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("nonDiff/dEdxTPCDG", "DG: TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}); + registry.add("nonDiff/dEdxTOFDG", "DG: TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); + registry.add("nonDiff/netChargeDG", "DG: Net charge; Net charge; DG collisions", {HistType::kTH1F, {{21, -10.5, 10.5}}}); + registry.add("nonDiff/IVMptSysDG", "DG: Invariant mass versus p_{T, system}; Invarian mass [GeV/c^2]; p_{T, system} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + registry.add("nonDiff/IVMptTrkDG", "DG: Invariant mass versus p_{T, tracks}; Invarian mass [GeV/c^2]; p_{T, tracks} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + // PYTHIA8 diffractive events + registry.add("MBRDiff/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{20, -0.5, 19.5}}}); + registry.add("MBRDiff/cleanFIT", "Statistics of collisions with empty FIT; Multiple of collision time resolution; FIT status; Collisions", {HistType::kTH2F, {{10, -0.5, 9.5}, {2, -0.5, 1.5}}}); + registry.add("MBRDiff/Tracks", "Number of tracks; Number of tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("MBRDiff/PVTracks", "Number of PV tracks; Number of PV tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("MBRDiff/globalTracks", "Number of global tracks; Number of global tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("MBRDiff/notPVTracks", "Not PV tracks; Track status bit; Not PV tracks", {HistType::kTH1F, {{17, -0.5, 16.5}}}); + registry.add("MBRDiff/tResvsrTOFTracks", "#tResvsrTOFTracks", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}); + registry.add("MBRDiff/PVposxy", "Vertex position in x and y direction; V_x; V_y; Collisions", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); + registry.add("MBRDiff/PVposz", "Vertex position in z direction; V_z; Collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("MBRDiff/etapt", "eta versus pT of all tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("MBRDiff/dEdxTPC", "TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}); + registry.add("MBRDiff/dEdxTOF", "TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); + registry.add("MBRDiff/PVposxyDG", "DG: Vertex position in x and y direction; V_x; V_y; Collisions", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); + registry.add("MBRDiff/PVposzDG", "DG: Vertex position in z direction; V_z; Collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("MBRDiff/etaptDG", "DG: eta versus pT of all tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("MBRDiff/dEdxTPCDG", "DG: TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}); + registry.add("MBRDiff/dEdxTOFDG", "DG: TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); + registry.add("MBRDiff/netChargeDG", "DG: Net charge; Net charge; DG collisions", {HistType::kTH1F, {{21, -10.5, 10.5}}}); + registry.add("MBRDiff/IVMptSysDG", "DG: Invariant mass versus p_{T, system}; Invarian mass [GeV/c^2]; p_{T, system} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + registry.add("MBRDiff/IVMptTrkDG", "DG: Invariant mass versus p_{T, tracks}; Invarian mass [GeV/c^2]; p_{T, tracks} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + // GRANIITTI diffractive events + registry.add("Diff/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{20, -0.5, 19.5}}}); + registry.add("Diff/cleanFIT", "Statistics of collisions with empty FIT; Multiple of collision time resolution; FIT status; Collisions", {HistType::kTH2F, {{10, -0.5, 9.5}, {2, -0.5, 1.5}}}); + registry.add("Diff/Tracks", "Number of tracks; Number of tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("Diff/PVTracks", "Number of PV tracks; Number of PV tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("Diff/globalTracks", "Number of global tracks; Number of global tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("Diff/notPVTracks", "Not PV tracks; Track status bit; Not PV tracks", {HistType::kTH1F, {{17, -0.5, 16.5}}}); + registry.add("Diff/tResvsrTOFTracks", "#tResvsrTOFTracks", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}); + registry.add("Diff/PVposxy", "Vertex position in x and y direction; V_x; V_y; Collisions", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); + registry.add("Diff/PVposz", "Vertex position in z direction; V_z; Collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("Diff/etapt", "eta versus pT of all tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("Diff/dEdxTPC", "TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}); + registry.add("Diff/dEdxTOF", "TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); + registry.add("Diff/PVposxyDG", "DG: Vertex position in x and y direction; V_x; V_y; Collisions", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); + registry.add("Diff/PVposzDG", "DG: Vertex position in z direction; V_z; Collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("Diff/etaptDG", "DG: eta versus pT of all tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("Diff/dEdxTPCDG", "DG: TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {3000, 0., 30000.}}}); + registry.add("Diff/dEdxTOFDG", "DG: TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); + registry.add("Diff/netChargeDG", "DG: Net charge; Net charge; DG collisions", {HistType::kTH1F, {{21, -10.5, 10.5}}}); + registry.add("Diff/IVMptSysDG", "DG: Invariant mass versus p_{T, system}; Invarian mass [GeV/c^2]; p_{T, system} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + registry.add("Diff/IVMptTrkDG", "DG: Invariant mass versus p_{T, tracks}; Invarian mass [GeV/c^2]; p_{T, tracks} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + } + + if (context.mOptions.get("processFV0")) { + registry.add("FV0/FV0Aamp", "Amplitudes of FV0A channels; FV0A amplitude; FV0A channel; Entries", {HistType::kTH2F, {{48, -0.5, 47.5}, {1000, 0., 1000.}}}); + } + + if (context.mOptions.get("processFT0")) { + registry.add("FT0/FT0Aamp", "Amplitudes of FT0A channels; FT0A amplitude; FT0A channel; Entries", {HistType::kTH2F, {{96, -0.5, 95.5}, {100, 0., 200.}}}); + registry.add("FT0/FT0Camp", "Amplitudes of FT0C channels; FT0C amplitude; FT0C channel; Entries", {HistType::kTH2F, {{112, -0.5, 111.5}, {100, 0., 200.}}}); + } + + if (context.mOptions.get("processFDD")) { + registry.add("FDD/FDDAamp", "Amplitudes of FDDA channels; FDDA amplitude; FDDA channel; Entries", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}); + registry.add("FDD/FDDCamp", "Amplitudes of FDDC channels; FDDC amplitude; FDDC channel; Entries", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}); + } + + if (context.mOptions.get("processZDC")) { + registry.add("ZdcEnergies", "Various ZDC energies; Energy type; ZDC energy; Entries", {HistType::kTH2F, {{22, -0.5, 21.5}, {100, 0., 1000.}}}); + } + + if (context.mOptions.get("processCalo")) { + registry.add("CaloCell", "Calorimeter cell; Calorimeter cell number; Entries", {HistType::kTH1I, {{18000, -0.5, 17999.5}}}); + registry.add("CaloAmplitude", "Calorimeter amplitudes; Calorimeter amplitude; Entries", {HistType::kTH1F, {{100, 0, 10.}}}); + } } void run(ProcessingContext& pc) @@ -230,12 +209,38 @@ struct DiffMCQA { Preslice perMcCollision = aod::mcparticle::mcCollisionId; - void process(CC const& collision, BCs const& bct0s, - TCs const& tracks, FWs const& fwdtracks, ATs const& ambtracks, AFTs const& ambfwdtracks, - aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds, - aod::Zdcs& zdcs, aod::Calos& calos, - aod::V0s const& v0s, aod::Cascades const& cascades, - aod::McCollisions const& McCols, aod::McParticles const& McParts) + // ............................................................................................................... + void processMCTruth(CCs const& collisions, aod::McCollisions const& mccollisions, aod::McParticles const& McParts) + { + LOGF(info, "Number of mc collisions %d collisions %d", mccollisions.size(), collisions.size()); + + // is this a central diffractive event? + // by default it is assumed to be a MB event + bool isPythiaDiff = false; + bool isGraniittiDiff = false; + for (auto mccol : mccollisions) { + registry.get(HIST("MCTruth/Stat"))->Fill(1., 1.); + auto MCPartSlice = McParts.sliceBy(perMcCollision, mccol.globalIndex()); + isPythiaDiff = udhelpers::isPythiaCDE(MCPartSlice); + isGraniittiDiff = udhelpers::isGraniittiCDE(MCPartSlice); + if (isPythiaDiff) { + registry.get(HIST("MCTruth/Stat"))->Fill(3., 1.); + } else if (isGraniittiDiff) { + registry.get(HIST("MCTruth/Stat"))->Fill(4., 1.); + } else { + registry.get(HIST("MCTruth/Stat"))->Fill(2., 1.); + } + } + } + PROCESS_SWITCH(DiffMCQA, processMCTruth, "Process NC truth", true); + + // ............................................................................................................... + void processMain(CC const& collision, BCs const& bct0s, + TCs const& tracks, FWs const& fwdtracks, ATs const& ambtracks, AFTs const& ambfwdtracks, + aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds, + aod::Zdcs& zdcs, aod::Calos& calos, + aod::V0s const& v0s, aod::Cascades const& cascades, + aod::McCollisions const& McCols, aod::McParticles const& McParts) { bool isDGcandidate = true; @@ -256,26 +261,26 @@ struct DiffMCQA { // update collision histograms if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(0., 1.); - registry.get(HIST("vtxPosxyDiff1"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPoszDiff1"))->Fill(collision.posZ()); - registry.get(HIST("TracksDiff1"))->Fill(tracks.size()); - registry.get(HIST("vtxTracksDiff1"))->Fill(collision.numContrib()); - registry.get(HIST("globalTracksDiff1"))->Fill(goodTracks.size()); + registry.get(HIST("MBRDiff/Stat"))->Fill(0., 1.); + registry.get(HIST("MBRDiff/PVposxy"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("MBRDiff/PVposz"))->Fill(collision.posZ()); + registry.get(HIST("MBRDiff/Tracks"))->Fill(tracks.size()); + registry.get(HIST("MBRDiff/PVTracks"))->Fill(collision.numContrib()); + registry.get(HIST("MBRDiff/globalTracks"))->Fill(goodTracks.size()); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(0., 1.); - registry.get(HIST("vtxPosxyDiff2"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPoszDiff2"))->Fill(collision.posZ()); - registry.get(HIST("TracksDiff2"))->Fill(tracks.size()); - registry.get(HIST("vtxTracksDiff2"))->Fill(collision.numContrib()); - registry.get(HIST("globalTracksDiff2"))->Fill(goodTracks.size()); + registry.get(HIST("Diff/Stat"))->Fill(0., 1.); + registry.get(HIST("Diff/PVposxy"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("Diff/PVposz"))->Fill(collision.posZ()); + registry.get(HIST("Diff/Tracks"))->Fill(tracks.size()); + registry.get(HIST("Diff/PVTracks"))->Fill(collision.numContrib()); + registry.get(HIST("Diff/globalTracks"))->Fill(goodTracks.size()); } else { - registry.get(HIST("Stat"))->Fill(0., 1.); - registry.get(HIST("vtxPosxy"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPosz"))->Fill(collision.posZ()); - registry.get(HIST("Tracks"))->Fill(tracks.size()); - registry.get(HIST("vtxTracks"))->Fill(collision.numContrib()); - registry.get(HIST("globalTracks"))->Fill(goodTracks.size()); + registry.get(HIST("nonDiff/Stat"))->Fill(0., 1.); + registry.get(HIST("nonDiff/PVposxy"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("nonDiff/PVposz"))->Fill(collision.posZ()); + registry.get(HIST("nonDiff/Tracks"))->Fill(tracks.size()); + registry.get(HIST("nonDiff/PVTracks"))->Fill(collision.numContrib()); + registry.get(HIST("nonDiff/globalTracks"))->Fill(goodTracks.size()); } // test influence of BCrange width @@ -286,27 +291,27 @@ struct DiffMCQA { isDGcandidate &= udhelpers::cleanFIT(bc, diffCuts.FITAmpLimits()); } if (isPythiaDiff) { - registry.get(HIST("cleanFITDiff1"))->Fill(NDtcoll, isDGcandidate * 1.); + registry.get(HIST("MBRDiff/cleanFIT"))->Fill(NDtcoll, isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("cleanFITDiff2"))->Fill(NDtcoll, isDGcandidate * 1.); + registry.get(HIST("Diff/cleanFIT"))->Fill(NDtcoll, isDGcandidate * 1.); } else { - registry.get(HIST("cleanFIT"))->Fill(NDtcoll, isDGcandidate * 1.); + registry.get(HIST("nonDiff/cleanFIT"))->Fill(NDtcoll, isDGcandidate * 1.); } } - // number of vertex tracks with TOF hit + // loop over all tracks float rgtrwTOF = 0.; for (auto const& track : tracks) { // update eta vs pt and dEdx histograms histogram if (isPythiaDiff) { - registry.get(HIST("etaptDiff1"))->Fill(track.eta(), track.pt()); - registry.get(HIST("dEdxTPCDiff1"))->Fill(track.p(), track.tpcSignal()); + registry.get(HIST("MBRDiff/etapt"))->Fill(track.eta(), track.pt()); + registry.get(HIST("MBRDiff/dEdxTPC"))->Fill(track.p(), track.tpcSignal()); } else if (isGraniittiDiff) { - registry.get(HIST("etaptDiff2"))->Fill(track.eta(), track.pt()); - registry.get(HIST("dEdxTPCDiff2"))->Fill(track.p(), track.tpcSignal()); + registry.get(HIST("Diff/etapt"))->Fill(track.eta(), track.pt()); + registry.get(HIST("Diff/dEdxTPC"))->Fill(track.p(), track.tpcSignal()); } else { - registry.get(HIST("etapt"))->Fill(track.eta(), track.pt()); - registry.get(HIST("dEdxTPC"))->Fill(track.p(), track.tpcSignal()); + registry.get(HIST("nonDiff/etapt"))->Fill(track.eta(), track.pt()); + registry.get(HIST("nonDiff/dEdxTPC"))->Fill(track.p(), track.tpcSignal()); } if (track.tpcSignal() > maxdEdxTPC) { maxdEdxTPC = track.tpcSignal(); @@ -316,11 +321,11 @@ struct DiffMCQA { // TOF hit? if (track.hasTOF()) { if (isPythiaDiff) { - registry.get(HIST("dEdxTOFDiff1"))->Fill(track.p(), track.tofSignal()); + registry.get(HIST("MBRDiff/dEdxTOF"))->Fill(track.p(), track.tofSignal()); } else if (isGraniittiDiff) { - registry.get(HIST("dEdxTOFDiff2"))->Fill(track.p(), track.tofSignal()); + registry.get(HIST("Diff/dEdxTOF"))->Fill(track.p(), track.tofSignal()); } else { - registry.get(HIST("dEdxTOF"))->Fill(track.p(), track.tofSignal()); + registry.get(HIST("nonDiff/dEdxTOF"))->Fill(track.p(), track.tofSignal()); } if (track.tofSignal() > maxdEdxTOF) { maxdEdxTOF = track.tofSignal(); @@ -338,11 +343,11 @@ struct DiffMCQA { } LOGF(debug, " Vertex tracks with TOF: %f [1]", rgtrwTOF); if (isPythiaDiff) { - registry.get(HIST("tResvsrTOFTracksDiff1"))->Fill(collision.collisionTimeRes(), rgtrwTOF); + registry.get(HIST("MBRDiff/tResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rgtrwTOF); } else if (isGraniittiDiff) { - registry.get(HIST("tResvsrTOFTracksDiff2"))->Fill(collision.collisionTimeRes(), rgtrwTOF); + registry.get(HIST("Diff/tResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rgtrwTOF); } else { - registry.get(HIST("tResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rgtrwTOF); + registry.get(HIST("nonDiff/tResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rgtrwTOF); } // is it a DG candidate? @@ -373,11 +378,11 @@ struct DiffMCQA { } } if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(1., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(1., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(1., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(1., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(1., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(1., isDGcandidate * 1.); } // no Zdc signal in bcSlice @@ -389,11 +394,11 @@ struct DiffMCQA { } } if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(2., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(2., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(2., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(2., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(2., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(2., isDGcandidate * 1.); } // no Calo signal in bcSlice @@ -404,41 +409,41 @@ struct DiffMCQA { } } if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(3., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(3., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(3., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(3., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(3., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(3., isDGcandidate * 1.); } // no V0s isDGcandidate &= (v0s.size() == 0); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(4., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(4., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(4., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(4., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(4., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(4., isDGcandidate * 1.); } // no Cascades isDGcandidate &= (cascades.size() == 0); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(5., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(5., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(5., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(5., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(5., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(5., isDGcandidate * 1.); } // number of forward tracks = 0 isDGcandidate &= (fwdtracks.size() == 0); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(6., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(6., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(6., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(6., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(6., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(6., isDGcandidate * 1.); } // no global tracks which are no vtx tracks @@ -453,14 +458,14 @@ struct DiffMCQA { } } if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(7., globalAndVtx * 1.); - registry.get(HIST("StatDiff1"))->Fill(8., vtxAndGlobal * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(7., globalAndVtx * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(8., vtxAndGlobal * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(7., globalAndVtx * 1.); - registry.get(HIST("StatDiff2"))->Fill(8., vtxAndGlobal * 1.); + registry.get(HIST("Diff/Stat"))->Fill(7., globalAndVtx * 1.); + registry.get(HIST("Diff/Stat"))->Fill(8., vtxAndGlobal * 1.); } else { - registry.get(HIST("Stat"))->Fill(7., globalAndVtx * 1.); - registry.get(HIST("Stat"))->Fill(8., vtxAndGlobal * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(7., globalAndVtx * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(8., vtxAndGlobal * 1.); } isDGcandidate &= globalAndVtx; if (diffCuts.globalTracksOnly()) { @@ -476,11 +481,11 @@ struct DiffMCQA { } } if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(9., noAmbTracks * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(9., noAmbTracks * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(9., noAmbTracks * 1.); + registry.get(HIST("Diff/Stat"))->Fill(9., noAmbTracks * 1.); } else { - registry.get(HIST("Stat"))->Fill(9., noAmbTracks * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(9., noAmbTracks * 1.); } // check a given bc for possible ambiguous FwdTracks @@ -492,39 +497,39 @@ struct DiffMCQA { } } if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(10., noAmbFwdTracks * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(10., noAmbFwdTracks * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(10., noAmbFwdTracks * 1.); + registry.get(HIST("Diff/Stat"))->Fill(10., noAmbFwdTracks * 1.); } else { - registry.get(HIST("Stat"))->Fill(10., noAmbFwdTracks * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(10., noAmbFwdTracks * 1.); } // at least one vtx track with TOF hit isDGcandidate &= (rgtrwTOF >= diffCuts.minRgtrwTOF()); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(11., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(11., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(11., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(11., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(11., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(11., isDGcandidate * 1.); } // number of vertex tracks <= n isDGcandidate &= (collision.numContrib() >= diffCuts.minNTracks()); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(12., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(12., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(12., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(12., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(12., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(12., isDGcandidate * 1.); } isDGcandidate &= (collision.numContrib() <= diffCuts.maxNTracks()); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(13., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(13., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(13., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(13., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(13., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(13., isDGcandidate * 1.); } // net charge and invariant mass @@ -544,62 +549,62 @@ struct DiffMCQA { // check also pt and eta of tracks for (auto& track : tracks) { - // update histogram rejectedTracks + // update histogram notPVTracks if (!track.isPVContributor()) { if (isPythiaDiff) { - registry.get(HIST("rejectedTracksDiff1"))->Fill(0., 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(1., track.isGlobalTrackSDD() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(2., track.passedTrackType() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(3., track.passedPtRange() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(4., track.passedEtaRange() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(5., track.passedTPCNCls() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(6., track.passedTPCCrossedRows() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(8., track.passedTPCChi2NDF() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(9., track.passedTPCRefit() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(10., track.passedITSNCls() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(11., track.passedITSChi2NDF() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(12., track.passedITSRefit() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(13., track.passedITSHits() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(14., track.passedGoldenChi2() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(15., track.passedDCAxy() * 1.); - registry.get(HIST("rejectedTracksDiff1"))->Fill(16., track.passedDCAz() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(0., 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(1., track.isGlobalTrackSDD() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(2., track.passedTrackType() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(3., track.passedPtRange() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(4., track.passedEtaRange() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(5., track.passedTPCNCls() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(6., track.passedTPCCrossedRows() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(8., track.passedTPCChi2NDF() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(9., track.passedTPCRefit() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(10., track.passedITSNCls() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(11., track.passedITSChi2NDF() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(12., track.passedITSRefit() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(13., track.passedITSHits() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(14., track.passedGoldenChi2() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(15., track.passedDCAxy() * 1.); + registry.get(HIST("MBRDiff/notPVTracks"))->Fill(16., track.passedDCAz() * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("rejectedTracksDiff2"))->Fill(0., 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(1., track.isGlobalTrackSDD() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(2., track.passedTrackType() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(3., track.passedPtRange() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(4., track.passedEtaRange() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(5., track.passedTPCNCls() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(6., track.passedTPCCrossedRows() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(8., track.passedTPCChi2NDF() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(9., track.passedTPCRefit() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(10., track.passedITSNCls() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(11., track.passedITSChi2NDF() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(12., track.passedITSRefit() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(13., track.passedITSHits() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(14., track.passedGoldenChi2() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(15., track.passedDCAxy() * 1.); - registry.get(HIST("rejectedTracksDiff2"))->Fill(16., track.passedDCAz() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(0., 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(1., track.isGlobalTrackSDD() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(2., track.passedTrackType() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(3., track.passedPtRange() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(4., track.passedEtaRange() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(5., track.passedTPCNCls() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(6., track.passedTPCCrossedRows() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(8., track.passedTPCChi2NDF() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(9., track.passedTPCRefit() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(10., track.passedITSNCls() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(11., track.passedITSChi2NDF() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(12., track.passedITSRefit() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(13., track.passedITSHits() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(14., track.passedGoldenChi2() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(15., track.passedDCAxy() * 1.); + registry.get(HIST("Diff/notPVTracks"))->Fill(16., track.passedDCAz() * 1.); } else { - registry.get(HIST("rejectedTracks"))->Fill(0., 1.); - registry.get(HIST("rejectedTracks"))->Fill(1., track.isGlobalTrackSDD() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(2., track.passedTrackType() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(3., track.passedPtRange() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(4., track.passedEtaRange() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(5., track.passedTPCNCls() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(6., track.passedTPCCrossedRows() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(8., track.passedTPCChi2NDF() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(9., track.passedTPCRefit() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(10., track.passedITSNCls() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(11., track.passedITSChi2NDF() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(12., track.passedITSRefit() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(13., track.passedITSHits() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(14., track.passedGoldenChi2() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(15., track.passedDCAxy() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(16., track.passedDCAz() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(0., 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(1., track.isGlobalTrackSDD() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(2., track.passedTrackType() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(3., track.passedPtRange() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(4., track.passedEtaRange() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(5., track.passedTPCNCls() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(6., track.passedTPCCrossedRows() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(8., track.passedTPCChi2NDF() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(9., track.passedTPCRefit() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(10., track.passedITSNCls() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(11., track.passedITSChi2NDF() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(12., track.passedITSRefit() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(13., track.passedITSHits() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(14., track.passedGoldenChi2() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(15., track.passedDCAxy() * 1.); + registry.get(HIST("nonDiff/notPVTracks"))->Fill(16., track.passedDCAz() * 1.); } continue; } @@ -618,19 +623,19 @@ struct DiffMCQA { } isDGcandidate &= goodpts; if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(14., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(14., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(14., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(14., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(14., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(14., isDGcandidate * 1.); } isDGcandidate &= goodetas; if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(15., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(15., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(15., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(15., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(15., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(15., isDGcandidate * 1.); } auto netChargeValues = diffCuts.netCharges(); if (std::find(netChargeValues.begin(), netChargeValues.end(), netCharge) == netChargeValues.end()) { @@ -638,150 +643,122 @@ struct DiffMCQA { } isDGcandidate &= goodnchs; if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(16., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(16., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(16., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(16., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(16., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(16., isDGcandidate * 1.); } if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(17., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(17., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(17., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(17., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(17., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(17., isDGcandidate * 1.); } isDGcandidate &= (ivm.M() >= diffCuts.minIVM()); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(18., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(18., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(18., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(18., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(18., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(18., isDGcandidate * 1.); } isDGcandidate &= (ivm.M() <= diffCuts.maxIVM()); if (isPythiaDiff) { - registry.get(HIST("StatDiff1"))->Fill(19., isDGcandidate * 1.); + registry.get(HIST("MBRDiff/Stat"))->Fill(19., isDGcandidate * 1.); } else if (isGraniittiDiff) { - registry.get(HIST("StatDiff2"))->Fill(19., isDGcandidate * 1.); + registry.get(HIST("Diff/Stat"))->Fill(19., isDGcandidate * 1.); } else { - registry.get(HIST("Stat"))->Fill(19., isDGcandidate * 1.); + registry.get(HIST("nonDiff/Stat"))->Fill(19., isDGcandidate * 1.); } // update some DG histograms if (isDGcandidate) { if (isPythiaDiff) { - registry.get(HIST("vtxPosxyDGDiff1"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPoszDGDiff1"))->Fill(collision.posZ()); - registry.get(HIST("netChargeDGDiff1"))->Fill(netCharge); - registry.get(HIST("IVMptSysDGDiff1"))->Fill(ivm.M(), ivm.Perp()); + registry.get(HIST("MBRDiff/PVposxyDG"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("MBRDiff/PVposzDG"))->Fill(collision.posZ()); + registry.get(HIST("MBRDiff/netChargeDG"))->Fill(netCharge); + registry.get(HIST("MBRDiff/IVMptSysDG"))->Fill(ivm.M(), ivm.Perp()); } else if (isGraniittiDiff) { - registry.get(HIST("vtxPosxyDGDiff2"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPoszDGDiff2"))->Fill(collision.posZ()); - registry.get(HIST("netChargeDGDiff2"))->Fill(netCharge); - registry.get(HIST("IVMptSysDGDiff2"))->Fill(ivm.M(), ivm.Perp()); + registry.get(HIST("Diff/PVposxyDG"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("Diff/PVposzDG"))->Fill(collision.posZ()); + registry.get(HIST("Diff/netChargeDG"))->Fill(netCharge); + registry.get(HIST("Diff/IVMptSysDG"))->Fill(ivm.M(), ivm.Perp()); } else { - registry.get(HIST("vtxPosxyDG"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPoszDG"))->Fill(collision.posZ()); - registry.get(HIST("netChargeDG"))->Fill(netCharge); - registry.get(HIST("IVMptSysDG"))->Fill(ivm.M(), ivm.Perp()); + registry.get(HIST("nonDiff/PVposxyDG"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("nonDiff/PVposzDG"))->Fill(collision.posZ()); + registry.get(HIST("nonDiff/netChargeDG"))->Fill(netCharge); + registry.get(HIST("nonDiff/IVMptSysDG"))->Fill(ivm.M(), ivm.Perp()); } // fill dEdx of DG event tracks for (auto& track : tracks) { if (track.isPVContributor()) { if (isPythiaDiff) { - registry.get(HIST("etaptDGDiff1"))->Fill(track.eta(), track.pt()); - registry.get(HIST("dEdxTPCDGDiff1"))->Fill(track.p(), track.tpcSignal()); - registry.get(HIST("IVMptTrkDGDiff1"))->Fill(ivm.M(), track.pt()); + registry.get(HIST("MBRDiff/etaptDG"))->Fill(track.eta(), track.pt()); + registry.get(HIST("MBRDiff/dEdxTPCDG"))->Fill(track.p(), track.tpcSignal()); + registry.get(HIST("MBRDiff/IVMptTrkDG"))->Fill(ivm.M(), track.pt()); } else if (isGraniittiDiff) { - registry.get(HIST("etaptDGDiff2"))->Fill(track.eta(), track.pt()); - registry.get(HIST("dEdxTPCDGDiff2"))->Fill(track.p(), track.tpcSignal()); - registry.get(HIST("IVMptTrkDGDiff2"))->Fill(ivm.M(), track.pt()); + registry.get(HIST("Diff/etaptDG"))->Fill(track.eta(), track.pt()); + registry.get(HIST("Diff/dEdxTPCDG"))->Fill(track.p(), track.tpcSignal()); + registry.get(HIST("Diff/IVMptTrkDG"))->Fill(ivm.M(), track.pt()); } else { - registry.get(HIST("etaptDG"))->Fill(track.eta(), track.pt()); - registry.get(HIST("dEdxTPCDG"))->Fill(track.p(), track.tpcSignal()); - registry.get(HIST("IVMptTrkDG"))->Fill(ivm.M(), track.pt()); + registry.get(HIST("nonDiff/etaptDG"))->Fill(track.eta(), track.pt()); + registry.get(HIST("nonDiff/dEdxTPCDG"))->Fill(track.p(), track.tpcSignal()); + registry.get(HIST("nonDiff/IVMptTrkDG"))->Fill(ivm.M(), track.pt()); } if (track.hasTOF()) { if (isPythiaDiff) { - registry.get(HIST("dEdxTOFDGDiff1"))->Fill(track.p(), track.tofSignal()); + registry.get(HIST("MBRDiff/dEdxTOFDG"))->Fill(track.p(), track.tofSignal()); } else if (isGraniittiDiff) { - registry.get(HIST("dEdxTOFDGDiff2"))->Fill(track.p(), track.tofSignal()); + registry.get(HIST("Diff/dEdxTOFDG"))->Fill(track.p(), track.tofSignal()); } else { - registry.get(HIST("dEdxTOFDG"))->Fill(track.p(), track.tofSignal()); + registry.get(HIST("nonDiff/dEdxTOFDG"))->Fill(track.p(), track.tofSignal()); } } } } } - }; -}; - -struct FV0Signals { - - HistogramRegistry registry{ - "registry", - { - {"FV0A", "#FV0A", {HistType::kTH2F, {{48, -0.5, 47.5}, {1000, 0., 1000.}}}}, - }}; + } + PROCESS_SWITCH(DiffMCQA, processMain, "Process Main", true); - void process(aod::FV0A const& fv0) + void processFV0(aod::FV0A const& fv0) { // side A for (size_t ind = 0; ind < fv0.channel().size(); ind++) { - registry.get(HIST("FV0A"))->Fill((fv0.channel())[ind], (fv0.amplitude())[ind]); + registry.get(HIST("FV0/FV0Aamp"))->Fill((fv0.channel())[ind], (fv0.amplitude())[ind]); } - }; -}; - -struct FT0Signals { - - HistogramRegistry registry{ - "registry", - { - {"FT0A", "#FT0A", {HistType::kTH2F, {{96, -0.5, 95.5}, {100, 0., 200.}}}}, - {"FT0C", "#FT0C", {HistType::kTH2F, {{112, -0.5, 111.5}, {100, 0., 200.}}}}, - }}; + } + PROCESS_SWITCH(DiffMCQA, processFV0, "Process FV0", true); - void process(aod::FT0 const& ft0) + void processFT0(aod::FT0 const& ft0) { // side A for (size_t ind = 0; ind < ft0.channelA().size(); ind++) { - registry.get(HIST("FT0A"))->Fill((ft0.channelA())[ind], (ft0.amplitudeA())[ind]); + registry.get(HIST("FT0/FT0Aamp"))->Fill((ft0.channelA())[ind], (ft0.amplitudeA())[ind]); } // side C for (size_t ind = 0; ind < ft0.channelC().size(); ind++) { - registry.get(HIST("FT0C"))->Fill((ft0.channelC())[ind], (ft0.amplitudeC())[ind]); + registry.get(HIST("FT0/FT0Camp"))->Fill((ft0.channelC())[ind], (ft0.amplitudeC())[ind]); } - }; -}; - -struct FDDSignals { - - HistogramRegistry registry{ - "registry", - { - {"FDDA", "#FDDA", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}}, - {"FDDC", "#FDDC", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}}, - }}; + } + PROCESS_SWITCH(DiffMCQA, processFT0, "Process FT0", true); - void process(aod::FDD const& fdd) + void processFDD(aod::FDD const& fdd) { // side A for (auto ind = 0; ind < 8; ind++) { - registry.get(HIST("FDDA"))->Fill(ind, (fdd.chargeA())[ind]); + registry.get(HIST("FDD/FDDAamp"))->Fill(ind, (fdd.chargeA())[ind]); } // side C for (auto ind = 0; ind < 8; ind++) { - registry.get(HIST("FDDC"))->Fill(ind, (fdd.chargeC())[ind]); + registry.get(HIST("FDD/FDDCamp"))->Fill(ind, (fdd.chargeC())[ind]); } - }; -}; - -struct ZDCSignals { + } + PROCESS_SWITCH(DiffMCQA, processFDD, "Process FDD", true); // energies: // 0: energyZEM1 @@ -806,13 +783,7 @@ struct ZDCSignals { // 19: energySectorZPC[1] // 20: energySectorZPC[2] // 21: energySectorZPC[3] - HistogramRegistry registry{ - "registry", - { - {"ZdcEnergies", "#ZdcEnergies", {HistType::kTH2F, {{22, -0.5, 21.5}, {100, 0., 1000.}}}}, - }}; - - void process(aod::Zdc const& zdc) + void processZDC(aod::Zdc const& zdc) { // Zdc energies registry.get(HIST("ZdcEnergies"))->Fill(0., zdc.energyZEM1()); @@ -837,36 +808,22 @@ struct ZDCSignals { registry.get(HIST("ZdcEnergies"))->Fill(19., (zdc.energySectorZPC())[1]); registry.get(HIST("ZdcEnergies"))->Fill(20., (zdc.energySectorZPC())[2]); registry.get(HIST("ZdcEnergies"))->Fill(21., (zdc.energySectorZPC())[3]); - }; -}; - -struct CaloSignals { - - HistogramRegistry registry{ - "registry", - { - {"CaloCell", "#CaloCell", {HistType::kTH1I, {{18000, -0.5, 17999.5}}}}, - {"CaloAmplitude", "#CaloAmplitude", {HistType::kTH1F, {{100, 0, 10.}}}}, - }}; + } + PROCESS_SWITCH(DiffMCQA, processZDC, "Process ZDC", true); - void process(aod::Calo const& calo) + void processCalo(aod::Calo const& calo) { // cell number registry.get(HIST("CaloCell"))->Fill(calo.cellNumber()); // amplitude registry.get(HIST("CaloAmplitude"))->Fill(calo.amplitude()); - }; + } + PROCESS_SWITCH(DiffMCQA, processCalo, "Process Calorimeter", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"diffmcqa"}), - adaptAnalysisTask(cfgc, TaskName{"fv0signals"}), - adaptAnalysisTask(cfgc, TaskName{"ft0signals"}), - adaptAnalysisTask(cfgc, TaskName{"fddsignals"}), - adaptAnalysisTask(cfgc, TaskName{"zdcsignals"}), - adaptAnalysisTask(cfgc, TaskName{"calosignals"}), - }; + adaptAnalysisTask(cfgc, TaskName{"diffmcqa"})}; } diff --git a/PWGUD/Tasks/diffQA.cxx b/PWGUD/Tasks/diffQA.cxx index a2c4da5b0a7..7534f347ea3 100644 --- a/PWGUD/Tasks/diffQA.cxx +++ b/PWGUD/Tasks/diffQA.cxx @@ -16,6 +16,7 @@ #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "TLorentzVector.h" +#include "CommonConstants/LHCConstants.h" #include "ReconstructionDataFormats/BCRange.h" #include "CommonConstants/PhysicsConstants.h" #include "PWGUD/Core/UDHelpers.h" @@ -30,15 +31,14 @@ struct DiffQA { Preslice perBCcalo = aod::calo::bcId; // constants - static const int nBCpOrbit = 3564; static const int ns = 20; // number of BCs to save (used in processF[V0, T0, DD]) static const int ncmin = 20; // minimum length of series of empty BCs (used in processF[V0, T0, DD]) - static constexpr std::string_view hFV0A[ns + 1] = {"fv0A00", "fv0A01", "fv0A02", "fv0A03", "fv0A04", "fv0A05", "fv0A06", "fv0A07", "fv0A08", "fv0A09", "fv0A10", "fv0A11", "fv0A12", "fv0A13", "fv0A14", "fv0A15", "fv0A16", "fv0A17", "fv0A18", "fv0A19", "fv0A20"}; - static constexpr std::string_view hFT0A[ns + 1] = {"ft0A00", "ft0A01", "ft0A02", "ft0A03", "ft0A04", "ft0A05", "ft0A06", "ft0A07", "ft0A08", "ft0A09", "ft0A10", "ft0A11", "ft0A12", "ft0A13", "ft0A14", "ft0A15", "ft0A16", "ft0A17", "ft0A18", "ft0A19", "ft0A20"}; - static constexpr std::string_view hFT0C[ns + 1] = {"ft0C00", "ft0C01", "ft0C02", "ft0C03", "ft0C04", "ft0C05", "ft0C06", "ft0C07", "ft0C08", "ft0C09", "ft0C10", "ft0C11", "ft0C12", "ft0C13", "ft0C14", "ft0C15", "ft0C16", "ft0C17", "ft0C18", "ft0C19", "ft0C20"}; - static constexpr std::string_view hFDDA[ns + 1] = {"fddA00", "fddA01", "fddA02", "fddA03", "fddA04", "fddA05", "fddA06", "fddA07", "fddA08", "fddA09", "fddA10", "fddA11", "fddA12", "fddA13", "fddA14", "fddA15", "fddA16", "fddA17", "fddA18", "fddA19", "fddA20"}; - static constexpr std::string_view hFDDC[ns + 1] = {"fddC00", "fddC01", "fddC02", "fddC03", "fddC04", "fddC05", "fddC06", "fddC07", "fddC08", "fddC09", "fddC10", "fddC11", "fddC12", "fddC13", "fddC14", "fddC15", "fddC16", "fddC17", "fddC18", "fddC19", "fddC20"}; + static constexpr std::string_view hFV0A[ns + 1] = {"FV0/A00", "FV0/A01", "FV0/A02", "FV0/A03", "FV0/A04", "FV0/A05", "FV0/A06", "FV0/A07", "FV0/A08", "FV0/A09", "FV0/A10", "FV0/A11", "FV0/A12", "FV0/A13", "FV0/A14", "FV0/A15", "FV0/A16", "FV0/A17", "FV0/A18", "FV0/A19", "FV0/A20"}; + static constexpr std::string_view hFT0A[ns + 1] = {"FT0/A00", "FT0/A01", "FT0/A02", "FT0/A03", "FT0/A04", "FT0/A05", "FT0/A06", "FT0/A07", "FT0/A08", "FT0/A09", "FT0/A10", "FT0/A11", "FT0/A12", "FT0/A13", "FT0/A14", "FT0/A15", "FT0/A16", "FT0/A17", "FT0/A18", "FT0/A19", "FT0/A20"}; + static constexpr std::string_view hFT0C[ns + 1] = {"FT0/C00", "FT0/C01", "FT0/C02", "FT0/C03", "FT0/C04", "FT0/C05", "FT0/C06", "FT0/C07", "FT0/C08", "FT0/C09", "FT0/C10", "FT0/C11", "FT0/C12", "FT0/C13", "FT0/C14", "FT0/C15", "FT0/C16", "FT0/C17", "FT0/C18", "FT0/C19", "FT0/C20"}; + static constexpr std::string_view hFDDA[ns + 1] = {"FDD/A00", "FDD/A01", "FDD/A02", "FDD/A03", "FDD/A04", "FDD/A05", "FDD/A06", "FDD/A07", "FDD/A08", "FDD/A09", "FDD/A10", "FDD/A11", "FDD/A12", "FDD/A13", "FDD/A14", "FDD/A15", "FDD/A16", "FDD/A17", "FDD/A18", "FDD/A19", "FDD/A20"}; + static constexpr std::string_view hFDDC[ns + 1] = {"FDD/C00", "FDD/C01", "FDD/C02", "FDD/C03", "FDD/C04", "FDD/C05", "FDD/C06", "FDD/C07", "FDD/C08", "FDD/C09", "FDD/C10", "FDD/C11", "FDD/C12", "FDD/C13", "FDD/C14", "FDD/C15", "FDD/C16", "FDD/C17", "FDD/C18", "FDD/C19", "FDD/C20"}; // global variables float maxdEdxTPC; @@ -78,57 +78,81 @@ struct DiffQA { // add histograms for the different process functions if (context.mOptions.get("processMain")) { - registry.add("Stat", "#Stat", {HistType::kTH1F, {{20, -0.5, 19.5}}}); - registry.add("Tracks", "#Tracks", {HistType::kTH1F, {{50, 0.5, 50.5}}}); - registry.add("vtxTracks", "#vtxTracks", {HistType::kTH1F, {{50, 0.5, 50.5}}}); - registry.add("globalTracks", "#globalTracks", {HistType::kTH1F, {{50, 0.5, 50.5}}}); - registry.add("rejectedTracks", "#rejectedTracks", {HistType::kTH1F, {{17, -0.5, 16.5}}}); - registry.add("tResvsrTOFTracks", "#tResvsrTOFTracks", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}); - registry.add("vtxPosxy", "#vtxPosxy", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); - registry.add("vtxPosz", "#vtxPosz", {HistType::kTH1F, {{1000, -100., 100.}}}); - registry.add("etapt", "#etapt", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); - registry.add("dEdxTPC", "#dEdxTPC", {HistType::kTH2F, {{120, -6., 6.}, {1000, 0., 1000.}}}); - registry.add("dEdxTOF", "#dEdxTOF", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); - registry.add("vtxPosxyDG", "#vtxPosxyDG", {HistType::kTH2F, {{200, -2., 2.}, {200, -2., 2.}}}); - registry.add("vtxPoszDG", "#vtxPoszDG", {HistType::kTH1F, {{1000, -100., 100.}}}); - registry.add("etaptDG", "#etaptDG", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); - registry.add("dEdxTPCDG", "#dEdxTPCDG", {HistType::kTH2F, {{120, -6., 6.0}, {1000, 0., 1000.}}}); - registry.add("dEdxTOFDG", "#dEdxTOFDG", {HistType::kTH2F, {{100, 0., 5.0}, {1000, 0., 500000.}}}); - registry.add("netChargeDG", "#netChargeDG", {HistType::kTH1F, {{21, -10.5, 10.5}}}); - registry.add("IVMptSysDG", "#IVMptSysDG", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); - registry.add("IVMptTrkDG", "#IVMptTrkDG", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + // collisions + registry.add("collisions/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{20, -0.5, 19.5}}}); + registry.add("collisions/PVposxy", "Vertex position in x and y direction; V_x; V_y; Collisions", {HistType::kTH2F, {{100, -0.5, 0.5}, {100, -0.5, 0.5}}}); + registry.add("collisions/PVposz", "Vertex position in z direction; V_z; Collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("collisions/Tracks", "Number of tracks; Number of tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("collisions/PVTracks", "Number of PV tracks; Number of PV tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("collisions/globalTracks", "Number of global tracks; Number of global tracks; Collisions", {HistType::kTH1F, {{50, 0.5, 50.5}}}); + registry.add("collisions/notPVTracks", "Not PV tracks; Track status bit; Not PV tracks", {HistType::kTH1F, {{17, -0.5, 16.5}}}); + registry.add("collisions/tResvsrTOFTracks", "Number of PV tracks with TOF hit versus collision time resolution; Collision time resolution [ns]; Fraction of PV tracks with TOF hit; Collisions", {HistType::kTH2F, {{1000, 0., 1.E3}, {101, -0.01, 1.01}}}); + + // tracks + registry.add("tracks/Stat", "Track bits as function of pT; Track pT; Track bit; Tracks", {HistType::kTH2F, {{100, 0., 5.}, {8, 0.5, 8.5}}}); + + registry.add("tracks/etapt", "eta versus pT of all tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("tracks/etapt2", "eta versus pT of all quality tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("tracks/etapt3", "eta versus pT of all global tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("tracks/etapt4", "eta versus pT of all PV tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("tracks/etapt5", "eta versus pT of all tracks with ITS; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("tracks/etapt6", "eta versus pT of all tracks with TPC; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("tracks/etapt7", "eta versus pT of all tracks with TRD; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("tracks/etapt8", "eta versus pT of all tracks with TOF; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + + registry.add("tracks/dEdxTPC", "TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{120, -6., 6.}, {1000, 0., 1000.}}}); + registry.add("tracks/dEdxTOF", "TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{120, -6., 6.}, {1000, 0., 30000.}}}); + + // DG + registry.add("DG/PVposxy", "DG: Vertex position in x and y direction; V_x [mm]; V_y [mm]; DG collisions", {HistType::kTH2F, {{100, -0.5, 0.5}, {100, -0.5, 0.5}}}); + registry.add("DG/PVposz", "DG: Vertex position in z direction; V_z; DG collisions", {HistType::kTH1F, {{1000, -100., 100.}}}); + registry.add("DG/netCharge", "DG: net charge; Net charge; DG collisions", {HistType::kTH1F, {{21, -10.5, 10.5}}}); + registry.add("DG/TrackStat", "Track bits as function of pT; Track pT; Track bit; Tracks", {HistType::kTH2F, {{100, 0., 5.}, {8, 0.5, 8.5}}}); + + registry.add("DG/etapt", "DG: eta versus pT of all tracks; eta of track; p_T of track [GeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("DG/etapt2", "DG: eta versus pT of all quality tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("DG/etapt3", "DG: eta versus pT of all global tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("DG/etapt4", "DG: eta versus pT of all PV tracks; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("DG/etapt5", "DG: eta versus pT of all tracks with ITS; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("DG/etapt6", "DG: eta versus pT of all tracks with TPC; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("DG/etapt7", "DG: eta versus pT of all tracks with TRD; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + registry.add("DG/etapt8", "DG: eta versus pT of all tracks with TOF; eta of track; p_T of track [MeV/c^2]; Tracks", {HistType::kTH2F, {{80, -2., 2.}, {100, 0., 5.}}}); + + registry.add("DG/dEdxTPC", "DG: TPC signal versus signed track momentum; Signed track momentum [GeV/c]; TPC signal; Tracks", {HistType::kTH2F, {{120, -6., 6.}, {1000, 0., 1000.}}}); + registry.add("DG/dEdxTOF", "DG: TOF signal versus signed track momentum; Signed track momentum [GeV/c]; TOF signal; Tracks", {HistType::kTH2F, {{120, -6., 6.}, {1000, 0., 30000.}}}); + registry.add("DG/IVMptSys", "DG: Invariant mass versus p_{T, system}; Invarian mass [GeV/c^2]; p_{T, system} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); + registry.add("DG/IVMptTrk", "DG: Invariant mass versus p_{T, tracks}; Invarian mass [GeV/c^2]; p_{T, tracks} [GeV/c]; DG collisions", {HistType::kTH2F, {{100, 0., 5.}, {350, 0., 3.5}}}); } - if (context.mOptions.get("processFewProng")) { - registry.add("fpStat", "#fpStat", {HistType::kTH1F, {{2, 0.5, 2.5}}}); - registry.add("allPVC", "#allPVC", {HistType::kTH1F, {{100, 0.5, 100.5}}}); - registry.add("fpPVC", "#fpPVC", {HistType::kTH1F, {{100, 0.5, 100.5}}}); + if (context.mOptions.get("processCleanFT0")) { + registry.add("cleanFT0/Stat", "Statistics of collisions with clean FT0; FT0 status; Collisions", {HistType::kTH1F, {{2, 0.5, 2.5}}}); + registry.add("cleanFT0/PVTracks", "Distribution of number of PV contributors for all collisions and for collisions with clean FT0; Number of PV contributors;", {HistType::kTH2F, {{100, 0.5, 100.5}, {2, 0.5, 2.5}}}); } if (context.mOptions.get("processCleanFIT1")) { - registry.add("cleanFIT1", "#cleanFIT1", {HistType::kTH2F, {{20, -0.5, 19.5}, {2, -0.5, 1.5}}}); - registry.add("cF1FV0Aamp", "#cF1FV0Aamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF1FT0Aamp", "#cF1FT0Aamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF1FT0Camp", "#cF1FT0Camp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF1FDDAamp", "#cF1FDDAamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF1FDDCamp", "#cF1FDDCamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT1/Stat", "Statistics of collisions with empty FT0; Multiple of collision time resolution; FT0 status; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {2, -0.5, 1.5}}}); + registry.add("cleanFIT1/FV0Aamp", "Amplitude of FV0A in collisions with empty FT0; Multiple of collision time resolution; FV0A amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT1/FT0Aamp", "Amplitude of FT0A in collisions with empty FT0; Multiple of collision time resolution; FT0A amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT1/FT0Camp", "Amplitude of FT0C in collisions with empty FT0; Multiple of collision time resolution; FT0C amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT1/FDDAamp", "Amplitude of FDDA in collisions with empty FT0; Multiple of collision time resolution; FDDA amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT1/FDDCamp", "Amplitude of FDDC in collisions with empty FT0; Multiple of collision time resolution; FDDC amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); } if (context.mOptions.get("processCleanFIT2")) { - registry.add("cleanFIT2", "#cleanFIT2", {HistType::kTH2F, {{20, -0.5, 19.5}, {2, -0.5, 1.5}}}); - registry.add("cF2FV0Aamp", "#cF2FV0Aamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF2FT0Aamp", "#cF2FT0Aamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF2FT0Camp", "#cF2FT0Camp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF2FDDAamp", "#cF2FDDAamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); - registry.add("cF2FDDCamp", "#cF2FDDCamp", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT2/Stat", "Statistics of collisions with empty FIT; Number of neighbouring BCs; FIT status; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {2, -0.5, 1.5}}}); + registry.add("cleanFIT2/FV0Aamp", "Amplitude of FV0A in collisions with empty FIT; Number of neighbouring BCs; FV0A amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT2/FT0Aamp", "Amplitude of FT0A in collisions with empty FIT; Number of neighbouring BCs; FT0A amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT2/FT0Camp", "Amplitude of FT0C in collisions with empty FIT; Number of neighbouring BCs; FT0C amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT2/FDDAamp", "Amplitude of FDDA in collisions with empty FIT; Number of neighbouring BCs; FDDA amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); + registry.add("cleanFIT2/FDDCamp", "Amplitude of FDDC in collisions with empty FIT; Number of neighbouring BCs; FDDC amplitude; Collisions", {HistType::kTH2F, {{20, -0.5, 19.5}, {1000, -0.5, 999.5}}}); } if (context.mOptions.get("processFV0")) { - registry.add("FV0A", "#FV0A", {HistType::kTH2F, {{48, -0.5, 47.5}, {2000, 0., 2000.}}}); - registry.add("FV0BCNUM", "#FV0BCNUM", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); + registry.add("FV0/FV0Aamp", "FV0A amplitudes; FV0A channel; FV0A amplitude; Entries", {HistType::kTH2F, {{48, -0.5, 47.5}, {2000, 0., 2000.}}}); + registry.add("FV0/emptyBCs", "Distribution of number of consecutive BCs with empty FV0A; Number of consecutive BCs with empty FV0A; Entries", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); } if (context.mOptions.get("processFT0")) { - registry.add("FT0A", "#FT0A", {HistType::kTH2F, {{96, -0.5, 95.5}, {400, 0., 400.}}}); - registry.add("FT0C", "#FT0C", {HistType::kTH2F, {{112, -0.5, 111.5}, {400, 0., 400.}}}); - registry.add("FT0ABCNUM", "#FT0ABCNUM", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); - registry.add("FT0CBCNUM", "#FT0CBCNUM", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); - registry.add("dFT0BCNUM", "#dFT0BCNUM", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); + registry.add("FT0/FT0Aamp", "FT0A amplitudes; FT0A channel; FT0A amplitude; Entries", {HistType::kTH2F, {{96, -0.5, 95.5}, {400, 0., 400.}}}); + registry.add("FT0/FT0Camp", "FT0C amplitudes; FT0C channel; FT0C amplitude; Entries", {HistType::kTH2F, {{112, -0.5, 111.5}, {400, 0., 400.}}}); + registry.add("FT0/AP2BC", "P2 BCs with FT0A signal; P2 BC; Entries", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); + registry.add("FT0/CP2BC", "P2 BCs with FT0C signal; P2 BC; Entries", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); + registry.add("FT0/emptyBCs", "Distribution of number of consecutive BCs with empty FT0; Number of consecutive BCs with empty FT0; Entries", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); // add amplitude histograms for (auto n{0}; n <= ns; n++) { @@ -137,12 +161,12 @@ struct DiffQA { } } if (context.mOptions.get("processFDD")) { - registry.add("FDDA", "#FDDA", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}); - registry.add("FDDC", "#FDDC", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}); - registry.add("FDDBCNUM", "#FDDBCNUM", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); + registry.add("FDD/FDDAamp", "FDDA amplitudes; FDDA channel; FDDA amplitude; Entries", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}); + registry.add("FDD/FDDCamp", "FDDA amplitudes; FDDA channel; FDDA amplitude; Entries", {HistType::kTH2F, {{8, -0.5, 7.5}, {100, 0., 100.}}}); + registry.add("FDD/emptyBCs", "Distribution of number of consecutive BCs with empty FDD; Number of consecutive BCs with empty FDD; Entries", {HistType::kTH1F, {{3564, -0.5, 3563.5}}}); } if (context.mOptions.get("processZDC")) { - registry.add("ZdcEnergies", "#ZdcEnergies", {HistType::kTH2F, {{22, -0.5, 21.5}, {100, 0., 1000.}}}); + registry.add("ZDC/Energies", "Registered energies in various ZDC channels; ZDC channel; Energy; Entries", {HistType::kTH2F, {{22, -0.5, 21.5}, {100, 0., 1000.}}}); } // if (context.mOptions.get("processCalo")) { // registry.add("CaloCell", "#CaloCell", {HistType::kTH1I, {{18000, -0.5, 17999.5}}}); @@ -200,32 +224,50 @@ struct DiffQA { aod::Zdcs& zdcs, aod::Calos& calos, aod::V0s const& v0s, aod::Cascades const& cascades) { - LOGF(debug, " Collision %d", collision.globalIndex()); LOGF(debug, " Start %i", abcrs.size()); bool isDGcandidate = true; - registry.get(HIST("Stat"))->Fill(0., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(0., isDGcandidate * 1.); // update collision histograms // vertex position - registry.get(HIST("vtxPosxy"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPosz"))->Fill(collision.posZ()); + registry.get(HIST("collisions/PVposxy"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("collisions/PVposz"))->Fill(collision.posZ()); // tracks - registry.get(HIST("Tracks"))->Fill(tracks.size()); + registry.get(HIST("collisions/Tracks"))->Fill(tracks.size()); // vertex tracks - registry.get(HIST("vtxTracks"))->Fill(collision.numContrib()); + registry.get(HIST("collisions/PVTracks"))->Fill(collision.numContrib()); // global tracks Partition goodTracks = requireGlobalTrackInFilter(); goodTracks.bindTable(tracks); - registry.get(HIST("globalTracks"))->Fill(goodTracks.size()); + registry.get(HIST("collisions/globalTracks"))->Fill(goodTracks.size()); - // number of vertex tracks with TOF hit + // loop over all tracks float rgtrwTOF = 0.; for (auto const& track : tracks) { - // update eta vs pt histogram - registry.get(HIST("etapt"))->Fill(track.eta(), track.pt()); + // update track stats + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 1., 1.); + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 2., track.isQualityTrack() * 1.); + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 3., track.isGlobalTrack() * 1.); + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 4., track.isPVContributor() * 1.); + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 5., track.hasITS() * 1.); + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 6., track.hasTPC() * 1.); + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 7., track.hasTRD() * 1.); + registry.get(HIST("tracks/Stat"))->Fill(track.pt(), 8., track.hasTOF() * 1.); + + // update eta vs pt histograms + registry.get(HIST("tracks/etapt"))->Fill(track.eta(), track.pt(), 1.); + registry.get(HIST("tracks/etapt2"))->Fill(track.eta(), track.pt(), track.isQualityTrack() * 1.); + registry.get(HIST("tracks/etapt3"))->Fill(track.eta(), track.pt(), track.isGlobalTrack() * 1.); + registry.get(HIST("tracks/etapt4"))->Fill(track.eta(), track.pt(), track.isPVContributor() * 1.); + registry.get(HIST("tracks/etapt5"))->Fill(track.eta(), track.pt(), track.hasITS() * 1.); + registry.get(HIST("tracks/etapt6"))->Fill(track.eta(), track.pt(), track.hasTPC() * 1.); + registry.get(HIST("tracks/etapt7"))->Fill(track.eta(), track.pt(), track.hasTRD() * 1.); + registry.get(HIST("tracks/etapt8"))->Fill(track.eta(), track.pt(), track.hasTOF() * 1.); + // update dEdx histograms - registry.get(HIST("dEdxTPC"))->Fill(track.p() * track.sign(), track.tpcSignal()); + registry.get(HIST("tracks/dEdxTPC"))->Fill(track.tpcInnerParam() * track.sign(), track.tpcSignal()); if (track.tpcSignal() > maxdEdxTPC) { maxdEdxTPC = track.tpcSignal(); LOGF(debug, " New maxdEdx TPC %f", maxdEdxTPC); @@ -233,10 +275,10 @@ struct DiffQA { // TOF hit? if (track.hasTOF()) { - registry.get(HIST("dEdxTOF"))->Fill(track.pt(), track.tofSignal()); + registry.get(HIST("tracks/dEdxTOF"))->Fill(track.p() * track.sign(), track.tofSignal()); if (track.tofSignal() > maxdEdxTOF) { maxdEdxTOF = track.tofSignal(); - LOGF(debug, " New maxdEdx tOF %f", maxdEdxTOF); + LOGF(debug, " New maxdEdx TOF %f", maxdEdxTOF); } // vertex track with TOF hit? @@ -245,11 +287,12 @@ struct DiffQA { } } } + // fraction of PV tracks with TOF hit if (collision.numContrib() > 0) { rgtrwTOF /= collision.numContrib(); } - LOGF(debug, " Vertex tracks with TOF: %f [1]", rgtrwTOF); - registry.get(HIST("tResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rgtrwTOF); + LOGF(debug, " PV tracks with TOF: %f [1]", rgtrwTOF); + registry.get(HIST("collisions/tResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rgtrwTOF); // is it a DG candidate? // DG = no FIT signal in compatible BCs @@ -278,7 +321,7 @@ struct DiffQA { isDGcandidate = false; } } - registry.get(HIST("Stat"))->Fill(1., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(1., isDGcandidate * 1.); // no Zdc signal in bcSlice std::vector lims(10, 0.); @@ -288,7 +331,7 @@ struct DiffQA { break; } } - registry.get(HIST("Stat"))->Fill(2., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(2., isDGcandidate * 1.); // no Calo signal in bcSlice for (auto const& bc : bcSlice) { @@ -297,19 +340,19 @@ struct DiffQA { break; } } - registry.get(HIST("Stat"))->Fill(3., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(3., isDGcandidate * 1.); // no V0s isDGcandidate &= (v0s.size() == 0); - registry.get(HIST("Stat"))->Fill(4., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(4., isDGcandidate * 1.); // no Cascades isDGcandidate &= (cascades.size() == 0); - registry.get(HIST("Stat"))->Fill(5., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(5., isDGcandidate * 1.); // number of forward tracks = 0 isDGcandidate &= (fwdtracks.size() == 0); - registry.get(HIST("Stat"))->Fill(6., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(6., isDGcandidate * 1.); // no global tracks which are no vtx tracks bool globalAndVtx = isDGcandidate; @@ -322,8 +365,8 @@ struct DiffQA { vtxAndGlobal = false; } } - registry.get(HIST("Stat"))->Fill(7., globalAndVtx * 1.); - registry.get(HIST("Stat"))->Fill(8., vtxAndGlobal * 1.); + registry.get(HIST("collisions/Stat"))->Fill(7., globalAndVtx * 1.); + registry.get(HIST("collisions/Stat"))->Fill(8., vtxAndGlobal * 1.); isDGcandidate &= globalAndVtx; if (diffCuts.globalTracksOnly()) { isDGcandidate &= vtxAndGlobal; @@ -337,7 +380,7 @@ struct DiffQA { break; } } - registry.get(HIST("Stat"))->Fill(9., noAmbTracks * 1.); + registry.get(HIST("collisions/Stat"))->Fill(9., noAmbTracks * 1.); // check a given bc for possible ambiguous FwdTracks auto noAmbFwdTracks = isDGcandidate; @@ -347,17 +390,17 @@ struct DiffQA { break; } } - registry.get(HIST("Stat"))->Fill(10., noAmbFwdTracks * 1.); + registry.get(HIST("collisions/Stat"))->Fill(10., noAmbFwdTracks * 1.); // fraction of PV tracks with TOF hit isDGcandidate &= (rgtrwTOF >= diffCuts.minRgtrwTOF()); - registry.get(HIST("Stat"))->Fill(11., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(11., isDGcandidate * 1.); // number of vertex tracks <= n isDGcandidate &= (collision.numContrib() >= diffCuts.minNTracks()); - registry.get(HIST("Stat"))->Fill(12., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(12., isDGcandidate * 1.); isDGcandidate &= (collision.numContrib() <= diffCuts.maxNTracks()); - registry.get(HIST("Stat"))->Fill(13., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(13., isDGcandidate * 1.); // net charge and invariant mass bool goodetas = true; @@ -376,25 +419,25 @@ struct DiffQA { // check also pt and eta of tracks for (auto const& track : tracks) { - // update histogram rejectedTracks + // update histogram notPVTracks if (!track.isPVContributor()) { - registry.get(HIST("rejectedTracks"))->Fill(0., 1.); - registry.get(HIST("rejectedTracks"))->Fill(1., track.isGlobalTrackSDD() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(2., track.passedTrackType() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(3., track.passedPtRange() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(4., track.passedEtaRange() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(5., track.passedTPCNCls() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(6., track.passedTPCCrossedRows() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(8., track.passedTPCChi2NDF() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(9., track.passedTPCRefit() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(10., track.passedITSNCls() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(11., track.passedITSChi2NDF() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(12., track.passedITSRefit() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(13., track.passedITSHits() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(14., track.passedGoldenChi2() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(15., track.passedDCAxy() * 1.); - registry.get(HIST("rejectedTracks"))->Fill(16., track.passedDCAz() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(0., 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(1., track.isGlobalTrackSDD() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(2., track.passedTrackType() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(3., track.passedPtRange() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(4., track.passedEtaRange() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(5., track.passedTPCNCls() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(6., track.passedTPCCrossedRows() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(7., track.passedTPCCrossedRowsOverNCls() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(8., track.passedTPCChi2NDF() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(9., track.passedTPCRefit() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(10., track.passedITSNCls() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(11., track.passedITSChi2NDF() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(12., track.passedITSRefit() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(13., track.passedITSHits() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(14., track.passedGoldenChi2() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(15., track.passedDCAxy() * 1.); + registry.get(HIST("collisions/notPVTracks"))->Fill(16., track.passedDCAz() * 1.); continue; } @@ -411,76 +454,72 @@ struct DiffQA { } } isDGcandidate &= goodpts; - registry.get(HIST("Stat"))->Fill(14., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(14., isDGcandidate * 1.); isDGcandidate &= goodetas; - registry.get(HIST("Stat"))->Fill(15., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(15., isDGcandidate * 1.); auto netChargeValues = diffCuts.netCharges(); if (std::find(netChargeValues.begin(), netChargeValues.end(), netCharge) == netChargeValues.end()) { goodnchs = false; } isDGcandidate &= goodnchs; - registry.get(HIST("Stat"))->Fill(16., isDGcandidate * 1.); - registry.get(HIST("Stat"))->Fill(17., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(16., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(17., isDGcandidate * 1.); isDGcandidate &= (ivm.M() >= diffCuts.minIVM()); - registry.get(HIST("Stat"))->Fill(18., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(18., isDGcandidate * 1.); isDGcandidate &= (ivm.M() <= diffCuts.maxIVM()); - registry.get(HIST("Stat"))->Fill(19., isDGcandidate * 1.); + registry.get(HIST("collisions/Stat"))->Fill(19., isDGcandidate * 1.); // update some DG histograms if (isDGcandidate) { // vertex position of DG events - registry.get(HIST("vtxPosxyDG"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vtxPoszDG"))->Fill(collision.posZ()); - registry.get(HIST("netChargeDG"))->Fill(netCharge); - registry.get(HIST("IVMptSysDG"))->Fill(ivm.M(), ivm.Perp()); + registry.get(HIST("DG/PVposxy"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("DG/PVposz"))->Fill(collision.posZ()); + registry.get(HIST("DG/netCharge"))->Fill(netCharge); + registry.get(HIST("DG/IVMptSys"))->Fill(ivm.M(), ivm.Perp()); - // fill dEdx of DG event tracks + // fill track status, eta vs pt, and dEdx of DG event tracks for (auto const& track : tracks) { if (track.isPVContributor()) { + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 1., 1.); + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 2., track.isQualityTrack() * 1.); + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 3., track.isGlobalTrack() * 1.); + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 4., 1.); + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 5., track.hasITS() * 1.); + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 6., track.hasTPC() * 1.); + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 7., track.hasTRD() * 1.); + registry.get(HIST("DG/TrackStat"))->Fill(track.pt(), 8., track.hasTOF() * 1.); + + registry.get(HIST("DG/etapt"))->Fill(track.eta(), track.pt(), 1.); + registry.get(HIST("DG/etapt2"))->Fill(track.eta(), track.pt(), track.isQualityTrack() * 1.); + registry.get(HIST("DG/etapt3"))->Fill(track.eta(), track.pt(), track.isGlobalTrack() * 1.); + registry.get(HIST("DG/etapt4"))->Fill(track.eta(), track.pt(), 1.); + registry.get(HIST("DG/etapt5"))->Fill(track.eta(), track.pt(), track.hasITS() * 1.); + registry.get(HIST("DG/etapt6"))->Fill(track.eta(), track.pt(), track.hasTPC() * 1.); + registry.get(HIST("DG/etapt7"))->Fill(track.eta(), track.pt(), track.hasTRD() * 1.); + registry.get(HIST("DG/etapt8"))->Fill(track.eta(), track.pt(), track.hasTOF() * 1.); + LOGF(debug, "dEdx TPC %f TOF %i %f", track.tpcSignal(), track.hasTOF(), track.hasTOF() ? track.tofSignal() : 0.); - registry.get(HIST("dEdxTPCDG"))->Fill(track.p() * track.sign(), track.tpcSignal()); - registry.get(HIST("etaptDG"))->Fill(track.eta(), track.pt()); - registry.get(HIST("IVMptTrkDG"))->Fill(ivm.M(), track.pt()); + registry.get(HIST("DG/dEdxTPC"))->Fill(track.tpcInnerParam() * track.sign(), track.tpcSignal()); if (track.hasTOF()) { - registry.get(HIST("dEdxTOFDG"))->Fill(track.pt(), track.tofSignal()); + registry.get(HIST("DG/dEdxTOF"))->Fill(track.p() * track.sign(), track.tofSignal()); } + registry.get(HIST("DG/IVMptTrk"))->Fill(ivm.M(), track.pt()); } } } } PROCESS_SWITCH(DiffQA, processMain, "Process Main", true); - // ............................................................................................................... - // Distribution of number of PV contributors for all collisions and those with empty FT0 - void processFewProng(CC const& collision, BCs const& bct0s, - aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds) - { - // count collisions - registry.get(HIST("fpStat"))->Fill(1., 1.); - registry.get(HIST("allPVC"))->Fill(collision.numContrib(), 1.); - - // check FT0 to be empty - auto bc = collision.foundBC_as(); - if (udhelpers::cleanFT0(bc, 0., 0.)) { - // only collisions with empty FT0 arrive here - registry.get(HIST("fpStat"))->Fill(2., 1.); - - // update #PV contributors in collisions with empty FT0 - registry.get(HIST("fpPVC"))->Fill(collision.numContrib(), 1.); - } - } - PROCESS_SWITCH(DiffQA, processFewProng, "Process FewProng", true); - // ............................................................................................................... // Fraction of collisions with empty FIT as function of NDtcoll void processCleanFIT1(CC const& collision, BCs const& bct0s, aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds) { - LOGF(debug, "(HIST("cleanFIT1"))->Fill(NDtcoll, isDGcandidate * 1.); + registry.get(HIST("cleanFIT1/Stat"))->Fill(NDtcoll, isDGcandidate * 1.); if (isDGcandidate) { - registry.get(HIST("cF1FV0Aamp"))->Fill(NDtcoll, ampFV0A); - registry.get(HIST("cF1FT0Aamp"))->Fill(NDtcoll, ampFT0A); - registry.get(HIST("cF1FT0Camp"))->Fill(NDtcoll, ampFT0C); - registry.get(HIST("cF1FDDAamp"))->Fill(NDtcoll, ampFDDA); - registry.get(HIST("cF1FDDCamp"))->Fill(NDtcoll, ampFDDC); + registry.get(HIST("cleanFIT1/FV0Aamp"))->Fill(NDtcoll, ampFV0A); + registry.get(HIST("cleanFIT1/FT0Aamp"))->Fill(NDtcoll, ampFT0A); + registry.get(HIST("cleanFIT1/FT0Camp"))->Fill(NDtcoll, ampFT0C); + registry.get(HIST("cleanFIT1/FDDAamp"))->Fill(NDtcoll, ampFDDA); + registry.get(HIST("cleanFIT1/FDDCamp"))->Fill(NDtcoll, ampFDDC); } } } - PROCESS_SWITCH(DiffQA, processCleanFIT1, "Process CleanFitTest1", true); + PROCESS_SWITCH(DiffQA, processCleanFIT1, "Process CleanFit1", true); // ............................................................................................................... void processCleanFIT2(CC const& collision, BCs const& bct0s, aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds) { - LOGF(debug, "(HIST("cleanFIT2"))->Fill(nMinBC, isDGcandidate * 1.); + registry.get(HIST("cleanFIT2/Stat"))->Fill(nMinBC, isDGcandidate * 1.); if (isDGcandidate) { - registry.get(HIST("cF2FV0Aamp"))->Fill(nMinBC, ampFV0A); - registry.get(HIST("cF2FT0Aamp"))->Fill(nMinBC, ampFT0A); - registry.get(HIST("cF2FT0Camp"))->Fill(nMinBC, ampFT0C); - registry.get(HIST("cF2FDDAamp"))->Fill(nMinBC, ampFDDA); - registry.get(HIST("cF2FDDCamp"))->Fill(nMinBC, ampFDDC); + registry.get(HIST("cleanFIT2/FV0Aamp"))->Fill(nMinBC, ampFV0A); + registry.get(HIST("cleanFIT2/FT0Aamp"))->Fill(nMinBC, ampFT0A); + registry.get(HIST("cleanFIT2/FT0Camp"))->Fill(nMinBC, ampFT0C); + registry.get(HIST("cleanFIT2/FDDAamp"))->Fill(nMinBC, ampFDDA); + registry.get(HIST("cleanFIT2/FDDCamp"))->Fill(nMinBC, ampFDDC); } } } - PROCESS_SWITCH(DiffQA, processCleanFIT2, "Process CleanFitTest2", true); + PROCESS_SWITCH(DiffQA, processCleanFIT2, "Process CleanFit2", true); + + // ............................................................................................................... + // Distribution of number of PV contributors for all collisions and those with empty FT0 + void processCleanFT0(CC const& collision, BCs const& bct0s, + aod::FT0s const& ft0s, aod::FV0As const& fv0as, aod::FDDs const& fdds) + { + // count collisions + registry.get(HIST("cleanFT0/Stat"))->Fill(1., 1.); + registry.get(HIST("cleanFT0/PVTracks"))->Fill(collision.numContrib(), 1., 1.); + + // check FT0 to be empty + auto bc = collision.foundBC_as(); + if (udhelpers::cleanFT0(bc, 0., 0.)) { + // only collisions with empty FT0 arrive here + registry.get(HIST("cleanFT0/Stat"))->Fill(2., 1.); + + // update #PV contributors in collisions with empty FT0 + registry.get(HIST("cleanFT0/PVTracks"))->Fill(collision.numContrib(), 2., 1.); + } + } + PROCESS_SWITCH(DiffQA, processCleanFT0, "Process CleanFIT1", true); // ............................................................................................................... void processFV0(aod::FV0As const& fv0s, BCs const&) @@ -559,25 +619,25 @@ struct DiffQA { if (fv0s.size() <= 0) { return; } - int64_t lastBCwFV0 = fv0s.begin().bc_as().globalBC(); - auto lastOrbit = lastBCwFV0 / nBCpOrbit; + auto lastOrbit = lastBCwFV0 / o2::constants::lhc::LHCMaxBunches; + for (auto fv0 : fv0s) { // side A for (size_t ind = 0; ind < fv0.channel().size(); ind++) { - registry.get(HIST("FV0A"))->Fill((fv0.channel())[ind], (fv0.amplitude())[ind]); + registry.get(HIST("FV0/FV0Aamp"))->Fill((fv0.channel())[ind], (fv0.amplitude())[ind]); } - // sequence of BCs + // length of series of empty BCs auto bc = fv0.bc_as(); int64_t aBC = bc.globalBC(); - auto aOrbit = aBC / nBCpOrbit; + auto aOrbit = aBC / o2::constants::lhc::LHCMaxBunches; auto ampA = udhelpers::FV0AmplitudeA(bc.foundFV0()); if (ampA > 0.) { // require both BCs to be in same orbit if (aOrbit == lastOrbit) - registry.get(HIST("FV0BCNUM"))->Fill((bc.globalBC() - lastBCwFV0)); + registry.get(HIST("FV0/emptyBCs"))->Fill((bc.globalBC() - lastBCwFV0)); lastBCwFV0 = aBC; lastOrbit = aOrbit; } @@ -592,47 +652,47 @@ struct DiffQA { int nc = 0; int64_t fBC = 0; // first BC with FIT activity int64_t aBC = 0; // actually processed BC - float minAmpA = 15., fAmpA = 0.; - float minAmpC = 15., fAmpC = 0.; + float minAmpA = 0., fAmpA = 0.; + float minAmpC = 0., fAmpC = 0.; int64_t lastBCwFT0 = ft0s.begin().bc_as().globalBC(); - auto lastOrbit = lastBCwFT0 / nBCpOrbit; + auto lastOrbit = lastBCwFT0 / o2::constants::lhc::LHCMaxBunches; for (auto ft0 : ft0s) { // side A for (size_t ind = 0; ind < ft0.channelA().size(); ind++) { - registry.get(HIST("FT0A"))->Fill((ft0.channelA())[ind], (ft0.amplitudeA())[ind]); + registry.get(HIST("FT0/FT0Aamp"))->Fill((ft0.channelA())[ind], (ft0.amplitudeA())[ind]); } // side C for (size_t ind = 0; ind < ft0.channelC().size(); ind++) { - registry.get(HIST("FT0C"))->Fill((ft0.channelC())[ind], (ft0.amplitudeC())[ind]); + registry.get(HIST("FT0/FT0Camp"))->Fill((ft0.channelC())[ind], (ft0.amplitudeC())[ind]); } // sequence of BCs auto bc = ft0.bc_as(); aBC = bc.globalBC(); - auto gBC = aBC % nBCpOrbit; + auto gBC = aBC % o2::constants::lhc::LHCMaxBunches; auto ampA = udhelpers::FT0AmplitudeA(bc.foundFT0()); auto ampC = udhelpers::FT0AmplitudeC(bc.foundFT0()); - // update FT0ABCNUM + // update AP2BC if (ampA > 0.) { - registry.get(HIST("FT0ABCNUM"))->Fill(gBC, 1.); + registry.get(HIST("FT0/AP2BC"))->Fill(gBC, 1.); } - // update FT0CBCNUM + // update FT0/CP2BC if (ampC > 0.) { - registry.get(HIST("FT0CBCNUM"))->Fill(gBC, 1.); + registry.get(HIST("FT0/CP2BC"))->Fill(gBC, 1.); } // update dFT0BCNUM // require both BCs to be in same orbit - auto aOrbit = aBC / nBCpOrbit; + auto aOrbit = aBC / o2::constants::lhc::LHCMaxBunches; LOGF(debug, "lastOrbit %d aOrbit %d", lastOrbit, aOrbit); if (ampA > 0. || ampC > 0.) { if (aOrbit == lastOrbit) - registry.get(HIST("dFT0BCNUM"))->Fill((aBC - lastBCwFT0)); + registry.get(HIST("FT0/emptyBCs"))->Fill((aBC - lastBCwFT0)); } else { continue; } @@ -645,7 +705,7 @@ struct DiffQA { fAmpC = ampC; } auto dBC = static_cast(aBC - fBC); - if (dBC <= ns) { + if (dBC >= 0 && dBC <= ns) { LOGF(debug, " dBC %d ampA %f ampC %f", dBC, ampA, ampC); switch (dBC) { case 0: @@ -787,29 +847,29 @@ struct DiffQA { LOGF(debug, " %d", fdds.size()); int64_t lastBCwFDD = fdds.begin().bc_as().globalBC(); - auto lastOrbit = lastBCwFDD / nBCpOrbit; + auto lastOrbit = lastBCwFDD / o2::constants::lhc::LHCMaxBunches; for (auto fdd : fdds) { // side A for (auto ind = 0; ind < 8; ind++) { - registry.get(HIST("FDDA"))->Fill(ind, (fdd.chargeA())[ind]); + registry.get(HIST("FDD/FDDAamp"))->Fill(ind, (fdd.chargeA())[ind]); } // side C for (auto ind = 0; ind < 8; ind++) { - registry.get(HIST("FDDC"))->Fill(ind, (fdd.chargeC())[ind]); + registry.get(HIST("FDD/FDDCamp"))->Fill(ind, (fdd.chargeC())[ind]); } // sequence of BCs auto bc = fdd.bc_as(); auto aBC = bc.globalBC(); - int64_t aOrbit = aBC / nBCpOrbit; + int64_t aOrbit = aBC / o2::constants::lhc::LHCMaxBunches; auto ampA = udhelpers::FDDAmplitudeA(bc.foundFDD()); auto ampC = udhelpers::FDDAmplitudeC(bc.foundFDD()); if (ampA > 0. || ampC > 0.) { // require both BCs to be in same orbit if (aOrbit == lastOrbit) - registry.get(HIST("FDDBCNUM"))->Fill((bc.globalBC() - lastBCwFDD)); + registry.get(HIST("FDD/emptyBCs"))->Fill((bc.globalBC() - lastBCwFDD)); lastBCwFDD = aBC; lastOrbit = aOrbit; } @@ -823,28 +883,28 @@ struct DiffQA { LOGF(debug, " %d", zdc.size()); // Zdc energies - registry.get(HIST("ZdcEnergies"))->Fill(0., zdc.energyZEM1()); - registry.get(HIST("ZdcEnergies"))->Fill(1., zdc.energyZEM2()); - registry.get(HIST("ZdcEnergies"))->Fill(2., zdc.energyCommonZNA()); - registry.get(HIST("ZdcEnergies"))->Fill(3., zdc.energyCommonZNC()); - registry.get(HIST("ZdcEnergies"))->Fill(4., zdc.energyCommonZPA()); - registry.get(HIST("ZdcEnergies"))->Fill(5., zdc.energyCommonZPC()); - registry.get(HIST("ZdcEnergies"))->Fill(6., (zdc.energySectorZNA())[0]); - registry.get(HIST("ZdcEnergies"))->Fill(7., (zdc.energySectorZNA())[1]); - registry.get(HIST("ZdcEnergies"))->Fill(8., (zdc.energySectorZNA())[2]); - registry.get(HIST("ZdcEnergies"))->Fill(9., (zdc.energySectorZNA())[3]); - registry.get(HIST("ZdcEnergies"))->Fill(10., (zdc.energySectorZNC())[0]); - registry.get(HIST("ZdcEnergies"))->Fill(11., (zdc.energySectorZNC())[1]); - registry.get(HIST("ZdcEnergies"))->Fill(12., (zdc.energySectorZNC())[2]); - registry.get(HIST("ZdcEnergies"))->Fill(13., (zdc.energySectorZNC())[2]); - registry.get(HIST("ZdcEnergies"))->Fill(14., (zdc.energySectorZPA())[0]); - registry.get(HIST("ZdcEnergies"))->Fill(15., (zdc.energySectorZPA())[1]); - registry.get(HIST("ZdcEnergies"))->Fill(16., (zdc.energySectorZPA())[2]); - registry.get(HIST("ZdcEnergies"))->Fill(17., (zdc.energySectorZPA())[3]); - registry.get(HIST("ZdcEnergies"))->Fill(18., (zdc.energySectorZPC())[0]); - registry.get(HIST("ZdcEnergies"))->Fill(19., (zdc.energySectorZPC())[1]); - registry.get(HIST("ZdcEnergies"))->Fill(20., (zdc.energySectorZPC())[2]); - registry.get(HIST("ZdcEnergies"))->Fill(21., (zdc.energySectorZPC())[3]); + registry.get(HIST("ZDC/Energies"))->Fill(0., zdc.energyZEM1()); + registry.get(HIST("ZDC/Energies"))->Fill(1., zdc.energyZEM2()); + registry.get(HIST("ZDC/Energies"))->Fill(2., zdc.energyCommonZNA()); + registry.get(HIST("ZDC/Energies"))->Fill(3., zdc.energyCommonZNC()); + registry.get(HIST("ZDC/Energies"))->Fill(4., zdc.energyCommonZPA()); + registry.get(HIST("ZDC/Energies"))->Fill(5., zdc.energyCommonZPC()); + registry.get(HIST("ZDC/Energies"))->Fill(6., (zdc.energySectorZNA())[0]); + registry.get(HIST("ZDC/Energies"))->Fill(7., (zdc.energySectorZNA())[1]); + registry.get(HIST("ZDC/Energies"))->Fill(8., (zdc.energySectorZNA())[2]); + registry.get(HIST("ZDC/Energies"))->Fill(9., (zdc.energySectorZNA())[3]); + registry.get(HIST("ZDC/Energies"))->Fill(10., (zdc.energySectorZNC())[0]); + registry.get(HIST("ZDC/Energies"))->Fill(11., (zdc.energySectorZNC())[1]); + registry.get(HIST("ZDC/Energies"))->Fill(12., (zdc.energySectorZNC())[2]); + registry.get(HIST("ZDC/Energies"))->Fill(13., (zdc.energySectorZNC())[2]); + registry.get(HIST("ZDC/Energies"))->Fill(14., (zdc.energySectorZPA())[0]); + registry.get(HIST("ZDC/Energies"))->Fill(15., (zdc.energySectorZPA())[1]); + registry.get(HIST("ZDC/Energies"))->Fill(16., (zdc.energySectorZPA())[2]); + registry.get(HIST("ZDC/Energies"))->Fill(17., (zdc.energySectorZPA())[3]); + registry.get(HIST("ZDC/Energies"))->Fill(18., (zdc.energySectorZPC())[0]); + registry.get(HIST("ZDC/Energies"))->Fill(19., (zdc.energySectorZPC())[1]); + registry.get(HIST("ZDC/Energies"))->Fill(20., (zdc.energySectorZPC())[2]); + registry.get(HIST("ZDC/Energies"))->Fill(21., (zdc.energySectorZPC())[3]); }; PROCESS_SWITCH(DiffQA, processZDC, "Process ZDC", true); From db025738e8d0936e3a710b437867a534df629596 Mon Sep 17 00:00:00 2001 From: Paul Date: Fri, 7 Apr 2023 17:27:32 +0200 Subject: [PATCH 2/4] avoid compilation warning because of ... comparison of integer expressions of different signedness --- EventFiltering/selectBCRange.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/EventFiltering/selectBCRange.cxx b/EventFiltering/selectBCRange.cxx index ff288657716..23793dbe56b 100644 --- a/EventFiltering/selectBCRange.cxx +++ b/EventFiltering/selectBCRange.cxx @@ -119,7 +119,7 @@ struct BCRangeSelector { auto localIter = bcIter; while (localIter.globalIndex() > 0) { --localIter; - if (localIter.globalBC() >= minBC.toLong()) { + if (localIter.globalBC() >= static_cast(minBC.toLong())) { minBCId = localIter.globalIndex(); } else { break; @@ -128,7 +128,7 @@ struct BCRangeSelector { localIter = bcIter; while (localIter.globalIndex() < bcs.size() - 1) { ++localIter; - if (localIter.globalBC() <= maxBC.toLong()) { + if (localIter.globalBC() <= static_cast(maxBC.toLong())) { maxBCId = localIter.globalIndex(); } else { break; @@ -197,7 +197,7 @@ struct BCRangeSelector { uint64_t second{bcs.iteratorAt(bcRanges[iR].second).globalBC()}; IR2.setFromLong(second); for (int i{0}; i < nToBeAddedPerRange && nToBeAdded > 0; ++i) { - if (bcRanges[iR].second < bcs.size() - 1) { + if (bcRanges[iR].second < static_cast(bcs.size() - 1)) { second = bcs.iteratorAt(bcRanges[iR].second + 1).globalBC(); IR1.setFromLong(second); if (IR1.differenceInBC(IR2) > o2::constants::lhc::LHCMaxBunches) { // protection against change of orbit in the DataFrame From f4aa31258956e40746456ae1aedbf3a6c190eb26 Mon Sep 17 00:00:00 2001 From: Paul Date: Fri, 7 Apr 2023 17:44:33 +0200 Subject: [PATCH 3/4] avoid compilation warning --- EventFiltering/selectBCRange.cxx | 2 +- PWGUD/Tasks/DGCandAnalyzer.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/EventFiltering/selectBCRange.cxx b/EventFiltering/selectBCRange.cxx index 23793dbe56b..f8f479d4fe5 100644 --- a/EventFiltering/selectBCRange.cxx +++ b/EventFiltering/selectBCRange.cxx @@ -173,7 +173,7 @@ struct BCRangeSelector { nBCselected += range.second - range.first + 1; } int nToBeAdded = std::ceil((bcs.size() - nBCselected) * fillFac); - int nToBeAddedPerRange = std::ceil(float(nToBeAdded) / bcRanges.size()); + int nToBeAddedPerRange = std::ceil(static_cast(nToBeAdded) / bcRanges.size()); LOGF(debug, "Extending ranges by %d BCs (%d per selected range)", nToBeAdded, nToBeAddedPerRange); InteractionRecord IR1, IR2; diff --git a/PWGUD/Tasks/DGCandAnalyzer.cxx b/PWGUD/Tasks/DGCandAnalyzer.cxx index a1e76cfb5d4..bc723d248dd 100644 --- a/PWGUD/Tasks/DGCandAnalyzer.cxx +++ b/PWGUD/Tasks/DGCandAnalyzer.cxx @@ -13,6 +13,7 @@ // \author Paul Buehler, paul.buehler@oeaw.ac.at // \since 06.06.2022 +#include #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" @@ -23,7 +24,6 @@ #include "PWGUD/Core/DGPIDSelector.h" #include "PWGUD/Core/UDGoodRunSelector.h" #include "PWGUD/Core/UDFSParser.h" -#include using namespace o2; using namespace o2::framework; From adb2c03f38b78b7da9cde291adff9ad02b80b91b Mon Sep 17 00:00:00 2001 From: Paul Date: Wed, 12 Apr 2023 23:33:09 +0200 Subject: [PATCH 4/4] moved UDTutorials to O2Physics/Tutorials/PWGUD --- PWGUD/Tasks/CMakeLists.txt | 15 ----------- Tutorials/CMakeLists.txt | 1 + Tutorials/PWGUD/CMakeLists.txt | 25 +++++++++++++++++++ .../PWGUD}/UDTutorial_01.cxx | 0 .../PWGUD}/UDTutorial_02a.cxx | 0 .../PWGUD}/UDTutorial_02b.cxx | 0 6 files changed, 26 insertions(+), 15 deletions(-) create mode 100644 Tutorials/PWGUD/CMakeLists.txt rename {PWGUD/Tasks => Tutorials/PWGUD}/UDTutorial_01.cxx (100%) rename {PWGUD/Tasks => Tutorials/PWGUD}/UDTutorial_02a.cxx (100%) rename {PWGUD/Tasks => Tutorials/PWGUD}/UDTutorial_02b.cxx (100%) diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index 5c242b7db90..1e0d98ed180 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -39,21 +39,6 @@ o2physics_add_dpl_workflow(dgcand-analyzer PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGCutparHolder O2Physics::UDGoodRunSelector O2Physics::DGPIDSelector O2Physics::UDFSParser COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(udtutorial-01 - SOURCES UDTutorial_01.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(udtutorial-02a - SOURCES UDTutorial_02a.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(udtutorial-02b - SOURCES UDTutorial_02b.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector - COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(upccand-analyzer SOURCES UPCCandidateAnalyzer.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase diff --git a/Tutorials/CMakeLists.txt b/Tutorials/CMakeLists.txt index d6c9440e57c..daecc4ed5c3 100644 --- a/Tutorials/CMakeLists.txt +++ b/Tutorials/CMakeLists.txt @@ -13,6 +13,7 @@ add_subdirectory(Skimming) add_subdirectory(OpenData) add_subdirectory(ML) add_subdirectory(PWGEM) +add_subdirectory(PWGUD) o2physics_add_dpl_workflow(histogram-track-selection SOURCES src/histogramTrackSelection.cxx diff --git a/Tutorials/PWGUD/CMakeLists.txt b/Tutorials/PWGUD/CMakeLists.txt new file mode 100644 index 00000000000..ca6f47a141a --- /dev/null +++ b/Tutorials/PWGUD/CMakeLists.txt @@ -0,0 +1,25 @@ +# 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. + +o2physics_add_dpl_workflow(udtutorial-01 + SOURCES UDTutorial_01.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(udtutorial-02a + SOURCES UDTutorial_02a.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(udtutorial-02b + SOURCES UDTutorial_02b.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector + COMPONENT_NAME Analysis) diff --git a/PWGUD/Tasks/UDTutorial_01.cxx b/Tutorials/PWGUD/UDTutorial_01.cxx similarity index 100% rename from PWGUD/Tasks/UDTutorial_01.cxx rename to Tutorials/PWGUD/UDTutorial_01.cxx diff --git a/PWGUD/Tasks/UDTutorial_02a.cxx b/Tutorials/PWGUD/UDTutorial_02a.cxx similarity index 100% rename from PWGUD/Tasks/UDTutorial_02a.cxx rename to Tutorials/PWGUD/UDTutorial_02a.cxx diff --git a/PWGUD/Tasks/UDTutorial_02b.cxx b/Tutorials/PWGUD/UDTutorial_02b.cxx similarity index 100% rename from PWGUD/Tasks/UDTutorial_02b.cxx rename to Tutorials/PWGUD/UDTutorial_02b.cxx