diff --git a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt index 96e49a62303..a7f4334d614 100644 --- a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt +++ b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt @@ -18,3 +18,8 @@ o2physics_add_dpl_workflow(gammaconversionstruthonlymc SOURCES gammaConversionsTruthOnlyMc.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(emc-pi0-qc + SOURCES emcalPi0QC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/PhotonMeson/Tasks/emcalPi0QC.cxx b/PWGEM/PhotonMeson/Tasks/emcalPi0QC.cxx new file mode 100644 index 00000000000..a190a071963 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/emcalPi0QC.cxx @@ -0,0 +1,385 @@ +// 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 +#include +#include +#include +#include +#include +#include +#include + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" +#include "Framework/HistogramRegistry.h" + +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Centrality.h" + +#include "EMCALBase/Geometry.h" +#include "PWGJE/DataModel/EMCALClusters.h" +#include "DataFormatsEMCAL/Cell.h" +#include "DataFormatsEMCAL/Constants.h" +#include "DataFormatsEMCAL/AnalysisCluster.h" + +#include "CommonDataFormat/InteractionRecord.h" + +#include "TLorentzVector.h" +#include "TVector3.h" + +// \struct Pi0QCTask +/// \brief Simple monitoring task for EMCal clusters +/// \author Joshua Koenig , Goethe University Frankfurt +/// \since 25.05.2022 +/// +/// This task is meant to be used for QC for the emcal using properties of the pi0 +/// - ... +/// Simple event selection using the flag doEventSel is provided, which selects INT7 events if set to 1 +/// For pilot beam data, instead of relying on the event selection, one can veto specific BC IDS using the flag +/// fDoVetoBCID. + +using namespace o2::framework; +using namespace o2::framework::expressions; +using collisionEvSelIt = o2::soa::Join::iterator; +using selectedClusters = o2::soa::Filtered; +using selectedCluster = o2::soa::Filtered; +using selectedAmbiguousClusters = o2::soa::Filtered; +using selectedAmbiguousCluster = o2::soa::Filtered; + +struct Photon { + Photon(float eta_tmp, float phi_tmp, float energy_tmp, int clusid = 0) + { + eta = eta_tmp; + phi = phi_tmp; + energy = energy_tmp; + theta = 2 * std::atan2(std::exp(-eta), 1); + px = energy * std::sin(theta) * std::cos(phi); + py = energy * std::sin(theta) * std::sin(phi); + pz = energy * std::cos(theta); + pt = std::sqrt(px * px + py * py); + photon.SetPxPyPzE(px, py, pz, energy); + id = clusid; + } + + TLorentzVector photon; + float pt; + float px; + float py; + float pz; + float eta; + float phi; + float energy; + float theta; + int id; +}; + +struct Meson { + Meson(Photon p1, Photon p2) : pgamma1(p1), + pgamma2(p2) + { + pMeson = p1.photon + p2.photon; + } + Photon pgamma1; + Photon pgamma2; + TLorentzVector pMeson; + + float getMass() const { return pMeson.M(); }; + float getPt() const { return pMeson.Pt(); }; + float getOpeningAngle() const { return pgamma1.photon.Angle(pgamma2.photon.Vect()); }; +}; + +struct Pi0QCTask { + HistogramRegistry mHistManager{"NeutralMesonHistograms"}; + o2::emcal::Geometry* mGeometry = nullptr; + + // configurable parameters + // TODO adapt mDoEventSel switch to also allow selection of other triggers (e.g. EMC7) + Configurable mDoEventSel{"doEventSel", 0, "demand kINT7"}; + Configurable mVetoBCID{"vetoBCID", "", "BC ID(s) to be excluded, this should be used as an alternative to the event selection"}; + Configurable mSelectBCID{"selectBCID", "all", "BC ID(s) to be included, this should be used as an alternative to the event selection"}; + Configurable mVertexCut{"vertexCut", -1, "apply z-vertex cut with value in cm"}; + Configurable mTimeMin{"TimeMinCut", -600, "apply min timing cut (in ns)"}; + Configurable mTimeMax{"TimeMaxCut", 900, "apply min timing cut (in ns)"}; + Configurable mClusterMinM02Cut{"MinM02Cut", 0.1, "apply min M02 cut"}; + Configurable mClusterMaxM02Cut{"MaxM02Cut", 0.7, "apply max M02 cut"}; + Configurable mMinEnergyCut{"MinEnergyCut", 0.7, "apply min cluster energy cut"}; + Configurable mMinNCellsCut{"MinNCellsCut", 1, "apply min cluster number of cell cut"}; + Configurable mClusterDefinition{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"}; + std::vector mVetoBCIDs; + std::vector mSelectBCIDs; + + // define cluster filter. It selects only those clusters which are of the type + // specified in the string mClusterDefinition,e.g. kV3Default, which is V3 clusterizer with default + // clusterization parameters + o2::aod::EMCALClusterDefinition clusDef = o2::aod::emcalcluster::getClusterDefinitionFromString(mClusterDefinition.value); + Filter clusterDefinitionSelection = o2::aod::emcalcluster::definition == static_cast(clusDef); + + // define container for photons + std::vector mPhotons; + + /// \brief Create output histograms and initialize geometry + void init(InitContext const&) + { + // create histograms + using o2HistType = HistType; + using o2Axis = AxisSpec; + + // load geometry just in case we need it + mGeometry = o2::emcal::Geometry::GetInstanceFromRunNumber(300000); + + // create common axes + LOG(info) << "Creating histograms"; + const o2Axis bcAxis{3501, -0.5, 3500.5}; + const o2Axis energyAxis{makeClusterBinning(), "#it{p}_{T} (GeV)"}; + + // event properties + mHistManager.add("eventsAll", "Number of events", o2HistType::kTH1F, {{1, 0.5, 1.5}}); + mHistManager.add("eventsSelected", "Number of events", o2HistType::kTH1F, {{1, 0.5, 1.5}}); + mHistManager.add("eventBCAll", "Bunch crossing ID of event (all events)", o2HistType::kTH1F, {bcAxis}); + mHistManager.add("eventBCSelected", "Bunch crossing ID of event (selected events)", o2HistType::kTH1F, {bcAxis}); + mHistManager.add("eventVertexZAll", "z-vertex of event (all events)", o2HistType::kTH1F, {{200, -20, 20}}); + mHistManager.add("eventVertexZSelected", "z-vertex of event (selected events)", o2HistType::kTH1F, {{200, -20, 20}}); + + // cluster properties + mHistManager.add("clusterE", "Energy of cluster", o2HistType::kTH1F, {energyAxis}); + mHistManager.add("clusterE_SimpleBinning", "Energy of cluster", o2HistType::kTH1F, {{400, 0, 100}}); + mHistManager.add("clusterEtaPhi", "Eta and phi of cluster", o2HistType::kTH2F, {{100, -1, 1}, {100, 0, 2 * TMath::Pi()}}); + mHistManager.add("clusterM02", "M02 of cluster", o2HistType::kTH1F, {{400, 0, 5}}); + mHistManager.add("clusterM20", "M20 of cluster", o2HistType::kTH1F, {{400, 0, 2.5}}); + mHistManager.add("clusterNLM", "Number of local maxima of cluster", o2HistType::kTH1I, {{10, 0, 10}}); + mHistManager.add("clusterNCells", "Number of cells in cluster", o2HistType::kTH1I, {{50, 0, 50}}); + mHistManager.add("clusterDistanceToBadChannel", "Distance to bad channel", o2HistType::kTH1F, {{100, 0, 100}}); + + // meson related histograms + mHistManager.add("invMassVsPt", "invariant mass and pT of meson candidates", o2HistType::kTH2F, {{400, 0, 0.8}, {energyAxis}}); + mHistManager.add("invMassVsPtBackground", "invariant mass and pT of background meson candidates", o2HistType::kTH2F, {{400, 0, 0.8}, {energyAxis}}); + + if (mVetoBCID->length()) { + std::stringstream parser(mVetoBCID.value); + std::string token; + int bcid; + while (std::getline(parser, token, ',')) { + bcid = std::stoi(token); + LOG(info) << "Veto BCID " << bcid; + mVetoBCIDs.push_back(bcid); + } + } + if (mSelectBCID.value != "all") { + std::stringstream parser(mSelectBCID.value); + std::string token; + int bcid; + while (std::getline(parser, token, ',')) { + bcid = std::stoi(token); + LOG(info) << "Select BCID " << bcid; + mSelectBCIDs.push_back(bcid); + } + } + } + /// \brief Process EMCAL clusters that are matched to a collisions + void processCollisions(collisionEvSelIt const& theCollision, selectedClusters const& clusters, o2::aod::BCs const& bcs) + { + mHistManager.fill(HIST("eventsAll"), 1); + + // do event selection if mDoEventSel is specified + // currently the event selection is hard coded to kINT7 + // but other selections are possible that are defined in TriggerAliases.h + if (mDoEventSel && (!theCollision.alias()[kINT7])) { + LOG(debug) << "Event not selected becaus it is not kINT7, skipping"; + return; + } + mHistManager.fill(HIST("eventVertexZAll"), theCollision.posZ()); + if (mVertexCut > 0 && TMath::Abs(theCollision.posZ()) > mVertexCut) { + LOG(debug) << "Event not selected because of z-vertex cut z= " << theCollision.posZ() << " > " << mVertexCut << " cm, skipping"; + return; + } + mHistManager.fill(HIST("eventsSelected"), 1); + mHistManager.fill(HIST("eventVertexZSelected"), theCollision.posZ()); + + ProcessClusters(theCollision, clusters, bcs); + ProcessMesons(theCollision, clusters, bcs); + } + PROCESS_SWITCH(Pi0QCTask, processCollisions, "Process clusters from collision", false); + + /// \brief Process EMCAL clusters that are not matched to a collision + /// This is not needed for most users + void processAmbiguous(o2::aod::BC const bc, selectedAmbiguousClusters const& clusters) + { + // loop over bc , if requested (mVetoBCID >= 0), reject everything from a certain BC + // this can be used as alternative to event selection (e.g. for pilot beam data) + // TODO: remove this loop and put it in separate process function that only takes care of ambiguous clusters + o2::InteractionRecord eventIR; + eventIR.setFromLong(bc.globalBC()); + mHistManager.fill(HIST("eventBCAll"), eventIR.bc); + if (std::find(mVetoBCIDs.begin(), mVetoBCIDs.end(), eventIR.bc) != mVetoBCIDs.end()) { + LOG(info) << "Event rejected because of veto BCID " << eventIR.bc; + return; + } + if (mSelectBCIDs.size() && (std::find(mSelectBCIDs.begin(), mSelectBCIDs.end(), eventIR.bc) == mSelectBCIDs.end())) { + return; + } + mHistManager.fill(HIST("eventBCSelected"), eventIR.bc); + + // ToDo: Add mode if collision is not found + // ProcessClusters(theCollision, clusters, bcs); + // ProcessMesons(theCollision, clusters, bcs); + } + PROCESS_SWITCH(Pi0QCTask, processAmbiguous, "Process Ambiguous clusters", false); + + /// \brief Process EMCAL clusters that are matched to a collisions + template + void ProcessClusters(collisionEvSelIt const& theCollision, Clusters const& clusters, o2::aod::BCs const& bcs) + { + // clear photon vector + mPhotons.clear(); + + // loop over all clusters from accepted collision + // auto eventClusters = clusters.select(o2::aod::emcalcluster::bcId == theCollision.bc().globalBC()); + for (const auto& cluster : clusters) { + // fill histograms of cluster properties + // in this implementation the cluster properties are directly + // loaded from the flat table, in the future one should + // consider using the AnalysisCluster object to work with + // after loading. + LOG(debug) << "Cluster energy: " << cluster.energy(); + LOG(debug) << "Cluster time: " << cluster.time(); + LOG(debug) << "Cluster M02: " << cluster.m02(); + mHistManager.fill(HIST("clusterE"), cluster.energy()); + mHistManager.fill(HIST("clusterE_SimpleBinning"), cluster.energy()); + mHistManager.fill(HIST("clusterEtaPhi"), cluster.eta(), cluster.phi()); + mHistManager.fill(HIST("clusterM02"), cluster.m02()); + mHistManager.fill(HIST("clusterM20"), cluster.m20()); + mHistManager.fill(HIST("clusterNLM"), cluster.nlm()); + mHistManager.fill(HIST("clusterNCells"), cluster.nCells()); + mHistManager.fill(HIST("clusterDistanceToBadChannel"), cluster.distanceToBadChannel()); + + // apply basic cluster cuts + if (cluster.energy() < mMinEnergyCut) { + LOG(debug) << "Cluster rejected because of energy cut"; + continue; + } + if (cluster.nCells() <= mMinNCellsCut) { + LOG(debug) << "Cluster rejected because of nCells cut"; + continue; + } + if (cluster.m02() < mClusterMinM02Cut || cluster.m02() > mClusterMaxM02Cut) { + LOG(debug) << "Cluster rejected because of m02 cut"; + continue; + } + if (cluster.time() < mTimeMin || cluster.time() > mTimeMax) { + LOG(debug) << "Cluster rejected because of time cut"; + continue; + } + + // put clusters in photon vector + // ToDo: At the moment, the eta and phi values are not corrected for a shift of the primary vertex! Should only be a small effect but has to be corrected + mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy(), cluster.id())); + } + } + + /// \brief Process meson candidates, caluclate invariant mass and pT and fill histograms + template + void ProcessMesons(collisionEvSelIt const& theCollision, Clusters const& clusters, o2::aod::BCs const& bcs) + { + // if less then 2 clusters are found, skip event + if (mPhotons.size() < 2) { + return; + } + + // loop over all photon combinations and build meson candidates + for (unsigned int ig1 = 0; ig1 < mPhotons.size(); ++ig1) { + for (unsigned int ig2 = ig1 + 1; ig2 < mPhotons.size(); ++ig2) { + + // build meson from photons + Meson meson(mPhotons[ig1], mPhotons[ig2]); + mHistManager.fill(HIST("invMassVsPt"), meson.getMass(), meson.getPt()); + + // calculate background candidates (rotation background) + CalculateBackground(meson, ig1, ig2); + } + } + } + + /// \brief Calculate background (using rotation background method) + void CalculateBackground(const Meson& meson, unsigned int ig1, unsigned int ig2) + { + // if less than 3 clusters are present, skip event + if (mPhotons.size() < 3) { + return; + } + const double rotationAngle = M_PI / 2.0; //0.78539816339; // rotaion angle 90° + + TLorentzVector lvRotationPhoton1; // photon candidates which get rotated + TLorentzVector lvRotationPhoton2; // photon candidates which get rotated + TVector3 lvRotationPion; // rotation axis + for (unsigned int ig3 = 0; ig3 < mPhotons.size(); ++ig3) { + // continue if photons are identical + if (ig3 == ig1 || ig3 == ig2) { + continue; + } + // calculate rotation axis + lvRotationPion = (meson.pMeson).Vect(); + + // initialize photons for rotation + lvRotationPhoton1.SetPxPyPzE(mPhotons[ig1].px, mPhotons[ig1].py, mPhotons[ig1].pz, mPhotons[ig1].energy); + lvRotationPhoton2.SetPxPyPzE(mPhotons[ig2].px, mPhotons[ig2].py, mPhotons[ig2].pz, mPhotons[ig2].energy); + + // rotate photons around rotation axis + lvRotationPhoton1.Rotate(rotationAngle, lvRotationPion); + lvRotationPhoton2.Rotate(rotationAngle, lvRotationPion); + + // initialize Photon objects for rotated photons + Photon rotPhoton1(lvRotationPhoton1.Eta(), lvRotationPhoton1.Phi(), lvRotationPhoton1.E(), mPhotons[ig1].id); + Photon rotPhoton2(lvRotationPhoton2.Eta(), lvRotationPhoton2.Phi(), lvRotationPhoton2.E(), mPhotons[ig2].id); + + // build meson from rotated photons + Meson mesonRotated1(rotPhoton1, mPhotons[ig3]); + Meson mesonRotated2(rotPhoton2, mPhotons[ig3]); + + // Fill histograms + mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated1.getMass(), mesonRotated1.getPt()); + mHistManager.fill(HIST("invMassVsPtBackground"), mesonRotated2.getMass(), mesonRotated2.getPt()); + } + } + + /// \brief Create binning for cluster energy/pT axis (variable bin size) + /// direct port from binning often used in AliPhysics for debugging + /// \return vector with bin limits + std::vector makeClusterBinning() const + { + + std::vector result; + int nBinsPt = 179; + double maxPt = 60; + for (Int_t i = 0; i < nBinsPt + 1; i++) { + if (i < 100) { + result.emplace_back(0.10 * i); + } else if (i < 140) { + result.emplace_back(10. + 0.25 * (i - 100)); + } else if (i < 180) { + result.emplace_back(20. + 1.00 * (i - 140)); + } else { + result.emplace_back(maxPt); + } + } + return result; + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, TaskName{"EMCPi0QCTask"}, SetDefaultProcesses{{{"processCollisions", true}, {"processAmbiguous", false}}}), + adaptAnalysisTask(cfgc, TaskName{"EMCPi0QCTaskAmbiguous"}, SetDefaultProcesses{{{"processCollisions", false}, {"processAmbiguous", true}}})}; + return workflow; +} \ No newline at end of file