diff --git a/EventFiltering/PWGUD/diffractionFilter.cxx b/EventFiltering/PWGUD/diffractionFilter.cxx index 42405d02aef..e3df65ff91d 100644 --- a/EventFiltering/PWGUD/diffractionFilter.cxx +++ b/EventFiltering/PWGUD/diffractionFilter.cxx @@ -39,55 +39,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -// The associations between collsisions and BCs can be ambiguous. -// By default a collision is associated with the BC closest in time. -// The collision time t_coll is determined by the tracks which are used to -// reconstruct the vertex. t_coll has an uncertainty dt_coll. -// Any BC with a BC time t_BC falling within a time window of +- ndt*dt_coll -// around t_coll could potentially be the true BC. ndt is typically 4. - -template -T compatibleBCs(soa::Join::iterator const& collision, int ndt, T const& bcs) -{ - LOGF(debug, "Collision time / resolution [ns]: %f / %f", collision.collisionTime(), collision.collisionTimeRes()); - - auto bcIter = collision.bc_as(); - - // due to the filling scheme the most probably BC may not be the one estimated from the collision time - uint64_t mostProbableBC = bcIter.globalBC(); - uint64_t meanBC = mostProbableBC - std::lround(collision.collisionTime() / o2::constants::lhc::LHCBunchSpacingNS); - int deltaBC = std::ceil(collision.collisionTimeRes() / o2::constants::lhc::LHCBunchSpacingNS * ndt); - int64_t minBC = meanBC - deltaBC; - uint64_t maxBC = meanBC + deltaBC; - if (minBC < 0) { - minBC = 0; - } - - // find slice of BCs table with BC in [minBC, maxBC] - int64_t maxBCId = bcIter.globalIndex(); - int moveCount = 0; // optimize to avoid to re-create the iterator - while (bcIter != bcs.end() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { - LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); - maxBCId = bcIter.globalIndex(); - ++bcIter; - ++moveCount; - } - - bcIter.moveByIndex(-moveCount); // Move back to original position - int64_t minBCId = collision.bcId(); - while (bcIter != bcs.begin() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { - LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); - minBCId = bcIter.globalIndex(); - --bcIter; - } - - LOGF(debug, " BC range: %i (%d) - %i (%d)", minBC, minBCId, maxBC, maxBCId); - - T slice{{bcs.asArrowTable()->Slice(minBCId, maxBCId - minBCId + 1)}, (uint64_t)minBCId}; - bcs.copyIndexBindings(slice); - return slice; -} - // Run 3 struct DGFilterRun3 { @@ -100,6 +51,12 @@ struct DGFilterRun3 { // DG selector DGSelector dgSelector; + HistogramRegistry registry{ + "registry", + { + {"aftercut", "#aftercut", {HistType::kTH1F, {{7, -0.5, 6.5}}}}, + }}; + // some general Collisions and Tracks filter using CCs = soa::Join; @@ -114,7 +71,7 @@ struct DGFilterRun3 { void process(CC const& collision, BCs const& bcs, TCs& tracks, - // MFs& mfttracks, + // MFs& mfttracks, FWs& fwdtracks, aod::Zdcs& zdcs, aod::FT0s& ft0s, @@ -127,28 +84,19 @@ struct DGFilterRun3 { // obtain slice of compatible BCs auto bcRange = compatibleBCs(collision, diffCuts->NDtcoll(), bcs); - LOGF(debug, " Number of compatible BCs in +- %i dtcoll: %i", diffCuts->NDtcoll(), bcRange.size()); - - // check that there are no FIT signals in any of the compatible BCs - // Double Gap (DG) condition - auto isDGEvent = true; - for (auto& bc : bcRange) { - if (bc.has_ft0() || bc.has_fv0a() || bc.has_fdd()) { - isDGEvent = false; - break; - } - } + LOGF(info, " Number of compatible BCs in +- %i dtcoll: %i", diffCuts->NDtcoll(), bcRange.size()); - // additional cuts - if (isDGEvent) { - isDGEvent = dgSelector.IsSelected(diffCuts, collision, bc, bcRange, tracks, fwdtracks); - } + // apply DG selection + auto isDGEvent = dgSelector.IsSelected(diffCuts, collision, bc, bcRange, tracks, fwdtracks); // fill filterTable - if (isDGEvent) { + if (isDGEvent == 0) { LOGF(info, "This collision is a DG candidate!"); } - filterTable(isDGEvent); + filterTable(isDGEvent == 0); + + // update histogram + registry.get(HIST("aftercut"))->Fill(isDGEvent); } }; diff --git a/EventFiltering/PWGUD/diffractionSelectors.h b/EventFiltering/PWGUD/diffractionSelectors.h index 28a973a8f2c..431fb115a17 100644 --- a/EventFiltering/PWGUD/diffractionSelectors.h +++ b/EventFiltering/PWGUD/diffractionSelectors.h @@ -24,6 +24,9 @@ using namespace o2::framework::expressions; template bool hasGoodPID(cutHolder diffCuts, TC track); +template +T compatibleBCs(soa::Join::iterator const& collision, int ndt, T const& bcs); + // add here Selectors for different types of diffractive events // Selector for Double Gap events struct DGSelector { @@ -32,19 +35,27 @@ struct DGSelector { // Function to check if collisions passes filter template - bool IsSelected(cutHolder diffCuts, CC const& collision, BC& bc, BCs& bcRange, TCs& tracks, FWs& fwdtracks) + int IsSelected(cutHolder diffCuts, CC const& collision, BC& bc, BCs& bcRange, TCs& tracks, FWs& fwdtracks) { LOGF(debug, "Collision %f BC %i", collision.collisionTime(), bc.globalBC()); LOGF(debug, "Number of close BCs: %i", bcRange.size()); - LOGF(debug, "Number of tracks: %i", tracks.size()); + + // check that there are no FIT signals in any of the compatible BCs + // Double Gap (DG) condition + for (auto& bc : bcRange) { + if (bc.has_ft0() || bc.has_fv0a() || bc.has_fdd()) { + return 1; + } + } // Number of tracks // All tracks in vertex + LOGF(debug, "Number of tracks: %i in vertex: %i", tracks.size(), collision.numContrib()); if (collision.numContrib() != tracks.size()) { - return false; + return 2; } if (tracks.size() < diffCuts.minNTracks() || tracks.size() > diffCuts.maxNTracks()) { - return false; + return 3; } // all tracks must be global tracks @@ -70,7 +81,7 @@ struct DGSelector { goodTracks &= (nTracksWithTOFHit >= diffCuts.minNTracksWithTOFHit()); goodTracks &= (netCharge >= diffCuts.minNetCharge() && netCharge <= diffCuts.maxNetCharge()); if (!goodTracks) { - return false; + return 4; } // only tracks with good PID @@ -79,7 +90,7 @@ struct DGSelector { goodPID &= hasGoodPID(diffCuts, track); } if (!goodPID) { - return false; + return 5; } // check no activity in muon arm @@ -88,14 +99,63 @@ struct DGSelector { LOGF(debug, " %i / %f / %f / %f", muon.trackType(), muon.eta(), muon.pt(), muon.p()); } if (fwdtracks.size() > 0) { - return false; + return 6; } - // ATTENTION: currently all events are selected - return true; + // if we arrive here then the event is good! + return 0; }; }; +// The associations between collsisions and BCs can be ambiguous. +// By default a collision is associated with the BC closest in time. +// The collision time t_coll is determined by the tracks which are used to +// reconstruct the vertex. t_coll has an uncertainty dt_coll. +// Any BC with a BC time t_BC falling within a time window of +- ndt*dt_coll +// around t_coll could potentially be the true BC. ndt is typically 4. + +template +T compatibleBCs(soa::Join::iterator const& collision, int ndt, T const& bcs) +{ + LOGF(debug, "Collision time / resolution [ns]: %f / %f", collision.collisionTime(), collision.collisionTimeRes()); + + auto bcIter = collision.bc_as(); + + // due to the filling scheme the most probably BC may not be the one estimated from the collision time + uint64_t mostProbableBC = bcIter.globalBC(); + uint64_t meanBC = mostProbableBC - std::lround(collision.collisionTime() / o2::constants::lhc::LHCBunchSpacingNS); + int deltaBC = std::ceil(collision.collisionTimeRes() / o2::constants::lhc::LHCBunchSpacingNS * ndt); + int64_t minBC = meanBC - deltaBC; + uint64_t maxBC = meanBC + deltaBC; + if (minBC < 0) { + minBC = 0; + } + + // find slice of BCs table with BC in [minBC, maxBC] + int64_t maxBCId = bcIter.globalIndex(); + int moveCount = 0; // optimize to avoid to re-create the iterator + while (bcIter != bcs.end() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { + LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); + maxBCId = bcIter.globalIndex(); + ++bcIter; + ++moveCount; + } + + bcIter.moveByIndex(-moveCount); // Move back to original position + int64_t minBCId = collision.bcId(); + while (bcIter != bcs.begin() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { + LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); + minBCId = bcIter.globalIndex(); + --bcIter; + } + + LOGF(debug, " BC range: %i (%d) - %i (%d)", minBC, minBCId, maxBC, maxBCId); + + T slice{{bcs.asArrowTable()->Slice(minBCId, maxBCId - minBCId + 1)}, (uint64_t)minBCId}; + bcs.copyIndexBindings(slice); + return slice; +} + // function to check if track provides good PID information // Checks the nSigma for any particle assumption to be within limits. template diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index 0ab2e94acac..9395a638e1e 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -23,3 +23,8 @@ o2physics_add_dpl_workflow(diff-mcdatascanner SOURCES diffMCDataScanner.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2::DetectorsVertexing COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(diff-mcqa + SOURCES diffMCQA.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2::DetectorsVertexing + COMPONENT_NAME Analysis) diff --git a/PWGUD/Tasks/diffMCDataScanner.cxx b/PWGUD/Tasks/diffMCDataScanner.cxx index 30fbc7c70c3..7fd5008e944 100644 --- a/PWGUD/Tasks/diffMCDataScanner.cxx +++ b/PWGUD/Tasks/diffMCDataScanner.cxx @@ -10,9 +10,6 @@ // or submit itself to any jurisdiction. /// /// \brief -/// ATTENTION Nov. 2021: -/// MFT is not implemented yet and can not be used -/// releated code is commented /// \author Paul Buehler, paul.buehler@oeaw.ac.at /// \since 01.10.2021 @@ -26,70 +23,20 @@ void customize(std::vector& workflowOptions) workflowOptions.push_back(ConfigParamSpec{"runCase", VariantType::Int, 0, {"runCase: 0 - histos, 1 - mcTruth, else - tree"}}); } -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" - -#include "PWGUD/DataModel/McPIDTable.h" - -#include "CommonConstants/LHCConstants.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/PID/PIDResponse.h" -#include "Common/Core/MC.h" +#include "PWGUD/Tasks/diffMCHelpers.h" using namespace o2::framework::expressions; -template -T getCompatibleBCs(soa::Join::iterator const& collision, T const& bcs) -{ - LOGF(debug, "Collision time / resolution [ns]: %f / %f", collision.collisionTime(), collision.collisionTimeRes()); - - auto bcIter = collision.bc_as(); - - // due to the filling scheme the most probably BC may not be the one estimated from the collision time - uint64_t mostProbableBC = bcIter.globalBC(); - uint64_t meanBC = mostProbableBC - std::lround(collision.collisionTime() / o2::constants::lhc::LHCBunchSpacingNS); - int deltaBC = std::ceil(collision.collisionTimeRes() / o2::constants::lhc::LHCBunchSpacingNS * 4); - int64_t minBC = meanBC - deltaBC; - uint64_t maxBC = meanBC + deltaBC; - if (minBC < 0) { - minBC = 0; - } - - // find slice of BCs table with BC in [minBC, maxBC] - int64_t maxBCId = bcIter.globalIndex(); - int moveCount = 0; // optimize to avoid to re-create the iterator - while (bcIter != bcs.end() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { - LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); - maxBCId = bcIter.globalIndex(); - ++bcIter; - ++moveCount; - } - - bcIter.moveByIndex(-moveCount); // Move back to original position - int64_t minBCId = collision.bcId(); - while (bcIter != bcs.begin() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { - LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); - minBCId = bcIter.globalIndex(); - --bcIter; - } - - LOGF(debug, " BC range: %i (%d) - %i (%d)", minBC, minBCId, maxBC, maxBCId); - - T slice{{bcs.asArrowTable()->Slice(minBCId, maxBCId - minBCId + 1)}, (uint64_t)minBCId}; - bcs.copyIndexBindings(slice); - return slice; -} - // Loop over collisions // find for each collision the number of compatible BCs +// runCase = 1 struct CompatibleBCs { using CCs = soa::Join; using CC = CCs::iterator; void process(CC const& collision, aod::BCs const& bcs) { - auto bcSlice = getCompatibleBCs(collision, bcs); + auto bcSlice = getMCCompatibleBCs(collision, 4, bcs); LOGF(debug, " Number of possible BCs: %i", bcSlice.size()); for (auto& bc : bcSlice) { LOGF(debug, " This collision may belong to BC %lld", bc.globalBC()); @@ -97,8 +44,8 @@ struct CompatibleBCs { } }; -// Loop over collisions // Fill histograms with collision and compatible BCs related information +// runCase = 1 struct collisionsInfo { int cnt = 0; HistogramRegistry registry{ @@ -107,7 +54,7 @@ struct collisionsInfo { {"timeResolution", "#timeResolution", {HistType::kTH1F, {{200, 0., 1.E3}}}}, {"numberBCs", "#numberBCs", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, {"DGCandidate", "#DGCandidate", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, - {"numberTracks", "#numberTracks", {HistType::kTH1F, {{301, -0.5, 300.5}}}}, + {"numberTracks", "#numberTracks", {HistType::kTH1F, {{5001, -0.5, 5000.5}}}}, {"numberVtxTracks", "#numberVtxTracks", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, {"numberGlobalTracks", "#numberGlobalTracks", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, {"netCharge", "#netCharge", {HistType::kTH1F, {{3, -1.5, 1.5}}}}, @@ -119,6 +66,8 @@ struct collisionsInfo { {"etaFWDDG", "#etaFWDDG", {HistType::kTH1F, {{100, -5.0, 5.0}}}}, // {"globalVsMFTAll", "#globalVsMFTAll", {HistType::kTH2F, {{21, -0.5, 20.5}, {21, -0.5, 20.5}}}}, // {"globalVsMFTDG", "#globalVsMFTDG", {HistType::kTH2F, {{21, -0.5, 20.5}, {21, -0.5, 20.5}}}}, + {"VtxvsTracks", "#VtxvsTracks", {HistType::kTH2F, {{101, -0.5, 100.5}, {5001, -0.5, 5000.5}}}}, + {"VtxvsGlobalTracks", "#VtxvsGlobalTracks", {HistType::kTH2F, {{101, -0.5, 100.5}, {101, -0.5, 100.5}}}}, }}; void init(o2::framework::InitContext&) @@ -139,6 +88,10 @@ struct collisionsInfo { // registry.get(HIST("globalVsMFTAll"))->GetYaxis()->SetTitle("Number of MFT tracks"); // registry.get(HIST("globalVsMFTDG"))->GetXaxis()->SetTitle("Number of global tracks"); // registry.get(HIST("globalVsMFTDG"))->GetYaxis()->SetTitle("Number of MFT tracks"); + registry.get(HIST("VtxvsTracks"))->GetXaxis()->SetTitle("Number of Vtx tracks"); + registry.get(HIST("VtxvsTracks"))->GetYaxis()->SetTitle("Number of tracks"); + registry.get(HIST("VtxvsGlobalTracks"))->GetXaxis()->SetTitle("Number of Vtx tracks"); + registry.get(HIST("VtxvsGlobalTracks"))->GetYaxis()->SetTitle("Number of global tracks"); } using CCs = soa::Join; @@ -154,7 +107,7 @@ struct collisionsInfo { { // obtain slice of compatible BCs - auto bcSlice = getCompatibleBCs(collision, bct0s); + auto bcSlice = getMCCompatibleBCs(collision, 4, bct0s); LOGF(debug, " Number of compatible BCs: %i", bcSlice.size()); registry.get(HIST("numberBCs"))->Fill(bcSlice.size()); @@ -166,19 +119,23 @@ struct collisionsInfo { break; } } + + // global tracks + Partition goodTracks = aod::track::isGlobalTrack > uint8_t(0); + goodTracks.bindTable(tracks); + int cntGlobal = goodTracks.size(); + // count tracks int cntAll = 0; - int cntGlobal = 0; int netCharge = 0; for (auto track : tracks) { cntAll++; netCharge += track.sign(); - if (track.isGlobalTrack()) { - cntGlobal++; - } } // netCharge /= abs(netCharge); + LOGF(info, "Global tracks: %i", cntGlobal); + if (isDGcandidate) { LOGF(info, " This is a DG candidate with %d tracks and %d net charge.", tracks.size(), netCharge); registry.get(HIST("DGCandidate"))->Fill(1.); @@ -196,6 +153,8 @@ struct collisionsInfo { // registry.get(HIST("numberMFTTracks"))->Fill(mfttracks.size()); registry.get(HIST("numberFWDTracks"))->Fill(fwdtracks.size()); // registry.get(HIST("globalVsMFTAll"))->Fill(cntGlobal, mfttracks.size()); + registry.get(HIST("VtxvsTracks"))->Fill(collision.numContrib(), cntAll); + registry.get(HIST("VtxvsGlobalTracks"))->Fill(collision.numContrib(), cntGlobal); // loop over MFT tracks // LOGF(debug, "MFT tracks: %i", mfttracks.size()); @@ -231,6 +190,7 @@ struct collisionsInfo { // Loop over BCs // check aliases, selection, and FIT signals per BC +// runCase = 1 struct BCInfo { int cnt = 0; HistogramRegistry registry{ @@ -303,6 +263,7 @@ struct BCInfo { // Loop over tracks // Make histograms with track type and time resolution +// runCase = 1 struct TrackTypes { HistogramRegistry registry{ "registry", @@ -375,6 +336,7 @@ struct TrackTypes { }; // MCTruth tracks +// runCase = 2 struct MCTracks { using CCs = soa::Join; @@ -431,6 +393,7 @@ struct MCTracks { }; // TPC nSigma +// runCase = 3 struct TPCnSigma { Produces nSigmas; diff --git a/PWGUD/Tasks/diffMCHelpers.h b/PWGUD/Tasks/diffMCHelpers.h new file mode 100644 index 00000000000..6bdd661616d --- /dev/null +++ b/PWGUD/Tasks/diffMCHelpers.h @@ -0,0 +1,94 @@ +// 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 +/// \author Paul Buehler, paul.buehler@oeaw.ac.at +/// \since 01.10.2021 + +#ifndef O2_ANALYSIS_DIFFRACTION_MCHELPER_H_ +#define O2_ANALYSIS_DIFFRACTION_MCHELPER_H_ + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "PWGUD/DataModel/McPIDTable.h" +#include "CommonConstants/LHCConstants.h" +#include "Common/CCDB/TriggerAliases.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/Core/PID/PIDResponse.h" +#include "Common/Core/MC.h" + +using namespace o2; +using namespace o2::framework; + +// The associations between collisions and BCs can be ambiguous. +// By default a collision is associated with the BC closest in time. +// The collision time t_coll is determined by the tracks which are used to +// reconstruct the vertex. t_coll has an uncertainty dt_coll. +// Any BC with a BC time t_BC falling within a time window of +- ndt*dt_coll +// around t_coll could potentially be the true BC. ndt is typically 4. + +template +T getMCCompatibleBCs(soa::Join::iterator const& collision, int ndt, T const& bcs) +{ + LOGF(debug, "Collision time / resolution [ns]: %f / %f", collision.collisionTime(), collision.collisionTimeRes()); + + auto bcIter = collision.bc_as(); + + // due to the filling scheme the most probably BC may not be the one estimated from the collision time + uint64_t mostProbableBC = bcIter.globalBC(); + uint64_t meanBC = mostProbableBC - std::lround(collision.collisionTime() / o2::constants::lhc::LHCBunchSpacingNS); + int deltaBC = std::ceil(collision.collisionTimeRes() / o2::constants::lhc::LHCBunchSpacingNS * ndt); + int64_t minBC = meanBC - deltaBC; + uint64_t maxBC = meanBC + deltaBC; + if (minBC < 0) { + minBC = 0; + } + + // find slice of BCs table with BC in [minBC, maxBC] + int64_t maxBCId = bcIter.globalIndex(); + int moveCount = 0; // optimize to avoid to re-create the iterator + while (bcIter != bcs.end() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { + LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); + maxBCId = bcIter.globalIndex(); + ++bcIter; + ++moveCount; + } + + bcIter.moveByIndex(-moveCount); // Move back to original position + int64_t minBCId = collision.bcId(); + while (bcIter != bcs.begin() && bcIter.globalBC() <= maxBC && (int64_t)bcIter.globalBC() >= minBC) { + LOGF(debug, "Table id %d BC %llu", bcIter.globalIndex(), bcIter.globalBC()); + minBCId = bcIter.globalIndex(); + --bcIter; + } + + LOGF(debug, " BC range: %i (%d) - %i (%d)", minBC, minBCId, maxBC, maxBCId); + + T slice{{bcs.asArrowTable()->Slice(minBCId, maxBCId - minBCId + 1)}, (uint64_t)minBCId}; + bcs.copyIndexBindings(slice); + return slice; +} + +// In PYTHIA a central diffractive produced (CD) particle hs the ID +// 9900110. Check the particles of a MC event to contain a CD particle. +template +bool isMcCDE(T MCparts) +{ + for (auto mcpart : MCparts) { + if (mcpart.pdgCode() == 9900110) { + return true; + } + } + return false; +} + +#endif // O2_ANALYSIS_DIFFRACTION_MCHELPER_H_ diff --git a/PWGUD/Tasks/diffMCQA.cxx b/PWGUD/Tasks/diffMCQA.cxx new file mode 100644 index 00000000000..de0ec8fb1a3 --- /dev/null +++ b/PWGUD/Tasks/diffMCQA.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 +/// \author Paul Buehler, paul.buehler@oeaw.ac.at +/// \since 01.10.2021 + +#include "PWGUD/Tasks/diffMCHelpers.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct DiffQA { + + // define histograms + HistogramRegistry registry{ + "registry", + { + {"timeResolution", "#timeResolution", {HistType::kTH1F, {{200, 0., 1.E3}}}}, + {"tResvsrTOFTracks", "#tResvsrTOFTracks", {HistType::kTH2F, {{200, 0., 1.E3}, {101, -0.01, 1.01}}}}, + {"numberBCs", "#numberBCs", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, + {"numberTracks", "#numberTracks", {HistType::kTH1F, {{3001, -0.5, 3000.5}}}}, + {"numberVtxTracks", "#numberVtxTracks", {HistType::kTH1F, {{151, -0.5, 150.5}}}}, + {"numberGlobalTracks", "#numberGlobalTracks", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, + {"numberFWDTracks", "#numberFWDTracks", {HistType::kTH1F, {{21, -0.5, 20.5}}}}, + {"VtxvsGlobalTracks", "#VtxvsGlobalTracks", {HistType::kTH2F, {{101, -0.5, 100.5}, {101, -0.5, 100.5}}}}, + {"CDEresolution", "#CDEresolution", {HistType::kTH2F, {{2, 0., 2.}, {2, 0., 2.}}}}, + {"CDEtResvsrTOFTracks", "#CDEtResvsrTOFTracks", {HistType::kTH2F, {{200, 0., 1.E3}, {101, -0.01, 1.01}}}}, + {"CDEnumberBCs", "#CDEnumberBCs", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, + {"CDEnumberTracks", "#CDEnumberTracks", {HistType::kTH1F, {{3001, -0.5, 3000.5}}}}, + {"CDEnumberVtxTracks", "#CDEnumberVtxTracks", {HistType::kTH1F, {{151, -0.5, 150.5}}}}, + {"CDEnumberGlobalTracks", "#CDEnumberGlobalTracks", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, + {"CDEnumberFWDTracks", "#CDEnumberFWDTracks", {HistType::kTH1F, {{21, -0.5, 20.5}}}}, + {"CDEVtxvsGlobalTracks", "#CDEVtxvsGlobalTracks", {HistType::kTH2F, {{101, -0.5, 100.5}, {101, -0.5, 100.5}}}}, + {"CDEtimeResolution", "#CDEtimeResolution", {HistType::kTH1F, {{200, 0., 1.E3}}}}, + }}; + + using CCs = soa::Join; + using CC = CCs::iterator; + using BCs = soa::Join; + using TCs = soa::Join; + using FWs = aod::FwdTracks; + + void process(CC const& collision, BCs const& bct0s, + TCs& tracks, FWs& fwdtracks, aod::FT0s& ft0s, aod::FV0As& fv0as, aod::FDDs& fdds, + aod::McCollisions& McCols, aod::McParticles& McParts) + { + + // is this a central diffractive events? + auto MCCol = collision.mcCollision(); + auto MCPartSlice = McParts.sliceBy(aod::mcparticle::mcCollisionId, MCCol.globalIndex()); + auto isCDE = isMcCDE(MCPartSlice); + + // obtain slice of compatible BCs + auto bcSlice = getMCCompatibleBCs(collision, 4, bct0s); + LOGF(debug, " Number of compatible BCs: %i", bcSlice.size()); + + // global tracks + Partition goodTracks = aod::track::isGlobalTrack > uint8_t(0); + goodTracks.bindTable(tracks); + + // number of gobal tracks with TOF hit + float rtrwTOF = 0.; + for (auto& track : goodTracks) { + if (track.hasTOF()) { + rtrwTOF += 1.; + } + } + rtrwTOF /= goodTracks.size(); + + // update histograms + if (isCDE) { + registry.get(HIST("CDEtimeResolution"))->Fill(collision.collisionTimeRes()); + registry.get(HIST("CDEtResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rtrwTOF); + registry.get(HIST("CDEnumberBCs"))->Fill(bcSlice.size()); + registry.get(HIST("CDEnumberTracks"))->Fill(tracks.size()); + registry.get(HIST("CDEnumberVtxTracks"))->Fill(collision.numContrib()); + registry.get(HIST("CDEnumberGlobalTracks"))->Fill(goodTracks.size()); + registry.get(HIST("CDEnumberFWDTracks"))->Fill(fwdtracks.size()); + registry.get(HIST("CDEVtxvsGlobalTracks"))->Fill(collision.numContrib(), goodTracks.size()); + } else { + registry.get(HIST("timeResolution"))->Fill(collision.collisionTimeRes()); + registry.get(HIST("tResvsrTOFTracks"))->Fill(collision.collisionTimeRes(), rtrwTOF); + registry.get(HIST("numberBCs"))->Fill(bcSlice.size()); + registry.get(HIST("numberTracks"))->Fill(tracks.size()); + registry.get(HIST("numberVtxTracks"))->Fill(collision.numContrib()); + registry.get(HIST("numberGlobalTracks"))->Fill(goodTracks.size()); + registry.get(HIST("numberFWDTracks"))->Fill(fwdtracks.size()); + registry.get(HIST("VtxvsGlobalTracks"))->Fill(collision.numContrib(), goodTracks.size()); + } + + // is it a DG candidate? + // DG = no FIT signal in compatible BCs + // & number of forward tracks = 0 + // & number of vertex tracks <= ntr + // & number of global tracks <= ntr + bool isDGcandidate = true; + int ntr = 2; + + // no FIT signal in compatible BCs + for (auto& bc : bcSlice) { + if (bc.has_ft0() || bc.has_fv0a() || bc.has_fdd()) { + isDGcandidate = false; + continue; + } + } + + // number of forward tracks = 0 + isDGcandidate &= (fwdtracks.size() == 0); + + // number of vertex tracks <= n + isDGcandidate &= (collision.numContrib() <= ntr); + + // number of global tracks <= ntr + isDGcandidate &= (goodTracks.size() <= ntr); + + if (isCDE) { + if (isDGcandidate) { + registry.get(HIST("CDEresolution"))->Fill(0.5, 0.5); + } else { + registry.get(HIST("CDEresolution"))->Fill(0.5, 1.5); + } + } else { + if (isDGcandidate) { + registry.get(HIST("CDEresolution"))->Fill(1.5, 0.5); + } else { + registry.get(HIST("CDEresolution"))->Fill(1.5, 1.5); + } + } + }; +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"diffqa"}), + }; +}