diff --git a/PWGCF/DataModel/DptDptFiltered.h b/PWGCF/DataModel/DptDptFiltered.h new file mode 100644 index 00000000000..0b1797b0440 --- /dev/null +++ b/PWGCF/DataModel/DptDptFiltered.h @@ -0,0 +1,68 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +#ifndef O2_ANALYSIS_DPTDPTFILTERED_H +#define O2_ANALYSIS_DPTDPTFILTERED_H + +#include "Framework/ASoA.h" +#include "Framework/AnalysisDataModel.h" + +namespace o2 +{ +namespace aod +{ +/* we have to change from int to bool when bool columns work properly */ +namespace dptdptfilter +{ +DECLARE_SOA_COLUMN(DptDptCFCollisionAccepted, collisionaccepted, bool); //! If the collision/event has been accepted or not +DECLARE_SOA_COLUMN(DptDptCFCollisionCentMult, centmult, float); //! The centrality/multiplicity pecentile +} // namespace dptdptfilter +DECLARE_SOA_TABLE(DptDptCFAcceptedCollisions, "AOD", "DPTDPTCFACCCOLL", //! Accepted reconstructed collisions/events filtered table + o2::soa::Index<>, + collision::BCId, + collision::PosZ, + dptdptfilter::DptDptCFCollisionAccepted, + dptdptfilter::DptDptCFCollisionCentMult); +using DptDptCFAcceptedCollision = DptDptCFAcceptedCollisions::iterator; +DECLARE_SOA_TABLE(DptDptCFAcceptedTrueCollisions, "AOD", "DPTCFACCGENCOLL", //! Accepted generated collisions/events filtered table + o2::soa::Index<>, + collision::BCId, + mccollision::PosZ, + dptdptfilter::DptDptCFCollisionAccepted, + dptdptfilter::DptDptCFCollisionCentMult); +using DptDptCFAcceptedTrueCollision = DptDptCFAcceptedTrueCollisions::iterator; +namespace dptdptfilter +{ +DECLARE_SOA_INDEX_COLUMN(DptDptCFAcceptedCollision, event); //! Reconstructed collision/event +DECLARE_SOA_INDEX_COLUMN(DptDptCFAcceptedTrueCollision, mcevent); //! Generated collision/event +DECLARE_SOA_COLUMN(TrackacceptedAsOne, trackacceptedasone, bool); //! Track accepted as type one +DECLARE_SOA_COLUMN(TrackacceptedAsTwo, trackacceptedastwo, bool); //! Track accepted as type two +DECLARE_SOA_COLUMN(Pt, pt, float); //! The track transverse momentum +DECLARE_SOA_COLUMN(Eta, eta, float); //! The track pseudorapidity +DECLARE_SOA_COLUMN(Phi, phi, float); //! The track azimuthal angle +} // namespace dptdptfilter +DECLARE_SOA_TABLE(ScannedTracks, "AOD", "SCANNEDTRACKS", //! The reconstructed tracks filtered table + dptdptfilter::DptDptCFAcceptedCollisionId, + dptdptfilter::TrackacceptedAsOne, + dptdptfilter::TrackacceptedAsTwo, + dptdptfilter::Pt, + dptdptfilter::Eta, + dptdptfilter::Phi); +DECLARE_SOA_TABLE(ScannedTrueTracks, "AOD", "SCANTRUETRACKS", //! The generated particles filtered table + dptdptfilter::DptDptCFAcceptedTrueCollisionId, + dptdptfilter::TrackacceptedAsOne, + dptdptfilter::TrackacceptedAsTwo, + dptdptfilter::Pt, + dptdptfilter::Eta, + dptdptfilter::Phi); +} // namespace aod +} // namespace o2 + +#endif // O2_ANALYSIS_DPTDPTFILTERED_H \ No newline at end of file diff --git a/PWGCF/TableProducer/CMakeLists.txt b/PWGCF/TableProducer/CMakeLists.txt index d43fd0f86b5..89e55e89743 100644 --- a/PWGCF/TableProducer/CMakeLists.txt +++ b/PWGCF/TableProducer/CMakeLists.txt @@ -13,3 +13,8 @@ o2physics_add_dpl_workflow(filter-correlations SOURCES filterCorrelations.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(dptdpt-filter + SOURCES dptdptfilter.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore + COMPONENT_NAME Analysis) diff --git a/PWGCF/TableProducer/dptdptfilter.cxx b/PWGCF/TableProducer/dptdptfilter.cxx new file mode 100644 index 00000000000..de46fb1e34c --- /dev/null +++ b/PWGCF/TableProducer/dptdptfilter.cxx @@ -0,0 +1,935 @@ +// 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. + +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Centrality.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/Core/MC.h" +#include "Common/Core/PID/PIDResponse.h" +#include "PWGCF/Core/AnalysisConfigurableCuts.h" +#include "PWGCF/DataModel/DptDptFiltered.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Framework/runDataProcessing.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "dptdptfilter.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::soa; +using namespace o2::framework::expressions; +using namespace o2::analysis; + +#define DPTDPTFILTERLOGCOLLISIONS debug +#define DPTDPTFILTERLOGTRACKS debug + +namespace o2::analysis::dptdptfilter +{ +using DptDptFullTracksPID = soa::Join; +using DptDptFullTracksPIDDetLevel = soa::Join; +using DptDptFullTracks = soa::Join; +using DptDptFullTracksDetLevel = soa::Join; + +/// \enum MatchRecoGenSpecies +/// \brief The species considered by the matching tast +enum MatchRecoGenSpecies { + kDptDptCharged = 0, ///< charged particle/track + kDptDptElectron, ///< electron + kDptDptMuon, ///< muon + kDptDptPion, ///< pion + kDptDptKaon, ///< kaon + kDptDptProton, ///< proton + kDptDptNoOfSpecies, ///< the number of considered species + kWrongSpecies = -1 +}; + +const char* speciesName[kDptDptNoOfSpecies] = {"h", "e", "mu", "pi", "ka", "p"}; + +const char* speciesTitle[kDptDptNoOfSpecies] = {"", "e", "#mu", "#pi", "K", "p"}; + +//============================================================================================ +// The DptDptFilter output objects +//============================================================================================ +std::string fTaskConfigurationString = "PendingToConfigure"; +TH1F* fhCentMultB = nullptr; +TH1F* fhCentMultA = nullptr; +TH1F* fhVertexZB = nullptr; +TH1F* fhVertexZA = nullptr; +TH1F* fhPB = nullptr; +TH1F* fhPA[kDptDptNoOfSpecies] = {nullptr}; +TH1F* fhPtB = nullptr; +TH1F* fhPtA[kDptDptNoOfSpecies] = {nullptr}; +TH1F* fhPtPosB = nullptr; +TH1F* fhPtPosA[kDptDptNoOfSpecies] = {nullptr}; +TH1F* fhPtNegB = nullptr; +TH1F* fhPtNegA[kDptDptNoOfSpecies] = {nullptr}; + +TH1F* fhEtaB = nullptr; +TH1F* fhEtaA = nullptr; + +TH1F* fhPhiB = nullptr; +TH1F* fhPhiA = nullptr; + +TH1F* fhDCAxyB = nullptr; +TH1F* fhDCAxyA = nullptr; +TH1F* fhFineDCAxyA = nullptr; +TH1F* fhDCAzB = nullptr; +TH1F* fhDCAzA = nullptr; +TH1F* fhFineDCAzA = nullptr; + +TH1F* fhTrueCentMultB = nullptr; +TH1F* fhTrueCentMultA = nullptr; +TH1F* fhTrueVertexZB = nullptr; +TH1F* fhTrueVertexZA = nullptr; +TH1F* fhTruePB = nullptr; +TH1F* fhTruePA[kDptDptNoOfSpecies] = {nullptr}; +TH1F* fhTruePtB = nullptr; +TH1F* fhTruePtA[kDptDptNoOfSpecies] = {nullptr}; +TH1F* fhTruePtPosB = nullptr; +TH1F* fhTruePtPosA[kDptDptNoOfSpecies] = {nullptr}; +TH1F* fhTruePtNegB = nullptr; +TH1F* fhTruePtNegA[kDptDptNoOfSpecies] = {nullptr}; + +TH1F* fhTrueEtaB = nullptr; +TH1F* fhTrueEtaA = nullptr; + +TH1F* fhTruePhiB = nullptr; +TH1F* fhTruePhiA = nullptr; + +TH1F* fhTrueDCAxyB = nullptr; +TH1F* fhTrueDCAxyA = nullptr; +TH1F* fhTrueDCAzB = nullptr; +TH1F* fhTrueDCAxyBid = nullptr; +TH1F* fhTrueDCAzA = nullptr; +} // namespace o2::analysis::dptdptfilter + +using namespace dptdptfilter; + +struct DptDptFilter { + Configurable cfgTrackType{"trktype", 1, "Type of selected tracks: 0 = no selection, 1 = global tracks FB96, 3 = Run3 tracks. Default 1"}; + Configurable cfgCentMultEstimator{"centmultestimator", "V0M", "Centrality/multiplicity estimator detector: V0M, NOCM: none. Default V0M"}; + Configurable cfgSystem{"syst", "PbPb", "System: pp, PbPb, Pbp, pPb, XeXe. Default PbPb"}; + Configurable cfgDataType{"datatype", "data", "Data type: data, MC, FastMC, OnTheFlyMC. Default data"}; + Configurable cfgBinning{"binning", + {28, -7.0, 7.0, 18, 0.2, 2.0, 16, -0.8, 0.8, 72, 0.5}, + "triplets - nbins, min, max - for z_vtx, pT, eta and phi, binning plus bin fraction of phi origin shift"}; + Configurable cfgTraceDCAOutliers{"trackdcaoutliers", {false, 0.0, 0.0}, "Track the generator level DCAxy outliers: false/true, low dcaxy, up dcaxy. Default {false,0.0,0.0}"}; + Configurable cfgTraceOutOfSpeciesParticles{"trackoutparticles", false, "Track the particles which are not e,mu,pi,K,p: false/true. Default false"}; + Configurable cfgRecoIdMethod{"recoidmethod", 0, "Method for identifying reconstructed tracks: 0 No PID, 1 PID, 2 mcparticle. Default 0"}; + Configurable cfgTrackSelection{"tracksel", {false, false, 0, 70, 0.8, 2.4, 3.2}, "Track selection: {useit: true/false, ongen: true/false, tpccls, tpcxrws, tpcxrfc, dcaxy, dcaz}. Default {false,0.70.0.8,2.4,3.2}"}; + Configurable cfgTraceCollId0{"tracecollid0", false, "Trace particles in collisions id 0. Default false"}; + + OutputObj fOutput{"DptDptFilterGlobalInfo", OutputObjHandlingPolicy::AnalysisObject}; + + Produces acceptedcollisions; + Produces scannedtracks; + Produces acceptedtrueevents; + Produces scannedtruetracks; + + template + void fillTrackHistosBeforeSelection(TrackObject const& track) + { + using namespace dptdptfilter; + + fhPB->Fill(track.p()); + fhPtB->Fill(track.pt()); + fhEtaB->Fill(track.eta()); + fhPhiB->Fill(track.phi()); + if (track.sign() > 0) { + fhPtPosB->Fill(track.pt()); + } else { + fhPtNegB->Fill(track.pt()); + } + fhDCAxyB->Fill(track.dcaXY()); + fhDCAzB->Fill(track.dcaZ()); + } + + template + void fillTrackHistosAfterSelection(TrackObject const& track, MatchRecoGenSpecies sp) + { + using namespace dptdptfilter; + + /* the charged species should have been called first so avoid double counting */ + if (sp == kDptDptCharged) { + fhEtaA->Fill(track.eta()); + fhPhiA->Fill(track.phi()); + fhDCAxyA->Fill(track.dcaXY()); + fhDCAzA->Fill(track.dcaZ()); + if (track.dcaXY() < 1.0) { + fhFineDCAxyA->Fill(track.dcaXY()); + } + if (track.dcaZ() < 1.0) { + fhFineDCAzA->Fill(track.dcaZ()); + } + } + fhPA[sp]->Fill(track.p()); + fhPtA[sp]->Fill(track.pt()); + if (track.sign() > 0) { + fhPtPosA[sp]->Fill(track.pt()); + } else { + fhPtNegA[sp]->Fill(track.pt()); + } + } + + template + void fillParticleHistosBeforeSelection(ParticleObject const& particle, MCCollisionObject const& collision, float charge) + { + using namespace dptdptfilter; + + fhTruePB->Fill(particle.p()); + fhTruePtB->Fill(particle.pt()); + fhTrueEtaB->Fill(particle.eta()); + fhTruePhiB->Fill(particle.phi()); + if (charge > 0) { + fhTruePtPosB->Fill(particle.pt()); + } else if (charge < 0) { + fhTruePtNegB->Fill(particle.pt()); + } + + float dcaxy = TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + + (particle.vy() - collision.posY()) * (particle.vy() - collision.posY())); + if (traceDCAOutliers.mDoIt and (traceDCAOutliers.mLowValue < dcaxy) and (dcaxy < traceDCAOutliers.mUpValue)) { + fhTrueDCAxyBid->Fill(TString::Format("%d", particle.pdgCode()).Data(), 1.0); + } + + fhTrueDCAxyB->Fill(TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + + (particle.vy() - collision.posY()) * (particle.vy() - collision.posY()))); + fhTrueDCAzB->Fill((particle.vz() - collision.posZ())); + } + + template + void fillParticleHistosAfterSelection(ParticleObject const& particle, MCCollisionObject const& collision, float charge, MatchRecoGenSpecies sp) + { + using namespace dptdptfilter; + + /* the charged species should have been called first so avoid double counting */ + if (sp == kDptDptCharged) { + fhTrueEtaA->Fill(particle.eta()); + fhTruePhiA->Fill(particle.phi()); + float dcaxy = TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + + (particle.vy() - collision.posY()) * (particle.vy() - collision.posY())); + if (traceDCAOutliers.mDoIt and (traceDCAOutliers.mLowValue < dcaxy) and (dcaxy < traceDCAOutliers.mUpValue)) { + LOGF(info, "DCAxy outlier: Particle with index %d and pdg code %d assigned to MC collision %d, pT: %f, phi: %f, eta: %f", + particle.globalIndex(), particle.pdgCode(), particle.mcCollisionId(), particle.pt(), particle.phi(), particle.eta()); + LOGF(info, " With status %d and flags %0X", particle.statusCode(), particle.flags()); + } + + fhTrueDCAxyA->Fill(TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + + (particle.vy() - collision.posY()) * (particle.vy() - collision.posY()))); + fhTrueDCAzA->Fill((particle.vz() - collision.posZ())); + } + fhTruePA[sp]->Fill(particle.p()); + fhTruePtA[sp]->Fill(particle.pt()); + if (charge > 0) { + fhTruePtPosA[sp]->Fill(particle.pt()); + } else { + fhTruePtNegA[sp]->Fill(particle.pt()); + } + } + + template + inline MatchRecoGenSpecies IdentifyTrack(TrackObject const& track) + { + using namespace o2::analysis::dptdptfilter; + + float nsigmas[kDptDptNoOfSpecies]; + if (track.p() < 0.8) { + nsigmas[kDptDptCharged] = 999.0f; + nsigmas[kDptDptElectron] = track.tpcNSigmaEl(); + nsigmas[kDptDptMuon] = track.tpcNSigmaMu(); + nsigmas[kDptDptPion] = track.tpcNSigmaPi(); + nsigmas[kDptDptKaon] = track.tpcNSigmaKa(); + nsigmas[kDptDptProton] = track.tpcNSigmaPr(); + } else { + /* introduce require TOF flag */ + if (track.hasTOF()) { + nsigmas[kDptDptCharged] = 999.0f; + nsigmas[kDptDptElectron] = sqrtf(track.tpcNSigmaEl() * track.tpcNSigmaEl() + track.tofNSigmaEl() * track.tofNSigmaEl()); + nsigmas[kDptDptMuon] = sqrtf(track.tpcNSigmaMu() * track.tpcNSigmaMu() + track.tofNSigmaMu() * track.tofNSigmaMu()); + nsigmas[kDptDptPion] = sqrtf(track.tpcNSigmaPi() * track.tpcNSigmaPi() + track.tofNSigmaPi() * track.tofNSigmaPi()); + nsigmas[kDptDptKaon] = sqrtf(track.tpcNSigmaKa() * track.tpcNSigmaKa() + track.tofNSigmaKa() * track.tofNSigmaKa()); + nsigmas[kDptDptProton] = sqrtf(track.tpcNSigmaPr() * track.tpcNSigmaPr() + track.tofNSigmaPr() * track.tofNSigmaPr()); + } else { + nsigmas[kDptDptCharged] = 999.0f; + nsigmas[kDptDptElectron] = track.tpcNSigmaEl(); + nsigmas[kDptDptMuon] = track.tpcNSigmaMu(); + nsigmas[kDptDptPion] = track.tpcNSigmaPi(); + nsigmas[kDptDptKaon] = track.tpcNSigmaKa(); + nsigmas[kDptDptProton] = track.tpcNSigmaPr(); + } + } + float min_nsigma = 999.0f; + MatchRecoGenSpecies sp_min_nsigma = kWrongSpecies; + for (int sp = 0; sp < kDptDptNoOfSpecies; ++sp) { + if (nsigmas[sp] < min_nsigma) { + min_nsigma = nsigmas[sp]; + sp_min_nsigma = MatchRecoGenSpecies(sp); + } + } + bool doublematch = false; + if (min_nsigma < 3.0) { + for (int sp = 0; (sp < kDptDptNoOfSpecies) and not doublematch; ++sp) { + if (sp != sp_min_nsigma) { + if (nsigmas[sp] < 3.0) { + doublematch = true; + } + } + } + if (doublematch) { + return kWrongSpecies; + } else { + return sp_min_nsigma; + } + } else { + return kWrongSpecies; + } + } + + template + inline MatchRecoGenSpecies IdentifyParticle(ParticleObject const& particle) + { + using namespace dptdptfilter; + + constexpr int pdgcodeEl = 11; + constexpr int pdgcodeMu = 13; + constexpr int pdgcodePi = 211; + constexpr int pdgcodeKa = 321; + constexpr int pdgcodePr = 2212; + + int pdgcode = abs(particle.pdgCode()); + + switch (pdgcode) { + case pdgcodeEl: + return kDptDptElectron; + break; + case pdgcodeMu: + return kDptDptMuon; + break; + case pdgcodePi: + return kDptDptPion; + break; + case pdgcodeKa: + return kDptDptKaon; + break; + case pdgcodePr: + return kDptDptProton; + break; + + default: + if (traceOutOfSpeciesParticles) { + LOGF(info, "Wrong particle passed selection cuts. PDG code: %d", pdgcode); + } + return kWrongSpecies; + break; + } + } + + template + MatchRecoGenSpecies trackIdentification(TrackObject const& track); + template + bool selectTrack(TrackObject const& track, int64_t colix); + template + void filterTracks(TrackListObject const& ftracks, int colix); + template + void filterParticles(ParticleListObject const& particles, MCCollisionObject const& mccollision, CollisionIndex colix); + + void init(InitContext const&) + { + using namespace dptdptfilter; + + LOGF(info, "DptDptFilterTask::init()"); + + /* update with the configurable values */ + /* the binning */ + ptbins = cfgBinning->mPTbins; + ptlow = cfgBinning->mPTmin; + ptup = cfgBinning->mPTmax; + etabins = cfgBinning->mEtabins; + etalow = cfgBinning->mEtamin; + etaup = cfgBinning->mEtamax; + zvtxbins = cfgBinning->mZVtxbins; + zvtxlow = cfgBinning->mZVtxmin; + zvtxup = cfgBinning->mZVtxmax; + /* the track types and combinations */ + tracktype = cfgTrackType.value; + /* the centrality/multiplicity estimation */ + fCentMultEstimator = getCentMultEstimator(cfgCentMultEstimator); + traceDCAOutliers = cfgTraceDCAOutliers; + traceOutOfSpeciesParticles = cfgTraceOutOfSpeciesParticles; + recoIdMethod = cfgRecoIdMethod; + if (cfgTrackSelection->mUseIt) { + useOwnTrackSelection = true; + if (cfgTrackSelection->mOnGen) { + useOwnParticleSelection = true; + particleMaxDCAxy = cfgTrackSelection->mDCAxy; + particleMaxDCAZ = cfgTrackSelection->mDCAz; + } + ownTrackSelection.SetMinNClustersTPC(cfgTrackSelection->mTPCclusters); + ownTrackSelection.SetMinNCrossedRowsTPC(cfgTrackSelection->mTPCxRows); + ownTrackSelection.SetMinNCrossedRowsOverFindableClustersTPC(cfgTrackSelection->mTPCXRoFClusters); + ownTrackSelection.SetMaxDcaXYPtDep(std::function{}); + ownTrackSelection.SetMaxDcaXY(cfgTrackSelection->mDCAxy); + ownTrackSelection.SetMaxDcaZ(cfgTrackSelection->mDCAz); + o2::aod::track::TrackTypeEnum ttype; + switch (tracktype) { + case 1: + ttype = o2::aod::track::Run2Track; + break; + case 3: + ttype = o2::aod::track::Track; + break; + default: + ttype = o2::aod::track::Track; + break; + } + ownTrackSelection.SetTrackType(ttype); + } else { + useOwnTrackSelection = false; + } + traceCollId0 = cfgTraceCollId0; + + /* if the system type is not known at this time, we have to put the initalization somewhere else */ + fSystem = getSystemType(cfgSystem); + fDataType = getDataType(cfgDataType); + fPDG = TDatabasePDG::Instance(); + + /* create the output list which will own the task histograms */ + TList* fOutputList = new TList(); + fOutputList->SetOwner(true); + fOutput.setObject(fOutputList); + + /* incorporate configuration parameters to the output */ + fOutputList->Add(new TParameter("TrackType", cfgTrackType, 'f')); + fOutputList->Add(new TParameter("TrackOneCharge", trackonecharge, 'f')); + fOutputList->Add(new TParameter("TrackTwoCharge", tracktwocharge, 'f')); + + if ((fDataType == kData) or (fDataType == kMC)) { + /* create the reconstructed data histograms */ + if (fSystem > kPbp) { + fhCentMultB = new TH1F("CentralityB", "Centrality before cut; centrality (%)", 100, 0, 100); + fhCentMultA = new TH1F("CentralityA", "Centrality; centrality (%)", 100, 0, 100); + } else { + /* for pp, pPb and Pbp systems use multiplicity instead */ + fhCentMultB = new TH1F("MultiplicityB", "Multiplicity before cut; multiplicity (%)", 100, 0, 100); + fhCentMultA = new TH1F("MultiplicityA", "Multiplicity; multiplicity (%)", 100, 0, 100); + } + + fhVertexZB = new TH1F("VertexZB", "Vertex Z; z_{vtx}", 60, -15, 15); + fhVertexZA = new TH1F("VertexZA", "Vertex Z; z_{vtx}", zvtxbins, zvtxlow, zvtxup); + + fhPB = new TH1F("fHistPB", "p distribution for reconstructed before;p (GeV/c);dN/dp (c/GeV)", 100, 0.0, 15.0); + fhPtB = new TH1F("fHistPtB", "p_{T} distribution for reconstructed before;p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); + fhPtPosB = new TH1F("fHistPtPosB", "P_{T} distribution for reconstructed (#plus) before;P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); + fhPtNegB = new TH1F("fHistPtNegB", "P_{T} distribution for reconstructed (#minus) before;P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); + fhEtaB = new TH1F("fHistEtaB", "#eta distribution for reconstructed before;#eta;counts", 40, -2.0, 2.0); + fhEtaA = new TH1F("fHistEtaA", "#eta distribution for reconstructed;#eta;counts", etabins, etalow, etaup); + fhPhiB = new TH1F("fHistPhiB", "#phi distribution for reconstructed before;#phi;counts", 360, 0.0, constants::math::TwoPI); + fhPhiA = new TH1F("fHistPhiA", "#phi distribution for reconstructed;#phi;counts", 360, 0.0, constants::math::TwoPI); + fhDCAxyB = new TH1F("DCAxyB", "DCA_{xy} distribution for reconstructed before;DCA_{xy} (cm);counts", 1000, -4.0, 4.0); + fhDCAxyA = new TH1F("DCAxyA", "DCA_{xy} distribution for reconstructed;DCA_{xy} (cm);counts", 1000, -4., 4.0); + fhFineDCAxyA = new TH1F("FineDCAxyA", "DCA_{xy} distribution for reconstructed;DCA_{xy} (cm);counts", 4000, -1.0, 1.0); + fhDCAzB = new TH1F("DCAzB", "DCA_{z} distribution for reconstructed before;DCA_{z} (cm);counts", 1000, -4.0, 4.0); + fhDCAzA = new TH1F("DCAzA", "DCA_{z} distribution for reconstructed;DCA_{z} (cm);counts", 1000, -4.0, 4.0); + fhFineDCAzA = new TH1F("FineDCAzA", "DCA_{z} distribution for reconstructed;DCA_{z} (cm);counts", 4000, -1.0, 1.0); + + for (int sp = 0; sp < kDptDptNoOfSpecies; ++sp) { + fhPA[sp] = new TH1F(TString::Format("fHistPA_%s", speciesName[sp]).Data(), + TString::Format("p distribution for reconstructed %s;p (GeV/c);dN/dp (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + fhPtA[sp] = new TH1F(TString::Format("fHistPtA_%s", speciesName[sp]), + TString::Format("p_{T} distribution for reconstructed %s;p_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + fhPtPosA[sp] = new TH1F(TString::Format("fHistPtPosA_%s", speciesName[sp]), + TString::Format("P_{T} distribution for reconstructed %s^{#plus};P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + fhPtNegA[sp] = new TH1F(TString::Format("fHistPtNegA_%s", speciesName[sp]), + TString::Format("P_{T} distribution for reconstructed %s^{#minus};P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + if (recoIdMethod == 0) { + /* no track identification, only charged */ + break; + } + } + + /* add the hstograms to the output list */ + fOutputList->Add(fhCentMultB); + fOutputList->Add(fhCentMultA); + fOutputList->Add(fhVertexZB); + fOutputList->Add(fhVertexZA); + fOutputList->Add(fhPB); + fOutputList->Add(fhPtB); + fOutputList->Add(fhPtPosB); + fOutputList->Add(fhPtNegB); + fOutputList->Add(fhEtaB); + fOutputList->Add(fhEtaA); + fOutputList->Add(fhPhiB); + fOutputList->Add(fhPhiA); + fOutputList->Add(fhDCAxyB); + fOutputList->Add(fhDCAxyA); + fOutputList->Add(fhFineDCAxyA); + fOutputList->Add(fhDCAzB); + fOutputList->Add(fhDCAzA); + fOutputList->Add(fhFineDCAzA); + + for (int sp = 0; sp < kDptDptNoOfSpecies; ++sp) { + fOutputList->Add(fhPA[sp]); + fOutputList->Add(fhPtA[sp]); + fOutputList->Add(fhPtPosA[sp]); + fOutputList->Add(fhPtNegA[sp]); + if (recoIdMethod == 0) { + /* no track identification, only charged */ + break; + } + } + } + + if (fDataType != kData) { + /* create the true data histograms */ + if (fSystem > kPbp) { + fhTrueCentMultB = new TH1F("TrueCentralityB", "Centrality before (truth); centrality (%)", 100, 0, 100); + fhTrueCentMultA = new TH1F("TrueCentralityA", "Centrality (truth); centrality (%)", 100, 0, 100); + } else { + /* for pp, pPb and Pbp systems use multiplicity instead */ + fhTrueCentMultB = new TH1F("TrueMultiplicityB", "Multiplicity before (truth); multiplicity (%)", 100, 0, 100); + fhTrueCentMultA = new TH1F("TrueMultiplicityA", "Multiplicity (truth); multiplicity (%)", 100, 0, 100); + } + + fhTrueVertexZB = new TH1F("TrueVertexZB", "Vertex Z before (truth); z_{vtx}", 60, -15, 15); + fhTrueVertexZA = new TH1F("TrueVertexZA", "Vertex Z (truth); z_{vtx}", zvtxbins, zvtxlow, zvtxup); + + fhTruePB = new TH1F("fTrueHistPB", "p distribution before (truth);p (GeV/c);dN/dp (c/GeV)", 100, 0.0, 15.0); + fhTruePtB = new TH1F("fTrueHistPtB", "p_{T} distribution before (truth);p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); + fhTruePtPosB = new TH1F("fTrueHistPtPosB", "P_{T} distribution (#plus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); + fhTruePtNegB = new TH1F("fTrueHistPtNegB", "P_{T} distribution (#minus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); + fhTrueEtaB = new TH1F("fTrueHistEtaB", "#eta distribution before (truth);#eta;counts", 40, -2.0, 2.0); + fhTrueEtaA = new TH1F("fTrueHistEtaA", "#eta distribution (truth);#eta;counts", etabins, etalow, etaup); + fhTruePhiB = new TH1F("fTrueHistPhiB", "#phi distribution before (truth);#phi;counts", 360, 0.0, constants::math::TwoPI); + fhTruePhiA = new TH1F("fTrueHistPhiA", "#phi distribution (truth);#phi;counts", 360, 0.0, constants::math::TwoPI); + fhTrueDCAxyB = new TH1F("TrueDCAxyB", "DCA_{xy} distribution for generated before;DCA_{xy} (cm);counts", 1000, -4.0, 4.0); + if (traceDCAOutliers.mDoIt) { + fhTrueDCAxyBid = new TH1F("PDGCodeDCAxyB", + TString::Format("PDG code within %.2f<|DCA_{#it{xy}}|<%.2f; PDG code", traceDCAOutliers.mLowValue, traceDCAOutliers.mUpValue).Data(), + 100, 0.5, 100.5); + } + fhTrueDCAxyA = new TH1F("TrueDCAxyA", "DCA_{xy} distribution for generated;DCA_{xy};counts (cm)", 1000, -4., 4.0); + fhTrueDCAzB = new TH1F("TrueDCAzB", "DCA_{z} distribution for generated before;DCA_{z} (cm);counts", 1000, -4.0, 4.0); + fhTrueDCAzA = new TH1F("TrueDCAzA", "DCA_{z} distribution for generated;DCA_{z} (cm);counts", 1000, -4.0, 4.0); + + for (int sp = 0; sp < kDptDptNoOfSpecies; ++sp) { + fhTruePA[sp] = new TH1F(TString::Format("fTrueHistPA_%s", speciesName[sp]).Data(), + TString::Format("p distribution %s (truth);p (GeV/c);dN/dp (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + fhTruePtA[sp] = new TH1F(TString::Format("fTrueHistPtA_%s", speciesName[sp]), + TString::Format("p_{T} distribution %s (truth);p_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + fhTruePtPosA[sp] = new TH1F(TString::Format("fTrueHistPtPosA_%s", speciesName[sp]), + TString::Format("P_{T} distribution %s^{#plus} (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + fhTruePtNegA[sp] = new TH1F(TString::Format("fTrueHistPtNegA_%s", speciesName[sp]), + TString::Format("P_{T} distribution %s^{#minus} (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), + ptbins, ptlow, ptup); + if (recoIdMethod == 0) { + /* no track identification, only charged */ + break; + } + } + + /* add the hstograms to the output list */ + fOutputList->Add(fhTrueCentMultB); + fOutputList->Add(fhTrueCentMultA); + fOutputList->Add(fhTrueVertexZB); + fOutputList->Add(fhTrueVertexZA); + fOutputList->Add(fhTruePB); + fOutputList->Add(fhTruePtB); + fOutputList->Add(fhTruePtPosB); + fOutputList->Add(fhTruePtNegB); + fOutputList->Add(fhTrueEtaB); + fOutputList->Add(fhTrueEtaA); + fOutputList->Add(fhTruePhiB); + fOutputList->Add(fhTruePhiA); + fOutputList->Add(fhTrueDCAxyB); + if (traceDCAOutliers.mDoIt) { + fOutputList->Add(fhTrueDCAxyBid); + } + fOutputList->Add(fhTrueDCAxyA); + fOutputList->Add(fhTrueDCAzB); + fOutputList->Add(fhTrueDCAzA); + + for (int sp = 0; sp < kDptDptNoOfSpecies; ++sp) { + fOutputList->Add(fhTruePA[sp]); + fOutputList->Add(fhTruePtA[sp]); + fOutputList->Add(fhTruePtPosA[sp]); + fOutputList->Add(fhTruePtNegA[sp]); + if (recoIdMethod == 0) { + /* no track identification, only charged */ + break; + } + } + } + } + + template + void processReconstructed(CollisionObject const& collision, TracksObject const& ftracks, float centormult); + + void processWithCent(aod::CollisionEvSelCent const& collision, DptDptFullTracks const& ftracks); + PROCESS_SWITCH(DptDptFilter, processWithCent, "Process reco with centrality", false); + + void processWithoutCent(aod::CollisionEvSel const& collision, DptDptFullTracks const& ftracks); + PROCESS_SWITCH(DptDptFilter, processWithoutCent, "Process reco without centrality", false); + + void processWithCentPID(aod::CollisionEvSelCent const& collision, DptDptFullTracksPID const& ftracks); + PROCESS_SWITCH(DptDptFilter, processWithCentPID, "Process PID reco with centrality", false); + + void processWithoutCentPID(aod::CollisionEvSel const& collision, DptDptFullTracksPID const& ftracks); + PROCESS_SWITCH(DptDptFilter, processWithoutCentPID, "Process PID reco without centrality", false); + + void processWithCentDetectorLevel(aod::CollisionEvSelCent const& collision, DptDptFullTracksDetLevel const& ftracks, aod::McParticles const&); + PROCESS_SWITCH(DptDptFilter, processWithCentDetectorLevel, "Process MC detector level with centrality", false); + + void processWithoutCentDetectorLevel(aod::CollisionEvSel const& collision, DptDptFullTracksDetLevel const& ftracks, aod::McParticles const&); + PROCESS_SWITCH(DptDptFilter, processWithoutCentDetectorLevel, "Process MC detector level without centrality", false); + + void processWithCentPIDDetectorLevel(aod::CollisionEvSelCent const& collision, DptDptFullTracksPIDDetLevel const& ftracks, aod::McParticles const&); + PROCESS_SWITCH(DptDptFilter, processWithCentPIDDetectorLevel, "Process PID MC detector level with centrality", false); + + void processWithoutCentPIDDetectorLevel(aod::CollisionEvSel const& collision, DptDptFullTracksPIDDetLevel const& ftracks, aod::McParticles const&); + PROCESS_SWITCH(DptDptFilter, processWithoutCentPIDDetectorLevel, "Process PID MC detector level without centrality", false); + + template + void processGenerated(CollisionObject const& mccollision, aod::McParticles const& mcparticles, float centormult); + + void processWithCentGeneratorLevel(aod::McCollision const& mccollision, + soa::Join const& collisions, + aod::McParticles const& mcparticles); + PROCESS_SWITCH(DptDptFilter, processWithCentGeneratorLevel, "Process generated with centrality", false); + + void processWithoutCentGeneratorLevel(aod::McCollision const& mccollision, + aod::McParticles const& mcparticles); + PROCESS_SWITCH(DptDptFilter, processWithoutCentGeneratorLevel, "Process generated without centrality", false); +}; + +template +MatchRecoGenSpecies DptDptFilter::trackIdentification(TrackObject const& track) +{ + using namespace dptdptfilter; + + MatchRecoGenSpecies sp = kWrongSpecies; + if (recoIdMethod == 0) { + sp = kDptDptCharged; + } else if (recoIdMethod == 1) { + if constexpr (framework::has_type_v) { + sp = IdentifyTrack(track); + } else { + LOGF(fatal, "Track identification required but PID information not present"); + } + } else if (recoIdMethod == 2) { + if constexpr (framework::has_type_v) { + sp = IdentifyParticle(track.mcParticle()); + } else { + LOGF(fatal, "Track identification required from MC particle but MC information not present"); + } + } + return sp; +} + +template +bool DptDptFilter::selectTrack(TrackObject const& track, int64_t colix) +{ + using namespace dptdptfilter; + + /* before track selection */ + fillTrackHistosBeforeSelection(track); + + /* track selection */ + /* tricky because the boolean columns issue */ + bool asone, astwo; + AcceptTrack(track, asone, astwo); + if (asone or astwo) { + /* the track has been accepted */ + fillTrackHistosAfterSelection(track, kDptDptCharged); + /* let's identify it */ + MatchRecoGenSpecies sp = trackIdentification(track); + if (sp != kWrongSpecies) { + /* fill the species histograms */ + fillTrackHistosAfterSelection(track, sp); + } + scannedtracks(colix, asone, astwo, track.pt(), track.eta(), track.phi()); + return true; + } + return false; +} + +template +void DptDptFilter::filterParticles(ParticleListObject const& particles, MCCollisionObject const& mccollision, CollisionIndex colix) +{ + using namespace dptdptfilter; + + int acceptedparticles = 0; + + for (auto& particle : particles) { + float charge = 0.0; + TParticlePDG* pdgparticle = fPDG->GetParticle(particle.pdgCode()); + if (pdgparticle != nullptr) { + charge = (pdgparticle->Charge() / 3 >= 1) ? 1.0 : ((pdgparticle->Charge() / 3 <= -1) ? -1.0 : 0.0); + } + + bool asone = false; + bool astwo = false; + if (charge != 0) { + /* before particle selection */ + fillParticleHistosBeforeSelection(particle, mccollision, charge); + + /* track selection */ + /* tricky because the boolean columns issue */ + AcceptParticle(particle, mccollision, asone, astwo); + if (asone or astwo) { + /* the track has been accepted */ + /* fill the charged particle histograms */ + fillParticleHistosAfterSelection(particle, mccollision, charge, kDptDptCharged); + /* let's identify the particle */ + MatchRecoGenSpecies sp = IdentifyParticle(particle); + if (sp != kWrongSpecies) { + fillParticleHistosAfterSelection(particle, mccollision, charge, sp); + } + acceptedparticles++; + scannedtruetracks(colix, asone, astwo, particle.pt(), particle.eta(), particle.phi()); + } + } else { + if ((particle.mcCollisionId() == 0) and traceCollId0) { + LOGF(info, "Particle %d with fractional charge or equal to zero", particle.globalIndex()); + } + } + } + LOGF(DPTDPTFILTERLOGCOLLISIONS, "Accepted %d generated particles", acceptedparticles); +} + +template +void DptDptFilter::filterTracks(TrackListObject const& ftracks, int colix) +{ + LOGF(fatal, "Track filtering not implemented for the passed track table"); +} + +template <> +void DptDptFilter::filterTracks(DptDptFullTracks const& ftracks, int colix) +{ + using namespace dptdptfilter; + + int acceptedtracks = 0; + + for (auto& track : ftracks) { + /* before track selection */ + if (selectTrack(track, colix)) { + acceptedtracks++; + } + } + LOGF(DPTDPTFILTERLOGCOLLISIONS, "Accepted %d reconstructed tracks", acceptedtracks); +} + +template <> +void DptDptFilter::filterTracks(DptDptFullTracksDetLevel const& ftracks, int colix) +{ + using namespace dptdptfilter; + + int acceptedtracks = 0; + + for (auto& track : ftracks) { + if (not(track.mcParticleId() < 0)) { + /* before track selection */ + if (selectTrack(track, colix)) { + acceptedtracks++; + } + } + } + LOGF(DPTDPTFILTERLOGCOLLISIONS, "Accepted %d reconstructed tracks", acceptedtracks); +} + +template <> +void DptDptFilter::filterTracks(DptDptFullTracksPID const& ftracks, int colix) +{ + using namespace dptdptfilter; + + int acceptedtracks = 0; + + for (auto& track : ftracks) { + /* before track selection */ + if (selectTrack(track, colix)) { + acceptedtracks++; + } + } + LOGF(DPTDPTFILTERLOGCOLLISIONS, "Accepted %d reconstructed tracks", acceptedtracks); +} + +template <> +void DptDptFilter::filterTracks(DptDptFullTracksPIDDetLevel const& ftracks, int colix) +{ + using namespace dptdptfilter; + + int acceptedtracks = 0; + + for (auto& track : ftracks) { + if (not(track.mcParticleId() < 0)) { + /* before track selection */ + if (selectTrack(track, colix)) { + acceptedtracks++; + } + } + } + LOGF(DPTDPTFILTERLOGCOLLISIONS, "Accepted %d reconstructed tracks", acceptedtracks); +} + +template +void DptDptFilter::processReconstructed(CollisionObject const& collision, TracksObject const& ftracks, float passedcent) +{ + using namespace dptdptfilter; + + LOGF(DPTDPTFILTERLOGCOLLISIONS, "DptDptFilterTask::processReconstructed(). New collision with %d tracks", ftracks.size()); + + fhCentMultB->Fill(passedcent); + fhVertexZB->Fill(collision.posZ()); + bool acceptedevent = false; + float centormult = passedcent; + if (IsEvtSelected(collision, centormult)) { + acceptedevent = true; + fhCentMultA->Fill(centormult); + fhVertexZA->Fill(collision.posZ()); + acceptedcollisions(collision.bcId(), collision.posZ(), acceptedevent, centormult); + + filterTracks(ftracks, acceptedcollisions.lastIndex()); + } else { + acceptedcollisions(collision.bcId(), collision.posZ(), acceptedevent, centormult); + for (auto& track : ftracks) { + scannedtracks(acceptedcollisions.lastIndex(), false, false, track.pt(), track.eta(), track.phi()); + } + } +} + +void DptDptFilter::processWithCent(aod::CollisionEvSelCent const& collision, DptDptFullTracks const& ftracks) +{ + processReconstructed(collision, ftracks, collision.centV0M()); +} + +void DptDptFilter::processWithoutCent(aod::CollisionEvSel const& collision, DptDptFullTracks const& ftracks) +{ + processReconstructed(collision, ftracks, 50.0); +} + +void DptDptFilter::processWithCentPID(aod::CollisionEvSelCent const& collision, DptDptFullTracksPID const& ftracks) +{ + processReconstructed(collision, ftracks, collision.centV0M()); +} + +void DptDptFilter::processWithoutCentPID(aod::CollisionEvSel const& collision, DptDptFullTracksPID const& ftracks) +{ + processReconstructed(collision, ftracks, 50.0); +} + +void DptDptFilter::processWithCentDetectorLevel(aod::CollisionEvSelCent const& collision, DptDptFullTracksDetLevel const& ftracks, aod::McParticles const&) +{ + processReconstructed(collision, ftracks, collision.centV0M()); +} + +void DptDptFilter::processWithoutCentDetectorLevel(aod::CollisionEvSel const& collision, DptDptFullTracksDetLevel const& ftracks, aod::McParticles const&) +{ + processReconstructed(collision, ftracks, 50.0); +} + +void DptDptFilter::processWithCentPIDDetectorLevel(aod::CollisionEvSelCent const& collision, DptDptFullTracksPIDDetLevel const& ftracks, aod::McParticles const&) +{ + processReconstructed(collision, ftracks, collision.centV0M()); +} + +void DptDptFilter::processWithoutCentPIDDetectorLevel(aod::CollisionEvSel const& collision, DptDptFullTracksPIDDetLevel const& ftracks, aod::McParticles const&) +{ + processReconstructed(collision, ftracks, 50.0); +} + +template +void DptDptFilter::processGenerated(CollisionObject const& mccollision, aod::McParticles const& mcparticles, float centormult) +{ + using namespace dptdptfilter; + + fhTrueCentMultB->Fill(centormult); + fhTrueVertexZB->Fill(mccollision.posZ()); + bool acceptedevent = false; + if (IsEvtSelected(mccollision, centormult)) { + acceptedevent = true; + fhTrueCentMultA->Fill(centormult); + fhTrueVertexZA->Fill(mccollision.posZ()); + acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), acceptedevent, centormult); + + filterParticles(mcparticles, mccollision, acceptedtrueevents.lastIndex()); + } else { + acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), acceptedevent, centormult); + for (auto& particle : mcparticles) { + scannedtruetracks(acceptedtrueevents.lastIndex(), false, false, particle.pt(), particle.eta(), particle.phi()); + } + } +} + +void DptDptFilter::processWithCentGeneratorLevel(aod::McCollision const& mccollision, + soa::Join const& collisions, + aod::McParticles const& mcparticles) +{ + using namespace dptdptfilter; + + LOGF(DPTDPTFILTERLOGCOLLISIONS, "DptDptFilterTask::processWithCentGeneratorLevel(). New generated collision %d reconstructed collisions and %d particles", collisions.size(), mcparticles.size()); + + /* TODO: in here we have to decide what to do in the following cases + - On the fly production -> clearly we will need a different process + - reconstructed collisions without generated associated -> we need a different process or a different signature + - multiplicity/centrality classes extracted from the reconstructed collision but then + - generated collision without associated reconstructed collision: how to extract mutliplicity/centrality classes? + - generated collision with several associated reconstructed collisions: from which to extract multiplicity/centrality classes? + */ + if (collisions.size() > 1) { + LOGF(error, "DptDptFilterTask::processWithCentGeneratorLevel(). Generated collision with more than one reconstructed collisions. Processing only the first for centrality/multiplicity classes extraction"); + } + + for (auto& collision : collisions) { + processGenerated(mccollision, mcparticles, collision.centV0M()); + break; /* TODO: only processing the first reconstructed collision for centrality/multiplicity class estimation */ + } +} + +void DptDptFilter::processWithoutCentGeneratorLevel(aod::McCollision const& mccollision, + aod::McParticles const& mcparticles) +{ + using namespace dptdptfilter; + + LOGF(DPTDPTFILTERLOGCOLLISIONS, "DptDptFilterTask::processWithoutCentGeneratorLevel(). New generated collision with %d particles", mcparticles.size()); + + /* the task does not have access to either centrality nor multiplicity + classes information, so it has to live without it. + For the time being we assign a value of 50% */ + processGenerated(mccollision, mcparticles, 50.0); +} + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, SetDefaultProcesses{{{"processWithoutCent", true}, {"processWithoutCentMC", true}}})}; + return workflow; +} diff --git a/PWGCF/TableProducer/dptdptfilter.h b/PWGCF/TableProducer/dptdptfilter.h new file mode 100644 index 00000000000..96b83896736 --- /dev/null +++ b/PWGCF/TableProducer/dptdptfilter.h @@ -0,0 +1,433 @@ +// 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. +#ifndef O2_ANALYSIS_DPTDPTFILTER_H +#define O2_ANALYSIS_DPTDPTFILTER_H + +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/Core/MC.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "PWGCF/Core/AnalysisConfigurableCuts.h" +#include + +namespace o2 +{ +namespace aod +{ +using CollisionsEvSelCent = soa::Join; +using CollisionEvSelCent = soa::Join::iterator; +using CollisionsEvSel = soa::Join; +using CollisionEvSel = soa::Join::iterator; +using TrackData = soa::Join::iterator; +} // namespace aod +namespace analysis +{ +namespace dptdptfilter +{ +/// \enum SystemType +/// \brief The type of the system under analysis +enum SystemType { + kNoSystem = 0, ///< no system defined + kpp, ///< **p-p** system + kpPb, ///< **p-Pb** system + kPbp, ///< **Pb-p** system + kPbPb, ///< **Pb-Pb** system + kXeXe, ///< **Xe-Xe** system + knSystems ///< number of handled systems +}; + +/// \enum DataType +/// \brief Which kind of data is the task addressing +enum DataType { + kData = 0, ///< actual data, not generated + kMC, ///< Generator level and detector level + kFastMC, ///< Gererator level but stored dataset + kOnTheFly, ///< On the fly generator level data + knGenData ///< number of different generator data types +}; + +/// \enum CentMultEstimatorType +/// \brief The detector used to estimate centrality/multiplicity +enum CentMultEstimatorType { + kNOCM = 0, ///< do not use centrality/multiplicity estimator + kV0M, ///< V0M centrality/multiplicity estimator + kV0A, ///< V0A centrality/multiplicity estimator + kV0C, ///< V0C centrality/multiplicity estimator + kCL0, ///< CL0 centrality/multiplicity estimator + kCL1, ///< CL1 centrality/multiplicity estimator + knCentMultEstimators ///< number of centrality/mutiplicity estimator +}; + +//============================================================================================ +// The DptDptFilter configuration objects +//============================================================================================ +int ptbins = 18; +float ptlow = 0.2, ptup = 2.0; +int etabins = 16; +float etalow = -0.8, etaup = 0.8; +int zvtxbins = 40; +float zvtxlow = -10.0, zvtxup = 10.0; +int phibins = 72; +float philow = 0.0; +float phiup = constants::math::TwoPI; + +int tracktype = 1; +int trackonecharge = 1; +int tracktwocharge = -1; + +SystemType fSystem = kNoSystem; +DataType fDataType = kData; +CentMultEstimatorType fCentMultEstimator = kV0M; + +/* adaptations for the pp nightly checks */ +analysis::CheckRangeCfg traceDCAOutliers; +bool traceOutOfSpeciesParticles = false; +int recoIdMethod = 0; +bool useOwnTrackSelection = false; +TrackSelection ownTrackSelection = getGlobalTrackSelection(); +bool useOwnParticleSelection = false; +bool particleMaxDCAxy = 999.9; +bool particleMaxDCAZ = 999.9; +bool traceCollId0 = false; + +TDatabasePDG* fPDG = nullptr; + +inline SystemType getSystemType(std::string const& sysstr) +{ + /* we have to figure out how extract the system type */ + if (sysstr.empty() or (sysstr == "PbPb")) { + return kPbPb; + } else if (sysstr == "pp") { + return kpp; + } else if (sysstr == "pPb") { + return kpPb; + } else if (sysstr == "Pbp") { + return kPbp; + } else if (sysstr == "pPb") { + return kpPb; + } else if (sysstr == "XeXe") { + return kXeXe; + } else { + LOGF(fatal, "DptDptCorrelations::getSystemType(). Wrong system type: %d", sysstr.c_str()); + } + return kPbPb; +} + +/// \brief Type of data according to the configuration string +/// \param datastr The data type configuration string +/// \return Internal code for the passed kind of data string +inline DataType getDataType(std::string const& datastr) +{ + /* we have to figure out how extract the type of data*/ + if (datastr.empty() or (datastr == "data")) { + return kData; + } else if (datastr == "MC") { + return kMC; + } else if (datastr == "FastMC") { + return kFastMC; + } else if (datastr == "OnTheFlyMC") { + return kOnTheFly; + } else { + LOGF(fatal, "DptDptCorrelations::getDataType(). Wrong type of dat: %d", datastr.c_str()); + } + return kData; +} + +inline CentMultEstimatorType getCentMultEstimator(std::string const& datastr) +{ + if (datastr == "V0M") { + return kV0M; + } else if (datastr == "NOCM") { + return kNOCM; + } else { + LOGF(fatal, "Centrality/Multiplicity estimator %s not supported yet", datastr.c_str()); + } + return kNOCM; +} + +////////////////////////////////////////////////////////////////////////////////// +/// Trigger selection +////////////////////////////////////////////////////////////////////////////////// + +/// \brief Trigger selection for reconstructed and detector level collision tables +template +inline bool triggerSelectionReco(CollisionObject const& collision) +{ + bool trigsel = false; + if (fDataType != kData) { + trigsel = true; + } else if (collision.alias()[kINT7]) { + if (collision.sel7()) { + trigsel = true; + } + } + return trigsel; +} + +/// \brief Trigger selection by default: unknow subscribed collision table +template +inline bool triggerSelection(CollisionObject const& collision) +{ + LOGF(fatal, "Trigger selection not implemented for this kind of collisions"); + return false; +} + +/// \brief Trigger selection for reconstructed collision tables without centrality/multiplicity +template <> +inline bool triggerSelection(aod::CollisionEvSel const& collision) +{ + return triggerSelectionReco(collision); +} + +/// \brief Trigger selection for reconstructed collision tables with centrality/multiplicity +template <> +inline bool triggerSelection(aod::CollisionEvSelCent const& collision) +{ + return triggerSelectionReco(collision); +} + +/// \brief Trigger selection for detector level collision tables without centrality/multiplicity +template <> +inline bool triggerSelection::iterator>(soa::Join::iterator const& collision) +{ + return triggerSelectionReco(collision); +} + +/// \brief Trigger selection for detector level collision tables with centrality/multiplicity +template <> +inline bool triggerSelection::iterator>(soa::Join::iterator const& collision) +{ + return triggerSelectionReco(collision); +} + +/// \brief Trigger selection for generator level collison table +template <> +inline bool triggerSelection(aod::McCollision const& collision) +{ + return true; +} + +////////////////////////////////////////////////////////////////////////////////// +/// Centrality selection +////////////////////////////////////////////////////////////////////////////////// + +/// \brief Centrality selection when there is centrality/multiplicity information +template +inline bool centralitySelectionMult(CollisionObject collision, float& centmult) +{ + bool centmultsel = false; + switch (fCentMultEstimator) { + case kV0M: + if (collision.centV0M() < 100 and 0 < collision.centV0M()) { + centmult = collision.centV0M(); + centmultsel = true; + } + break; + default: + break; + } + return centmultsel; +} + +/// \brief Centrality selection when there is not centrality/multiplicity information +template +inline bool centralitySelectionNoMult(CollisionObject collision, float& centmult) +{ + bool centmultsel = false; + switch (fCentMultEstimator) { + case kNOCM: + centmult = 50.0; + centmultsel = true; + break; + default: + break; + } + return centmultsel; +} + +/// \brief Centrality selection by default: unknown subscribed collision table +template +inline bool centralitySelection(CollisionObject const& collision, float& centmult) +{ + LOGF(fatal, "Centrality selection not implemented for this kind of collisions"); + return false; +} + +/// \brief Centrality selection for reconstructed and detector level collision tables with centrality/multiplicity information +template <> +inline bool centralitySelection(aod::CollisionEvSelCent const& collision, float& centmult) +{ + return centralitySelectionMult(collision, centmult); +} + +/// \brief Centrality selection for reconstructed and detector level collision tables without centrality/multiplicity information +template <> +inline bool centralitySelection(aod::CollisionEvSel const& collision, float& centmult) +{ + return centralitySelectionNoMult(collision, centmult); +} + +/// \brief Centrality selection for detector level collision tables without centrality/multiplicity +template <> +inline bool centralitySelection::iterator>(soa::Join::iterator const& collision, float& centmult) +{ + return centralitySelectionNoMult(collision, centmult); +} + +/// \brief Centrality selection for detector level collision tables with centrality/multiplicity +template <> +inline bool centralitySelection::iterator>(soa::Join::iterator const& collision, float& centmult) +{ + return centralitySelectionMult(collision, centmult); +} + +/// \brief Centrality selection for generator level collision table +template <> +inline bool centralitySelection(aod::McCollision const& collision, float& centmult) +{ + if (centmult < 100 and 0 < centmult) { + return true; + } else { + return false; + } +} + +////////////////////////////////////////////////////////////////////////////////// +/// Event selection +////////////////////////////////////////////////////////////////////////////////// + +template +inline bool IsEvtSelected(CollisionObject const& collision, float& centormult) +{ + bool trigsel = triggerSelection(collision); + + bool zvtxsel = false; + /* TODO: vertex quality checks */ + if (zvtxlow < collision.posZ() and collision.posZ() < zvtxup) { + zvtxsel = true; + } + + bool centmultsel = centralitySelection(collision, centormult); + + return trigsel and zvtxsel and centmultsel; +} + +////////////////////////////////////////////////////////////////////////////////// +/// Track selection +////////////////////////////////////////////////////////////////////////////////// + +template +inline bool matchTrackType(TrackObject const& track) +{ + if (useOwnTrackSelection) { + return ownTrackSelection.IsSelected(track); + } else { + switch (tracktype) { + case 1: + if (track.isGlobalTrack() != 0 || track.isGlobalTrackSDD() != 0) { + return true; + } else { + return false; + } + break; + case 3: /* Run3 track */ + if (track.isGlobalTrack() != 0 || track.isGlobalTrackSDD() != 0) { + return true; + } else { + return false; + } + break; + default: + return false; + } + } +} + +template +inline void AcceptTrack(TrackObject const& track, bool& asone, bool& astwo) +{ + asone = false; + astwo = false; + + /* TODO: incorporate a mask in the scanned tracks table for the rejecting track reason */ + if (matchTrackType(track)) { + if (ptlow < track.pt() and track.pt() < ptup and etalow < track.eta() and track.eta() < etaup) { + if (((track.sign() > 0) and (trackonecharge > 0)) or ((track.sign() < 0) and (trackonecharge < 0))) { + asone = true; + } + if (((track.sign() > 0) and (tracktwocharge > 0)) or ((track.sign() < 0) and (tracktwocharge < 0))) { + astwo = true; + } + } + } +} + +template +inline void AcceptParticle(ParticleObject& particle, MCCollisionObject const& collision, bool& asone, bool& astwo) +{ + asone = false; + astwo = false; + + float charge = (fPDG->GetParticle(particle.pdgCode())->Charge() / 3 >= 1) ? 1.0 : ((fPDG->GetParticle(particle.pdgCode())->Charge() / 3 <= -1) ? -1.0 : 0.0); + + if (MC::isPhysicalPrimary(particle)) { + if ((particle.mcCollisionId() == 0) and traceCollId0) { + LOGF(info, "Particle %d passed isPhysicalPrimary", particle.globalIndex()); + } + if (useOwnParticleSelection) { + float dcaxy = TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + + (particle.vy() - collision.posY()) * (particle.vy() - collision.posY())); + float dcaz = TMath::Abs(particle.vz() - collision.posZ()); + if (not((dcaxy < particleMaxDCAxy) and (dcaz < particleMaxDCAZ))) { + if ((particle.mcCollisionId() == 0) and traceCollId0) { + auto currparticle = particle; + LOGF(info, "Rejecting particle with dcaxy: %.2f and dcaz: %.2f", dcaxy, dcaz); + LOGF(info, " assigned collision Id: %d, looping on collision Id: %d", currparticle.mcCollisionId(), collision.globalIndex()); + LOGF(info, " Collision x: %.5f, y: %.5f, z: %.5f", collision.posX(), collision.posY(), collision.posZ()); + LOGF(info, " Particle x: %.5f, y: %.5f, z: %.5f", particle.vx(), particle.vy(), particle.vz()); + LOGF(info, " index: %d, pdg code: %d", currparticle.globalIndex(), currparticle.pdgCode()); + while (currparticle.has_mother0()) { + LOGF(info, " mother0 index: %d, mother1 index: %d", currparticle.mother0Id(), currparticle.mother1Id()); + LOGF(info, " Tracking back mother0 index"); + auto newcurrparticle = currparticle.template mother0_as(); + LOGF(info, " assigned collision Id: %d, looping on collision Id: %d", newcurrparticle.mcCollisionId(), collision.globalIndex()); + LOGF(info, " index: %d, pdg code: %d", newcurrparticle.globalIndex(), newcurrparticle.pdgCode()); + LOGF(info, " Passed isPhysicalPrimary(): %s", MC::isPhysicalPrimary(newcurrparticle) ? "YES" : "NO"); + currparticle = newcurrparticle; + } + } + return; + } + } + if (ptlow < particle.pt() and particle.pt() < ptup and etalow < particle.eta() and particle.eta() < etaup) { + if (((charge > 0) and (trackonecharge > 0)) or ((charge < 0) and (trackonecharge < 0))) { + asone = true; + } + if (((charge > 0) and (tracktwocharge > 0)) or ((charge < 0) and (tracktwocharge < 0))) { + astwo = true; + } + } + } else { + if ((particle.mcCollisionId() == 0) and traceCollId0) { + LOGF(info, "Particle %d NOT passed isPhysicalPrimary", particle.globalIndex()); + } + } +} + +} // namespace dptdptfilter +} // namespace analysis +} // namespace o2 + +#endif // O2_ANALYSIS_DPTDPTFILTER_H diff --git a/PWGCF/Tasks/CMakeLists.txt b/PWGCF/Tasks/CMakeLists.txt index a3c4b13b4d1..7dff6b4ec4e 100644 --- a/PWGCF/Tasks/CMakeLists.txt +++ b/PWGCF/Tasks/CMakeLists.txt @@ -14,6 +14,11 @@ o2physics_add_dpl_workflow(dptdptcorrelations PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(dptdpt-filter-qa + SOURCES dptdptfilterqa.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(match-reco-gen SOURCES match-reco-gen.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore diff --git a/PWGCF/Tasks/dptdptcorrelations.cxx b/PWGCF/Tasks/dptdptcorrelations.cxx index 1aa544ee2a3..88c427858b3 100644 --- a/PWGCF/Tasks/dptdptcorrelations.cxx +++ b/PWGCF/Tasks/dptdptcorrelations.cxx @@ -17,7 +17,10 @@ #include "Common/Core/TrackSelection.h" #include "Common/Core/MC.h" #include "PWGCF/Core/AnalysisConfigurableCuts.h" +#include "PWGCF/DataModel/DptDptFiltered.h" +#include "PWGCF/TableProducer/dptdptfilter.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "Framework/runDataProcessing.h" #include #include #include @@ -39,79 +42,18 @@ using namespace o2::framework::expressions; #define DPTDPTLOGCOLLISIONS debug #define DPTDPTLOGTRACKS debug -#include "Framework/runDataProcessing.h" - namespace o2 { namespace aod { -/* we have to change from int to bool when bool columns work properly */ -namespace dptdptcorrelations -{ -DECLARE_SOA_COLUMN(EventAccepted, eventaccepted, uint8_t); //! If the collision/event has been accepted or not -DECLARE_SOA_COLUMN(EventCentMult, centmult, float); //! The centrality/multiplicity pecentile -} // namespace dptdptcorrelations -DECLARE_SOA_TABLE(AcceptedEvents, "AOD", "ACCEPTEDEVENTS", //! Accepted reconstructed collisions/events filtered table - o2::soa::Index<>, - collision::BCId, - collision::PosZ, - dptdptcorrelations::EventAccepted, - dptdptcorrelations::EventCentMult); -using AcceptedEvent = AcceptedEvents::iterator; -DECLARE_SOA_TABLE(AcceptedTrueEvents, "AOD", "ACCTRUEEVENTS", //! Accepted generated collisions/events filtered table - o2::soa::Index<>, - collision::BCId, - mccollision::PosZ, - dptdptcorrelations::EventAccepted, - dptdptcorrelations::EventCentMult); -using AcceptedTrueEvent = AcceptedTrueEvents::iterator; -namespace dptdptcorrelations -{ -DECLARE_SOA_INDEX_COLUMN(AcceptedEvent, event); //! Reconstructed collision/event -DECLARE_SOA_INDEX_COLUMN(AcceptedTrueEvent, mcevent); //! Generated collision/event -DECLARE_SOA_COLUMN(TrackacceptedAsOne, trackacceptedasone, uint8_t); //! Track accepted as type one -DECLARE_SOA_COLUMN(TrackacceptedAsTwo, trackacceptedastwo, uint8_t); //! Track accepted as type two -DECLARE_SOA_COLUMN(Pt, pt, float); //! The track transverse momentum -DECLARE_SOA_COLUMN(Eta, eta, float); //! The track pseudorapidity -DECLARE_SOA_COLUMN(Phi, phi, float); //! The track azimuthal angle -} // namespace dptdptcorrelations -DECLARE_SOA_TABLE(ScannedTracks, "AOD", "SCANNEDTRACKS", //! The reconstructed tracks filtered table - dptdptcorrelations::AcceptedEventId, - dptdptcorrelations::TrackacceptedAsOne, - dptdptcorrelations::TrackacceptedAsTwo, - dptdptcorrelations::Pt, - dptdptcorrelations::Eta, - dptdptcorrelations::Phi); -DECLARE_SOA_TABLE(ScannedTrueTracks, "AOD", "SCANTRUETRACKS", //! The generated particles filtered table - dptdptcorrelations::AcceptedTrueEventId, - dptdptcorrelations::TrackacceptedAsOne, - dptdptcorrelations::TrackacceptedAsTwo, - dptdptcorrelations::Pt, - dptdptcorrelations::Eta, - dptdptcorrelations::Phi); - -using CollisionsEvSelCent = soa::Join; -using CollisionEvSelCent = soa::Join::iterator; -using CollisionsEvSel = soa::Join; -using CollisionEvSel = soa::Join::iterator; -using TrackData = soa::Join::iterator; using FilteredTracks = soa::Filtered>; using FilteredTrackData = Partition::filtered_iterator; } // namespace aod } // namespace o2 -namespace dptdptcorrelations +namespace correlationstask { -/* all this is made configurable */ -int ptbins = 18; -float ptlow = 0.2, ptup = 2.0; -int etabins = 16; -float etalow = -0.8, etaup = 0.8; -int zvtxbins = 40; -float zvtxlow = -10.0, zvtxup = 10.0; -int phibins = 72; -float philow = 0.0; -float phiup = M_PI * 2; +using namespace o2::analysis::dptdptfilter; float phibinshift = 0.5; float etabinwidth = (etaup - etalow) / float(etabins); float phibinwidth = (phiup - philow) / float(phibins); @@ -119,134 +61,13 @@ int deltaetabins = etabins * 2 - 1; float deltaetalow = etalow - etaup, deltaetaup = etaup - etalow; float deltaetabinwidth = (deltaetaup - deltaetalow) / float(deltaetabins); int deltaphibins = phibins; -float deltaphibinwidth = M_PI * 2 / deltaphibins; +float deltaphibinwidth = constants::math::TwoPI / deltaphibins; float deltaphilow = 0.0 - deltaphibinwidth / 2.0; -float deltaphiup = M_PI * 2 - deltaphibinwidth / 2.0; +float deltaphiup = constants::math::TwoPI - deltaphibinwidth / 2.0; -int tracktype = 1; -int trackonecharge = 1; -int tracktwocharge = -1; bool processpairs = false; std::string fTaskConfigurationString = "PendingToConfigure"; -/// \enum SystemType -/// \brief The type of the system under analysis -enum SystemType { - kNoSystem = 0, ///< no system defined - kpp, ///< **p-p** system - kpPb, ///< **p-Pb** system - kPbp, ///< **Pb-p** system - kPbPb, ///< **Pb-Pb** system - kXeXe, ///< **Xe-Xe** system - knSystems ///< number of handled systems -}; - -/// \enum GeneratorType -/// \brief Which kid of generator data is the task addressing -enum GenType { - kData = 0, ///< actual data, not generated - kMC, ///< Generator level and detector level - kFastMC, ///< Gererator level but stored dataset - kOnTheFly, ///< On the fly generator level data - knGenData ///< number of different generator data types -}; - -/// \enum CentMultEstimatorType -/// \brief The detector used to estimate centrality/multiplicity -enum CentMultEstimatorType { - kNOCM = 0, ///< do not use centrality/multiplicity estimator - kV0M, ///< V0M centrality/multiplicity estimator - kV0A, ///< V0A centrality/multiplicity estimator - kV0C, ///< V0C centrality/multiplicity estimator - kCL0, ///< CL0 centrality/multiplicity estimator - kCL1, ///< CL1 centrality/multiplicity estimator - knCentMultEstimators ///< number of centrality/mutiplicity estimator -}; - -namespace filteranalysistask -{ -//============================================================================================ -// The DptDptCorrelationsFilterAnalysisTask output objects -//============================================================================================ -SystemType fSystem = kNoSystem; -GenType fDataType = kData; -CentMultEstimatorType fCentMultEstimator = kV0M; -TDatabasePDG* fPDG = nullptr; -TH1F* fhCentMultB = nullptr; -TH1F* fhCentMultA = nullptr; -TH1F* fhVertexZB = nullptr; -TH1F* fhVertexZA = nullptr; -TH1F* fhPtB = nullptr; -TH1F* fhPtA = nullptr; -TH1F* fhPtPosB = nullptr; -TH1F* fhPtPosA = nullptr; -TH1F* fhPtNegB = nullptr; -TH1F* fhPtNegA = nullptr; - -TH1F* fhEtaB = nullptr; -TH1F* fhEtaA = nullptr; - -TH1F* fhPhiB = nullptr; -TH1F* fhPhiA = nullptr; - -TH2F* fhEtaVsPhiB = nullptr; -TH2F* fhEtaVsPhiA = nullptr; - -TH2F* fhPtVsEtaB = nullptr; -TH2F* fhPtVsEtaA = nullptr; - -TH1F* fhDCAxyB = nullptr; -TH1F* fhDCAxyA = nullptr; -TH1F* fhFineDCAxyA = nullptr; -TH1F* fhDCAzB = nullptr; -TH1F* fhDCAzA = nullptr; -TH1F* fhFineDCAzA = nullptr; - -TH1F* fhTrueCentMultB = nullptr; -TH1F* fhTrueCentMultA = nullptr; -TH1F* fhTrueVertexZB = nullptr; -TH1F* fhTrueVertexZA = nullptr; -TH1F* fhTruePtB = nullptr; -TH1F* fhTruePtA = nullptr; -TH1F* fhTruePtPosB = nullptr; -TH1F* fhTruePtPosA = nullptr; -TH1F* fhTruePtNegB = nullptr; -TH1F* fhTruePtNegA = nullptr; - -TH1F* fhTrueEtaB = nullptr; -TH1F* fhTrueEtaA = nullptr; - -TH1F* fhTruePhiB = nullptr; -TH1F* fhTruePhiA = nullptr; - -TH2F* fhTrueEtaVsPhiB = nullptr; -TH2F* fhTrueEtaVsPhiA = nullptr; - -TH2F* fhTruePtVsEtaB = nullptr; -TH2F* fhTruePtVsEtaA = nullptr; - -TH1F* fhTrueDCAxyB = nullptr; -TH1F* fhTrueDCAxyA = nullptr; -TH1F* fhTrueDCAzB = nullptr; -TH1F* fhTrueDCAzA = nullptr; -} // namespace filteranalysistask - -namespace filteranalysistaskqa -{ -TH1F* fhTracksOne = nullptr; -TH1F* fhTracksTwo = nullptr; -TH1F* fhTracksOneAndTwo = nullptr; -TH1F* fhTracksNone = nullptr; -TH1F* fhTracksOneUnsel = nullptr; -TH1F* fhTracksTwoUnsel = nullptr; -TH1F* fhTracksOneAndTwoUnsel = nullptr; -TH1F* fhTracksNoneUnsel = nullptr; -TH1F* fhSelectedEvents = nullptr; -} // namespace filteranalysistaskqa - -namespace correlationstask -{ - /// \enum TrackPairs /// \brief The track combinations hadled by the class enum TrackPairs { @@ -258,758 +79,6 @@ enum TrackPairs { }; } // namespace correlationstask -/// \brief System type according to configuration string -/// \param sysstr The system configuration string -/// \return The internal code for the passed system string -SystemType getSystemType(std::string const& sysstr) -{ - /* we have to figure out how extract the system type */ - if (sysstr.empty() or (sysstr == "PbPb")) { - return kPbPb; - } else if (sysstr == "pp") { - return kpp; - } else if (sysstr == "pPb") { - return kpPb; - } else if (sysstr == "Pbp") { - return kPbp; - } else if (sysstr == "pPb") { - return kpPb; - } else if (sysstr == "XeXe") { - return kXeXe; - } else { - LOGF(fatal, "DptDptCorrelations::getSystemType(). Wrong system type: %d", sysstr.c_str()); - } - return kPbPb; -} - -/// \brief Type of data according to the configuration string -/// \param datastr The data type configuration string -/// \return Internal code for the passed kind of data string -GenType getGenType(std::string const& datastr) -{ - /* we have to figure out how extract the type of data*/ - if (datastr.empty() or (datastr == "data")) { - return kData; - } else if (datastr == "MC") { - return kMC; - } else if (datastr == "FastMC") { - return kFastMC; - } else if (datastr == "OnTheFlyMC") { - return kOnTheFly; - } else { - LOGF(fatal, "DptDptCorrelations::getGenType(). Wrong type of dat: %d", datastr.c_str()); - } - return kData; -} - -template -bool IsEvtSelected(CollisionObject const& collision, float& centormult) -{ - using namespace filteranalysistask; - using namespace dptdptcorrelations; - - bool trigsel = false; - if (fDataType != kData) { - trigsel = true; - } else if (collision.alias()[kINT7]) { - if (collision.sel7()) { - trigsel = true; - } - } - - bool zvtxsel = false; - /* TODO: vertex quality checks */ - if (zvtxlow < collision.posZ() and collision.posZ() < zvtxup) { - zvtxsel = true; - } - - bool centmultsel = false; - switch (fCentMultEstimator) { - case kV0M: - if (collision.centV0M() < 100 and 0 < collision.centV0M()) { - centormult = collision.centV0M(); - centmultsel = true; - } - break; - default: - break; - } - return trigsel and zvtxsel and centmultsel; -} - -template -bool IsEvtSelectedNoCentMult(CollisionObject const& collision, float& centormult) -{ - using namespace filteranalysistask; - using namespace dptdptcorrelations; - - bool trigsel = false; - if (fDataType != kData) { - trigsel = true; - } else if (collision.alias()[kINT7]) { - if (collision.sel7() or collision.sel8()) { - trigsel = true; - } - } - - bool zvtxsel = false; - /* TODO: vertex quality checks */ - if (zvtxlow < collision.posZ() and collision.posZ() < zvtxup) { - zvtxsel = true; - } - - bool centmultsel = false; - switch (fCentMultEstimator) { - case kNOCM: - centormult = 50.0; - centmultsel = true; - break; - default: - break; - } - return trigsel and zvtxsel and centmultsel; -} - -template -bool IsTrueEvtSelected(CollisionObject const& collision, float centormult) -{ - using namespace filteranalysistask; - using namespace dptdptcorrelations; - - bool zvtxsel = false; - /* TODO: vertex quality checks */ - if (zvtxlow < collision.posZ() and collision.posZ() < zvtxup) { - zvtxsel = true; - } - - bool centmultsel = false; - if (centormult < 100 and 0 < centormult) { - centmultsel = true; - } - - return zvtxsel and centmultsel; -} - -template -bool matchTrackType(TrackObject const& track) -{ - using namespace filteranalysistask; - - switch (tracktype) { - case 1: - if (track.isGlobalTrack() != 0 || track.isGlobalTrackSDD() != 0) { - return true; - } else { - return false; - } - break; - case 3: /* Run3 track */ - if (track.isGlobalTrack() != 0 || track.isGlobalTrackSDD() != 0) { - return true; - } else { - return false; - } - break; - default: - return false; - } -} - -template -inline void AcceptTrack(TrackObject const& track, bool& asone, bool& astwo) -{ - using namespace filteranalysistask; - - asone = false; - astwo = false; - - /* TODO: incorporate a mask in the scanned tracks table for the rejecting track reason */ - if (matchTrackType(track)) { - if (ptlow < track.pt() and track.pt() < ptup and etalow < track.eta() and track.eta() < etaup) { - if (((track.sign() > 0) and (trackonecharge > 0)) or ((track.sign() < 0) and (trackonecharge < 0))) { - asone = true; - } - if (((track.sign() > 0) and (tracktwocharge > 0)) or ((track.sign() < 0) and (tracktwocharge < 0))) { - astwo = true; - } - } - } -} - -template -inline void AcceptParticle(ParticleObject& particle, ParticleListObject& particles, bool& asone, bool& astwo) -{ - using namespace filteranalysistask; - - asone = false; - astwo = false; - - float charge = (fPDG->GetParticle(particle.pdgCode())->Charge() / 3 >= 1) ? 1.0 : ((fPDG->GetParticle(particle.pdgCode())->Charge() / 3 <= -1) ? -1.0 : 0.0); - - /* TODO: matchTrackType will not work. We need at least is physical primary */ - if (MC::isPhysicalPrimary(particle)) { - if (ptlow < particle.pt() and particle.pt() < ptup and etalow < particle.eta() and particle.eta() < etaup) { - if (((charge > 0) and (trackonecharge > 0)) or ((charge < 0) and (trackonecharge < 0))) { - asone = true; - } - if (((charge > 0) and (tracktwocharge > 0)) or ((charge < 0) and (tracktwocharge < 0))) { - astwo = true; - } - } - } -} - -template -void fillTrackHistosBeforeSelection(TrackObject const& track) -{ - using namespace filteranalysistask; - - fhPtB->Fill(track.pt()); - fhEtaB->Fill(track.eta()); - fhPhiB->Fill(track.phi()); - fhEtaVsPhiB->Fill(track.phi(), track.eta()); - fhPtVsEtaB->Fill(track.eta(), track.pt()); - if (track.sign() > 0) { - fhPtPosB->Fill(track.pt()); - } else { - fhPtNegB->Fill(track.pt()); - } - fhDCAxyB->Fill(track.dcaXY()); - fhDCAzB->Fill(track.dcaZ()); -} - -template -void fillTrackHistosAfterSelection(TrackObject const& track) -{ - using namespace filteranalysistask; - - fhPtA->Fill(track.pt()); - fhEtaA->Fill(track.eta()); - fhPhiA->Fill(track.phi()); - fhEtaVsPhiA->Fill(track.phi(), track.eta()); - fhPtVsEtaA->Fill(track.eta(), track.pt()); - if (track.sign() > 0) { - fhPtPosA->Fill(track.pt()); - } else { - fhPtNegA->Fill(track.pt()); - } - fhDCAxyA->Fill(track.dcaXY()); - fhDCAzA->Fill(track.dcaZ()); - if (track.dcaXY() < 1.0) { - fhFineDCAxyA->Fill(track.dcaXY()); - } - if (track.dcaZ() < 1.0) { - fhFineDCAzA->Fill(track.dcaZ()); - } -} - -template -void fillParticleHistosBeforeSelection(ParticleObject const& particle, MCCollisionObject const& collision, float charge) -{ - using namespace filteranalysistask; - - fhTruePtB->Fill(particle.pt()); - fhTrueEtaB->Fill(particle.eta()); - fhTruePhiB->Fill(particle.phi()); - fhTrueEtaVsPhiB->Fill(particle.phi(), particle.eta()); - fhTruePtVsEtaB->Fill(particle.eta(), particle.pt()); - if (charge > 0) { - fhTruePtPosB->Fill(particle.pt()); - } else if (charge < 0) { - fhTruePtNegB->Fill(particle.pt()); - } - - fhTrueDCAxyB->Fill(TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + - (particle.vy() - collision.posY()) * (particle.vy() - collision.posY()))); - fhTrueDCAzB->Fill((particle.vz() - collision.posZ())); -} - -template -void fillParticleHistosAfterSelection(ParticleObject const& particle, MCCollisionObject const& collision, float charge) -{ - using namespace filteranalysistask; - - fhTruePtA->Fill(particle.pt()); - fhTrueEtaA->Fill(particle.eta()); - fhTruePhiA->Fill(particle.phi()); - fhTrueEtaVsPhiA->Fill(particle.phi(), particle.eta()); - fhTruePtVsEtaA->Fill(particle.eta(), particle.pt()); - if (charge > 0) { - fhTruePtPosA->Fill(particle.pt()); - } else { - fhTruePtNegA->Fill(particle.pt()); - } - - fhTrueDCAxyA->Fill(TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + - (particle.vy() - collision.posY()) * (particle.vy() - collision.posY()))); - fhTrueDCAzA->Fill((particle.vz() - collision.posZ())); -} - -} /* end namespace dptdptcorrelations */ - -// Task for correlations analysis -// FIXME: this should really inherit from AnalysisTask but -// we need GCC 7.4+ for that - -using namespace dptdptcorrelations; - -struct DptDptCorrelationsFilterAnalysisTask { - Configurable cfgTrackType{"trktype", 1, "Type of selected tracks: 0 = no selection, 1 = global tracks FB96"}; - Configurable cfgProcessPairs{"processpairs", true, "Process pairs: false = no, just singles, true = yes, process pairs"}; - Configurable cfgCentMultEstimator{"centmultestimator", "V0M", "Centrality/multiplicity estimator detector: V0M, NOCM: none. Default V0M"}; - Configurable cfgDataType{"datatype", "data", "Data type: data, MC, FastMC, OnTheFlyMC. Default data"}; - Configurable cfgSystem{"syst", "PbPb", "System: pp, PbPb, Pbp, pPb, XeXe. Default PbPb"}; - - Configurable cfgBinning{"binning", - {28, -7.0, 7.0, 18, 0.2, 2.0, 16, -0.8, 0.8, 72, 0.5}, - "triplets - nbins, min, max - for z_vtx, pT, eta and phi, binning plus bin fraction of phi origin shift"}; - - OutputObj fOutput{"DptDptCorrelationsGlobalInfo", OutputObjHandlingPolicy::AnalysisObject}; - - Produces acceptedevents; - Produces scannedtracks; - Produces acceptedtrueevents; - Produces scannedtruetracks; - - template - void filterTracks(TrackListObject const& ftracks, CollisionIndex colix) - { - using namespace filteranalysistask; - - int acceptedtracks = 0; - - for (auto& track : ftracks) { - /* before track selection */ - fillTrackHistosBeforeSelection(track); - - /* track selection */ - /* tricky because the boolean columns issue */ - bool asone, astwo; - AcceptTrack(track, asone, astwo); - if (asone or astwo) { - /* the track has been accepted */ - fillTrackHistosAfterSelection(track); - acceptedtracks++; - scannedtracks(colix, (uint8_t)asone, (uint8_t)astwo, track.pt(), track.eta(), track.phi()); - } - } - LOGF(DPTDPTLOGCOLLISIONS, "Accepted %d reconstructed tracks", acceptedtracks); - } - - template - void filterDetectorLevelTracks(TrackListObject const& ftracks, CollisionIndex colix) - { - using namespace filteranalysistask; - - int acceptedtracks = 0; - - for (auto& track : ftracks) { - bool asone = false; - bool astwo = false; - if (not(track.mcParticleId() < 0)) { - /* correctly reconstructed track */ - /* before track selection */ - fillTrackHistosBeforeSelection(track); - - /* track selection */ - /* tricky because the boolean columns issue */ - AcceptTrack(track, asone, astwo); - if (asone or astwo) { - /* the track has been accepted */ - fillTrackHistosAfterSelection(track); - acceptedtracks++; - scannedtracks(colix, (uint8_t)asone, (uint8_t)astwo, track.pt(), track.eta(), track.phi()); - } - } - } - LOGF(DPTDPTLOGCOLLISIONS, "Accepted %d reconstructed tracks", acceptedtracks); - } - - template - void filterParticles(ParticleListObject const& particles, MCCollisionObject const& mccollision, CollisionIndex colix) - { - using namespace filteranalysistask; - - int acceptedparticles = 0; - - for (auto& particle : particles) { - float charge = 0.0; - TParticlePDG* pdgparticle = fPDG->GetParticle(particle.pdgCode()); - if (pdgparticle != nullptr) { - charge = (pdgparticle->Charge() / 3 >= 1) ? 1.0 : ((pdgparticle->Charge() / 3 <= -1) ? -1.0 : 0.0); - } - - bool asone = false; - bool astwo = false; - if (charge != 0) { - /* before particle selection */ - fillParticleHistosBeforeSelection(particle, mccollision, charge); - - /* track selection */ - /* tricky because the boolean columns issue */ - AcceptParticle(particle, particles, asone, astwo); - if (asone or astwo) { - /* the track has been accepted */ - fillParticleHistosAfterSelection(particle, mccollision, charge); - acceptedparticles++; - scannedtruetracks(colix, (uint8_t)asone, (uint8_t)astwo, particle.pt(), particle.eta(), particle.phi()); - } - } - } - LOGF(DPTDPTLOGCOLLISIONS, "Accepted %d generated particles", acceptedparticles); - } - - void init(InitContext const&) - { - using namespace filteranalysistask; - - LOGF(info, "FilterAnalysisTask::init()"); - - /* update with the configurable values */ - /* the binning */ - ptbins = cfgBinning->mPTbins; - ptlow = cfgBinning->mPTmin; - ptup = cfgBinning->mPTmax; - etabins = cfgBinning->mEtabins; - etalow = cfgBinning->mEtamin; - etaup = cfgBinning->mEtamax; - zvtxbins = cfgBinning->mZVtxbins; - zvtxlow = cfgBinning->mZVtxmin; - zvtxup = cfgBinning->mZVtxmax; - /* the track types and combinations */ - tracktype = cfgTrackType.value; - /* the centrality/multiplicity estimation */ - if (cfgCentMultEstimator->compare("V0M") == 0) { - fCentMultEstimator = kV0M; - } else if (cfgCentMultEstimator->compare("NOCM") == 0) { - fCentMultEstimator = kNOCM; - } else { - LOGF(fatal, "Centrality/Multiplicity estimator %s not supported yet", cfgCentMultEstimator->c_str()); - } - - /* if the system type is not known at this time, we have to put the initalization somewhere else */ - fSystem = getSystemType(cfgSystem); - fDataType = getGenType(cfgDataType); - fPDG = TDatabasePDG::Instance(); - - /* create the output list which will own the task histograms */ - TList* fOutputList = new TList(); - fOutputList->SetOwner(true); - fOutput.setObject(fOutputList); - - /* incorporate configuration parameters to the output */ - fOutputList->Add(new TParameter("TrackType", cfgTrackType, 'f')); - fOutputList->Add(new TParameter("TrackOneCharge", trackonecharge, 'f')); - fOutputList->Add(new TParameter("TrackTwoCharge", tracktwocharge, 'f')); - - if ((fDataType == kData) or (fDataType == kMC)) { - /* create the reconstructed data histograms */ - if (fSystem > kPbp) { - fhCentMultB = new TH1F("CentralityB", "Centrality before cut; centrality (%)", 100, 0, 100); - fhCentMultA = new TH1F("CentralityA", "Centrality; centrality (%)", 100, 0, 100); - } else { - /* for pp, pPb and Pbp systems use multiplicity instead */ - fhCentMultB = new TH1F("MultiplicityB", "Multiplicity before cut; multiplicity (%)", 100, 0, 100); - fhCentMultA = new TH1F("MultiplicityA", "Multiplicity; multiplicity (%)", 100, 0, 100); - } - - fhVertexZB = new TH1F("VertexZB", "Vertex Z; z_{vtx}", 60, -15, 15); - fhVertexZA = new TH1F("VertexZA", "Vertex Z; z_{vtx}", zvtxbins, zvtxlow, zvtxup); - - fhPtB = new TH1F("fHistPtB", "p_{T} distribution for reconstructed before;p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhPtA = new TH1F("fHistPtA", "p_{T} distribution for reconstructed;p_{T} (GeV/c);dN/dP_{T} (c/GeV)", ptbins, ptlow, ptup); - fhPtPosB = new TH1F("fHistPtPosB", "P_{T} distribution for reconstructed (#plus) before;P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhPtPosA = new TH1F("fHistPtPosA", "P_{T} distribution for reconstructed (#plus);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", ptbins, ptlow, ptup); - fhPtNegB = new TH1F("fHistPtNegB", "P_{T} distribution for reconstructed (#minus) before;P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhPtNegA = new TH1F("fHistPtNegA", "P_{T} distribution for reconstructed (#minus);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", ptbins, ptlow, ptup); - fhEtaB = new TH1F("fHistEtaB", "#eta distribution for reconstructed before;#eta;counts", 40, -2.0, 2.0); - fhEtaA = new TH1F("fHistEtaA", "#eta distribution for reconstructed;#eta;counts", etabins, etalow, etaup); - fhPhiB = new TH1F("fHistPhiB", "#phi distribution for reconstructed before;#phi;counts", 360, 0.0, 2 * M_PI); - fhPhiA = new TH1F("fHistPhiA", "#phi distribution for reconstructed;#phi;counts", 360, 0.0, 2 * M_PI); - fhEtaVsPhiB = new TH2F(TString::Format("CSTaskEtaVsPhiB_%s", fTaskConfigurationString.c_str()), "#eta vs #phi before;#phi;#eta", 360, 0.0, 2 * M_PI, 100, -2.0, 2.0); - fhEtaVsPhiA = new TH2F(TString::Format("CSTaskEtaVsPhiA_%s", fTaskConfigurationString.c_str()), "#eta vs #phi;#phi;#eta", 360, 0.0, 2 * M_PI, etabins, etalow, etaup); - fhPtVsEtaB = new TH2F(TString::Format("fhPtVsEtaB_%s", fTaskConfigurationString.c_str()), "p_{T} vs #eta before;#eta;p_{T} (GeV/c)", etabins, etalow, etaup, 100, 0.0, 15.0); - fhPtVsEtaA = new TH2F(TString::Format("fhPtVsEtaA_%s", fTaskConfigurationString.c_str()), "p_{T} vs #eta;#eta;p_{T} (GeV/c)", etabins, etalow, etaup, ptbins, ptlow, ptup); - fhDCAxyB = new TH1F("DCAxyB", "DCA_{xy} distribution for reconstructed before;DCA_{xy} (cm);counts", 1000, -4.0, 4.0); - fhDCAxyA = new TH1F("DCAxyA", "DCA_{xy} distribution for reconstructed;DCA_{xy} (cm);counts", 1000, -4., 4.0); - fhFineDCAxyA = new TH1F("FineDCAxyA", "DCA_{xy} distribution for reconstructed;DCA_{xy} (cm);counts", 4000, -1.0, 1.0); - fhDCAzB = new TH1F("DCAzB", "DCA_{z} distribution for reconstructed before;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - fhDCAzA = new TH1F("DCAzA", "DCA_{z} distribution for reconstructed;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - fhFineDCAzA = new TH1F("FineDCAzA", "DCA_{z} distribution for reconstructed;DCA_{z} (cm);counts", 4000, -1.0, 1.0); - - /* add the hstograms to the output list */ - fOutputList->Add(fhCentMultB); - fOutputList->Add(fhCentMultA); - fOutputList->Add(fhVertexZB); - fOutputList->Add(fhVertexZA); - fOutputList->Add(fhPtB); - fOutputList->Add(fhPtA); - fOutputList->Add(fhPtPosB); - fOutputList->Add(fhPtPosA); - fOutputList->Add(fhPtNegB); - fOutputList->Add(fhPtNegA); - fOutputList->Add(fhEtaB); - fOutputList->Add(fhEtaA); - fOutputList->Add(fhPhiB); - fOutputList->Add(fhPhiA); - fOutputList->Add(fhEtaVsPhiB); - fOutputList->Add(fhEtaVsPhiA); - fOutputList->Add(fhPtVsEtaB); - fOutputList->Add(fhPtVsEtaA); - fOutputList->Add(fhDCAxyB); - fOutputList->Add(fhDCAxyA); - fOutputList->Add(fhFineDCAxyA); - fOutputList->Add(fhDCAzB); - fOutputList->Add(fhDCAzA); - fOutputList->Add(fhFineDCAzA); - } - - if (fDataType != kData) { - /* create the true data histograms */ - if (fSystem > kPbp) { - fhTrueCentMultB = new TH1F("TrueCentralityB", "Centrality before (truth); centrality (%)", 100, 0, 100); - fhTrueCentMultA = new TH1F("TrueCentralityA", "Centrality (truth); centrality (%)", 100, 0, 100); - } else { - /* for pp, pPb and Pbp systems use multiplicity instead */ - fhTrueCentMultB = new TH1F("TrueMultiplicityB", "Multiplicity before (truth); multiplicity (%)", 100, 0, 100); - fhTrueCentMultA = new TH1F("TrueMultiplicityA", "Multiplicity (truth); multiplicity (%)", 100, 0, 100); - } - - fhTrueVertexZB = new TH1F("TrueVertexZB", "Vertex Z before (truth); z_{vtx}", 60, -15, 15); - fhTrueVertexZA = new TH1F("TrueVertexZA", "Vertex Z (truth); z_{vtx}", zvtxbins, zvtxlow, zvtxup); - - fhTruePtB = new TH1F("fTrueHistPtB", "p_{T} distribution before (truth);p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhTruePtA = new TH1F("fTrueHistPtA", "p_{T} distribution (truth);p_{T} (GeV/c);dN/dP_{T} (c/GeV)", ptbins, ptlow, ptup); - fhTruePtPosB = new TH1F("fTrueHistPtPosB", "P_{T} distribution (#plus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhTruePtPosA = new TH1F("fTrueHistPtPosA", "P_{T} distribution (#plus) (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", ptbins, ptlow, ptup); - fhTruePtNegB = new TH1F("fTrueHistPtNegB", "P_{T} distribution (#minus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhTruePtNegA = new TH1F("fTrueHistPtNegA", "P_{T} distribution (#minus) (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", ptbins, ptlow, ptup); - fhTrueEtaB = new TH1F("fTrueHistEtaB", "#eta distribution before (truth);#eta;counts", 40, -2.0, 2.0); - fhTrueEtaA = new TH1F("fTrueHistEtaA", "#eta distribution (truth);#eta;counts", etabins, etalow, etaup); - fhTruePhiB = new TH1F("fTrueHistPhiB", "#phi distribution before (truth);#phi;counts", 360, 0.0, 2 * M_PI); - fhTruePhiA = new TH1F("fTrueHistPhiA", "#phi distribution (truth);#phi;counts", 360, 0.0, 2 * M_PI); - fhTrueEtaVsPhiB = new TH2F(TString::Format("CSTaskTrueEtaVsPhiB_%s", fTaskConfigurationString.c_str()), "#eta vs #phi before (truth);#phi;#eta", 360, 0.0, 2 * M_PI, 100, -2.0, 2.0); - fhTrueEtaVsPhiA = new TH2F(TString::Format("CSTaskTrueEtaVsPhiA_%s", fTaskConfigurationString.c_str()), "#eta vs #phi (truth);#phi;#eta", 360, 0.0, 2 * M_PI, etabins, etalow, etaup); - fhTruePtVsEtaB = new TH2F(TString::Format("fhTruePtVsEtaB_%s", fTaskConfigurationString.c_str()), "p_{T} vs #eta before (truth);#eta;p_{T} (GeV/c)", etabins, etalow, etaup, 100, 0.0, 15.0); - fhTruePtVsEtaA = new TH2F(TString::Format("fhTruePtVsEtaA_%s", fTaskConfigurationString.c_str()), "p_{T} vs #eta (truth);#eta;p_{T} (GeV/c)", etabins, etalow, etaup, ptbins, ptlow, ptup); - fhTrueDCAxyB = new TH1F("TrueDCAxyB", "DCA_{xy} distribution for generated before;DCA_{xy} (cm);counts", 1000, -4.0, 4.0); - fhTrueDCAxyA = new TH1F("TrueDCAxyA", "DCA_{xy} distribution for generated;DCA_{xy};counts (cm)", 1000, -4., 4.0); - fhTrueDCAzB = new TH1F("TrueDCAzB", "DCA_{z} distribution for generated before;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - fhTrueDCAzA = new TH1F("TrueDCAzA", "DCA_{z} distribution for generated;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - - /* add the hstograms to the output list */ - fOutputList->Add(fhTrueCentMultB); - fOutputList->Add(fhTrueCentMultA); - fOutputList->Add(fhTrueVertexZB); - fOutputList->Add(fhTrueVertexZA); - fOutputList->Add(fhTruePtB); - fOutputList->Add(fhTruePtA); - fOutputList->Add(fhTruePtPosB); - fOutputList->Add(fhTruePtPosA); - fOutputList->Add(fhTruePtNegB); - fOutputList->Add(fhTruePtNegA); - fOutputList->Add(fhTrueEtaB); - fOutputList->Add(fhTrueEtaA); - fOutputList->Add(fhTruePhiB); - fOutputList->Add(fhTruePhiA); - fOutputList->Add(fhTrueEtaVsPhiB); - fOutputList->Add(fhTrueEtaVsPhiA); - fOutputList->Add(fhTruePtVsEtaB); - fOutputList->Add(fhTruePtVsEtaA); - fOutputList->Add(fhTrueDCAxyB); - fOutputList->Add(fhTrueDCAxyA); - fOutputList->Add(fhTrueDCAzB); - fOutputList->Add(fhTrueDCAzA); - } - } - - void processWithCent(aod::CollisionEvSelCent const& collision, soa::Join const& ftracks) - { - using namespace filteranalysistask; - - LOGF(DPTDPTLOGCOLLISIONS, "FilterAnalysisTask::processWithCent(). New collision with %d tracks", ftracks.size()); - - fhCentMultB->Fill(collision.centV0M()); - fhVertexZB->Fill(collision.posZ()); - bool acceptedevent = false; - float centormult = -100.0; - if (IsEvtSelected(collision, centormult)) { - acceptedevent = true; - fhCentMultA->Fill(collision.centV0M()); - fhVertexZA->Fill(collision.posZ()); - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - - filterTracks(ftracks, acceptedevents.lastIndex()); - } else { - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& track : ftracks) { - scannedtracks(acceptedevents.lastIndex(), (uint8_t) false, (uint8_t) false, track.pt(), track.eta(), track.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithCent, "Process reco with centrality", false); - - void processWithoutCent(aod::CollisionEvSel const& collision, soa::Join const& ftracks) - { - using namespace filteranalysistask; - - LOGF(DPTDPTLOGCOLLISIONS, "FilterAnalysisTask::processWithoutCent(). New collision with collision id %d and with %d tracks", collision.bcId(), ftracks.size()); - - /* the task does not have access to either centrality nor multiplicity - classes information, so it has to live without it. - For the time being we assign a value of 50% */ - fhCentMultB->Fill(50.0); - fhVertexZB->Fill(collision.posZ()); - bool acceptedevent = false; - float centormult = -100.0; - if (IsEvtSelectedNoCentMult(collision, centormult)) { - acceptedevent = true; - fhCentMultA->Fill(50.0); - fhVertexZA->Fill(collision.posZ()); - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - - filterTracks(ftracks, acceptedevents.lastIndex()); - } else { - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& track : ftracks) { - scannedtracks(acceptedevents.lastIndex(), (uint8_t) false, (uint8_t) false, track.pt(), track.eta(), track.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithoutCent, "Process reco without centrality", false); - - void processWithCentDetectorLevel(aod::CollisionEvSelCent const& collision, - soa::Join const& ftracks) - { - using namespace filteranalysistask; - - LOGF(DPTDPTLOGCOLLISIONS, "FilterAnalysisTask::processWithCentDetectorLevel(). New collision with %d tracks", ftracks.size()); - - fhCentMultB->Fill(collision.centV0M()); - fhVertexZB->Fill(collision.posZ()); - bool acceptedevent = false; - float centormult = -100.0; - if (IsEvtSelected(collision, centormult)) { - acceptedevent = true; - fhCentMultA->Fill(collision.centV0M()); - fhVertexZA->Fill(collision.posZ()); - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - - filterDetectorLevelTracks(ftracks, acceptedevents.lastIndex()); - } else { - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& track : ftracks) { - scannedtracks(acceptedevents.lastIndex(), (uint8_t) false, (uint8_t) false, track.pt(), track.eta(), track.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithCentDetectorLevel, "Process MC detector level with centrality", false); - - void processWithoutCentDetectorLevel(aod::CollisionEvSel const& collision, - soa::Join const& ftracks) - { - using namespace filteranalysistask; - - LOGF(DPTDPTLOGCOLLISIONS, "FilterAnalysisTask::processWithoutCentDetectorLevel(). New collision with collision id %d and with %d tracks", collision.bcId(), ftracks.size()); - - /* the task does not have access to either centrality nor multiplicity - classes information, so it has to live without it. - For the time being we assign a value of 50% */ - fhCentMultB->Fill(50.0); - fhVertexZB->Fill(collision.posZ()); - bool acceptedevent = false; - float centormult = -100.0; - if (IsEvtSelectedNoCentMult(collision, centormult)) { - acceptedevent = true; - fhCentMultA->Fill(50.0); - fhVertexZA->Fill(collision.posZ()); - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - - filterDetectorLevelTracks(ftracks, acceptedevents.lastIndex()); - } else { - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& track : ftracks) { - scannedtracks(acceptedevents.lastIndex(), (uint8_t) false, (uint8_t) false, track.pt(), track.eta(), track.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithoutCentDetectorLevel, "Process MC detector level without centrality", false); - - void processWithCentGeneratorLevel(aod::McCollision const& mccollision, - soa::Join const& collisions, - aod::McParticles const& mcparticles) - { - using namespace filteranalysistask; - - LOGF(DPTDPTLOGCOLLISIONS, "FilterAnalysisTask::processWithCentGeneratorLevel(). New generated collision %d reconstructed collisions and %d particles", collisions.size(), mcparticles.size()); - - /* TODO: in here we have to decide what to do in the following cases - - On the fly production -> clearly we will need a different process - - reconstructed collisions without generated associated -> we need a different process or a different signature - - multiplicity/centrality classes extracted from the reconstructed collision but then - - generated collision without associated reconstructed collision: how to extract mutliplicity/centrality classes? - - generated collision with several associated reconstructed collisions: from which to extract multiplicity/centrality classes? - */ - if (collisions.size() > 1) { - LOGF(error, "FilterAnalysisTask::processWithCentGeneratorLevel(). Generated collision with more than one reconstructed collisions. Processing only the first for centrality/multiplicity classes extraction"); - } - - for (auto& collision : collisions) { - float cent = collision.centV0M(); - fhTrueCentMultB->Fill(cent); - fhTrueVertexZB->Fill(mccollision.posZ()); - - bool acceptedevent = false; - if (IsTrueEvtSelected(mccollision, cent)) { - acceptedevent = true; - fhTrueCentMultA->Fill(cent); - fhTrueVertexZA->Fill(mccollision.posZ()); - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, cent); - - filterParticles(mcparticles, collision, acceptedtrueevents.lastIndex()); - } else { - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, cent); - for (auto& particle : mcparticles) { - scannedtruetracks(acceptedtrueevents.lastIndex(), (uint8_t) false, (uint8_t) false, particle.pt(), particle.eta(), particle.phi()); - } - } - break; /* TODO: only processing the first reconstructed collision for centrality/multiplicity class estimation */ - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithCentGeneratorLevel, "Process generated with centrality", false); - - void processWithoutCentGeneratorLevel(aod::McCollision const& mccollision, - aod::McParticles const& mcparticles) - { - using namespace filteranalysistask; - - LOGF(DPTDPTLOGCOLLISIONS, "FilterAnalysisTask::processWithoutCentGeneratorLevel(). New generated collision with %d particles", mcparticles.size()); - - /* the task does not have access to either centrality nor multiplicity - classes information, so it has to live without it. - For the time being we assign a value of 50% */ - fhTrueCentMultB->Fill(50.0); - fhTrueVertexZB->Fill(mccollision.posZ()); - - bool acceptedevent = false; - float centormult = 50.0; - if (IsTrueEvtSelected(mccollision, centormult)) { - acceptedevent = true; - fhTrueCentMultA->Fill(centormult); - fhTrueVertexZA->Fill(mccollision.posZ()); - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, centormult); - - filterParticles(mcparticles, mccollision, acceptedtrueevents.lastIndex()); - } else { - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& particle : mcparticles) { - scannedtruetracks(acceptedtrueevents.lastIndex(), (uint8_t) false, (uint8_t) false, particle.pt(), particle.eta(), particle.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithoutCentGeneratorLevel, "Process generated without centrality", false); -}; - // Task for building correlations struct DptDptCorrelationsTask { @@ -1048,8 +117,10 @@ struct DptDptCorrelationsTask { /// \return the track phi origin shifted azimuthal angle float GetShiftedPhi(float phi) { + using namespace correlationstask; + using namespace o2::analysis::dptdptfilter; if (not(phi < phiup)) { - return phi - M_PI * 2; + return phi - constants::math::TwoPI; } else { return phi; } @@ -1070,6 +141,9 @@ struct DptDptCorrelationsTask { template int GetEtaPhiIndex(TrackObject const& t) { + using namespace correlationstask; + using namespace o2::analysis::dptdptfilter; + int etaix = int((t.eta() - etalow) / etabinwidth); /* consider a potential phi origin shift */ float phi = GetShiftedPhi(t.phi()); @@ -1090,6 +164,7 @@ struct DptDptCorrelationsTask { int GetDEtaDPhiGlobalIndex(TrackObject const& t1, TrackObject const& t2) { using namespace correlationstask; + using namespace o2::analysis::dptdptfilter; /* rule: ix are always zero based while bins are always one based */ int etaix_1 = int((t1.eta() - etalow) / etabinwidth); @@ -1236,6 +311,7 @@ struct DptDptCorrelationsTask { void init(TList* fOutputList) { using namespace correlationstask; + using namespace o2::analysis::dptdptfilter; /* create the histograms */ Bool_t oldstatus = TH1::AddDirectoryStatus(); @@ -1369,6 +445,7 @@ struct DptDptCorrelationsTask { void init(InitContext const&) { using namespace correlationstask; + using namespace o2::analysis::dptdptfilter; /* update with the configurable values */ ptbins = cfgBinning->mPTbins; @@ -1382,7 +459,7 @@ struct DptDptCorrelationsTask { zvtxup = cfgBinning->mZVtxmax; phibins = cfgBinning->mPhibins; philow = 0.0f; - phiup = M_PI * 2; + phiup = constants::math::TwoPI; phibinshift = cfgBinning->mPhibinshift; processpairs = cfgProcessPairs.value; /* update the potential binning change */ @@ -1394,9 +471,9 @@ struct DptDptCorrelationsTask { deltaetalow = etalow - etaup, deltaetaup = etaup - etalow; deltaetabinwidth = (deltaetaup - deltaetalow) / float(deltaetabins); deltaphibins = phibins; - deltaphibinwidth = M_PI * 2 / deltaphibins; + deltaphibinwidth = constants::math::TwoPI / deltaphibins; deltaphilow = 0.0 - deltaphibinwidth / 2.0; - deltaphiup = M_PI * 2 - deltaphibinwidth / 2.0; + deltaphiup = constants::math::TwoPI - deltaphibinwidth / 2.0; /* create the output directory which will own the task output */ TList* fGlobalOutputList = new TList(); @@ -1470,18 +547,18 @@ struct DptDptCorrelationsTask { return ixDCE; } - Filter onlyacceptedevents = (aod::dptdptcorrelations::eventaccepted == (uint8_t) true); - Filter onlyacceptedtracks = ((aod::dptdptcorrelations::trackacceptedasone == (uint8_t) true) or (aod::dptdptcorrelations::trackacceptedastwo == (uint8_t) true)); + Filter onlyacceptedcollisions = (aod::dptdptfilter::collisionaccepted == true); + Filter onlyacceptedtracks = ((aod::dptdptfilter::trackacceptedasone == true) or (aod::dptdptfilter::trackacceptedastwo == true)); - void processRecLevel(soa::Filtered::iterator const& collision, soa::Filtered& tracks) + void processRecLevel(soa::Filtered::iterator const& collision, soa::Filtered& tracks) { using namespace correlationstask; /* locate the data collecting engine for the collision centrality/multiplicity */ int ixDCE = getDCEindex(collision); if (not(ixDCE < 0)) { - Partition TracksOne = aod::dptdptcorrelations::trackacceptedasone == (uint8_t) true; - Partition TracksTwo = aod::dptdptcorrelations::trackacceptedastwo == (uint8_t) true; + Partition TracksOne = aod::dptdptfilter::trackacceptedasone == true; + Partition TracksTwo = aod::dptdptfilter::trackacceptedastwo == true; TracksOne.bindTable(tracks); TracksTwo.bindTable(tracks); @@ -1492,15 +569,15 @@ struct DptDptCorrelationsTask { } PROCESS_SWITCH(DptDptCorrelationsTask, processRecLevel, "Process reco level correlations", false); - void processGenLevel(soa::Filtered::iterator const& collision, soa::Filtered& tracks) + void processGenLevel(soa::Filtered::iterator const& collision, soa::Filtered& tracks) { using namespace correlationstask; /* locate the data collecting engine for the collision centrality/multiplicity */ int ixDCE = getDCEindex(collision); if (not(ixDCE < 0)) { - Partition TracksOne = aod::dptdptcorrelations::trackacceptedasone == (uint8_t) true; - Partition TracksTwo = aod::dptdptcorrelations::trackacceptedastwo == (uint8_t) true; + Partition TracksOne = aod::dptdptfilter::trackacceptedasone == true; + Partition TracksTwo = aod::dptdptfilter::trackacceptedastwo == true; TracksOne.bindTable(tracks); TracksTwo.bindTable(tracks); @@ -1520,138 +597,9 @@ struct DptDptCorrelationsTask { PROCESS_SWITCH(DptDptCorrelationsTask, processCleaner, "Cleaner process for not used output", false); }; -// Checking the filtered tables -/* it seems we cannot use a base class task */ -// struct TracksAndEventClassificationQABase { - -void initQATask(InitContext const&, TList* outlst) -{ - using namespace filteranalysistaskqa; - - fhTracksOne = new TH1F("TracksOne", "Tracks as track one;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksTwo = new TH1F("TracksTwo", "Tracks as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksOneAndTwo = new TH1F("TracksOneAndTwo", "Tracks as track one and as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksNone = new TH1F("TracksNone", "Not selected tracks;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksOneUnsel = new TH1F("TracksOneUnsel", "Tracks as track one;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksTwoUnsel = new TH1F("TracksTwoUnsel", "Tracks as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksOneAndTwoUnsel = new TH1F("TracksOneAndTwoUnsel", "Tracks as track one and as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksNoneUnsel = new TH1F("TracksNoneUnsel", "Not selected tracks;number of tracks;events", 1500, 0.0, 1500.0); - fhSelectedEvents = new TH1F("SelectedEvents", "Selected events;;events", 2, 0.0, 2.0); - fhSelectedEvents->GetXaxis()->SetBinLabel(1, "Not selected events"); - fhSelectedEvents->GetXaxis()->SetBinLabel(2, "Selected events"); - - outlst->Add(fhTracksOne); - outlst->Add(fhTracksTwo); - outlst->Add(fhTracksOneAndTwo); - outlst->Add(fhTracksNone); - outlst->Add(fhTracksOneUnsel); - outlst->Add(fhTracksTwoUnsel); - outlst->Add(fhTracksOneAndTwoUnsel); - outlst->Add(fhTracksNoneUnsel); - outlst->Add(fhSelectedEvents); -} - -template -void processQATask(FilteredCollision const& collision, - FilteredTracks const& tracks) -{ - using namespace filteranalysistaskqa; - - if (collision.eventaccepted() != (uint8_t) true) { - fhSelectedEvents->Fill(0.5); - } else { - fhSelectedEvents->Fill(1.5); - } - - int ntracks_one = 0; - int ntracks_two = 0; - int ntracks_one_and_two = 0; - int ntracks_none = 0; - for (auto& track : tracks) { - if ((track.trackacceptedasone() != (uint8_t) true) and (track.trackacceptedastwo() != (uint8_t) true)) { - ntracks_none++; - } - if ((track.trackacceptedasone() == (uint8_t) true) and (track.trackacceptedastwo() == (uint8_t) true)) { - ntracks_one_and_two++; - } - if (track.trackacceptedasone() == (uint8_t) true) { - ntracks_one++; - } - if (track.trackacceptedastwo() == (uint8_t) true) { - ntracks_two++; - } - } - if (collision.eventaccepted() != (uint8_t) true) { - /* control for non selected events */ - fhTracksOneUnsel->Fill(ntracks_one); - fhTracksTwoUnsel->Fill(ntracks_two); - fhTracksNoneUnsel->Fill(ntracks_none); - fhTracksOneAndTwoUnsel->Fill(ntracks_one_and_two); - } else { - fhTracksOne->Fill(ntracks_one); - fhTracksTwo->Fill(ntracks_two); - fhTracksNone->Fill(ntracks_none); - fhTracksOneAndTwo->Fill(ntracks_one_and_two); - } -} -// }; - -/* it seems we cannot use a base class task */ -// struct TracksAndEventClassificationQARec : TracksAndEventClassificationQABase { -struct TracksAndEventClassificationQARec { - OutputObj fOutput{"FliterTaskRecoQA", OutputObjHandlingPolicy::AnalysisObject}; - - void init(InitContext const& context) - { - TList* fOutputList = new TList(); - fOutputList->SetName("FilterTaskRecoQA"); - fOutputList->SetOwner(true); - fOutput.setObject(fOutputList); - - initQATask(context, fOutputList); - } - - Filter onlyacceptedevents = (aod::dptdptcorrelations::eventaccepted == (uint8_t) true); - Filter onlyacceptedtracks = ((aod::dptdptcorrelations::trackacceptedasone == (uint8_t) true) or (aod::dptdptcorrelations::trackacceptedastwo == (uint8_t) true)); - - void process(soa::Filtered::iterator const& collision, soa::Filtered const& tracks) - { - LOGF(DPTDPTLOGCOLLISIONS, "New filtered collision with BC id %d and with %d accepted tracks", collision.bcId(), tracks.size()); - processQATask(collision, tracks); - } -}; - -/* it seems we cannot use a base class task */ -//struct TracksAndEventClassificationQAGen : TracksAndEventClassificationQABase { -struct TracksAndEventClassificationQAGen { - OutputObj fOutput{"FliterTaskGenQA", OutputObjHandlingPolicy::AnalysisObject}; - - void init(InitContext const& context) - { - TList* fOutputList = new TList(); - fOutputList->SetName("FilterTaskGenQA"); - fOutputList->SetOwner(true); - fOutput.setObject(fOutputList); - - initQATask(context, fOutputList); - } - - Filter onlyacceptedevents = (aod::dptdptcorrelations::eventaccepted == (uint8_t) true); - Filter onlyacceptedtracks = ((aod::dptdptcorrelations::trackacceptedasone == (uint8_t) true) or (aod::dptdptcorrelations::trackacceptedastwo == (uint8_t) true)); - - void process(soa::Filtered::iterator const& collision, soa::Filtered const& tracks) - { - LOGF(DPTDPTLOGCOLLISIONS, "New filtered generated collision with BC id %d and with %d accepted tracks", collision.bcId(), tracks.size()); - processQATask(collision, tracks); - } -}; - WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { WorkflowSpec workflow{ - adaptAnalysisTask(cfgc, SetDefaultProcesses{{{"processWithoutCent", true}, {"processWithoutCentMC", true}}}), - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc, TaskName{"DptDptCorrelationsTaskRec"}, SetDefaultProcesses{{{"processRecLevel", true}, {"processCleaner", false}}}), adaptAnalysisTask(cfgc, TaskName{"DptDptCorrelationsTaskGen"}, SetDefaultProcesses{{{"processGenLevel", true}, {"processCleaner", false}}})}; return workflow; diff --git a/PWGCF/Tasks/dptdptfilterqa.cxx b/PWGCF/Tasks/dptdptfilterqa.cxx new file mode 100644 index 00000000000..52c2b5352c9 --- /dev/null +++ b/PWGCF/Tasks/dptdptfilterqa.cxx @@ -0,0 +1,161 @@ +// 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. + +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/runDataProcessing.h" +#include "PWGCF/DataModel/DptDptFiltered.h" +#include "PWGCF/TableProducer/dptdptfilter.h" + +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::soa; +using namespace o2::framework::expressions; + +#define DPTDPTFILTERLOGCOLLISIONS debug +#define DPTDPTFILTERLOGTRACKS debug + +namespace o2 +{ +namespace aod +{ +using FilteredTracks = soa::Filtered>; +using FilteredTrackData = Partition::filtered_iterator; +} // namespace aod +} // namespace o2 + +namespace o2::analysis::dptdptfilterqa +{ +typedef enum { kRECO = 0, + kGEN } innerdatatype; +static constexpr std::string_view dirname[] = {"reconstructed/", "generated/"}; +} // namespace o2::analysis::dptdptfilterqa + +// Checking the filtered tables +struct DptDptFilterQA { + Configurable cfgDataType{"datatype", "data", "Data type: data, MC, FastMC, OnTheFlyMC. Default data"}; + HistogramRegistry histos{"DptDptFilterQA", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + o2::analysis::dptdptfilter::DataType datatype; + + template + void createHistograms() + { + + using namespace o2::analysis::dptdptfilterqa; + + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksOne").Data(), "Tracks as track one", kTH1F, {{1500, 0.0, 1500.0, "number of tracks"}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksTwo").Data(), "Tracks as track two", kTH1F, {{1500, 0.0, 1500.0, "number of tracks"}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksOneAndTwo").Data(), "Tracks as track one and as track two", kTH1F, {{1500, 0.0, 1500.0, "number of tracks"}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksNone").Data(), "Not selected tracks", kTH1F, {{1500, 0.0, 1500.0, "number of tracks"}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksOneUnsel").Data(), "Tracks as track one", kTH1F, {{1500, 0.0, 1500.0, "number of tracks"}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksTwoUnsel").Data(), "Tracks as track two", kTH1F, {{1500, 0.0, 1500.0, "number of tracks"}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksOneAndTwoUnsel").Data(), "Tracks as track one and as track two", kTH1F, {{1500, 0.0, 1500.0, "number of tracks"}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "TracksNoneUnsel").Data(), "Not selected tracks", kTH1F, {{1500, 0.0, 1500.0}}); + histos.add(TString::Format("%s%s", dirname[dir].data(), "SelectedEvents").Data(), "Selected events", kTH1F, {{2, 0.0, 2.0}}); + histos.get(HIST(dirname[dir]) + HIST("SelectedEvents"))->GetXaxis()->SetBinLabel(1, "Not selected events"); + histos.get(HIST(dirname[dir]) + HIST("SelectedEvents"))->GetXaxis()->SetBinLabel(2, "Selected events"); + }; + + void init(InitContext const& context) + { + using namespace o2::analysis::dptdptfilter; + using namespace o2::analysis::dptdptfilterqa; + + switch (getDataType(cfgDataType)) { + case kData: + createHistograms(); + break; + case kMC: + createHistograms(); + createHistograms(); + break; + case kFastMC: + case kOnTheFly: + createHistograms(); + break; + default: + LOGF(fatal, "Data type %s not supported", cfgDataType); + break; + } + } + + template + void processQATask(FilteredCollision const& collision, + FilteredTracks const& tracks) + { + using namespace o2::analysis::dptdptfilterqa; + + if (!collision.collisionaccepted()) { + histos.fill(HIST(dirname[dir]) + HIST("SelectedEvents"), 0.5); + } else { + histos.fill(HIST(dirname[dir]) + HIST("SelectedEvents"), 1.5); + } + + int ntracks_one = 0; + int ntracks_two = 0; + int ntracks_one_and_two = 0; + int ntracks_none = 0; + for (auto& track : tracks) { + if (!track.trackacceptedasone() and !track.trackacceptedastwo()) { + ntracks_none++; + } + if (track.trackacceptedasone() and track.trackacceptedastwo()) { + ntracks_one_and_two++; + } + if (track.trackacceptedasone()) { + ntracks_one++; + } + if (track.trackacceptedastwo()) { + ntracks_two++; + } + } + if (!collision.collisionaccepted()) { + /* control for non selected events */ + histos.fill(HIST(dirname[dir]) + HIST("TracksOneUnsel"), ntracks_one); + histos.fill(HIST(dirname[dir]) + HIST("TracksTwoUnsel"), ntracks_two); + histos.fill(HIST(dirname[dir]) + HIST("TracksNoneUnsel"), ntracks_none); + histos.fill(HIST(dirname[dir]) + HIST("TracksOneAndTwoUnsel"), ntracks_one_and_two); + } else { + histos.fill(HIST(dirname[dir]) + HIST("TracksOne"), ntracks_one); + histos.fill(HIST(dirname[dir]) + HIST("TracksTwo"), ntracks_two); + histos.fill(HIST(dirname[dir]) + HIST("TracksNone"), ntracks_none); + histos.fill(HIST(dirname[dir]) + HIST("TracksOneAndTwo"), ntracks_one_and_two); + } + } + + Filter onlyacceptedcollisions = (aod::dptdptfilter::collisionaccepted == true); + Filter onlyacceptedtracks = ((aod::dptdptfilter::trackacceptedasone == true) or (aod::dptdptfilter::trackacceptedastwo == true)); + + void processGeneratorLevel(soa::Filtered::iterator const& collision, soa::Filtered const& tracks) + { + using namespace o2::analysis::dptdptfilterqa; + LOGF(DPTDPTFILTERLOGCOLLISIONS, "New filtered generated collision with BC id %d and with %d accepted tracks", collision.bcId(), tracks.size()); + processQATask(collision, tracks); + } + PROCESS_SWITCH(DptDptFilterQA, processGeneratorLevel, "Process generator level filter task QA", true); + + void processDetectorLevel(soa::Filtered::iterator const& collision, soa::Filtered const& tracks) + { + using namespace o2::analysis::dptdptfilterqa; + LOGF(DPTDPTFILTERLOGCOLLISIONS, "New filtered collision with BC id %d and with %d accepted tracks", collision.bcId(), tracks.size()); + processQATask(collision, tracks); + } + PROCESS_SWITCH(DptDptFilterQA, processDetectorLevel, "Process detector level filter task QA", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; + return workflow; +} diff --git a/PWGCF/Tasks/match-reco-gen.cxx b/PWGCF/Tasks/match-reco-gen.cxx index 9f0b5fb74b2..4f2965697cd 100644 --- a/PWGCF/Tasks/match-reco-gen.cxx +++ b/PWGCF/Tasks/match-reco-gen.cxx @@ -19,7 +19,10 @@ #include "Common/Core/MC.h" #include "Common/Core/PID/PIDResponse.h" #include "PWGCF/Core/AnalysisConfigurableCuts.h" +#include "PWGCF/DataModel/DptDptFiltered.h" +#include "PWGCF/TableProducer/dptdptfilter.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "Framework/runDataProcessing.h" #include #include #include @@ -41,688 +44,26 @@ using namespace o2::framework::expressions; #define MATCHRECGENLOGCOLLISIONS debug #define MATCHRECGENLOGTRACKS debug -#include "Framework/runDataProcessing.h" - namespace o2 { namespace aod { -/* we have to change from int to bool when bool columns work properly */ -namespace dptdptcorrelations -{ -DECLARE_SOA_COLUMN(EventAccepted, eventaccepted, uint8_t); //! If the collision/event has been accepted or not -DECLARE_SOA_COLUMN(EventCentMult, centmult, float); //! The centrality/multiplicity pecentile -} // namespace dptdptcorrelations -DECLARE_SOA_TABLE(AcceptedEvents, "AOD", "ACCEPTEDEVENTS", //! Accepted reconstructed collisions/events filtered table - o2::soa::Index<>, - collision::BCId, - collision::PosZ, - dptdptcorrelations::EventAccepted, - dptdptcorrelations::EventCentMult); -using AcceptedEvent = AcceptedEvents::iterator; -DECLARE_SOA_TABLE(AcceptedTrueEvents, "AOD", "ACCTRUEEVENTS", //! Accepted generated collisions/events filtered table - o2::soa::Index<>, - collision::BCId, - mccollision::PosZ, - dptdptcorrelations::EventAccepted, - dptdptcorrelations::EventCentMult); -using AcceptedTrueEvent = AcceptedTrueEvents::iterator; -namespace dptdptcorrelations -{ -DECLARE_SOA_INDEX_COLUMN(AcceptedEvent, event); //! Reconstructed collision/event -DECLARE_SOA_INDEX_COLUMN(AcceptedTrueEvent, mcevent); //! Generated collision/event -DECLARE_SOA_COLUMN(TrackacceptedAsOne, trackacceptedasone, uint8_t); //! Track accepted as type one -DECLARE_SOA_COLUMN(TrackacceptedAsTwo, trackacceptedastwo, uint8_t); //! Track accepted as type two -DECLARE_SOA_COLUMN(Pt, pt, float); //! The track transverse momentum -DECLARE_SOA_COLUMN(Eta, eta, float); //! The track pseudorapidity -DECLARE_SOA_COLUMN(Phi, phi, float); //! The track azimuthal angle -} // namespace dptdptcorrelations -DECLARE_SOA_TABLE(ScannedTracks, "AOD", "SCANNEDTRACKS", //! The reconstructed tracks filtered table - dptdptcorrelations::AcceptedEventId, - dptdptcorrelations::TrackacceptedAsOne, - dptdptcorrelations::TrackacceptedAsTwo, - dptdptcorrelations::Pt, - dptdptcorrelations::Eta, - dptdptcorrelations::Phi); -DECLARE_SOA_TABLE(ScannedTrueTracks, "AOD", "SCANTRUETRACKS", //! The generated particles filtered table - dptdptcorrelations::AcceptedTrueEventId, - dptdptcorrelations::TrackacceptedAsOne, - dptdptcorrelations::TrackacceptedAsTwo, - dptdptcorrelations::Pt, - dptdptcorrelations::Eta, - dptdptcorrelations::Phi); - -using CollisionsEvSelCent = soa::Join; -using CollisionEvSelCent = soa::Join::iterator; -using CollisionsEvSel = soa::Join; -using CollisionEvSel = soa::Join::iterator; -using TrackData = soa::Join::iterator; using FilteredTracks = soa::Filtered>; using FilteredTrackData = Partition::filtered_iterator; } // namespace aod } // namespace o2 -namespace dptdptcorrelations -{ -/* all this is made configurable */ -int ptbins = 18; -float ptlow = 0.2, ptup = 2.0; -int etabins = 16; -float etalow = -0.8, etaup = 0.8; -int zvtxbins = 40; -float zvtxlow = -10.0, zvtxup = 10.0; -int phibins = 72; -float philow = 0.0; -float phiup = M_PI * 2; -float phibinshift = 0.5; -float etabinwidth = (etaup - etalow) / float(etabins); -float phibinwidth = (phiup - philow) / float(phibins); -int deltaetabins = etabins * 2 - 1; -float deltaetalow = etalow - etaup, deltaetaup = etaup - etalow; -float deltaetabinwidth = (deltaetaup - deltaetalow) / float(deltaetabins); -int deltaphibins = phibins; -float deltaphibinwidth = M_PI * 2 / deltaphibins; -float deltaphilow = 0.0 - deltaphibinwidth / 2.0; -float deltaphiup = M_PI * 2 - deltaphibinwidth / 2.0; - -int tracktype = 1; -int trackonecharge = 1; -int tracktwocharge = -1; -bool processpairs = false; -std::string fTaskConfigurationString = "PendingToConfigure"; - -/// \enum SystemType -/// \brief The type of the system under analysis -enum SystemType { - kNoSystem = 0, ///< no system defined - kpp, ///< **p-p** system - kpPb, ///< **p-Pb** system - kPbp, ///< **Pb-p** system - kPbPb, ///< **Pb-Pb** system - kXeXe, ///< **Xe-Xe** system - knSystems ///< number of handled systems -}; - -/// \enum GeneratorType -/// \brief Which kid of generator data is the task addressing -enum GenType { - kData = 0, ///< actual data, not generated - kMC, ///< Generator level and detector level - kFastMC, ///< Gererator level but stored dataset - kOnTheFly, ///< On the fly generator level data - knGenData ///< number of different generator data types -}; - -/// \enum CentMultEstimatorType -/// \brief The detector used to estimate centrality/multiplicity -enum CentMultEstimatorType { - kNOCM = 0, ///< do not use centrality/multiplicity estimator - kV0M, ///< V0M centrality/multiplicity estimator - kV0A, ///< V0A centrality/multiplicity estimator - kV0C, ///< V0C centrality/multiplicity estimator - kCL0, ///< CL0 centrality/multiplicity estimator - kCL1, ///< CL1 centrality/multiplicity estimator - knCentMultEstimators ///< number of centrality/mutiplicity estimator -}; - -/// \enum MatchRecoGenSpecies -/// \brief The species considered by the matching tast -enum MatchRecoGenSpecies { - kCharged = 0, ///< charged particle/track - kElectron, ///< electron - kMuon, ///< muon - kPion, ///< pion - kKaon, ///< kaon - kProton, ///< proton - kNoOfSpecies, ///< the number of considered species - kWrongSpecies = -1 -}; - -const char* speciesName[kNoOfSpecies] = {"h", "e", "mu", "pi", "ka", "p"}; - -const char* speciesTitle[kNoOfSpecies] = {"", "e", "#mu", "#pi", "K", "p"}; - -namespace filteranalysistask -{ -//============================================================================================ -// The DptDptCorrelationsFilterAnalysisTask output objects -//============================================================================================ -SystemType fSystem = kNoSystem; -GenType fDataType = kData; -CentMultEstimatorType fCentMultEstimator = kV0M; -analysis::CheckRangeCfg traceDCAOutliers; -bool traceOutOfSpeciesParticles = false; -int recoIdMethod = 0; -bool useOwnTrackSelection = false; -TrackSelection ownTrackSelection = getGlobalTrackSelection(); -bool useOwnParticleSelection = false; -bool particleMaxDCAxy = 999.9; -bool particleMaxDCAZ = 999.9; -bool traceCollId0 = false; - -TDatabasePDG* fPDG = nullptr; -TH1F* fhCentMultB = nullptr; -TH1F* fhCentMultA = nullptr; -TH1F* fhVertexZB = nullptr; -TH1F* fhVertexZA = nullptr; -TH1F* fhPB = nullptr; -TH1F* fhPA[kNoOfSpecies] = {nullptr}; -TH1F* fhPtB = nullptr; -TH1F* fhPtA[kNoOfSpecies] = {nullptr}; -TH1F* fhPtPosB = nullptr; -TH1F* fhPtPosA[kNoOfSpecies] = {nullptr}; -TH1F* fhPtNegB = nullptr; -TH1F* fhPtNegA[kNoOfSpecies] = {nullptr}; - -TH1F* fhEtaB = nullptr; -TH1F* fhEtaA = nullptr; - -TH1F* fhPhiB = nullptr; -TH1F* fhPhiA = nullptr; - -TH1F* fhDCAxyB = nullptr; -TH1F* fhDCAxyA = nullptr; -TH1F* fhFineDCAxyA = nullptr; -TH1F* fhDCAzB = nullptr; -TH1F* fhDCAzA = nullptr; -TH1F* fhFineDCAzA = nullptr; - -TH1F* fhTrueCentMultB = nullptr; -TH1F* fhTrueCentMultA = nullptr; -TH1F* fhTrueVertexZB = nullptr; -TH1F* fhTrueVertexZA = nullptr; -TH1F* fhTruePB = nullptr; -TH1F* fhTruePA[kNoOfSpecies] = {nullptr}; -TH1F* fhTruePtB = nullptr; -TH1F* fhTruePtA[kNoOfSpecies] = {nullptr}; -TH1F* fhTruePtPosB = nullptr; -TH1F* fhTruePtPosA[kNoOfSpecies] = {nullptr}; -TH1F* fhTruePtNegB = nullptr; -TH1F* fhTruePtNegA[kNoOfSpecies] = {nullptr}; - -TH1F* fhTrueEtaB = nullptr; -TH1F* fhTrueEtaA = nullptr; - -TH1F* fhTruePhiB = nullptr; -TH1F* fhTruePhiA = nullptr; - -TH1F* fhTrueDCAxyB = nullptr; -TH1F* fhTrueDCAxyA = nullptr; -TH1F* fhTrueDCAzB = nullptr; -TH1F* fhTrueDCAxyBid = nullptr; -TH1F* fhTrueDCAzA = nullptr; -} // namespace filteranalysistask - -namespace filteranalysistaskqa -{ -TH1F* fhTracksOne = nullptr; -TH1F* fhTracksTwo = nullptr; -TH1F* fhTracksOneAndTwo = nullptr; -TH1F* fhTracksNone = nullptr; -TH1F* fhTracksOneUnsel = nullptr; -TH1F* fhTracksTwoUnsel = nullptr; -TH1F* fhTracksOneAndTwoUnsel = nullptr; -TH1F* fhTracksNoneUnsel = nullptr; -TH1F* fhSelectedEvents = nullptr; -} // namespace filteranalysistaskqa - -namespace correlationstask -{ - -/// \enum TrackPairs -/// \brief The track combinations hadled by the class -enum TrackPairs { - kOO = 0, ///< one-one pairs - kOT, ///< one-two pairs - kTO, ///< two-one pairs - kTT, ///< two-two pairs - nTrackPairs ///< the number of track pairs -}; -} // namespace correlationstask - -/// \brief System type according to configuration string -/// \param sysstr The system configuration string -/// \return The internal code for the passed system string -SystemType getSystemType(std::string const& sysstr) -{ - /* we have to figure out how extract the system type */ - if (sysstr.empty() or (sysstr == "PbPb")) { - return kPbPb; - } else if (sysstr == "pp") { - return kpp; - } else if (sysstr == "pPb") { - return kpPb; - } else if (sysstr == "Pbp") { - return kPbp; - } else if (sysstr == "pPb") { - return kpPb; - } else if (sysstr == "XeXe") { - return kXeXe; - } else { - LOGF(fatal, "DptDptCorrelations::getSystemType(). Wrong system type: %d", sysstr.c_str()); - } - return kPbPb; -} - -/// \brief Type of data according to the configuration string -/// \param datastr The data type configuration string -/// \return Internal code for the passed kind of data string -GenType getGenType(std::string const& datastr) -{ - /* we have to figure out how extract the type of data*/ - if (datastr.empty() or (datastr == "data")) { - return kData; - } else if (datastr == "MC") { - return kMC; - } else if (datastr == "FastMC") { - return kFastMC; - } else if (datastr == "OnTheFlyMC") { - return kOnTheFly; - } else { - LOGF(fatal, "DptDptCorrelations::getGenType(). Wrong type of dat: %d", datastr.c_str()); - } - return kData; -} - -template -bool IsEvtSelected(CollisionObject const& collision, float& centormult) -{ - using namespace filteranalysistask; - using namespace dptdptcorrelations; - - bool trigsel = false; - if (fDataType != kData) { - trigsel = true; - } else if (collision.alias()[kINT7]) { - if (collision.sel7()) { - trigsel = true; - } - } - - bool zvtxsel = false; - /* TODO: vertex quality checks */ - if (zvtxlow < collision.posZ() and collision.posZ() < zvtxup) { - zvtxsel = true; - } - - bool centmultsel = false; - switch (fCentMultEstimator) { - case kV0M: - if (collision.centV0M() < 100 and 0 < collision.centV0M()) { - centormult = collision.centV0M(); - centmultsel = true; - } - break; - default: - break; - } - return trigsel and zvtxsel and centmultsel; -} - -template -bool IsEvtSelectedNoCentMult(CollisionObject const& collision, float& centormult) -{ - using namespace filteranalysistask; - using namespace dptdptcorrelations; - - bool trigsel = false; - if (fDataType != kData) { - trigsel = true; - } else if (collision.alias()[kINT7]) { - if (collision.sel7() or collision.sel8()) { - trigsel = true; - } - } - - bool zvtxsel = false; - /* TODO: vertex quality checks */ - if (zvtxlow < collision.posZ() and collision.posZ() < zvtxup) { - zvtxsel = true; - } - - bool centmultsel = false; - switch (fCentMultEstimator) { - case kNOCM: - centormult = 50.0; - centmultsel = true; - break; - default: - break; - } - return trigsel and zvtxsel and centmultsel; -} - -template -bool IsTrueEvtSelected(CollisionObject const& collision, float centormult) -{ - using namespace filteranalysistask; - using namespace dptdptcorrelations; - - bool zvtxsel = false; - /* TODO: vertex quality checks */ - if (zvtxlow < collision.posZ() and collision.posZ() < zvtxup) { - zvtxsel = true; - } - - bool centmultsel = false; - if (centormult < 100 and 0 < centormult) { - centmultsel = true; - } - - return zvtxsel and centmultsel; -} - -template -bool matchTrackType(TrackObject const& track) -{ - using namespace filteranalysistask; - - if (useOwnTrackSelection) { - return ownTrackSelection.IsSelected(track); - } else { - switch (tracktype) { - case 1: - if (track.isGlobalTrack() != 0 || track.isGlobalTrackSDD() != 0) { - return true; - } else { - return false; - } - break; - case 3: /* Run3 track */ - if (track.isGlobalTrack() != 0 || track.isGlobalTrackSDD() != 0) { - return true; - } else { - return false; - } - break; - default: - return false; - } - } -} - -template -inline void AcceptTrack(TrackObject const& track, bool& asone, bool& astwo) -{ - using namespace filteranalysistask; - - asone = false; - astwo = false; - - /* TODO: incorporate a mask in the scanned tracks table for the rejecting track reason */ - if (matchTrackType(track)) { - if (ptlow < track.pt() and track.pt() < ptup and etalow < track.eta() and track.eta() < etaup) { - if (((track.sign() > 0) and (trackonecharge > 0)) or ((track.sign() < 0) and (trackonecharge < 0))) { - asone = true; - } - if (((track.sign() > 0) and (tracktwocharge > 0)) or ((track.sign() < 0) and (tracktwocharge < 0))) { - astwo = true; - } - } - } -} - -template -inline void AcceptParticle(ParticleObject& particle, MCCollisionObject const& collision, bool& asone, bool& astwo) -{ - using namespace filteranalysistask; - - asone = false; - astwo = false; - - float charge = (fPDG->GetParticle(particle.pdgCode())->Charge() / 3 >= 1) ? 1.0 : ((fPDG->GetParticle(particle.pdgCode())->Charge() / 3 <= -1) ? -1.0 : 0.0); - - if (MC::isPhysicalPrimary(particle)) { - if ((particle.mcCollisionId() == 0) and traceCollId0) { - LOGF(info, "Particle %d passed isPhysicalPrimary", particle.globalIndex()); - } - if (useOwnParticleSelection) { - float dcaxy = TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + - (particle.vy() - collision.posY()) * (particle.vy() - collision.posY())); - float dcaz = TMath::Abs(particle.vz() - collision.posZ()); - if (not((dcaxy < particleMaxDCAxy) and (dcaz < particleMaxDCAZ))) { - if ((particle.mcCollisionId() == 0) and traceCollId0) { - auto currparticle = particle; - LOGF(info, "Rejecting particle with dcaxy: %.2f and dcaz: %.2f", dcaxy, dcaz); - LOGF(info, " assigned collision Id: %d, looping on collision Id: %d", currparticle.mcCollisionId(), collision.globalIndex()); - LOGF(info, " Collision x: %.5f, y: %.5f, z: %.5f", collision.posX(), collision.posY(), collision.posZ()); - LOGF(info, " Particle x: %.5f, y: %.5f, z: %.5f", particle.vx(), particle.vy(), particle.vz()); - LOGF(info, " index: %d, pdg code: %d", currparticle.globalIndex(), currparticle.pdgCode()); - while (currparticle.has_mother0()) { - LOGF(info, " mother0 index: %d, mother1 index: %d", currparticle.mother0Id(), currparticle.mother1Id()); - LOGF(info, " Tracking back mother0 index"); - auto newcurrparticle = currparticle.template mother0_as(); - LOGF(info, " assigned collision Id: %d, looping on collision Id: %d", newcurrparticle.mcCollisionId(), collision.globalIndex()); - LOGF(info, " index: %d, pdg code: %d", newcurrparticle.globalIndex(), newcurrparticle.pdgCode()); - LOGF(info, " Passed isPhysicalPrimary(): %s", MC::isPhysicalPrimary(newcurrparticle) ? "YES" : "NO"); - currparticle = newcurrparticle; - } - } - return; - } - } - if (ptlow < particle.pt() and particle.pt() < ptup and etalow < particle.eta() and particle.eta() < etaup) { - if (((charge > 0) and (trackonecharge > 0)) or ((charge < 0) and (trackonecharge < 0))) { - asone = true; - } - if (((charge > 0) and (tracktwocharge > 0)) or ((charge < 0) and (tracktwocharge < 0))) { - astwo = true; - } - } - } else { - if ((particle.mcCollisionId() == 0) and traceCollId0) { - LOGF(info, "Particle %d NOT passed isPhysicalPrimary", particle.globalIndex()); - } - } -} - -template -void fillTrackHistosBeforeSelection(TrackObject const& track) -{ - using namespace filteranalysistask; - - fhPB->Fill(track.p()); - fhPtB->Fill(track.pt()); - fhEtaB->Fill(track.eta()); - fhPhiB->Fill(track.phi()); - if (track.sign() > 0) { - fhPtPosB->Fill(track.pt()); - } else { - fhPtNegB->Fill(track.pt()); - } - fhDCAxyB->Fill(track.dcaXY()); - fhDCAzB->Fill(track.dcaZ()); -} - -template -void fillTrackHistosAfterSelection(TrackObject const& track, MatchRecoGenSpecies sp) -{ - using namespace filteranalysistask; - - /* the charged species should have been called first so avoid double counting */ - if (sp == kCharged) { - fhEtaA->Fill(track.eta()); - fhPhiA->Fill(track.phi()); - fhDCAxyA->Fill(track.dcaXY()); - fhDCAzA->Fill(track.dcaZ()); - if (track.dcaXY() < 1.0) { - fhFineDCAxyA->Fill(track.dcaXY()); - } - if (track.dcaZ() < 1.0) { - fhFineDCAzA->Fill(track.dcaZ()); - } - } - fhPA[sp]->Fill(track.p()); - fhPtA[sp]->Fill(track.pt()); - if (track.sign() > 0) { - fhPtPosA[sp]->Fill(track.pt()); - } else { - fhPtNegA[sp]->Fill(track.pt()); - } -} - -template -void fillParticleHistosBeforeSelection(ParticleObject const& particle, MCCollisionObject const& collision, float charge) -{ - using namespace filteranalysistask; - - fhTruePB->Fill(particle.p()); - fhTruePtB->Fill(particle.pt()); - fhTrueEtaB->Fill(particle.eta()); - fhTruePhiB->Fill(particle.phi()); - if (charge > 0) { - fhTruePtPosB->Fill(particle.pt()); - } else if (charge < 0) { - fhTruePtNegB->Fill(particle.pt()); - } - - float dcaxy = TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + - (particle.vy() - collision.posY()) * (particle.vy() - collision.posY())); - if (traceDCAOutliers.mDoIt and (traceDCAOutliers.mLowValue < dcaxy) and (dcaxy < traceDCAOutliers.mUpValue)) { - fhTrueDCAxyBid->Fill(TString::Format("%d", particle.pdgCode()).Data(), 1.0); - } - - fhTrueDCAxyB->Fill(TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + - (particle.vy() - collision.posY()) * (particle.vy() - collision.posY()))); - fhTrueDCAzB->Fill((particle.vz() - collision.posZ())); -} - -template -void fillParticleHistosAfterSelection(ParticleObject const& particle, MCCollisionObject const& collision, float charge, MatchRecoGenSpecies sp) +namespace o2::analysis::recogenmap { - using namespace filteranalysistask; - - /* the charged species should have been called first so avoid double counting */ - if (sp == kCharged) { - fhTrueEtaA->Fill(particle.eta()); - fhTruePhiA->Fill(particle.phi()); - float dcaxy = TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + - (particle.vy() - collision.posY()) * (particle.vy() - collision.posY())); - if (traceDCAOutliers.mDoIt and (traceDCAOutliers.mLowValue < dcaxy) and (dcaxy < traceDCAOutliers.mUpValue)) { - LOGF(info, "DCAxy outlier: Particle with index %d and pdg code %d assigned to MC collision %d, pT: %f, phi: %f, eta: %f", - particle.globalIndex(), particle.pdgCode(), particle.mcCollisionId(), particle.pt(), particle.phi(), particle.eta()); - LOGF(info, " With status %d and flags %0X", particle.statusCode(), particle.flags()); - } - - fhTrueDCAxyA->Fill(TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) + - (particle.vy() - collision.posY()) * (particle.vy() - collision.posY()))); - fhTrueDCAzA->Fill((particle.vz() - collision.posZ())); - } - fhTruePA[sp]->Fill(particle.p()); - fhTruePtA[sp]->Fill(particle.pt()); - if (charge > 0) { - fhTruePtPosA[sp]->Fill(particle.pt()); - } else { - fhTruePtNegA[sp]->Fill(particle.pt()); - } -} - -template -inline MatchRecoGenSpecies IdentifyTrack(TrackObject const& track) -{ - float nsigmas[kNoOfSpecies]; - if (track.p() < 0.8) { - nsigmas[kCharged] = 999.0f; - nsigmas[kElectron] = track.tpcNSigmaEl(); - nsigmas[kMuon] = track.tpcNSigmaMu(); - nsigmas[kPion] = track.tpcNSigmaPi(); - nsigmas[kKaon] = track.tpcNSigmaKa(); - nsigmas[kProton] = track.tpcNSigmaPr(); - } else { - /* introduce require TOF flag */ - if (track.hasTOF()) { - nsigmas[kCharged] = 999.0f; - nsigmas[kElectron] = sqrtf(track.tpcNSigmaEl() * track.tpcNSigmaEl() + track.tofNSigmaEl() * track.tofNSigmaEl()); - nsigmas[kMuon] = sqrtf(track.tpcNSigmaMu() * track.tpcNSigmaMu() + track.tofNSigmaMu() * track.tofNSigmaMu()); - nsigmas[kPion] = sqrtf(track.tpcNSigmaPi() * track.tpcNSigmaPi() + track.tofNSigmaPi() * track.tofNSigmaPi()); - nsigmas[kKaon] = sqrtf(track.tpcNSigmaKa() * track.tpcNSigmaKa() + track.tofNSigmaKa() * track.tofNSigmaKa()); - nsigmas[kProton] = sqrtf(track.tpcNSigmaPr() * track.tpcNSigmaPr() + track.tofNSigmaPr() * track.tofNSigmaPr()); - } else { - nsigmas[kCharged] = 999.0f; - nsigmas[kElectron] = track.tpcNSigmaEl(); - nsigmas[kMuon] = track.tpcNSigmaMu(); - nsigmas[kPion] = track.tpcNSigmaPi(); - nsigmas[kKaon] = track.tpcNSigmaKa(); - nsigmas[kProton] = track.tpcNSigmaPr(); - } - } - float min_nsigma = 999.0f; - MatchRecoGenSpecies sp_min_nsigma = kWrongSpecies; - for (int sp = 0; sp < kNoOfSpecies; ++sp) { - if (nsigmas[sp] < min_nsigma) { - min_nsigma = nsigmas[sp]; - sp_min_nsigma = MatchRecoGenSpecies(sp); - } - } - bool doublematch = false; - if (min_nsigma < 3.0) { - for (int sp = 0; (sp < kNoOfSpecies) and not doublematch; ++sp) { - if (sp != sp_min_nsigma) { - if (nsigmas[sp] < 3.0) { - doublematch = true; - } - } - } - if (doublematch) { - return kWrongSpecies; - } else { - return sp_min_nsigma; - } - } else { - return kWrongSpecies; - } -} - -template -inline MatchRecoGenSpecies IdentifyParticle(ParticleObject const& particle) -{ - using namespace filteranalysistask; - constexpr int pdgcodeEl = 11L; - constexpr int pdgcodeMu = 13L; - constexpr int pdgcodePi = 211L; - constexpr int pdgcodeKa = 321L; - constexpr int pdgcodePr = 2212L; - - int pdgcode = abs(particle.pdgCode()); - - switch (pdgcode) { - case pdgcodeEl: - return kElectron; - break; - case pdgcodeMu: - return kMuon; - break; - case pdgcodePi: - return kPion; - break; - case pdgcodeKa: - return kKaon; - break; - case pdgcodePr: - return kProton; - break; - - default: - if (traceOutOfSpeciesParticles) { - LOGF(info, "Wrong particle passed selection cuts. PDG code: %d", pdgcode); - } - return kWrongSpecies; - break; - } -} - -} /* end namespace dptdptcorrelations */ - -// Task for correlations analysis -// FIXME: this should really inherit from AnalysisTask but -// we need GCC 7.4+ for that - -using namespace dptdptcorrelations; +std::vector> mclabelpos[2]; +std::vector> mclabelneg[2]; +} // namespace o2::analysis::recogenmap -struct DptDptCorrelationsFilterAnalysisTask { - Configurable cfgTrackType{"trktype", 1, "Type of selected tracks: 0 = no selection, 1 = global tracks FB96, 3 = Run3 tracks. Default 1"}; - Configurable cfgCentMultEstimator{"centmultestimator", "V0M", "Centrality/multiplicity estimator detector: V0M, NOCM: none. Default V0M"}; - Configurable cfgSystem{"syst", "PbPb", "System: pp, PbPb, Pbp, pPb, XeXe. Default PbPb"}; +/// \brief Checks the correspondence generator level <=> detector level +struct CheckGeneratorLevelVsDetectorLevel { + Configurable cfgTrackType{"trktype", 1, "Type of selected tracks: 0 = no selection, 1 = global tracks FB96"}; + Configurable cfgCentMultEstimator{"centmultestimator", "V0M", "Centrality/multiplicity estimator detector: V0M, NOCM: none. Default V0M"}; + Configurable cfgDataType{"datatype", "data", "Data type: data, MC, FastMC, OnTheFlyMC. Default data"}; Configurable cfgBinning{"binning", {28, -7.0, 7.0, 18, 0.2, 2.0, 16, -0.8, 0.8, 72, 0.5}, "triplets - nbins, min, max - for z_vtx, pT, eta and phi, binning plus bin fraction of phi origin shift"}; @@ -731,105 +72,21 @@ struct DptDptCorrelationsFilterAnalysisTask { Configurable cfgRecoIdMethod{"recoidmethod", 0, "Method for identifying reconstructed tracks: 0 PID, 1 mcparticle. Default 0"}; Configurable cfgTrackSelection{"tracksel", {false, false, 0, 70, 0.8, 2.4, 3.2}, "Track selection: {useit: true/false, ongen: true/false, tpccls, tpcxrws, tpcxrfc, dcaxy, dcaz}. Default {false,0.70.0.8,2.4,3.2}"}; Configurable cfgTraceCollId0{"tracecollid0", false, "Trace particles in collisions id 0. Default false"}; + Configurable cfgTrackMultiRec{"trackmultirec", false, "Track muli-reconstructed particles: true, false. Default false"}; + Configurable cfgTrackCollAssoc{"trackcollassoc", false, "Track collision id association, track-mcparticle-mccollision vs. track-collision-mccollision: true, false. Default false"}; - OutputObj fOutput{"MatchingRecoGenGlobalInfo", OutputObjHandlingPolicy::AnalysisObject}; - - Produces acceptedevents; - Produces scannedtracks; - Produces acceptedtrueevents; - Produces scannedtruetracks; - - template - void filterDetectorLevelTracks(TrackListObject const& ftracks, CollisionIndex colix, aod::McParticles const&) - { - using namespace filteranalysistask; - - int acceptedtracks = 0; - - for (auto& track : ftracks) { - bool asone = false; - bool astwo = false; - if (not(track.mcParticleId() < 0)) { - /* correctly reconstructed track */ - /* before track selection */ - fillTrackHistosBeforeSelection(track); - - /* track selection */ - /* tricky because the boolean columns issue */ - AcceptTrack(track, asone, astwo); - if (asone or astwo) { - /* the track has been accepted */ - /* fill the charged tracks histograms */ - fillTrackHistosAfterSelection(track, kCharged); - /* let's identify it */ - MatchRecoGenSpecies sp = kWrongSpecies; - if (cfgRecoIdMethod == 0) { - sp = IdentifyTrack(track); - } else if (cfgRecoIdMethod == 1) { - sp = IdentifyParticle(track.mcParticle()); - } - if (sp != kWrongSpecies) { - /* fill the species histograms */ - fillTrackHistosAfterSelection(track, sp); - } - acceptedtracks++; - scannedtracks(colix, (uint8_t)asone, (uint8_t)astwo, track.pt(), track.eta(), track.phi()); - LOGF(MATCHRECGENLOGTRACKS, "Accepted track with global Id %d and with assigned collision Id %d", track.globalIndex(), track.collisionId()); - } - } - } - LOGF(MATCHRECGENLOGCOLLISIONS, "Accepted %d reconstructed tracks", acceptedtracks); - } - - template - void filterParticles(ParticleListObject const& particles, MCCollisionObject const& mccollision, CollisionIndex colix) - { - using namespace filteranalysistask; - - int acceptedparticles = 0; - - for (auto& particle : particles) { - float charge = 0.0; - TParticlePDG* pdgparticle = fPDG->GetParticle(particle.pdgCode()); - if (pdgparticle != nullptr) { - charge = (pdgparticle->Charge() / 3 >= 1) ? 1.0 : ((pdgparticle->Charge() / 3 <= -1) ? -1.0 : 0.0); - } - - bool asone = false; - bool astwo = false; - if (charge != 0) { - /* before particle selection */ - fillParticleHistosBeforeSelection(particle, mccollision, charge); - - /* track selection */ - /* tricky because the boolean columns issue */ - AcceptParticle(particle, mccollision, asone, astwo); - if (asone or astwo) { - /* the track has been accepted */ - /* fill the charged particle histograms */ - fillParticleHistosAfterSelection(particle, mccollision, charge, kCharged); - /* let's identify the particle */ - MatchRecoGenSpecies sp = IdentifyParticle(particle); - if (sp != kWrongSpecies) { - fillParticleHistosAfterSelection(particle, mccollision, charge, sp); - } - acceptedparticles++; - scannedtruetracks(colix, (uint8_t)asone, (uint8_t)astwo, particle.pt(), particle.eta(), particle.phi()); - } - } else { - if ((particle.mcCollisionId() == 0) and traceCollId0) { - LOGF(info, "Particle %d with fractional charge or equal to zero", particle.globalIndex()); - } - } - } - LOGF(MATCHRECGENLOGCOLLISIONS, "Accepted %d generated particles", acceptedparticles); - } + HistogramRegistry histos{"RecoGenHistograms", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + typedef enum { kBEFORE = 0, + kAFTER } beforeafterselection; + typedef enum { kPOSITIVE = 0, + kNEGATIVE } colllabelsign; + enum { kMATCH = 0, + kDONTMATCH }; - void init(InitContext const&) + void init(InitContext const& context) { - using namespace filteranalysistask; - - LOGF(info, "FilterAnalysisTask::init()"); + using namespace o2::analysis::recogenmap; + using namespace o2::analysis::dptdptfilter; /* update with the configurable values */ /* the binning */ @@ -845,13 +102,7 @@ struct DptDptCorrelationsFilterAnalysisTask { /* the track types and combinations */ tracktype = cfgTrackType.value; /* the centrality/multiplicity estimation */ - if (cfgCentMultEstimator->compare("V0M") == 0) { - fCentMultEstimator = kV0M; - } else if (cfgCentMultEstimator->compare("NOCM") == 0) { - fCentMultEstimator = kNOCM; - } else { - LOGF(fatal, "Centrality/Multiplicity estimator %s not supported yet", cfgCentMultEstimator->c_str()); - } + fCentMultEstimator = getCentMultEstimator(cfgCentMultEstimator); traceDCAOutliers = cfgTraceDCAOutliers; traceOutOfSpeciesParticles = cfgTraceOutOfSpeciesParticles; recoIdMethod = cfgRecoIdMethod; @@ -886,488 +137,11 @@ struct DptDptCorrelationsFilterAnalysisTask { } traceCollId0 = cfgTraceCollId0; - /* if the system type is not known at this time, we have to put the initalization somewhere else */ - fSystem = getSystemType(cfgSystem); - fDataType = kMC; - fPDG = TDatabasePDG::Instance(); - - /* create the output list which will own the task histograms */ - TList* fOutputList = new TList(); - fOutputList->SetOwner(true); - fOutput.setObject(fOutputList); - - /* incorporate configuration parameters to the output */ - fOutputList->Add(new TParameter("TrackType", cfgTrackType, 'f')); - fOutputList->Add(new TParameter("TrackOneCharge", trackonecharge, 'f')); - fOutputList->Add(new TParameter("TrackTwoCharge", tracktwocharge, 'f')); - - if ((fDataType == kData) or (fDataType == kMC)) { - /* create the reconstructed data histograms */ - if (fSystem > kPbp) { - fhCentMultB = new TH1F("CentralityB", "Centrality before cut; centrality (%)", 100, 0, 100); - fhCentMultA = new TH1F("CentralityA", "Centrality; centrality (%)", 100, 0, 100); - } else { - /* for pp, pPb and Pbp systems use multiplicity instead */ - fhCentMultB = new TH1F("MultiplicityB", "Multiplicity before cut; multiplicity (%)", 100, 0, 100); - fhCentMultA = new TH1F("MultiplicityA", "Multiplicity; multiplicity (%)", 100, 0, 100); - } - - fhVertexZB = new TH1F("VertexZB", "Vertex Z; z_{vtx}", 60, -15, 15); - fhVertexZA = new TH1F("VertexZA", "Vertex Z; z_{vtx}", zvtxbins, zvtxlow, zvtxup); - - fhPB = new TH1F("fHistPB", "p distribution for reconstructed before;p (GeV/c);dN/dp (c/GeV)", 100, 0.0, 15.0); - fhPtB = new TH1F("fHistPtB", "p_{T} distribution for reconstructed before;p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhPtPosB = new TH1F("fHistPtPosB", "P_{T} distribution for reconstructed (#plus) before;P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhPtNegB = new TH1F("fHistPtNegB", "P_{T} distribution for reconstructed (#minus) before;P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhEtaB = new TH1F("fHistEtaB", "#eta distribution for reconstructed before;#eta;counts", 40, -2.0, 2.0); - fhEtaA = new TH1F("fHistEtaA", "#eta distribution for reconstructed;#eta;counts", etabins, etalow, etaup); - fhPhiB = new TH1F("fHistPhiB", "#phi distribution for reconstructed before;#phi;counts", 360, 0.0, 2 * M_PI); - fhPhiA = new TH1F("fHistPhiA", "#phi distribution for reconstructed;#phi;counts", 360, 0.0, 2 * M_PI); - fhDCAxyB = new TH1F("DCAxyB", "DCA_{xy} distribution for reconstructed before;DCA_{xy} (cm);counts", 1000, -4.0, 4.0); - fhDCAxyA = new TH1F("DCAxyA", "DCA_{xy} distribution for reconstructed;DCA_{xy} (cm);counts", 1000, -4., 4.0); - fhFineDCAxyA = new TH1F("FineDCAxyA", "DCA_{xy} distribution for reconstructed;DCA_{xy} (cm);counts", 4000, -1.0, 1.0); - fhDCAzB = new TH1F("DCAzB", "DCA_{z} distribution for reconstructed before;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - fhDCAzA = new TH1F("DCAzA", "DCA_{z} distribution for reconstructed;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - fhFineDCAzA = new TH1F("FineDCAzA", "DCA_{z} distribution for reconstructed;DCA_{z} (cm);counts", 4000, -1.0, 1.0); - - for (int sp = 0; sp < kNoOfSpecies; ++sp) { - fhPA[sp] = new TH1F(TString::Format("fHistPA_%s", speciesName[sp]).Data(), - TString::Format("p distribution for reconstructed %s;p (GeV/c);dN/dp (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - fhPtA[sp] = new TH1F(TString::Format("fHistPtA_%s", speciesName[sp]), - TString::Format("p_{T} distribution for reconstructed %s;p_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - fhPtPosA[sp] = new TH1F(TString::Format("fHistPtPosA_%s", speciesName[sp]), - TString::Format("P_{T} distribution for reconstructed %s^{#plus};P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - fhPtNegA[sp] = new TH1F(TString::Format("fHistPtNegA_%s", speciesName[sp]), - TString::Format("P_{T} distribution for reconstructed %s^{#minus};P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - } - - /* add the hstograms to the output list */ - fOutputList->Add(fhCentMultB); - fOutputList->Add(fhCentMultA); - fOutputList->Add(fhVertexZB); - fOutputList->Add(fhVertexZA); - fOutputList->Add(fhPB); - fOutputList->Add(fhPtB); - fOutputList->Add(fhPtPosB); - fOutputList->Add(fhPtNegB); - fOutputList->Add(fhEtaB); - fOutputList->Add(fhEtaA); - fOutputList->Add(fhPhiB); - fOutputList->Add(fhPhiA); - fOutputList->Add(fhDCAxyB); - fOutputList->Add(fhDCAxyA); - fOutputList->Add(fhFineDCAxyA); - fOutputList->Add(fhDCAzB); - fOutputList->Add(fhDCAzA); - fOutputList->Add(fhFineDCAzA); - - for (int sp = 0; sp < kNoOfSpecies; ++sp) { - fOutputList->Add(fhPA[sp]); - fOutputList->Add(fhPtA[sp]); - fOutputList->Add(fhPtPosA[sp]); - fOutputList->Add(fhPtNegA[sp]); - } - } - - if (fDataType != kData) { - /* create the true data histograms */ - if (fSystem > kPbp) { - fhTrueCentMultB = new TH1F("TrueCentralityB", "Centrality before (truth); centrality (%)", 100, 0, 100); - fhTrueCentMultA = new TH1F("TrueCentralityA", "Centrality (truth); centrality (%)", 100, 0, 100); - } else { - /* for pp, pPb and Pbp systems use multiplicity instead */ - fhTrueCentMultB = new TH1F("TrueMultiplicityB", "Multiplicity before (truth); multiplicity (%)", 100, 0, 100); - fhTrueCentMultA = new TH1F("TrueMultiplicityA", "Multiplicity (truth); multiplicity (%)", 100, 0, 100); - } - - fhTrueVertexZB = new TH1F("TrueVertexZB", "Vertex Z before (truth); z_{vtx}", 60, -15, 15); - fhTrueVertexZA = new TH1F("TrueVertexZA", "Vertex Z (truth); z_{vtx}", zvtxbins, zvtxlow, zvtxup); - - fhTruePB = new TH1F("fTrueHistPB", "p distribution before (truth);p (GeV/c);dN/dp (c/GeV)", 100, 0.0, 15.0); - fhTruePtB = new TH1F("fTrueHistPtB", "p_{T} distribution before (truth);p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhTruePtPosB = new TH1F("fTrueHistPtPosB", "P_{T} distribution (#plus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhTruePtNegB = new TH1F("fTrueHistPtNegB", "P_{T} distribution (#minus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0); - fhTrueEtaB = new TH1F("fTrueHistEtaB", "#eta distribution before (truth);#eta;counts", 40, -2.0, 2.0); - fhTrueEtaA = new TH1F("fTrueHistEtaA", "#eta distribution (truth);#eta;counts", etabins, etalow, etaup); - fhTruePhiB = new TH1F("fTrueHistPhiB", "#phi distribution before (truth);#phi;counts", 360, 0.0, 2 * M_PI); - fhTruePhiA = new TH1F("fTrueHistPhiA", "#phi distribution (truth);#phi;counts", 360, 0.0, 2 * M_PI); - fhTrueDCAxyB = new TH1F("TrueDCAxyB", "DCA_{xy} distribution for generated before;DCA_{xy} (cm);counts", 1000, -4.0, 4.0); - if (traceDCAOutliers.mDoIt) { - fhTrueDCAxyBid = new TH1F("PDGCodeDCAxyB", - TString::Format("PDG code within %.2f<|DCA_{#it{xy}}|<%.2f; PDG code", traceDCAOutliers.mLowValue, traceDCAOutliers.mUpValue).Data(), - 100, 0.5, 100.5); - } - fhTrueDCAxyA = new TH1F("TrueDCAxyA", "DCA_{xy} distribution for generated;DCA_{xy};counts (cm)", 1000, -4., 4.0); - fhTrueDCAzB = new TH1F("TrueDCAzB", "DCA_{z} distribution for generated before;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - fhTrueDCAzA = new TH1F("TrueDCAzA", "DCA_{z} distribution for generated;DCA_{z} (cm);counts", 1000, -4.0, 4.0); - - for (int sp = 0; sp < kNoOfSpecies; ++sp) { - fhTruePA[sp] = new TH1F(TString::Format("fTrueHistPA_%s", speciesName[sp]).Data(), - TString::Format("p distribution %s (truth);p (GeV/c);dN/dp (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - fhTruePtA[sp] = new TH1F(TString::Format("fTrueHistPtA_%s", speciesName[sp]), - TString::Format("p_{T} distribution %s (truth);p_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - fhTruePtPosA[sp] = new TH1F(TString::Format("fTrueHistPtPosA_%s", speciesName[sp]), - TString::Format("P_{T} distribution %s^{#plus} (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - fhTruePtNegA[sp] = new TH1F(TString::Format("fTrueHistPtNegA_%s", speciesName[sp]), - TString::Format("P_{T} distribution %s^{#minus} (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", speciesTitle[sp]).Data(), - ptbins, ptlow, ptup); - } - - /* add the hstograms to the output list */ - fOutputList->Add(fhTrueCentMultB); - fOutputList->Add(fhTrueCentMultA); - fOutputList->Add(fhTrueVertexZB); - fOutputList->Add(fhTrueVertexZA); - fOutputList->Add(fhTruePB); - fOutputList->Add(fhTruePtB); - fOutputList->Add(fhTruePtPosB); - fOutputList->Add(fhTruePtNegB); - fOutputList->Add(fhTrueEtaB); - fOutputList->Add(fhTrueEtaA); - fOutputList->Add(fhTruePhiB); - fOutputList->Add(fhTruePhiA); - fOutputList->Add(fhTrueDCAxyB); - if (traceDCAOutliers.mDoIt) { - fOutputList->Add(fhTrueDCAxyBid); - } - fOutputList->Add(fhTrueDCAxyA); - fOutputList->Add(fhTrueDCAzB); - fOutputList->Add(fhTrueDCAzA); - - for (int sp = 0; sp < kNoOfSpecies; ++sp) { - fOutputList->Add(fhTruePA[sp]); - fOutputList->Add(fhTruePtA[sp]); - fOutputList->Add(fhTruePtPosA[sp]); - fOutputList->Add(fhTruePtNegA[sp]); - } - } - } - - using FullTracksPID = soa::Join; - - void processWithCentDetectorLevel(aod::CollisionEvSelCent const& collision, FullTracksPID const& ftracks, aod::McParticles const& particles) - { - using namespace filteranalysistask; - - LOGF(MATCHRECGENLOGCOLLISIONS, "FilterAnalysisTask::processWithCentDetectorLevel(). New collision with %d tracks", ftracks.size()); - - fhCentMultB->Fill(collision.centV0M()); - fhVertexZB->Fill(collision.posZ()); - bool acceptedevent = false; - float centormult = -100.0; - if (IsEvtSelected(collision, centormult)) { - acceptedevent = true; - fhCentMultA->Fill(collision.centV0M()); - fhVertexZA->Fill(collision.posZ()); - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - - filterDetectorLevelTracks(ftracks, acceptedevents.lastIndex(), particles); - } else { - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& track : ftracks) { - scannedtracks(acceptedevents.lastIndex(), (uint8_t) false, (uint8_t) false, track.pt(), track.eta(), track.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithCentDetectorLevel, "Process MC detector level with centrality", false); - - void processWithoutCentDetectorLevel(aod::CollisionEvSel const& collision, FullTracksPID const& ftracks, aod::McParticles const& particles) - { - using namespace filteranalysistask; - - LOGF(MATCHRECGENLOGCOLLISIONS, "FilterAnalysisTask::processWithoutCentDetectorLevel(). New collision with collision id %d and with %d tracks", collision.bcId(), ftracks.size()); - - /* the task does not have access to either centrality nor multiplicity - classes information, so it has to live without it. - For the time being we assign a value of 50% */ - fhCentMultB->Fill(50.0); - fhVertexZB->Fill(collision.posZ()); - bool acceptedevent = false; - float centormult = -100.0; - if (IsEvtSelectedNoCentMult(collision, centormult)) { - acceptedevent = true; - fhCentMultA->Fill(50.0); - fhVertexZA->Fill(collision.posZ()); - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - - LOGF(MATCHRECGENLOGCOLLISIONS, "Accepted collision with BC id %d and collision Id %d", collision.bcId(), collision.globalIndex()); - filterDetectorLevelTracks(ftracks, acceptedevents.lastIndex(), particles); - } else { - acceptedevents(collision.bcId(), collision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& track : ftracks) { - scannedtracks(acceptedevents.lastIndex(), (uint8_t) false, (uint8_t) false, track.pt(), track.eta(), track.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithoutCentDetectorLevel, "Process MC detector level without centrality", false); - - void processWithCentGeneratorLevel(aod::McCollision const& mccollision, - soa::Join const& collisions, - aod::McParticles const& mcparticles) - { - using namespace filteranalysistask; - - LOGF(MATCHRECGENLOGCOLLISIONS, "FilterAnalysisTask::processWithCentGeneratorLevel(). New generated collision %d reconstructed collisions and %d particles", collisions.size(), mcparticles.size()); - - /* TODO: in here we have to decide what to do in the following cases - - On the fly production -> clearly we will need a different process - - reconstructed collisions without generated associated -> we need a different process or a different signature - - multiplicity/centrality classes extracted from the reconstructed collision but then - - generated collision without associated reconstructed collision: how to extract mutliplicity/centrality classes? - - generated collision with several associated reconstructed collisions: from which to extract multiplicity/centrality classes? - */ - if (collisions.size() > 1) { - LOGF(error, "FilterAnalysisTask::processWithCentGeneratorLevel(). Generated collision with more than one reconstructed collisions. Processing only the first for centrality/multiplicity classes extraction"); - } - - for (auto& collision : collisions) { - float cent = collision.centV0M(); - fhTrueCentMultB->Fill(cent); - fhTrueVertexZB->Fill(mccollision.posZ()); - - bool acceptedevent = false; - if (IsTrueEvtSelected(mccollision, cent)) { - acceptedevent = true; - fhTrueCentMultA->Fill(cent); - fhTrueVertexZA->Fill(mccollision.posZ()); - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, cent); - - filterParticles(mcparticles, collision, acceptedtrueevents.lastIndex()); - } else { - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, cent); - for (auto& particle : mcparticles) { - scannedtruetracks(acceptedtrueevents.lastIndex(), (uint8_t) false, (uint8_t) false, particle.pt(), particle.eta(), particle.phi()); - } - } - break; /* TODO: only processing the first reconstructed collision for centrality/multiplicity class estimation */ - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithCentGeneratorLevel, "Process generated with centrality", false); - - void processWithoutCentGeneratorLevel(aod::McCollision const& mccollision, - aod::McParticles const& mcparticles) - { - using namespace filteranalysistask; - - LOGF(MATCHRECGENLOGCOLLISIONS, "FilterAnalysisTask::processWithoutCentGeneratorLevel(). New generated collision with %d particles", mcparticles.size()); - - /* the task does not have access to either centrality nor multiplicity - classes information, so it has to live without it. - For the time being we assign a value of 50% */ - fhTrueCentMultB->Fill(50.0); - fhTrueVertexZB->Fill(mccollision.posZ()); - - bool acceptedevent = false; - float centormult = 50.0; - if (IsTrueEvtSelected(mccollision, centormult)) { - acceptedevent = true; - fhTrueCentMultA->Fill(centormult); - fhTrueVertexZA->Fill(mccollision.posZ()); - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, centormult); - - filterParticles(mcparticles, mccollision, acceptedtrueevents.lastIndex()); - } else { - acceptedtrueevents(mccollision.bcId(), mccollision.posZ(), (uint8_t)acceptedevent, centormult); - for (auto& particle : mcparticles) { - scannedtruetracks(acceptedtrueevents.lastIndex(), (uint8_t) false, (uint8_t) false, particle.pt(), particle.eta(), particle.phi()); - } - } - } - PROCESS_SWITCH(DptDptCorrelationsFilterAnalysisTask, processWithoutCentGeneratorLevel, "Process generated without centrality", false); -}; - -// Checking the filtered tables -/* it seems we cannot use a base class task */ -// struct TracksAndEventClassificationQABase { - -void initQATask(InitContext const&, TList* outlst) -{ - using namespace filteranalysistaskqa; - - fhTracksOne = new TH1F("TracksOne", "Tracks as track one;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksTwo = new TH1F("TracksTwo", "Tracks as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksOneAndTwo = new TH1F("TracksOneAndTwo", "Tracks as track one and as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksNone = new TH1F("TracksNone", "Not selected tracks;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksOneUnsel = new TH1F("TracksOneUnsel", "Tracks as track one;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksTwoUnsel = new TH1F("TracksTwoUnsel", "Tracks as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksOneAndTwoUnsel = new TH1F("TracksOneAndTwoUnsel", "Tracks as track one and as track two;number of tracks;events", 1500, 0.0, 1500.0); - fhTracksNoneUnsel = new TH1F("TracksNoneUnsel", "Not selected tracks;number of tracks;events", 1500, 0.0, 1500.0); - fhSelectedEvents = new TH1F("SelectedEvents", "Selected events;;events", 2, 0.0, 2.0); - fhSelectedEvents->GetXaxis()->SetBinLabel(1, "Not selected events"); - fhSelectedEvents->GetXaxis()->SetBinLabel(2, "Selected events"); - - outlst->Add(fhTracksOne); - outlst->Add(fhTracksTwo); - outlst->Add(fhTracksOneAndTwo); - outlst->Add(fhTracksNone); - outlst->Add(fhTracksOneUnsel); - outlst->Add(fhTracksTwoUnsel); - outlst->Add(fhTracksOneAndTwoUnsel); - outlst->Add(fhTracksNoneUnsel); - outlst->Add(fhSelectedEvents); -} - -template -void processQATask(FilteredCollision const& collision, - FilteredTracks const& tracks) -{ - using namespace filteranalysistaskqa; - - if (collision.eventaccepted() != (uint8_t) true) { - fhSelectedEvents->Fill(0.5); - } else { - fhSelectedEvents->Fill(1.5); - } - - int ntracks_one = 0; - int ntracks_two = 0; - int ntracks_one_and_two = 0; - int ntracks_none = 0; - for (auto& track : tracks) { - if ((track.trackacceptedasone() != (uint8_t) true) and (track.trackacceptedastwo() != (uint8_t) true)) { - ntracks_none++; - } - if ((track.trackacceptedasone() == (uint8_t) true) and (track.trackacceptedastwo() == (uint8_t) true)) { - ntracks_one_and_two++; - } - if (track.trackacceptedasone() == (uint8_t) true) { - ntracks_one++; - } - if (track.trackacceptedastwo() == (uint8_t) true) { - ntracks_two++; - } - } - if (collision.eventaccepted() != (uint8_t) true) { - /* control for non selected events */ - fhTracksOneUnsel->Fill(ntracks_one); - fhTracksTwoUnsel->Fill(ntracks_two); - fhTracksNoneUnsel->Fill(ntracks_none); - fhTracksOneAndTwoUnsel->Fill(ntracks_one_and_two); - } else { - fhTracksOne->Fill(ntracks_one); - fhTracksTwo->Fill(ntracks_two); - fhTracksNone->Fill(ntracks_none); - fhTracksOneAndTwo->Fill(ntracks_one_and_two); - } -} -// }; - -/* it seems we cannot use a base class task */ -// struct TracksAndEventClassificationQARec : TracksAndEventClassificationQABase { -struct TracksAndEventClassificationQARec { - OutputObj fOutput{"FliterTaskRecoQA", OutputObjHandlingPolicy::AnalysisObject}; - - void init(InitContext const& context) - { - TList* fOutputList = new TList(); - fOutputList->SetName("FilterTaskRecoQA"); - fOutputList->SetOwner(true); - fOutput.setObject(fOutputList); - - initQATask(context, fOutputList); - } - - Filter onlyacceptedevents = (aod::dptdptcorrelations::eventaccepted == (uint8_t) true); - Filter onlyacceptedtracks = ((aod::dptdptcorrelations::trackacceptedasone == (uint8_t) true) or (aod::dptdptcorrelations::trackacceptedastwo == (uint8_t) true)); - - void process(soa::Filtered::iterator const& collision, soa::Filtered const& tracks) - { - LOGF(MATCHRECGENLOGCOLLISIONS, "New filtered collision with BC id %d and with %d accepted tracks", collision.bcId(), tracks.size()); - processQATask(collision, tracks); - } -}; - -/* it seems we cannot use a base class task */ -// struct TracksAndEventClassificationQAGen : TracksAndEventClassificationQABase { -struct TracksAndEventClassificationQAGen { - OutputObj fOutput{"FliterTaskGenQA", OutputObjHandlingPolicy::AnalysisObject}; - - void init(InitContext const& context) - { - TList* fOutputList = new TList(); - fOutputList->SetName("FilterTaskGenQA"); - fOutputList->SetOwner(true); - fOutput.setObject(fOutputList); - - initQATask(context, fOutputList); - } - - Filter onlyacceptedevents = (aod::dptdptcorrelations::eventaccepted == (uint8_t) true); - Filter onlyacceptedtracks = ((aod::dptdptcorrelations::trackacceptedasone == (uint8_t) true) or (aod::dptdptcorrelations::trackacceptedastwo == (uint8_t) true)); - - void process(soa::Filtered::iterator const& collision, soa::Filtered const& tracks) - { - LOGF(MATCHRECGENLOGCOLLISIONS, "New filtered generated collision with BC id %d and with %d accepted tracks", collision.bcId(), tracks.size()); - processQATask(collision, tracks); - } -}; - -namespace recogenmap -{ -std::vector> mclabelpos[2]; -std::vector> mclabelneg[2]; -} // namespace recogenmap - -/// \brief Checks the correspondence generator level <=> detector level -struct CheckGeneratorLevelVsDetectorLevel { - Configurable cfgTrackType{"trktype", 1, "Type of selected tracks: 0 = no selection, 1 = global tracks FB96"}; - Configurable cfgCentMultEstimator{"centmultestimator", "V0M", "Centrality/multiplicity estimator detector: V0M, NOCM: none. Default V0M"}; - Configurable cfgDataType{"datatype", "data", "Data type: data, MC, FastMC, OnTheFlyMC. Default data"}; - Configurable cfgBinning{"binning", - {28, -7.0, 7.0, 18, 0.2, 2.0, 16, -0.8, 0.8, 72, 0.5}, - "triplets - nbins, min, max - for z_vtx, pT, eta and phi, binning plus bin fraction of phi origin shift"}; - Configurable cfgTrackMultiRec{"trackmultirec", false, "Track muli-reconstructed particles: true, false. Default false"}; - Configurable cfgTrackCollAssoc{"trackcollassoc", false, "Track collision id association, track-mcparticle-mccollision vs. track-collision-mccollision: true, false. Default false"}; - - HistogramRegistry histos{"RecoGenHistograms", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; - TDatabasePDG* fPDG = nullptr; - typedef enum { kBEFORE = 0, - kAFTER } beforeafterselection; - typedef enum { kPOSITIVE = 0, - kNEGATIVE } colllabelsign; - enum { kMATCH = 0, - kDONTMATCH }; - - void init(InitContext const& context) - { - using namespace dptdptcorrelations; - using namespace dptdptcorrelations::filteranalysistask; - - /* update with the configurable values */ - /* the binning */ - ptbins = cfgBinning->mPTbins; - ptlow = cfgBinning->mPTmin; - ptup = cfgBinning->mPTmax; - etabins = cfgBinning->mEtabins; - etalow = cfgBinning->mEtamin; - etaup = cfgBinning->mEtamax; - zvtxbins = cfgBinning->mZVtxbins; - zvtxlow = cfgBinning->mZVtxmin; - zvtxup = cfgBinning->mZVtxmax; - /* the track types and combinations */ - tracktype = cfgTrackType.value; - /* the centrality/multiplicity estimation */ - if (cfgCentMultEstimator->compare("V0M") == 0) { - fCentMultEstimator = kV0M; - } else if (cfgCentMultEstimator->compare("NOCM") == 0) { - fCentMultEstimator = kNOCM; - } else { - LOGF(fatal, "Centrality/Multiplicity estimator %s not supported yet", cfgCentMultEstimator->c_str()); - } - fDataType = getGenType(cfgDataType); - - constexpr float TWOPI = 2.0f * static_cast(M_PI); + fDataType = getDataType(cfgDataType); fPDG = TDatabasePDG::Instance(); AxisSpec deltaEta = {100, -2, 2, "#Delta#eta"}; - AxisSpec deltaPhi = {100, 0, TWOPI, "#Delta#varphi (rad)"}; + AxisSpec deltaPhi = {100, 0, constants::math::TwoPI, "#Delta#varphi (rad)"}; AxisSpec deltaPt = {1000, 0, 4, "#Delta#it{p}_{T} (GeV/#it{c})"}; AxisSpec mrectimes = {11, -0.5f, 10.5f, "##/particle"}; AxisSpec detectors = {32, -0.5, 31.5, "Detectors"}; @@ -1382,15 +156,15 @@ struct CheckGeneratorLevelVsDetectorLevel { histos.add("before/positivecolid/mrDeltaPt", "#Delta#it{p}_{T} multirec tracks", kTH1F, {deltaPt}); histos.add("before/positivecolid/multirec", "Multiple reconstruction", kTH1F, {mrectimes}); histos.add("before/positivecolid/genrecoeta", "#eta Generated vs reconstructed", kTH2F, {{100, -1.0, 1.0, "#eta reco"}, {100, -1.0, 1.0, "#eta gen"}}); - histos.add("before/positivecolid/genrecophi", "#varphi Generated vs reconstructed", kTH2F, {{100, 0, TWOPI, "#varphi (rad) reco"}, {100, 0, TWOPI, "#varphi (rad) gen"}}); + histos.add("before/positivecolid/genrecophi", "#varphi Generated vs reconstructed", kTH2F, {{100, 0, constants::math::TwoPI, "#varphi (rad) reco"}, {100, 0, constants::math::TwoPI, "#varphi (rad) gen"}}); histos.add("before/positivecolid/genrecopt", "#it{p}_{T} Generated vs reconstructed", kTH2F, {{1000, 0, 10.0, "#it{p}_{T} (GeV/#it{c}) reco"}, {1000, 0, 10.0, "#it{p}_{T} (GeV/#it{c}) gen"}}); histos.add("before/positivecolid/detectormap", "Active detectors", kTH1F, {detectors}); histos.add("before/positivecolid/matchcollid", "particle MC coll Id <=> track coll MC coll Id", kTH1F, {{2, 0.0, 2.0}}); histos.add("before/positivecolid/genrecomreta", "#eta Generated vs reconstructed (mr)", kTH2F, {{100, -1.0, 1.0, "#eta reco"}, {100, -1.0, 1.0, "#eta gen"}}); - histos.add("before/positivecolid/genrecomrphi", "#varphi Generated vs reconstructed (mr)", kTH2F, {{100, 0, TWOPI, "#varphi (rad) reco"}, {100, 0, TWOPI, "#varphi (rad) gen"}}); + histos.add("before/positivecolid/genrecomrphi", "#varphi Generated vs reconstructed (mr)", kTH2F, {{100, 0, constants::math::TwoPI, "#varphi (rad) reco"}, {100, 0, constants::math::TwoPI, "#varphi (rad) gen"}}); histos.add("before/positivecolid/genrecomrpt", "#it{p}_{T} Generated vs reconstructed (mr)", kTH2F, {{1000, 0, 10.0, "#it{p}_{T} (GeV/#it{c}) reco"}, {1000, 0, 10.0, "#it{p}_{T} (GeV/#it{c}) gen"}}); histos.add("before/positivecolid/recomreta", "#eta Reconstructed (mr)", kTH1F, {{100, -1.0, 1.0, "#eta"}}); - histos.add("before/positivecolid/recomrphi", "#varphi Reconstructed (mr)", kTH1F, {{100, 0, TWOPI, "#varphi (rad)"}}); + histos.add("before/positivecolid/recomrphi", "#varphi Reconstructed (mr)", kTH1F, {{100, 0, constants::math::TwoPI, "#varphi (rad)"}}); histos.add("before/positivecolid/recomrpt", "#it{p}_{T} Reconstructed (mr)", kTH1F, {{1000, 0, 10.0, "#it{p}_{T} (GeV/#it{c})"}}); histos.add("before/positivecolid/detectormapmr", "Active detectors (mr)", kTH1F, {detectors}); histos.add("before/positivecolid/matchcollidmr", "particle MC coll Id <=> track coll MC coll Id (mr)", kTH1F, {{2, 0.0, 2.0}}); @@ -1421,11 +195,10 @@ struct CheckGeneratorLevelVsDetectorLevel { template void collectData(TracskListObject const& tracks, ParticlesListObject const& mcParticles, CollisionsListObject const& colls) { - using namespace recogenmap; + using namespace o2::analysis::recogenmap; static constexpr std::string_view dir[] = {"before/", "after/"}; static constexpr std::string_view colldir[] = {"positivecolid/", "negativecolid/"}; - constexpr float TWOPI = 2.0F * static_cast(M_PI); int nrec_poslabel = 0; int nrec_neglabel = 0; @@ -1482,10 +255,10 @@ struct CheckGeneratorLevelVsDetectorLevel { float deltaeta = track1.eta() - track2.eta(); float deltaphi = track1.phi() - track2.phi(); if (deltaphi < 0) { - deltaphi += TWOPI; + deltaphi += constants::math::TwoPI; } - if (deltaphi > TWOPI) { - deltaphi -= TWOPI; + if (deltaphi > constants::math::TwoPI) { + deltaphi -= constants::math::TwoPI; } float deltapt = (track1.pt() > track2.pt()) ? track1.pt() - track2.pt() : track2.pt() - track1.pt(); @@ -1552,12 +325,11 @@ struct CheckGeneratorLevelVsDetectorLevel { } } - void - processMapChecksBeforeCuts(soa::Join const& tracks, - soa::Join const& collisions, - aod::McParticles const& mcParticles) + template + void processMapChecksBeforeCuts(TracksObject const& tracks, CollisionsObject const& collisions, aod::McParticles const& mcParticles) { - using namespace recogenmap; + using namespace o2::analysis::recogenmap; + using namespace o2::analysis::dptdptfilter; for (int i = 0; i < 2; ++i) { mclabelpos[i].clear(); @@ -1606,68 +378,12 @@ struct CheckGeneratorLevelVsDetectorLevel { collectData(tracks, mcParticles, collisions); collectData(tracks, mcParticles, collisions); } - PROCESS_SWITCH(CheckGeneratorLevelVsDetectorLevel, processMapChecksBeforeCuts, "Process detector <=> generator levels mapping checks before selection cuts", false); - void processMapChecksCutsWithCent(soa::Join const& tracks, - soa::Join const& collisions, - aod::McParticles const& mcParticles) + template + void processMapChecksAfterCuts(TracksObject const& tracks, CollisionsObject const& collisions, aod::McParticles const& mcParticles) { - using namespace recogenmap; - - for (int i = 0; i < 2; ++i) { - mclabelpos[i].clear(); - mclabelneg[i].clear(); - mclabelpos[i].resize(mcParticles.size()); - mclabelneg[i].resize(mcParticles.size()); - } - - size_t nreco = 0; - size_t ngen = 0; - - for (auto& part : mcParticles) { - auto pdgpart = fPDG->GetParticle(part.pdgCode()); - if (pdgpart != nullptr) { - float charge = (pdgpart->Charge() >= 3) ? 1.0 : ((pdgpart->Charge() <= -3) ? -1.0 : 0.0); - if (charge != 0.0) { - ngen++; - } - } - } - - // Let's go through the reco-gen mapping to detect multi-reconstructed particles - for (auto& track : tracks) { - int64_t recix = track.globalIndex(); - int32_t label = track.mcParticleId(); - if (not(track.collisionId() < 0)) { - auto coll = collisions.iteratorAt(track.collisionId()); - float centormult = -100.0f; - if (IsEvtSelected(coll, centormult)) { - if (not(label < 0)) { - bool asone = false; - bool astwo = false; - - AcceptTrack(track, asone, astwo); - if (asone or astwo) { - /* the track has been accepted */ - nreco++; - LOGF(MATCHRECGENLOGTRACKS, "Accepted track with global Id %d and collision Id %d has label %d associated to MC collision %d", recix, track.collisionId(), label, track.mcParticle().mcCollisionId()); - mclabelpos[kPOSITIVE][label].push_back(recix); - } - } - } - } - } - LOGF(info, "New dataframe (DF) with %d generated charged particles and %d reconstructed accepted tracks", ngen, nreco); - - collectData(tracks, mcParticles, collisions); - } - PROCESS_SWITCH(CheckGeneratorLevelVsDetectorLevel, processMapChecksCutsWithCent, "Process detector <=> generator levels mapping checks after selection cuts", false); - - void processMapChecksCutsWithoutCent(soa::Join const& tracks, - soa::Join const& collisions, - aod::McParticles const& mcParticles) - { - using namespace recogenmap; + using namespace o2::analysis::recogenmap; + using namespace o2::analysis::dptdptfilter; for (int i = 0; i < 2; ++i) { mclabelpos[i].clear(); @@ -1695,9 +411,9 @@ struct CheckGeneratorLevelVsDetectorLevel { int32_t label = track.mcParticleId(); if (not(label < 0)) { if (not(track.collisionId() < 0)) { - auto coll = collisions.iteratorAt(track.collisionId()); + typename CollisionsObject::iterator coll = collisions.iteratorAt(track.collisionId()); float centormult = -100.0f; - if (IsEvtSelectedNoCentMult(coll, centormult)) { + if (IsEvtSelected(coll, centormult)) { bool asone = false; bool astwo = false; @@ -1716,21 +432,28 @@ struct CheckGeneratorLevelVsDetectorLevel { collectData(tracks, mcParticles, collisions); } - PROCESS_SWITCH(CheckGeneratorLevelVsDetectorLevel, processMapChecksCutsWithoutCent, "Process detector <=> generator levels mapping checks after selection cuts", false); - void processDummy(aod::Collisions const& colls) + void processMapChecksWithCent(soa::Join const& tracks, + soa::Join const& collisions, + aod::McParticles const& mcParticles) + { + processMapChecksBeforeCuts(tracks, collisions, mcParticles); + processMapChecksAfterCuts(tracks, collisions, mcParticles); + } + PROCESS_SWITCH(CheckGeneratorLevelVsDetectorLevel, processMapChecksWithCent, "Process detector <=> generator levels with centrality/multiplicity information", false); + + void processMapChecksWithoutCent(soa::Join const& tracks, + soa::Join const& collisions, + aod::McParticles const& mcParticles) { - LOGF(MATCHRECGENLOGCOLLISIONS, "Got a new set of %d collisions", colls.size()); + processMapChecksBeforeCuts(tracks, collisions, mcParticles); + processMapChecksAfterCuts(tracks, collisions, mcParticles); } - PROCESS_SWITCH(CheckGeneratorLevelVsDetectorLevel, processDummy, "Dummy process of detector <=> generator levels mapping checks", true); + PROCESS_SWITCH(CheckGeneratorLevelVsDetectorLevel, processMapChecksWithoutCent, "Process detector <=> generator levels without centrality/multiplicity information", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - WorkflowSpec workflow{ - adaptAnalysisTask(cfgc, SetDefaultProcesses{{{"processWithoutCent", true}, {"processWithoutCentMC", true}}}), - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc, SetDefaultProcesses{{{"processMapChecksCutsWithCent", false}}})}; + WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; return workflow; }