diff --git a/Detectors/GlobalTracking/CMakeLists.txt b/Detectors/GlobalTracking/CMakeLists.txt index 197f7730a3ea4..34bfe394032b7 100644 --- a/Detectors/GlobalTracking/CMakeLists.txt +++ b/Detectors/GlobalTracking/CMakeLists.txt @@ -13,6 +13,7 @@ o2_add_library(GlobalTracking TARGETVARNAME targetName SOURCES src/MatchTPCITS.cxx src/MatchTOF.cxx + src/MatchHMP.cxx src/MatchTPCITSParams.cxx src/MatchCosmics.cxx src/MatchCosmicsParams.cxx @@ -28,16 +29,18 @@ o2_add_library(GlobalTracking O2::DataFormatsFT0 O2::DataFormatsTOF O2::DataFormatsTRD + O2::DataFormatsHMP O2::ITSReconstruction O2::FT0Reconstruction O2::TPCFastTransformation O2::GPUO2Interface O2::TPCBase O2::TPCReconstruction - O2::TPCCalibration O2::TOFBase O2::TOFCalibration - O2::TOFWorkflowUtils + O2::HMPIDBase + #O2::HMPIDCalibration + O2::HMPIDWorkflow O2::SimConfig O2::DataFormatsFT0 O2::DataFormatsGlobalTracking @@ -54,6 +57,7 @@ o2_target_root_dictionary(GlobalTracking include/GlobalTracking/MatchGlobalFwdParam.h include/GlobalTracking/MatchGlobalFwdAssessment.h include/GlobalTracking/MatchTOF.h + include/GlobalTracking/MatchHMP.h include/GlobalTracking/MatchCosmics.h include/GlobalTracking/MatchCosmicsParams.h include/GlobalTracking/MatchITSTPCQC.h diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h new file mode 100644 index 0000000000000..18eff074c0a8a --- /dev/null +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h @@ -0,0 +1,315 @@ +// 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 ALICEO2_GLOBTRACKING_MATCHHMP_ +#define ALICEO2_GLOBTRACKING_MATCHHMP_ + +#include +#include +#include +#include +#include +#include +#include "ReconstructionDataFormats/Track.h" +#include "ReconstructionDataFormats/TrackTPCITS.h" +#include "ReconstructionDataFormats/TrackTPCTOF.h" +#include "ReconstructionDataFormats/MatchInfoTOFReco.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "CommonDataFormat/EvIndex.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "CommonUtils/TreeStreamRedirector.h" + +#include "SimulationDataFormat/MCTruthContainer.h" +#include "DetectorsBase/Propagator.h" +#include "MathUtils/Cartesian.h" +#include "MathUtils/Utils.h" +#include "CommonConstants/MathConstants.h" +#include "CommonConstants/PhysicsConstants.h" +#include "CommonConstants/GeomConstants.h" +#include "DetectorsBase/GeometryManager.h" + +#include "DataFormatsHMP/Cluster.h" +#include "GlobalTracking/MatchTPCITS.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsTRD/TrackTRD.h" +#include "ReconstructionDataFormats/PID.h" +#include "TPCFastTransform.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "ReconstructionDataFormats/MatchInfoHMP.h" + +#include "HMPIDBase/Geo.h" +#include "DataFormatsHMP/Cluster.h" + +namespace o2 +{ + +namespace globaltracking +{ +class RecoContainer; +} + +namespace dataformats +{ +template +class MCTruthContainer; +} + +namespace globaltracking +{ +/* enum from TOF::cluster +enum { kUpLeft = 0, // 2^0, 1st bit + kUp = 1, // 2^1, 2nd bit + kUpRight = 2, // 2^2, 3rd bit + kRight = 3, // 2^3, 4th bit + kDownRight = 4, // 2^4, 5th bit + kDown = 5, // 2^5, 6th bit + kDownLeft = 6, // 2^6, 7th bit + kLeft = 7 }; // 2^7, 8th bit */ + +/* ef: I can not get it to compile w/o commenting out this struct +///< original track in the currently loaded TPC-ITS reco output +struct TrackLocTPCITS : public o2::track::TrackParCov { + o2::dataformats::EvIndex source; ///< track origin id + o2::math_utils::Bracket timeBins; ///< bracketing time-bins + float zMin = 0; // min possible Z of this track + float zMax = 0; // max possible Z of this track + int matchID = MinusOne; ///< entry (none if MinusOne) of TOF matchTOF struct in the mMatchesTOF + TrackLocTPCITS(const o2::track::TrackParCov& src, int tch, int tid) : o2::track::TrackParCov(src), source(tch, tid) {} + TrackLocTPCITS() = default; + ClassDefNV(TrackLocTPCITS, 1); // RS TODO: is this class needed? +};// */ + +class MatchHMP +{ + + using Geo = o2::hmpid::Geo; + using Cluster = o2::hmpid::Cluster; // ef: not hmp + using evGIdx = o2::dataformats::EvIndex; + using evIdx = o2::dataformats::EvIndex; + using timeEst = o2::dataformats::TimeStampWithError; + using matchTrack = std::pair; + using trkType = o2::dataformats::MatchInfoTOFReco::TrackType; + + public: + // MatchHMP() = default; + // ~MatchHMP() = default; + ///< perform matching for provided input + void run(const o2::globaltracking::RecoContainer& inp); + + ///< print settings + void print() const; + void printCandidatesHMP() const; + + ///< set time tolerance on track-HMP times comparison + void setTimeTolerance(float val) { mTimeTolerance = val; } + ///< get tolerance on track-HMP times comparison + float getTimeTolerance() const { return mTimeTolerance; } + + ///< set space tolerance on track-HMP times comparison + void setSpaceTolerance(float val) { mSpaceTolerance = val; } + ///< get tolerance on track-HMP times comparison + float getSpaceTolerance() const { return mSpaceTolerance; } + + ///< set number of sigma used to do the matching + void setSigmaTimeCut(float val) { mSigmaTimeCut = val; } + ///< get number of sigma used to do the matching + float getSigmaTimeCut() const { return mSigmaTimeCut; } + + enum DebugFlagTypes : UInt_t { + MatchTreeAll = 0x1 << 1, // ///< produce matching candidates tree for all candidates + }; + ///< check if partucular flags are set + bool isDebugFlag(UInt_t flags) const { return mDBGFlags & flags; } + + ///< get debug trees flags + UInt_t getDebugFlags() const { return mDBGFlags; } + + ///< set or unset debug stream flag + void setDebugFlag(UInt_t flag, bool on = true); + + std::vector& getMatchedTrackVector(trkType index) + { + return mMatchedTracks[index]; + } + + std::vector& getMatchedHMPLabelsVector(trkType index) { return mOutHMPLabels[index]; } ///< get vector of HMP label of matched tracks + + ///< set input TPC tracks cluster indices + void setTPCTrackClusIdxInp(const gsl::span inp) + { + mTPCTrackClusIdx = inp; + } + + ///< set input TPC cluster sharing map + void setTPCClustersSharingMap(const gsl::span inp) + { + mTPCRefitterShMap = inp; + } + + ///< set input TPC clusters + void setTPCClustersInp(const o2::tpc::ClusterNativeAccess* inp) + { + mTPCClusterIdxStruct = inp; + } + + void setFIT(bool value = true) { mIsFIT = value; } + int findFITIndex(int bc); + + // bool makeConstrainedTPCTrack(int matchedID, o2::dataformats::TrackTPCHMP& trConstr); has not been declared + bool makeConstrainedTPCTrack(int matchedID, o2::dataformats::TrackTPCTOF& trConstr); + + ///< populate externally provided container by TOF-time-constrained TPC tracks + template + + /* void makeConstrainedTPCTracks(V& container) + { + //checkRefitter(); + int nmatched = mMatchedTracks[trkType::TPC].size(), nconstrained = 0; + container.resize(nmatched); + for (unsigned i = 0; i < nmatched; i++) { + if (makeConstrainedTPCTrack(i, container[nconstrained])) { + nconstrained++; + } + } + container.resize(nconstrained); + } + */ + + void setTS(unsigned long creationTime) + { + mTimestamp = creationTime; + } + unsigned long getTS() const { return mTimestamp; } + + private: + // bool prepareFITData(); + int prepareInteractionTimes(); + bool prepareTracks(); + bool prepareHMPClusters(); + void doMatching(int sec); + + bool propagateToRefX(o2::track::TrackParCov& trc, float xRef /*in cm*/, float stepInCm /*in cm*/, o2::track::TrackLTIntegral& intLT); + bool propagateToRefXWithoutCov(o2::track::TrackParCov& trc, float xRef /*in cm*/, float stepInCm /*in cm*/, float bz); + + void updateTimeDependentParams(); + + /* + void addTPCSeed(const o2::tpc::TrackTPC& _tr, o2::dataformats::GlobalTrackID srcGID); + void addITSTPCSeed(const o2::dataformats::TrackTPCITS& _tr, o2::dataformats::GlobalTrackID srcGID); + void addTRDSeed(const o2::trd::TrackTRD& _tr, o2::dataformats::GlobalTrackID srcGID, float time0, float terr); + void addConstrainedSeed(o2::track::TrackParCov& trc, o2::dataformats::GlobalTrackID srcGID, o2::track::TrackLTIntegral intLT0, timeEst timeMUS); + void addTPCTRDSeed(const o2::track::TrackParCov& _tr, o2::dataformats::GlobalTrackID srcGID, int tpcID); + void addITSTPCTRDSeed(const o2::track::TrackParCov& _tr, o2::dataformats::GlobalTrackID srcGID, int tpcID); + */ + + //================================================================ + + // Data members + const o2::globaltracking::RecoContainer* mRecoCont = nullptr; + + o2::InteractionRecord mStartIR{0, 0}; ///< IR corresponding to the start of the TF + + // for derived class + int mCurrTracksTreeEntry = 0; ///< current tracks tree entry loaded to memory + + float mXRef = 500.; ///< reference radius to propage tracks for matching + + bool mMCTruthON = false; ///< flag availability of MC truth + + ///========== Parameters to be set externally, e.g. from CCDB ==================== + float mBz = 0; ///< nominal Bz + float mMaxInvPt = 999.; ///< derived from nominal Bz + + // to be done later + float mTPCTBinMUS = 0.; ///< TPC time bin duration in microseconds + float mTPCTBinMUSInv = 0.; ///< inverse TPC time bin duration in microseconds + float mTPCBin2Z = 0.; ///< conversion coeff from TPC time-bin to Z + + bool mIsCosmics = false; ///< switch on to reconstruct cosmics and match with TPC + float mTimeTolerance = 1e3; ///< tolerance in ns for track-TOF time bracket matching + float mSpaceTolerance = 10; ///< tolerance in cm for track-TOF time bracket matching + int mSigmaTimeCut = 30.; ///< number of sigmas to cut on time when matching the track to the TOF cluster + + bool mIsFIT = false; + bool mIsITSTPCused = false; + bool mIsTPCused = false; + bool mIsTPCTRDused = false; + bool mIsITSTPCTRDused = false; + bool mSetHighPurity = false; + + unsigned long mTimestamp = 0; ///< in ms + + // from ruben + gsl::span mTPCTracksArray; ///< input TPC tracks span + + ///>>>------ these are input arrays which should not be modified by the matching code + // since this info is provided by external device + std::vector mITSTPCTracksArrayInp; ///< input tracks + gsl::span mHMPClustersArrayInp; ///< input HMPID clusters + + /// data needed for refit of time-constrained TPC tracks + gsl::span mTPCTrackClusIdx; ///< input TPC track cluster indices span + gsl::span mTPCRefitterShMap; ///< externally set TPC clusters sharing map + const o2::tpc::ClusterNativeAccess* mTPCClusterIdxStruct = nullptr; ///< struct holding the TPC cluster indices + std::unique_ptr mTPCTransform; ///< TPC cluster transformation + std::unique_ptr mTPCRefitter; ///< TPC refitter used for TPC tracks refit during the reconstruction + + const o2::dataformats::MCTruthContainer* mHMPClusLabels; ///< input HMP clusters MC labels (pointer to read from tree) + + int mNotPropagatedToHMP[trkType::SIZE]; ///< number of tracks failing in propagation + + ///< working copy of the input tracks + std::vector mTracksWork[trkType::SIZE]; ///< track params prepared for matching + time value + std::vector mTracksLblWork[trkType::SIZE]; ///< TPCITS track labels + std::vector mLTinfos[trkType::SIZE]; ///< expected times and others + std::vector mTrackGid[trkType::SIZE]; ///< expected times and others + ///< per sector indices of track entry in mTracksWork + std::array, o2::constants::math::NSectors> mTracksSectIndexCache[trkType::SIZE]; + + std::vector mExtraTPCFwdTime; ///< track extra params for TPC tracks: Fws Max time + std::vector mHMPClusWork; ///< track params prepared for matching + + std::vector mSideTPC; ///< track side for TPC tracks + + ///< per sector indices of HMP cluster entry in mHMPClusWork + std::array, o2::constants::math::NSectors> mHMPClusSectIndexCache; + + ///< array of track-HMPCluster pairs from the matching + std::vector mMatchedTracksPairs; + + ///< array of matched HMPCluster with matching information (residuals, expected times...) with the corresponding vector of indices + std::vector mMatchedTracks[trkType::SIZEALL]; // this is the output of the matching -> UNCONS, CONSTR + std::vector mOutHMPLabels[trkType::SIZEALL]; ///< TOF label of matched tracks + + std::vector mMatchedTracksIndex[trkType::SIZE]; // vector of indexes of the tracks to be matched + + // ef : move to cxx + // int mNumOfClusters; // number of clusters to be matched + // std::unique_ptr mMatchedClustersIndex; //[mNumOfClusters] ef : change to smart-pointer + // int* mMatchedClustersIndex = nullptr; //[mNumOfClusters] + + std::unique_ptr mDBGOut; + UInt_t mDBGFlags = 0; + + ///----------- aux stuff --------------/// + static constexpr float MAXSNP = 0.85; // max snp of ITS or TPC track at xRef to be matched + + TStopwatch mTimerTot; + TStopwatch mTimerMatchITSTPC; + TStopwatch mTimerMatchTPC; + TStopwatch mTimerDBG; + + ClassDef(MatchHMP, 1); +}; +} // namespace globaltracking +} // namespace o2 + +#endif diff --git a/Detectors/GlobalTracking/src/MatchHMP.cxx b/Detectors/GlobalTracking/src/MatchHMP.cxx new file mode 100644 index 0000000000000..930c3c33efe6d --- /dev/null +++ b/Detectors/GlobalTracking/src/MatchHMP.cxx @@ -0,0 +1,596 @@ +// 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 "FairLogger.h" +#include "Field/MagneticField.h" +#include "Field/MagFieldFast.h" +#include "TOFBase/Geo.h" + +#include "SimulationDataFormat/MCTruthContainer.h" + +#include "DetectorsBase/Propagator.h" + +#include "MathUtils/Cartesian.h" +#include "MathUtils/Utils.h" +#include "CommonConstants/MathConstants.h" +#include "CommonConstants/PhysicsConstants.h" +#include "CommonConstants/GeomConstants.h" +#include "DetectorsBase/GeometryManager.h" + +#include +#include +#include +#include +#include "DataFormatsParameters/GRPObject.h" +#include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/TrackLTIntegral.h" + +#include "GlobalTracking/MatchHMP.h" + +#include "TPCBase/ParameterGas.h" +#include "TPCBase/ParameterElectronics.h" +#include "TPCReconstruction/TPCFastTransformHelperO2.h" + +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "DataFormatsGlobalTracking/RecoContainerCreateTracksVariadic.h" +#include "HMPIDBase/Param.h" + +using namespace o2::globaltracking; +using evGIdx = o2::dataformats::EvIndex; +using Cluster = o2::hmpid::Cluster; +using GTrackID = o2::dataformats::GlobalTrackID; +using timeEst = o2::dataformats::TimeStampWithError; + +ClassImp(MatchHMP); + +//______________________________________________ +void MatchHMP::run(const o2::globaltracking::RecoContainer& inp) +{ + ///< running the matching + mRecoCont = &inp; + mRecoCont = &inp; + mStartIR = inp.startIR; + + for (int i = 0; i < trkType::SIZEALL; i++) { + mMatchedTracks[i].clear(); + mOutHMPLabels[i].clear(); + } + + for (int i = 0; i < trkType::SIZE; i++) { + mTracksWork[i].clear(); + mTrackGid[i].clear(); + } + for (int it = 0; it < trkType::SIZE; it++) { + mMatchedTracksIndex[it].clear(); + mLTinfos[it].clear(); + if (mMCTruthON) { + mTracksLblWork[it].clear(); + } + for (int sec = o2::constants::math::NSectors; sec--;) { + mTracksSectIndexCache[it][sec].clear(); + } + } + + mSideTPC.clear(); + mExtraTPCFwdTime.clear(); + + // mTimerTot.Start(); + bool isPrepareHMPClusters = prepareHMPClusters(); + mTimerTot.Stop(); + LOGF(info, "Timing prepareTOFCluster: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1); + + if (!isPrepareHMPClusters) { // check cluster before of tracks to see also if MC is required + return; + } + + mTimerTot.Start(); + if (!prepareTracks()) { + return; + } + mTimerTot.Stop(); + LOGF(info, "Timing prepare TPC tracks: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1); + + mTimerTot.Start(); + + // if (!prepareFITData()) { ef : removed in HMP.h + // return; + // } + + mTimerTot.Stop(); + LOGF(info, "Timing prepare FIT data: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1); + + mTimerTot.Start(); + for (int sec = o2::constants::math::NSectors; sec--;) { + mMatchedTracksPairs.clear(); // sector + LOG(debug) << "Doing matching for sector " << sec << "..."; + if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) { + mTimerMatchITSTPC.Start(sec == o2::constants::math::NSectors - 1); + doMatching(sec); + mTimerMatchITSTPC.Stop(); + } + + /*if (mIsTPCused) { + mTimerMatchTPC.Start(sec == o2::constants::math::NSectors - 1); + // doMatchingForTPC(sec); ef : removed in HMP.h + mTimerMatchTPC.Stop(); + }*/ + + LOG(debug) << "...done. Now check the best matches"; + // selectBestMatches(); ef : removed in HMP + } + + // re-arrange outputs from constrained/unconstrained to the 4 cases (TPC, ITS-TPC, TPC-TRD, ITS-TPC-TRD) to be implemented as soon as TPC-TRD and ITS-TPC-TRD tracks available + + mIsTPCused = false; + mIsITSTPCused = false; + mIsTPCTRDused = false; + mIsITSTPCTRDused = false; + + mTimerTot.Stop(); + LOGF(info, "Timing Do Matching: Cpu: %.3e s Real: %.3e s in %d slots", mTimerTot.CpuTime(), mTimerTot.RealTime(), mTimerTot.Counter() - 1); + LOGF(info, "Timing Do Matching Constrained: Cpu: %.3e s Real: %.3e s in %d slots", mTimerMatchITSTPC.CpuTime(), mTimerMatchITSTPC.RealTime(), mTimerMatchITSTPC.Counter() - 1); + LOGF(info, "Timing Do Matching TPC : Cpu: %.3e s Real: %.3e s in %d slots", mTimerMatchTPC.CpuTime(), mTimerMatchTPC.RealTime(), mTimerMatchTPC.Counter() - 1); +} + +//______________________________________________ +bool MatchHMP::prepareTracks() +{ + mNotPropagatedToHMP[trkType::UNCONS] = 0; + mNotPropagatedToHMP[trkType::CONSTR] = 0; + + auto creator = [this](auto& trk, GTrackID gid, float time0, float terr) { + const int nclustersMin = 0; + if constexpr (isTPCTrack()) { + if (trk.getNClusters() < nclustersMin) { + return true; + } + + if (std::abs(trk.getQ2Pt()) > mMaxInvPt) { + return true; + } + } + if constexpr (isTPCITSTrack()) { + if (trk.getParamOut().getX() < o2::constants::geom::XTPCOuterRef - 1.) { + return true; + } + } + + return true; + }; + + mRecoCont->createTracksVariadic(creator); + + for (int it = 0; it < trkType::SIZE; it++) { + mMatchedTracksIndex[it].resize(mTracksWork[it].size()); + std::fill(mMatchedTracksIndex[it].begin(), mMatchedTracksIndex[it].end(), -1); // initializing all to -1 + } + + if (mIsTPCused) { + LOG(debug) << "Number of UNCONSTRAINED tracks that failed to be propagated to TOF = " << mNotPropagatedToHMP[trkType::UNCONS]; + + // sort tracks in each sector according to their time (increasing in time) + for (int sec = o2::constants::math::NSectors; sec--;) { + auto& indexCache = mTracksSectIndexCache[trkType::UNCONS][sec]; + LOG(debug) << "Sorting sector" << sec << " | " << indexCache.size() << " tracks"; + if (!indexCache.size()) { + continue; + } + std::sort(indexCache.begin(), indexCache.end(), [this](int a, int b) { + auto& trcA = mTracksWork[trkType::UNCONS][a].second; + auto& trcB = mTracksWork[trkType::UNCONS][b].second; + return ((trcA.getTimeStamp() - trcA.getTimeStampError()) - (trcB.getTimeStamp() - trcB.getTimeStampError()) < 0.); + }); + } // loop over tracks of single sector + } + if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused) { + LOG(debug) << "Number of CONSTRAINED tracks that failed to be propagated to TOF = " << mNotPropagatedToHMP[trkType::CONSTR]; // was mNotPropagatedToTOF + + // sort tracks in each sector according to their time (increasing in time) + for (int sec = o2::constants::math::NSectors; sec--;) { + auto& indexCache = mTracksSectIndexCache[trkType::CONSTR][sec]; + LOG(debug) << "Sorting sector" << sec << " | " << indexCache.size() << " tracks"; + if (!indexCache.size()) { + continue; + } + std::sort(indexCache.begin(), indexCache.end(), [this](int a, int b) { + auto& trcA = mTracksWork[trkType::CONSTR][a].second; + auto& trcB = mTracksWork[trkType::CONSTR][b].second; + return ((trcA.getTimeStamp() - mSigmaTimeCut * trcA.getTimeStampError()) - (trcB.getTimeStamp() - mSigmaTimeCut * trcB.getTimeStampError()) < 0.); + }); + } // loop over tracks of single sector + } + + return true; +} +//______________________________________________ +bool MatchHMP::prepareHMPClusters() +{ + mHMPClustersArrayInp = mRecoCont->getHMPClusters(); + mHMPClusLabels = mRecoCont->getHMPClustersMCLabels(); + mMCTruthON = mHMPClusLabels && mHMPClusLabels->getNElements(); + + ///< prepare the tracks that we want to match to HMPID + + // copy the track params, propagate to reference X and build sector tables + mHMPClusWork.clear(); + // mHMPClusWork.reserve(mNumOfClusters); // we cannot do this, we don't have mNumOfClusters yet + // if (mMCTruthON) { + // mTOFClusLblWork.clear(); + // mTOFClusLblWork.reserve(mNumOfClusters); + // } + + for (int sec = o2::constants::math::NSectors; sec--;) { + mHMPClusSectIndexCache[sec].clear(); + // mHMPClusSectIndexCache[sec].reserve(100 + 1.2 * mNumOfClusters / o2::constants::math::NSectors); // we cannot do this, we don't have mNumOfClusters yet + } + + int mNumOfClusters = 0; + + int nClusterInCurrentChunk = mHMPClustersArrayInp.size(); + LOG(debug) << "nClusterInCurrentChunk = " << nClusterInCurrentChunk; + mNumOfClusters += nClusterInCurrentChunk; + mHMPClusWork.reserve(mHMPClusWork.size() + mNumOfClusters); + for (int it = 0; it < nClusterInCurrentChunk; it++) { + const Cluster& clOrig = mHMPClustersArrayInp[it]; + // create working copy of track param + mHMPClusWork.emplace_back(clOrig); + auto& cl = mHMPClusWork.back(); + // cache work track index + // mHMPClusSectIndexCache[cl.getSector()].push_back(mHMPClusWork.size() - 1); fix + } // ef: o2::hmpid::Cluster has no type getSector(); I changed cluster to hmpid::cluster instead of tof + // + // sort clusters in each sector according to their time (increasing in time) + for (int sec = o2::constants::math::NSectors; sec--;) { + auto& indexCache = mHMPClusSectIndexCache[sec]; + LOG(debug) << "Sorting sector" << sec << " | " << indexCache.size() << " TOF clusters"; + if (!indexCache.size()) { + continue; + } + std::sort(indexCache.begin(), indexCache.end(), [this](int a, int b) { + auto& clA = mHMPClusWork[a]; + auto& clB = mHMPClusWork[b]; + // return (clA.getTime() - clB.getTime()) < 0.; fix ef: original statement + return 0; //(clA.getTime() - clB.getTime()) < 0.; //fix ef: o2::hmpid::Cluster + }); // has no type getTime() + } // loop over TOF clusters of single sector + + std::unique_ptr mMatchedClustersIndex; + mMatchedClustersIndex.reset(new int[mNumOfClusters]); + + /* ef : changed to smart/pointer + if (mMatchedClustersIndex) { + mMatchedClustersIndex; + } */ + //= std::unique_ptr(new int[mNumOfClusters]); // ef : change to smart-pointer + // std::fill_n(mMatchedClustersIndex, mNumOfClusters, -1); // initializing all to -1 + + // ef : want to avoid raw-pointers, + // are any of the following valid? : + + /* + for(auto& f : mMatchedClustersIndex){ + f = -1; + } */ + + for (int n = 0; n < mNumOfClusters; n++) { + mMatchedClustersIndex[n] = -1; + } + + return true; +} +//______________________________________________ +void MatchHMP::doMatching(int cham) +{ + trkType type = trkType::CONSTR; + + ///< do the real matching per chmaber + auto& cacheHMP = mHMPClusSectIndexCache[cham]; // array of cached HMP cluster indices for this chmaber; reminder: they are ordered in time! + auto& cacheTrk = mTracksSectIndexCache[type][cham]; // array of cached tracks indices for this chamber; reminder: they are ordered in time! + int nTracks = cacheTrk.size(), nHMPCls = cacheHMP.size(); + + // ef: changed here from sec 2 cham + LOG(debug) << "Matching sector " << cham << ": number of tracks: " << nTracks << ", number of HMP clusters: " << nHMPCls; + if (!nTracks || !nHMPCls) { + return; + } + int ihmp0 = 0; // starting index in HMP clusters for matching of the track + int detId[2][5]; // at maximum one track can fall in 2 strips during the propagation; the second dimention of the array is the TOF det index + float deltaPos[2][3]; // at maximum one track can fall in 2 strips during the propagation; the second dimention of the array is the residuals + o2::track::TrackLTIntegral trkLTInt[2]; // Here we store the integrated track length and time for the (max 2) matched strips + int nStepsInsideSameStrip[2] = {0, 0}; // number of propagation steps in the same strip (since we have maximum 2 strips, it has dimention = 2) + float deltaPosTemp[3]; + std::array pos; + std::array posBeforeProp; + float posFloat[3]; + + // prematching for TPC only tracks (identify BC candidate to correct z for TPC track accordingly to v_drift) + + LOG(debug) << "Trying to match %d tracks" << cacheTrk.size(); + for (int itrk = 0; itrk < cacheTrk.size(); itrk++) { + for (int ii = 0; ii < 2; ii++) { + detId[ii][2] = -1; // before trying to match, we need to inizialize the detId corresponding to the strip number to -1; this is the array that we will use to save the det id of the maximum 2 strips matched + nStepsInsideSameStrip[ii] = 0; + } + int nStripsCrossedInPropagation = 0; // how many strips were hit during the propagation + auto& trackWork = mTracksWork[type][cacheTrk[itrk]]; + auto& trefTrk = trackWork.first; + auto& intLT = mLTinfos[type][cacheTrk[itrk]]; + + // Printf("intLT (before doing anything): length = %f, time (Pion) = %f", intLT.getL(), intLT.getTOF(o2::track::PID::Pion)); + float minTrkTime = (trackWork.second.getTimeStamp() - mSigmaTimeCut * trackWork.second.getTimeStampError()) * 1.E6; // minimum time in ps + float maxTrkTime = (trackWork.second.getTimeStamp() + mSigmaTimeCut * trackWork.second.getTimeStampError()) * 1.E6; // maximum time in ps + int istep = 1; // number of steps + float step = 1.0; // step size in cm + + // uncomment for local debug + /* + //trefTrk.getXYZGlo(posBeforeProp); + //float posBeforeProp[3] = {trefTrk.getX(), trefTrk.getY(), trefTrk.getZ()}; // in local ref system + //printf("Global coordinates: posBeforeProp[0] = %f, posBeforeProp[1] = %f, posBeforeProp[2] = %f\n", posBeforeProp[0], posBeforeProp[1], posBeforeProp[2]); + //Printf("Radius xy = %f", TMath::Sqrt(posBeforeProp[0]*posBeforeProp[0] + posBeforeProp[1]*posBeforeProp[1])); + //Printf("Radius xyz = %f", TMath::Sqrt(posBeforeProp[0]*posBeforeProp[0] + posBeforeProp[1]*posBeforeProp[1] + posBeforeProp[2]*posBeforeProp[2])); + */ + + // initializing + for (int ii = 0; ii < 2; ii++) { + for (int iii = 0; iii < 5; iii++) { + detId[ii][iii] = -1; + } + } + + int detIdTemp[5] = {-1, -1, -1, -1, -1}; // HMP detector id at the current propagation point + + double reachedPoint = mXRef + istep * step; + + // modify this lines for the tracks propagation to the HMPID chambers + + /* while (propagateToRefX(trefTrk, reachedPoint, step, intLT) && nStripsCrossedInPropagation <= 2 && reachedPoint < Geo::RMAX) { + // while (o2::base::Propagator::Instance()->PropagateToXBxByBz(trefTrk, mXRef + istep * step, MAXSNP, step, 1, &intLT) && nStripsCrossedInPropagation <= 2 && mXRef + istep * step < Geo::RMAX) { + + trefTrk.getXYZGlo(pos); + for (int ii = 0; ii < 3; ii++) { // we need to change the type... + posFloat[ii] = pos[ii]; + }*/ + + // uncomment below only for local debug; this will produce A LOT of output - one print per propagation step + /* + Printf("posFloat[0] = %f, posFloat[1] = %f, posFloat[2] = %f", posFloat[0], posFloat[1], posFloat[2]); + Printf("radius xy = %f", TMath::Sqrt(posFloat[0]*posFloat[0] + posFloat[1]*posFloat[1])); + Printf("radius xyz = %f", TMath::Sqrt(posFloat[0]*posFloat[0] + posFloat[1]*posFloat[1] + posFloat[2]*posFloat[2])); + */ + + for (int idet = 0; idet < 5; idet++) { + detIdTemp[idet] = -1; + } + // cham is passed as input, but for tof it is sec? ef + // Geo::getPadDxDyDz(posFloat, detIdTemp, deltaPosTemp, cham); + // Geo::getPadDxDyDz(posFloat, detIdTemp, deltaPosTemp, sec); + reachedPoint += step; + + if (detIdTemp[2] == -1) { + continue; + } + + if (nStripsCrossedInPropagation == 0 || // we are crossing a strip for the first time... + (nStripsCrossedInPropagation >= 1 && (detId[nStripsCrossedInPropagation - 1][0] != detIdTemp[0] || detId[nStripsCrossedInPropagation - 1][1] != detIdTemp[1] || detId[nStripsCrossedInPropagation - 1][2] != detIdTemp[2]))) { // ...or we are crossing a new strip + if (nStripsCrossedInPropagation == 0) { + LOG(debug) << "We cross a strip for the first time"; + } + if (nStripsCrossedInPropagation == 2) { + break; // we have already matched 2 strips, we cannot match more + } + nStripsCrossedInPropagation++; + } + } + + /* for (Int_t imatch = 0; imatch < nStripsCrossedInPropagation; imatch++) { + // we take as residual the average of the residuals along the propagation in the same strip + deltaPos[imatch][0] /= nStepsInsideSameStrip[imatch]; + deltaPos[imatch][1] /= nStepsInsideSameStrip[imatch]; + deltaPos[imatch][2] /= nStepsInsideSameStrip[imatch]; + } + + if (nStripsCrossedInPropagation == 0) { + continue; // the track never hit a TOF strip during the propagation + }*/ + + bool foundCluster = false; + + for (auto ihmp = ihmp0; ihmp < nHMPCls; ihmp++) { + // printf("ihmp = %d\n", ihmp); + auto& trefTOF = mHMPClusWork[cacheHMP[ihmp]]; + // compare the times of the track and the TOF clusters - remember that they both are ordered in time! + // Printf("trefTOF.getTime() = %f, maxTrkTime = %f, minTrkTime = %f", trefTOF.getTime(), maxTrkTime, minTrkTime); + + /* // ef start comment if :no function getTime in hmp-clusters + if (trefTOF.getTime() < minTrkTime) { // this cluster has a time that is too small for the current track, we will get to the next one + //Printf("In trefTOF.getTime() < minTrkTime"); + ihmp0 = ihmp + 1; // but for the next track that we will check, we will ignore this cluster (the time is anyway too small) + continue; + } */ + // ef end comment if + + // ef no method getTime in hmpid: + // if (trefTOF.getTime() > maxTrkTime) { // no more TOF clusters can be matched to this track + // break; + //} + + // int mainChannel = trefTOF.getMainContributingChannel(); // ef: o2::hmpid::Cluster + // has no type getMainContributingChannel + + int indices[5]; + // Geo::getVolumeIndices(mainChannel, indices); + + // compute fine correction using cluster position instead of pad center + // this because in case of multiple-hit cluster position is averaged on all pads contributing to the cluster (then error position matrix can be used for Chi2 if nedeed) + int ndigits = 1; + float posCorr[3] = {0, 0, 0}; + + // ef: o2::hmpid::Cluster + // has no type isBitSet + /* ef: uncomment later + if (trefTOF.isBitSet(/*o2::tof::Cluster::*kLeft)) { // ef : temporary workaround bc + posCorr[0] += Geo::XPAD, ndigits++; // hmpid does not have enum-type + } // for bit + if (trefTOF.isBitSet(kUpLeft)) { + posCorr[0] += Geo::XPAD, posCorr[2] -= Geo::ZPAD, ndigits++; + } + if (trefTOF.isBitSet(kDownLeft)) { + posCorr[0] += Geo::XPAD, posCorr[2] += Geo::ZPAD, ndigits++; + } + if (trefTOF.isBitSet(kUp)) { + posCorr[2] -= Geo::ZPAD, ndigits++; + } + if (trefTOF.isBitSet(kDown)) { + posCorr[2] += Geo::ZPAD, ndigits++; + } + if (trefTOF.isBitSet(kRight)) { + posCorr[0] -= Geo::XPAD, ndigits++; + } + if (trefTOF.isBitSet(kUpRight)) { + posCorr[0] -= Geo::XPAD, posCorr[2] -= Geo::ZPAD, ndigits++; + } + if (trefTOF.isBitSet(kDownRight)) { + posCorr[0] -= Geo::XPAD, posCorr[2] += Geo::ZPAD, ndigits++; + } */ + // ef uncomment later + // ef: o2::hmpid::Cluster has no type isBitSet // + + float ndifInv = 1. / ndigits; + if (ndigits > 1) { + posCorr[0] *= ndifInv; + posCorr[1] *= ndifInv; + posCorr[2] *= ndifInv; + } + + int trackIdTOF; + int eventIdTOF; + int sourceIdTOF; + + /*for (auto iPropagation = 0; iPropagation < nStripsCrossedInPropagation; iPropagation++) { + LOG(debug) << "TOF Cluster [" << ihmp << ", " << cacheHMP[ihmp] << "]: indices = " << indices[0] << ", " << indices[1] << ", " << indices[2] << ", " << indices[3] << ", " << indices[4]; + LOG(debug) << "Propagated Track [" << itrk << "]: detId[" << iPropagation << "] = " << detId[iPropagation][0] << ", " << detId[iPropagation][1] << ", " << detId[iPropagation][2] << ", " << detId[iPropagation][3] << ", " << detId[iPropagation][4]; + float resX = deltaPos[iPropagation][0] - (indices[4] - detId[iPropagation][4]) * Geo::XPAD + posCorr[0]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster + float resZ = deltaPos[iPropagation][2] - (indices[3] - detId[iPropagation][3]) * Geo::ZPAD + posCorr[2]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster + float res = TMath::Sqrt(resX * resX + resZ * resZ); + + LOG(debug) << "resX = " << resX << ", resZ = " << resZ << ", res = " << res; + if (indices[0] != detId[iPropagation][0]) { + continue; + } + if (indices[1] != detId[iPropagation][1]) { + continue; + } + if (indices[2] != detId[iPropagation][2]) { + continue; + } + float chi2 = res; // TODO: take into account also the time! + + if (res < mSpaceTolerance) { // matching ok! + LOG(debug) << "MATCHING FOUND: We have a match! between track " << mTracksSectIndexCache[type][cham][itrk] << " and TOF cluster " << mHMPClusSectIndexCache[indices[0]][ihmp]; // sec2cham + foundCluster = true; + // set event indexes (to be checked) + int eventIndexTOFCluster = mHMPClusSectIndexCache[indices[0]][ihmp]; + + // ef commented out next line, no method named getTime() + //mMatchedTracksPairs.emplace_back(cacheTrk[itrk], eventIndexTOFCluster, mHMPClusWork[cacheHMP[ihmp]].getTime(), chi2, trkLTInt[iPropagation], mTrackGid[type][cacheTrk[itrk]], type, (trefTOF.getTime() - (minTrkTime + maxTrkTime) * 0.5) * 1E-6, 0., resX, resZ); // TODO: check if this is correct! + } + } + } */ + // end for-loop + } + return; +} + +//______________________________________________ +bool MatchHMP::propagateToRefX(o2::track::TrackParCov& trc, float xRef, float stepInCm, o2::track::TrackLTIntegral& intLT) +{ + // propagate track to matching reference X + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; // material correction method + static const float tanHalfSector = tan(o2::constants::math::SectorSpanRad / 2); + bool refReached = false; + float xStart = trc.getX(); + // the first propagation will be from 2m, if the track is not at least at 2m + if (xStart < 50.) { + xStart = 50.; + } + int istep = 1; + bool hasPropagated = o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, xStart + istep * stepInCm, MAXSNP, stepInCm, matCorr, &intLT); + while (hasPropagated) { + if (trc.getX() > xRef) { + refReached = true; // we reached the 371cm reference + } + istep++; + if (fabs(trc.getY()) > trc.getX() * tanHalfSector) { // we are still in the same sector + // we need to rotate the track to go to the new sector + // Printf("propagateToRefX: changing sector"); + auto alphaNew = o2::math_utils::angle2Alpha(trc.getPhiPos()); + if (!trc.rotate(alphaNew) != 0) { + // Printf("propagateToRefX: failed to rotate"); + break; // failed (this line is taken from MatchTPCITS and the following comment too: RS: check effect on matching tracks to neighbouring sector) + } + } + if (refReached) { + break; + } + hasPropagated = o2::base::Propagator::Instance()->PropagateToXBxByBz(trc, xStart + istep * stepInCm, MAXSNP, stepInCm, matCorr, &intLT); + } + + // if (std::abs(trc.getSnp()) > MAXSNP) Printf("propagateToRefX: condition on snp not ok, returning false"); + // Printf("propagateToRefX: snp of teh track is %f (--> %f grad)", trc.getSnp(), TMath::ASin(trc.getSnp())*TMath::RadToDeg()); + return refReached && std::abs(trc.getSnp()) < 0.95; // Here we need to put MAXSNP +} + +//______________________________________________ +bool MatchHMP::propagateToRefXWithoutCov(o2::track::TrackParCov& trc, float xRef, float stepInCm, float bzField) +{ + // propagate track to matching reference X without using the covariance matrix + // we create the copy of the track in a TrackPar object (no cov matrix) + o2::track::TrackPar trcNoCov(trc); + const float tanHalfSector = tan(o2::constants::math::SectorSpanRad / 2); + bool refReached = false; + float xStart = trcNoCov.getX(); + // the first propagation will be from 2m, if the track is not at least at 2m + if (xStart < 50.) { + xStart = 50.; + } + int istep = 1; + bool hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField); + while (hasPropagated) { + if (trcNoCov.getX() > xRef) { + refReached = true; // we reached the 371cm reference + } + istep++; + if (fabs(trcNoCov.getY()) > trcNoCov.getX() * tanHalfSector) { // we are still in the same sector + // we need to rotate the track to go to the new sector + // Printf("propagateToRefX: changing sector"); + auto alphaNew = o2::math_utils::angle2Alpha(trcNoCov.getPhiPos()); + if (!trcNoCov.rotateParam(alphaNew) != 0) { + // Printf("propagateToRefX: failed to rotate"); + break; // failed (this line is taken from MatchTPCITS and the following comment too: RS: check effect on matching tracks to neighbouring sector) + } + } + if (refReached) { + break; + } + hasPropagated = trcNoCov.propagateParamTo(xStart + istep * stepInCm, bzField); + } + // if (std::abs(trc.getSnp()) > MAXSNP) Printf("propagateToRefX: condition on snp not ok, returning false"); + // Printf("propagateToRefX: snp of teh track is %f (--> %f grad)", trcNoCov.getSnp(), TMath::ASin(trcNoCov.getSnp())*TMath::RadToDeg()); + + // return refReached && std::abs(trcNoCov.getSnp()) < 0.95 && TMath::Abs(trcNoCov.getZ()) < Geo::MAXHZTOF; // Here we need to put MAXSNP + + return 999.; +} + +//______________________________________________ diff --git a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt index 0079fbc415b68..d3c1cab08e82d 100644 --- a/Detectors/GlobalTrackingWorkflow/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/CMakeLists.txt @@ -23,6 +23,7 @@ o2_add_library(GlobalTrackingWorkflow src/CosmicsMatchingSpec.cxx src/TrackCosmicsWriterSpec.cxx src/TOFMatcherSpec.cxx + src/HMPMatcherSpec.cxx src/TOFMatchChecker.cxx src/TOFEventTimeChecker.cxx src/GlobalFwdTrackWriterSpec.cxx @@ -36,6 +37,7 @@ o2_add_library(GlobalTrackingWorkflow O2::ITSWorkflow O2::MFTTracking O2::MFTWorkflow + O2::HMPIDWorkflow O2::TPCWorkflow O2::FT0Workflow O2::ITSMFTWorkflow diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h new file mode 100644 index 0000000000000..d853e7bd39a05 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h @@ -0,0 +1,33 @@ +// 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. + +/// @file HMPMatcherSpec.h // ef ; change to hmp + +#ifndef O2_HMP_MATCHER_SPEC // hmp +#define O2_HMP_MATCHER_SPEC // + +#include "Framework/DataProcessorSpec.h" +#include "ReconstructionDataFormats/MatchInfoTOFReco.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace globaltracking +{ + +/// create a processor spec +framework::DataProcessorSpec getHMPMatcherSpec(o2::dataformats::GlobalTrackID::mask_t src, bool useMC, bool useFIT, bool tpcRefit, bool strict); + +} // namespace globaltracking +} // namespace o2 + +#endif /* O2_HMP_MATCHER_SPEC */ diff --git a/Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx new file mode 100644 index 0000000000000..c517132650ffa --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx @@ -0,0 +1,225 @@ +// 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. + +/// @file HMPMatcherSpec.cxx + +#include +#include +#include "TStopwatch.h" +#include "Framework/ConfigParamRegistry.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "CommonUtils/NameConf.h" +#include "DataFormatsParameters/GRPObject.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "Framework/Task.h" +#include "Framework/DataProcessorSpec.h" + +// from Tracks +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "ReconstructionDataFormats/GlobalTrackAccessor.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DataFormatsTPC/TrackTPC.h" +#include "DataFormatsITS/TrackITS.h" +#include "ReconstructionDataFormats/TrackTPCITS.h" +#include "ReconstructionDataFormats/TrackTPCTOF.h" + +// from HMPID +#include "DataFormatsHMP/Cluster.h" +#include "GlobalTracking/MatchHMP.h" +#include "GlobalTrackingWorkflow/HMPMatcherSpec.h" + +using namespace o2::framework; +// using MCLabelsTr = gsl::span; +// using GID = o2::dataformats::GlobalTrackID; +// using DetID = o2::detectors::DetID; + +using evIdx = o2::dataformats::EvIndex; +using MatchOutputType = std::vector; +using GID = o2::dataformats::GlobalTrackID; + +namespace o2 +{ +namespace globaltracking +{ + +class HMPMatcherSpec : public Task +{ + public: + HMPMatcherSpec(std::shared_ptr dr, bool useMC, bool useFIT, bool tpcRefit, bool strict) : mDataRequest(dr), mUseMC(useMC), mUseFIT(useFIT), mDoTPCRefit(tpcRefit), mStrict(strict) {} + ~HMPMatcherSpec() override = default; + void init(InitContext& ic) final; + void run(ProcessingContext& pc) final; + void endOfStream(framework::EndOfStreamContext& ec) final; + + private: + std::shared_ptr mDataRequest; + bool mUseMC = true; + bool mUseFIT = false; + bool mDoTPCRefit = false; + bool mStrict = false; + MatchHMP mMatcher; ///< Cluster finder + TStopwatch mTimer; +}; + +void HMPMatcherSpec::init(InitContext& ic) +{ + mTimer.Stop(); + mTimer.Reset(); + //-------- init geometry and field --------// + o2::base::GeometryManager::loadGeometry(); + o2::base::Propagator::initFieldFromGRP(); + std::unique_ptr grp{o2::parameters::GRPObject::loadFrom()}; + + // this is a hack to provide Mat.LUT from the local file, in general will be provided by the framework from CCDB + std::string matLUTPath = ic.options().get("material-lut-path"); + std::string matLUTFile = o2::base::NameConf::getMatLUTFileName(matLUTPath); + if (o2::utils::Str::pathExists(matLUTFile)) { + auto* lut = o2::base::MatLayerCylSet::loadFromFile(matLUTFile); + o2::base::Propagator::Instance()->setMatLUT(lut); + LOG(debug) << "Loaded material LUT from " << matLUTFile; + } else { + LOG(debug) << "Material LUT " << matLUTFile << " file is absent, only TGeo can be used"; + } + if (mStrict) { + // mMatcher.setHighPurity(); ef: function is commented out in MatchHMP.h + } +} + +void HMPMatcherSpec::run(ProcessingContext& pc) +{ + mTimer.Start(false); + + RecoContainer recoData; + recoData.collectData(pc, *mDataRequest.get()); + + auto creationTime = DataRefUtils::getHeader(pc.inputs().getFirstValid(true))->creation; + + LOG(debug) << "isTrackSourceLoaded: TPC -> " << recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::TPC); + LOG(debug) << "isTrackSourceLoaded: ITSTPC -> " << recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::ITSTPC); + LOG(debug) << "isTrackSourceLoaded: TPCTRD -> " << recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::TPCTRD); + LOG(debug) << "isTrackSourceLoaded: ITSTPCTRD -> " << recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::ITSTPCTRD); + + bool isTPCused = recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::TPC); + bool isITSTPCused = recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::ITSTPC); + bool isTPCTRDused = recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::TPCTRD); + bool isITSTPCTRDused = recoData.isTrackSourceLoaded(o2::dataformats::GlobalTrackID::Source::ITSTPCTRD); + uint32_t ss = o2::globaltracking::getSubSpec(mStrict ? o2::globaltracking::MatchingType::Strict : o2::globaltracking::MatchingType::Standard); + + // mMatcher.setFIT(mUseFIT); + + // mMatcher.setTS(creationTime); + + mMatcher.run(recoData); + + if (isTPCused) { + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MTC_TPC", ss, Lifetime::Timeframe}, mMatcher.getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC)); + if (mUseMC) { + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MCMTC_TPC", ss, Lifetime::Timeframe}, mMatcher.getMatchedHMPLabelsVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC)); + // ef : what should the namescope of TrackType be? + } // ef changed to getMatchedHMPLabelsVector + + auto nmatch = mMatcher.getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC).size(); + if (mDoTPCRefit) { + LOG(debug) << "Refitting " << nmatch << " matched TPC tracks with TOF time info"; + } else { + LOG(debug) << "Shifting Z for " << nmatch << " matched TPC tracks according to TOF time info"; + } + auto& tracksTPCTOF = pc.outputs().make>(OutputRef{"tpctofTracks", ss}, nmatch); + // mMatcher.makeConstrainedTPCTracks(tracksTPCTOF); + } + + if (isITSTPCused) { + pc.outputs().snapshot(Output{o2::header::gDataOriginHMP, "MTC_ITSTPC", 0, Lifetime::Timeframe}, mMatcher.getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::ITSTPC)); + if (mUseMC) { + pc.outputs().snapshot(Output{o2::header::gDataOriginHMP, "MCMTC_ITSTPC", 0, Lifetime::Timeframe}, mMatcher.getMatchedHMPLabelsVector(o2::dataformats::MatchInfoTOFReco::TrackType::ITSTPC)); + } + } + + if (isTPCTRDused) { + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MTC_TPCTRD", ss, Lifetime::Timeframe}, mMatcher.getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPCTRD)); + if (mUseMC) { + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MCMTC_TPCTRD", ss, Lifetime::Timeframe}, mMatcher.getMatchedHMPLabelsVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPCTRD)); + } + } + + if (isITSTPCTRDused) { + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MTC_ITSTPCTRD", 0, Lifetime::Timeframe}, mMatcher.getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::ITSTPCTRD)); + if (mUseMC) { + pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MCMTC_ITSTPCTRD", 0, Lifetime::Timeframe}, mMatcher.getMatchedHMPLabelsVector(o2::dataformats::MatchInfoTOFReco::TrackType::ITSTPCTRD)); + } + } + + // TODO: TRD-matched tracks + // pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "CALIBDATA", 0, Lifetime::Timeframe}, mMatcher.getCalibVector()); + + mTimer.Stop(); +} + +void HMPMatcherSpec::endOfStream(EndOfStreamContext& ec) +{ + LOGF(debug, "HMP matching total timing: Cpu: %.3e Real: %.3e s in %d slots", + mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); +} + +DataProcessorSpec getHMPMatcherSpec(GID::mask_t src, bool useMC, bool useFIT, bool tpcRefit, bool strict) +{ + uint32_t ss = o2::globaltracking::getSubSpec(strict ? o2::globaltracking::MatchingType::Strict : o2::globaltracking::MatchingType::Standard); + auto dataRequest = std::make_shared(); + if (strict) { + dataRequest->setMatchingInputStrict(); + } + dataRequest->requestTracks(src, useMC); + dataRequest->requestClusters(GID::getSourceMask(GID::TOF), useMC); + if (useFIT) { + dataRequest->requestClusters(GID::getSourceMask(GID::FT0), false); + } + + std::vector outputs; + if (GID::includesSource(GID::TPC, src)) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MTC_TPC", ss, Lifetime::Timeframe); + // outputs.emplace_back(OutputLabel{"tpctofTracks"}, o2::header::gDataOriginTOF, "TOFTRACKS_TPC", ss, Lifetime::Timeframe); + if (useMC) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MCMTC_TPC", ss, Lifetime::Timeframe); + } + } + if (GID::includesSource(GID::ITSTPC, src)) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MTC_ITSTPC", 0, Lifetime::Timeframe); + if (useMC) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MCMTC_ITSTPC", 0, Lifetime::Timeframe); + } + } + if (GID::includesSource(GID::ITSTPCTRD, src)) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MTC_ITSTPCTRD", 0, Lifetime::Timeframe); + if (useMC) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MCMTC_ITSTPCTRD", 0, Lifetime::Timeframe); + } + } + if (GID::includesSource(GID::TPCTRD, src)) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MTC_TPCTRD", ss, Lifetime::Timeframe); + if (useMC) { + outputs.emplace_back(o2::header::gDataOriginHMP, "MCMTC_TPCTRD", ss, Lifetime::Timeframe); + } + } + // outputs.emplace_back(o2::header::gDataOriginTOF, "CALIBDATA", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "hmp-matcher", + dataRequest->inputs, + outputs, + AlgorithmSpec{adaptFromTask(dataRequest, useMC, useFIT, tpcRefit, strict)}, + Options{ + {"material-lut-path", VariantType::String, "", {"Path of the material LUT file"}}}}; +} + +} // namespace globaltracking +} // namespace o2 diff --git a/Detectors/HMPID/base/CMakeLists.txt b/Detectors/HMPID/base/CMakeLists.txt index da9b9e77b73d0..dd57151bee537 100644 --- a/Detectors/HMPID/base/CMakeLists.txt +++ b/Detectors/HMPID/base/CMakeLists.txt @@ -10,6 +10,7 @@ # or submit itself to any jurisdiction. o2_add_library(HMPIDBase + TARGETVARNAME targetName SOURCES src/Param.cxx src/Geo.cxx PUBLIC_LINK_LIBRARIES O2::CommonDataFormat O2::SimulationDataFormat diff --git a/Detectors/HMPID/base/include/HMPIDBase/Param.h b/Detectors/HMPID/base/include/HMPIDBase/Param.h index 9d5b044e155cf..69dd471f409ee 100644 --- a/Detectors/HMPID/base/include/HMPIDBase/Param.h +++ b/Detectors/HMPID/base/include/HMPIDBase/Param.h @@ -16,11 +16,25 @@ #include #include //base class #include //Instance() -#include //Lors2Mars() Mars2Lors() + +// ef : to use XYZVector +#include //fields +#include "Math/Vector3D.h" //fields + +//#include //Lors2Mars() Mars2Lors() + +/* are these necessary?: +#include //ef +#include "Math/GenVector/RotationX.h" //ef +#include "Math/GenVector/RotationY.h" //ef +#include "Math/GenVector/RotationZ.h" //ef +*/ class TGeoVolume; class TGeoHMatrix; +using XYZVector = ROOT::Math::XYZVector; + namespace o2 { namespace hmpid @@ -100,8 +114,33 @@ class Param bool getInstType() const { return fgInstanceType; } // return if the instance is from geom or ideal - inline static bool isInDead(float x, float y); // is the point in dead area? - inline static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch); // is a dead pad? + // is the point in dead area? + static bool isInDead(float x, float y) + { // ef : moved function-definition from cxx + // Check is the current point is outside of sensitive area or in dead zones + // Arguments: x,y -position + // Returns: 1 if not in sensitive zone + for (Int_t iPc = 0; iPc < 6; iPc++) { + if (x >= fgkMinPcX[iPc] && x <= fgkMaxPcX[iPc] && y >= fgkMinPcY[iPc] && y <= fgkMaxPcY[iPc]) { + return kFALSE; // in current pc + } + } + return kTRUE; + } + + static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch) // ef : moved function-definition from cxx + { + + // Check is the current pad is active or not + // Arguments: padx,pady pad integer coord + // Returns: kTRUE if dead, kFALSE if active + + if (fgMapPad[padx - 1][pady - 1][ch]) { + return kFALSE; // current pad active + } + + return kTRUE; + } inline void setChStatus(Int_t ch, bool status = kTRUE); inline void setSectStatus(Int_t ch, Int_t sect, bool status); @@ -188,12 +227,17 @@ class Param double l[3] = {x - mX, y - mY, z}; mM[c]->LocalToMaster(l, m); } - TVector3 lors2Mars(Int_t c, double x, double y, Int_t pl = kPc) const + + // template + XYZVector lors2Mars(Int_t c, double x, double y, Int_t pl = kPc) const { double m[3]; lors2Mars(c, x, y, m, pl); - return TVector3(m); + + return XYZVector(m[0], m[1], m[2]); // TVector3(m); + } // MRS->LRS + void mars2Lors(Int_t c, double* m, double& x, double& y) const { double l[3]; @@ -210,12 +254,13 @@ class Param ph = TMath::ATan2(l[1], l[0]); } void lors2MarsVec(Int_t c, double* m, double* l) const { mM[c]->LocalToMasterVect(m, l); } // LRS->MRS - TVector3 norm(Int_t c) const + + XYZVector norm(Int_t c) const // TVector3 { double n[3]; norm(c, n); - return TVector3(n); - } // norm + return XYZVector(n[0], n[1], n[2]); // TVector3(n); + } // norm void norm(Int_t c, double* n) const { double l[3] = {0, 0, 1}; diff --git a/Detectors/HMPID/base/src/Geo.cxx b/Detectors/HMPID/base/src/Geo.cxx index 6615bdaae3e74..ed065661d86fc 100644 --- a/Detectors/HMPID/base/src/Geo.cxx +++ b/Detectors/HMPID/base/src/Geo.cxx @@ -19,7 +19,7 @@ #include "HMPIDBase/Param.h" #include "TGeoManager.h" #include "TMath.h" -#include +#include "FairLogger.h" #include "DetectorsBase/GeometryManager.h" ClassImp(o2::hmpid::Geo); diff --git a/Detectors/HMPID/base/src/Param.cxx b/Detectors/HMPID/base/src/Param.cxx index 9696067801fd6..61e4ef6d8358c 100644 --- a/Detectors/HMPID/base/src/Param.cxx +++ b/Detectors/HMPID/base/src/Param.cxx @@ -14,13 +14,22 @@ #include //TestTrans() #include //TestTrans() #include //TestTrans() -#include + +// #include // ef: this is not used? +// TGeoHMatrix is used; could not find any indication of that being dis-advised + #include //Stack() #include //ctor #include #include //ctor #include +// ef : rotations +#include //ef +#include "Math/GenVector/RotationX.h" //ef +#include "Math/GenVector/RotationY.h" //ef +#include "Math/GenVector/RotationZ.h" //ef + using namespace o2::hmpid; ClassImp(o2::hmpid::Param); @@ -57,14 +66,14 @@ float Param::fgAllY = 0; bool Param::fgInstanceType = kTRUE; -Param* Param::fgInstance = nullptr; //singleton pointer +Param* Param::fgInstance = nullptr; // singleton pointer Int_t Param::fgNSigmas = 4; Int_t Param::fgThreshold = 4; //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Param::Param(bool noGeo) : mX(0), mY(0), mRefIdx(1.28947), mPhotEMean(6.675), mTemp(25) -//just set a refractive index for C6F14 at ephot=6.675 eV @ T=25 C +// just set a refractive index for C6F14 at ephot=6.675 eV @ T=25 C { // Here all the intitializition is taken place when Param::Instance() is invoked for the first time. // In particular, matrices to be used for LORS<->MARS trasnformations are initialized from TGeo structure. @@ -90,12 +99,12 @@ Param::Param(bool noGeo) : mX(0), mY(0), mRefIdx(1.28947), mPhotEMean(6.675), mT } } */ - mRefIdx = meanIdxRad(); //initialization of the running ref. index of freon + mRefIdx = meanIdxRad(); // initialization of the running ref. index of freon float dead = 2.6; // cm of the dead zones between PCs-> See 2CRC2099P1 if (noGeo == kTRUE) { - fgInstanceType = kFALSE; //instance from ideal geometry, no actual geom is present + fgInstanceType = kFALSE; // instance from ideal geometry, no actual geom is present } if (noGeo == kFALSE && !gGeoManager) { @@ -148,7 +157,7 @@ Param::Param(bool noGeo) : mX(0), mY(0), mRefIdx(1.28947), mPhotEMean(6.675), mT for (Int_t ich = kMinCh; ich <= kMaxCh; ich++) { for (Int_t padx = 0; padx < 160; padx++) { for (Int_t pady = 0; pady < 144; pady++) { - fgMapPad[padx][pady][ich] = kTRUE; //init all the pads are active at the beginning.... + fgMapPad[padx][pady][ich] = kTRUE; // init all the pads are active at the beginning.... } } } @@ -157,7 +166,7 @@ Param::Param(bool noGeo) : mX(0), mY(0), mRefIdx(1.28947), mPhotEMean(6.675), mT if (gGeoManager && gGeoManager->IsClosed()) { TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(Form("/HMPID/Chamber%i", i)); if (!pne) { - //AliErrorClass(Form("The symbolic volume %s does not correspond to any physical entry!",Form("HMPID_%i",i))); + // AliErrorClass(Form("The symbolic volume %s does not correspond to any physical entry!",Form("HMPID_%i",i))); mM[i] = new TGeoHMatrix; idealPosition(i, mM[i]); } else { @@ -175,7 +184,7 @@ Param::Param(bool noGeo) : mX(0), mY(0), mRefIdx(1.28947), mPhotEMean(6.675), mT } } fgInstance = this; -} //ctor +} // ctor //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Param::print(Option_t* opt) const { @@ -184,7 +193,7 @@ void Param::print(Option_t* opt) const for (Int_t i = 0; i < 7; i++) { mM[i]->Print(opt); } -} //Print() +} // Print() //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Param::idealPosition(Int_t iCh, TGeoHMatrix* pMatrix) { @@ -201,46 +210,46 @@ void Param::idealPosition(Int_t iCh, TGeoHMatrix* pMatrix) case 0: pMatrix->RotateY(kAngHor); pMatrix->RotateZ(-kAngVer); - break; //right and down + break; // right and down case 1: pMatrix->RotateZ(-kAngVer); - break; //down + break; // down case 2: pMatrix->RotateY(kAngHor); - break; //right + break; // right case 3: - break; //no rotation + break; // no rotation case 4: pMatrix->RotateY(-kAngHor); - break; //left + break; // left case 5: pMatrix->RotateZ(kAngVer); - break; //up + break; // up case 6: pMatrix->RotateY(-kAngHor); pMatrix->RotateZ(kAngVer); - break; //left and up + break; // left and up } - pMatrix->RotateZ(kAngCom); //apply common rotation in XY plane + pMatrix->RotateZ(kAngCom); // apply common rotation in XY plane } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /*Int_t Param::Stack(Int_t evt,Int_t tid) { // Prints some useful info from stack // Arguments: evt - event number. if not -1 print info only for that event -// tid - track id. if not -1 then print it and all it's mothers if any +// tid - track id. if not -1 then print it and all it's mothers if any // Returns: mother tid of the given tid if any - AliRunLoader *pAL=AliRunLoader::Open(); + AliRunLoader *pAL=AliRunLoader::Open(); if(pAL->LoadHeader()) return -1; if(pAL->LoadKinematics()) return -1; - + Int_t mtid=-1; Int_t iNevt=pAL->GetNumberOfEvents(); - + for(Int_t iEvt=0;iEvtGetEvent(iEvt); - AliStack *pStack=pAL->Stack(); + pAL->GetEvent(iEvt); + AliStack *pStack=pAL->Stack(); if(tid==-1){ //print all tids for this event for(Int_t i=0;iGetNtrack();i++) pStack->Particle(i)->Print(); Printf("totally %i tracks including %i primaries for event %i out of %i event(s)", @@ -253,8 +262,8 @@ void Param::idealPosition(Int_t iCh, TGeoHMatrix* pMatrix) while((tid=pTrack->GetFirstMother()) >= 0){ pTrack=pStack->Particle(tid); str+=" from ";str+=pTrack->GetName(); - } - }//if(tid==-1) + } + }//if(tid==-1) }//events loop pAL->UnloadHeader(); pAL->UnloadKinematics(); return mtid; @@ -263,15 +272,15 @@ void Param::idealPosition(Int_t iCh, TGeoHMatrix* pMatrix) /*Int_t Param::StackCount(Int_t pid,Int_t evt) { // Counts total number of particles of given sort (including secondary) for a given event - AliRunLoader *pAL=AliRunLoader::Open(); - pAL->GetEvent(evt); + AliRunLoader *pAL=AliRunLoader::Open(); + pAL->GetEvent(evt); if(pAL->LoadHeader()) return 0; if(pAL->LoadKinematics()) return 0; AliStack *pStack=pAL->Stack(); - + Int_t iCnt=0; for(Int_t i=0;iGetNtrack();i++) if(pStack->Particle(i)->GetPdgCode()==pid) iCnt++; - + pAL->UnloadHeader(); pAL->UnloadKinematics(); return iCnt; }*/ @@ -286,16 +295,21 @@ double Param::sigma2(double trkTheta, double trkPhi, double ckovTh, double ckovP // MIP beta // Returns: absolute error on Cerenkov angle, [radians] - TVector3 v(-999, -999, -999); + // TVector3 v(-999, -999, -999); : ef : changed to : + XYZVector v(-999, -999, -999); + double trkBeta = 1. / (TMath::Cos(ckovTh) * getRefIdx()); if (trkBeta > 1) { - trkBeta = 1; //protection against bad measured thetaCer + trkBeta = 1; // protection against bad measured thetaCer } if (trkBeta < 0) { trkBeta = 0.0001; // } + // ef : why not initialize directly to these values? //math_utils::Vector3D v(sigLoc, sigGeom, sigCrom); + + // ef : following methods are valid for DisplacementVector3D v.SetX(sigLoc(trkTheta, trkPhi, ckovTh, ckovPh, trkBeta)); v.SetY(sigGeom(trkTheta, trkPhi, ckovTh, ckovPh, trkBeta)); v.SetZ(sigCrom(trkTheta, trkPhi, ckovTh, ckovPh, trkBeta)); @@ -337,7 +351,7 @@ double Param::sigLoc(double trkTheta, double trkPhi, double thetaC, double phiC, // formula (7) pag.4 double dtdyc = kk * (k * (cosfd * sinf + cost * sinfd * cosf) + (alpha * e / (betaM * betaM)) * sint * sinfd); - double errX = 0.2, errY = 0.25; //end of page 7 + double errX = 0.2, errY = 0.25; // end of page 7 return TMath::Sqrt(errX * errX * dtdxc * dtdxc + errY * errY * dtdyc * dtdyc); } //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -365,7 +379,7 @@ double Param::sigCrom(double trkTheta, double trkPhi, double thetaC, double phiC double f = 0.0172 * (7.75 - 5.635) / TMath::Sqrt(24.); return f * dtdn; -} //SigCrom() +} // SigCrom() //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double Param::sigGeom(double trkTheta, double trkPhi, double thetaC, double phiC, double betaM) { @@ -404,7 +418,7 @@ double Param::sigGeom(double trkTheta, double trkPhi, double thetaC, double phiC double trErr = radThick() / (TMath::Sqrt(12.) * cost); return trErr * dtdT; -} //SigGeom() +} // SigGeom() //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ double Param::sigmaCorrFact(Int_t iPart, double occupancy) { @@ -437,10 +451,10 @@ Param* Param::instance() // Arguments: none // Returns: pointer to the instance of AliHMPIDParam or 0 if no geometry if (!fgInstance) { - new Param(kFALSE); //default setting for reconstruction, if no geometry.root -> AliFatal + new Param(kFALSE); // default setting for reconstruction, if no geometry.root -> AliFatal } return fgInstance; -} //Instance() +} // Instance() //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Param* Param::instanceNoGeo() @@ -449,12 +463,14 @@ Param* Param::instanceNoGeo() // Arguments: none // Returns: pointer to the instance of AliHMPIDParam or 0 if no geometry if (!fgInstance) { - new Param(kTRUE); //to avoid AliFatal, for MOOD and displays, use ideal geometry parameters + new Param(kTRUE); // to avoid AliFatal, for MOOD and displays, use ideal geometry parameters } return fgInstance; -} //Instance() +} // Instance() //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -bool Param::isInDead(float x, float y) + +/* ef : moved to header +inline bool Param::isInDead(float x, float y) // ef : bool -> inline bool { // Check is the current point is outside of sensitive area or in dead zones // Arguments: x,y -position @@ -466,9 +482,11 @@ bool Param::isInDead(float x, float y) } return kTRUE; -} +} */ //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -bool Param::isDeadPad(Int_t padx, Int_t pady, Int_t ch) + +/*ef : moved to header +inline bool Param::isDeadPad(Int_t padx, Int_t pady, Int_t ch) // ef : bool -> inline bool { // Check is the current pad is active or not // Arguments: padx,pady pad integer coord @@ -479,8 +497,9 @@ bool Param::isDeadPad(Int_t padx, Int_t pady, Int_t ch) } return kTRUE; -} +} */ //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + void Param::lors2Pad(float x, float y, Int_t& pc, Int_t& px, Int_t& py) { // Check the pad of given position @@ -490,25 +509,25 @@ void Param::lors2Pad(float x, float y, Int_t& pc, Int_t& px, Int_t& py) if (x > fgkMinPcX[0] && x < fgkMaxPcX[0]) { pc = 0; px = Int_t(x / sizePadX()); - } //PC 0 or 2 or 4 + } // PC 0 or 2 or 4 else if (x > fgkMinPcX[1] && x < fgkMaxPcX[1]) { pc = 1; px = Int_t((x - fgkMinPcX[1]) / sizePadX()); - } //PC 1 or 3 or 5 + } // PC 1 or 3 or 5 else { return; } if (y > fgkMinPcY[0] && y < fgkMaxPcY[0]) { py = Int_t(y / sizePadY()); - } //PC 0 or 1 + } // PC 0 or 1 else if (y > fgkMinPcY[2] && y < fgkMaxPcY[2]) { pc += 2; py = Int_t((y - fgkMinPcY[2]) / sizePadY()); - } //PC 2 or 3 + } // PC 2 or 3 else if (y > fgkMinPcY[4] && y < fgkMaxPcY[4]) { pc += 4; py = Int_t((y - fgkMinPcY[4]) / sizePadY()); - } //PC 4 or 5 + } // PC 4 or 5 else { return; } @@ -516,9 +535,9 @@ void Param::lors2Pad(float x, float y, Int_t& pc, Int_t& px, Int_t& py) //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Int_t Param::inHVSector(float y) { - //Calculate the HV sector corresponding to the cluster position - //Arguments: y - //Returns the HV sector in the single module + // Calculate the HV sector corresponding to the cluster position + // Arguments: y + // Returns the HV sector in the single module Int_t hvsec = -1; Int_t pc, px, py; @@ -535,15 +554,15 @@ Int_t Param::inHVSector(float y) double Param::findTemp(double tLow, double tHigh, double y) { // Model for gradient in temperature - double yRad = hinRad(y); //height in a given radiator + double yRad = hinRad(y); // height in a given radiator if (tHigh < tLow) { - tHigh = tLow; //if Tout < Tin consider just Tin as reference... + tHigh = tLow; // if Tout < Tin consider just Tin as reference... } if (yRad < 0) { - yRad = 0; //protection against fake y values + yRad = 0; // protection against fake y values } if (yRad > sizePcY()) { - yRad = sizePcY(); //protection against fake y values + yRad = sizePcY(); // protection against fake y values } double gradT = (tHigh - tLow) / sizePcY(); // linear gradient @@ -552,9 +571,9 @@ double Param::findTemp(double tLow, double tHigh, double y) //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Param::setChStatus(Int_t ch, bool status) { - //Set a chamber on or off depending on the status - //Arguments: ch=chamber,status=kTRUE = active, kFALSE=off - //Returns: none + // Set a chamber on or off depending on the status + // Arguments: ch=chamber,status=kTRUE = active, kFALSE=off + // Returns: none for (Int_t padx = 0; padx < kMaxPcx + 1; padx++) { for (Int_t pady = 0; pady < kMaxPcy + 1; pady++) { fgMapPad[padx][pady][ch] = status; @@ -564,10 +583,10 @@ void Param::setChStatus(Int_t ch, bool status) //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Param::setSectStatus(Int_t ch, Int_t sect, bool status) { - //Set a given sector sect for a chamber ch on or off depending on the status - //Sector=0,5 (6 sectors) - //Arguments: ch=chamber,sect=sector,status: kTRUE = active, kFALSE=off - //Returns: none + // Set a given sector sect for a chamber ch on or off depending on the status + // Sector=0,5 (6 sectors) + // Arguments: ch=chamber,sect=sector,status: kTRUE = active, kFALSE=off + // Returns: none Int_t npadsect = (kMaxPcy + 1) / 6; Int_t padSectMin = npadsect * sect; @@ -582,9 +601,9 @@ void Param::setSectStatus(Int_t ch, Int_t sect, bool status) //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Param::setPcStatus(Int_t ch, Int_t pc, bool status) { - //Set a given PC pc for a chamber ch on or off depending on the status - //Arguments: ch=chamber,pc=PC,status: kTRUE = active, kFALSE=off - //Returns: none + // Set a given PC pc for a chamber ch on or off depending on the status + // Arguments: ch=chamber,pc=PC,status: kTRUE = active, kFALSE=off + // Returns: none Int_t deltaX = pc % 2; Int_t deltaY = pc / 2; @@ -602,9 +621,9 @@ void Param::setPcStatus(Int_t ch, Int_t pc, bool status) //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Param::printChStatus(Int_t ch) { - //Print the map status of a chamber on or off depending on the status - //Arguments: ch=chamber - //Returns: none + // Print the map status of a chamber on or off depending on the status + // Arguments: ch=chamber + // Returns: none Printf(" "); Printf(" --------- C H A M B E R %d ---------------", ch); for (Int_t pady = kMaxPcy; pady >= 0; pady--) { @@ -624,9 +643,9 @@ void Param::printChStatus(Int_t ch) //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Param::setGeomAccept() { - //Set the real acceptance of the modules, due to ineficciency or hardware problems (up tp 1/6/2010) - //Arguments: none - //Returns: none + // Set the real acceptance of the modules, due to ineficciency or hardware problems (up tp 1/6/2010) + // Arguments: none + // Returns: none setSectStatus(0, 3, kFALSE); setSectStatus(4, 0, kFALSE); setSectStatus(5, 1, kFALSE); diff --git a/Detectors/HMPID/reconstruction/CMakeLists.txt b/Detectors/HMPID/reconstruction/CMakeLists.txt index c56e087ba555d..f833f2d50545b 100644 --- a/Detectors/HMPID/reconstruction/CMakeLists.txt +++ b/Detectors/HMPID/reconstruction/CMakeLists.txt @@ -10,16 +10,20 @@ # or submit itself to any jurisdiction. o2_add_library(HMPIDReconstruction + TARGETVARNAME targetName SOURCES src/Clusterer.cxx src/HmpidDecoder2.cxx src/HmpidEquipment.cxx src/CTFHelper.cxx + src/Recon.cxx src/CTFCoder.cxx PUBLIC_LINK_LIBRARIES O2::HMPIDBase O2::DataFormatsHMP + ROOT::Physics O2::HMPIDSimulation) - + o2_target_root_dictionary(HMPIDReconstruction HEADERS include/HMPIDReconstruction/Clusterer.h include/HMPIDReconstruction/HmpidDecoder2.h + include/HMPIDReconstruction/Recon.h include/HMPIDReconstruction/HmpidEquipment.h) diff --git a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/CTFHelper.h b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/CTFHelper.h index 094aa7f367b3f..a4e12e880e0e0 100644 --- a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/CTFHelper.h +++ b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/CTFHelper.h @@ -75,50 +75,12 @@ class CTFHelper return (I&)(*this); } - const I operator++(int) - { - auto res = *this; - ++mIndex; - return res; - } - const I& operator--() { mIndex--; return (I&)(*this); } - const I operator--(int) - { - auto res = *this; - --mIndex; - return res; - } - - const I& operator+=(difference_type i) - { - mIndex += i; - return (I&)(*this); - } - - const I operator+=(difference_type i) const - { - auto tmp = *const_cast(this); - return tmp += i; - } - - const I& operator-=(difference_type i) - { - mIndex -= i; - return (I&)(*this); - } - - const I operator-=(difference_type i) const - { - auto tmp = *const_cast(this); - return tmp -= i; - } - difference_type operator-(const I& other) const { return mIndex - other.mIndex; } difference_type operator-(size_t idx) const { return mIndex - idx; } @@ -133,8 +95,6 @@ class CTFHelper bool operator==(const I& other) const { return mIndex == other.mIndex; } bool operator>(const I& other) const { return mIndex > other.mIndex; } bool operator<(const I& other) const { return mIndex < other.mIndex; } - bool operator>=(const I& other) const { return mIndex >= other.mIndex; } - bool operator<=(const I& other) const { return mIndex <= other.mIndex; } protected: gsl::span mData{}; @@ -159,18 +119,6 @@ class CTFHelper } return 0; } - value_type operator[](difference_type i) const - { - size_t id = mIndex + i; - if (id) { - if (mData[id].getOrbit() == mData[id - 1].getOrbit()) { - return mData[id].getBc() - mData[id - 1].getBc(); - } else { - return mData[id].getBc(); - } - } - return 0; - } }; //_______________________________________________ @@ -180,11 +128,6 @@ class CTFHelper public: using _Iter::_Iter; value_type operator*() const { return mIndex ? mData[mIndex].getOrbit() - mData[mIndex - 1].getOrbit() : 0; } - value_type operator[](difference_type i) const - { - size_t id = mIndex + i; - return id ? mData[id].getOrbit() - mData[id - 1].getOrbit() : 0; - } }; //_______________________________________________ @@ -194,7 +137,6 @@ class CTFHelper public: using _Iter::_Iter; value_type operator*() const { return mData[mIndex].getNumberOfObjects(); } - value_type operator[](difference_type i) const { return mData[mIndex + i].getNumberOfObjects(); } }; //_______________________________________________ @@ -213,11 +155,6 @@ class CTFHelper { return (*mTrigStart)[mIndex] ? mData[mIndex].getCh() : mData[mIndex].getCh() - mData[mIndex - 1].getCh(); } - value_type operator[](difference_type i) const - { - size_t id = mIndex + i; - return (*mTrigStart)[id] ? mData[id].getCh() : mData[id].getCh() - mData[id - 1].getCh(); - } }; //_______________________________________________ @@ -226,7 +163,6 @@ class CTFHelper public: using _Iter::_Iter; value_type operator*() const { return mData[mIndex].getQ(); } - value_type operator[](difference_type i) const { return mData[mIndex + i].getQ(); } }; //_______________________________________________ @@ -235,7 +171,6 @@ class CTFHelper public: using _Iter::_Iter; value_type operator*() const { return mData[mIndex].getPh(); } - value_type operator[](difference_type i) const { return mData[mIndex + i].getPh(); } }; //_______________________________________________ @@ -244,7 +179,6 @@ class CTFHelper public: using _Iter::_Iter; value_type operator*() const { return mData[mIndex].getX(); } - value_type operator[](difference_type i) const { return mData[mIndex + i].getX(); } }; //_______________________________________________ @@ -253,7 +187,6 @@ class CTFHelper public: using _Iter::_Iter; value_type operator*() const { return mData[mIndex].getY(); } - value_type operator[](difference_type i) const { return mData[mIndex + i].getY(); } }; //<<< =========================== ITERATORS ======================================== diff --git a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h index 35ee3e311f363..fb2fabfd9ce8c 100644 --- a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h +++ b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Clusterer.h @@ -22,6 +22,8 @@ #include "SimulationDataFormat/MCTruthContainer.h" #include "SimulationDataFormat/MCCompLabel.h" +#include "TMatrixF.h" // ef: added + namespace o2 { @@ -47,7 +49,7 @@ class Clusterer static void Dig2Clu(gsl::span digs, std::vector& clus, float* pUserCut, bool isUnfold = kTRUE); // digits->clusters static void FormClu(Cluster& pClu, int pDig, gsl::span digs, TMatrixF& pDigMap); // cluster formation recursive algorithm static int UseDig(int padX, int padY, TMatrixF& pDigMap); // use this pad's digit to form a cluster - inline bool IsDigSurvive(Digit* pDig) const; // check for sigma cut + inline bool IsDigSurvive(Digit* pDig) const; // check for sigma cut private: // void processChamber(std::vector& clusters, MCLabelContainer const* digitMCTruth); diff --git a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder.h b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder.h index 9392b2f61e8ff..be4b1d427b3fd 100644 --- a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder.h +++ b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder.h @@ -23,7 +23,7 @@ #include #include -#include +#include "FairLogger.h" #include "CommonDataFormat/InteractionRecord.h" #include "HMPIDReconstruction/HmpidEquipment.h" diff --git a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder2.h b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder2.h index 6dbb2cbd7fbbc..843fb4aafc98c 100644 --- a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder2.h +++ b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/HmpidDecoder2.h @@ -30,7 +30,7 @@ #include "DataFormatsHMP/Digit.h" -#include +#include "FairLogger.h" #include "HMPIDReconstruction/HmpidEquipment.h" diff --git a/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Recon.h b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Recon.h new file mode 100644 index 0000000000000..627f0ce98d3a5 --- /dev/null +++ b/Detectors/HMPID/reconstruction/include/HMPIDReconstruction/Recon.h @@ -0,0 +1,201 @@ +// CR statements + +#ifndef ALICEO2_HMPID_RECON_H +#define ALICEO2_HMPID_RECON_H + +/* +Changed from TCloneArrays of Cluster-pointers to vectors of clusters +changed par-names of cluster-pointers from pClu to cluster (name of cluster-object) +Changed name of clusters from pCluLst (TCloneArrays) to clusters (vector) +Changed raw-pointers to smart-pointers +Changed to cammelCase convention for AliceO2 member functions and variables +Changed from legacy physics classes TVector2 and TVector3 to Math/Vector3D +Changed from legacy physics classes TRotaiton to Rotation 3D +*/ + +#include //base class + + + +// ef : change from TRotation legacy class +#include +#include "Math/GenVector/RotationX.h" +#include "Math/GenVector/RotationY.h" +#include "Math/GenVector/RotationZ.h" + +// ef : new includes to replace TVector2/3 +#include +#include "Math/Vector3D.h" + +// ef : cartiesian Vector3D class XYZVector; // dirCkovTRS, dirCkovLORS + pc, rad in findPhotCkov + +// ef : polar Vector3D class Polar3DVector; // fTrkDir, dirTRS, dirLORS, dirCkov + refract -> dir + +#include +// ef: using vectors instead of TClonesArray + +// class AliESDtrack; //CkovAngle() ef : commented out +// ef: what is eq in O2? + +#include "HMPIDBase/Param.h" +class Param; +#include "DataFormatsHMP/Cluster.h" +#include "ReconstructionDataFormats/Track.h" + + +using Polar3DVector = ROOT::Math::Polar3DVector; + +namespace o2 +{ +namespace hmpid +{ +class Recon : public TNamed +{ + public: + using TrackParCov = o2::track::TrackParCov; + + // ef : moved Recon(): ctor from .cxx file + // Recon() = default; + + Recon() : TNamed("RichRec", "RichPat"), // ef : moved from cxx + fPhotCnt(-1), + fPhotFlag(0x0), + fPhotClusIndex(0x0), + fPhotCkov(0x0), + fPhotPhi(0x0), + fPhotWei(0x0), + fCkovSigma2(0), + fIsWEIGHT(kFALSE), + fDTheta(0.001), + fWindowWidth(0.045), + fRingArea(0), + fRingAcc(0), + fTrkDir(0, 0, 1), // Just for test + fTrkPos(30, 40), // Just for test + fMipPos(0, 0), + fPc(0, 0), + // fParam(o2::hmpid::Param::instance()) + fParam(Param::instance()) + { + //.. + // init of data members + //.. + + fParam->setRefIdx(fParam->meanIdxRad()); // initialization of ref index to a default one // ef:ok + } + + ~Recon() = default; + // virtual ~Recon() {;} //dtor + + // ef : methods in these classeS? + void initVars(int n); // init space for variables + + // ef : commented out: no need to delete variables when smart-pointer + // void deleteVars() const; // delete variables + + // void CkovAngle (AliESDtrack *pTrk,TClonesArray *pCluLst,int index,double nmean,float xRa,float yRa ); + void cKovAngle(TrackParCov trackParCov, const std::vector& clusters, int index, double nmean, float xRa, float yRa); // reconstructed Theta Cerenkov + + template + bool findPhotCkov(double cluX, double cluY, double& thetaCer, double& phiCer); // find ckov angle for single photon candidate + double findRingCkov(int iNclus); // best ckov for ring formed by found photon candidates + void findRingGeom(double ckovAng, int level = 1); // estimated area of ring in cm^2 and portion accepted by geometry + + template + const o2::math_utils::Vector2D intWithEdge(o2::math_utils::Vector2D p1, o2::math_utils::Vector2D p2); // find intercection between plane and lines of 2 thetaC + + // int flagPhot (double ckov,TClonesArray *pCluLst,AliESDtrack *pTrk ); + int flagPhot(double ckov, const std::vector& clusters, TrackParCov trackParCov); // is photon ckov near most probable track ckov + + double houghResponse(); // most probable track ckov angle + // template + void propagate(const Polar3DVector& dir, math_utils::Vector3D& pos, double z) const; // propagate photon alogn the line + // void refract(math_utils::Vector3D& dir, double n1, double n2) const; // refract photon on the boundary + + template + void refract(Polar3DVector& dir, double n1, double n2) const; // + + template + o2::math_utils::Vector2D tracePhot(double ckovTh, double ckovPh) const; // trace photon created by track to PC + + // ef : commented out addObjectToFriends + // void addObjectToFriends(TClonesArray *pCluLst, int photonIndex, AliESDtrack *pTrk ); // Add AliHMPIDCluster object to ESD friends + template + o2::math_utils::Vector2D traceForward(Polar3DVector dirCkov) const; // tracing forward a photon from (x,y) to PC + void lors2Trs(Polar3DVector dirCkov, double& thetaCer, double& phiCer) const; // LORS to TRS + void trs2Lors(math_utils::Vector3D dirCkov, double& thetaCer, double& phiCer) const; // TRS to LORS + + math_utils::Vector2D getMip() const + { + return fMipPos; + } // mip coordinates + + double getRingArea() const + { + return fRingArea; + } // area of the current ring in cm^2 + + double getRingAcc() const + { + return fRingAcc; + } // portion of the ring ([0,1]) accepted by geometry.To scale n. of photons + double findRingExt(double ckov, int ch, double xPc, double yPc, double thRa, double phRa); // find ring acceptance by external parameters + + void setTrack(double xRad, double yRad, double theta, double phi) + { + fTrkDir.SetCoordinates(1, theta, phi); + fTrkPos.SetCoordinates(xRad, yRad); + } // set track parameter at RAD + void setImpPC(double xPc, double yPc) + { + fPc.SetCoordinates(xPc, yPc); + } // set track impact to PC + void setMip(double xmip, double ymip) + { + fMipPos.SetCoordinates(xmip, ymip); + } // set track impact to PC + enum eTrackingFlags { kNotPerformed = -20, + kMipDistCut = -9, + kMipQdcCut = -5, + kNoPhotAccept = -11, + kNoRad = -22 }; + // + protected: + int fPhotCnt; // counter of photons candidate + + // ef : changed to smart-pointer arrays + std::unique_ptr fPhotFlag; // flags of photon candidates + std::unique_ptr fPhotClusIndex; // cluster index of photon candidates + std::unique_ptr fPhotCkov; // Ckov angles of photon candidates, [rad] + std::unique_ptr fPhotPhi; // phis of photons candidates, [rad] + std::unique_ptr fPhotWei; // weigths of photon candidates + + // int *fPhotClusIndex; // cluster index of photon candidates + + double fCkovSigma2; // sigma2 of the reconstructed ring + + bool fIsWEIGHT; // flag to consider weight procedure + float fDTheta; // Step for sliding window + float fWindowWidth; // Hough width of sliding window + + double fRingArea; // area of a given ring + double fRingAcc; // fraction of the ring accepted by geometry + + math_utils::Vector3D fTrkDir; // track direction in LORS at RAD + + math_utils::Vector2D fTrkPos; // track positon in LORS at RAD // XY mag + math_utils::Vector2D fMipPos; // mip positon for a given trackf // XY + math_utils::Vector2D fPc; // track position at PC // XY + + std::unique_ptr fParam; // Pointer to HMPIDParam + + private: + Recon(const Recon& r); // dummy copy constructor + Recon& operator=(const Recon& r); // dummy assignment operator + // + ClassDef(Recon, 3) +}; + +} // namespace hmpid +} // namespace o2 +#endif diff --git a/Detectors/HMPID/reconstruction/src/Clusterer.cxx b/Detectors/HMPID/reconstruction/src/Clusterer.cxx index 9ddacbb3de005..dc170369bd690 100644 --- a/Detectors/HMPID/reconstruction/src/Clusterer.cxx +++ b/Detectors/HMPID/reconstruction/src/Clusterer.cxx @@ -12,7 +12,7 @@ /// \file Clusterer.cxx /// \brief Implementation of the HMPID cluster finder #include -#include // for LOG +#include "FairLogger.h" // for LOG #include "Framework/Logger.h" #include "HMPIDBase/Param.h" #include "DataFormatsHMP/Cluster.h" @@ -63,7 +63,7 @@ void Clusterer::Dig2Clu(gsl::span digs, std::vector // for LOG +#include "FairLogger.h" // for LOG #include "Framework/Logger.h" #include "HMPIDReconstruction/HmpidDecodeRawFile.h" diff --git a/Detectors/HMPID/reconstruction/src/HmpidDecodeRawMem.cxx b/Detectors/HMPID/reconstruction/src/HmpidDecodeRawMem.cxx index 5a4f2acbfd97b..a1b4f6329afb0 100644 --- a/Detectors/HMPID/reconstruction/src/HmpidDecodeRawMem.cxx +++ b/Detectors/HMPID/reconstruction/src/HmpidDecodeRawMem.cxx @@ -18,7 +18,7 @@ /* ------ HISTORY --------- */ -#include // for LOG +#include "FairLogger.h" // for LOG #include "Framework/Logger.h" #include "DataFormatsHMP/Digit.h" diff --git a/Detectors/HMPID/reconstruction/src/HmpidDecoder.cxx b/Detectors/HMPID/reconstruction/src/HmpidDecoder.cxx index 93713811174b9..f680d00477908 100644 --- a/Detectors/HMPID/reconstruction/src/HmpidDecoder.cxx +++ b/Detectors/HMPID/reconstruction/src/HmpidDecoder.cxx @@ -19,7 +19,7 @@ /* ------ HISTORY --------- */ -#include // for LOG +#include "FairLogger.h" // for LOG #include "Framework/Logger.h" #include "Headers/RAWDataHeader.h" #include "HMPIDReconstruction/HmpidDecoder.h" diff --git a/Detectors/HMPID/reconstruction/src/HmpidDecoder2.cxx b/Detectors/HMPID/reconstruction/src/HmpidDecoder2.cxx index 6bbfb552f5154..0fa3a043ef4bc 100644 --- a/Detectors/HMPID/reconstruction/src/HmpidDecoder2.cxx +++ b/Detectors/HMPID/reconstruction/src/HmpidDecoder2.cxx @@ -19,7 +19,7 @@ /* ------ HISTORY --------- */ -#include // for LOG +#include "FairLogger.h" // for LOG #include "Framework/Logger.h" #include "Headers/RAWDataHeader.h" #include "HMPIDReconstruction/HmpidEquipment.h" @@ -1098,7 +1098,6 @@ void HmpidDecoder2::dumpHmpidError(int ErrorField) std::cout << "HMPID Decoder2 : [ERROR] " << "HMPID Error field = " << ErrorField << " : " << printbuf << std::endl; } - LOG(warn) << "HMPID Header Error Field =" << ErrorField << " [" << printbuf << "]"; } return; } diff --git a/Detectors/HMPID/reconstruction/src/README.m b/Detectors/HMPID/reconstruction/src/README.m new file mode 100644 index 0000000000000..a795c0d7a3127 --- /dev/null +++ b/Detectors/HMPID/reconstruction/src/README.m @@ -0,0 +1,38 @@ + + // in h : + // TVector3 posCkov(fTrkPos.X(), fTrkPos.Y(), zRad); + + + +Methods using SetPhi / SetMag / SetTheta must be Polar3D or CylindricalEta3D +// pc, rad in findPhotCkov are XYZ + +// fTrkDir, dirCkovTRS, dirCkovLORS XYZ + +// dirTRS, dirLORS, dirCkov + refract -> dir need Polar3D + + +dir should be ok now; but input to refract? (her er dir brukt) + + // theta: + // findPhotCkov -> fTrkDir + // lors2Trs -> dirCkovTRS + // trs2Lors -> dirCkovLORS + // findRingCkov -> fTrkDir + + // setTheta + // refract -> dir + + // Phi + // findPhotCkov -> pc, rad : Ikke inputs + // lors2Trs -> fTrkDir, dirCkovTRS : Ikke inputs + // trs2Lors -> fTrkDir, dirCkovLORS : Ikke inputs + // intWithEdge -> fTrkDir : Ikke input + + //SetMagThetaPhi -> SetR; SetTheta; SetPhi -> Only Polar3D + // dirTRS, dirLORS, dirCkov + + // in h : + // dir.SetXYZ(-999, -999, -999); + // dir.SetTheta(TMath::ASin(sinref)); + // Theta, Phi, SetMagThetaPhi diff --git a/Detectors/HMPID/reconstruction/src/Recon.cxx b/Detectors/HMPID/reconstruction/src/Recon.cxx new file mode 100644 index 0000000000000..7d4c28953d360 --- /dev/null +++ b/Detectors/HMPID/reconstruction/src/Recon.cxx @@ -0,0 +1,681 @@ + +////////////////////////////////////////////////////////////////////////// +// // +// HMPIDRecon // +// // +// HMPID class to perfom pattern recognition based on Hough transfrom // +// for single chamber // +////////////////////////////////////////////////////////////////////////// +#include "HMPIDBase/Param.h" +#include "HMPIDReconstruction/Recon.h" //class header + +#include //TracePhot() +#include //HoughResponse() +//#include //CkovAngle() ef : changed to std::vector + +//#include //CkovAngle() ef:? +//#include //CkovAngle() ef:? + +#include "ReconstructionDataFormats/Track.h" + +/* ef : + Changed from TCloneArrays of Cluster-pointers to vectors of clusters + changed par-names of cluster-pointers from pClu to cluster (name of cluster-object) + Changed name of clusters from pCluLst (TCloneArrays) to clusters (vector) +*/ + +// ef : moved isInDead and isDeadPad from Param.cxx to h +// because they are inline-static + +// ef : changed all functions to cammelcase convention per coding-guidelines + +// changed to smart-pointers +// not totally sure whether the initialization and deletion of the vars in initialize() and delete()-function is the best way to do it + +// ef : commented out all usage of AliESDtrack pTrk; +// changed AliESDtrack to TrackParCov (not sure if valid) + +// commented out addObjectToFriends + +// commented out deleteVars; not necessary to delete smart-pointers + +using namespace o2::hmpid; +// ClassImp(o2::hmpid::Recon); +ClassImp(o2::hmpid::Param); +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +// ef: moved to .h +/* +Recon::Recon(): + TNamed("RichRec","RichPat"), + fPhotCnt(-1), + fPhotFlag(0x0), + fPhotClusIndex(0x0), + fPhotCkov(0x0), + fPhotPhi(0x0), + fPhotWei(0x0), + fCkovSigma2(0), + fIsWEIGHT(kFALSE), + fDTheta(0.001), + fWindowWidth(0.045), + fRingArea(0), + fRingAcc(0), + fTrkDir(0,0,1), // Just for test + fTrkPos(30,40), // Just for test + fMipPos(0,0), + fPc(0,0), + fParam(o2::hmpid::Param::instance()) +{ +//.. +//init of data members +//.. + + fParam->setRefIdx(fParam->meanIdxRad()); // initialization of ref index to a default one +} */ +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Recon::initVars(int n) +{ + //.. + // Init some variables + //.. + if (n <= 0) { + return; + } + + // ef : changed to smart-pointer Array + // fPhotFlag = new int[n]; + fPhotFlag = std::unique_ptr(new int[n]); + fPhotClusIndex = std::unique_ptr(new int[n]); + + fPhotCkov = std::unique_ptr(new double[n]); + fPhotPhi = std::unique_ptr(new double[n]); + fPhotWei = std::unique_ptr(new double[n]); + // +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +// ef : commented out: no need to delete variables when smart-pointer +/*void Recon::deleteVars() const +{ + // ef: should not be done using this method? + // delete [] fPhotFlag; fPhotClusIndex; fPhotCkov; fPhotPhi;fPhotWei; +} */ + +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// void Recon::cKovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,int index,double nmean,float xRa,float yRa) +void Recon::cKovAngle(TrackParCov trackParCov, const std::vector& clusters, int index, double nmean, float xRa, float yRa) // ef : Cammelcase convention +{ + // Pattern recognition method based on Hough transform + // Arguments: pTrk - track for which Ckov angle is to be found + // pCluLst - list of clusters for this chamber + // Returns: - track ckov angle, [rad], + + const int nMinPhotAcc = 3; // Minimum number of photons required to perform the pattern recognition + + int nClusTot = clusters.size(); + + initVars(nClusTot); + + float xPc, yPc, th, ph; + + // ef : commented out: + // pTrk->GetHMPIDtrk(xPc,yPc,th,ph); //initialize this track: th and ph angles at middle of RAD + // ef : AliESDtrack::GetHMPIDtrk + + setTrack(xRa, yRa, th, ph); + + fParam->setRefIdx(nmean); + + float mipX = -1, mipY = -1; + int chId = -1, mipQ = -1, sizeClu = -1; + + fPhotCnt = 0; + + int nPads = 0; + + for (int iClu = 0; iClu < clusters.size(); iClu++) { // clusters loop + + o2::hmpid::Cluster cluster = clusters[iClu]; + nPads += clusters.size(); + if (iClu == index) { // this is the MIP! not a photon candidate: just store mip info + mipX = cluster.x(); + mipY = cluster.y(); + mipQ = (int)cluster.q(); + sizeClu = cluster.size(); + continue; + } + chId = cluster.ch(); + if (cluster.q() > 2 * fParam->qCut()) + continue; + double thetaCer, phiCer; + if (findPhotCkov(cluster.x(), cluster.y(), thetaCer, phiCer)) { // find ckov angle for this photon candidate + fPhotCkov[fPhotCnt] = thetaCer; // actual theta Cerenkov (in TRS) + fPhotPhi[fPhotCnt] = phiCer; + fPhotClusIndex[fPhotCnt] = iClu; // actual phi Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z) + fPhotCnt++; // increment counter of photon candidates + } + } // clusters loop + + // pTrk->SetHMPIDmip(mipX,mipY,mipQ,fPhotCnt); //store mip info in any case + // pTrk->SetHMPIDcluIdx(chId,index+1000*sizeClu); //set index of cluster + + if (fPhotCnt < nMinPhotAcc) { // no reconstruction with <=3 photon candidates + /* ef : commented out : + // pTrk->SetHMPIDsignal(kNoPhotAccept); //set the appropriate flag + */ + return; + } + + fMipPos.SetCoordinates(mipX, mipY); // ef: ok if Tvector2 can be used + + // PATTERN RECOGNITION STARTED: + if (fPhotCnt > fParam->multCut()) + fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction + else + fIsWEIGHT = kFALSE; + + int iNrec = flagPhot(houghResponse(), clusters, trackParCov); // flag photons according to individual theta ckov with respect to most probable + + /* ef : commented out : + // pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec); //store mip info + */ + if (iNrec < nMinPhotAcc) { // ef:ok + /* ef : commented out : + // pTrk->SetHMPIDsignal(kNoPhotAccept); //no photon candidates are accepted + */ + return; + } + + int occupancy = (int)(1000 * (nPads / (6. * 80. * 48.))); + + double thetaC = findRingCkov(clusters.size()); // find the best reconstructed theta Cherenkov + // FindRingGeom(thetaC,2); + + /* ef : commented out : + // pTrk->SetHMPIDsignal(thetaC+occupancy); //store theta Cherenkov and chmaber occupancy + // pTrk->SetHMPIDchi2(fCkovSigma2); //store experimental ring angular resolution squared + */ + + // deleteVars(); ef : in case of smart-pointers, should not be necessary? +} // CkovAngle() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\ + +template +bool Recon::findPhotCkov(double cluX, double cluY, double& thetaCer, double& phiCer) +{ + // Finds Cerenkov angle for this photon candidate + // Arguments: cluX,cluY - position of cadidate's cluster + // Returns: Cerenkov angle + + Polar3DVector dirCkov; // ef : TVector3->Polar3D + + double zRad = -0.5 * fParam->radThick() - 0.5 * fParam->winThick(); // z position of middle of RAD + + // ef : TVector3->XYZVector + math_utils::Vector3D rad(fTrkPos.X(), fTrkPos.Y(), zRad); // impact point at middle of RAD + math_utils::Vector3D pc(cluX, cluY, 0.5 * fParam->winThick() + fParam->gapIdx()); // mip at PC + + double cluR = TMath::Sqrt((cluX - fTrkPos.X()) * (cluX - fTrkPos.X()) + + (cluY - fTrkPos.Y()) * (cluY - fTrkPos.Y())); // ref. distance impact RAD-CLUSTER + + double phi = (pc - rad).Phi(); // phi of photon + + double ckov1 = 0; + double ckov2 = 0.75 + fTrkDir.Theta(); // start to find theta cerenkov in DRS + const double kTol = 0.01; + int iIterCnt = 0; + + while (1) { + if (iIterCnt >= 50) + return kFALSE; + double ckov = 0.5 * (ckov1 + ckov2); + + // ef SetMagThetaPhi -> SetCoordinates + dirCkov.SetCoordinates(1, ckov, phi); + o2::math_utils::Vector2D posC = traceForward(dirCkov); // trace photon with actual angles + double dist = cluR - (posC - fTrkPos).Mag2(); // get distance between trial point and cluster position // ef mod->Mag + // .Mod() Tvector2 + + if (posC.X() == -999) + dist = -999; // total reflection problem + iIterCnt++; // counter step + if (dist > kTol) + ckov1 = ckov; // cluster @ larger ckov + else if (dist < -kTol) + ckov2 = ckov; // cluster @ smaller ckov + else { // precision achived: ckov in DRS found + + // ef SetMagThetaPhi -> SetCoordinates + dirCkov.SetCoordinates(1, ckov, phi); // SetMagThetaPhi in TVecto3 + lors2Trs(dirCkov, thetaCer, phiCer); // find ckov (in TRS:the effective Cherenkov angle!) + return kTRUE; + } + } +} // FindPhotTheta() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +/*EF : COMMENTED traceForward */ +template // typename +o2::math_utils::Vector2D Recon::traceForward(Polar3DVector dirCkov) const +{ + // Trace forward a photon from (x,y) up to PC + // Arguments: dirCkov photon vector in LORS + // Returns: pos of traced photon at PC + + math_utils::Vector2D pos(-999, -999); + double thetaCer = dirCkov.Theta(); + if (thetaCer > TMath::ASin(1. / fParam->getRefIdx())) + return pos; // total refraction on WIN-GAP boundary + double zRad = -0.5 * fParam->radThick() - 0.5 * fParam->winThick(); // z position of middle of RAD + + // ef + math_utils::Vector3D posCkov(fTrkPos.X(), fTrkPos.Y(), zRad); + // TVector3 posCkov(fTrkPos.X(), fTrkPos.Y(), zRad); // RAD: photon position is track position @ middle of RAD + + propagate(dirCkov, posCkov, -0.5 * fParam->winThick()); // go to RAD-WIN boundary + refract(dirCkov, fParam->getRefIdx(), fParam->winIdx()); // RAD-WIN refraction + propagate(dirCkov, posCkov, 0.5 * fParam->winThick()); // go to WIN-GAP boundary + refract(dirCkov, fParam->winIdx(), fParam->gapIdx()); // WIN-GAP refraction + propagate(dirCkov, posCkov, 0.5 * fParam->winThick() + fParam->gapThick()); // go to PC + pos.SetCoordinates(posCkov.X(), posCkov.Y()); + return pos; +} // TraceForward() + +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Recon::lors2Trs(Polar3DVector dirCkov, double& thetaCer, double& phiCer) const +{ + // Theta Cerenkov reconstruction + // Arguments: dirCkov photon vector in LORS + // Returns: thetaCer of photon in TRS + // phiCer of photon in TRS + // TVector3 dirTrk; + // dirTrk.SetMagThetaPhi(1,fTrkDir.Theta(),fTrkDir.Phi()); -> dirTrk.SetCoordinates(1,fTrkDir.Theta(),fTrkDir.Phi()) + // double thetaCer = TMath::ACos(dirCkov*dirTrk); + + // TRotation mtheta; // TRotation-> Rotaiton3D : matrix 3x3 + // mtheta.RotateY(-fTrkDir.Theta()); ef : change to : + ROOT::Math::Rotation3D mtheta(ROOT::Math::RotationY(-fTrkDir.Theta())); + + // math_utils::Rotation3D mphi; + // mphi.RotateZ(-fTrkDir.Phi()); ef : change to : + + ROOT::Math::Rotation3D mphi(ROOT::Math::RotationZ(-fTrkDir.Phi())); + + ROOT::Math::Rotation3D mrot = mtheta * mphi; + + // ef : TVector3->Polar3D + math_utils::Vector3D dirCkovTRS; + dirCkovTRS = mrot * dirCkov; + phiCer = dirCkovTRS.Phi(); // actual value of the phi of the photon + thetaCer = dirCkovTRS.Theta(); // actual value of thetaCerenkov of the photon +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Recon::trs2Lors(math_utils::Vector3D dirCkov, double& thetaCer, double& phiCer) const +{ + // Theta Cerenkov reconstruction + // Arguments: dirCkov photon vector in TRS + // Returns: thetaCer of photon in LORS + // phiCer of photon in LORS + + // TRotation mtheta; + // mtheta.RotateY(fTrkDir.Theta()); ef : changed to : + + ROOT::Math::Rotation3D mtheta(ROOT::Math::RotationY(fTrkDir.Theta())); + + // TRotation mphi; + // mphi.RotateZ(fTrkDir.Phi()); ef : changed to : + + ROOT::Math::Rotation3D mphi(ROOT::Math::RotationZ(fTrkDir.Phi())); + + ROOT::Math::Rotation3D mrot = mphi * mtheta; + + // ef : TVector3->Polar3D + math_utils::Vector3D dirCkovLORS; + dirCkovLORS = mrot * dirCkov; + phiCer = dirCkovLORS.Phi(); // actual value of the phi of the photon + thetaCer = dirCkovLORS.Theta(); // actual value of thetaCerenkov of the photon +} +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void Recon::findRingGeom(double ckovAng, int level) +{ + // Find area covered in the PC acceptance + // Arguments: ckovAng - cerenkov angle + // level - precision in finding area and portion of ring accepted (multiple of 50) + // Returns: area of the ring in cm^2 for given theta ckov + + int kN = 50 * level; + int nPoints = 0; + double area = 0; + + bool first = kFALSE; + + // this needs to be changed? + // TVector2 + math_utils::Vector2D pos1; + + for (int i = 0; i < kN; i++) { + if (!first) { + pos1 = tracePhot(ckovAng, double(TMath::TwoPi() * (i + 1) / kN)); // find a good trace for the first photon + if (pos1.X() == -999) + continue; // no area: open ring + + if (!fParam->isInside(pos1.X(), pos1.Y(), 0)) { + pos1 = intWithEdge(fMipPos, pos1); // find the very first intersection... + } else { + if (!Param::isInDead(1.0f, 1.0f)) // ef : moved method from Param.cxx to h + nPoints++; // photon is accepted if not in dead zone + } + first = kTRUE; + continue; + } + math_utils::Vector2D pos2 = tracePhot(ckovAng, double(TMath::TwoPi() * (i + 1) / kN)); // trace the next photon + if (pos2.X() == -999) + continue; // no area: open ring + if (!fParam->isInside(pos2.X(), pos2.Y(), 0)) { + pos2 = intWithEdge(fMipPos, pos2); + } else { + if (!Param::isInDead(pos2.X(), pos2.Y())) + nPoints++; // photon is accepted if not in dead zone + } + + area += TMath::Abs((pos1 - fMipPos).X() * (pos2 - fMipPos).Y() - (pos1 - fMipPos).Y() * (pos2 - fMipPos).X()); // add area of the triangle... + pos1 = pos2; + } + //--- find area and length of the ring; + fRingAcc = (double)nPoints / (double)kN; + area *= 0.5; + fRingArea = area; +} // FindRingGeom() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +template // typename +const o2::math_utils::Vector2D Recon::intWithEdge(o2::math_utils::Vector2D p1, math_utils::Vector2D p2) +{ + // It finds the intersection of the line for 2 points traced as photons + // and the edge of a given PC + // Arguments: 2 points obtained tracing the photons + // Returns: intersection point with detector (PC) edges + + double xmin = (p1.X() < p2.X()) ? p1.X() : p2.X(); + double xmax = (p1.X() < p2.X()) ? p2.X() : p1.X(); + double ymin = (p1.Y() < p2.Y()) ? p1.Y() : p2.Y(); + double ymax = (p1.Y() < p2.Y()) ? p2.Y() : p1.Y(); + + double m = TMath::Tan((p2 - p1).Phi()); + math_utils::Vector2D pint; + // intersection with low X + pint.SetCoordinates((double)(p1.X() + (0 - p1.Y()) / m), 0.); + if (pint.X() >= 0 && pint.X() <= fParam->sizeAllX() && + pint.X() >= xmin && pint.X() <= xmax && + pint.Y() >= ymin && pint.Y() <= ymax) + return pint; + // intersection with high X + pint.SetCoordinates((double)(p1.X() + (fParam->sizeAllY() - p1.Y()) / m), (double)(fParam->sizeAllY())); + if (pint.X() >= 0 && pint.X() <= fParam->sizeAllX() && + pint.X() >= xmin && pint.X() <= xmax && + pint.Y() >= ymin && pint.Y() <= ymax) + return pint; + // intersection with left Y + pint.SetCoordinates(0., (double)(p1.Y() + m * (0 - p1.X()))); + if (pint.Y() >= 0 && pint.Y() <= fParam->sizeAllY() && + pint.Y() >= ymin && pint.Y() <= ymax && + pint.X() >= xmin && pint.X() <= xmax) + return pint; + // intersection with righ Y + pint.SetCoordinates((double)(fParam->sizeAllX()), (double)(p1.Y() + m * (fParam->sizeAllX() - p1.X()))); // ef: Set->SetCoordinates + if (pint.Y() >= 0 && pint.Y() <= fParam->sizeAllY() && + pint.Y() >= ymin && pint.Y() <= ymax && + pint.X() >= xmin && pint.X() <= xmax) + return pint; + return p1; +} // IntWithEdge() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +double Recon::findRingCkov(int) +{ + // Loops on all Ckov candidates and estimates the best Theta Ckov for a ring formed by those candidates. Also estimates an error for that Theat Ckov + // collecting errors for all single Ckov candidates thetas. (Assuming they are independent) + // Arguments: iNclus- total number of clusters in chamber for background estimation + // Return: best estimation of track Theta ckov + + double wei = 0.; + double weightThetaCerenkov = 0.; + + double ckovMin = 9999., ckovMax = 0.; + double sigma2 = 0; // to collect error squared for this ring + + for (int i = 0; i < fPhotCnt; i++) { // candidates loop + if (fPhotFlag[i] == 2) { + if (fPhotCkov[i] < ckovMin) + ckovMin = fPhotCkov[i]; // find max and min Theta ckov from all candidates within probable window + if (fPhotCkov[i] > ckovMax) + ckovMax = fPhotCkov[i]; + weightThetaCerenkov += fPhotCkov[i] * fPhotWei[i]; + wei += fPhotWei[i]; // collect weight as sum of all candidate weghts + + sigma2 += 1. / fParam->sigma2(fTrkDir.Theta(), fTrkDir.Phi(), fPhotCkov[i], fPhotPhi[i]); + } + } // candidates loop + + if (sigma2 > 0) + fCkovSigma2 = 1. / sigma2; + else + fCkovSigma2 = 1e10; + + if (wei != 0.) + weightThetaCerenkov /= wei; + else + weightThetaCerenkov = 0.; + return weightThetaCerenkov; +} // FindCkovRing() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// int Recon::flagPhot(double ckov,TClonesArray *pCluLst, AliESDtrack *pTrk) +int Recon::flagPhot(double ckov, const std::vector& clusters, TrackParCov trackParCov) +{ + // Flag photon candidates if their individual ckov angle is inside the window around ckov angle returned by HoughResponse() + // Arguments: ckov- value of most probable ckov angle for track as returned by HoughResponse() + // Returns: number of photon candidates happened to be inside the window + + // Photon Flag: Flag = 0 initial set; + // Flag = 1 good candidate (charge compatible with photon); + // Flag = 2 photon used for the ring; + + int steps = (int)((ckov) / fDTheta); // how many times we need to have fDTheta to fill the distance between 0 and thetaCkovHough + + double tmin = (double)(steps - 1) * fDTheta; + double tmax = (double)(steps)*fDTheta; + double tavg = 0.5 * (tmin + tmax); + + tmin = tavg - 0.5 * fWindowWidth; + tmax = tavg + 0.5 * fWindowWidth; + + int iInsideCnt = 0; // count photons which Theta ckov inside the window + for (int i = 0; i < fPhotCnt; i++) { // photon candidates loop + fPhotFlag[i] = 0; + if (fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax) { + fPhotFlag[i] = 2; + // ef: comment + // addObjectToFriends(clusters, i, trackParCov); + iInsideCnt++; + } + } + + return iInsideCnt; + +} // FlagPhot() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +// ef : commented out addObjectToFriends +// void Recon::addObjectToFriends(TClonesArray *pCluLst, int photonIndex, AliESDtrack *pTrk) +/* +{ + // Add AliHMPIDcluster object to ESD friends + + o2::hmpid::Cluster cluster = clusters[fPhotClusIndex[photonIndex]]; + + // o2::hmpid::Cluster *pClu=(o2::hmpid::Cluster*)pCluLst->UncheckedAt(fPhotClusIndex[photonIndex]); + + // o2::hmpid::Cluster *pClus = new o2::hmpid::Cluster(*pClu); // ef : old + cluster.setChi2(fPhotCkov[photonIndex]); + // pTrk->AddCalibObject(pClus); // AliESDtrack +} */ +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +template // typename +o2::math_utils::Vector2D Recon::tracePhot(double ckovThe, double ckovPhi) const +{ + // Trace a single Ckov photon from emission point somewhere in radiator up to photocathode taking into account ref indexes of materials it travereses + // Arguments: ckovThe,ckovPhi- photon ckov angles in TRS, [rad] + // Returns: distance between photon point on PC and track projection + + double theta, phi; + math_utils::Vector3D dirTRS; // ef TVector3 -> Polar3D + + Polar3DVector dirLORS; + + // ef SetMagThetaPhi->SetCoordinates + + dirTRS.SetCoordinates(1, ckovThe, ckovPhi); // photon in TRS + trs2Lors(dirTRS, theta, phi); + dirLORS.SetCoordinates(1, theta, phi); // photon in LORS + return traceForward(dirLORS); // now foward tracing +} // TracePhot() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +// template +void Recon::propagate(const Polar3DVector& dir, math_utils::Vector3D& pos, double z) const +{ + // Finds an intersection point between a line and XY plane shifted along Z. + // Arguments: dir,pos - vector along the line and any point of the line + // z - z coordinate of plain + // Returns: none + // On exit: pos is the position if this intesection if any + static math_utils::Vector3D nrm(0, 0, 1); + math_utils::Vector3D pnt(0, 0, z); + + math_utils::Vector3D diff = pnt - pos; + double sint = 0; //(nrm * diff) / (nrm * dir); + pos += sint * dir; +} // Propagate() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +template // typename +void Recon::refract(Polar3DVector& dir, double n1, double n2) const +{ + // Refract direction vector according to Snell law + // Arguments: + // n1 - ref idx of first substance + // n2 - ref idx of second substance + // Returns: none + // On exit: dir is new direction + double sinref = (n1 / n2) * TMath::Sin(dir.Theta()); + if (TMath::Abs(sinref) > 1.) + dir.SetXYZ(-999, -999, -999); // dette er ok! + else + dir.SetTheta(TMath::ASin(sinref)); // dette er ok! +} + +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +double Recon::houghResponse() +{ + // fIdxMip = mipId; + + double kThetaMax = 0.75; + int nChannels = (int)(kThetaMax / fDTheta + 0.5); + + // ef : change to smart-pointer + + std::unique_ptr phots, photsw, resultw; + phots.reset(new TH1D("Rphot", "phots", nChannels, 0, kThetaMax)); + photsw.reset(new TH1D("RphotWeighted", "photsw", nChannels, 0, kThetaMax)); + resultw.reset(new TH1D("resultw", "resultw", nChannels, 0, kThetaMax)); + + /* ef : changed from this: + // TH1D *resultw = new TH1D("resultw","resultw" ,nChannels,0,kThetaMax); + // TH1D *phots = new TH1D("Rphot" ,"phots" ,nChannels,0,kThetaMax); + // TH1D *photsw = new TH1D("RphotWeighted" ,"photsw" ,nChannels,0,kThetaMax); */ + + int nBin = (int)(kThetaMax / fDTheta); + int nCorrBand = (int)(fWindowWidth / (2 * fDTheta)); + + for (int i = 0; i < fPhotCnt; i++) { // photon cadidates loop + double angle = fPhotCkov[i]; + if (angle < 0 || angle > kThetaMax) + continue; + phots->Fill(angle); + int bin = (int)(0.5 + angle / (fDTheta)); + double weight = 1.; + if (fIsWEIGHT) { + double lowerlimit = ((double)bin) * fDTheta - 0.5 * fDTheta; + double upperlimit = ((double)bin) * fDTheta + 0.5 * fDTheta; + findRingGeom(lowerlimit); + double areaLow = getRingArea(); + findRingGeom(upperlimit); + double areaHigh = getRingArea(); + double diffArea = areaHigh - areaLow; + if (diffArea > 0) + weight = 1. / diffArea; + } + photsw->Fill(angle, weight); + fPhotWei[i] = weight; + } // photon candidates loop + + for (int i = 1; i <= nBin; i++) { + int bin1 = i - nCorrBand; + int bin2 = i + nCorrBand; + if (bin1 < 1) + bin1 = 1; + if (bin2 > nBin) + bin2 = nBin; + double sumPhots = phots->Integral(bin1, bin2); + if (sumPhots < 3) + continue; // if less then 3 photons don't trust to this ring + double sumPhotsw = photsw->Integral(bin1, bin2); + if ((double)((i + 0.5) * fDTheta) > 0.7) + continue; + resultw->Fill((double)((i + 0.5) * fDTheta), sumPhotsw); + } + // evaluate the "BEST" theta ckov as the maximum value of histogramm + + // ef : get() method should not be used to create new pointers for raw-pointers from smart-pointers, + // does this apply to the GetArray-method too? + double* pVec = resultw->GetArray(); + int locMax = TMath::LocMax(nBin, pVec); + + // ef: not this method, raw-pointers should not be used with new/delete-keywords + // smart-pointers are deleted when the fcuntion exits scope : + // delete phots;delete photsw;delete resultw; // Reset and delete objects + + return (double)(locMax * fDTheta + 0.5 * fDTheta); // final most probable track theta ckov +} // HoughResponse() +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +double Recon::findRingExt(double ckov, Int_t ch, double xPc, double yPc, double thRa, double phRa) +{ + // To find the acceptance of the ring even from external inputs. + // + // + double xRa = xPc - (fParam->radThick() + fParam->winThick() + fParam->gapThick()) * TMath::Cos(phRa) * TMath::Tan(thRa); // just linear extrapolation back to RAD + double yRa = yPc - (fParam->radThick() + fParam->winThick() + fParam->gapThick()) * TMath::Sin(phRa) * TMath::Tan(thRa); + + int nStep = 500; + int nPhi = 0; + + Int_t ipc, ipadx, ipady; + + if (ckov > 0) { + setTrack(xRa, yRa, thRa, phRa); + for (int j = 0; j < nStep; j++) { + math_utils::Vector2D pos; + pos = tracePhot(ckov, j * TMath::TwoPi() / (double)(nStep - 1)); + if (Param::isInDead(pos.X(), pos.Y())) // ef : moved method from Param.cxx to h + continue; // ef + fParam->lors2Pad(pos.X(), pos.Y(), ipc, ipadx, ipady); + ipadx += (ipc % 2) * fParam->kPadPcX; + ipady += (ipc / 2) * fParam->kPadPcY; + if (ipadx < 0 || ipady > 160 || ipady < 0 || ipady > 144 || ch < 0 || ch > 6) + continue; + if (Param::isDeadPad(ipadx, ipady, ch)) // ef : moved method from Param.cxx to h + continue; + nPhi++; + } // point loop + return ((double)nPhi / (double)nStep); + } // if + return -1; +} diff --git a/Detectors/HMPID/simulation/src/Detector.cxx b/Detectors/HMPID/simulation/src/Detector.cxx index 53a7656237938..fed7f9bd541b4 100644 --- a/Detectors/HMPID/simulation/src/Detector.cxx +++ b/Detectors/HMPID/simulation/src/Detector.cxx @@ -73,20 +73,20 @@ bool Detector::ProcessHits(FairVolume* v) Int_t volID = fMC->CurrentVolID(copy); auto stack = (o2::data::Stack*)fMC->GetStack(); - //Treat photons - //photon (Ckov or feedback) hits on module PC (Hpad) + // Treat photons + // photon (Ckov or feedback) hits on module PC (Hpad) if ((fMC->TrackPid() == 50000050 || fMC->TrackPid() == 50000051) && (volID == mHpad0VolID || volID == mHpad1VolID || volID == mHpad2VolID || volID == mHpad3VolID || volID == mHpad4VolID || volID == mHpad5VolID || volID == mHpad6VolID)) { - if (fMC->Edep() > 0) { //photon survided QE test i.e. produces electron + if (fMC->Edep() > 0) { // photon survided QE test i.e. produces electron if (IsLostByFresnel()) { fMC->StopTrack(); return false; - } //photon lost due to fersnel reflection on PC - Int_t tid = fMC->GetStack()->GetCurrentTrackNumber(); //take TID - Int_t pid = fMC->TrackPid(); //take PID - Float_t etot = fMC->Etot(); //total hpoton energy, [GeV] + } // photon lost due to fersnel reflection on PC + Int_t tid = fMC->GetStack()->GetCurrentTrackNumber(); // take TID + Int_t pid = fMC->TrackPid(); // take PID + Float_t etot = fMC->Etot(); // total hpoton energy, [GeV] Double_t x[3]; - fMC->TrackPosition(x[0], x[1], x[2]); //take MARS position at entrance to PC - Float_t hitTime = (Float_t)fMC->TrackTime(); //hit formation time + fMC->TrackPosition(x[0], x[1], x[2]); // take MARS position at entrance to PC + Float_t hitTime = (Float_t)fMC->TrackTime(); // hit formation time Int_t idch; // chamber number if (volID == mHpad0VolID) { idch = 0; @@ -104,20 +104,20 @@ bool Detector::ProcessHits(FairVolume* v) idch = 6; } Double_t xl, yl; - o2::hmpid::Param::instance()->mars2Lors(idch, x, xl, yl); //take LORS position - AddHit(x[0], x[1], x[2], hitTime, etot, tid, idch); //HIT for photon, position at P, etot will be set to Q - GenFee(etot); //generate feedback photons etot is modified in hit ctor to Q of hit + o2::hmpid::Param::instance()->mars2Lors(idch, x, xl, yl); // take LORS position + AddHit(x[0], x[1], x[2], hitTime, etot, tid, idch); // HIT for photon, position at P, etot will be set to Q + GenFee(etot); // generate feedback photons etot is modified in hit ctor to Q of hit stack->addHit(GetDetId()); - } //photon hit PC and DE >0 + } // photon hit PC and DE >0 return kTRUE; - } //photon hit PC + } // photon hit PC - //Treat charged particles - static Float_t eloss; //need to store mip parameters between different steps + // Treat charged particles + static Float_t eloss; // need to store mip parameters between different steps static Double_t in[3]; if (fMC->IsTrackEntering() && fMC->TrackCharge() && (volID == mHpad0VolID || volID == mHpad1VolID || volID == mHpad2VolID || volID == mHpad3VolID || volID == mHpad4VolID || volID == mHpad5VolID || volID == mHpad6VolID)) { - //Trackref stored when entering in the pad volume + // Trackref stored when entering in the pad volume o2::TrackReference tr(*fMC, GetDetId()); tr.setTrackID(stack->GetCurrentTrackNumber()); // tr.setUserId(lay); @@ -126,18 +126,18 @@ bool Detector::ProcessHits(FairVolume* v) if (fMC->TrackCharge() && (volID == mHcel0VolID || volID == mHcel1VolID || volID == mHcel2VolID || volID == mHcel3VolID || volID == mHcel4VolID || volID == mHcel5VolID || volID == mHcel6VolID)) { // charged particle in amplification gap (Hcel) - if (fMC->IsTrackEntering() || fMC->IsNewTrack()) { //entering or newly created - eloss = 0; //reset Eloss collector - fMC->TrackPosition(in[0], in[1], in[2]); //take position at the entrance - } else if (fMC->IsTrackExiting() || fMC->IsTrackStop() || fMC->IsTrackDisappeared()) { //exiting or disappeared - eloss += fMC->Edep(); //take into account last step Eloss - Int_t tid = fMC->GetStack()->GetCurrentTrackNumber(); //take TID - Int_t pid = fMC->TrackPid(); //take PID + if (fMC->IsTrackEntering() || fMC->IsNewTrack()) { // entering or newly created + eloss = 0; // reset Eloss collector + fMC->TrackPosition(in[0], in[1], in[2]); // take position at the entrance + } else if (fMC->IsTrackExiting() || fMC->IsTrackStop() || fMC->IsTrackDisappeared()) { // exiting or disappeared + eloss += fMC->Edep(); // take into account last step Eloss + Int_t tid = fMC->GetStack()->GetCurrentTrackNumber(); // take TID + Int_t pid = fMC->TrackPid(); // take PID Double_t out[3]; - fMC->TrackPosition(out[0], out[1], out[2]); //take MARS position at exit - Float_t hitTime = (Float_t)fMC->TrackTime(); //hit formation time + fMC->TrackPosition(out[0], out[1], out[2]); // take MARS position at exit + Float_t hitTime = (Float_t)fMC->TrackTime(); // hit formation time out[0] = 0.5 * (out[0] + in[0]); // - out[1] = 0.5 * (out[1] + in[1]); //take hit position at the anod plane + out[1] = 0.5 * (out[1] + in[1]); // take hit position at the anod plane out[2] = 0.5 * (out[2] + in[2]); Int_t idch; // chamber number if (volID == mHcel0VolID) { @@ -156,20 +156,20 @@ bool Detector::ProcessHits(FairVolume* v) idch = 6; } Double_t xl, yl; - o2::hmpid::Param::instance()->mars2Lors(idch, out, xl, yl); //take LORS position + o2::hmpid::Param::instance()->mars2Lors(idch, out, xl, yl); // take LORS position if (eloss > 0) { // HIT for MIP, position near anod plane, eloss will be set to Q AddHit(out[0], out[1], out[2], hitTime, eloss, tid, idch); - GenFee(eloss); //generate feedback photons + GenFee(eloss); // generate feedback photons stack->addHit(GetDetId()); eloss = 0; } } else { - //just going inside - eloss += fMC->Edep(); //collect this step eloss + // just going inside + eloss += fMC->Edep(); // collect this step eloss } return kTRUE; - } //MIP in GAP + } // MIP in GAP // later on return true if there was a hit! return false; @@ -185,16 +185,21 @@ void Detector::GenFee(Float_t qtot) { // Generate FeedBack photons for the current particle. To be invoked from StepManager(). // eloss=0 means photon so only pulse height distribution is to be analysed. + + // ef: how to implement bc of TVirtualMC* ? + // the method used (TrackPosition) takes TLorentzvector as input, but not Lorentzvector + + // ef: should use ROOT::Math::LorentzVector x4; TLorentzVector x4; fMC->TrackPosition(x4); Int_t iNphotons = fMC->GetRandom()->Poisson(0.02 * qtot); //# of feedback photons is proportional to the charge of hit - //AliDebug(1,Form("N photons=%i",iNphotons)); + // AliDebug(1,Form("N photons=%i",iNphotons)); Int_t j; Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop, dir[3], phi, pol[3], mom[4]; - //Generate photons - for (Int_t i = 0; i < iNphotons; i++) { //feedbacks loop + // Generate photons + for (Int_t i = 0; i < iNphotons; i++) { // feedbacks loop Double_t ranf[2]; - fMC->GetRandom()->RndmArray(2, ranf); //Sample direction + fMC->GetRandom()->RndmArray(2, ranf); // Sample direction cthf = ranf[0] * 2 - 1.0; if (cthf < 0) { continue; @@ -276,12 +281,13 @@ void Detector::GenFee(Float_t qtot) } fMC->Gdtom(pol, pol, 2); Int_t outputNtracksStored; - } //feedbacks loop -} //GenerateFeedbacks() + } // feedbacks loop +} // GenerateFeedbacks() //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Bool_t Detector::IsLostByFresnel() { // Calculate probability for the photon to be lost by Fresnel reflection. + // ROOT::Math::LorentzVector p4; TLorentzVector p4; Double_t mom[3], localMom[3]; fMC->TrackMomentum(p4); @@ -301,7 +307,7 @@ Bool_t Detector::IsLostByFresnel() } else { return kFALSE; } -} //IsLostByFresnel() +} // IsLostByFresnel() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Float_t Detector::Fresnel(Float_t ene, Float_t pdoti, Bool_t pola) { @@ -325,8 +331,8 @@ Float_t Detector::Fresnel(Float_t ene, Float_t pdoti, Bool_t pola) Float_t cn = csin[j] + ((csin[j + 1] - csin[j]) / 0.1) * (xe - en[j]); Float_t ck = csik[j] + ((csik[j + 1] - csik[j]) / 0.1) * (xe - en[j]); - //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR - //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197 + // FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR + // W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197 Float_t sinin = TMath::Sqrt((1. - pdoti) * (1. + pdoti)); Float_t tanin = sinin / pdoti; @@ -340,8 +346,8 @@ Float_t Detector::Fresnel(Float_t ene, Float_t pdoti, Bool_t pola) Float_t rp = rs * ((aO - sinin * tanin) * (aO - sinin * tanin) + b2) / ((aO + sinin * tanin) * (aO + sinin * tanin) + b2); - //CORRECTION FACTOR FOR SURFACE ROUGHNESS - //B.J. STAGG APPLIED OPTICS, 30(1991),4113 + // CORRECTION FACTOR FOR SURFACE ROUGHNESS + // B.J. STAGG APPLIED OPTICS, 30(1991),4113 Float_t sigraf = 18.; Float_t lamb = 1240 / ene; @@ -351,7 +357,7 @@ Float_t Detector::Fresnel(Float_t ene, Float_t pdoti, Bool_t pola) (4 * TMath::Pi() * pdoti * sigraf / lamb)); if (pola) { - Float_t pdotr = 0.8; //DEGREE OF POLARIZATION : 1->P , -1->S + Float_t pdotr = 0.8; // DEGREE OF POLARIZATION : 1->P , -1->S fresn = 0.5 * (rp * (1 + pdotr) + rs * (1 - pdotr)); } else { fresn = 0.5 * (rp + rs); @@ -359,7 +365,7 @@ Float_t Detector::Fresnel(Float_t ene, Float_t pdoti, Bool_t pola) fresn = fresn * rO; return fresn; -} //Fresnel() +} // Fresnel() //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ void Detector::Register() { FairRootManager::Instance()->RegisterAny(addNameTo("Hit").data(), mHits, true); } diff --git a/Detectors/HMPID/workflow/CMakeLists.txt b/Detectors/HMPID/workflow/CMakeLists.txt index 4599e692330e9..395d1f5001c5d 100644 --- a/Detectors/HMPID/workflow/CMakeLists.txt +++ b/Detectors/HMPID/workflow/CMakeLists.txt @@ -15,7 +15,8 @@ o2_add_library(HMPIDWorkflow src/DigitsToRawSpec.cxx src/DigitsToClustersSpec.cxx src/DigitsToRootSpec.cxx - src/ClustersToRootSpec.cxx + src/ClustersToRootSpec2.h # ef : use MakeTreeRootWriter + #src/ClustersToRootSpec.cxx src/DumpDigitsSpec.cxx src/PedestalsCalculationSpec.cxx src/RawToDigitsSpec.cxx @@ -79,10 +80,12 @@ o2_add_executable(digits-to-root-workflow SOURCES src/digits-to-root-workflow.cxx PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) -o2_add_executable(clusters-to-root-workflow - COMPONENT_NAME hmpid - SOURCES src/clusters-to-root-workflow.cxx - PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) + +# ef : not used anymore, digits-to-cluster performs writing to file or stream +#o2_add_executable(clusters-to-root-workflow +# COMPONENT_NAME hmpid +# SOURCES src/clusters-to-root-workflow.cxx +# PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) o2_add_executable(digits-to-raw-stream-workflow COMPONENT_NAME hmpid @@ -93,3 +96,4 @@ o2_add_executable(entropy-encoder-workflow COMPONENT_NAME hmpid SOURCES src/entropy-encoder-workflow.cxx PUBLIC_LINK_LIBRARIES O2::HMPIDWorkflow) + diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClustersToRootSpec.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClustersToRootSpec.h deleted file mode 100644 index 0e6d8d7ad0e5a..0000000000000 --- a/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClustersToRootSpec.h +++ /dev/null @@ -1,63 +0,0 @@ -// 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. - -/// -/// \file ClustersToRootSpec.h -/// \author Antonio Franco -/// -/// \brief Definition of a data processor to write Root File from Clusters stream -/// - -#ifndef DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_CLUSTERSTOROOTSPEC_H_ -#define DETECTORS_HMPID_WORKFLOW_INCLUDE_HMPIDWORKFLOW_CLUSTERSTOROOTSPEC_H_ - -#include "TTree.h" -#include "TFile.h" - -#include "Framework/DataProcessorSpec.h" -#include "Framework/Task.h" - -#include "HMPIDBase/Common.h" -#include "DataFormatsHMP/Cluster.h" -#include "DataFormatsHMP/Trigger.h" -#include "CommonDataFormat/InteractionRecord.h" - -#include "Framework/WorkflowSpec.h" - -namespace o2 -{ -namespace hmpid -{ - -class ClustersToRootTask : public framework::Task -{ - public: - ClustersToRootTask() = default; - ~ClustersToRootTask() override = default; - void init(framework::InitContext& ic) final; - void run(framework::ProcessingContext& pc) final; - void endOfStream(framework::EndOfStreamContext& ec) override; - - private: - ExecutionTimer mExTimer; - std::vector mTriggers; - std::vector mClusters; - TTree* mTheTree; - std::string mOutRootFileName; - TFile* mfileOut; -}; - -o2::framework::DataProcessorSpec getClustersToRootSpec(std::string inputSpec = "HMP/CLUSTERS"); - -} // end namespace hmpid -} // end namespace o2 - -#endif diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitReaderSpec.h_notused.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitReaderSpec.h_notused.h deleted file mode 100644 index eea9b134bd911..0000000000000 --- a/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitReaderSpec.h_notused.h +++ /dev/null @@ -1,53 +0,0 @@ -// 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. - -/// @file DigitReader.h - -#ifndef O2_HMPID_DIGITREADER -#define O2_HMPID_DIGITREADER - -#include "TFile.h" -#include "Framework/DataProcessorSpec.h" -#include "Framework/Task.h" -#include "DataFormatsHMP/Digit.h" -#include "SimulationDataFormat/MCCompLabel.h" -#include "SimulationDataFormat/MCTruthContainer.h" - -namespace o2 -{ -namespace hmpid -{ - -class DigitReader : public o2::framework::Task -{ - public: - DigitReader(bool useMC) : mUseMC(useMC) {} - ~DigitReader() override = default; - void init(o2::framework::InitContext& ic) final; - void run(o2::framework::ProcessingContext& pc) final; - - private: - int mState = 0; - bool mUseMC = true; - std::unique_ptr mFile = nullptr; - - std::vector mDigits, *mPdigits = &mDigits; - - o2::dataformats::MCTruthContainer mLabels, *mPlabels = &mLabels; -}; - -/// read simulated HMPID digits from a root file -framework::DataProcessorSpec getDigitReaderSpec(bool useMC); - -} // namespace hmpid -} // namespace o2 - -#endif /* O2_HMPID_DIGITREADER */ diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToClustersSpec.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToClustersSpec.h index 1f6473c95e7c2..6b0a2d3d3afe0 100644 --- a/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToClustersSpec.h +++ b/Detectors/HMPID/workflow/include/HMPIDWorkflow/DigitsToClustersSpec.h @@ -1,6 +1,6 @@ // 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. +// 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". @@ -22,9 +22,13 @@ #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" +#include "Framework/WorkflowSpec.h" #include "HMPIDBase/Common.h" #include "HMPIDReconstruction/Clusterer.h" -#include "Framework/WorkflowSpec.h" +#include "DataFormatsHMP/Trigger.h" + +#include "TFile.h" +#include "TTree.h" namespace o2 { @@ -34,25 +38,39 @@ namespace hmpid class DigitsToClustersTask : public framework::Task { public: - DigitsToClustersTask() = default; + DigitsToClustersTask(bool readFile) + : mReadFile(readFile) {} ~DigitsToClustersTask() override = default; void init(framework::InitContext& ic) final; void run(framework::ProcessingContext& pc) final; void endOfStream(framework::EndOfStreamContext& ec) override; private: + bool mReadFile = false; std::string mSigmaCutPar; float mSigmaCut[7] = {4, 4, 4, 4, 4, 4, 4}; - o2::hmpid::Clusterer* mRec; + std::unique_ptr mFile; ///< input file containin the tree + std::unique_ptr mTree; ///< input tree + + std::vector* mDigitsFromFile; + std::vector* mTriggersFromFile; + + std::unique_ptr mRec; // ef: changed to smart-pointer long mDigitsReceived; + void initFileIn(const std::string& fileName); + ExecutionTimer mExTimer; - void strToFloatsSplit(std::string s, std::string delimiter, float* res, int maxElem = 7); + void strToFloatsSplit(std::string s, std::string delimiter, float* res, + int maxElem = 7); }; -o2::framework::DataProcessorSpec getDigitsToClustersSpec(std::string inputSpec = "HMP/DIGITS"); +// ef : read from stream by default: +o2::framework::DataProcessorSpec + getDigitsToClustersSpec(std::string inputSpec = "HMP/DIGITS", bool readFile = false); +o2::framework::DataProcessorSpec getClustersToRootWriter(); } // end namespace hmpid } // end namespace o2 diff --git a/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClusterizerSpec.h_notused.h b/Detectors/HMPID/workflow/include/HMPIDWorkflow/HMPIDDigitizerSpec.h similarity index 71% rename from Detectors/HMPID/workflow/include/HMPIDWorkflow/ClusterizerSpec.h_notused.h rename to Detectors/HMPID/workflow/include/HMPIDWorkflow/HMPIDDigitizerSpec.h index 6102ec481c97c..c5992a26c5ff1 100644 --- a/Detectors/HMPID/workflow/include/HMPIDWorkflow/ClusterizerSpec.h_notused.h +++ b/Detectors/HMPID/workflow/include/HMPIDWorkflow/HMPIDDigitizerSpec.h @@ -9,8 +9,8 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#ifndef STEER_DIGITIZERWORKFLOW_HMPIDCLUSTERIZER_H_ -#define STEER_DIGITIZERWORKFLOW_HMPIDCLUSTERIZER_H_ +#ifndef STEER_DIGITIZERWORKFLOW_SRC_HMPIDDIGITIZERSPEC_H_ +#define STEER_DIGITIZERWORKFLOW_SRC_HMPIDDIGITIZERSPEC_H_ #include "Framework/DataProcessorSpec.h" @@ -19,9 +19,9 @@ namespace o2 namespace hmpid { -o2::framework::DataProcessorSpec getHMPIDClusterizerSpec(bool useMC); +o2::framework::DataProcessorSpec getHMPIDDigitizerSpec(int channel, bool mctruth = true); } // end namespace hmpid } // end namespace o2 -#endif /* STEER_DIGITIZERWORKFLOW_HMPIDCLUSTERIZERSPEC_H_ */ +#endif /* STEER_DIGITIZERWORKFLOW_SRC_HMPIDDIGITIZERSPEC_H_ */ diff --git a/Detectors/HMPID/workflow/src/ClustersToRootSpec.cxx b/Detectors/HMPID/workflow/src/ClustersToRootSpec.cxx deleted file mode 100644 index 6b39d74b41f87..0000000000000 --- a/Detectors/HMPID/workflow/src/ClustersToRootSpec.cxx +++ /dev/null @@ -1,153 +0,0 @@ -// 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. - -/// \file ClustersToRootSpec.cxx -/// \author Antonio Franco - INFN Bari -/// \version 1.0 -/// \date 20 nov 2021 -/// \brief Implementation of a data processor to produce Root File from Clusters Stream -/// - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "DPLUtils/DPLRawParser.h" -#include "DPLUtils/MakeRootTreeWriterSpec.h" - -#include "TTree.h" -#include "TFile.h" - -#include - -#include "Framework/DataRefUtils.h" -#include "Framework/InputSpec.h" -#include "Framework/CallbackService.h" -#include "Framework/ConfigParamRegistry.h" -#include "Framework/ControlService.h" -#include "Framework/DataProcessorSpec.h" -#include "Framework/Lifetime.h" -#include "Framework/Output.h" -#include "Framework/Task.h" -#include "Framework/WorkflowSpec.h" -#include "Framework/Logger.h" -#include "Framework/InputRecordWalker.h" - -#include "Headers/RAWDataHeader.h" -#include "DetectorsRaw/RDHUtils.h" - -#include "HMPIDBase/Geo.h" -#include "HMPIDWorkflow/ClustersToRootSpec.h" - -namespace o2 -{ -namespace hmpid -{ - -using namespace o2; -using namespace o2::framework; -using RDH = o2::header::RDHAny; - -//======================= -// Data decoder -void ClustersToRootTask::init(framework::InitContext& ic) -{ - LOG(info) << "[HMPID Write Root File From Clusters stream - init()]"; - - // get line command options - mOutRootFileName = ic.options().get("out-file"); - - mTriggers.clear(); - mClusters.clear(); - - TString filename; - TString tit; - - filename = TString::Format("%s", mOutRootFileName.c_str()); - LOG(info) << "Create the ROOT file " << filename.Data(); - mfileOut = new TFile(TString::Format("%s", filename.Data()), "RECREATE"); - tit = TString::Format("HMPID Clusters File Decoding"); - mTheTree = new TTree("o2hmp", tit); - mTheTree->Branch("InteractionRecords", &mTriggers); - mTheTree->Branch("HMPIDClusters", &mClusters); - - mExTimer.start(); - return; -} - -void ClustersToRootTask::run(framework::ProcessingContext& pc) -{ - std::vector triggers; - std::vector clusters; - - for (auto const& ref : InputRecordWalker(pc.inputs())) { - if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{header::gDataOriginHMP, "INTRECORDS1"}})) { - triggers = pc.inputs().get>(ref); - // LOG(info) << "We receive triggers =" << triggers.size(); - } - if (DataRefUtils::match(ref, {"check", ConcreteDataTypeMatcher{header::gDataOriginHMP, "CLUSTERS"}})) { - clusters = pc.inputs().get>(ref); - // LOG(info) << "The size of the vector =" << clusters.size(); - } - - for (int i = 0; i < triggers.size(); i++) { - // LOG(info) << "Trigger Event Orbit = " << triggers[i].getOrbit() << " BC = " << triggers[i].getBc(); - int startClusterIndex = mClusters.size(); - int numberOfClusters = 0; - for (int j = triggers[i].getFirstEntry(); j <= triggers[i].getLastEntry(); j++) { - mClusters.push_back(clusters[j]); // append the cluster - numberOfClusters++; - } - mTriggers.push_back(triggers[i]); - mTriggers.back().setDataRange(startClusterIndex, numberOfClusters); - } - } - mExTimer.stop(); - return; -} - -void ClustersToRootTask::endOfStream(framework::EndOfStreamContext& ec) -{ - mExTimer.logMes("Received an End Of Stream !"); - for (int i = 0; i < mClusters.size(); i++) { - mClusters.at(i).print(); - } - mTheTree->Fill(); - mTheTree->Write(); - mfileOut->Close(); - mExTimer.logMes("Register Tree ! "); - return; -} - -//_________________________________________________________________________________________________ -o2::framework::DataProcessorSpec getClustersToRootSpec(std::string inputSpec) -{ - std::vector inputs; - inputs.emplace_back("clusters", o2::header::gDataOriginHMP, "CLUSTERS", 0, Lifetime::Timeframe); - inputs.emplace_back("intrecord", o2::header::gDataOriginHMP, "INTRECORDS1", 0, Lifetime::Timeframe); - - std::vector outputs; - - return DataProcessorSpec{ - "HMP-ClustersToRoot", - inputs, - outputs, - AlgorithmSpec{adaptFromTask()}, - Options{{"out-file", VariantType::String, "hmpClusters.root", {"name of the output file"}}}}; -} - -} // namespace hmpid -} // end namespace o2 diff --git a/Detectors/HMPID/workflow/src/ClustersToRootSpec2.h b/Detectors/HMPID/workflow/src/ClustersToRootSpec2.h new file mode 100644 index 0000000000000..794d65725d1f8 --- /dev/null +++ b/Detectors/HMPID/workflow/src/ClustersToRootSpec2.h @@ -0,0 +1,48 @@ +// 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 STEER_DIGITIZERWORKFLOW_SRC_HMPCLUSTERWRITERSPEC_H_ +#define STEER_DIGITIZERWORKFLOW_SRC_HMPCLUSTERWRITERSPEC_H_ + +#include "Framework/DataProcessorSpec.h" +#include "DPLUtils/MakeRootTreeWriterSpec.h" +#include "Framework/InputSpec.h" +#include "DataFormatsHMP/Digit.h" +#include "DataFormatsHMP/Trigger.h" +#include "DataFormatsHMP/Cluster.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" + +namespace o2 +{ +namespace hmpid +{ + +template +using BranchDefinition = framework::MakeRootTreeWriterSpec::BranchDefinition; + +o2::framework::DataProcessorSpec getClustersToRootWriter() +{ + using InputSpec = framework::InputSpec; + using MakeRootTreeWriterSpec = framework::MakeRootTreeWriterSpec; + + return MakeRootTreeWriterSpec("HMPClustersWriter", + "hmpidclusters.root", + "o2sim", + 1, + BranchDefinition>{InputSpec{"hmpclusterinput", "HMP", "CLUSTERS"}, "HMPIDclusters"}, + BranchDefinition>{InputSpec{"hmpinteractionrecords", "HMP", "INTRECORDS1"}, "InteractionRecords"})(); +} + +} // end namespace hmpid +} // end namespace o2 + +#endif /* STEER_DIGITIZERWORKFLOW_SRC_HMPIDDIGITWRITERSPEC_H_ */ diff --git a/Detectors/HMPID/workflow/src/DigitsToClustersSpec.cxx b/Detectors/HMPID/workflow/src/DigitsToClustersSpec.cxx index 07da5f50a2ee4..e8263f57fffea 100644 --- a/Detectors/HMPID/workflow/src/DigitsToClustersSpec.cxx +++ b/Detectors/HMPID/workflow/src/DigitsToClustersSpec.cxx @@ -1,6 +1,6 @@ // 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. +// 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". @@ -16,36 +16,38 @@ /// \brief Implementation of a data processor to Clusterize the Digits /// -#include -#include -#include -#include #include +#include #include +#include +#include +#include #include -#include "HMPIDWorkflow/DigitsToClustersSpec.h" #include "Framework/CallbackService.h" #include "Framework/ConfigParamRegistry.h" #include "Framework/ControlService.h" #include "Framework/DataProcessorSpec.h" +#include "Framework/DataRefUtils.h" +#include "Framework/InputRecordWalker.h" #include "Framework/Lifetime.h" +#include "Framework/Logger.h" #include "Framework/Output.h" #include "Framework/Task.h" #include "Framework/WorkflowSpec.h" -#include "Framework/Logger.h" -#include "Framework/DataRefUtils.h" -#include "Framework/InputRecordWalker.h" +#include "HMPIDWorkflow/DigitsToClustersSpec.h" -#include "Headers/RAWDataHeader.h" -#include "DetectorsRaw/RDHUtils.h" #include "DPLUtils/DPLRawParser.h" +#include "DetectorsRaw/RDHUtils.h" +#include "Headers/RAWDataHeader.h" +#include "DataFormatsHMP/Cluster.h" #include "DataFormatsHMP/Digit.h" #include "DataFormatsHMP/Trigger.h" -#include "DataFormatsHMP/Cluster.h" #include "HMPIDBase/Geo.h" +#include "CommonUtils/NameConf.h" // ef : o2::utils::Str + namespace o2 { namespace hmpid @@ -56,8 +58,11 @@ using namespace o2::header; using namespace o2::framework; using RDH = o2::header::RDHAny; -// Splits a string in float array for string delimiter, TODO: Move this in a HMPID common library -void DigitsToClustersTask::strToFloatsSplit(std::string s, std::string delimiter, float* res, int maxElem) +// Splits a string in float array for string delimiter, TODO: Move this in a +// HMPID common library +void DigitsToClustersTask::strToFloatsSplit(std::string s, + std::string delimiter, float* res, + int maxElem) { int index = 0; size_t pos_start = 0; @@ -73,71 +78,187 @@ void DigitsToClustersTask::strToFloatsSplit(std::string s, std::string delimiter } } res[index++] = (std::stof(s.substr(pos_start))); + return; } //======================= // void DigitsToClustersTask::init(framework::InitContext& ic) { - LOG(info) << "[HMPID Clusterization - init() ] "; + LOG(info) << "[HMPID Clusterization - init() ; mReadFile = ] " + << mReadFile; mSigmaCutPar = ic.options().get("sigma-cut"); + if (mSigmaCutPar != "") { strToFloatsSplit(mSigmaCutPar, ",", mSigmaCut, 7); } mDigitsReceived = 0; - mRec = new o2::hmpid::Clusterer(); + mRec.reset(new o2::hmpid::Clusterer()); // ef: changed to smart-pointer + mExTimer.start(); + + // specify location and filename for output in case of writing to file + if (mReadFile) { + // Build the file name + const auto filename = o2::utils::Str::concat_string( + o2::utils::Str::rectifyDirectory( + ic.options().get("input-dir")), + ic.options().get("hmpid-digit-infile")); + initFileIn(filename); + } } void DigitsToClustersTask::run(framework::ProcessingContext& pc) { - LOG(debug) << "[HMPID DClusterization - run() ] Enter ..."; - auto triggers = pc.inputs().get>("intrecord"); - auto digits = pc.inputs().get>("digits"); + // outputs std::vector clusters; std::vector clusterTriggers; + LOG(info) << "[HMPID DClusterization - run() ] Enter ..."; + clusters.clear(); + clusterTriggers.clear(); - for (const auto& trig : triggers) { - if (trig.getNumberOfObjects()) { - gsl::span trigDigits{digits.data() + trig.getFirstEntry(), size_t(trig.getNumberOfObjects())}; - size_t clStart = clusters.size(); - mRec->Dig2Clu(trigDigits, clusters, mSigmaCut, true); - clusterTriggers.emplace_back(trig.getIr(), clStart, clusters.size() - clStart); + //===============mReadFromFile============================================= + if (mReadFile) { + LOG(info) << "[HMPID DClusterization - run() ] Entries = " << mTree->GetEntries(); + // check if more entries in tree + + if (mTree->GetReadEntry() + 1 >= mTree->GetEntries()) { + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + mExTimer.stop(); + mExTimer.logMes("End Clusterization ! digits = " + + std::to_string(mDigitsReceived)); } - } - LOGP(info, "Received {} triggers with {} digits -> {} triggers with {} clusters", triggers.size(), digits.size(), clusterTriggers.size(), clusters.size()); - mDigitsReceived += digits.size(); - pc.outputs().snapshot(o2::framework::Output{"HMP", "CLUSTERS", 0, o2::framework::Lifetime::Timeframe}, clusters); - pc.outputs().snapshot(o2::framework::Output{"HMP", "INTRECORDS1", 0, o2::framework::Lifetime::Timeframe}, clusterTriggers); + // there are more entries in file + else { + // =============== read digits and digit-triggers ===================== + // iterate through TTree + + auto entry = mTree->GetReadEntry() + 1; + assert(entry < mTree->GetEntries()); + mTree->GetEntry(0); + + // =============== create clusters ===================== + for (const auto& trig : *mTriggersFromFile) { + if (trig.getNumberOfObjects()) { + gsl::span trigDigits{ + mDigitsFromFile->data() + trig.getFirstEntry(), + size_t(trig.getNumberOfObjects())}; + size_t clStart = clusters.size(); + mRec->Dig2Clu(trigDigits, clusters, mSigmaCut, true); + clusterTriggers.emplace_back(trig.getIr(), clStart, + clusters.size() - clStart); + } + } + LOGP(info, "Received {} triggers with {} digits -> {} triggers with {} clusters", + mTriggersFromFile->size(), mDigitsFromFile->size(), clusterTriggers.size(), + clusters.size()); + mDigitsReceived += mDigitsFromFile->size(); + } // + } //=============== + + else { // ========= if readfromStream============================================== + + auto triggers = pc.inputs().get>("intrecord"); + auto digits = pc.inputs().get>("digits"); + + // =============== create clusters ===================== + for (const auto& trig : triggers) { + if (trig.getNumberOfObjects()) { + gsl::span trigDigits{ + digits.data() + trig.getFirstEntry(), + size_t(trig.getNumberOfObjects())}; + size_t clStart = clusters.size(); + mRec->Dig2Clu(trigDigits, clusters, mSigmaCut, true); + clusterTriggers.emplace_back(trig.getIr(), clStart, + clusters.size() - clStart); + } + } + } //========= + //===================================================================================== + + pc.outputs().snapshot( + o2::framework::Output{"HMP", "CLUSTERS", 0, + o2::framework::Lifetime::Timeframe}, + clusters); + pc.outputs().snapshot( + o2::framework::Output{"HMP", "INTRECORDS1", 0, + o2::framework::Lifetime::Timeframe}, + clusterTriggers); - mExTimer.elapseMes("Clusterization of Digits received = " + std::to_string(mDigitsReceived)); + mExTimer.elapseMes("Clusterization of Digits received = " + + std::to_string(mDigitsReceived)); } void DigitsToClustersTask::endOfStream(framework::EndOfStreamContext& ec) { + mExTimer.stop(); - mExTimer.logMes("End Clusterization ! digits = " + std::to_string(mDigitsReceived)); + mExTimer.logMes("End Clusterization ! digits = " + + std::to_string(mDigitsReceived)); +} + +void DigitsToClustersTask::initFileIn(const std::string& filename) +{ + // Create the TFIle + mFile = std::make_unique(filename.c_str(), "OLD"); + assert(mFile && !mFile->IsZombie()); + mTree.reset((TTree*)mFile->Get("o2sim")); + + if (!mTree) { + LOG(error) + << "HMPID DigitToClusterSpec::init() : Did not find o2sim tree in " + << filename.c_str(); + throw std::runtime_error( + "HMPID DigitToClusterSpec::init() : Did not find " + "o2sim file in digits tree"); + } + + mTree->SetBranchAddress("HMPDigit", &mDigitsFromFile); + mTree->SetBranchAddress("InteractionRecords", &mTriggersFromFile); + mTree->Print("toponly"); } //_________________________________________________________________________________________________ -o2::framework::DataProcessorSpec getDigitsToClustersSpec(std::string inputSpec) +o2::framework::DataProcessorSpec + getDigitsToClustersSpec(std::string inputSpec, bool readFile) { + + // define inputs if reading from stream: std::vector inputs; - inputs.emplace_back("digits", o2::header::gDataOriginHMP, "DIGITS", 0, Lifetime::Timeframe); - inputs.emplace_back("intrecord", o2::header::gDataOriginHMP, "INTRECORDS", 0, Lifetime::Timeframe); + if (!readFile) { + inputs.emplace_back("digits", o2::header::gDataOriginHMP, "DIGITS", 0, + Lifetime::Timeframe); + inputs.emplace_back("intrecord", o2::header::gDataOriginHMP, "INTRECORDS", + 0, Lifetime::Timeframe); + } + + // define outputs + + // outputs are streamed, and optionally stored in a root-file if the --write-to-file + // option in digits-to-clusters-workflow.cxx is passed std::vector outputs; - outputs.emplace_back("HMP", "CLUSTERS", 0, o2::framework::Lifetime::Timeframe); - outputs.emplace_back("HMP", "INTRECORDS1", 0, o2::framework::Lifetime::Timeframe); + + outputs.emplace_back("HMP", "CLUSTERS", 0, + o2::framework::Lifetime::Timeframe); + outputs.emplace_back("HMP", "INTRECORDS1", 0, + o2::framework::Lifetime::Timeframe); return DataProcessorSpec{ - "HMP-Clusterization", - inputs, - outputs, - AlgorithmSpec{adaptFromTask()}, - Options{{"sigma-cut", VariantType::String, "", {"sigmas as comma separated list"}}}}; + "HMP-Clusterization", inputs, outputs, + AlgorithmSpec{adaptFromTask(readFile)}, + Options{{"sigma-cut", + VariantType::String, + "", + {"sigmas as comma separated list"}}, + {"hmpid-digit-infile", + VariantType::String, + "hmpiddigits.root", + {"Name of the input file"}}, + {"input-dir", VariantType::String, "./", {"Input directory"}}}}; } } // namespace hmpid diff --git a/Detectors/HMPID/workflow/src/DigitsToRawSpec.cxx b/Detectors/HMPID/workflow/src/DigitsToRawSpec.cxx index bf6a8e7f34f11..21a5a0e5c7057 100644 --- a/Detectors/HMPID/workflow/src/DigitsToRawSpec.cxx +++ b/Detectors/HMPID/workflow/src/DigitsToRawSpec.cxx @@ -35,7 +35,7 @@ #include "Framework/Task.h" #include "Framework/WorkflowSpec.h" -#include // for LOG +#include "FairLogger.h" // for LOG #include "Framework/Logger.h" #include "Framework/InputRecordWalker.h" #include "DataFormatsParameters/GRPObject.h" @@ -197,7 +197,7 @@ o2::framework::DataProcessorSpec getDigitsToRawSpec() outputs, AlgorithmSpec{adaptFromTask()}, Options{{"outdir", VariantType::String, "./", {"base dir for output file"}}, - {"file-for", VariantType::String, "all", {"single file per: all,flp,link,crorcendpoint"}}, + {"file-for", VariantType::String, "all", {"single file per: all,flp,link,cru"}}, {"outfile", VariantType::String, "HMP", {"base name for output file"}}, {"in-file", VariantType::String, "hmpiddigits.root", {"name of the input sim root file"}}, {"dump-digits", VariantType::Bool, false, {"out the digits file in /tmp/hmpDumpDigits.dat"}}, diff --git a/Detectors/HMPID/workflow/src/EntropyEncoderSpec.cxx b/Detectors/HMPID/workflow/src/EntropyEncoderSpec.cxx index d9f032ec5ef7f..348a28cef0061 100644 --- a/Detectors/HMPID/workflow/src/EntropyEncoderSpec.cxx +++ b/Detectors/HMPID/workflow/src/EntropyEncoderSpec.cxx @@ -70,7 +70,7 @@ void EntropyEncoderSpec::run(ProcessingContext& pc) auto triggers = pc.inputs().get>("triggers"); auto digits = pc.inputs().get>("digits"); if (mSelIR) { - mCTFCoder.setSelectedIRFrames(pc.inputs().get>("selIRFrames")); + mCTFCoder.getIRFramesSelector().setSelectedIRFrames(pc.inputs().get>("selIRFrames")); } auto& buffer = pc.outputs().make>(Output{"HMP", "CTFDATA", 0, Lifetime::Timeframe}); auto iosize = mCTFCoder.encode(buffer, triggers, digits); @@ -104,8 +104,6 @@ DataProcessorSpec getEntropyEncoderSpec(bool selIR) {{"ctfrep"}, "HMP", "CTFENCREP", 0, Lifetime::Timeframe}}, AlgorithmSpec{adaptFromTask(selIR)}, Options{{"ctf-dict", VariantType::String, "ccdb", {"CTF dictionary: empty or ccdb=CCDB, none=no external dictionary otherwise: local filename"}}, - {"irframe-margin-bwd", VariantType::UInt32, 0u, {"margin in BC to add to the IRFrame lower boundary when selection is requested"}}, - {"irframe-margin-fwd", VariantType::UInt32, 0u, {"margin in BC to add to the IRFrame upper boundary when selection is requested"}}, {"mem-factor", VariantType::Float, 1.f, {"Memory allocation margin factor"}}}}; } diff --git a/Detectors/HMPID/workflow/src/HMPIDDigitizerSpec.cxx b/Detectors/HMPID/workflow/src/HMPIDDigitizerSpec.cxx new file mode 100644 index 0000000000000..50703fe558030 --- /dev/null +++ b/Detectors/HMPID/workflow/src/HMPIDDigitizerSpec.cxx @@ -0,0 +1,175 @@ +// 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 "HMPIDWorkflow/HMPIDDigitizerSpec.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/DataRefUtils.h" +#include "Framework/Lifetime.h" +#include "Headers/DataHeader.h" +#include "TStopwatch.h" +#include "Steer/HitProcessingManager.h" // for DigitizationContext +#include "TChain.h" +#include +#include +#include "Framework/Task.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsHMP/Digit.h" +#include "DataFormatsHMP/Trigger.h" +#include "HMPIDSimulation/HMPIDDigitizer.h" +#include "HMPIDSimulation/Detector.h" +#include "DetectorsBase/BaseDPLDigitizer.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include + +using namespace o2::framework; +using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; + + +namespace o2 +{ +namespace hmpid +{ + +class HMPIDDPLDigitizerTask : public o2::base::BaseDPLDigitizer +{ + public: + HMPIDDPLDigitizerTask() : o2::base::BaseDPLDigitizer(o2::base::InitServices::GEOM) {} + + void initDigitizerTask(framework::InitContext& ic) override + { + } + + void run(framework::ProcessingContext& pc) + { + static bool finished = false; + if (finished) { + return; + } + LOG(info) << "Doing HMPID digitization"; + + // read collision context from input + auto context = pc.inputs().get("collisioncontext"); + + context->initSimChains(o2::detectors::DetID::HMP, mSimChains); + + auto& irecords = context->getEventRecords(); + for (auto& record : irecords) { + LOG(info) << "HMPID TIME RECEIVED " << record.getTimeNS(); + } + + auto& eventParts = context->getEventParts(); + std::vector digitsAccum; // accumulator for digits + o2::dataformats::MCTruthContainer labelAccum; // timeframe accumulator for labels + mIntRecord.clear(); + + auto flushDigitsAndLabels = [this, &digitsAccum, &labelAccum]() { + // flush previous buffer + mDigits.clear(); + mLabels.clear(); + mDigitizer.flush(mDigits); + LOG(info) << "HMPID flushed " << mDigits.size() << " digits at this time "; + LOG(info) << "NUMBER OF LABEL OBTAINED " << mLabels.getNElements(); + int32_t first = digitsAccum.size(); // this is the first + std::copy(mDigits.begin(), mDigits.end(), std::back_inserter(digitsAccum)); + labelAccum.mergeAtBack(mLabels); + + // save info for the triggers accepted + LOG(info) << "Trigger Orbit :" << mDigitizer.getOrbit() << " BC:" << mDigitizer.getBc(); + mIntRecord.push_back(o2::hmpid::Trigger(o2::InteractionRecord(mDigitizer.getBc(), mDigitizer.getOrbit()), first, digitsAccum.size() - first)); + }; + + // loop over all composite collisions given from context + // (aka loop over all the interaction records) + for (int collID = 0; collID < irecords.size(); ++collID) { + // try to start new readout cycle by setting the trigger time + auto triggeraccepted = mDigitizer.setTriggerTime(irecords[collID].getTimeNS()); + if (triggeraccepted) { + flushDigitsAndLabels(); // flush previous readout cycle + } + auto withinactivetime = mDigitizer.setEventTime(irecords[collID].getTimeNS()); + if (withinactivetime) { + // for each collision, loop over the constituents event and source IDs + // (background signal merging is basically taking place here) + for (auto& part : eventParts[collID]) { + mDigitizer.setEventID(part.entryID); + mDigitizer.setSrcID(part.sourceID); + + // get the hits for this event and this source + std::vector hits; + context->retrieveHits(mSimChains, "HMPHit", part.sourceID, part.entryID, &hits); + LOG(info) << "For collision " << collID << " eventID " << part.entryID << " found HMP " << hits.size() << " hits "; + + mDigitizer.setLabelContainer(&mLabels); + mLabels.clear(); + mDigits.clear(); + + mDigitizer.process(hits, mDigits); + } + + } else { + LOG(info) << "COLLISION " << collID << "FALLS WITHIN A DEAD TIME"; + } + } + // final flushing step; getting everything not yet written out + flushDigitsAndLabels(); + + // send out to next stage + pc.outputs().snapshot(Output{"HMP", "DIGITS", 0, Lifetime::Timeframe}, digitsAccum); + pc.outputs().snapshot(Output{"HMP", "INTRECORDS", 0, Lifetime::Timeframe}, mIntRecord); + if (pc.outputs().isAllowed({"HMP", "DIGITLBL", 0})) { + pc.outputs().snapshot(Output{"HMP", "DIGITLBL", 0, Lifetime::Timeframe}, labelAccum); + } + LOG(info) << "HMP: Sending ROMode= " << mROMode << " to GRPUpdater"; + pc.outputs().snapshot(Output{"HMP", "ROMode", 0, Lifetime::Timeframe}, mROMode); + + // we should be only called once; tell DPL that this process is ready to exit + pc.services().get().readyToQuit(QuitRequest::Me); + finished = true; + } + + private: + HMPIDDigitizer mDigitizer; + std::vector mSimChains; + std::vector mDigits; + o2::dataformats::MCTruthContainer mLabels; // labels which get filled + std::vector mIntRecord; + + // RS: at the moment using hardcoded flag for continuous readout + o2::parameters::GRPObject::ROMode mROMode = o2::parameters::GRPObject::PRESENT; // readout mode +}; + +o2::framework::DataProcessorSpec getHMPIDDigitizerSpec(int channel, bool mctruth) +{ + // create the full data processor spec using + // a name identifier + // input description + // algorithmic description (here a lambda getting called once to setup the actual processing function) + // options that can be used for this processor (here: input file names where to take the hits) + std::vector outputs; + outputs.emplace_back("HMP", "DIGITS", 0, Lifetime::Timeframe); + outputs.emplace_back("HMP", "INTRECORDS", 0, Lifetime::Timeframe); + if (mctruth) { + outputs.emplace_back("HMP", "DIGITLBL", 0, Lifetime::Timeframe); + } + outputs.emplace_back("HMP", "ROMode", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "HMPIDDigitizer", + Inputs{InputSpec{"collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast(channel), Lifetime::Timeframe}}, + + outputs, + AlgorithmSpec{adaptFromTask()}, Options{}}; +} + +} // end namespace hmpid +} // end namespace o2 diff --git a/Detectors/HMPID/workflow/src/clusters-to-root-workflow.cxx b/Detectors/HMPID/workflow/src/clusters-to-root-workflow.cxx deleted file mode 100644 index 8da14f051dffd..0000000000000 --- a/Detectors/HMPID/workflow/src/clusters-to-root-workflow.cxx +++ /dev/null @@ -1,71 +0,0 @@ -// 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. - -/// \file clusters-to-root-workflow.cxx -/// \author Antonio Franco - INFN Bari -/// \version 1.0 -/// \date 22 nov 2021 -/// - -#include "Framework/WorkflowSpec.h" -#include "Framework/DataSpecUtils.h" -#include "Framework/CallbackService.h" -#include "Framework/ControlService.h" -#include "Framework/Task.h" -#include "Framework/CompletionPolicy.h" -#include "Framework/CompletionPolicyHelpers.h" -#include "Framework/DispatchPolicy.h" -#include "Framework/ConfigParamSpec.h" -#include "Framework/Variant.h" -#include "CommonUtils/ConfigurableParam.h" -#include "CommonUtils/NameConf.h" -#include "DetectorsRaw/HBFUtilsInitializer.h" -#include "Framework/CallbacksPolicy.h" - -void customize(std::vector& policies) -{ - o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); -} - -// customize the completion policy -void customize(std::vector& policies) -{ - using o2::framework::CompletionPolicy; - using o2::framework::CompletionPolicyHelpers; - policies.push_back(o2::framework::CompletionPolicyHelpers::defineByName("clusters-hmpid-root", CompletionPolicy::CompletionOp::Consume)); -} - -// we need to add workflow options before including Framework/runDataProcessing -void customize(std::vector& workflowOptions) -{ - std::string keyvaluehelp("Semicolon separated key=value strings ..."); - workflowOptions.push_back(o2::framework::ConfigParamSpec{"configKeyValues", o2::framework::VariantType::String, "", {keyvaluehelp}}); - o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); -} - -#include "Framework/runDataProcessing.h" -#include "HMPIDWorkflow/ClustersToRootSpec.h" - -using namespace o2; -using namespace o2::framework; - -WorkflowSpec defineDataProcessing(const ConfigContext& configcontext) -{ - WorkflowSpec specs; - o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); - DataProcessorSpec consumer = o2::hmpid::getClustersToRootSpec(); - specs.push_back(consumer); - - // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit - o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); - - return specs; -} diff --git a/Detectors/HMPID/workflow/src/digits-to-clusters-workflow.cxx b/Detectors/HMPID/workflow/src/digits-to-clusters-workflow.cxx index 513b6adcae70f..c8bdc2f9a2a77 100644 --- a/Detectors/HMPID/workflow/src/digits-to-clusters-workflow.cxx +++ b/Detectors/HMPID/workflow/src/digits-to-clusters-workflow.cxx @@ -1,6 +1,6 @@ // 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. +// 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". @@ -15,20 +15,26 @@ /// \date 22 nov 2021 /// -#include "Framework/WorkflowSpec.h" -#include "Framework/DataSpecUtils.h" +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/NameConf.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" #include "Framework/CallbackService.h" -#include "Framework/ControlService.h" -#include "Framework/Task.h" +#include "Framework/CallbacksPolicy.h" #include "Framework/CompletionPolicy.h" #include "Framework/CompletionPolicyHelpers.h" -#include "Framework/DispatchPolicy.h" #include "Framework/ConfigParamSpec.h" +#include "Framework/ControlService.h" +#include "Framework/DataSpecUtils.h" +#include "Framework/DispatchPolicy.h" +#include "Framework/Task.h" #include "Framework/Variant.h" -#include "CommonUtils/ConfigurableParam.h" -#include "CommonUtils/NameConf.h" -#include "DetectorsRaw/HBFUtilsInitializer.h" -#include "Framework/CallbacksPolicy.h" +#include "Framework/WorkflowSpec.h" + +/* + ef : perform clusterization: + either based on simulated data from a file, or on real data trhough a stream + The executable reads upstream and writes upstream by default. +*/ void customize(std::vector& policies) { @@ -40,19 +46,39 @@ void customize(std::vector& policies) { using o2::framework::CompletionPolicy; using o2::framework::CompletionPolicyHelpers; - policies.push_back(o2::framework::CompletionPolicyHelpers::defineByName("clusters-hmpid-write", CompletionPolicy::CompletionOp::Consume)); + policies.push_back(o2::framework::CompletionPolicyHelpers::defineByName( + "clusters-hmpid-write", CompletionPolicy::CompletionOp::Consume)); } // we need to add workflow options before including Framework/runDataProcessing void customize(std::vector& workflowOptions) { std::string keyvaluehelp("Semicolon separated key=value strings ..."); - workflowOptions.push_back(o2::framework::ConfigParamSpec{"configKeyValues", o2::framework::VariantType::String, "", {keyvaluehelp}}); + workflowOptions.push_back( + o2::framework::ConfigParamSpec{"configKeyValues", + o2::framework::VariantType::String, + "", + {keyvaluehelp}}); + workflowOptions.push_back( + o2::framework::ConfigParamSpec{"read-from-file", + o2::framework::VariantType::Bool, + false, + {"read upstream by default"}}); + workflowOptions.push_back( + o2::framework::ConfigParamSpec{"write-to-file", + o2::framework::VariantType::Bool, + false, + {"read upstream by default"}}); + o2::raw::HBFUtilsInitializer::addConfigOption(workflowOptions); } #include "Framework/runDataProcessing.h" #include "HMPIDWorkflow/DigitsToClustersSpec.h" +//#include "HMPIDWorkflow/ClustersToRootSpec.h" ef :no longer used, choice of +// writing to root/stream is defined here +#include "ClustersToRootSpec2.h" +//#include "HMPIDWorkflow/HMPIDDigitizerSpec.h" using namespace o2; using namespace o2::framework; @@ -60,8 +86,29 @@ using namespace o2::framework; WorkflowSpec defineDataProcessing(const ConfigContext& configcontext) { WorkflowSpec specs; - o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); - DataProcessorSpec consumer = o2::hmpid::getDigitsToClustersSpec(); + o2::conf::ConfigurableParam::updateFromString( + configcontext.options().get("configKeyValues")); + auto mFromFile = configcontext.options().get( + "read-from-file"); // read upstream by default + auto mToFile = configcontext.options().get( + "write-to-file"); // write upstream by default + + // HMPIDDigitizerSpec + /* + if(!mFromFile){ + DataProcessorSpec consumerFromFile = + o2::hmpid::getHMPIDDigitizerSpec(36, true); // int channel, bool mctruth + specs.push_back(consumerFromFile); + } */ + + DataProcessorSpec consumer = + o2::hmpid::getDigitsToClustersSpec("HMP/DIGITS", mFromFile); specs.push_back(consumer); + + if (mToFile) { // Write to File + DataProcessorSpec consumerClusterToRoot = o2::hmpid::getClustersToRootWriter(); + specs.push_back(consumerClusterToRoot); + } + return specs; }